Phase Retrieval via Matrix Completion

Introduction

This paper considers the fundamental problem of recovering a general signal, an image for example, from the magnitude of its Fourier transform. This problem, also known as phase retrieval, arises in many applications and has challenged engineers, physicists, and mathematicians for decades. Its origin comes from the fact that detectors can often times only record the squared modulus of the Fresnel or Fraunhofer diffraction pattern of the radiation that is scattered from an object. In such settings, one cannot measure the phase of the optical wave reaching the detector and, therefore, much information about the scattered object or the optical field is lost since, as is well known, the phase encodes a lot of the structural content of the image we wish to form.

Historically, the first application of phase retrieval is X-ray crystallography , and today this may still very well be the most important application. Over the last century or so, this field has developed a wide array of techniques to recover Bragg peaks from missing-phase data. Of course, the phase retrieval problem permeates many other areas of imaging science, and other applications include diffraction imaging , optics , astronomical imaging , microscopy , to name just a few. In particular, X-ray tomography has become an invaluable tool in biomedical imaging to generate quantitative 3D density maps of extended specimens on the nanoscale . Other subjects where phase retrieval plays an important role are quantum mechanics and even differential geometry . We note that phase retrieval has seen a resurgence in activity in recent years, fueled on the one hand by the desire to image individual molecules and other nano-particles, and on the other hand by new imaging capabilities: one such recent modality is the availability of new X-ray synchrotron sources that provide extraordinary X-ray fluxes, see for example . References and various instances of the phase retrieval problem as well as some theoretical and numerical solutions can be found in .

2 Main approaches to phase retrieval

Holographic techniques are among the more popular methods that have been proposed to measure the phase of the optical wave. While holographic techniques have been successfully applied in certain areas of optical imaging, they are generally difficult to implement in practice . Hence, the development of algorithms for signal recovery from magnitude measurements is still a very active field of research. Existing methods for phase retrieval rely on all kinds of a priori information about the signal, such as positivity, atomicity, support constraints, real-valuedness, and so on . Direct methods are limited in their applicability to small-scale problems due to their large computational complexity.

Oversampling in the Fourier domain has been proposed as a means to mitigate the non-uniqueness of the phase retrieval problem. While oversampling offers no benefit for most one-dimensional signals, the situation is more favorable for multidimensional signals, where it has been shown that twofold oversampling in each dimension almost always yields uniqueness for finitely supported, real-valued and non-negative signals . In other words, a digital image of the form $x=\{x[t_{1},t_{2}]\}$ with $0\leq t_{1}<n_{1}$ and $0\leq t_{2}<n_{2}$, whose Fourier transform is given by

is usually uniquely determined from the values of $|{\hat{x}}[\omega_{1},\omega_{2}]|$ on the oversampled grid $\omega=(\omega_{1},\omega_{2})\in\Omega=\Omega_{1}\times\Omega_{2}$ in which $\Omega_{i}=\{0,1/2,1,3/2,\ldots,n_{i}+1/2\}$. This holds provided $x$ has proper spatial support, is real valued and non-negative.

As pointed out in , these uniqueness results do not say anything about how a signal can be recovered from its intensity measurements, or about the robustness and stability of commonly used reconstruction algorithms — a fact we shall make very clear in the sequel. In fact, theoretical uniqueness conditions do not readily translate into numerical methods and as a result, the algorithmical and practical aspects of the phase retrieval problem (from noisy data) still pose significant challenges.

By and large, the most popular methods for phase retrieval from oversampled data are alternating projection algorithms pioneered by Gerchberg and Saxton and Fienup . These methods often require careful exploitation of signal constraints and delicate parameter selection to increase the likelihood of convergence to a correct solution . We describe the simplest realization of a widely used approach based on alternating projections , which assumes support constraints in the spatial domain and oversampled measurements in the frequency domain. With $T$ being a known subset containing the support of the signal $x$ ($\text{supp}(x)\subset T$) and Fourier magnitude measurements $\{y[\omega]\}_{\omega\in\Omega}$ with $y[\omega]=|\hat{x}[\omega]|$, the method works as follows:

Initialization: Choose an initial guess $x_{0}$ and set $z_{0}[\omega]=y[\omega]\frac{\hat{x}_{0}[\omega]}{|\hat{x}_{0}[\omega]|}$ for $\omega\in\Omega$.

Loop: For $k=1,2,\dots$ inductively define

While this algorithm is simple to implement and amenable to additional constraints such as the positivity of $x$, its convergence remains problematic. Projection algorithms onto convex sets are well understood . However, the set $\{z:|\hat{z}[\omega]|=|\hat{x}[\omega]|\}$ is not convex and, therefore, the algorithm is not known to converge in general or even to give a reasonable solution . Good results have been reported in certain settings but they appear to be nevertheless somewhat problematic in light of our numerical experiments from Section 4. Moreover, as discussed in , one of the most stringent limitations of these methods is the need for isolated objects (the support constraint). Finally, points out that oversampling is not always practically feasible as certain experimental geometries allow only for sub-Nyquist sampling; an example is the Bragg sampling from periodic crystalline structures.

In a different direction, a frame-theoretic approach to phase retrieval has been proposed in , where the authors derive various necessary and sufficient conditions for the uniqueness of the solution, as well as various numerical algorithms. While theoretically appealing, the practical applicability of these results is limited by the fact that very specific types of measurements are required, which cannot be realized in most applications of interest.

To summarize our discussion, we have seen many methods which all represent some important attempts to find efficient algorithms, and work well in certain situations. However, these techniques do not always provide a consistent and robust result.

3 PhaseLift – a novel methodology

This paper develops a novel methodology for phase retrieval based on a rigorous and flexible numerical framework. Whereas most of the existing methods seek to overcome nonuniqueness by imposing additional constraints on the signal, we pursue a different direction by assuming no constraints at all on the signal. There are two main components to our approach.

Multiple structured illuminations. We suggest collecting several diffraction patterns providing ‘different views’ of the sample or specimen. This can be accomplished in a number of ways: for instance, by modulating the light beam falling onto the sample or by placing a mask right after the sample, see Section 2 for details. Taking multiple diffraction patterns usually yields uniqueness as discussed in Section 3.

The concept of using multiple measurements as an attempt to resolve the phase ambiguity for diffraction imaging is of course not new, and was suggested in . Since then, a variety of methods have been proposed to carry out these multiple measurements; depending on the particular application, these may include the use of various gratings and/or of masks, the rotation of the axial position of the sample, and the use of defocusing implemented in a spatial light modulator, see for details and references. Other approaches include ptychography, an exciting field of research, where one records several diffraction patterns from overlapping areas of the sample, see and references therein.

Formulation of phase recovery as a matrix completion problem. We suggest (1) lifting up the problem of recovering a vector from quadratic constraints into that of a recovering of a rank-one matrix from affine constraints, and (2) relaxing the combinatorial problem into a convenient convex program. Since the lifting step is fundamental to our approach, we will refer to the proposed numerical framework as PhaseLift. The price we pay for trading the nonconvex quadratic constraints into convex constraints is that we must deal with a highly underdetermined problem. However, recent advances in the areas of compressive sensing and matrix completion have shown that such convex approximations are often exact.

Although our algorithmic framework appears to be novel for phase retrieval, the idea of solving problems involving nonconvex quadratic constraints by semidefinite relaxations has a long history in optimization, see and references therein, and Section 1.4 below for more discussion.

The goal of this paper is to demonstrate that taken together, multiple coded illuminations and convex programming (trace-norm minimization) provide a powerful new approach to phase retrieval. Further, a significant aspect of our methodology is that our systematic optimization framework offers a principled way of dealing with noise, and makes it easy to handle various statistical noise models. This is important because in practice, measurements are always noisy. In fact, our framework can be understood as an elaborate regularized maximum likelihood method. Lastly, our framework can also include a priori knowledge about the signal that can be formulated or relaxed as convex constraints.

4 Precedents

The idea of linearizing the phase retrieval problem by reformulating it as a problem of recovering a matrix from linear measurements can be found in . While this reference also contains some intriguing numerical recovery algorithms, their practical relevance for most applications is limited by the fact that the proposed measurement matrices either require a very specific algebraic structure which does not seem to be compatible with the physical properties of diffraction, or the number of measurements is proportional to the square of the signal dimension, which is not feasible in most applications.

In terms of framework, the closest approach is the paper , in which the authors use a matrix completion approach for array imaging from intensity measurements. Although this paper executes a similar relaxation as ours, there are some differences. We present a “noise-aware” framework, which makes it possible to account for a variety of noise models in a systematic way. Moreover, our emphasis is on a novel combination of structured illuminations and convex programming, which seems to bear great potential.

Methodology

Suppose $x=\{x[t]\}$ is the object of interest ($t$ may be a one- or multi-dimensional index). In this paper, we shall discuss illumination schemes collecting the diffraction pattern of the modulated object $w[t]x[t]$, where the waveforms or patterns $w[t]$ may be selected by the user. There are many ways in which this can be implemented in practice, and we discuss just a few of those.

Masking. One possibility is to modify the phase front after the sample by inserting a mask or a phase plate, see for example. A schematic layout is shown in Figure 1. In , the sample is scanned by shifting the phase plate as in ptychography (discussed below); the difference is that one scans the known phase plate rather than the object being imaged.

Optical grating. Another standard approach would be to change the profile or modulate the illuminating beam, which can easily be accomplished by the use of optical gratings . A simplified representation would look similar to the scheme depicted in Figure 1, with a grating instead of the mask (the grating could be placed before or after the sample).

Ptychography. Here, one measures multiple diffraction patterns by scanning a finite illumination on an extended specimen . In this setup, it is common to maintain a substantial overlap between adjacent illumination positions.

Oblique illuminations. One can use illuminating beams hitting the sample at user specified angle , see Figure 2 for a schematic illustration of this approach. One can also imagine having multiple simultaneous oblique illuminations.

As is clear, there is no shortage of options and one might prefer solutions which require generating as few diffraction patterns as possible for stable recovery.

2 Lifting

In the setting where we would collect the diffraction pattern of $w[t]x_{0}[t]$ as discussed earlier, then the waveform $a_{k}[t]$ can be written as

here, $\omega_{k}$ is a frequency value so that $a_{k}[t]$ is a patterned complex sinusoid. One can assume for convenience that the normalizing constant is such that $a_{k}$ is unit normed, i.e. $\|a_{k}\|_{2}^{2}=\sum_{t}|a_{k}[t]|^{2}=1$. Phase retrieval is then the feasibility problem

As is well known, quadratic measurements can be lifted up and interpreted as linear measurements about the rank-one matrix $X=xx^{*}$. Indeed,

where $A_{k}$ is the rank-one matrix $a_{k}a_{k}^{*}$. In what follows, we will let $\mathcal{A}$ be the linear operator mapping positive semidefinite matrices into $\{\operatorname{Tr}(A_{k}X):k=1,\ldots,m\}$. Hence, the phase retrieval problem is equivalent to

Upon solving the left-hand side of (2.3), we would factorize the rank-one solution $X$ as $xx^{*}$, hence finding solutions to the phase-retrieval problem. Note that the equivalence between the left- and right-hand side of (2.3) is straightforward since by definition, there exists a rank-one solution. Therefore, our problem is a rank minimization problem over an affine slice of the positive semidefinite cone. As such, it falls in the realm of low-rank matrix completion or matrix recovery, a class of optimization problems that has gained tremendous attention in recent years, see e.g. . Just as in matrix completion, the linear system $\mathcal{A}(X)=b$, with unknown in the positive semidefinite cone, is highly underdetermined. For instance suppose our signal $x_{0}$ has $n$ complex unknowns. Then we may imagine collecting six diffraction patterns with $n$ measurements for each (no oversampling). Thus $m=6n$ whereas the dimension of the space of $n\times n$ Hermitian matrices over the reals is $n^{2}$, which is obviously much larger.

We are of course interested in low-rank solutions and this makes the search feasible. This also raises an important question: what is the minimal number of diffraction patterns needed to recover $x$, whatever $x$ may be? Since each pattern yields $n$ real-valued coefficients and $x$ has $n$ complex-valued unknowns, the answer is at least two. Further, in the context of quantum state tomography, Theorem II in shows one needs at least $3n-2$ intensity measurements to guarantee uniqueness, hence suggesting an absolute minimum of three diffraction patterns. Are three patterns sufficient? For some answers to this question, see Section 3.

3 Recovery via convex programming

The rank minimization problem (2.3) is NP hard. We propose using the trace norm as a convex surrogate for the rank functional, giving the familiar SDP (and a crucial component of PhaseLift),

This problem is convex and there exists a wide array of numerical solvers including the popular Nesterov’s accelerated first order method . As far as the relationship between (2.3) and (2.4) is concerned, the matrix $\mathcal{A}$ in most diffraction imaging applications is not known to obey any of the conditions derived in the literature that would guarantee a formal equivalence between the two programs. Nevertheless, the formulation (2.4) enjoys great empirical performance as demonstrated in Section 4.

with optimization variables $\mu$ and $X$ (in other words, find $X\succeq 0$ such that $-\log p(b;\mathcal{A}(X))+\lambda\operatorname{Tr}(X)$ is minimum). Above, $\lambda$ is a positive scalar and, hence, our approach is a penalized or regularized maximum likelihood method, which trades off between goodness and complexity of the fit. When the likelihood is log-concave, problem (2.6) is convex and solvable. We give two examples for concreteness:

Poisson data. Suppose that $\{b_{k}\}$ is a sequence of independent samples from the Poisson distributions $\text{Poi}(\mu_{k})$. The Poisson log-likelihood for independent samples has the form $\sum_{k}b_{k}\log\mu_{k}-\mu_{k}$ (up to an additive constant factor) and thus, our problem becomes

Gaussian data. Suppose that $\{b_{k}\}$ is a sequence of independent samples from the Gaussian distribution with mean $\mu_{k}$ and variance $\sigma_{k}^{2}$ (or is well approximated by Gaussian variables). Then our problem becomes

If $\Sigma$ is a diagonal matrix with diagonal elements $\sigma_{k}^{2}$, this can be written as

Both formulations are of course convex and in both cases, one recovers the noiseless trace-minimization problem (2.4) as $\lambda\rightarrow 0^{+}$.

In addition, it is straightforward to include further constraints frequently discussed in the phase retrieval literature such as real-valuedness, positivity, atomicity and so on. Suppose the support of $x$ is known to be included in a set $T$ known a priori. Then we would add the linear constraint

(Algorithmically, one would simply work with a reduced-size matrix.) Suppose we would like to enforce real-valuedness, then we simply assume that $X$ is real valued and positive semidefinite. Finally positivity can be expressed as linear inequalities

Of course, many other types of constraints can be incorporated in this framework, which provides appreciable flexibility.

4 PhaseLift with reweighting

The trace norm promotes low-rank solutions and this is why it is often used as a convex proxy for the rank. However, it is possible to further promote low-rank solutions by solving a sequence of weighted trace-norm problems, a technique which has been shown to provide even more accurate solutions . The reweighting scheme works like this: choose $\varepsilon>0$; start with $W_{0}=I$ and for $k=0,1,\ldots,$ inductively define $X_{k}$ as the optimal solution to

The algorithm terminates on convergence or when the iteration count $k$ attains a specified maximum number of iterations $k_{\text{max}}$. One can see that the first step of this procedure is precisely (2.4); after this initial step, the algorithm proceeds in solving a sequence of trace-norm problems in which the matrix weights $W_{k}$ are roughly the inverse of the current guess.

As explained in the literature , this reweighting scheme can be viewed as attempting to solve

by minimizing the tangent approximation to $f$ at each iterate; that is to say, at step $k$, (2.7) is equivalent to minimizing $f(X_{k-1})+\langle\nabla f(X_{k-1}),X-X_{k-1}\rangle$ over the feasible set. This can also be applied to noise-aware variants where one would simply replace the objective functional in (2.6) with

at each step, and update $W_{k}$ in exactly the same way as before.

The log-det functional is closer to the rank functional than the trace norm. In fact, minimizing this functional solves the phase retrieval problem as incorporated in the following theorem.

Since the reweighting algorithm is a good heuristic for solving (2.8), we potentially have an interesting and tractable method for phase retrieval. It is not a perfect heuristic, however, as we cannot expect this procedure to always find the global minimum since the objective functional is concave.

The assumption that the identity matrix is in the span of the $A_{k}$’s holds whenever the modulus of the Fourier transform of the sample is measured. Indeed, if $|Fx|^{2}$ is observed, then letting $\{f_{k}^{*}\}$ be the rows of $F$, we have $\sum_{k}f_{k}f_{k}^{*}=I$.

Proof Our assumption implies that for any feasible $X$, $\operatorname{Tr}(X)=\sum_{k}h_{k}\operatorname{Tr}(A_{k}X)=\sum_{k}h_{k}b_{k}$ is fixed. Assume without loss of generality that feasible points obey $\operatorname{Tr}(X)=1$ (if $\operatorname{Tr}(X)=0$, then the unique solution is $X=0$). If $x_{0}$ is the unique solution to phase retrieval (up to global phase), then $X_{0}=x_{0}x_{0}^{*}$ is the only rank-one feasible point. We thus need to show that any feasible $X$ with $\operatorname{rank}(X)>1$, obeys $f(X)>f(X_{0})$, a fact which follows from the strong concavity of $f$ (of the logarithm). Let $X=\sum_{j}\lambda_{j}u_{j}u_{j}^{*}$ be any eigenvalue decomposition of a feasible point. Then

and it follows from the strict concavity of the log that

The first strict inequality holds unless $X$ is rank one, in which case, we have equality. The equality follows from $\sum_{j}\lambda_{j}=\operatorname{Tr}(X)=1$. Since the right-hand side is nothing else than $f(X_{0})$, the theorem is established.

For the second part, set $B_{k}=R^{-1}A_{k}R^{-1}$ and consider a new data problem with constraints $\{X:\operatorname{Tr}(B_{k}X)=b_{k},X\succeq 0\}$. Now $X$ is feasible for our problem if and only if $RXR$ is feasible for this new problem. (This is because the mapping $X\mapsto RXR$ preserves the positive semidefinite cone.) Now suppose $x_{0}$ is the solution to phase retrieval and set $X_{0}=x_{0}x_{0}^{*}$ as before. Since $I$ is in the span of the sensing matrices $B_{k}$, we have just learned that for all for all $RXR\neq RX_{0}R$ and $X$ feasible for our problem,

Theory

Our PhaseLift framework poses two main theoretical questions:

When do multiple diffracted images imply unicity of the solution?

When does our convex heuristic succeed in recovering the unique solution to the phase-retrieval problem?

Developing comprehensive answers to these questions constitutes a whole research program, which clearly is beyond the scope of this work. In this paper, we shall limit ourselves to introducing some theoretical results showing simple ways of designing diffraction patterns, which give unicity. Our focus is on getting uniqueness from a very limited number of diffraction patterns. For example, we shall demonstrate that in some cases three diffraction images are sufficient for perfect recovery. Thus, we give below partial answers to the first question and will address the second in a later publication.

A frequently discussed approach to retrieve phase information uses a technique from holography. Roughly speaking, the idea is to let the signal of interest $x$ interfere with a known reference beam $y$. One typically measures $|x+y|^{2}$ and $|x-iy|^{2}$ and precise knowledge of $y$ allows, in principle, to recover the amplitude and phase of $x$. Holographic techniques are hard to implement in practice. Instead, we propose using a modulated version of the signal itself as a reference beam which in some cases may be easier to implement.

These measurements can be obtained by illuminating the sample with the three light fields $1$, $e^{i2\pi st/n}$ and $e^{i2\pi(st/n-1/4)}$. We show below that these $3n$ measurements are generally sufficient for perfect recovery.

Conversely, if the DFT vanishes at two frequency points $k$ and $k^{\prime}$ obeying $k-k^{\prime}\neq s\text{ mod }n$, then recovery is not possible. from the $3n$ real numbers (3.1).

The proof of this theorem is constructive and we give a simple algorithm that achieves perfect reconstruction. Further, one can use masks to scramble the Fourier transform as to make sure it does not vanish. Suppose for instance that we collect

where the $z[t]$’s are iid $\mathcal{N}(0,1)$. Then since the Fourier transform of $z[t]x[t]$ does not vanish with probability one, we have the following corollary.

Of course, one could derive similar results by scrambling the Fourier transform with the aid of other types of masks, e.g. binary masks. We do not pursue such calculations.

$\mathcal{F}_{n_{1}\times n_{2}}$ is the 2D unitary Fourier transform defined by (1.2) in which the frequencies belong to the 2D grid $\{0,1,\ldots,n_{1}-1\}\times\{0,1,\ldots,n_{2}-1\}$, and $\mathcal{D}^{s}$ is the modulation

With these definitions, we have the following result:

Suppose the DFT of $x\in C^{n_{1}\times n_{2}}$ does not vanish. Then $x$ can be recovered up to global phase from the $3n_{1}n_{2}$ real numbers (3.2) if and only if $s_{1}$ is prime with $n_{1}$, $s_{2}$ is prime with $n_{2}$ and $n_{1}$ is prime with $n_{2}$. Under these assumptions, if the trace-minimization program (2.4) or the iteratively reweighted algorithm return a rank-1 solution, then this solution is exact.

Again, one can apply a random mask to turn this statement into a probabilistic statement holding either with probability one or with very large probability depending upon the mask that is used.

One can always choose $s_{1}$ and $s_{2}$ such that they be prime with $n_{1}$ and $n_{2}$ respectively. The last condition may be less friendly but one can decide to pad one dimension with zeros to guarantee primality. This is equivalent to a slight oversampling of the DFT along one direction. An alternative is to take $5n_{1}n_{2}$ measurements in which we modulate the signal horizontally and then vertically; that is to say, we modulate with $s=(s_{1},0)$ and then with $s=(0,s_{2})$. These $5n_{1}n_{2}$ measurements guarantee recovery if $s_{1}$ is prime with $n_{1}$ and $s_{2}$ is prime with $n_{2}$ for all sizes $n_{1}$ and $n_{2}$, see Section 3.3 for details.

for all $k\in\{0,1,\ldots,n-1\}$ (above, $k-s$ is understood mod $n$). Write $\hat{x}[k]=|\hat{x}[k]|e^{i\phi[k]}$ so that $\phi[k]$ is the missing phase, and observe that

Hence, if $\hat{x}[k]\neq 0$ for all $k\in\{0,1,\ldots,n-1\}$, our data gives us knowledge of all phase shifts of the form

We can, therefore, initialize $\phi$ to be zero and then get the values of $\phi[-s]$, $\phi[-2s]$ and so on.

For the second part of the theorem, assume without loss of generality, that $s=-1$ and that $(k,k^{\prime})=(0,k_{0})$ ($1<k_{0}<n-1$). For simplicity suppose these are the only zeros of the DFT. This creates two disjoint sets of frequency indices: those for which $0<k<k_{0}$ and those for which $k_{0}<k\leq n-1$. We are given no information about the phase difference between elements of these two subgroups, and hence recovery is not possible. This argument extends to situations where the DFT vanishes more often, in which case, we have even more indeterminacy.

2 Proof of Theorem 3.3

Let $\hat{x}=\{\hat{x}[k_{1},k_{2}]\}$, where $(k_{1},k_{2})\in\{0,1,\ldots n_{1}-1\}\times\{0,1,\ldots,n_{2}-1\}$ be the DFT of $x$. Then we have knowledge of

for all $(k_{1}k_{2})$. With the same notations as before, this gives knowledge of all phase shifts of the form

3 Extensions

It is clear from our analysis that if we were to collect $|\mathcal{F}_{n_{1}\times n_{2}}x|^{2}$ together with

A simple instance consists in choosing one modulation pattern to be $(s_{1},0)$ and another to be $(0,s_{2})$. If $s_{1}$ is prime is $n_{1}$ and $s_{2}$ with $n_{2}$, these two modulations generate the whole group regardless of the relationship between $n_{1}$ and $n_{2}$. An algorithmic way to see this is as follows. Initialize $\phi(0,0)$. Then by using horizontal modulations, one recovers all phases of the form $\phi(k_{1},0)$. Further, by using vertical modulations (starting with $\phi(k_{1},0)$), one can recover all phases of the form $\phi(k_{1},k_{2})$ by moving upward.

Numerical Experiments

This section introduces numerical simulations to illustrate and study the effectiveness of PhaseLift.

All numerical algorithms were implemented in Matlab using TFOCS as well as modifications of TFOCS template files. TFOCS is a library of Matlab-files designed to facilitate the construction of first-order methods for a variety of convex optimization problems, which include those we consider.

In a nutshell, suppose we wish to solve the problem

where $\{t_{k}\}$ is a sequence of stepsize rules and $\mathcal{P}$ is the projection onto the positive semidefinite cone. (Various stepsize rules are typically considered including fixed stepsizes, backtracking line search, exact line search and so on.)

TFOCS implements a variety of accelerated first-ordered methods pioneered by Nesterov, see and references therein. One variant works as follows. Choose $X_{0}$, set $Y_{0}=X_{0}$ and $\theta_{0}=1$, and inductively define

where $\{t_{k}\}$ is a sequence of stepsize rules as before. The sequence $\{\theta_{k}\}$ is usually referred to as a sequence of accelerated parameters, and $\{Y_{k}\}$ is an auxiliary sequence at which the gradient is to be evaluated. The advantage of this approach is that the computational work per iteration is as in the projected gradient method but the number of iterations needed to reach a certain accuracy is usually much lower . TFOCS implements such iterations and others like it but with various improvements.

For large problems, e.g. images with a large number $N$ of pixels, it is costly to hold the $N\times N$ optimization variable $X$ in memory. To overcome this issue, our computational approach maintains a low-rank factorization of $X$. This is achieved by substituting the projection onto the semidefinite cone (the expensive step) with a proxy. Whereas $\mathcal{P}$ dumps the negative eigenvalues as in

where $\sum_{i}\lambda_{i}u_{i}u_{i}^{*}$ ($\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{N}$) is any eigenvalue decomposition of $X$, our proxy only keeps the $k$ largest eigenvalues in the expansion as in

For small values of $k$—we use $k$ between 10 and 20—this can be efficiently computed since we only need to compute the top eigenvectors of a low-rank matrix at each step. Although this approximation gives good empirical results, convergence is no longer guaranteed.

2 Setup

To measure performance, we will use the mean-square error (MSE). However, since a solution $x_{0}$ is only unique up to global phase, it makes no sense to compute the squared distance between $x_{0}$ and the recovery $\hat{x}_{0}$. Rather, we compute the distance to the solution space, i.e. we are interested in the relative MSE defined as

This is the definition we will adopt throughout the paper;Alternatively, we could use $\|x_{0}x_{0}^{*}-\hat{x}_{0}\hat{x}_{0}^{*}\|_{F}/\|x_{0}x_{0}^{*}\|_{F}$, which gives very similar values. the Signal-to-Noise Ratio (SNR) is analogous, namely, $\text{SNR}=\log_{10}(\text{rel.~{}MSE})$.

3 1-D simulations

Phase retrieval for one-dimensional signals arises in fiber optics , speech recognition , but also in the determination of concentration profiles in diffraction imaging . We evaluate PhaseLift for noiseless and noisy data using a variety of different ‘illuminations’ and test signals.

In the first set of experiments we demonstrate the recovery of two very different signals from noiseless data. Both test signals are of length $n=128$. The first signal, shown in Figure 3(a)) is a linear combination of a few sinusoids and represents a typical transfer function one might encounter in optics. The second signal is a complex signal, with independent Gaussian complex entries (each entry is of the form $a+ib$ where $a$ and $b$ are independent $\mathcal{N}(0,1)$ variables) so that the real and imaginary parts are independent white noise sequences; the real part of the signal is shown in Figure 3(b).

Four random binary masks are used to perform the structured illumination so that we measure $|Ax|^{2}$ in which

where each $W_{i}$ is diagonal with either or $1$ on the diagonal, resulting in a total of 512 intensity measurements. We work with the objective functional $\frac{1}{2}\|b-\mathcal{A}(X)\|_{2}^{2}+\lambda\operatorname{Tr}(X)$ and the constraint $X\succeq 0$ to recover the signal, in which we use a small value for $\lambda$ such as $0.05$ since we are dealing with noisefree data. We apply the reweighting scheme as discussed in Section 2.4. (To achieve perfect reconstruction, one would have to let $\lambda\to 0$ as the iteration count increases.) The algorithm is terminated when the relative residual error is less than a fixed tolerance, namely, $\|\mathcal{A}(\hat{x}_{0}\hat{x}_{0}^{*})-b\|_{2}\leq 10^{-6}\|b\|_{2}$, where $\hat{x}_{0}$ is the reconstructed signal just as before. The original and recovered signals are plotted in Figure 3(a) and (b). The SNR is 87.3dB in the first case and 90.5dB in the second.

We have repeated these experiments with the same test signals and the same algorithm, but using Gaussian masks instead of binary masks. In other words, the $W_{i}$’s have Gaussian entries on the diagonal. It turns out that in this case, three illuminations—instead of four—were sufficient to obtain similar performance. This seems to be empirical support for a long-standing conjecture in quantum mechanics due to Ron Wright (see e.g. the concluding section of ). The conjecture is that there exist three unitary operators $U_{1},U_{2},U_{3}$ such that the phase of the (finite-dimensional) signal $x$ is uniquely determined by the measurements $|U_{1}x|,|U_{2}x|,|U_{3}x|$. Our simulations suggest that one can choose $U_{1}=F$, $U_{2}=FW$, and $U_{3}=FW^{\prime}$, where $W,W^{\prime}$ are diagonal matrices with i.i.d. complex normal random variables as diagonal entries. The choice $U_{1}=F$, $U_{2}=I$, $U_{3}=FW$ was equally successful in our experiments, and is a bit closer to the quantum mechanical setting. Furthermore, we point out that no reweighting was needed, when we used six or more Gaussian masks. Expressed differently, plain trace-norm minimization succeeds with $6n$ or more intensity measurements of this kind.

3.2 Noisy measurements

In the next set of experiments, we consider the case when the measurements are contaminated with Poisson noise. The test signal is again a complex random signal as above. Eight illuminations with binary masks are used. We add random Poisson noise to the measurements for five different SNR levels, ranging from about 16dB to about 52dB. Since the solution is known, we have calculated reconstructions for various values of the parameter $\lambda$ balancing the negative log-likelihood and the trace norm, and report results for that $\lambda$ giving the lowest MSE. We implemented this strategy via the standard Golden Section Search . In practice one would have to find the best $\lambda$ via a strategy like cross validation (CV) or generalized cross validation (GCV). For each SNR level we repeated the experiment ten times with different random noise and different binary masks.

Figure 4 shows the average relative MSE in dB (the values of $10\log_{10}(\text{rel.~{}MSE})$ are plotted) versus the SNR. The error curves show clearly the desirable linear behavior between SNR and MSE with respect to the log-log scale. The performance degrades very gracefully with decreasing SNR. Furthermore, the difference of about 5dB between the error curve associated with four measurement and the error curve associated with eight measurements corresponds to an almost twofold error reduction, which is about as much improvement as one can hope to gain by doubling the number of measurements.

3.3 Multiple measurements via oversampling

It is well-known that for one-dimensional signals, oversampling does not result in unique solutions to the phase problem . This might mainly be a theoretical issue without real practical consequences. Our numerical simulations confirm that this is not the case and demonstrate very clearly that oversampling is not useful for most one-dimensional signals. We apply one of Fienup’s algorithms as discussed in Section 4.A in , and PhaseLift with reweighting to real-valued non-negative random signals of length 128. We use oversampling rates ranging from 2 to 5, and stop the algorithms when the relative residual error is less than $10^{-3}$.

The results are displayed in Table 1. As can be seen, both algorithms find nearly perfect fits to the data and yet, they are very far from the original solution. Hence, we have a problem with multiple solutions. The average SNR of the “reconstructions” obtained via Fienup’s method is around 4dB, which is just barely better than if we had used a random guess as a solution. The average SNR for the solutions computed by PhaseLift is only marginally better.

4 2-D simulations

We consider a stylized version of a setup one encounters in X-ray crystallography or diffraction imaging. The test image, shown in Figure 5(a) (magnitude), is a complex-valued image of size $256\times 256$, whose pixel values correspond to the complex transmission coefficients of a collection of gold balls embedded in a medium.

In the first experiment, we demonstrate the recovery of the image shown in Figure 5(a) from noiseless measurements. We consider two different types of illuminations. The first uses Gaussian random masks in which the coefficients on the diagonal of $W_{k}$ are independent real-valued standard normal variables. We use three illuminations, one being constant, i.e. $W_{1}=I$, and the other two Gaussian. Again, we choose a small value of $\lambda$ set to $0.05$ in $\frac{1}{2}\|b-\mathcal{A}(X)\|_{2}^{2}+\lambda\operatorname{Tr}(X)$ since we have no noise, and stop the reweighting iterations as soon as the residual error obeys $\|\mathcal{A}(\hat{x}_{0}\hat{x}_{0}^{*})-b\|_{2}\leq 10^{-3}\|b\|_{2}$. The reconstruction, shown in Figure 5(b), is visually indistinguishable from the original. Since the original image and the reconstruction are complex-valued, we only display the absolute value of each image throughout this and the next subsection.

Gaussian random masks may not be realizable in practice. Our second example uses simple random binary masks, where the entries are either or $1$ with equal probability. In this case, a larger number of illuminations as well as a larger number of reweighting steps are required to achieve a reconstruction of comparable quality. The result for eight binary illuminations is shown in Figure 5(c).

4.2 Noisy measurements

In the second set of experiments we consider the same test image as before, but now with noisy measurements. In the first experiment the SNR is 20dB, in the second experiment the SNR is 60dB. We use 32 Gaussian random masks in each case. The resulting reconstructions are depicted in Figure 6(a) (20dB case) and Figure 6(b) (60dB case). The SNR in the 20dB case is 11.83dB. While the reconstructed image appears slightly more “fuzzy” than the original image, all features of the image are clearly visible. In the 60dB case the SNR is 47.96dB, and the reconstruction is virtually indistinguishable from the original image.

4.3 Multiple measurements via oversampling

Oversampling of two-dimensional signals is widely used to overcome the nonuniqueness of the phase retrieval problem. We now explore whether this approach is viable.

Here, we consider signals with real, non-negative values as test images, a case frequently considered in the literature, see e.g. . These images are of size $128\times 128$. We take noiseless measurements and apply PhaseLift as well as Fienup’s iterative algorithms [3, Section 4.A]. For each method, we terminate the iterations if the relative residual error is less than $10^{-3}$ or if the relative error between two successive iterates is less than $10^{-6}$.

The simulations show that PhaseLift yields reconstructions that fit the measured data well, yielding a small relative residual error $\|\mathcal{A}(X)-y\|_{2}/\|y\|_{2}$, yet the reconstructions are far away from the true signal. This behavior is indicative of an ill-conditioned problem.

The iterates of Fienup’s algorithm stagnated most of the time without converging to a solution. At other times it did yield reconstructions that fit the measured data well, but in either case the reconstruction was always very different from the true signal. Moreover, the reconstructions vary widely depending on the initial (random) guess.

Table 2 displays the results of PhaseLift as well as one of Fienup’s algorithms, namely, the Error Reduction Algorithm [3, Section 4. A] (the other versions yield comparable results). The setup is this: we oversample the signal in each dimension by a factor of $r$, where $r=2,3,4,5$. For each oversampling rate, we run ten experiments using a different test signal each time. The table shows the average residual errors over ten runs as well as the average relative MSE. The ill-posedness of the problem is evident from the disconnect between small residual error and large reconstruction error; that is to say, we fit the data very well and yet observe a large reconstruction error. Thus, in stark contrast to what is widely believed, our simulations indicate that oversampling by itself is not a viable strategy for phase retrieval even for non-negative, real-valued images.

Discussion

This paper introduces a novel framework for phase retrieval, combining multiple illuminations with tools from convex optimization, which been shown to work very well in practice and bears great potential. Having said this, our work also calls for improved theory, improved algorithms and a physical implementation of these ideas. Regarding this last point, it would be interesting to design physical experiments to test our methodology on real data, and we hope to report on this in a future publication. For now, we would like to bring up important open problems.

At the theoretical level, we need to understand for which families of physically implementable structured illuminations does the trace-norm heuristic succeed? How many diffraction patterns are provably sufficient for our convex programming approach to work? Also, we have shown that our approach is robust to noise in the sense that the performance degrades very gracefully as the SNR decreases. Can this be made rigorous? Can we prove that our proposed framework is indeed robust to noise? Here, it is very likely that the tools and ideas developed in the theories of compressed sensing and of matrix completion will play a key role in addressing these fundamental issues.

At the algorithmic level, we need to address the fact that the lifting creates optimization problems of potentially enormous sizes. A tantalizing prospect is whether or not it is possible to use knowledge that the solution has low-rank, e. g. rank one, to design algorithms which do not need to assemble or store very large matrices. If so, how can this be done? Here, randomized algorithms holding up a sketch of the full matrix may prove very helpful.

Finally, we would like to close by returning to another finding of this paper. Namely, oversampling the Fourier transform—this is the same as assuming finite support of the specimen—appears extremely problematic in practice, even for real-valued nonnegative signals. To be sure, we have demonstrated that there typically exist very distinct 2D signals whose modulus of the Fourier transform nearly coincide, whatever the degree of oversampling. In light of this extreme ill posedness, we have trouble understanding why this technique is used so heavily when it does not produce useful results in the absence of very specific a priori information about the image. Moreover, our concern is compounded by the additional fact that algorithms in common use tend to return solutions that depend on an initial guess so that different runs return widely different solutions. In the spirit of reproducible research, this calls for documented results making publicly available both data sets and software so that researchers can reproduce published results or results yet to be published. Of course, one might also be willing to rely on image priors far stronger than finite spatial extent, real-valuedness and positivity; they would, however, need to be specified.

E. C. is partially supported by NSF via grant CCF-0963835 and the 2006 Waterman Award. Y. E. is partially supported by the Israel Science Foundation under Grant no. 170/10. T. S. acknowledges partial support from the NSF via grants DMS 0811169 and DMS-1042939, and from DARPA via grant N66001-11-1-4090. V. V. is supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program. We are indebted to Stefano Marchesini for inspiring and helpful discussions on the phase problem in X-ray crystallography as well as for providing us with the gold balls data set used in Section 4. We want to thank Philippe Jaming for bringing Wright’s conjecture in to our attention.

References