Phase Retrieval via Wirtinger Flow: Theory and Algorithms

Introduction

We are interested in solving quadratic equations of the form

On the theoretical side, several combinatorial optimization problems—optimization programs with discrete design variables which take on integer or Boolean values—can be cast as solving quadratic equations or as minimizing a linear objective subject to quadratic inequalities. In their most general form these problems are known to be notoriously difficult (NP-hard) [52, Section 4.3]. Nevertheless, many heuristics have been developed for addressing such problems.For a partial review of some of these heuristics as well as some recent theoretical advances in related problems we refer to our companion paper [16, Section 1.6] and references therein . One popular heuristic is based on a class of convex relaxations known as Shor’s relaxations [52, Section 4.3.1] which can be solved using tractable semi-definite programming (SDP). For certain random models, some recent SDP relaxations such as PhaseLift are known to provide exact solutions (up to global phase) to the generalized phase retrieval problem using a near minimal number of sampling vectors . While in principle SDP based relaxations offer tractable solutions, they become computationally prohibitive as the dimension of the signal increases. Indeed, for a large number of unknowns in the tens of thousands, say, the memory requirements are far out of reach of desktop computers so that these SDP relaxations are de facto impractical.

Algorithm: Wirtinger Flow

This paper introduces an approach to phase retrieval based on non-convex optimization as well as a solution algorithm, which has two components: (1) a careful initialization obtained by means of a spectral method, and (2) a series of updates refining this initial estimate by iteratively applying a novel update rule, much like in a gradient descent scheme. We refer to the combination of these two steps, introduced in reverse order below, as the Wirtinger flow (WF) algorithm.

Our approach to (2.1) is simply stated: start with an initialization $\bm{z}_{0}$, and for $\tau=0,1,2,\ldots$, inductively define

2 Initialization via a spectral method

Our main result states that for a certain random model, if the initialization $\bm{z}_{0}$ is sufficiently accurate, then the sequence $\{\bm{z}_{\tau}\}$ will converge toward a solution to the generalized phase problem (1.1). In this paper, we propose computing the initial guess $\bm{z}_{0}$ via a spectral method, detailed in Algorithm 1. In words, $\bm{z}_{0}$ is the leading eigenvector of the positive semidefinite Hermitian matrix $\sum_{r}y_{r}\bm{a}_{r}\bm{a}_{r}^{*}$ constructed from the knowledge of the sampling vectors and observations. (As usual, $\bm{a}_{r}^{*}$ is the adjoint of $\bm{a}_{r}$.) Letting $\bm{A}$ be the $m\times n$ matrix whose $r$th row is $\bm{a}_{r}^{*}$ so that with obvious notation $\bm{y}=|\bm{A}\bm{x}|^{2}$, $\bm{z}_{0}$ is the leading eigenvector of $\bm{A}^{*}\,\operatorname{diag}\{\bm{y}\}\,\bm{A}$ and can be computed via the power method by repeatedly applying $\bm{A}$, entrywise multiplication by $\bm{y}$ and $\bm{A}^{*}$. In the theoretical framework we study below, a constant number of power iterations would give machine accuracy because of an eigenvalue gap between the top two eigenvalues, please see Appendix B for additional information.

3 Wirtinger flow as a stochastic gradient scheme

We would like to motivate the Wirtinger flow algorithm and provide some insight as to why we expect it to work in a model where the sampling vectors are random. First, we emphasize that our statements in this section are heuristic in nature; as it will become clear in the proof Section 7, a correct mathematical formalization of these ideas is far more complicated than our heuristic development here may suggest. Second, although our ideas are broadly applicable, it makes sense to begin understanding the algorithm in a setting where everything is real valued, and in which the vectors $a_{r}$ are i.i.d. $\mathcal{N}(\bm{0},\bm{I})$. Also without any loss in generality and to simplify exposition in this section we shall assume $\left\|\bm{x}\right\|=1$.

Let $\bm{x}$ be a solution to (1.1) so that $y_{r}=|\langle\bm{a}_{r},\bm{x}\rangle|^{2}$, and consider the initialization step first. In the Gaussian model, a simple moment calculation gives

We now turn our attention to the gradient-update (2.2) and define

where here and below, $\bm{x}$ is once again our planted solution. The first term ensures that the direction of $\bm{z}$ matches the direction of $\bm{x}$ and the second term penalizes the deviation of the Euclidean norm of $\bm{z}$ from that of $\bm{x}$. Obviously, the minimizers of this function are $\pm\bm{x}$. Now consider the gradient scheme

In Section 7.9, we show that if $\text{min }\left\|\bm{z}_{0}\pm\bm{x}\right\|\leq 1/8\left\|\bm{x}\right\|$, then $\{\bm{z}_{\tau}\}$ converges to $\bm{x}$ up to a global sign. However, this is all ideal as we would need knowledge of $\bm{x}$ itself to compute the gradient of $F$; we simply cannot run this algorithm.

Hence, the average WF update is the same as that in (2.3) so that we can interpret the Wirtinger flow algorithm as a stochastic gradient scheme in which we only get to observe an unbiased estimate $\nabla f(\bm{z})$ of the “true” gradient $\nabla F(\bm{z})$.

Regarding WF as a stochastic gradient scheme helps us in choosing the learning parameter or step size $\mu_{\tau}$. Lemma 7.7 asserts that

holds with high probability. Looking at the right-hand side, this says that the uncertainty about the gradient estimate depends on how far we are from the actual solution $\bm{x}$. The further away, the larger the uncertainty or the noisier the estimate. This suggests that in the early iterations we should use a small learning parameter as the noise is large since we are not yet close to the solution. However, as the iterations count increases and we make progress, the size of the noise also decreases and we can pick larger values for the learning parameter. This heuristic together with experimentation lead us to consider

shown in Figure 1. Values of $\tau_{0}$ around $330$ and of $\mu_{\text{max}}$ around $0.4$ worked well in our simulations. This makes sure that $\mu_{\tau}$ is rather small at the beginning (e.g. $\mu_{1}\approx 0.003$ but quickly increases and reaches a maximum value of about $0.4$ after 200 iterations or so.

Main Results

Our main result establishes the correctness of the Wirtinger flow algorithm in the Gaussian model defined below. Later in Section 5, we shall also develop exact recovery results for a physically inspired diffraction model.

We also need to define the distance to the solution set.

with probability at least $1-10e^{-\gamma n}-{8}/{n^{2}}$ ($\gamma$ is a fixed positive numerical constant). Further, take a constant learning parameter sequence, $\mu_{\tau}=\mu$ for all $\tau=1,2,\ldots$ and assume $\mu\leq c_{1}/n$ for some fixed numerical constant $c_{1}$. Then there is an event of probability at least $1-13e^{-\gamma n}-me^{-1.5m}-{8}/{n^{2}}$, such that on this event, starting from any initial solution $\bm{z}_{0}$ obeying (3.1), we have

Setting $\mu=c_{1}/n$ yields $\epsilon$ accuracy in a relative sense, namely, $\operatorname{dist}(\bm{z},\bm{x})\leq\epsilon\left\|\bm{x}\right\|$, in $\mathcal{O}(n\,\log 1/\epsilon)$ iterations. The computational work at each iteration is dominated by two matrix-vector products of the form $\bm{A}\bm{z}$ and $\bm{A}^{*}\bm{v}$, see Appendix B. It follows that the overall computational complexity of the WF algorithm is $\mathcal{O}(mn^{2}\log{1}/{\epsilon})$. Later in the paper, we will exhibit a modification to the WF algorithm of mere theoretical interest, which also yields exact recovery under the same sampling complexity and an $\mathcal{O}(mn\log{1}/{\epsilon})$ computational complexity; that is to say, the computational workload is now just linear in the problem size.

2 Comparison with other non-convex schemes

We now pause to comment on a few other non-convex schemes in the literature. Other comparisons may be found in our companion paper .

Earlier, we discussed the Gerchberg-Saxton and Fienup algorithms. These formulations assume that $\bm{A}$ is a Fourier transform and can be described as follows: suppose $\bm{z}_{\tau}$ is the current guess, then one computes the image of $\bm{z}_{\tau}$ through $\bm{A}$ and adjust its modulus so that it matches that of the observed data vector: with obvious notation,

where $\odot$ is elementwise multiplication, and $\bm{b}=|\bm{A}\bm{x}|$ so that $b_{r}^{2}=y_{r}$ for all $r=1,\ldots,m$. Then

(In the case of Fourier data, the step (3.2)–(3.3) essentially adjusts the modulus of the Fourier transform of the current guess so that it fits the measured data.) Finally, if we know that the solution belongs to a convex set $\mathcal{C}$ (as in the case where the signal is known to be real-valued, possibly non-negative and of finite support), then the next iterate is

where $P_{\mathcal{C}}$ is the projection onto the convex set $\mathcal{C}$. If no such information is available, then $\bm{z}_{\tau+1}=\bm{v}_{\tau+1}$. The first step (3.3) is not a projection onto a convex set and, therefore, it is in general completely unclear whether the Gerchberg-Saxton algorithm actually converges. (And if it were to converge, at what speed?) It is also unclear how the procedure should be initialized to yield accurate final estimates. This is in contrast to the Wirtinger flow algorithm, which in the Gaussian model is shown to exhibit geometric convergence to the solution to the phase retrieval problem. Another benefit is that the Wirtinger flow algorithm does not require solving a least-squares problem (3.3) at each iteration; each step enjoys a reduced computational complexity.

A recent contribution related to ours is the interesting paper , which proposes an alternating minimization scheme named AltMinPhase for the general phase retrieval problem. AltMinPhase is inspired by the Gerchberg-Saxton update (3.2)–(3.3) as well as other established alternating projection heuristics . We describe the algorithm in the setup of Theorem 3.3 for which gives theoretical guarantees. To begin with, AltMinPhase partitions the sampling vectors $\bm{a}_{r}$ (the rows of the matrix $\bm{A}$) and corresponding observations $y_{r}$ into $B+1$ disjoint blocks $(\bm{y}^{(0)},\bm{A}^{(0)})$, $(\bm{y}^{(1)},\bm{A}^{(1)})$, $\ldots$, $(\bm{y}^{(B)},\bm{A}^{(B)})$ of roughly equal size. Hence, distinct blocks are stochastically independent from each other. The first block $(\bm{y}^{(0)},\bm{A}^{(0)})$ is used to compute an initial estimate $\bm{z}_{0}$. After initialization, AltMinPhase goes through a series of iterations of the form (3.2)–(3.3), however, with the key difference that each iteration uses a fresh set of sampling vectors and observations: in details,

The main point is that in practice, it is not realistic to imagine (1) that we will divide the samples in distinct blocks (how many blocks should we form a priori? of which sizes?) and (2) that we will use measured data only once. With respect to the latter, observe that the Gerchberg-Saxton procedure (3.2)–(3.3) uses all the samples at each iteration. This is the reason why AltMinPhase is of little practical value, and of theoretical interest only. As a matter of fact, its design and study seem merely to stem from analytical considerations: since one uses an independent set of measurements at each iteration, $A^{(\tau+1)}$ and $z_{\tau}$ are stochastically independent, a fact which considerably simplifies the convergence analysis. In stark contrast, the WF iterate uses all the samples at each iteration and thus introduces some dependencies, which makes for some delicate analysis. Overcoming these difficulties is crucial because the community is preoccupied with convergence properties of algorithms one actually runs, like Gerchberg-Saxton (3.2)–(3.3), or would actually want to run. Interestingly, it may be possible to use some of the ideas developed in this paper to develop a rigorous theory of convergence for algorithms in the style of Gerchberg-Saxton and Fienup, please see .

In a recent paper , which appeared on the arXiv preprint server as the final version of this paper was under preparation, the authors explore necessary and sufficient conditions for the global convergence of an alternative minimization scheme with generic sampling vectors. The issue is that we do not know when these conditions hold. Further, even when the algorithm converges, it does not come with an explicit convergence rate so that is is not known whether the algorithm converges in polynomial time. As before, some of our methods as well as those from our companion paper may have some bearing upon the analysis of this algorithm. Similarly, another class of nonconvex algorithms that have recently been proposed in the literature are iterative algorithms based on Generalized Approximate Message Passing (GAMP), see and as well as for some background literature on AMP. In , the authors demonstrate a favorable runtime for an algorithm of this nature. However, this does not come with any theoretical guarantees.

Moving away from the phase retrieval problem, we would like to mention some very interesting work on the matrix completion problem using non-convex schemes by Montanari and coauthors , see also . Although the problems and models are quite different, there are some general similarities between the algorithm named OptSpace in and ours. Indeed, OptSpace operates by computing an initial guess of the solution to a low-rank matrix completion problem by means of a spectral method. It then sets up a nonconvex problem, and proposes an iterative algorithm for solving it. Under suitable assumptions, demonstrates the correctness of this method in the sense that OptSpace will eventually converge to a low-rank solution, although it is not shown to converge in polynomial time.

Numerical Experiments

We present some numerical experiments to assess the empirical performance of the Wirtinger flow algorithm. Here, we mostly consider a model of coded diffraction patterns reviewed below.

We consider an acquisition model, where we collect data of the form

We emphasize that these models are physically realizable in optical applications specially those that arise in microscopy. However, we should note that phase retrieval has many different applications and in some cases other models may be more convenient. We refer to our companion paper Section 2.2 for a discussion of other practically relevant models.

2 The Gaussian and coded diffraction models

Random low-pass signals. Here, $\bm{x}$ is given by

with $M=n/8$ and $X_{k}$ and $Y_{k}$ are i.i.d. $\mathcal{N}(0,1)$.

where $X_{k}$ and $Y_{k}$ are are i.i.d. $\mathcal{N}(0,1/8)$ so that the low-pass model is a ‘bandlimited’ version of this high-pass random model (variances are adjusted so that the expected signal power is the same).

Below, we set $n=128$, and generate one signal of each type which will be used in all the experiments.

The initialization step of the Wirtinger flow algorithm is run by applying $50$ iterations of the power method outlined in Algorithm 3 from Appendix B. In the iteration (2.2), we use the parameter value $\mu_{\tau}=\min(1-\exp(-\tau/\tau_{0}),0.2)$ where $\tau_{0}\approx 330$. We stop after $2,500$ iterations, and report the empirical probability of success for the two different signal models. The empirical probability of succcess is an average over $100$ trials, where in each instance, we generate new random sampling vectors according to the Gaussian or CDP models. We declare a trial successful if the relative error of the reconstruction $\operatorname{dist}(\hat{\bm{x}},\bm{x})/\left\|\bm{x}\right\|$ falls below $10^{-5}$.

Figure 2 shows that around $4.5n$ Gaussian phaseless measurements suffice for exact recovery with high probability via the Wirtinger flow algorithm. We also see that about six octanary patterns are sufficient.

3 Performance on natural images

We move on to testing the Wirtinger flow algorithm on various images of different sizes; these are photographs of the Naqsh-e Jahan Square in the central Iranian city of Esfahan, the Stanford main quad, and the Milky Way galaxy. Since each photograph is in color, we run the WF algorithm on each of the three RGB images that make up the photograph. Color images are viewed as $n_{1}\times n_{2}\times 3$ arrays, where the first two indices encode the pixel location, and the last the color band.

We generate $L=20$ random octanary patterns and gather the coded diffraction patterns for each color band using these 20 samples. As before, we run 50 iterations of the power method as the initialization step. The updates use the sequence $\mu_{\tau}=\min(1-\exp({-\tau/\tau_{0}}),0.4)$ where $\tau_{0}\approx 330$ as before. In all cases we run $300$ iterations and record the relative recovery error as well as the running time. If $\bm{x}$ and $\hat{\bm{x}}$ are the original and recovered images, the relative error is equal to $\left\|\hat{\bm{x}}-\bm{x}\right\|/\left\|\bm{x}\right\|$, where $\left\|\cdot\right\|$ is the Euclidean norm $\left\|\bm{x}\right\|^{2}=\sum_{i,j,k}|x(i,j,k)|^{2}$. The computational time we report is the the computational time averaged over the three RGB images. All experiments were carried out on a MacBook Pro with a 2.4 GHz Intel Core i7 Processor and 8 GB 1600 MHz DDR3 memory.

Figure 3 shows the images recovered via the Wirtinger flow algorithm. In all cases, WF gets 12 or 13 digits of precision in a matter of minutes. To convey an idea of timing that is platform-independent, we also report time in units of FFTs; one FFT unit is the amount of time it takes to perform a single FFT on an image of the same size. Now all the workload is dominated by matrix vector products of the form $\bm{A}\bm{z}$ and $\bm{A}^{*}\bm{v}$. In details, each iteration of the power method in the initialization step, or each update (2.2) requires $40$ FFTs; the factor of 40 comes from the fact that we have 20 patterns and that each iteration involves one FFT and one adjoint or inverse FFT. Hence, the total number of FFTs is equal to

Another way to state this is that the total workload of our algorithm is roughly equal to $350$ applications of the sensing matrix $\bm{A}$ and its adjoint $\bm{A}^{*}$. For about 13 digits of accuracy (relative error of about $10^{-13}$), Figure 3 shows that we need between 21,000 and 42,000 FFT units. This is within a factor between 1.5 and 3 of the optimal number computed above. This increase has to do with the fact that in our implementation, certain variables are copied into other temporary variables and these types of operations cause some overhead. This overhead is non-linear and becomes more prominent as the size of the signal increases.

For comparison, SDP based solutions such as PhaseLift and PhaseCut would be prohibitive on a laptop computer as the lifted signal would not fit into memory. In the SDP approach an $n$ pixel image become an $n^{2}/2$ array, which in the first example already takes storing the lifted signal even for the smallest image requires $(189\times 768)^{2}\times 1/2\times 8$ Bytes, which is approximately 85 GB of space. (For the image of the Milky Way, storage would be about 17 TB.) These large memory requirements prevent the application of full-blown SDP solvers on desktop computers.

4 3D molecules

Understanding molecular structure is a great contemporary scientific challenge, and several techniques are currently employed to produce 3D images of molecules; these include electron microscopy and X-ray imaging. In X-ray imaging, for instance, the experimentalist illuminates an object of interest, e.g. a molecule, and then collects the intensity of the diffracted rays, please see Figure 4 for an illustrative setup. Figures 5 and 6 show the schematic representation and the corresponding electron density maps for the Caffeine and Nicotine molecules: the density map $\rho(x_{1},x_{2},x_{3})$ is the 3D object we seek to infer. In this paper, we do not go as far 3D reconstruction but demonstrate the performance of the Wirtinger flow algorithm for recovering projections of 3D molecule density maps from simulated data. For related simulations using convex schemes we refer the reader to .

To derive signal equations, consider an experimental apparatus as in Figure 4. If we imagine that light propagates in the direction of the $x_{3}$-axis, an approximate model for the collected data reads

Therefore, by imputing the missing phase using phase retrieval algorithms, one can recover a slice of the 3D Fourier transfom of the electron density map, i.e. $R(f_{1},f_{2},0)$. Viewing the object from different angles or directions gives us different slices. In a second step we do not perform in this paper, one can presumably recover the 3D Fourier transform of the electron density map from all these slices (this is the tomography or blind tomography problem depending upon whether or not the projection angles are known) and, in turn, the 3D electron density map.

We now generate $51$ observation planes by rotating the $x_{1}x_{2}$-plane around the $x_{1}$-axis by equally spaced angles in the interval $[0,2\pi]$. Each of these planes is associated with a 2D projection of size $1024\times 1024$, giving us 20 coded diffraction octanary patterns (we use the same patterns for all 51 projections). We run the Wirtinger flow algorithm with exactly the same parameters as in the previous section, and stop after 150 gradient iterations. Figure 8 reports the average relative error over the 51 projections and the total computational time required for reconstructing all 51 images.

Theory for the Coded Diffraction Model

We complement our study with theoretical results applying to the model of coded diffraction patterns. These results concern a variation of the Wirtinger flow algorithm: whereas the iterations are exactly the same as (2.2), the initialization applies an iterative scheme which uses fresh sets of sample at each iteration. This is described in Algorithm 2. In the CDP model, the partitioning assigns to the same group all the observations and sampling vectors corresponding to a given realization of the random code. This is equivalent to partitioning the random patterns into $B+1$ groups. As a result, sampling vectors in distinct groups are stochastically independent.

with probability at least $1-(4L+2)/{n^{3}}$. Further, take a constant learning parameter sequence, $\mu_{\tau}=\mu$ for all $\tau=1,2,\ldots$ and assume $\mu\leq c_{1}$ for some fixed numerical constant $c_{1}$. Then there is an event of probability at least $1-(2L+1)/{n^{3}}-{1}/{n^{2}}$, such that on this event, starting from any initial solution $\bm{z}_{0}$ obeying (5.1), we have

In the Gaussian model, both statements also hold with high probability provided that $m\geq c_{0}\cdot n\,(\log n)^{2}$, where $c_{0}$ is a sufficiently large numerical constant.

Hence, we achieve perfect recovery from on the order of $n(\log n)^{4}$ samples arising from a coded diffraction experiment. Our companion paper established that PhaseLift—the SDP relaxation—is also exact with a sampling complexity on the order of $n(\log n)^{4}$ (this has recently been improved to $n(\log n)^{2}$ ). We believe that the sampling complexity of both approaches (WF and SDP) can be further reduced to $n\log n$ (or even $n$ for certain kind of patterns). We leave this to future research.

Setting $\mu=c_{1}$ yields $\epsilon$ accuracy in $\mathcal{O}(\log 1/\epsilon)$ iterations. As the computational work at each iteration is dominated by two matrix-vector products of the form $\bm{A}\bm{z}$ and $\bm{A}^{*}\bm{v}$, it follows that the overall computational is at most $\mathcal{O}(nL\log n\log{1}/{\epsilon})$. In particular, this approach yields a near-linear time algorithm in the CDP model (linear in the dimension of the signal $n$). In the Gaussian model, the complexity scales like $\mathcal{O}(mn\log{1}/{\epsilon})$.

Wirtinger Derivatives

of several complex variables can be written in the form $f(\bm{z},\bar{\bm{z}})$, where $f$ is holomorphic in $\bm{z}=\bm{x}+i\bm{y}$ for fixed $\bar{\bm{z}}$ and holomorphic in $\bar{\bm{z}}=\bm{x}-i\bm{y}$ for fixed $\bm{z}$. This holds as long as a the real-valued functions $u$ and $v$ are differentiable as functions of the real variables $\bm{x}$ and $\bm{y}$. As an example, consider

This fact underlies the development of the Wirtinger calculus. In essence, the conjugate coordinates

Our definitions follow standard notation from multivariate calculus so that derivatives are row vectors and gradients are column vectors. In this new coordinate system the complex gradient is given by

In this coordinate system the complex Hessian is given by

If we were to run gradient descent in this new coordinate system, the iterates would be

Note that when $f$ is a real-valued function (as in this paper) we have

Therefore, the second set of updates in (6.1) is just the conjugate of the first. Thus, it is sufficient to keep track of the first update, namely,

For real valued functions of complex variables, setting

Proofs

We first note that in the CDP model with admissible CDPs $\left\|\bm{a}_{r}\right\|\leq\sqrt{6n}$ for all $r=1,2,\ldots,m$, as in our CDP model $\left|d\right|\leq\sqrt{3}<\sqrt{6}$. In the Gaussian model the measurements vectors also obey $\left\|\bm{a}_{r}\right\|\leq\sqrt{6n}$ for all $r=1,2,\ldots,m$ with probability at least $1-me^{-1.5n}$. Throughout the proofs, we assume we are on this event without explicitly mentioning it each time.

Before we begin with the proofs we should mention that we will prove our result using the update

Since $\left|\left\|\bm{z}_{0}\right\|^{2}-\left\|\bm{x}\right\|^{2}\right|\leq\frac{1}{64}\left\|\bm{x}\right\|^{2}$ holds with high probability as proven in Section 7.8, we have

Therefore, the results for the update (7.1) automatically carry over to the WF update with a simple rescaling of the upper bound on the learning parameter. More precisely, if we prove that the update (7.1) converges to a global optimum as long as $\mu\leq\mu_{0}$, then the convergence of the WF update to a global optimum is guaranteed as long as $\mu_{\text{WF}}\leq\frac{63}{64}\mu_{0}$. Also, the update in (7.1) is invariant to the Euclidean norm of $\bm{x}$. Therefore, without loss of generality we will assume throughout the proofs that $\left\|\bm{x}\right\|=1$.

We remind the reader that throughout $\bm{x}$ is a solution to our quadratic equations, i.e. obeys $y=|\bm{A}\bm{x}|^{2}$ and that the sampling vectors are independent from $\bm{x}$. Define

to be the set of all vectors that differ from the planted solution $\bm{x}$ only by a global phase factor. We also introduce the set of all points that are close to $P$,

2 Formulas for the complex gradient and Hessian

We gather some useful gradient and Hessian calculations that will be used repeatedly. Starting with

3 Expectation and concentration

This section gathers some useful intermediate results whose proofs are deferred to Appendix A. The first two lemmas establish the expectation of the Hessian, gradient and a related random variable in both the Gaussian and admissible CDP models.In the CDP model the expectation is with respect to the random modulation pattern.

Recall that $\bm{x}$ is a solution obeying $\left\|\bm{x}\right\|=1$, which is independent from the sampling vectors. Furthermore, assume the sampling vectors $\bm{a}_{r}$ are distributed according to either the Gaussian or admissible CDP model. Then

The next lemma gathers some useful identities in the Gaussian model.

The next lemma establishes the concentration of the Hessian around its mean for both the Gaussian and the CDP model.

In the setup of Lemma 7.1, assume the vectors $\bm{a}_{r}$ are distributed according to either the Gaussian or admissible CDP model with a sufficiently large number of measurements. This means that the number of samples obeys $m\geq c(\delta)\cdot n\log n$ in the Gaussian model and the number of patterns obeys $L\geq c(\delta)\cdot\log^{3}n$ in the CDP model. Then

holds with probability at least $1-10e^{-\gamma n}-{8}/{n^{2}}$ and $1-(2L+1)/{n^{3}}$ for the Gaussian and CDP models, respectively.

We will also make use of the two results below, which are corollaries of the three lemmas above. These corollaries are also proven in Appendix A.

The next lemma establishes the concentration of the gradient around its mean for both Gaussian and admissible CDP models.

holds with probability at least $1-20e^{-\gamma m}-{4m}/{n^{4}}$ in the Gaussian model and $1-(4L+2)/{n^{3}}$ in the CDP model.

We finish with a result concerning the concentration of the sample covariance matrix.

holds with probability at least $1-2e^{-\gamma m}$ for the Gaussian model and $1-{1}/{n^{2}}$ in the CDP model. On this event,

4 General convergence analysis

We will assume that the function $f$ satisfies a regularity condition on $E(\epsilon)$, which essentially states that the gradient of the function is well behaved. We remind the reader that $E(\epsilon)$, as defined in (7.4), is the set of points that are close to the path of global minimizers.

We say that the function $f$ satisfies the regularity condition or $RC(\alpha,\beta,\epsilon)$ if for all vectors $\bm{z}\in E(\epsilon)$ we have

In the lemma below we show that as long as the regularity condition holds on $E(\epsilon)$ then Wirtinger Flow starting from an initial solution in $E(\epsilon)$ converges to a global optimizer at a geometric rate. Subsequent sections shall establish that this property holds.

Assume that $f$ obeys RC$(\alpha,\beta,\epsilon)$ for all $\bm{z}\in E(\epsilon)$. Furthermore, suppose $\bm{z}_{0}\in E$, and assume $0<\mu\leq{2}/{\beta}$. Consider the following update

Then for all $\tau$ we have $\bm{z}_{\tau}\in E(\epsilon)$ and

We note that for $\alpha\beta<4$, (7.10) holds with the direction of the inequality reversed.One can see this by applying Cauchy-Schwarz and calculating the determinant of the resulting quadratic form. Thus, if $RC(\alpha,\beta,\epsilon)$ holds, $\alpha$ and $\beta$ must obey $\alpha\beta\geq 4$. As a result, under the stated assumptions of Lemma 7.10 above, the factor $1-2\mu/\alpha\geq 1-4/(\alpha\beta)$ is non-negative.

Proof The proof follows a structure similar to related results in the convex optimization literature e.g. [53, Theorem 2.1.15]. However, unlike these classical results where the goal is often to prove convergence to a unique global optimum, the objective function $f$ does not have a unique global optimum. Indeed, in our problem, if $\bm{x}$ is solution, then $e^{i\phi}\bm{x}$ is also solution. Hence, proper modification is required to prove convergence results.

We prove that if $\bm{z}\in E(\epsilon)$ then for all $0<\mu\leq{2}/{\beta}$

Therefore, if $\bm{z}\in E(\epsilon)$ then we also have $\bm{z}_{+}\in E(\epsilon)$. The lemma follows by inductively applying (7.11). Now simple algebraic manipulations together with the regularity condition (7.10) give

where the last line follows from $\mu\leq{2}/{\beta}$. The definition of $\phi(\bm{z}_{+})$ gives

5 Proof of the regularity condition

For any $\bm{z}\in E(\epsilon)$, we need to show that

We prove this with $\delta=0.01$ by establishing that our gradient satisfies the local smoothness and local curvature conditions defined below. Combining both these two properties gives (7.12).

We say that the function $f$ satisfies the local curvature condition or $LCC(\alpha,\epsilon,\delta)$ if for all vectors $\bm{z}\in E(\epsilon)$,

This condition essentially states that the function curves sufficiently upwards (along most directions) near the curve of global optimizers.

We say that the function $f$ satisfies the local smoothness condition or $LSC(\beta,\epsilon,\delta)$ if for all vectors $\bm{z}\in E(\epsilon)$ we have

This condition essentially states that the gradient of the function is well behaved (the function does not vary too much) near the curve of global optimizers.

6 Proof of the local curvature condition

For any $\bm{z}\in E(\epsilon)$, we want to prove the local curvature condition (7.13). Recall that

and define $\bm{h}:=e^{-i\phi(\bm{z})}\bm{z}-\bm{x}$. To establish (7.13) it suffices to prove that

holds for all $\bm{h}$ satisfying $\operatorname{Im}(\bm{h}^{*}\bm{x})=0$, $\|\bm{h}\|_{2}\leq\epsilon$. Equivalently, we only need to prove that for all $\bm{h}$ satisfying $\operatorname{Im}(\bm{h}^{*}\bm{x})=0$, $\|\bm{h}\|_{2}=1$ and for all $s$ with $0\leq s\leq\epsilon$,

holds for all $\bm{h}$ obeying $\left\|\bm{h}\right\|=1$. Therefore, to establish the local curvature condition (7.13) it suffices to show that

We will establish (7.17) for different measurement models and different values of $\epsilon$. Below, it shall be convenient to use the shorthand

Set $\epsilon={1}/{8\sqrt{n}}$. We show that with high probability, (7.17) holds for all $\bm{h}$ satisfying Im$(\bm{h}^{*}\bm{x})=0$, $\|\bm{h}\|_{2}=1$, $0\leq s\leq\epsilon$, $\delta\leq 0.01$, and $\alpha\geq 30$. First, note that by Cauchy-Schwarz inequality,

The last inequality follows from $(a-b)^{2}\geq\frac{a^{2}}{2}-b^{2}$. By Corollary 7.5,

holds with high probability for all $\bm{h}$ obeying $\left\|\bm{h}\right\|=1$. Furthermore, by applying Lemma 7.8,

holds with high probability. Plugging (7.19) and (7.20) in (7.6.1) yields

(7.17) follows by using $\alpha\geq 30$, $\epsilon=\frac{1}{8\sqrt{n}}$ and $\delta=0.01$.

6.2 Proof of (7.17) with ϵ=1/8italic-ϵ18\epsilon={1}/{8} in the Gaussian model

Set $\epsilon={1}/{8}$. We show that with high probability, (7.17) holds for all $\bm{h}$ satisfying $\operatorname{Im}(\bm{h}^{*}\bm{x})=0$, $\|\bm{h}\|_{2}=1$, $0\leq s\leq\epsilon$, $\delta\leq 2$, and $\alpha\geq 8$. To this end, we first state a result about the tail of a sum of i.i.d. random variables. Below, $\Phi$ is the cumulative distribution function of a standard normal variable.

To establish (7.17) we first prove it for a fixed $\bm{h}$, and then use a covering argument. Observe that

compare (7.5) and (7.6). Therefore, using $s\leq\frac{1}{8}$,

We have all the elements to apply Lemma 7.13 with $\sigma^{2}=m\max(9^{2},210)=210m$ and $y={m}/{4}$:

with $\gamma=1/1260$. Therefore, with probability at least $1-e^{-3\gamma m}$, we have

holds with probability at least $1-e^{\gamma m}$ as long as $m\geq c\cdot n\log n$. This concludes the proof of (7.17) with $\alpha\geq 8$.

7 Proof of the local smoothness condition

Define $\bm{h}:=e^{-\phi(\bm{z})}\bm{z}-\bm{x}$ and $\bm{w}:=e^{-i\phi(\bm{z})}\bm{u}$, to establish (7.14) it suffices to prove that

holds for all $\bm{h}$ and $\bm{w}$ satisfying Im$(\bm{h}^{*}\bm{x})=0$, $\left\|\bm{h}\right\|\leq\epsilon$ and $\left\|\bm{w}\right\|=1$. Equivalently, we only need to prove for all $\bm{h}$ and $\bm{w}$ satisfying Im$(\bm{h}^{*}\bm{x})=0$, $\left\|\bm{h}\right\|=\left\|\bm{w}\right\|=1$ and $\forall s:0\leq s\leq\epsilon$,

Note that since $(a+b+c)^{2}\leq 3(a^{2}+b^{2}+c^{2})$

We now bound each of the terms on the right-hand side. For the first term we use Cauchy-Schwarz and Corollary 7.6, which give

Finally, for the third term we use the Cauchy-Schwarz inequality together with Lemma 7.8 (inequality) (7.9) to derive

We now plug these inequalities into (7.7) and get

which completes the proof of (7.26) and, in turn, establishes the local smoothness condition in (7.14). However, the last line of (7.7) holds as long as

In our theorems we use two different values of $\epsilon=\frac{1}{8\sqrt{n}}$ and $\epsilon=\frac{1}{8}$. Using $\delta\leq 0.01$ in (7.32) we conclude that the local smoothness condition (7.7) holds as long as

8 Wirtinger flow initialization

In this section, we prove that the WF initialization obeys (3.1) from Theorem 3.3. Recall that

Also, since $\lambda_{0}$ is the top eigenvalue of $\bm{Y}$, and $\left\|\bm{x}\right\|=1$, we have

Combining the above two inequalities together, we have

9 Initialization via resampled Wirtinger Flow

In this section, we prove that that the output of Algorithm 2 obeys (5.1) from Theorem 5.1. Introduce

We prove that $F$ obeys a regularization condition in $E(1/8)$, namely,

Lemma 7.7 implies that for a fixed vector $\bm{z}$,

holds with high probability. The last inequality follows from (7.33). Applying Lemma 7.7, we also have

It only remains to establish (7.33). First, without loss of generality, we can assume that $\phi(\bm{z})=0$, which implies $\operatorname{Re}(\bm{z}^{*}\bm{x})=\left|\bm{z}^{*}\bm{x}\right|$ and use $\left\|\bm{z}-\bm{x}\right\|$ in lieu of dist$(\bm{z},\bm{x})$. Set $\bm{h}:=\bm{z}-\bm{x}$ so that $\operatorname{Im}(\bm{x}^{*}\bm{h})=0$. This implies

where the last inequality is due to $\left\|\bm{h}\right\|\leq\epsilon\leq{1}/{8}$. Furthermore,

where the last inequality also holds because $\left\|\bm{h}\right\|\leq\epsilon\leq{1}/{8}$. Finally, (7.35) and (7.9) imply

where $\alpha^{\prime}=8$ and $\beta^{\prime}=200$.

Acknowledgements

E. C. is partially supported by AFOSR under grant FA9550-09-1-0643 and by gift from the Broadcom Foundation. X. L. is supported by the Wharton Dean’s Fund for Post-Doctoral Research and by funding from the National Institutes of Health. M. S. is supported by a Benchmark Stanford Graduate Fellowship and AFOSR grant FA9550-09-1-0643. M. S. would like to thank Andrea Montanari for helpful discussions and for the class which inspired him to explore provable algorithms based on non-convex schemes. He would also like to thank Alexandre d’Aspremont, Fajwel Fogel and Filipe Maia for sharing some useful code regarding 3D molecule reconstructions, Ben Recht for introducing him to reference , and Amit Singer for bringing the paper to his attention.

References

Appendix A Expectations and deviations

The proof for admissible coded diffraction patterns follows from Lemmas 3.1 and 3.2 in . For the Gaussian model, it is a consequence of the two lemmas below, whose proofs are ommitted.

A.2 Proof of Lemma 7.2

Thus by applying Lemma 3.1 in (for the CDP model) and Lemma A.1 above (for the Gaussian model) we have

A.3 Proof of Lemma 7.3

The identity (7.7) follows from standard normal moment calculations.

A.4 Proof of Lemma 7.4

where $R$ is a positive scalar whose value will be determined shortly, and define the events

A.4.2 The Gaussian model

By unitary invariance, we may take $\bm{x}=\bm{e}_{1}$. Letting $\bm{z}(1)$, be the first coordinate of a vector $\bm{z}$, to prove Lemma 7.4 for the Gaussian model it suffices to prove the two inequalities,

For any $\epsilon>0$, there is a constant $C>0$ with the property that $m\geq C\cdot n$ implies

with probability at least $1-3n^{-2}$; this is a consquence of Chebyshev’s inequality. Moreover a union bound gives

We now turn our attention to the last two terms of (A.4.2). For the second term, the ordinary Hoeffding’s inequality (Proposition 10 in ) gives that for any constants $\delta_{0}$ and $\gamma$, there exists a constant $C(\delta_{0},\gamma)$, such that for $m\geq C(\delta_{0},\gamma)\sqrt{n\left(\sum_{r=1}^{m}\left|\bm{a}_{r}(1)\right|^{6}\right)}$,

holds with probability at least $1-3e^{-2\gamma n}$. To control the final term, we apply the Bernstein-type inequality (Proposition 16 in ) to assert the following: for any positive constants $\delta_{0}$ and $\gamma$, there exists a constant $C(\delta_{0},\gamma)$, such that for $m\geq C(\delta_{0},\gamma)(\sqrt{n(\sum_{r=1}^{m}\left|\bm{a}_{r}(1)\right|^{4})}+n\max_{r=1}^{m}|\bm{a}_{r}(1)|^{2})$,

holds with probability at least $1-2e^{-2\gamma n}$.

Therefore, for any unit norm vector $\bm{y}$,

holds with probability at least $1-5e^{-2\gamma n}$. By Lemma 5.4 in , we can bound the operator norm via an $\epsilon$-net argument:

The proof of (A.3) is similar. The only difference is that the random matrix is not Hermitian, so we work with

where $\bm{u}$ and $\bm{v}$ are unit vectors.

A.5 Proof of Corollary 7.5

The other inequality is established in a similar fashion.

A.6 Proof of Corollary 7.6

In the proof of Lemma 7.4, we established that with high probability,

This concludes the proof of one side. The other side is similar.

A.7 Proof of Lemma 7.7

Define $\bm{h}:=e^{-\phi(\bm{z})}\bm{z}-\bm{x}$ and $\bm{w}:=e^{-i\phi(\bm{z})}\bm{u}$, we have

Combining (A.7) and (A.8) together with the triangular inequality and Lemma 7.4 give

A.8 Proof of Lemma 7.8

The result for the CDP model follows from Lemma 3.3 in . For the Gaussian model, it is a consequence of standard results, e.g. Theorem 5.39 in , concerning the deviation of the sample covariance matrix from its mean.

Appendix B The Power Method

We use the power method (Algorithm 3) with a random initialization to compute the first eigenvector of $\bm{Y}=\bm{A}\,\text{diag}\{\bm{y}\}\,\bm{A}^{*}$. Since, each iteration of the power method asks to compute the matrix-vector product

we simply need to apply $\bm{A}$ and $\bm{A}^{*}$ to an arbitrary vector. In the Gaussian model, this costs $2mn$ multiplications while in the CDP model the cost is that of $2L$ $n$-point FFTs. We now turn our attention to the number of iterations required to achieve a sufficiently accurate initialization.

Standard results from numerical linear algebra show that after $k$ iterations of the power method, the accuracy of the eigenvector is $\mathcal{O}(\tan\theta_{0}(\lambda_{2}/\lambda_{1})^{k})$, where $\lambda_{1}$ and $\lambda_{2}$ are the top two eigenvalues of the positive semidefinite matrix $\bm{Y}$, and $\theta_{0}$ is the angle between the initial guess and the top eigenvector. Hence, we would need on the order of $\log(n/\epsilon)\,/\,\log(\lambda_{1}/\lambda_{2})$ for $\epsilon$ accuracy. Under the stated assumptions, Lemma 7.4 bounds below the eigenvalue gap by a numerical constant so that we can see that few iterations of the power method would yield accurate estimates.