PhaseMax: Convex Phase Retrieval via Basis Pursuit

I Introduction

Here, $\langle\mathbf{x},\hat{}\mathbf{x}\rangle_{\Re}$ denotes the real-part of the inner product between the vectors $\mathbf{x}$ and $\hat{}\mathbf{x}$. The main idea behind PhaseMax is to find the vector $\mathbf{x}$ that is most aligned with the approximation vector $\hat{}\mathbf{x}$ and satisfies a convex relaxation of the measurement constraints in (1).

Our main goal is to develop sharp lower bounds on the probability with which PhaseMax succeeds in recovering the true signal $\mathbf{x}^{0}$, up to an arbitrary phase ambiguity that does not affect the measurement constraints in (1). By assuming noiseless measurements, one of our main results is as follows.

be the angle between the true vector $\mathbf{x}^{0}$ and the approximation $\hat{}\mathbf{x}$, and define the constant

In words, if $m>4n/\!\alpha$ and $\alpha>0$, then PhaseMax will succeed with non-zero probability. Furthermore, for a fixed signal dimension $n$ and an arbitrary approximation vector $\hat{}\mathbf{x}$ that satisfies $\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x})<\frac{\pi}{2}$, i.e., one that is not orthogonal to the vector $\mathbf{x}^{0}$, we can make the success probability of PhaseMax arbitrarily close to one by increasing the number of measurements $m$. As we shall see, our recovery guarantees are sharp and accurately predict the performance of PhaseMax in practice.

We emphasize that the convex formulation (PM) has been studied by several other authors, including Bahmani and Romberg , whose work appeared shortly before our own. Related work will be discussed in detail in Section I-C.

It is quite intriguing that the following Basis Pursuit problem

with $\operatorname*{phase}(z)=z/|z|$ for $z\neq 0$ and $\operatorname*{phase}(0)=1.$ This observation not only reveals a fundamental connection between phase retrieval and sparse signal recovery, but also implies that Basis Pursuit solvers can be used to recover the signal from the phase-less measurements in (1).

I-B A Brief History of Phase Retrieval

Phase retrieval is a well-studied problem with a long history and enjoys widespread use in applications such as X-ray crystallography , microscopy , imaging , and many more . Early algorithms, such as the Gerchberg-Saxton or Fienup algorithms, rely on alternating projection to recover complex-valued signals from magnitude-only measurements. The papers sparked new interest in the phase retrieval problem by showing that it can be relaxed to a semidefinite program. Prominent convex relaxations include PhaseLift and PhaseCut . These methods come with recovery guarantees, but require the problem to be lifted to a higher dimensional space, which prevents their use for large-scale problems. Several authors recently addressed this problem by presenting efficient methods that can solve large-scale semidefinite programs without the complexity of representing the full-scale (or lifted) matrix of unknowns. These approaches include gauge duality methods and sketching methods .

More recently, a number of non-convex algorithms have been proposed (see e.g., ) that directly operate in the original signal dimension and exhibit excellent empirical performance. The algorithms in come with recovery guarantees that mainly rely on accurate initializers, such as the (truncated) spectral initializer , the Null initializer , the orthogonality-promoting method , or more recent initializers that guarantee optimal performance (see Section VI for additional details). These initializers enable non-convex phase retrieval algorithms to succeed, given a sufficiently large number of measurements; see for more details on the geometry of such non-convex problems.

I-C Related Work on Non-Lifting Convex Methods

Because of the extreme computational burden of traditional convex relaxations (which require “lifting” to a higher dimension), a number of authors have recently been interested in non-lifting convex relaxations for phase retrieval. Shortly before the appearance of this work, the formulation (PM) was studied by Bahmani and Romberg . Using methods from machine learning theory, the authors derived bounds on the recovery of signals with and without noise. The results in are stronger than those presented here in that they are uniform with respect to the approximation vector $\hat{}\mathbf{x},$ but weaker in the sense that they require significantly more measurements.

Several authors have studied PhaseMax after the initial appearance of this work. An alternative proof of accurate signal recovery was derived using measure concentration bounds in . Compressed sensing methods for sparse signals , and corruption-robust methods for noisy signals have also been presented. The closely related non-lifting relaxation BranchHull has also been proposed for recovering signals from entry-wise product measurements .

Finally, we note that recent works have proved tight asymptotic bounds for PhaseMax in the asymptotic limit, i.e., where $\beta=m/n$ is constant and $m\to\infty$. The authors of derive an exact asymptotic bound of the performance of PhaseMax, and a related non-rigorous analysis was given in . It was also shown that better signal recovery guarantees can be obtained by iteratively applying PhaseMax. The resulting method, called PhaseLamp, was analyzed in .

I-D Contributions and Paper Outline

In contrast to algorithms relying on semidefinite relaxation or non-convex problem formulations, we propose PhaseMax, a novel, convex method for phase retrieval that directly operates in the original signal dimension. In Section II, we establish a deterministic condition that guarantees uniqueness of the solution to the (PM) problem. We borrow methods from geometric probability to derive sharp lower bounds on the success probability for real- and complex-valued systems in Section III. Section V generalizes our results to a broader range of random measurement ensembles and to systems with measurement noise. We show in Section VI that randomly chosen approximation vectors are sufficient to ensure faithful recovery, given a sufficiently large number of measurements. We numerically demonstrate the sharpness of our recovery guarantees and showcase the practical limits of PhaseMax in Section VII. We conclude in Section VIII.

I-E Notation

II Uniqueness Condition

The measurement constraints in (1) do not uniquely define a vector. If $\mathbf{x}$ is a vector that satisfies (1), then any vector $\mathbf{x}^{\prime}=e^{j\phi}\mathbf{x}$ for $\phi\in[0,2\pi)$ also satisfies the constraints. In contrast, if $\mathbf{x}$ is a solution to (PM), then $e^{j\phi}\mathbf{x}$ with $\phi\neq 0$ will not be another solution. In fact, consider any vector $\mathbf{x}$ in the feasible set of (PM) with $\langle\mathbf{x},\hat{}\mathbf{x}\rangle_{\Im}\neq 0$. By choosing $\omega=\operatorname*{phase}(\langle\mathbf{x},\hat{}\mathbf{x}\rangle)$, we have

which implies that given such a vector $\mathbf{x},$ one can always increase the objective function of (PM) simply by aligning $\mathbf{x}$ with the approximation $\hat{}\mathbf{x}$ (i.e., modifying its phase so that $\langle\mathbf{x},\hat{}\mathbf{x}\rangle$ is real valued). The following definition makes this observation rigorous.

A vector $\mathbf{x}$ is said to be aligned with another vector $\hat{}\mathbf{x}$, if the inner product $\langle\mathbf{x},\hat{}\mathbf{x}\rangle$ is real-valued and non-negative.

From all the vectors that satisfy the measurement constraints in (1), there is only one that is a candidate solution to the convex problem (PM), which is also the solution that is aligned with $\hat{}\mathbf{x}.$ For this reason, we adopt the following important convention throughout the rest of this paper.

The true vector $\mathbf{x}^{0}$ denotes a solution to (1) that is aligned with the approximation vector $\hat{}\mathbf{x}.$

There is an interesting relation between the convex formulation of PhaseMax and the semidefinite relaxation method PhaseLift . Recall that the set of solutions to any convex problem is always convex. However, the solution set of the measurement constraints (1) is invariant under phase rotations, and thus non-convex (provided it is non-zero). It is therefore impossible to design a convex problem that yields this set of solutions. PhaseMax and PhaseLift differ in how they remove the phase ambiguity from the problem to enable a convex formulation. Rather than trying to identify the true vector $\mathbf{x}^{0},$ PhaseLift reformulates the problem in terms of the quantity $\mathbf{x}^{0}(\mathbf{x}^{0})^{H},$ which is unaffected by phase rotations in $\mathbf{x}^{0}.$ Hence, PhaseLift removes the rotation symmetry from the solution set, yielding a problem with a convex set of solutions. PhaseMax does something much simpler: it pins down the phase of the solution to an arbitrary quantity, thus removing the phase ambiguity and restoring convexity to the solution set. This arbitrary phase choice is made when selecting the phase of the approximation $\hat{}\mathbf{x}.$

We are now ready to state a deterministic condition under which PhaseMax succeeds in recovering the true vector $\mathbf{x}^{0}.$ The result applies to the noiseless case, i.e., $\eta_{i}=0$, $i=1,2,\ldots,m$. In this case, all inequality constraints in (PM) are active at $\mathbf{x}^{0}.$ The noisy case will be discussed in Section V-B.

The true vector $\mathbf{x}^{0}$ is the unique maximizer of (PM) if, for any unit vector $\boldsymbol{\delta}\in\mathcal{H}^{n}$ that is aligned with the approximation $\hat{}\mathbf{x},$

Suppose the conditions of this theorem hold, and consider some candidate solution $\mathbf{x}^{\prime}$ in the feasible set for (PM) with $\langle\mathbf{x}^{\prime},\hat{}\mathbf{x}\rangle\geq\langle\mathbf{x}^{0},\hat{}\mathbf{x}\rangle.$ Without loss of generality, we assume $\mathbf{x}^{\prime}$ to be aligned with $\hat{}\mathbf{x}.$ Then, the vector $\mathbf{\Delta}=\mathbf{x}^{\prime}-\mathbf{x}^{0}$ is also aligned with $\hat{}\mathbf{x}$, and satisfies

Since $\mathbf{x}^{\prime}$ is a feasible solution for (PM), we have

But $|\mathbf{a}_{i}^{T}\mathbf{x}^{0}|^{2}=b_{i}^{2},$ and so

Now, if $\|\mathbf{\Delta}\|_{2}>0,$ then the unit-length vector $\boldsymbol{\delta}=\mathbf{\Delta}/\|\mathbf{\Delta}\|_{2}$ satisfies $[\langle\mathbf{a}_{i},\mathbf{x}^{0}\rangle^{*}\langle\mathbf{a}_{i},\boldsymbol{\delta}\rangle]_{\Re}\leq 0$ for all $i,$ which contradicts the hypothesis of the theorem. It follows that $\|\mathbf{\Delta}\|_{2}=0$ and $\mathbf{x}^{\prime}=\mathbf{x}^{0}.$ ∎

Theorem 2 has an intuitive geometrical interpretation. If $\mathbf{x}^{0}$ is an optimal point and $\boldsymbol{\delta}$ is an ascent direction, then one cannot move in the direction of $\boldsymbol{\delta}$ starting at $\mathbf{x}^{0}$ without leaving the feasible set. This condition is met if there is an $\mathbf{a}_{i}$ such that $\mathbf{x}^{0}$ and $\boldsymbol{\delta}$ both lie on the same side of the plane through the origin orthogonal to the measurement vector $\mathbf{a}_{i}.$

III Preliminaries: classical sphere covering problems and geometric probability

To derive sharp conditions on the success probability of PhaseMax, we require a set of tools from geometric probability. Many classical problems in geometric probability involve calculating the likelihood of a sphere being covered by random “caps,” or semi-spheres, which we define below.

Consider the set $\mathcal{S}^{n-1}_{\mathcal{H}}=\{\mathbf{x}\in\mathcal{H}^{n}\,|\,\|\mathbf{x}\|_{2}=1\},$ the unit sphere embedded in $\mathcal{H}^{n}.$ Given a vector $\mathbf{a}\in\mathcal{H}^{n},$ the cap centered at $\mathbf{a}$ with central angle $\theta$ is defined as

This cap contains all vectors that form an angle with $\mathbf{a}$ of less than $\theta$ radians. When $\theta=\pi/2,$ we have a semisphere centered at $\mathbf{a},$ which is simply denoted by

We say that a collection of caps covers the entire sphere if the sphere is contained in the union of the caps. Before we can say anything useful about when a collection of caps covers the sphere, we will need the following classical result, which is often attributed to Schläfli . Proofs that use simple induction methods can be found in . For completeness, we briefly include a short proof in Appendix A.

Classical results in geometric probability study the likelihood of a sphere being covered by random caps with centers chosen independently and uniformly from the sphere’s surface. For our purposes, we need to study the more specific case in which caps are only chosen from a subset of the sphere. While calculating this probability is hard in general, it is quite simple when the set obeys the following symmetry condition.

We say that the set $\mathcal{A}$ is symmetric if, for all $\mathbf{x}\in\mathcal{A},$ we also have $-\mathbf{x}\in\mathcal{A}.$

We are now ready to prove a general result that states when the sphere is covered by random caps.

This is the probability of turning up $n$ or more heads when flipping $m_{\mathcal{A}}-1$ fair coins.

Consider the following two-step process for constructing the set $\{\mathbf{a}_{i}\}.$ First, we sample $m_{\mathcal{A}}$ vectors $\{\mathbf{a}^{\prime}_{i}\}$ independently and uniformly from $\mathcal{A}.$ Note that, because $\mathcal{A}$ has positive measure, it holds with probability 1 that any subset of $n$ vectors will be linearly independent.

Second, we define $\mathbf{a}_{i}=c_{i}\mathbf{a}^{\prime}_{i},$ where $\{c_{i}\}$ are i.i.d. Bernoulli variables that take value $+1$ or $-1$ with probability $\nicefrac{{1}}{{2}}.$ We can think of this second step as randomly “flipping” a subset of uniform random vectors. Since $\mathcal{A}$ is symmetric and $\{\mathbf{a}^{\prime}_{i}\}$ is sampled independently and uniformly, the random vectors $\{\mathbf{a}_{i}\}$ also have an independent and uniform distribution over $\mathcal{A}.$ This construction may seem superfluous since both $\{\mathbf{a}_{i}\}$ and $\{\mathbf{a}^{\prime}_{i}\}$ have the same distribution, but we will see below that this becomes useful.

Several papers have studied the probability of covering the sphere using points independently and uniformly chosen over the entire sphere. The only aspect that is unusual about Lemma 2 is the observation that this probability remains the same if we restrict our choices to the set $\mathcal{A},$ provided $\mathcal{A}$ is symmetric. We note that this result was observed by Gilbert in the case $n=3,$ and we generalize it to any $n>1$ using a similar argument.

We now present a somewhat more complicated covering theorem. The next result considers the case where the measurement vectors are drawn only from a semisphere. We consider the question of whether these vectors cover enough area to contain not only their home semisphere, but another nearby semisphere as well.

Because the mapping $\mathbf{a}_{i}\to\mathbf{a}^{\prime}_{i}$ is onto and (piecewise) isometric, $\{\mathbf{a}^{\prime}_{i}\}$ will be uniformly distributed over the half sphere $\{\mathbf{a}\,|\,\langle\mathbf{a},\mathbf{x}^{0}\rangle>0\}$ whenever $\{\mathbf{a}_{i}\}$ are independently and uniformly distributed over the entire sphere.

Consider the “hourglass” shaped, symmetric set

As noted in Lemma 2, the expression on the right is the chance of turning up more than $n$ heads when flipping $m_{\mathcal{A}}$ fair coins.

The probability $p_{\text{cover}}(m,n;\mathbf{x},\mathbf{y})$ is then given by

The expression on the right hand side is the probability of getting more than $n$ heads when one fair coin is flipped for every measurement $\mathbf{a}_{i}$ that lies in $\mathcal{A}$.

Using Hoeffding’s inequality, we obtain the following lower bound

Lemma 3 contains most of the machinery needed for the proofs that follow. In the sequel, we prove a number of exact reconstruction theorems for (PM). Most of the results rely on short arguments followed by the invocation of Lemma 3.

We finally state a result that bounds from below the probability of covering the sphere with caps of small central angle. The following Lemma is a direct corollary of the results of Burgisser, Cucker, and Lotz in . A derivation that uses their results is given in Appendix B.

IV Recovery Guarantees

Using the uniqueness condition provided by Theorem 2 and the tools derived in Section III, we now develop sharp lower bounds on the success probability of PhaseMax for noiseless real- and complex-valued systems. The noisy case will be discussed in Section V-B.

The true vector $\mathbf{x}^{0}$ will be the unique maximizer of (PM) if

Using this observation, we can develop the following lower bound on the success probability of PhaseMax for real-valued systems.

where $\alpha=1-\frac{2}{\pi}\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x})$ and $m>2n/\alpha.$

IV-B The Complex Case

We now prove Theorem 1 given in Section I, which characterizes the success probability of PhaseMax for phase retrieval in complex-valued systems. For clarity, we restate our result in shorter form.

where $\alpha=1-\frac{2}{\pi}\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x})$ and $m>4n/\alpha.$

By Theorem 2, the true signal $\mathbf{x}^{0}$ will be the unique maximizer of (PM) if

Let us bound the probability of this event. Consider the set $\mathcal{B}=\{\boldsymbol{\delta}\,|\,\langle\boldsymbol{\delta},\mathbf{x}^{0}\rangle_{\Im}=0\}.$ We now claim that (11) holds whenever

Theorems 1 and 3 guarantee exact recovery for a sufficiently large number of measurements $m$ provided that $\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x})<\frac{\pi}{2}.$ In the case $\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x})>\frac{\pi}{2},$ our theorems guarantee convergence to $-\mathbf{x}^{0}$ (which is also a valid solution) for sufficiently large $m$. Our theorems only fail for large $m$ if $\arccos(\mathbf{x}^{0},\hat{}\mathbf{x})=\pi/2,$ which happens with probability zero when the approximation vector $\hat{}\mathbf{x}$ is generated at random. See Section VI for more details.

V Generalizations

Our theory thus far addressed the idealistic case in which the measurement vectors are independently and uniformly sampled from a unit sphere and for noiseless measurements. We now extend our results to more general random measurement ensembles and to noisy measurements.

The theorems of Section IV require the measurement vectors $\{\mathbf{a}_{i}\}$ to be drawn independently and uniformly from the surface of the unit sphere. This condition can easily be generalized to other sampling ensembles. In particular, our results still hold for all rotationally symmetric distributions. A distribution $D$ is rotationally symmetric if the distribution of $\mathbf{a}/\|\mathbf{a}\|_{2}$ is uniform over the sphere when $\mathbf{a}\sim D.$ For such a distribution, one can make the change of variables $\mathbf{a}\leftarrow\mathbf{a}/\|\mathbf{a}\|_{2},$ and then apply Theorems 1 and 3 to the resulting problem. Note that this change of variables does not change the feasible set for (PM), and thus does not change the solution. Consequently, the same recovery guarantees apply to the original problem without explicitly implementing this change of variables. We thus have the following simple corollary.

The results of Theorem 1 and Theorem 3 still hold if the samples $\{\mathbf{a}_{i}\}$ are drawn from a multivariate Gaussian distribution with independent and identically distributed (i.i.d.) entries.

A multivariate Gaussian distribution with i.i.d. entries is rotationally symmetric, and thus the change of variables $\mathbf{a}\leftarrow\mathbf{a}/\|\mathbf{a}\|_{2}$ yields an equivalent problem with measurements sampled uniformly from the unit sphere. ∎

What happens when the distribution is not spherically symmetric? In this case, we can still guarantee recovery, but we require a larger number of measurements. The following result is, analogous to Theorem 1, for the noiseless complex case.

The ratio of the non-uniform density (15) to the uniform density (14) is

V-B Noisy Measurements

We now analyze the sensitivity of PhaseLift to the measurement noise $\{\eta_{i}\}.$ For brevity, we focus only on the case of complex-valued signals. To analyze the impact of noise, we re-write the problem (PM) in the following equivalent form:

Here, $\hat{b}_{i}^{2}=|\langle\mathbf{a}_{i},\mathbf{x}^{0}\rangle|^{2}$ is the (unknown) true magnitude measurement and $b_{i}^{2}=\hat{b}_{i}^{2}+\eta_{i}.$ We are interested in bounding the impact that these measurement errors have on the solution to (PM). Note that the severity of a noise perturbation of size $\eta_{i}$ depends on the (arbitrary) magnitude of the measurement vector $\mathbf{a}_{i}.$ For this reason, we assume the vectors $\{\mathbf{a}_{i}\}$ have unit norm throughout this section.

We will begin by proving results only for the case of non-negative noise. We will then generalize our analysis to the case of arbitrary bounded noise. The following result gives a geometric characterization of the reconstruction error.

Our goal is to put a bound on the magnitude of the recovery error $\mathbf{\Delta}.$ We start by reformulating the constraints in (24) to get

Since $\eta_{i}$ is non-negative, we have

Now suppose that $\|\mathbf{\Delta}\|_{2}>\varepsilon.$ From (26) we have

Using this result, we can bound the reconstruction error in the noisy case. For brevity, we present results only for the complex-valued case.

We now claim that the conditions of Theorem 5 hold whenever

We now consider the case of noise that takes on both positive and negative values. In this case, we bound the error by converting the problem into an equivalent problem with non-negative noise, and then apply Theorem 6.

Suppose the vectors $\{\mathbf{a}_{i}\}$ in (19) are normalized to have unit length, and that $\eta_{i}>-\hat{b}_{i}^{2}$ for all $i.$This assumption is required so that $b_{i}^{2}=\hat{b}_{i}^{2}+\eta_{i}>0$. Define the following measures of the noise

Choose some error bound $\varepsilon>r/2,$ and define the angle $\phi=\arccos(r/2\varepsilon)-\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x}).$ Then, we have the bound

Consider the “shrunk” version of problem (19)

for some real-valued “shrink factor” $s>0.$ Clearly, if $\mathbf{x}^{0}$ is aligned with $\hat{}\mathbf{x}$ and satisfies $|\langle\mathbf{a}_{i},\mathbf{x}^{0}\rangle|=b_{i}$ for all $i,$ then $s\mathbf{x}^{0}$ is aligned with $\hat{}\mathbf{x}$ and satisfies $|\langle\mathbf{a}_{i},s\mathbf{x}^{0}\rangle|=sb_{i}.$ We can now transform the noisy problem (19) into an equivalent problem with non-negative noise by choosing

We then have $(sb_{i})^{2}+\zeta_{i}=b_{i}^{2}+\eta_{i}^{2},$ and so problem (30) is equivalent to problem (19). However, the noise $\zeta_{i}$ in problem (30) is non-negative, and thus we can apply Theorem 6. This theorem requires the constant $r$ for the shrunken problem, which is now

The solution to the shrunk problem (30) satisfies $\|\mathbf{x}^{\star}-s\mathbf{x}^{0}\|_{2}\leq\epsilon,$ with probability $p_{\text{cover}}(m,2n,\phi),$ where $\phi=\arccos(r_{\text{shrunk}}/2\epsilon)-\operatorname*{angle}(\mathbf{x}^{0},\hat{}\mathbf{x}).$ If this condition is fulfilled, then we have

Theorem 7 requires the noise to be sufficiently small so that the measurements are non-negative, and the PhaseMax formulation is feasible. A natural extension that avoids this caveat is to enforce constraints with a hinge penalty rather than a hard constraint. In addition, it is possible to achieve better noise robustness in this case. This direction was studied in .

VI How to Compute Approximation Vectors?

There exist a variety of algorithms that compute approximation vectorsApproximation vectors are also known as initialization or anchor vectors., such as the (truncated) spectral initializer or corresponding optimized variants , the Null initializer , the orthogonality-promoting method , or least-squares methods . We now show that even randomly generated approximation vectors guarantee the success of PhaseMax with high probability given a sufficiently large number of measurements. We then show that more sophisticated methods guarantee success with high probability if the number of measurements depends linearly on $n.$

Consider the angle $\beta=\operatorname*{angle}(\mathbf{x},\mathbf{y})$ between two random vectors $\mathbf{x},\mathbf{y}\in\mathcal{S}^{{n}-1}_{\mathcal{H}}$ sampled independently and uniformly from the unit sphere. Then, the expected magnitude of the cosine distance between the two random vectors satisfies

We first consider the real case. The quantity $\cos(\beta)=\langle\mathbf{x},\mathbf{y}\rangle/(\|\mathbf{x}\|_{2}\|\mathbf{y}\|_{2})$ is simply the sample correlation between two random vectors, whose distribution function is given by

where $B$ is the beta function and $z\in[-1,+1]$. Hence, the expectation of the magnitude of the inner product is given by

The integral in the numerator was studied in [50, Eq. 31] and evaluates to $\frac{1}{2}B(1,\frac{n-1}{2}).$ Plugging this expression into (33), and using the identity $\Gamma(a)/\Gamma(a/2)=2^{a-1}\Gamma(\frac{a+1}{2})/\sqrt{\pi},$ we get

for the real case. Plugging this bound on the expected value for $\alpha$ into Theorem 3, we see that, for an average randomly-generated approximation vector (one with $\alpha>\sqrt{8/(\pi^{3}n)}$), the probability of exact reconstruction goes to $1$ rapidly as $n$ goes to infinity, provided that the number of measurements satisfies $m>cn^{3/2}$ for any $c>\sqrt{\pi^{3}/2}.$ For complex-valued signals and measurement matrices, this becomes $c>2\sqrt{\pi^{3}}.$

Our results indicate that the use of random approximation vectors requires $O(n^{3/2})$ measurements rather than the $O(n)$ required by other phase retrieval algorithms with other initialization methods, e.g., the ones proposed in (see also Section VII). Hence, it may be more practical for PhaseMax to use approximation vectors obtained from more sophisticated initialization algorithms.

VI-B Spectral Initializers

Spectral initializers were first used in phase retrieval by , and enable the computation of an approximation vector $\hat{}\mathbf{x}$ that exhibits strong theoretical properties . This class of initializers was analyzed in detail in , and developed an optimal spectral initializer for a class of problems including phase retrieval.

Fix some $\delta>1,$ and consider a measurement process with $m=\delta n$ measurement vectors sampled from a Gaussian distribution. Then, the spectral initializers in are known to deliver an initializer $\hat{}\mathbf{x}$ with $\lim_{n\to\infty}\operatorname*{angle}(\hat{}\mathbf{x},\mathbf{x}^{0})>\epsilon$ for some positive $\epsilon.$ Equivalently, there is some accuracy parameter $\alpha_{\star}=\lim_{n\to\infty}1-\frac{2}{\pi}\operatorname*{angle}(\hat{}\mathbf{x},\mathbf{x}^{0})>0$. By combining this result with Theorem 1, we see that spectral initializers enable PhaseMax to succeed for large $n$ with high probability provided that $m>n/\alpha_{*}$ measurements are used for signal recoveryOur results assume the initial vector is independent of the measurement vectors. The total number of measurements needed is thus $m>n/\alpha_{*}+\delta n$ for initialization and recovery combined.; a number of measurements that is linear in $n$. Note that similar non-asymptotic guarantees can be achieved for finite $n$ using other results on spectral initializers (e.g., Prop. 8 in ), although with less sharp bounds.

In Section VII-B, we investigate the numerical value of $\alpha_{*}$ and make quantitative statements about how many measurements are required when using a spectral initializer. As we will see there, PhaseMax requires roughly $m=5.5n$ noiseless Gaussian measurements in theoryThe constants that appear in the bounds were approximated using stochastic numerical methods, and were not obtained exactly using analysis., and somewhat fewer measurements in practice.

VII Discussion

This section briefly compares our theoretical results to that of existing algorithms. We furthermore demonstrate the sharpness of our recovery guarantees and show the practical limits of PhaseMax.

Table I compares our noiseless recovery guarantees in a complex system to that of PhaseLift , truncated Wirtinger flow (TWF) , and truncated amplitude flow (TAW) . We also compare to the recovery guarantee provided for PhaseMax using classical machine learning methods in . Since AltMinPhase requires an online measurement model that differs significantly from the other algorithms considered here, we omit a comparison. We see that PhaseMax requires the same sample complexity (number of required measurements) as compared to PhaseLift, TWF, and TAW, when used together with the truncated spectral initializer . While the constants $c_{0}$, $c_{1}$, and $c_{2}$ in the recovery guarantees for all of the other methods are generally very large, our recovery guarantees contain no unspecified constants, explicitly depend on the approximation factor $\alpha$, and are surprisingly sharp. We next demonstrate the accuracy of our results via numerical simulations.

VII-B Tightness of Recovery Guarantees

We now investigate the tightness of the recovery guarantees of PhaseMax using both experiments and theory. First, we compare the empirical success probability of PhaseMax in a noiseless and complex-valued scenario with measurement vectors taken independently and uniformly from the unit sphere. All experiments were obtained by using the implementations provided in the software library PhasePack . This software library provides efficient implementations of phase retrieval methods using fast adaptive gradient solvers . We declare signal recovery a success whenever the relative reconstruction error satisfies

We compare empirical rates of success to the theoretical lower bound in Theorem 1. Figure 2 shows results for $n=100$ and $n=500$ measurements, where we artificially generate an approximation $\hat{}\mathbf{x}$ for different angles $\beta=\operatorname*{angle}(\hat{}\mathbf{x},\mathbf{x}^{0})$ measured in degrees. Clearly, our theoretical lower bound accurately predicts the real-world performance of PhaseMax. For large $n$ and large $\beta$, the gap between theory and practice becomes extremely tight. We furthermore observe a sharp phase transition between failure and success, with the transition getting progressively sharper for larger dimensions $n$.

Next, we investigate the tightness of the recovery guarantees when an initial vector is chosen using the spectral initialization method . Using analytical formulas, the authors of calculate the accuracy of their proposed spectral initializer for real-valued signal estimation from Gaussian measurements with large $n$. They find that, when $2n$ measurements are used for spectral estimation, the initializer and signal have a squared cosine similarity of over $0.6.$ This corresponds to an accuracy parameter of $\alpha>0.55,$ and PhaseMax can recover the signal exactly with $m=n/\alpha\approx 3.5n$ measurements. This is within a factor of 2 of the information-theoretic lower bound for real-valued signal recovery, which is $2n-2$ measurements.

Note that our analysis of PhaseMax assumes that the measurements used for recovery are statistically independent from those used by the initializer. For this reason, our theory does not allow us to perform signal recovery by “recycling” the measurements used for initialization. While we do not empirically observe any change in behavior of the method when the initialization measurements are used for recovery, the above reconstruction bounds formally require $5.5n$ measurements ($2n$ for initialization, plus $3.5n$ for recovery). The bounds in are uniform with respect to the initializer, and thus enable measurement recycling (although this does not result in tighter bounds for the overall number of measurements because the required constants are larger).

Finally, we would like to mention that asymptotically exact performance bounds for PhaseMax have recently been derived in , and tighter bounds for non-lifting phase retrieval have been shown for the method PhaseLamp, which repeatedly uses PhaseMax within an iterative process.

VII-C Experimental Results

We compare PhaseMax to other algorithms using synthetic and empirical datasets. The implementations used here are publicly available as part of the PhasePack software library . This software library contains scripts for comparing different phase retrieval methods, including scripts for reproducing the experiments shown here.

We briefly compare PhaseMax to a select set of phase retrieval algorithms in terms of relative reconstruction error with Gaussian measurements. We emphasize that this comparison is by no means intended to be exhaustive and serves the sole purpose of demonstrating the efficacy and limits of PhaseMax (see, e.g., for more extensive phase retrieval algorithm comparisons). We compare the Gerchberg-Saxton algorithm , the Fienup algorithm , truncated Wirtinger flow , and PhaseMax—all of these methods use the truncated spectral initializer . We also run simulations using the semidefinite relaxation (SDR)-based method PhaseLift implemented via FASTA ; this is, together with PhaseCut , the only convex alternative to PhaseMax, but lifts the problem to a higher dimension.

Figure 3 reveals that PhaseMax requires larger oversampling ratios $m/n$ to enable faithful signal recovery compared to non-convex phase-retrieval algorithms that operate in the original signal dimension. This is because the truncated spectral initializer requires oversampling ratios of about six or higher to yield sufficiently accurate approximation vectors $\hat{}\mathbf{x}$ that enable PhaseMax to succeed.

VII-C2 Comparisons Using Empirical Data

The PhasePack library contains scripts for comparing phase retrieval algorithms using synthetic measurements as well as publicly available real datasets. The datasets contain measurements obtained using the experimental setup described in , in which a binary mask is imaged through a diffusive medium.

Figure 4 shows reconstructions from measurements obtained using two different test images. Both images were acquired using phaseless measurements from the same measurement operator $\mathbf{A}$. Reconstructions are shown for the Gerchberg-Saxton , Wirtinger flow , and truncated Wirtinger flow methods. We also show results of the two convex methods PhaseLift (which lifts the problem dimension) and PhaseMax (the proposed non-lifting relaxation). For each algorithm, Figure 4 shows the recovered images and reports the relative measurement error $\||\mathbf{A}\mathbf{x}|-\mathbf{b}\|_{2}/\|\mathbf{b}\|_{2}.$ All phase retrieval algorithms were initialized using a variant of the spectral method proposed in .

We find that the Gerchberg-Saxton algorithm outperforms all other considered methods by a small margin (in terms of measurement error) followed by the Wirtinger flow. On this dataset, measurements seem to be quite valuable, and the truncated Wirtinger flow method (with default truncation parameters as described in ) performs slightly worse than the original Wirtinger flow, which uses the full dataset. PhaseMax produces results that are visually comparable to that of Wirtinger flow, but with slightly higher relative measurement error.

VII-D Advantages of PhaseMax

While PhaseMax does not achieve exact reconstruction with the lowest number of measurements, it is convex, operates in the original signal dimension, can be implemented via efficient solvers for Basis Pursuit, and comes with sharp performance guarantees that do not sweep constants under the rug (cf. Figure 2). The convexity of PhaseMax enables a natural extension to sparse phase retrieval or other signal priors (e.g., total variation or bounded infinity norm) that can be formulated with convex functions. Such non-differentiable priors cannot be efficiently minimized using simple gradient descent methods (which form the basis of Wirtinger or amplitude flow, and many other methods), but can potentially be solved using standard convex solvers when combined with the PhaseMax formulation.

VIII Conclusions

We have proposed a novel, convex phase retrieval algorithm, which we call PhaseMax. We have provided accurate bounds on the success probability that depend on the signal dimension, the number of measurements, and the angle between the approximation vector and the true vector. Our analysis covers a broad range of random measurement ensembles and characterizes the impact of general measurement noise on the solution accuracy. We have demonstrated the sharpness of our recovery guarantees and studied the practical limits of PhaseMax via simulations.

There are many avenues for future work. Developing a computationally efficient algorithm for solving PhaseMax (or the related PhaseLamp procedure) accurately and for large problems is a challenging problem. An accurate analysis of more general noise models is left for future work.

Appendix A Proof of Lemma 1

We leave it to the reader to verify that (5) satisfies this recurrence relation and base case.

Appendix B Proof of Lemma 4

In this section, we prove Lemma 4. This Lemma is a direct corollary of the following result of Burgisser, Cucker, and Lotz . For a complete proof of this result, see Theorem 1.1 of , and the upper bound on the constant “C” given in Proposition 5.5.

While Theorem 8 provides a bound on $p_{\text{cover}}(m,n,\phi),$ the formulation of this bound does not provide any intuition of the scaling of $p_{\text{cover}}(m,n,\phi)$ or its dependence on $m$ and $n.$ For this reason, we derive Lemma 4, which is a weaker but more intuitive result. We restate Lemma 4 here for clarity.

Let us simplify the result of Theorem 8. If we assume $m>2n,$ then Hoeffding’s inequality yields

Next, we derive a lower bound as follows:

We have used the fact that $\sqrt{(n/2-1)/\pi}\frac{1}{n-1}>\frac{1}{\sqrt{8n}}$ for $n\geq 4,$ and also the “Wallis ratio” bound $\frac{\Gamma(n/2)}{\Gamma((n-1)/2)}\geq\sqrt{n/2-1}$ . Finally, we plug in the inequality ${{m}\choose{n}}\leq\frac{(em)^{n}}{n^{n}}$. We now have

Now we simplify the integral. Using the identity $(1-x)^{a}<e^{-ax},$ which holds for $x\leq 1,$ we can convert each term in the integrand into an exponential. We do this first with $x=t^{2}$ and then with $x=\frac{1}{\sqrt{8n}}(1-t^{2})^{(n-1)/2}$ to obtain

We then apply the Cauchy-Schwarz inequality to get

Replacing the integral with this bound and using the definition $\epsilon=\cos(\phi)$ yields the result.

References