GESPAR: Efficient Phase Retrieval of Sparse Signals

I Introduction

Recovery of a signal from the magnitude of its Fourier transform, also known as phase retrieval, is of great interest in applications such as optical imaging , crystallography , and more . Due to the loss of Fourier phase information, the problem (in 1D) is generally ill-posed. A common approach to overcome this ill-posedeness is to exploit prior information on the signal. A variety of methods have been developed that use such prior information, which may be the signal’s support (region in which the signal is nonzero), non-negativity, or the signal’s magnitude , .

A popular class of algorithms is based on the use of alternate projections between the different constraints. In order to increase the probability of correct recovery, these methods require the prior information to be very precise, for example, exact/or “almost” exact knowledge of the support set. Since the projections are generally onto non-convex sets, convergence to a correct recovery is not guaranteed . A more recent approach is to use matrix-lifting of the problem which allows to recast phase retrieval as a semi-definite programming (SDP) problem . The algorithm developed in does not require prior information about the signal but instead uses multiple signal measurements (e.g., using different illumination settings, in an optical setup).

In order to obtain more robust recovery without requiring multiple measurements, we develop a method that exploits signal sparsity. Existing approaches aimed at recovering sparse signals from their Fourier magnitude belong to two main categories: SDP-based techniques ,,, and algorithms that use alternate projections (Fienup-type methods) . Phase retrieval of sparse signals can be viewed as a special case of the more general quadratic compressed sensing (QCS) problem considered in . Specifically, QCS treats recovery of sparse vectors from quadratic measurements of the form $y_{i}=\mathbf{x}^{T}\mathbf{A}_{i}\mathbf{x},\,\,\,i=1,\ldots,N$, where $\mathbf{x}$ is the unknown sparse vector to be recovered, $y_{i}$ are the measurements, and $\mathbf{A}_{i}$ are known matrices. In (discrete) phase retrieval, $\mathbf{A}_{i}={\bf F}_{i}^{*}{\bf F}_{i}$ where ${\bf F}_{i}$ is the $i$th row of the discrete Fourier transform (DFT) matrix. QCS is encountered, for example, when imaging a sparse object using partially spatially-incoherent illumination .

A general approach to QCS was developed in based on matrix lifting. More specifically, the quadratic constraints where lifted to a higher dimension by defining a matrix variable ${\bf X}={\bf x}{\bf x}^{T}$. The problem was then recast as an SDP involving minimization of the rank of the lifted matrix subject to the recovery constraints as well as row sparsity constraints on ${\bf X}$. An iterative thresholding algorithm based on a sequence of SDPs was then proposed to recover a sparse solution. Similar SDP-type ideas were recently used in the context of phase retrieval ,. However, due to the increase in dimension created by the matrix lifting procedure, the SDP approach is not suitable for large-scale problems.

Another approach for phase retrieval of sparse signals is adding a sparsity constraint to the well-known iterative error reduction algorithm of Fienup . In general, Fienup-type approaches are known to suffer from convergence issues and often do not lead to correct recovery especially in 1D problems; simulation results show that even with the additional information that the input is sparse, convergence is still problematic and the algorithm often recovers erroneous solutions.

In this paper we propose an efficient method for phase retrieval which also leads to good recovery performance. Our approach is based on a fast 2-opt local search method (see for an excellent introduction to such techniques) applied to a sparsity constrained non-linear optimization formulation of the problem. We refer to the resulting algorithm as GESPAR: GrEedy Sparse PhAse Retrieval. Sparsity constrained nonlinear optimization problems have been considered recently in ; the method derived in this paper is motivated – although different in many aspects – by the local search-type techniques of . In essence, GESPAR is a local-search method, where the support of the sought signal is updated iteratively, according to selection rules described in detail in Section III. A local minimum of the objective function is then found given the current support using the damped Gauss Newton algorithm. Theorem 1 establishes convergence of the iterations to a stationary point of the objective under suitable conditions.

We demonstrate through numerical simulations that GESPAR is both efficient and more accurate than current techniques. Several other aspects of the algorithm are explored via simulations such as robustness to noise, and scalability for larger dimensions. In the simulations performed we found that the number of measurements needed for reliable recovery from Fourier magnitudes seems to scale like $s^{3}$, where $s$ is the sparsity level.

GESPAR is applicable to recovery of a sparse vector from general quadratic measurements, and is not restricted to Fourier magnitude measurements. Nonetheless, when the measurements are obtained in the Fourier domain, the algorithm can be implemented efficiently by exploiting the fast Fourier transform, as we discuss in Section IV.

The remainder of the paper is organized as follows. We formulate the problem in Section II. Section III describes our proposed algorithm in detail and establishes convergence of the local iterations. Implementation details for Fourier-based problems are provided in Section IV. Extensive numerical experiments illustrating the empirical performance of GESPAR are presented in Section V.

II Problem Formulation

The mathematical formulation of the problem that we consider consists of minimizing the sum of squared errors subject to the sparsity constraint:

where ${\bf F}_{i}$ is the $i$th row of the DFT matrix ${\bf F}$, and $\|\cdot\|_{0}$ stands for the zero-“norm”, that is, the number of nonzero elements. Note that the unknown vector ${\bf x}$ can only be found up to trivial degeneracies that are the result of the loss of Fourier phase information: circular shift, global phase, and signal “mirroring”.

To aid in solving the phase retrieval problem, we can rely on the fact that the autocorrelation sequence of $\bar{{\bf x}}$ (the first $n$ components of ${\bf x}$) may be determined from ${\bf y}$ if $N\geq 2n-1$. Specifically, let

denote the correlation sequence of length $2n-1$. If we choose $N\geq 2n-1$, then $\{g_{m}\}$ can be obtained by taking the inverse DFT of ${\bf y}$.

Determining $g_{m}$ requires oversampling, or zero-padding of $\mathbf{x}$. While this additional information improves the recovery performance, as is demonstrated in the simulations section, it is not actually needed for GESPAR to work. Nevertheless, when this information is available, GESPAR exploits it, in the following way. First of all, we assume that no support cancelations occur in $\{g_{m}\}$, namely, if $x_{i}\neq 0$ and $x_{j}\neq 0$ for some $i,j$, then $g_{|i-j|}\neq 0$. When the values of ${\bf x}$ are random, this is true with probability 1. This fact can be used in GESPAR in order to obtain initial information on the support of ${\bf x}$, which we capture by two sets $J_{1}$ and $J_{2}$.

Denote by $J_{1}$ the set of indices known in advance to be in the support. To derive the set $J_{1}$, note that due to the existing degree of freedom relating to shift-invariance of ${\bf x}$, the index 1 can be assumed to be in the support, thereby removing this degree of freedom; as a consequence, the index corresponding to the last nonzero element in the autocorrelation sequence is also in the support, i.e.

which will be the formulation to be studied.

Note that even with knowledge of the exact support of $\mathbf{x}$ there is no guarantee for uniqueness beyond the aforementioned trivial degeneracies. Consider for example the two vectors $\mathbf{u}=(1,0,-2,0,-2)$ and $\mathbf{v=}(1-\sqrt{3},0,1,0,1+\sqrt{3})$. Both of these vectors are $s=3$ sparse, and they have the same autocorrelation function $g_{m}=(-2,0,2,0,9,0,2,0,-2)$. This ambiguity therefore cannot be resolved using any method that uses sparsity (even exact support information) and autocorrelation (or Fourier magnitude) measurements alone.

Finally, when the measurements are noisy, the autocorrelation information is not very useful for support estimation, since very small (noise level) values in the autocorrelation sequence cannot be treated as zero. For this reason, the autocorrelation-derived support information is not used in GESPAR at all in the noisy case. Formally, ignoring this information is equivalent to setting $J_{1}=\{1\}$ and $J_{2}=\{1,2,\ldots,n\}$.

II-B Sparse Phase Retrieval: General Measurements

Although the problem formulation above assumes Fourier measurements and sparsity of $\mathbf{\bar{x}}$, we show below that our approach applies to arbitrary quadratic measurements of $\mathbf{\bar{x}}$. This includes the case in which $\bar{\mathbf{x}}$ is sparse in a basis other than the identity basis. In fact, in this general case, the formulation given in (4) remains the same, with the only change being the definition of the matrices $\mathbf{A}_{i}$.

Consider the phase retrieval problem with respect to arbitrary linear measurements, so that

In the next section, we propose GESPAR—an iterative local-search based algorithm for solving (4). We note that although in the context of phase retrieval the parameters ${\bf A}_{i},J_{1},J_{2}$ have special properties (e.g., ${\bf A}_{i}$ is positive semidefinite of at most rank 2, $|J_{1}|=2$), we will not use these properties in GESPAR. Therefore, our approach is capable of handling general instances of (4) with the sole assumption that ${\bf A}_{i}$ is symmetric for any $i=1,2,\ldots,N$. In the Fourier case, the algorithm can be implemented more efficiently, as we discuss in Section IV.

III GrEedy Sparse PhAse Retrieval (GESPAR)

Before describing our algorithm, we begin by presenting a variant of the damped Gauss-Newton (DGN) method , that is in fact the core step of our approach. The DGN method is invoked in order to solve the problem of minimizing the objective function $f$ over a given support $S\subseteq\{1,2,\ldots,n\}\;(|S|=s)$:

The minimization in (7) is a nonlinear least-squares problem. A natural approach for tackling it is via the DGN method. This algorithm begins with an arbitrary vector ${\bf z}_{0}$. In our simulations, we choose it as a white random Gaussian vector with zero mean and unit variance. At each iteration, all the terms inside the squares in $g({\bf z})$ are linearized around the previous guess. Namely, we write $g({\bf z})$ from (7) as:

with $h_{i}({\bf z})={\bf z}^{T}{\bf B}_{i}{\bf z}-y_{i}$, and ${\bf B}_{i}={\bf U}_{S}^{T}{\bf A}_{i}{\bf U}_{S}$. At each step we replace $h_{i}$ by its linear approximation around ${\bf z}_{k-1}$:

We then choose ${\bf z}_{k}$ to be the solution of the problem

Problem (10) can be written as a linear least-squares problem

The following theorem establishes the rate of convergence of the norm of the gradient of the objective function to zero, and consequently proves that the limit points of the sequence are stationary points.

Let $\{{\bf z}_{k}\}$ be the sequence generated by the DGN method. Assume that $\sum_{i=1}^{N}{\bf B}_{i}\succ{\bf 0}$ and that there exists $\underline{\lambda}>0$ such that for all $k$

Then $\nabla g({\bf z}_{k})\rightarrow{\bf 0}$ as $k\rightarrow\infty$ and there exists a constant $C>0$ such that

Moreover, each limit point of the sequence is a stationary point of $g$.

Note that the proof requires $J({\bf z}_{k})$ to have full column rank, and in fact that the minimum eigenvalues of $J({\bf z}_{k})^{T}J({\bf z}_{k})$ are uniformly bounded below. In the vast majority of our runs this assumption held true; however, we did encounter in our numerical experiments a few cases in which this condition was not valid. In these situations, our implementation chose one of the optimal solutions of the corresponding least-squares problem. We noticed that these cases had negligible effect on the results.

III-B The 2-opt Local Search Method

The GESPAR method consists of repeatedly invoking a local-search method on an initial random support set. In this section we describe the local search procedure. At the beginning, the support is chosen to be a set of $s$ random indices chosen to satisfy the support constraints $J_{1}\subseteq S\subseteq J_{2}$. Then, at each iteration a swap between a support and an off-support index is performed such that the resulting solution via the DGN method improves the objective function. Since at each iteration only two elements are changed (one in the support and one in the off-support), this is a so-called “2-opt” method (see ). The swaps are always chosen to be between the index corresponding to components in the current iterate ${\bf x}_{k-1}$ with the smallest absolute value and the off-support index corresponding to the component of $\nabla f({\bf x}_{k-1})=4\sum_{i}(\mathbf{x}_{k-1}^{T}\mathbf{A}_{i}\mathbf{x}_{k-1}-\mathbf{c}_{i})\mathbf{A}_{i}\mathbf{x}_{k-1}$ with the largest absolute value. This process continues as long as the objective function decreases and stops when no improvement can be made. A detailed description of the method is given in Algorithm 2.

III-C The GESPAR Algorithm

The 2-opt method can have the tendency to get stuck at local optima points. Therefore, our final algorithm, which we call GESPAR, is a restarted version of 2-opt. The 2-opt method is repeatedly invoked with different initial random support sets until the resulting objective function value is smaller than a certain threshold (success) or the number of maximum allowed total number of swaps was passed (failure). A detailed description of the method is given in Algorithm 3. One element of our specific implementation that is not described in Algorithm 3 is the incorporation of random weights added to the objective function, giving randomly different weights to the different measurements. Namely, the objective function used is actually chosen as $f({\bf x})=\sum_{i=1}^{N}w_{i}({\bf x}^{T}{\bf A}_{i}{\bf x}-y_{i})^{2}$ with $w_{i}=1$ or $2$ with equal probability. The random generation of weights is done each time the DGN procedure is invoked. We observed that this modification reduced the probability of the 2-opt procedure to get stuck in non-optimal points.

IV Fourier Implementation Details

In principle, GESPAR may be used to find sparse solutions to any system of quadratic equations, i.e. problems of the form:

However, when the matrices ${\bf A}_{i}$ correspond to transforms that can be implemented efficiently, GESPAR takes on a particularly simple form.

For example, consider the case in which $\{{\bf A}_{i}\}$ represent Fourier measurements. In this case, the creation and storing of the matrices ${\bf A}_{i}$ defined in Section II, can be avoided in the implementation, by using the FFT. Specifically, to calculate the weighted objective function, we note that

where $\hat{x}_{i}$ is the $i$th DFT component of ${\bf x}$, which can be computed via the FFT. Clearly, $J({\bf z})$, which is used in the DGN procedure (Algorithm 1) can also be computed efficiently since ${\bf B}_{i}={\bf U}_{S}^{T}{\bf A}_{i}{\bf U}_{S}$ only involves a small ($s$) number of columns of the Fourier matrix $\mathbf{F}$.

The FFT can also be used in the calculation of the gradient $\nabla f({\bf x})$, used in the 2-opt stage 2:

Consequently, in no step of the algorithm is it necessary to calculate the set of matrices $\mathbf{A}_{i}$ explicitly.

This fact is even more important in the 2D Fourier phase retrieval problem, as the relevant vector sizes become very large. Since a major advantage of GESPAR over other methods (e.g. SDP based) is its low computational cost, GESPAR may be used to find a sparse solution to the 2D Fourier phase retrieval - or phase retrieval of images. The only adjustments needed in the algorithm are in the implementation, for example, using FFT2 instead of storing the large matrices ${\bf A}_{i}$.

Figure 1 shows a recovery example of a sparse $195\times 195$ pixel image, comprised of $s=15$ circles at random locations and random values on a grid containing $225$ points, recovered from its $38,025$ 2D-Fourier magnitude measurements, using GESPAR. The dictionary used in this example contains 225 elements consisting of non-overlapping circles located on a $15\times 15$ point cartesian grid, each with a 13 pixel diameter. The solution took 80 seconds. Solving the same problem using the sparse Fienup algorithm did not yield a successful reconstruction, and using the SDP method is not practical due to the large matrix sizes.

Further investigation of the algorithm’s performance in the 2D case is presented in Section V.

V Numerical Simulations

In order to demonstrate the performance of GESPAR, we conduct several numerical simulations. The algorithm is compared to other existing methods, and is evaluated in terms of signal-recovery accuracy, computational efficiency, and robustness to noise.

In this subsection we examine the recovery success rate of GESPAR as a function of the number of non-zero elements in the signal. A runtime comparison of the tested methods is also performed.

We choose $\bar{{\bf x}}$ as a random vector of length $n$. The vector contains uniformly distributed values in the range $\cup$ in $s$ randomly chosen elements. The $N$ point DFT of the signal is calculated, and its magnitude-square is taken as ${\bf y}$, the vector of measurements. The $2n-1$ point correlation is also calculated. In order to recover the unknown vector ${\bf x}$, the GESPAR algorithm is used with $\tau=10^{-4}$ and $ITER=6400$. We also test two other algorithms for comparison purposes: An SDP based algorithm (Algorithm 2, ), and an iterative Fienup algorithm with a sparsity constraint . In our simulation $n=64$ and $N=128$. The Sparse-Fienup algorithm is run using $100$ random initial points, out of which the chosen solution is the one that best matches the measurements. Namely, $\hat{\mathbf{x}}$ is selected as the $s$ sparse output of the Sparse-Fienup algorithm with the minimal cost $f({\bf x})=\sum_{i=1}^{N}(|{\bf F}_{i}{\bf x}|^{2}-y_{i})^{2}$ out of the $100$ runs.

Signal recovery results of the numerical simulation are shown in Fig. 2, where the probability of successful recovery is plotted for different sparsity levels. The success probability is defined as the ratio of correctly recovered signals ${\bf x}$ out of $100$ simulations. In each simulation both the support and the signal values are randomly selected. The three algorithms (GESPAR, SDP and Sparse-Fienup) are compared. The results clearly show that GESPAR outperforms the other methods in terms of probability of successful recovery - over 90% successful recovery up to $s=15$, vs. $s=8$ and $s=7$ in the other two techniques.

Average runtime comparison of the three algorithms is shown in Table I for $n=64$ and $N=128$. The runtime is averaged over all successful recoveries. The computer used has an intel i5 CPU and 4GB of RAM. As seen in the table, the SDP based algorithm is significantly slower than the other two methods. Fienup iterations are slightly slower than GESPAR and lead to a much lower success rate. In these simulations, GESPAR is both fast and more accurate than its competitors.

V-B Sensitivity to exact sparsity knowledge

Since the exact value of the signal’s sparsity $s$ may not be known, the performance of GESPAR is examined when only an upper limit on $s$ is given. To this end we run GESPAR twice: Once with $s$ known exactly at each realization, and once with only an upper limit on $s$. The upper limit is taken as $25$. The other simulation settings are the same as in SectionV-A.

Figure 3 shows the probability for successful recovery of the two simulations. The rather loose upper limit on $s$ does not seem to affect the results significantly— in fact, the performance is somewhat improved when allowing more nonzero elements during the iterations.

V-C Effect of the number of allowed swaps

One of the stopping criteria for the GESPAR algorithm is when the total number of swaps exceeds a predefined parameter (the input parameter $ITER$ in Algorithm 3). Naturally, increasing the allowed number of index swaps will increase the probability of finding a correct solution, but at the cost of increased computation time. It is therefore important to quantify this effect, which is the purpose of the current simulation.

We run GESPAR with the same parameters as in SectionV-A several times, where in each simulation a different value for the parameter $ITER$ is used, in the range $$. Figure 4 shows the results. As expected, increasing the number of possible swaps increases the recovery probability. Note that increasing the value of$ITER$above$6400$demonstrated no improvement in the recovery results - for the unsuccessful recoveries, increasing the number of swaps even up to$ITER=25600$did not help. This means that for these simulation values (e.g.$N=128,\,s<25$), using a value of$ITER$larger than$6400$ only increases computation time without improving the results.

V-D Effect of oversampling and support information

Here we examine the effect of oversampling and of autocorrelation-derived support information. GESPAR is run on random vectors $\mathbf{x}$ of length $n=64$, with a varying amount of noiseless Fourier magnitude measurements, obtained by the $N$ point DFT of $\mathbf{x}$ with $N=64,128,256$. In these cases, no support information was used - i.e. $J_{1}=\{1\}$ and $J_{2}=\{1,2,\ldots,n\}$. In addition, in order to investigate the effect of support information, we run GESPAR with $n=64,N=128$ (i.e. oversampling by a factor of 2), and use the support information derived from the autocorrelation sequence. The results, shown in Fig. 5, clearly show that both oversampling and support information improve performance.

V-E Robustness to noise

We now evaluate GESPAR as a function of SNR, and compare it with sparse Fienup . The SDP based method presented in is not designed to deal with noise and therefore we did not apply it here. The SDP approach of considers random measurements, and does not produce comparable results from direct Fourier measurements.

As in Section V-A, we choose $\mathbf{\bar{x}}$ as a vector of length $n$, with $s$ randomly chosen elements containing uniformly distributed values, and evaluate its $N$ point Fourier magnitude-square. White-gaussian noise $\mathbf{v}$ is added to the measurements, at different SNR values, defined as: $SNR=20\text{log}\frac{\|\mathbf{y\|}}{\|\mathbf{v\|}}$. In order to recover the unknown vector ${\bf x}$, the GESPAR algorithm is used with $\tau=10^{-4}$ and $ITER=10000$, as well as the sparse-Fienup algorithm, for comparison purposes. In our simulation $n=64$ and $N=128$. The sparse-Fienup algorithm is run with a maximum of $1000$ iterations, and with $100$ random initial points.

Note that even with little noise, the information on the support obtained by the zeros of the autocorrelation is no longer available. This is due to the fact that in the presence of noise, there will be no true zeros in the measured (or calculated) autocorrelation. In this case, one might try to threshold the autocorrelation values, rendering small autocorrelation values as zeros. However, this might result in zeroing of small (yet non-zero) values of the true autocorrelation function. Therefore, in the noisy case, we do not use support information obtained by the autocorrelation function in GESPAR, namely $J_{2}=\{1,2,\ldots,n\}$.

Figure 6 shows the normalized mean squared reconstruction error (NMSE), defined as $NMSE=\frac{\|\mathbf{x}-\mathbf{\hat{x}}\|_{2}}{\|\mathbf{x}\|_{2}}$, as a function of sparsity, for different SNR values. Each point represents an average over 100 different random realizations. The performance under different SNR values is plotted for GESPAR (full lines), and for sparse-Fienup (dashed-lines). The performance of GESPAR naturally improves as SNR increases, and it clearly outperforms sparse-Fienup in terms of noise-robustness.

V-F Scalability

As one of the main advantages of GESPAR over SDP based methods is its ability to solve large problems efficiently, we now examine its performance for different vector sizes.

We simulate GESPAR for various values of $n\in$. In all cases $N=2n$. The other simulation parameters are as in Section V-A. The recovery probability vs. sparsity $s$ for different vector lengths is shown in Fig. 7. The maximal sparsity $s$ allowing successful recovery is shown to increase with vector length $n$, and seems to scale like $n^{1/3}$, which is consistent with the same scaling observation presented in . The mean reconstruction time for a signal with $n=512,\,s=35$ from $N=1024$ measurements, allowing $ITER=6400$ replacements, is $33.5$ seconds. For comparison, a corresponding plot representing the scalability of the sparse-Fienup algorithm is presented in Fig. 8. Plotting a similar scalability plot for the SDP based method is not possible due to the high computational cost which under our simulation conditions limits the application of this method to around $n\sim 400$.

V-G Computation Time

The most time consuming part of GESPAR is the matrix inversion process in the DGN segment of the algorithm. Therefore, computation time scales approximately linearly with the number of swaps - as each swap corresponds to a single DGN iteration. The approximately linear dependence of runtime in the number of swaps is displayed in Fig. 9. Each point in the plot represents the mean time it took GESPAR to run $ITER$ iterations, averaged over 50 random input signals with $N=128,\,n=64,s=10$.

A major factor that determines the computation time is the number of index swaps required to find a solution. The mean number of swaps as a function of $s,n$ is shown in Fig. 10. Beyond the successful recovery region (the white region in Fig. 7), the maximal number of swaps ($6400$) is used, without yielding a correct solution.

V-H Two-Dimensional Fourier Phase Retrieval

In this section we apply GESPAR to 2D Fourier phase retrieval problems, showing its ability to solve large scale problems efficiently.

We generate random $s-$sparse 2D signals of sizes $\sqrt{n}\times\sqrt{n}$, with varying values for $s$ and $n$, in the ranges $s\in[2:82]$ and $n\in[256:6400]$. Each signal is recovered from the noiseless magnitude of its 2D DFT, with no oversampling, using GESPAR. Similarly to the 1D noisy simulation, no autocorrelation-derived support information was used here. The parameter $ITER$ is taken as $6400$. The recovery probability vs. sparsity $s$ for different vector lengths is shown in Fig. 11. Similarly to the 1D case, the maximal sparsity allowing successful recovery increases with $n$. For comparison, Fig. 12 shows the result of a sparse-Fienup scalability simulation for the 2D case, under the same conditions, with 200 initial points per signal (increasing this parameter did not affect the results significantly). GESPAR is shown to outperform the sparse-Fienup method in the 2D case as well. As in the 1D case, a comparison to SDP based methods is not included here, since applying the SDP based method on the 2D case is very difficult due to memory limitations.

A comparison between GESPAR and the sparse-Fienup method is shown in Fig. 13. The comparison shows the average time a successful recovery in the simulation took, as a function of vector size $n$. Sparse-Fienup is seen to be faster, however comparing Fig. 11 to Fig. 12 shows that GESPAR can recover signals up with a higher value of $s$: For example, when $n=6400$, GESPAR recovers with very high probability signals up to sparsity $s=57$, while sparse Fienup only recovers up to $s=42$.

VI Conclusion

We proposed and demonstrated GESPAR - a fast algorithm for recovering a sparse vector from its Fourier magnitude, or more generally, from quadratic measurements. We showed via simulations that GESPAR outperforms alternative approaches suggested for this problem in terms of complexity and success probability. The algorithm does not require matrix-lifting, and therefore is potentially suitable for large scale problems such as 2D images. The simulations demonstrated robustness of GESPAR to noise and other inexact knowledge, as well as its ability to successfully treat a variety of phase retrieval problems in one and two dimensions.

References

Appendix A Proof of Theorem

Define the vector-valued function ${\bf h}$ by

with $h_{i}({\bf z})={\bf z}^{T}{\bf B}_{i}{\bf z}-y_{i}$ With this notation, the vector ${\bf b}_{k}$ can be written as

and the solution of the least-squares problem is

From (16) it follows that $-{\bf d}_{k}$ is a descent direction since

We now show that the sequence generated by the DGN method is bounded. Indeed, since $-{\bf d}_{k}$ is a descent direction,

where the second inequality is due to Cauchy-Schwarz and the last inequality is a result of the fact that $\sum_{i=1}^{N}{\bf B}_{i}\succ 0$ and $y_{i}\geq 0$. Therefore,

proving that $\{{\bf z}_{k}\}\subseteq B[{\bf 0},\sqrt{\alpha}]=\{{\bf z}:\|{\bf z}\|\leq\sqrt{\alpha}\}.$

Since $g$ is twice continuously differentiable, and $J({\bf z})$ is continuous, it follows that there exists $M>0$ and $\Lambda>0$ such that $\lambda_{\max}(\nabla^{2}g({\bf z}))\leq M$ and $\lambda_{\max}(J({\bf z})^{T}J({\bf z}))\leq\Lambda$ for any ${\bf z}\in B[{\bf 0},2\sqrt{\alpha}]$. In addition, since $\nabla g$ is continuous over $B[{\bf 0},2\sqrt{\alpha}]$, there exist $\beta>0$ such that $\|\nabla g({\bf z})\|\leq\beta$ for all ${\bf z}\in B[0,2\sqrt{\alpha}]$. Therefore, by (17) it follows that

The fact that $\lambda_{\max}(\nabla^{2}g({\bf z}))\leq M$ for all ${\bf z}\in B[{\bf 0},2\sqrt{\alpha}]$ implies that $\nabla g$ is Lipschitz continuous over $B[{\bf 0},2\sqrt{\alpha}]$ with parameter $M>0$. Hence, by the descent lemma ,

for any ${\bf x},{\bf y}\in B[{\bf 0},2\sqrt{\alpha}]$.

From $\|{\bf z}_{k-1}\|\leq\sqrt{\alpha}$ and $\|{\bf d}_{k}\|\leq\beta/(2\underline{\lambda})$, it follows that ${\bf z}_{k-1}-t{\bf d}_{k}\in B[{\bf 0},2\sqrt{\alpha}]$ whenever $t\leq\frac{2\underline{\lambda}\sqrt{\alpha}}{\beta}$. Therefore, we can plug ${\bf y}={\bf z}_{k-1}-t{\bf d}_{k}$ and ${\bf x}={\bf z}_{k-1}$ into (20) to obtain

Therefore, if $t\leq\min\left\{\frac{2\underline{\lambda}}{M},\frac{2\underline{\lambda}\sqrt{\alpha}}{\beta}\right\}$, then

By the way the backtracking procedure is defined, we have that either $t_{k}=1$ or $2t_{k}>\min\left\{\frac{2\underline{\lambda}}{M},\frac{2\underline{\lambda}\sqrt{\alpha}}{\beta}\right\}$ and hence $t_{k}\geq\min\left\{1,\frac{\underline{\lambda}}{M},\frac{\underline{\lambda}\sqrt{\alpha}}{\beta}\right\}$. Together with (21) this results in the inequality

where $C=\min\left\{\frac{1}{2\Lambda},\frac{\underline{\lambda}}{2M\Lambda},\frac{\underline{\lambda}\sqrt{\alpha}}{2\beta\Lambda}\right\}$. Noting that $\{g({\bf z}_{k})\}$ is a bounded below and nonincreasing sequence, it follows that it converges. The left-hand side of (22) therefore converges to zero and we obtain the result that $\nabla g({\bf z}_{k})$ converges to zero as $k$ tends to infinity. This fact also readily implies that all accumulation points of the sequence are stationary. Summing the inequality (22) over $p=1,2,\ldots,k+1$ we obtain that

and consequently, (also using the fact that $g({\bf z}_{k+1})\geq 0$),

from which the inequality (12) follows. $\Box$