Fundamental Limits of PhaseMax for Phase Retrieval: A Replica Analysis

I Introduction

we are interested in reconstructing $\boldsymbol{\xi}$ up to a global sign change. This is the real-valued version of the classical phase retrieval problem , which has attracted much renewed interests in the signal processing community in recent years (see, e.g., ). The main challenge of the phase retrieval problem comes from the nonconvex nature of the constraints (1). Recently, a simple yet very effective convex relaxation was independently proposed by two groups of authors . Following , we shall refer to it as the PhaseMax method, which seeks to estimate $\boldsymbol{\xi}$ via a linear programming problem:

Here, the nonconvex equality constraints in (1) have been relaxed to convex inequality constraints. The vector $\boldsymbol{x}_{\text{init}}$ is an initial guess of the target vector $\boldsymbol{\xi}$. In practice, $\boldsymbol{x}_{\text{init}}$ can be obtained if we have additional prior knowledge about $\boldsymbol{\xi}$ (e.g., nonnegativity) or by using a simple spectral method .

The performance of the PhaseMax method has been investigated in (see also ), where the authors provide sufficient conditions for PhaseMax to successfully recover the target vector $\boldsymbol{\xi}$. In this paper, we present an exact performance analysis of the method in the high-dimensional ($n\rightarrow\infty$) limit. In particular, we show that a sharp phase transition phenomenon takes place, with a simple analytical formula characterizing the phase transition boundary.

We shall quantify the performance of PhaseMax in terms of the normalized mean squared error (NMSE), defined as $\text{NMSE}_{n}\overset{\text{def}}{=}{\min\{\mathinner{\!\left\lVert\boldsymbol{\xi}-\widehat{}\boldsymbol{x}\right\rVert}_{2}^{2},\mathinner{\!\left\lVert\boldsymbol{\xi}+\widehat{}\boldsymbol{x}\right\rVert}_{2}^{2}\}}/{\mathinner{\!\left\lVert\boldsymbol{\xi}\right\rVert}_{2}^{2}}$. The NMSE depends on two parameters: the oversampling ratio $\alpha\overset{\text{def}}{=}m/n$, and the quality of the initial guess $\boldsymbol{x}_{\text{init}}$, measured via the input cosine similarity

Taking values between and $1$, the parameter $\rho_{\text{init}}$ assesses the degree of alignment between the true signal vector $\boldsymbol{\xi}$ and the initial guess $\boldsymbol{x}_{\text{init}}$.

As the main contribution of our work, we derive the following exact asymptotic characterization of PhaseMax, under the assumption that the sensing vectors are drawn from the normal distribution:

where $s(\rho_{\text{init}},\alpha)$ is a positive function that can be computed by solving a fixed point equation (see (13), (14) and (16) in Section II-C.) The above expression characterizes the fundamental limits of PhaseMax: for any fixed input cosine similarity $\rho_{\text{init}}$, there is a critical threshold $\alpha_{c}(\rho_{\text{init}})$ such that PhaseMax perfectly recovers $\boldsymbol{\xi}$ if $\alpha>\alpha_{c}(\rho_{\text{init}})$, and that it fails to recover $\boldsymbol{\xi}$ if $\alpha<\alpha_{c}(\rho_{\text{init}})$.

Figure 1 illustrates our asymptotic characterization and compares it with results from numerical simulations. Specifically, the red curve in the figure shows the phase transition boundary $\alpha_{c}(\rho_{\text{init}})$ as a function of the input cosine similarity $\rho_{\text{init}}$, which can be seen to have excellent agreement with the actual performance of the algorithm. In , the authors show that PhaseMax is successful with high probability if

This sufficient condition is plotted as the blue curve in Figure 1. We note that our theoretical prediction significantly reduces the required oversampling ratio as given in (5) for any considered quality of the initial guess vector.

Our analysis is based on the powerful replica method from statistical mechanics. Although certain key steps of the replica method have not yet been mathematically proven, the method has been successful in the analysis of a wide-range of high-dimensional inference problems in signal and information processing (see, e.g., ). Some of its sharp predictions have later been proven through alternative mathematical approaches (e.g., ). In this work, we use the replica method to derive our asymptotic predictions, and corroborate these analytical results—rigorously speaking, conjectures—via numerical simulations.

The rest of this paper is organized as follows. After precisely laying out the various technical assumptions, we present the main results of this work in Section II. Additional numerical results are provided in Section III to validate our theoretical predictions. Section IV concludes the paper. For readers interested in our replica calculations, we present some of our key derivations in the appendix, and leave the full technical details to a follow-up paper.

II Main Results

In what follows, we first state the assumptions under which we derive our analytical predictions.

The sensing vectors $\left\{\boldsymbol{a}_{i}\right\}_{1\leq i\leq m}$ are independent random vectors whose entries are i.i.d. standard normal random variables.

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

Both the target vector $\boldsymbol{\xi}$ and the initial guess $\boldsymbol{x}_{\text{init}}$ are independent from the sensing vectors.

The target vector $\boldsymbol{\xi}$ has a positive cosine with the initial guess $\boldsymbol{x}_{\text{init}}$.

$\mathinner{\!\left\lVert\boldsymbol{\xi}\right\rVert}_{2}=\mathinner{\!\left\lVert\boldsymbol{x}_{\text{init}}\right\rVert}_{2}=\sqrt{n}$.

Note that the last two assumptions can be made without loss of generality, since $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are both valid targets and thanks to the scale invariant nature of the convex optimization problem in (2), respectively.

II-B The Boltzmann Distribution

The first step of our replica analysis is to “soften” the optimization problem (2) via a probability distribution. To that end, we introduce the following function

where $U(x)$ represents the unit-step function, i.e., $U(x)=1$ if $x\geq 0$ and $U(x)=0$ otherwise. Clearly, the convex optimization problem (2) is equivalent to minimizing the function $\mathcal{H}$ over the variable $\boldsymbol{x}$. Now consider the following probability distribution

where in reaching (11) we have used the assumption that $\mathinner{\!\left\lVert\boldsymbol{\xi}\right\rVert}_{2}=\sqrt{n}$. Thus, the task of analyzing the asymptotic performance of PhaseMax boils down to calculating the values of $q^{\ast}$ and $\vartheta^{\ast}$, which we do next by using the replica method.

II-C Asymptotic Predictions via the Replica Method

The challenge here is to compute the partition function $Z_{n}(\beta)$, which involves a high-dimensional integration. And this is where the replica method comes in. Using this method, we can calculate $f(\beta)$ for all $\beta>0$. In particular, its limit as $\beta\to\infty$ can be derived as

To solve the extremization problem in (II-C), we set the gradient with respect to the variables to zero, which leads to a set of nonlinear saddle point equations. After some further simplifications, we can eliminate the variables $\vartheta$, $\chi$, $\widehat{\vartheta}$, $\widehat{q}$ and $\widehat{\chi}$ and just need to study a simple fixed-point equation:

where $c=\tan(\pi/\alpha)/2$ is a constant,

The solution to the above equations then gives us the key parameters of interest $\vartheta^{\ast}$ and $q^{\ast}$ as defined in (10), from which we can compute the asymptotic NMSE by using (11).

We observe that $q=1$ is always a fixed point of (13). However, for any fixed $\rho_{\text{init}}$ and when we reduce the oversampling ratio $\alpha$ to below a threshold, a second fixed point emerges and the original solution $q=1$ becomes unstable. This is indeed the origin of the phase transition. To locate the phase transition boundary, we study the stability of the solution $q=1$. Specifically, by the definition of $h(q)$, we can verify that \dfrac{\operatorname{d\!}{}h}{\operatorname{d\!}{q}}\mathinner{\bigr{\rvert}}_{q=1}\equiv 1 and

Thus, the solution $q=1$ becomes unstable (i.e. a phase transition happens) when

which is exactly when \dfrac{\operatorname{d\!}{{}^{2}}h}{\operatorname{d\!}{q^{2}}}\mathinner{\bigr{\rvert}}_{q=1} changes its sign. Finally, the function $s(\rho_{\text{init}},\alpha)$ in (4) can be obtained as

where $q^{\ast}$ is the stable solution of $\eqref{eq:q_fixed_point}$ and $\vartheta^{\ast}$ is given by (14).

III Numerical Results

In this section, we present additional numerical results to verify our analytical predictions found through the replica method. In all of our experiments, we solve the convex optimization problem (2) using the approach presented in where the signal dimension is set to $n=1000$. The results are also averaged over $50$ independent Monte Carlo trials.

Our first simulation example, shown in Figure 2, studies the performance of our analytical prediction of the NMSE [see (16)] as a function of the input cosine similarity $\rho_{\text{init}}$ for two different values of the oversampling ratio: $\alpha=2.5$ and $\alpha=3.5$, respectively. As seen from the figure, the theoretical prediction obtained by the replica method is in excellent agreement with the experimental results obtained by numerically solving the convex optimization problem (2). The results also validate our theoretical prediction of the phase transition points: the critical input cosine similarity corresponding to each value of $\alpha$ is $\rho_{\text{init}}(\alpha=2.5)=0.769$ and $\rho_{\text{init}}(\alpha=3.5)=0.533$.

A different example is shown in Figure 2, where we examine the performance of our analytical prediction of the NMSE as a function of the oversampling ratio $\alpha$ for two different values of the input cosine similarity: $\rho_{\text{init}}=0.25$ and $\rho_{\text{init}}=0.35$, respectively. Again, as seen from the figure, our theoretical results can accurately predict the actual performance of the algorithm.

IV Conclusion

We presented in this paper an exact characterization of the performance of the PhaseMax method for phase retrieval. Our replica analysis leads to an analytical formula for the asymptotic normalized MSE of the estimate given by PhaseMax in the high-dimensional limit. It also reveals a sharp phase transition phenomenon: for PhaseMax to succeed, the oversampling ratio must be above a critical threshold, given as a function of the input cosine similarity. Simulation results confirm the validity of our theoretical predictions. They also show that our theoretical results significantly reduce the required oversampling ratio given by an existing sufficient condition in the literature.

Appendix A Technical Details

This appendix provides a sketch of our derivations leading to (II-C). To start, we write the partition function $Z_{n}(\beta)$ as

Using the replica trick , the free energy density for any given parameter $\beta$ can be expressed as follows

where $\boldsymbol{x}_{a}$ denotes the $a$th replica signal and where the expectation is over the random vector $\boldsymbol{a}$ which is normally distributed with zero mean and covariance matrix ${\bf I}_{n}$.

where $\alpha=\frac{m}{n}$ denotes the oversampling ratio and where the function $\mathcal{G}$ can be expressed as follows

with $u_{0}=\boldsymbol{a}^{T}\boldsymbol{\xi}/\sqrt{n}$ and $u_{a}=\boldsymbol{a}^{T}\boldsymbol{x}_{a}/\sqrt{n}$, for all $a$. Furthermore, using a result in large deviation theory known as the Gartner-Ellis theorem , the rate function $\mathcal{I}$ can be expressed as the Fenchel–Legendre transform of a cumulant generating function. Specifically, the rate function $\mathcal{I}$ can be expressed as follows

where the function $\mathcal{\lambda}$ represents the cumulant generating function and is given by

where $s_{0}$, $t$ and $\{s_{a},1\leq a\leq n\}$ are i.i.d. Gaussian random variables with zero mean and unit variance. Using the introduced representation of the random variables $u_{0}$ and $u_{a}$, the free energy density can be rewritten as follows

References