Phase Retrieval via Linear Programming: Fundamental Limits and Algorithmic Improvements

I Introduction

Among the most well-established methods are those based on semidefinite relaxation (e.g., .) Such convex optimization methods operate by lifting the original $n$-dimensional natural parameter space to a higher dimensional matrix space. Unfortunately, the increase in the dimensionality introduces challenges in computational complexity and memory requirement for the resulting algorithms. Subsequent works suggest going around this issue by developing nonconvex formulations of the phase retrieval problem and solution algorithms that start with a careful spectral initialization , which is then iteratively refined by a gradient-descent-like scheme of low computational complexity.

More recently, an alternative convex formulation of the phase retrieval problem in the original $n$-dimensional parameter space was independently proposed by two groups of authors . The resulting method, referred to as PhaseMax in , relaxes the nonconvex equality constraints in (1) to convex inequality constraints, and solves the following linear program:

Here, $\boldsymbol{x}_{\text{init}}$ represents an initial guess (or “anchor vector”) that is correlated with the target vector $\boldsymbol{\xi}$.

Despite its simple formulation, PhaseMax has strong theoretical performance guarantees. Existing analysis shows that PhaseMax achieves exact signal recovery from a nearly optimal number of random measurements. Specifically, in the case when the sensing vectors are drawn from the Gaussian distribution, the required number of measurements for perfect reconstruction is shown to be linear with respect to the underlying dimension, i.e., $m=c\,n$ for some constant $c$ that depends on the quality of the initial vector $\boldsymbol{x}_{\text{init}}$. The analysis in gives various upper bounds on the constant $c$. The exact value of $c$, namely the sharp phase transition threshold, is predicted in a recent work by a subset of the authors of the current paper, but the analysis in uses the non-rigorous replica method from statistical physics.

I-B Contributions

Our main contributions in this paper are two-fold.

1. Exact recovery guarantees. We present an exact performance analysis of the PhaseMax method for the (real-valued) phase retrieval problem with Gaussian sensing vectors in the large system limit. When $m,n\rightarrow\infty$ at a proportional ratio $\alpha=m/n$, we rigorously establish the exact phase transition threshold. Furthermore, in the regime where perfect recovery is not feasible, we derive asymptotically exact formulas for the normalized mean squared error (NMSE), defined as

Our formulas reveal the precise dependence of the NMSE on the oversampling ratio $\alpha$ and on the quality of the initial guess $\boldsymbol{x}_{\text{init}}$ as measured via the input cosine similarity

Our main results can be summarized by the following asymptotic characterization of the NMSE:

and $f(\rho_{\text{init}},\alpha)$ is explicitly determined by solving a one-dimensional deterministic fixed point equation [see (8) and Theorem 1.]

Our analysis builds upon the recently developed convex Gaussian min-max theorem (CGMT) , which involves a tight version of a classical Gaussian comparison inequality. The CGMT framework has been successfully applied to derive precise performance guarantees for structured signal recovery under (noisy) linear Gaussian measurements, e.g., . In , the CGMT is used to study signal recovery from a class of non-linear measurements. However, this excludes magnitude-only or quadratic measurements that are relevant for the phase retrieval problem considered here.

The precise nature of our results serves to tighten up the previously known performance bounds of PhaseMax . They also exactly match and thus rigorously verify the predictions in obtained from the non-rigorous replica method from statistical physics. In fact, to the best of our knowledge, this is the first exact performance analysis of any of the existing solution methods for the phase retrieval problem.

2. From precise analysis to algorithmic improvements. As the second contribution of this paper, we propose a new nonconvex formulation and an efficient iterative algorithm for the phase retrieval problem. Our new formulation is inspired by PhaseMax, the key idea of which is to relax the nonconvex equality constraints in (1) to convex inequality constraints. The intersections of all these inequality constraints form a high-dimensional (random) polytope. Our analysis of the PhaseMax method provides useful insights on the exact high-dimensional geometry of that random polytope. These insights then lead us to a novel nonconvex formulation of the phase retrieval problem, as follows:

Note that (6) is indeed a nonconvex problem, as we aim to maximize a convex function over a convex domain. We devise an efficient iterative method, which we call PhaseLamp, to solve (6). The name comes from the fact that the algorithm is based on the idea of successive linearization and maximization over a polytope, where in each step we solve a PhaseMax problem with the initialization given by the estimate from the previous iteration.

We prove that the proposed PhaseLamp method has (strictly) superior recovery performance over PhaseMax. Specifically, we show that a sufficient condition for PhaseLamp to perfectly recover the target signal $\boldsymbol{\xi}$ is

where $\rho_{s}(\alpha)$ is determined explicitly by solving a one-dimensional deterministic fixed point equation (see (17) and Theorem 3.) In particular, $\rho_{s}(\alpha)$ is strictly smaller than $\rho_{c}(\alpha)$ as defined in (5). This is illustrated through a numerical example shown in Figure 1. We can see that the proposed PhaseLamp method significantly improves the recovery performance of the PhaseMax method, especially in the more challenging, and arguably the more practically relevant regime of small input cosine similarities $\rho_{\text{init}}$. Moreover, although (7) is only a sufficient condition, it nevertheless provides a good estimate of the actual performance of the algorithm.

II Precise Performance Analysis of PhaseMax

The asymptotic predictions derived in this paper are based on the following assumptions.

The sensing vectors $\left\{\boldsymbol{a}_{i}\right\}_{1\leq i\leq m}$ are known and drawn independently from a Gaussian distribution with zero mean and covariance matrix $\boldsymbol{I}_{n}$.

$m=m(n)$ with $\alpha_{n}=m(n)/n\rightarrow\alpha>0$ as $n\rightarrow\infty$.

Both the target signal vector $\boldsymbol{\xi}$ and the initial guess vector $\boldsymbol{x}_{\text{init}}$ are independent from the sensing vectors $\left\{\boldsymbol{a}_{i}\right\}_{1\leq i\leq m}$.

For convenience, we shall also assume that the initial guess vector $\boldsymbol{x}_{\text{init}}$ has a positive correlation with the target signal vector $\boldsymbol{\xi}$, i.e., $\boldsymbol{\xi}^{T}\boldsymbol{x}_{\text{init}}>0$. This can be made without loss of generality since the vectors $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are both valid target vectors.

II-B Fundamental Limits of PhaseMax

In this section, we characterize the asymptotic NMSE of the PhaseMax method under the stated assumptions. In particular, our results point out necessary and sufficient conditions on the oversampling ratio $\alpha$ and on the cosine similarity $\rho_{\text{init}}$ for perfect recovery.

In order to state our results we need a few definitions. For any fixed cosine similarity $\rho_{\text{init}}$ and fixed oversampling ratio $\alpha>2$, define $s^{*}$ as follows:

where the function $g_{\alpha}(s)\mathrel{\mathop{\mathchar 58\relax}}\rightarrow(0,\infty)$ (parametrized by $\alpha$) is given by

with $c_{\alpha}=1/\text{tan}\left(\pi/\alpha\right)$ and

For any fixed input cosine similarity $\rho_{\text{init}}>0$ and any fixed oversampling ratio $\alpha>2$, let $s^{\ast},r^{\ast}$ be defined as in (8) and (11), respectively. Then, under the assumptions in Section II-A, the NMSE of the PhaseMax method converges in probability as follows:

The proof of Theorem 1 is based on the CGMT. To streamline our presentation, we postpone a sketch of the proof to the appendix. Theorem 1 accurately predicts the NMSE of PhaseMax in the large system limit. The prediction is expressed in terms of $s^{\ast}$, the solution to the one-dimensional deterministic maximization problem in (8). It can be shown that this optimization problem is concave and that $s^{\ast}$ can be uniquely determined by a fixed point equation.

A sketch of the proof of Theorem 2 can be found at the end of the appendix. The theorem establishes a precise phase transition behavior on the performance of PhaseMax: for any fixed oversampling ratio $\alpha>2$, there is a critical cosine similarity $\rho_{c}(\alpha)$ such that the PhaseMax method perfectly recovers the target signal vector $\boldsymbol{\xi}$ if and only if $\rho_{\text{init}}>\rho_{c}(\alpha)$.

II-C Numerical Simulations

In this section, we present simulation results that verify the validity of our predictions given in Theorems 1 and 2. We solve the convex optimization problem (2) using the technique introduced in . The signal dimension is set to $n=1000$.

In the first example, we investigate the performance of our asymptotic predictions for PhaseMax given in (12). Specifically, we compare the asymptotic predictions against simulation results for different values of the input cosine similarity $\rho_{\text{init}}$ and the oversampling ratio $\alpha$. Figure 2 illustrates the NMSE of PhaseMax as a function of the input cosine similarity given in (3), for two different values of the oversampling ratio. It can be noticed that our asymptotic predictions of the PhaseMax performance obtained using the CGMT perfectly match the asymptotic predictions derived using the non-rigorous replica method . The asymptotic predictions are also in excellent agreement with the actual performance of the PhaseMax method in finite dimensions. Figure 2 further shows that the critical cosine similarity for the considered oversampling ratio $\alpha$ is given by $\rho_{\text{init}}(\alpha=3)\approx 0.63$ and $\rho_{\text{init}}(\alpha=5)\approx 0.37$, respectively, which perfectly matches the theoretical predictions given in Theorem 2.

Figure 2 provides an additional example, where we plot the NMSE of the PhaseMax method as a function of the oversampling ratio, for two different values of the input cosine similarity. Again, the asymptotic performance obtained in our analysis perfectly matches the actual performance of the algorithm.

III Algorithmic Improvements

This section proposes an efficient iterative algorithm to solve the norm maximization problem formulated in (6). The optimization problem (6) consists of maximizing a convex function over a convex feasibility set. Hence, it is nonconcave where the cost function can be written as a difference of concave functions. To solve this problem, we propose the following scheme, named PhaseLamp, based on the idea of successive linearization and maximization over a polytope:

In essence, at each iteration we approximate (i.e., linearize) the cost function of (6) via

There are several ways to interpret the proposed PhaseLamp algorithm. First, it can be viewed as an iterative and bootstrapped version of the PhaseMax method (2) where at each iteration the previous optimal solution is used as an (improved) initial guess of the target signal vector $\boldsymbol{\xi}$. Second, PhaseLamp is a special case of a minorize-maximization (MM) algorithm . To see this, we note from the convexity of the cost function $\mathinner{\!\left\lVert\boldsymbol{x}\right\rVert}_{2}^{2}$ that

The PhaseLamp procedure consists of iteratively maximizing the lower bound in (15) over the convex feasibility set in (6). One particular property of the MM procedure is that it guarantees that the objective value of the optimization problem (6) is nondecreasing, i.e.,

Due to the nonconcavity of the maximization problem (6), the proposed iterative algorithm is not guaranteed to converge to the global optimal solution of (6). One particular property of the fixed points of the optimization problem (14) is that it is an extreme point of the feasibility set given in (6).

III-B Performance Guarantees for PhaseLamp

Using the analysis strategy that leads to Theorems 1 and 2, we are further able to derive a sufficient condition for PhaseLamp to perfectly recover the target signal $\boldsymbol{\xi}$.

Again, we first need a few definitions. One can show that, for any $\alpha>2$, the equation

has a unique solution in the interval $\theta\in(0,\pi/2)$. We denote that solution by $\theta^{\ast}_{\alpha}$. Let

where $c_{\alpha}=1/\text{tan}\left(\pi/\alpha\right)$, and

where $g_{\alpha}$ is the function defined in (9). We are now ready to state the main theorem of this section.

The proof of Theorem 3 is based on CGMT and the properties of the fixed points of the optimization problem (14). Due to space constraint, we defer the proof of this theorem to the long version of the current paper. Unlike the asymptotically exact characterization given in Theorem 2, the condition given in Theorem 3 is sufficient but not necessary. However, as shown in Figure 1 and the additional simulation results given in the next section, the condition in (20) provides a reasonably tight bound on the actual performance of the PhaseLamp algorithm.

III-C Numerical Results

We present some numerical results to illustrate the performance of the proposed PhaseLamp method and our theoretical predictions given in Theorem 3. In our experiments, the signal dimension is set to $n=1000$. Figure 3 plots the NMSE as a function of the oversampling ratio for two different values of the input cosine similarity ($\rho_{\text{init}}=0.1$ and $\rho_{\text{init}}=0.3$, respectively.)

We observe that the proposed PhaseLamp method indeed outperforms the original PhaseMax method, and the amount of improvement is greater when the input cosine similarity is smaller. Specifically, for $\rho_{\text{init}}=0.1$, the empirical minimum sampling ratio for PhaseLamp to perfectly recover $\boldsymbol{\xi}$ is at $\alpha\approx 3.3$, whereas PhaseMax requires $\alpha\approx 18.2$ . [The latter point is not shown in Figure 3.] Moreover, Figure 3 also shows that the sufficient condition developed in Theorem 3 for PhaseLamp provides a good estimate of the actual performance of the proposed algorithm. For example, Figure 3 demonstrates that the actual critical oversampling ratio of PhaseLamp is at $\alpha_{c}\approx 2.9$, whereas the sufficient oversampling ratio as given in Theorem 3 is $\alpha_{s}\approx 3.3$.

IV Conclusion

We presented in this paper an asymptotically exact characterization of the performance of the PhaseMax method for phase retrieval. Specifically, our analysis reveals a sharp phase transition behavior in the performance of the method as one varies the oversampling ratio and the input cosine similarity. Our analysis is based on the CGMT, and the results match previous predictions derived from the non-rigorous replica method. Moreover, we also presented a new nonconvex formulation of the phase retrieval problem and PhaseLamp, an iterative algorithm based on linearization and maximization over a polytope. We provided a sufficient condition for PhaseLamp to perfectly retrieve the target vector. Simulation results confirm the validity of our theoretical predictions. They also show that the proposed iterative algorithm significantly improves the recovery performance of the original PhaseMax method.

In this appendix we provide proof sketches for Theorems 1 and 2.

Similarly, define $\eta_{1}$ and $\widetilde{\boldsymbol{\eta}}$ such that

Finally, for a vector $\boldsymbol{c}$, we let $\mathinner{\!\left\lvert\boldsymbol{c}\right\rvert}$ and $\text{sign}(\boldsymbol{c})$ to denote its component-wise absolute value and sign, respectively. Also, we let $\min(\boldsymbol{c})$ return the minimum value in the vector, and, $\boldsymbol{z}=\boldsymbol{c}\wedge\mathbf{0}$ be a vector such that $\boldsymbol{z}_{i}=\min(\boldsymbol{c}_{i},0)$.

-B Convex Gaussian Min-Max Theorem (CGMT)

The proof follows the CGMT framework. For ease of reference we summarize here the essential ideas of the framework; please see [16, Section 6] for the formal statement of the theorem and further details. The CGMT associates with a primary optimization (PO) problem a simplified auxiliary optimization (AO) problem from which we can tightly infer properties of the original (PO), including the optimal cost and the optimal solution. The two problems are of the following form:

In words, concentration of the optimal cost of the AO problem around $\mu$ implies concentration of the optimal cost of the corresponding PO problem around the same value $\mu$. Moreover, starting from (23) and under strict convexity conditions, the CGMT shows that concentration of the optimal solution of the AO problem implies concentration of the optimal solution of the PO to the same value. For example, if minimizers of (22) satisfy $\mathinner{\!\left\lVert\boldsymbol{w}^{\ast}(\boldsymbol{g},\boldsymbol{h})\right\rVert}_{2}\to\zeta^{\ast}$ for some $\zeta^{\ast}>0$, then, the same holds true for the minimizers of (21): $\mathinner{\!\left\lVert\boldsymbol{w}^{\ast}(\boldsymbol{C})\right\rVert}_{2}\to\zeta^{\ast}$ [16, Theorem 6.1(iii)]. Thus, one can analyze the AO to infer corresponding properties of the PO, the premise being of course that the former is simpler to handle than the latter.

-C CGMT for the PhaseMax Method

We apply the CGMT to characterize the asymptotic NMSE of the PhaseMax optimization in (2) as in (12) and (13). In this section, we write the PhaseMax optimization in the form of a PO as in (21), which in turn leads to a corresponding AO optimization problem. For these problems, we can show that the conditions of the CGMT on convexity and compactness are satisfied.

First, we appropriately write the linear program in (2) as a minmax program. Start with its dual:

Since $m>n$ (recall: $\alpha>2$) both the primal and the dual are bounded feasible with probability one, and strong duality holds . Therefore, (24) is equivalent to the following

Further note that the constraint sets are convex compact and $\psi$ is concave-convex on $\mathcal{S}_{\boldsymbol{x}}\times(\mathcal{S}_{\boldsymbol{\lambda}}\times\mathcal{S}_{\boldsymbol{\mu}})$, where $\boldsymbol{x}=[x_{1}~{}\widetilde{\boldsymbol{x}}^{T}]^{T}$ and $\mathcal{S}_{\boldsymbol{x}}=\mathcal{S}_{x_{1}}\times\mathcal{S}_{\widetilde{\boldsymbol{x}}}$.

We are now ready to formulate the corresponding AO problem:

Following the CGMT framework we proceed onwards with analyzing (-C).

-D Analysis of the Auxiliary Optimization Problem

Consider the following change of variables: $\boldsymbol{v}=\boldsymbol{\lambda}-\boldsymbol{\mu}$ and $\boldsymbol{b}=\boldsymbol{\lambda}+\boldsymbol{\mu}$. To respect the nonnegativity of $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$, it must be that $\boldsymbol{b}\geq|\boldsymbol{v}|$. In fact, it can be checked that the optimal solution for $\boldsymbol{b}$ is $\boldsymbol{b}=|\boldsymbol{v}|$. Thus, the optimization problem (-C) can be reduced to the following:

Next, observe that if we fix $\mathinner{\!\left\lvert\boldsymbol{v}\right\rvert}$, then the optimal $\boldsymbol{v}$ satisfies $\text{sign}(\boldsymbol{v})=-\text{sign}\left(\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}\right\rVert}_{2}\boldsymbol{g}+\boldsymbol{q}x_{1}\right)$ which simplifies the optimization to the following

Next, in the optimization above one can fix the norm of $\boldsymbol{v}$ and optimize over its direction. Omitting some details, the optimization becomes

The final step in simplifying the AO problem is as follows. For fixed value of $x_{1}$ (say $x_{1}=s>0$), and for fixed norm of $\widetilde{\boldsymbol{x}}$ (say, $\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}\right\rVert}_{2}=r$), we optimize over the direction of $\widetilde{\boldsymbol{x}}$. It can be shown that this optimization further reduces (29) to the following two-dimensional optimization problem:

-D2 Convergence Analysis

Now that we have simplified the AO to a maximization problem over only two scalar variables as in (30), we are ready to study its asymptotic behavior in the regime $m,n\rightarrow\infty,m/n\rightarrow\alpha$. Specifically, it can be shown that the optimization problem (30) converges to the following deterministic optimization problem

The full technical details of obtaining the convergence result in (31) are deferred to the full version of the paper. In short, pointwise convergence of the objective function of (30) to (31) for fixed $s$ and $r$ follows easily from the weak law of large numbers. The corresponding convergence of the optimal costs requires proof of uniform convergence, which follows by pointwise convergence and concavity of the objective function [26, Lemma 7.75]. We call the deterministic two-dimensional optimization problem in (31) as the scalar performance optimization (SPO); according to the CGMT solving the SPO allows us to conclude on the asymptotic performance of the PhaseMax problem (cc. the PO).

-D3 Solving the scalar performance optimization

Recall that the SPO in (31) is the converging limit of the AO in (-C). Specifically, the optimization variables $s$ and $r$ in (31) correspond exactly to $x_{1}$ and $\|\widetilde{}\boldsymbol{x}\|_{2}$ in (-C). From this and uniform convergence discussed previously, the optimal values of $s$ and $r$ are the converging limits of $x_{1}$ and of $\|\widetilde{}\boldsymbol{x}\|_{2}$, respectively. In what follows, we solve the SPO problem for the optimal $s$ and $r$. First, using the assumption of the theorem that $\alpha>2$, it can be shown that the feasible set of (31) is nonempty iff $\mathinner{\!\left\lvert s\right\rvert}\leq 1$. Second, for fixed $|s|\leq 1$, the problem

At this point, note that we can always find a large enough constant $\widetilde{B}>0$ such that $r^{\ast}(s)<\widetilde{B}$ for all $|s|<1$. Therefore, choosing $B$ in (31) such that $B=\widetilde{B}$ guarantees that the optimal value of $r$ in (31) is given by (34). Substituting this value back in (31), we can now optimize over $s$ by solving the following:

A few algebra manipulations show that (35) is equivalent to (8) in the statement of the theorem. To show the equivalence, further note that $\eta_{1}$ and $\widetilde{\boldsymbol{\eta}}$ in (35) are related to the input cosine similarity $\rho_{\text{init}}$, defined in (3), as follows (recall: $\boldsymbol{\xi}=\boldsymbol{e}_{1}$.),

Finally, note that the optimization in (35) [eqv., in (8)] inherits the concavity of (31), i.e., it is a concave program.

-D4 Phase transition calculations

In this section, we compute the phase transition boundary of the PhaseMax method. Our goal is to find necessary and sufficient conditions under which the solution $\hat{}\boldsymbol{x}$ of the PhaseMax is, with high probability, equal to $\boldsymbol{\xi}=\boldsymbol{e}_{1}$. Mapping this to the SPO in (35) [eqv., see (8)], we seek conditions under which $s^{\ast}=1$ and $r^{\ast}=0$. From concavity, this happens if and only if the derivative of the cost function of the optimization problem (8) at $s=1$ is nonnegative. Hence, the necessary and sufficient condition for perfect recovery of the PhaseMax method is given by

for $\alpha>2$. Equivalently, the oversampling ratio $\alpha$ and the input cosine similarity given in (3) must satisfy

This then gives us the statement of Theorem 2.

References