Non-Convex Phase Retrieval from STFT Measurements

I Introduction

The problem of recovering a signal from its Fourier transform magnitude arises in many areas in engineering and science, such as optics, X-ray crystallography, speech recognition, blind channel estimation, alignment and astronomy . This problem is called Fourier phase retrieval and can be viewed as a special case of a quadratic system of equations. The latter area received considerable attention recently, partially due to its strong connections with the fields of compressed sensing and matrix completion; see for instance . Contemporary surveys of the phase retrieval problem from a signal processing point of view can be found in .

Phase retrieval for one-dimensional (1D) signals is an ill-posed problem unless the signal has the minimum phase property . In this special case, the signal can be recovered by several tractable algorithms (see for instance Section 2.6 of ). Particularly, in it was shown that a semidefinite program (SDP) relaxation achieves the optimal solution in the least-squares (LS) sense. For general signals, two main approaches are typically suggested. The first builds upon prior knowledge on the signal’s support, such as sparsity or a portion of the underlying signal . An alternative strategy makes use of additional measurements. Such measurements can be obtained by structured illuminations and masks or by measuring the magnitude of the short-time Fourier transform (STFT) . In , it was demonstrated that for the same number of measurements, the STFT magnitude leads to better performance than an over-sampled discrete Fourier transform (DFT).

where $k=0,\dots,N-1$, $m=0,\dots,\left\lceil\frac{N}{L}\right\rceil-1$ and $L$ determines the separation in time between adjacent sections. The pseudo-inverse of the STFT is given by

The problem of recovering a signal from its STFT magnitude $|\mathbf{X}[m,k]|^{2}$, frequently called spectrogram, arises in several applications in optics and speech processing . Particularly, it serves as the model for a popular variant of an ultra-short laser pulse characterization technique called Frequency-Resolved Optical Gating (referred to as X-FROG) . Another important application is ptychography in which a moving probe is used to sense multiple diffraction measurements .

Several algorithms were suggested to recover a signal from the magnitude of its STFT. The classical method, proposed by Griffin and Lim , is a modification of the alternating projection (or reduction error) algorithms of Gerchberg and Saxton and Fienup . The properties of this method are not well understood (for analysis of alternating projection algorithms in phase retrieval, see ). In , the authors prove that a non-vanishing signal can be recovered by an SDP with maximal overlap between adjacent windows ($L=1$). They also demonstrate empirically that the algorithm works well with less restrictive requirements on the window and is robust to noise. Despite the appealing numerical performance, solving an SDP requires high computational resources. Recently, an interesting recovery approach was proposed in . This paper suggests a multi-stage method, based on spectral clustering and phase synchronization. It is shown that the algorithm achieves stable estimation (and exact in the noise-free setting) with only $\mathcal{O}(N\log N)$ phaseless STFT measurements. However, this technique requires a random window of length $W=N$, while in most applications it is common to work with shorter windows. Another line of works suggest applying a phase synchronization framework . It was shown that even for short windows, the sought signal can be recovered exactly and efficiently by spectral and greedy techniques. These methods are accompanied by stability guarantees. Their main drawback is that they rely on reliable estimates of the temporal magnitude, which do not always exist.

Here, we take a different approach and propose a data-driven initialization technique, followed by non-convex gradient algorithms. We begin by taking the 1D DFT of the acquired data with respect to the frequency variable (the second variable of the STFT). This transformation reveals the underlying structure of the data and greatly simplifies the analysis. As a direct consequence, we show that for $L=1$ and sufficiently long windows $W\geq\left\lceil\frac{N+1}{2}\right\rceil$ (and some mild additional conditions), one can recover the signal by extracting the principal eigenvector of a designed matrix, constructed as the solution of a simple linear LS problem. We refer to this matrix as the approximation matrix since it approximates the correlation matrix $\mathbf{X}:=\mathbf{x}\mathbf{x}^{*}$.

When the conditions for a closed-form solution are not met, we propose using the principal eigenvector of the approximation matrix to initialize two non-convex algorithms. The first is based on minimizing a standard quadratic loss function, frequently called the empirical risk (ER). Inspired by the phasecut method , we also propose a new phase retrieval algorithm, called Non-Convex PhaseCut (NCPC), that maximizes a quadratic function over the set of phases. Each step of the algorithm follows the component of the gradient which agrees with the phase constraints. As will be shown, the ER technique is more stable in the low signal–to–noise ratio (SNR) regimes, while NCPC is superior in high SNR environments and for short windows. Our approach deviates in two important aspects from the recent line of work in non-convex phase retrieval . First, all these papers focus their attention on the setup of phase retrieval with random sensing vectors and rely heavily on probabilistic considerations. In this case, efficient algorithms were designed to estimate the signal from $\mathcal{O}(N)$ measurements. In contrast, we consider a deterministic framework. Second, we construct our approximation matrix by the solution of a LS problem, whereas the aforementioned papers take a superposition of the measurements to approximate $\mathbf{X}$.

The properties of non-convex algorithms depend heavily on the initialization method and the geometry of the loss functions. For $L=1$, we estimate the distance between the proposed initialization and the target signal, which decays to zero as $W$ tends to $\frac{N+1}{2}$. If the signal has unit modulus entries, then a slight modification of our initialization recovers the signal exactly for $W\geq 2$. In the later case, we also prove the existence of a basin of attraction around the global minimum of the ER loss function and estimate its size. In the basin of attraction, the algorithm is guaranteed to converge to a global minimum at a geometric rate. We note that while the theoretical guarantees of the algorithms are limited, their experimental performance is significantly better. Particularly, the algorithms perform well with small redundancy in the measurements and are robust in the presence of noise.

The paper is organized as follows. We begin in Section II by formulating mathematically the problem of phase retrieval from STFT magnitude measurements. In Section III we discuss the uniqueness of the solution and present conditions under which it has a closed-form LS expression. Additionally, we present a method that recovers signals with unit modulus entries under mild conditions. Section IV presents the two non-convex algorithms with the proposed initialization. Section V shows numerical results and Section VI presents our theoretical findings regarding the proposed initialization and the ER loss function. Proofs are provided in Section VII. Section VIII concludes the paper, discusses its main implications and draws potential future research directions.

II Problem Formulation

The distance between two vectors is defined as

If $d\left(\mathbf{z},\mathbf{x}\right)=0$ then we say that $\mathbf{x}$ and $\mathbf{z}$ are equal up to global phase. The phase $\phi\in[0,2\pi)$ attaining the minimum is denoted by $\phi(\mathbf{z})$, i.e.,

In the sequel, we make use of the notion of non-vanishing signals, defined as follows:

Instead of treating the measurements (II.1) directly, we often consider the acquired data in a transformed domain by taking its 1D DFT with respect to the frequency variable (normalized by $1/N$). Then, our measurement model reads

and $\mathbf{f}_{k}^{*}$ is the $k$th row of the DFT matrix.

Our problem of recovering $\mathbf{x}$ from the measurements (II.1) can therefore be posed as a constrained LS problem:

III Uniqueness and Basic Algorithms

A fundamental question in phase retrieval problems is whether the quadratic measurement operator of (II.1), or equivalently the non-convex problem (II.8), determines the underlying signal $\mathbf{x}$ uniquely (up to global phase, see Definition II.1). In other words, one wants to know the conditions on the window $\mathbf{g}$ and the signal $\mathbf{x}$ such that the non-linear transformation that maps $\mathbf{x}$ to $\mathbf{Z}$ is injective. Before treating this question, we introduce some basic window definitions:

A window $\mathbf{g}$ is called a rectangular window of length $W$ if $\mathbf{g}[n]=1$ for all $n=0,\dots,W-1$ and zero elsewhere. It is a non-vanishing window of length $W$ if $\mathbf{g}[n]\neq 0$ for all $n=0,\dots,W-1$ and zero elsewhere.

An important example for an admissible window is a rectangular window. Specifically, we have the following lemma:

A rectangular window $\mathbf{g}$ of length $2\leq W\leq N/2$ is an admissible window of length $W$ if $\alpha$ and $N$ are co-prime numbers for all $\alpha=2,\dots,W$. This holds trivially when $N$ is a prime number.

We repeated this process 100 times for several values of $W$. As can be seen in Table I, $|\lambda_{\min}|$ is bounded away from zero, implying that the windows are indeed admissible.

We now analyze the uniqueness of the measurement operator for the case $L=1$. Uniqueness results for $L>1$ are discussed in . Our results are constructive in the sense that their proofs provide an explicit scheme to recover the signal.

Our first result concerns non-vanishing signals. In this case, the magnitude of the STFT determines the underlying signal uniquely under mild conditions. This conclusion was already derived in based on different considerations. Nevertheless, the following proposition comes with an explicit recovery scheme as presented in Appendix -A.

A similar uniqueness result was derived in . There, it is required that the DFT of $|\mathbf{g}[n]|^{2}$ is non-vanishing, $N\geq 2W-1$ and $N$ and $W-1$ are co-prime numbers.

In the special case in which the signal is known to have unit modulus entries, the signal can be recovered as the principal eigenvector of a matrix designed as follows:

where $P:=\{n\thinspace:\thinspace(\mathbf{G}_{0}^{\dagger}\mathbf{y}_{0})[n]>0\}$. If $\mathbf{G}_{0}$ is invertible then

The following proposition shows that Algorithm 1 recovers the underlying signal for $L=1$ if the window is sufficiently long and satisfies some additional technical conditions. In , an equivalent uniqueness result was derived but without providing an algorithm. Algorithm 1 is equivalent to the discretized version of Wigner deconvolution that was suggested previously without theoretical analysis in .

Let $L=1$ and suppose that $\mathbf{g}$ is an admissible window of length $W\geq\left\lceil\frac{N+1}{2}\right\rceil$ (see Definition III.2). Then, Algorithm 1 recovers any complex signal uniquely up to global phase.

IV Local Non-Convex Algorithms

In this section we present our main algorithmic approach to recover a signal from its STFT magnitude (II.1). First, we propose two non-convex gradient algorithms to estimate the signal. As the problem is inherently non-convex, we then suggest a systematic, data–driven, technique for initialization. This non-convex approach for STFT phase retrieval is summarized in Algorithm 2. The code for all algorithms is publicly available at http://webee.technion.ac.il/Sites/People/YoninaEldar.

The equality between the two loss functions is proven in Appendix -D. In the sequel, we use both formulations.

Figure IV.1 presents the two-dimensional (first two variables) plane of the loss function (IV.1) for the signal $\mathbf{x}=[0.2,0.2,0,0,0]$ (i.e., $N=5$) with $L=1$ and a rectangular window of length $W=2$. The function has no sharp transitions and contains two saddle points and two global minima (as a result of the global phase ambiguity). Accordingly, in this specific case and bearing in mind that our view is restricted to two of the five dimensions only, it seems that a gradient descent algorithm will converge to a global minimum from almost any initialization (see also ). While this phenomenon does not occur for any arbitrary parameter selection, this example motivates applying a gradient algorithm directly on the non-convex loss function (for a similar demonstration of the loss function with random sensing vectors, see ).

One way to minimize the ER loss function (IV.1) or (IV.2) is by employing a gradient algorithm, where the $k$th iteration takes on the form

for step size $\mu$. For real signals, direct computation of the gradient in (IV.1) gives

Similar computations can be performed for (IV.2). If the signal is complex, then one can use the elegant formulation of Wirtinger derivatives, see . The loss functions (IV.1) or (IV.2) may be minimized by many other methods. For instance, in Section V we employ a trust-region algorithm.

IV-B Non-Convex PhaseCut (NCPC)

When minimizing the empirical risk (IV.1) or (IV.2), the unknown signal itself is the optimization variable. Alternatively, we may take the point of view that the unknowns are the phases of the STFT measurements. Indeed, if these phases were known, then one could recover the signal by applying (I.2). We may therefore rework the problem into one where only the phases are variables .

Since $\mathbf{I}-\operatorname{STFT}\circ\operatorname{STFT}^{\dagger}$ is an orthogonal projector, this further simplifies into the following non-convex optimization problem over complex phases:

The term involving the identity operator $\mathbf{I}$ is constant under the constraints, so that the problem is equivalent to the following maximization problem:

Notice that $\operatorname{STFT}\circ\operatorname{STFT}^{\dagger}$ is the orthogonal projector onto the subspace of matrices which are the STFT of some signal. As a result, applying $\operatorname{STFT}\circ\operatorname{STFT}^{\dagger}$ to the matrix $\mathbf{Z}^{1/2}\odot\mathbf{U}$ produces the matrix which, in the LS sense, is closest to being the STFT of a signal. Thus, the cost function in (IV.5) favors phases $\mathbf{U}$ such that $\mathbf{Z}^{1/2}\odot\mathbf{U}$ is as close as possible to an STFT. We recall that this projection operator can be computed efficiently by applying (I.1) and (I.2) using FFT.

Problem (IV.5) resembles the phase synchronization problem . In , the authors pursue a convex relaxation of (IV.5) named phasecut. Here, following , we use the Manopt toolbox to run local optimization of (IV.5) over the manifold of phases . In its simplest form, the algorithm follows the gradient’s component which is consistent with the feasible set of solutions (see details below). To initialize the local optimization algorithm, we set $\mathbf{U}_{0}$ to be the phases of $\operatorname{STFT}\left(\mathbf{x}_{0}\right)$, where $\mathbf{x}_{0}$ is the initialization used by Algorithm 2. This approach is summarized in Algorithm 3.

For completeness, we provide a brief overview of step 2 of Algorithm 3, that is, optimization of the phases. We restrict attention to a simple Riemannian optimization algorithm, namely, the gradient ascent algorithm. See for details about the more sophisticated Riemannian trust-region method (RTR), which we use in practice.

The variable $\mathbf{U}$ lives on a smooth manifold, namely, the set of phases

which is a Cartesian product of unit circles in the complex plane (a torus). This smooth nonlinear space can be linearized about every point $\mathbf{U}$ by differentiating the constraints. This yields a linear subspace known as the tangent space to $\mathcal{M}$ at $\mathbf{U}$:

That is, it subtracts from each entry $\mathbf{V}[m,k]$ its component aligned with $\mathbf{U}[m,k]$. Explicitly, the Riemannian gradient is then given by

Now that we are equipped with a notion of gradient on the manifold, the only missing ingredient to implement a gradient ascent optimization algorithm is a means of moving away from a point (a current iterate) along a chosen tangent direction (here, the gradient vector), while remaining on the manifold $\mathcal{M}$. The standard tool to achieve this is known as a retraction [2, §4.1]. An obvious retraction for $\mathcal{M}$ is

The gradient ascent algorithm takes the form

As will be shown next, this is also the stagnation point of Fienup’s algorithm. This approach is summarized in Algorithm 4.

We stress that this approach is different from a projected gradient method. Indeed, in a projected gradient method, one would alternate between following the classical gradient $\nabla f(\mathbf{U})$ and projecting to $\mathcal{M}$ with the $\operatorname{phase}$ operator. That is, each iteration resembles (IV.8) with $\nabla f$ instead of $\operatorname{grad}f$. In contrast, the Riemannian gradient method follows the tangent part of the gradient, $\operatorname{grad}f$ (IV.7) and then projects onto $\mathcal{M}$. One advantage is that, close to convergence, the Riemannian gradient has small norm (as expected), whereas the classical gradient may still be large.

IV-B2 Relation to Fienup’s Algorithm

Our method can be compared with the classical Fienup algorithm for the STFT case, also called Griffin–Lim algorithm , as follows. One approach to optimize (IV.5), instead of the Riemannian gradient iterations that we describe in Algorithm 4, is an iterative technique called projected power method (PPM), or generalized power method . This algorithm iterates as the power method, with the difference that, at each iteration, it keeps only the phases of the current iterate. Specifically, the $k$th iteration is of the form

Similarly to NCPC, the algorithm stops when (IV.9) is satisfied. On the other hand, each iteration of Fienup’s algorithm takes on the form

Applying the operator $\operatorname{phase}\circ\operatorname{STFT}$ on the iterations of (IV.11) shows that it is equivalent to PPM through the mapping $\mathbf{U}_{k}=\operatorname{phase}\circ\operatorname{STFT}(\mathbf{x}_{k})$. In this sense, one can understand Fienup’s algorithm as a particular iterative method to solve the optimization problem (IV.5).

According to [14, Lemma 15], all fixed points of (IV.10)–and hence of (IV.11)–map to critical points of the optimization problem (IV.5), that is, they map to points $\mathbf{U}_{k}$ where the Riemannian gradient is zero. These are only the first-order necessary optimality conditions. Numerical experiments (not displayed here) show that some of the stable fixed points of (IV.11) map to critical points which do not satisfy the second-order necessary optimality conditions (their Riemannian Hessian admits a positive eigenvalue) and are therefore suboptimal. In contrast, such points would be unstable fixed points for any reasonable Riemannian optimization algorithm as confirmed in the same experiments. This distinction at least partially explains why the empirical performance of the NCPC algorithm is superior to that of Fienup’s algorithm, as demonstrated in Section V.

IV-C Initialization

In practical applications, a variety of approaches are used to initialize the refinement techniques. While the specific initialization method is application-dependent, these approaches can be broadly classified into two categories. The first is based on the structure of the expected signal. For instance, in some applications it is common to use a Gaussian pulse with random phases as an initial point . This, however, may lead to a phenomenon called model bias in which the estimate tends to capture characteristics of the model rather than the true signal. An alternative strategy, also used by commercial software, is based on random initialization. This is very different from our initialization which exploits the acquired data.

IV-C2 Initialization for L>1𝐿1L>1

The following lemma states that in this case, no information is lost by choosing $L=L_{BW}$ compared to taking maximal overlap $L=1$. Moreover, it suggests to upsample the measurement vector by expansion and low-pass interpolation. Our technique resembles standard upsampling arguments in digital signal processing (DSP) (see for instance Section 4.6 of ).

and $\mathbf{F}_{p}$ is a partial Fourier matrix consisting of the first $N/L$ rows of the DFT matrix $\mathbf{F}$.

While Lemma IV.1 shows that no information is lost using an ideal low-pass window with integer bandwidth $N/L$, in practice we do not use these windows. Instead, we approximate the low-pass interpolation of $\mathbf{F}_{p}^{*}\mathbf{F}_{p}$ as suggested in Lemma IV.1 by a simple smooth interpolation. This leads to better numerical results and reduces the computational complexity. In Section V we show simulations with both linear and cubic interpolations.

Following the upsampling stage, the algorithm proceeds as for $L=1$ by extracting the principal eigenvector (with the appropriate normalization) of an approximation matrix. This initialization is summarized in Algorithm 5.

V Numerical Results

This section is devoted to numerical experiments examining the proposed non-convex algorithms. In all experiments, the underlying signal was drawn from $\mathbf{x}\sim\mathcal{N}\left(0,\mathbf{I}\right)$, where $\mathbf{I}$ is the identity matrix. The measurements $\mathbf{Z}$ (II.1) were contaminated with either i.i.d. additive Gaussian noise or Poisson noise. The recovery error is computed by $\frac{d\left(\mathbf{x},\hat{\mathbf{x}}\right)}{\left\|\mathbf{x}\right\|_{2}}$, where $\hat{\mathbf{x}}$ is the estimated signal and the distance function $d\left(\cdot,\cdot\right)$ is defined in Definition II.1. We optimize both the empirical risk loss function (IV.1) and the non-convex phasecut (NCPC) objective function by a trust-region algorithm using the Manopt toolbox .

The first experiment examines the estimation quality of the initialization method described in Algorithm 5. Figure V.1 presents the initialization error as a function of the window’s length. We considered a Gaussian window defined by $\mathbf{g}[n]=e^{\frac{-n^{2}}{2\sigma^{2}}}$ and cubic and linear interpolations. For $n>3\sigma$, we set the entries of the window to be zero so that $W=3\sigma$. The results demonstrate the effectiveness of the smooth interpolation technique. For low values of $L$, it seems that the two interpolations achieve similar performance. For larger $L$, namely, fewer measurements, cubic interpolation outperforms linear interpolation. In the following experiments we use cubic interpolation.

The next experiment aims to estimate the basin of attraction of the loss function (IV.1) or (IV.2). That is to say, the area in which a local optimization method will converge to a global minimum. To do that, we set the initialization vector to be $\mathbf{x}_{0}=\mathbf{x}+\mathbf{z}$, where $\mathbf{x}\sim\mathcal{N}\left(0,\mathbf{I}\right)$ is the underlying signal. The perturbation vector $\mathbf{z}$ takes on the values $\pm\sigma$ (with random signs) for some $\sigma>0$ so that $d(\mathbf{x}_{0},\mathbf{x})\leq\sqrt{N}\sigma$. Then, we applied the trust-region algorithm and checked whether the algorithm converges to $\mathbf{x}$. As can be seen in Figure V.2, the algorithm converges to the global minimum as long as $\sigma\leq 0.3$ for $L=1,2$ (the case of $L=1$ is not presented in the figure) and $\sigma\leq 0.25$ for $L=4$. These experimental results indicate that the actual basin of attraction is larger than our theoretical estimation in Section VI and Theorem VI.2.

Figure V.3 shows a representative example of the performance of Algorithm 2 where we minimized the empirical risk loss function (IV.1) or (IV.2). The experiment was conducted on a signal of length $N=23$ with a rectangular window in a noisy environment of SNR$=20$ dB.

Figure V.4 presents the success rate of the algorithms as a function of the window’s length in a noise-free environment. As can be seen, NCPC achieves the highest success rate, implying that it requires less redundancy in the data. Figure V.5 presents the recovery error for different noise models. Figure V.5a shows the error when the measurements are contaminated with normal noise as a function of the SNR level. The proposed algorithms are compared with Fienup’s method that iterates according to (IV.11). In the low SNR regime, minimizing the ER loss function seems to be better. Figure V.5b shows the error with Poisson noise as a function of $W$. For short windows, NCPC works best. The performance for longer windows is comparable for all algorithms. Figure V.6 presents the same experiments with low-pass data. This reflects a phenomenon that typically occurs in optical applications in which the fine details of the data are blurred by the measurement process. Estimating a signal from its low-resolution measurements, when the phases are available, has been investigated thoroughly in the last years, see for instance . Accordingly, we assume that we can acquire the data $\mathbf{Z}[m,k]$ for all $m$ but only for $k=-K_{\max},\dots,K_{\max}$ for some cut-off frequency $K_{\max}$. Particularly, in Figure V.6 we consider $N=53$ and $K_{\max}=18$ (i.e., $70\%$ of the spectral content) for the two proposed algorithms. In this case, if the SNR is not too bad, then NCPC works significantly better than ER in both cases. As in Figure V.5a, in the low SNR regime, minimizing the ER loss function achieves superior performance for Gaussian noise.

VI Theory

This section presents the theoretical contribution of this work, focusing on the case of maximum overlap between adjacent windows $(L=1)$. As explained and demonstrated numerically, the non-convex approaches also tend to work well for $L>1$ and when the high-frequencies of the data are suppressed.

In our first theoretical result, Theorem VI.1, we analyze the initialization algorithm presented in Algorithm 1 and estimate the distance between the initialization vector and the ground truth. Next, we study the geometry of the loss function (IV.2), which controls the behavior of our ER minimization algorithm. To this end, suppose we minimize the ER loss function (IV.2) using gradient descent followed by a thresholding step that can be used if the signal is bounded. This scheme is presented in Algorithm 6. In Theorem VI.2 we establish the existence of a basin of attraction of size $\frac{1}{8\sqrt{N}W^{2}}$ around the global minimum for signals with unit modulus entries. In the basin of attraction, a gradient algorithm is guaranteed to converge to a global minimum at a geometric rate. This result is true for any gradient scheme with a thresholding step as in Algorithm 6. We stress that the theoretical contribution of this result is limited. As presented in Corollary VI.3, the estimated basin of attraction is small so that theoretically Algorithm 6 converges in the same area in which the problem has a closed linear LS solution. To the best of our knowledge, this is the first result quantifying the size of the basin of attraction of a gradient algorithm in a deterministic phase retrieval setup. This is in contrast to the basin of attraction of random phase retrieval setups which is quite well–understood.

A crucial condition for the success of gradient algorithms is that its initialization is sufficiently close to the global minimum. The following result quantifies the estimation error of the proposed initialization presented in Algorithm 1 for bounded signals and $L=1$. The error reduces to zero as $W$ approaches $\frac{N+1}{2}$. The case of $L>1$ is discussed briefly in Section IV. The result is stated for a normalized signal. The norm of the signal can be estimated easily from the main diagonal of $\mathbf{x}\mathbf{x}^{*}$ as explained in Section III.

Suppose that $L=1$, $\|\mathbf{x}\|_{2}=1$, $\mathbf{g}$ is an admissible window of length $W\geq 2$ and that $\|\mathbf{x}\|_{\infty}\leq\sqrt{\frac{B}{N}}$ for some $0<B\leq\frac{N}{2\left(N-2W+1\right)}$. Then under the measurement model of (II.1), the initialization vector given in Algorithm 1 satisfies

The properties of the gradient algorithm minimizing the ER rely on the geometry of the loss function (IV.2) near the global minimum. The following result quantifies the size of the basin of attraction of the loss function (IV.2), namely, the area in which a gradient algorithm is guaranteed to converge to a global minimum at a geometric rate. As demonstrated in Figure V.2, in practice the basin of attraction is quite large for a broad family of signals. The proof relies on a geometric analysis of the loss function as presented in Lemmas VII.3 and VII.4.

where $\alpha\geq\frac{4N}{W}$ and $\beta\geq 256N^{2}W^{3}$.

Combining Theorems VI.1 and VI.2 leads to the following corollary:

Then, under the measurement model of (II.1), Algorithm 6, initialized by Algorithm 1, with thresholding parameter $B=\frac{1}{\sqrt{N}}$ and step size $0<\mu\leq 2/\beta$ achieves the following geometric convergence:

where $\alpha\geq\frac{4N}{W}$ and $\beta\geq 256N^{2}W^{3}$.

We mention that the result of Corollary VI.3 is good merely for long windows. However, in practice we observe that the algorithm works well also for short windows. As we discuss in Section VIII, bridging this theoretical gap is an important direction for future research.

VII Proofs

The same bound holds for $\left\|\mathbf{E}\right\|_{1}:=\max_{j}\sum_{i}\left|\mathbf{E}[i,j]\right|$ and therefore by Hölder’s inequality we get

In order to complete the proof, we still need to show that if $\left\|\mathbf{X}-\mathbf{X}_{0}\right\|_{2}$ is small, then $d\left(\mathbf{x},\mathbf{x}_{0}\right)$ is small as well, where $\mathbf{{x}}_{0}$ is the principal eigenvector of $\mathbf{X}_{0}$ with appropriate normalization. To show that, we follow the outline of Section 7.8 in . Observe that as $\mathbf{G}_{0}$ is invertible by assumption, the norm of $\mathbf{x}$ is known by

Accordingly, we assume hereinafter without loss of generality that $\mathbf{x}$ and $\mathbf{x}_{0}$ have unit norm. Let $\lambda_{0}$ be the top eigenvalue of $\mathbf{X}_{0}$, associated with $\mathbf{{x}}_{0}$. We observe that

Furthermore, as $\left\|\mathbf{x}\right\|_{2}=1$ we also have

Combining the last two inequalities we get

where the term in the square root is positive by assumption.

VII-B Proof of Theorem VI.2

We say that a function $f$ satisfies the regularity condition in $\mathcal{E}$ if for all vectors $\mathbf{z}\in\mathcal{E}$ we have

for some positive constants $\alpha,\beta$.

The following lemma states that if the regularity condition is met, then the gradient step converges to a global minimum at a geometric rate.

Assume that $f$ satisfies the regularity condition for all $\mathbf{z}\in\mathcal{E}$. Consider the following update rule

In order to show that the regularity condition of Definition VII.1 is met, we present two lemmas for signals with unit modulus entries. The first result shows that the gradient of the loss function (IV.2), given explicitly in (IV.3), is bounded near its global minimum. This implies that the loss function is smooth. We consider here only the case of a rectangular window $\mathbf{g}$ of length $W$. The extension to non-vanishing windows of length $W$ is straightforward (see remark in Appendix -F):

The second lemma shows that the inner product between the gradient and the vector $\mathbf{z}-\mathbf{x}e^{j\phi(\mathbf{z})}$ is positive if $d\left(\mathbf{x},\mathbf{z}\right)\leq\frac{1}{8\sqrt{N}W^{2}}$. This result implies that $-\nabla f(\mathbf{z})$ points approximately towards $\mathbf{x}$. As in Lemma VII.3, we consider for simplicity rectangular windows of length $W$. Yet, the analysis can be extended to non-vanishing windows of length $W$. In this case, the bounds are dependent on the dynamic range of $\mathbf{g}$ (for details, see remark in Appendix Remark).

where $\nabla f(\mathbf{z})$ is given in (IV.3).

We notice that the thresholding stage of Algorithm 6 cannot increase the error as the signal is assumed to be bounded. The proof of Theorem VI.2 is then completed by directly leveraging lemmas VII.3 and VII.4 and seeing that Definition VII.1 holds in our case with constants $\alpha\geq\frac{4N}{W}$ and $\beta\geq 256N^{2}W^{3}$.

VII-C Proof of Corollary VI.3

As $N$ is a prime number, $\mathbf{g}$ is an admissible window of length $W$ (see Lemma III.3). According to Theorem VI.2, we merely need to show that the initialization point is within the basin of attraction, namely, $d\left(\mathbf{x},\mathbf{x}_{0}\right)\leq\frac{1}{8\sqrt{N}W^{2}}$. From Lemma VI.1, we know that the initialization obeys

Using the fact that $a\leq\sqrt{a}$ for all $0\leq a\leq 1$ and some standard algebraic calculations, we conclude that the initialization of Algorithm 1 is within the basin of attraction as long as

VIII Discussion

This paper explores practical, efficient, non-convex phase retrieval algorithms with some deterministic theoretical guarantees. Particularly, we propose two local optimization methods based on minimizing the ER loss function and optimizing on the manifold of phases. The latter is a new phase retrieval algorithm that takes into account the special geometry of the phase retrieval problem.

Since the optimization problems are non-convex, we also propose an initialization method. The method is based on the insight that, for sufficiently long windows, the signal can be recovered as the solution of a linear LS problem. While this may not be true for shorter windows, we use the LS solution to construct a special matrix and initialize the local optimization algorithms with the principal eigenvector of this matrix. Similar initialization approaches were suggested recently for phase retrieval problems. However, they are mainly focused on random setups and based on probabilistic considerations. For $L=1$, we estimate the distance between the initialization point and the ground truth. The case of $L>1$ raises some interesting questions. As a heuristic, we suggested to smoothly interpolate the missing entries. This practice works quite well since the window acts as an averaging operator. Clearly, the interpolation method depends on the window shape. A main challenge for future research is analyzing the setting of $L>1$.

For signals with unit modulus entries, we prove in Theorem VI.2 that the ER loss function has a basin of attraction. We show numerically that the actual basin of attraction is larger than the theoretical bound and exists for a broader family of signals. The gap between the actual size of the basin of attraction and the theoretical result is the bottleneck that prevents a full theoretical understanding of the proposed algorithms. Specifically, improving Lemma VII.4 will lead directly to tighter estimation of the size of the basin of attraction. Ideally, this would lead to the conclusion that the proposed initial guess lies in the basin.

acknowledgment

We would like to thank Mahdi Soltanolkotabi, Iréne Waldspurger, Pavel Sidorenko and Laura Waller for their remarks on an initial draft of this paper.

References

-A Proof of Proposition III.4

-B Proof of Proposition III.5

Then, $\mathbf{x}$ is a principal eigenvectors of $\mathbf{X}_{0}$ (up to global phase).

Based on the special structure of $\mathbf{X_{0}}$, the following calculation shows that $\mathbf{x}$ is an eigenvector of $\mathbf{X}_{0}$ with $\frac{2}{N}$ as the associated eigenvalue:

We still need to show that $\mathbf{x}$ is a principal eigenvector of $\mathbf{X}_{0}$. Since each column and row of $\mathbf{X}_{0}$ is composed of two non-zero values, it is evident that

-C Proof of Proposition III.6

-D Proof of the equality between the loss functions (IV.2) and (IV.1)

Let $\mathbf{U}$ be a unitary matrix. Since unitary matrices do not change the length of a vector, we have

By choosing $\mathbf{U}$ to be the DFT matrix and normalize, we get exactly the loss function in (IV.2).

-E Proof of Lemma IV.1

where $\mathbf{F}_{L}$ consists of the $\left\{jL\thinspace:\thinspace j=0,\dots,N/L-1\right\}$ columns of $\mathbf{F}$ (notice the difference between $\mathbf{F}_{L}$ and $\mathbf{F}_{p}$). We aim at showing that expanding and interpolating $\mathbf{y}_{L}$ as explained in Lemma IV.1 results in $\mathbf{y}$. Direct computation shows that the expansion stage as described in (IV.12) is equivalent to multiplying both sides by $\mathbf{F^{*}F}_{L}$:

Let us denote $\mathbf{T:=F}_{L}\mathbf{F}_{L}^{*}$, which is a Toeplitz matrix with $L$ on the $\frac{jN}{L}$ diagonals for $j=0,\dots,N/L-1$ and zero otherwise. Because of the structure of $\mathbf{\Sigma}$ we can then write

Comparing (E.2) with (E.1) completes the proof.

-F Proof of Lemma VII.3

For convenience, let us denote $d\left(\mathbf{x},\mathbf{z}\right)=\frac{\varepsilon}{\sqrt{N}}$ for some $\varepsilon\leq 1$ and therefore $|\mathbf{z}[n]|\geq\frac{1-\varepsilon}{\sqrt{N}}$ for all $n$. Accordingly, for any $(n,k)$,

Since $\mathbf{x}[n]$ and $\mathbf{z}[n]|$ have the same sign pattern, we have

In case of a non-vanishing window of length $W$, one can easily bound the gradient using the same technique, while taking into account $\max_{n}\left|\mathbf{g}[n]\right|$ in the inequalities.

-G Proof of Lemma VII.4

Clearly, if $\mathbf{z}=\mathbf{x}e^{j\phi(\mathbf{z})}$ then $\left\langle\nabla f(\mathbf{z}),\mathbf{z}-\mathbf{x}e^{j\phi(\mathbf{z})}\right\rangle=0$. Otherwise, the first term of (G.1) is strictly positive. Hence, in order to achieve a lower bound on (G.1), we first derive an upper bound on the second term and then bound the first term from below.

Next, we aim to bound the first term of (G.1) from below as follows:

where the last inequality is true since $\mathbf{x}^{2}[n]\geq\mathbf{z}^{2}[n]$ and for any positive (or negative) sequence $\left\{a_{i}\right\}$ we have $\left(\sum_{i}a_{i}\right)^{2}\geq\sum_{i}a_{i}^{2}$. Furthermore, since $|\mathbf{z}[n]|=\frac{1-\varepsilon_{n}}{\sqrt{N}}$ we have

Therefore, since $\varepsilon_{n}\leq 1$ for all $n$ we conclude that

Plugging (G.2) and (G.4) into (G.1) yields

where the last inequality holds for $d(\mathbf{x},\mathbf{z})\leq\frac{1}{8\sqrt{N}W^{2}}.$

Observe that the analysis for non-vanishing windows of length $W$ requires only a small modification. In this case, one should use the maximal and the minimal values of the window in the above inequalities. For instance, one would need to take $\mbox{g}_{\mbox{min}}:=\min_{n=0,\dots,W-1}\left|\mathbf{g}[n]\right|$ into account in (G.3).