Phase Retrieval via Polytope Optimization: Geometry, Phase Transitions, and New Algorithms

I Introduction

Among the most well-established methods are those based on semidefinite relaxation (e.g., ), which operate by lifting the original $n$-dimensional natural parameter space to a higher dimensional matrix space. Despite the strong theoretical performance guarantees enjoyed by these convex-relaxation methods, the aforementioned lifting step significantly increases the computational complexity and memory requirement for the resulting algorithms. To address these challenges, recent work studies algorithms that directly solve the nonconvex formulations of the phase retrieval problem. Typically, such nonconvex methods follow a two-step approach, combining a careful initialization step with further local refinement such as iterative gradient descent .

Taking a different approach, two groups of authors independently proposed a simple yet highly effective scheme that is based on convex programming in the original $n$-dimensional signal space. 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:

I-B Contributions

In this paper, we present an exact performance analysis of the PhaseMax method in the high-dimensional ($n\to\infty$) limit. In particular, we show that a phase transition phenomenon takes place, with a simple analytical formula characterizing the phase transition boundary. Moreover, we extend the idea of PhaseMax by proposing a new nonconvex formulation of the phase retrieval problem and an accompanying iterative algorithm. We show that this new algorithm, which we call PhaseLamp, has provably superior recovery guarantees over the original PhaseMax method. In what follows, we highlight our main results with more technical details.

1. Exact performance analysis of PhaseMax. We quantify the performance of PhaseMax in terms of the normalized mean squared error (NMSE), defined as

The NMSE depends on two parameters: the oversampling ratio

and the quality of the initial guess $\boldsymbol{x}_{\text{init}}$, measured via the input cosine similarity

Note that the parameter $\rho_{\text{init}}$ quantifies the degree of alignment between the target vector $\boldsymbol{\xi}$ and the initial guess $\boldsymbol{x}_{\text{init}}$.

As one of the main contributions of our work, we derive the following asymptotically exact characterization of PhaseMax, under the assumption that the sensing vectors are drawn from the normal distribution: as $m,n\rightarrow\infty$ with their ratio $\alpha$ fixed,

and $f(\rho_{\text{init}},\alpha)$ is a positive function that can be explicitly determined by solving a one-dimensional deterministic fixed point equation (see Theorem 2). We note that the asymptotic characterization in (4) establishes an exact phase transition boundary on the minimum required number of measurements for PhaseMax to be successful: for any fixed sampling ratio $\alpha$, there is a critical threshold $\rho_{\text{c}}(\alpha)$ such that PhaseMax perfectly recovers $\boldsymbol{\xi}$ if and only if the input cosine similarity $\rho_{\text{init}}>\rho_{c}(\alpha)$.

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 $\rho_{c}(\alpha)$, 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

which is plotted as the blue curve in Figure 1. We can see that our theoretical prediction serves to tighten the sufficient condition given in (6).

2. Nonconvex formulation and new algorithms. The insights gained from the exact analysis of PhaseMax lead us to a new nonconvex formulation of the phase retrieval problem:

Note that (7) is indeed a nonconvex problem, as we aim to maximize a convex function over a convex domain. We propose an efficient iterative method, which we call PhaseLamp, to solve (7). 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 complement PhaseLamp with performance guarantees. Due to the iterative nature of PhaseLamp, the analysis here is more challenging than that of PhaseMax. By carefully characterizing the stationary points of (7), we prove that a sufficient condition for PhaseLamp to perfectly recover the target signal $\boldsymbol{\xi}$ (or, $-\boldsymbol{\xi}$) is

where $\rho_{s}(\alpha)$ is determined explicitly by solving a one-dimensional deterministic equation (see (4) and Theorem 5.) Importantly, $\rho_{s}(\alpha)$ is strictly smaller than $\rho_{c}(\alpha)$ as defined in (5). Therefore, the proposed PhaseLamp method has (strictly) superior recovery performance over PhaseMax with respect to the minimum number of measurements needed to guarantee perfect solution of (1).

We illustrate this improvement in Figure 1, where it is shown that PhaseLamp has significantly better recovery performance, especially in the more challenging, and arguably the more practically relevant regime of small input cosine similarities $\rho_{\text{init}}$. Moreover, the numerical simulations shown at the same figure suggest that, although (8) is only a sufficient condition, it nevertheless provides a good estimate of the actual performance of the algorithm. Finally, as yet another variation on the theme of performing phase retrieval via polytope optimization, we propose in Section IV-C a weighted version of PhaseLamp. This new version is empirically shown to further outperform PhaseLamp.

Although our theoretical analysis is carried out for generic Gaussian measurements, the proposed PhaseLamp algorithm and its weighted version perform well under more realistic measurement models that arise in phase retrieval applications. In Figure 2, we compare the performance of PhaseLamp to PhaseMax and three other leading methods in the literature, where the measurement model corresponds to coded-diffraction patterns . In this experiment, PhaseLamp successfully recovers the underlying image and outperforms the other competing methods. More details about the setup of this experiment as well as additional numerical results can be found in Section V.

I-C Related Work

The performance of PhaseMax has been previously investigated in the literature. Existing analysis shows that PhaseMax can achieve exact signal recovery from a nearly optimal number of random linear measurements. Specifically, in the case where 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 guess $\boldsymbol{x}_{\text{init}}$. The analysis in gives various upper bounds on the constant $c$. In a more recent work , a subset of the authors of the current paper were able to pinpoint the exact value of $c$, but the analysis in uses the nonrigorous replica method from statistical physics. Therefore, the precise nature of the results of our paper serves to(a) tighten up the previously known performance bounds of PhaseMax as given in ; and (b) rigorously verify the predictions based on the replica method given in . Moreover, our novel theoretical analysis builds upon an exact characterization of the geometry of the feasibility set of the PhaseMax problem in (2). This geometric insight plays a key role in both the formulation and the analysis of the improved PhaseLamp method proposed in this paper.

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 nonlinear measurements. However, this excludes magnitude-only or quadratic measurements that are relevant for the phase retrieval problem considered here.

This paper is a significantly extended version of our earlier conference paper , which announced our results with proof sketches. A limitation of our work is that we only consider the real-valued version of the phase retrieval problem. Very recently, our analysis techniques have been extended by the authors of to the complex-valued case. As another limitation, we assume that we have access to noiseless measurements as in (1). However, we believe that our technical approaches based on CGMT can be generalized to study the case of noisy measurements as well as robust versions of PhaseMax (see, e.g., ).

I-D Paper Outline

The rest of the paper is organized as follows. Central to our work is an exact characterization of the geometric properties of the feasibility polytope of the optimization in (2). Thus, we start by presenting in Section II a rigorous high dimensional analysis of this polytope. Section III focuses on PhaseMax, where we establish accurate performance guarantees for this method in the high-dimensional limit. The new nonconvex formulation (7) and the accompanying PhaseLamp algorithm are introduced in Section IV. We also provide sufficient conditions for PhaseLamp to achieve perfect recovery. Additional simulation results are shown in Section V, comparing PhaseLamp (and its weighted variation) with several other existing algorithms for the phase retrieval problem. Section VI collects the proofs and technical details of all the results introduced in the previous sections. We conclude the paper in Section VII.

I-E Technical Assumptions and Notation

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 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$, where $\alpha>1$.

The initial guess $\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$.

The assumption in (A.3) can be made without loss of generality, as both $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are valid target signals. Similarly, the assumptions made in (A.4) only serve to simplify the notation but they are not restrictive either, thanks to the rotational invariance of the Gaussian distribution and since the optimization problem (2) is scale invariant.

For any set $\mathcal{A}$ in a finite-dimensional Euclidean space, we define its $L_{2}$ norm as

II Polytope geometry

In this section, we study the geometry of the feasibility set of PhaseMax in (2), which is given as follows:

Under the assumption of Gaussian sensing vectors, ${\mathcal{C}}_{\rm feas}$ forms a high-dimensional random polytope. It is essential, both for the analysis of PhaseMax and also for motivating PhaseLamp, to understand the exact structure of the above polytope.

Before we delve into the details of our analysis, we provide a visualization of ${\mathcal{C}}_{\rm feas}$ via a simulation example, which aims to explain intuitively why PhaseMax is expected to succeed at recovering the unknown signal as the number of measurements increases. Specifically, Figure 3 shows a projection of the random polytope ${\mathcal{C}}_{\rm feas}$ as a function of the number of measurements $m$.

Note that as the number $m$ of constraints increases, the feasibility set ${\mathcal{C}}_{\rm feas}$ looks more and more like a needle pointed towards the target signal vectors $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$. This observation suggests the existence of a phase transition behavior in the performance of PhaseMax method. In particular, the target signal vector $\boldsymbol{\xi}$ has the highest correlation with the initial guess vector $\boldsymbol{x}_{\text{init}}$ among all the feasible vectors as long as $m$ is sufficiently large (as a function of the correlation of $\boldsymbol{x}_{\text{init}}$ with $\boldsymbol{\xi}$). Of course, if we hope to make this observation rigorous, we need to develop formal analytic results regarding the properties of the random high-dimensional set ${\mathcal{C}}_{\rm feas}$. Despite the challenge of the task at hand, we show in the next sections that this is possible.

II-B The Sufficient Feasible Set

Note that what determines the error in a solution $\boldsymbol{x}$ of either (2) or (7) are the magnitudes of $|\boldsymbol{x}_{1}-1|$ and of $\|\widetilde{\boldsymbol{x}}\|_{2}$. This essentially simplifies our task to that of understanding the geometry of the two dimensional projection of ${\mathcal{C}}_{\rm feas}$:

Assume that the oversampling ratio $\alpha>1$. Define the deterministic set ${\mathcal{D}}_{\rm feas}^{\epsilon}$ as follows

where $\epsilon>0$. Then, for all $\epsilon>0$ it holds that

The take away message of Theorem 1, whose proof is detailed in Appendix A-A, is that the random feasibility set $\mathcal{S}_{\text{feas}}$ is essentially a subset (with high-probability in the large system limit) of any $\epsilon$-perturbation of the following deterministic set

Hence, in order to understand the properties of $\mathcal{S}_{\text{feas}}$, it is essential to study the properties of ${\mathcal{D}}_{\rm feas}$.

We start with a visualization of ${\mathcal{D}}_{\rm feas}$ for different values of the oversampling ratio $\alpha$ in Figure 4. Observe that, for sufficiently large oversampling ratio $\alpha$, ${\mathcal{D}}_{\rm feas}$ looks like a needle pointed towards the target signal vectors $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$. Recall, that this is consistent with the observations of Section II-A. Furthermore, note that ${\mathcal{D}}_{\rm feas}$ is always convex and bounded.

The following lemma, which is proved in Appendix B-A, formalizes these observations.

The deterministic set ${\mathcal{D}}_{\rm feas}$ satisfies the following properties:

For $\alpha\geq 2$ and $s\in[-1,~{}1]$, the maximum radius $r_{\text{max}}(s)$ of the set ${\mathcal{D}}_{\rm feas}$ satisfies $c_{d}(s,r_{\text{max}}(s))=r_{\text{max}}(s)/\alpha$. Moreover, for $s=0$, $r_{\text{max}}(0)>0$.

For any $\alpha>1$ and $\epsilon>0$, the set ${\mathcal{D}}^{\epsilon}_{\rm feas}$ is compact. Moreover, it satisfies $\mathcal{D}_{\text{feas}}^{\epsilon_{1}}\subseteq\mathcal{D}_{\text{feas}}^{\epsilon_{2}}$ for any $0<\epsilon_{1}\leq\epsilon_{2}$ and we have

II-C Sufficient Condition for PhaseMax

With Theorem 1 and Lemma 1 at hand, we have established an exact characterization of the high-dimensional geometry of the feasibility set of PhaseMax. Naturally, this leads to a sufficient condition under which its solution ${\widehat{\boldsymbol{x}}}$ is the true unknown vector $\boldsymbol{\xi}$.

All we need in addition to Theorem 1 is the following simple observation regarding $\widehat{\boldsymbol{x}}$, which follows directly by its optimality in the optimization problem in (2).

The optimal solution set of PhaseMax is a subset of the following deterministic set

Without loss of generality, we can assume (due to symmetry) that $\rho_{\text{init}}\geq 0$. Let $\widehat{\boldsymbol{x}}$ be an optimal solution of (2) and partition it as follows $\widehat{\boldsymbol{x}}^{T}={[s~{}~{}{{\widetilde{\boldsymbol{x}}}}]}^{T}.$ From optimality it holds:

This implies that $\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}r\geq\eta_{1}-\eta_{1}s$ with $r=\mathinner{\!\left\lVert{\widetilde{\boldsymbol{x}}}\right\rVert}_{2}$. Recalling that

and rearranging terms, completes the proof of the lemma. ∎

With these at hand, we have shown that in the high dimensional limit the solution of PhaseMax belongs to the intersection of the sets ${\mathcal{D}}_{\text{fp}}(\rho_{\text{init}})$ and ${\mathcal{D}}_{\rm feas}$. Therefore, a natural sufficient condition for perfect recovery is that this intersection only contains the desired points $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$. Proposition 1 below formalizes this geometric condition and Figure 5 serves as a numerical illustration of it.

Note that in the high dimensional limit the solution of PhaseMax should belong to the intersection of the sets ${\mathcal{D}}_{\text{fp}}(\rho_{\text{init}})$ and ${\mathcal{D}}_{\rm feas}$. This fact leads us to a sufficient condition for perfect recovery of PhaseMax as stated in the following Proposition.

Selecting the input cosine similarity in this way guarantees that for any $\epsilon>0$ there exists $\delta>0$ such that

which means that for any $k\geq k_{0}$ and $0\leq\delta\leq\delta_{0}/2$, we have $\mathcal{D}_{\text{feas}}^{\epsilon_{k}}\cap{\mathcal{D}}_{\rm fp}(\rho_{\text{init}})\not\subset\mathcal{B}^{\delta}$. Now, define the following set

III Precise Analysis of PhaseMax

In Section II-C we derived a sufficient condition for perfect recovery of PhaseMax. In this section, we establish a tight such result by further assuming that the initial guess vector $\boldsymbol{x}_{\text{init}}$ is independent of the sensing vectors $\left\{\boldsymbol{a}_{i}\right\}_{1\leq i\leq m}$ and the target signal $\boldsymbol{\xi}$. In particular, we precisely characterize the minimum required number of measurements as a function of the input cosine similarity $\rho_{\text{init}}$ so that PhaseMax finds the true vector $\boldsymbol{\xi}$. Moreover, when this is not possible we precisely quantify the NMSE.

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^{\ast}$ as follows:

where the function $g_{\alpha}$ (parametrized by $\alpha$) is given by

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

We are now ready to state the main result of this section. Its proof uses the recently developed CGMT framework and is deferred to Section VI-C.

Assume that $\boldsymbol{x}_{\text{init}}$ is independent of the sensing vectors $\left\{\boldsymbol{a}_{i}\right\}_{1\leq i\leq m}$ and of the target signal $\boldsymbol{\xi}$. 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 (25) and (28), respectively. Then, the NMSE of the PhaseMax method converges in probability as follows:

Moreover, the optimal cost and the optimal solution $\widehat{\boldsymbol{x}}$ of PhaseMax satisfy the following :

where $\widehat{\boldsymbol{x}}=[s(\widehat{\boldsymbol{x}})~{}~{}\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})^{T}]^{T}$ and $r(\widehat{\boldsymbol{x}})=\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})\right\rVert}_{2}$.

Theorem 2 accurately predicts the NMSE of PhaseMax in the large system limit. The formulae involve solving the one-dimensional deterministic maximization problem in (25). In Section VI-C2 we show that this optimization is strictly concave, thus, $s^{\ast}$ is unique and can be efficiently determined by solving a fixed point equation.

Clearly, we can use the formula on the NMSE given by Theorem 2 to quantify necessary and sufficient conditions under which zero error is achieved. This is the content of the next theorem, which we prove in Section VI-C3.

Theorem 3 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 algorithm perfectly recovers the target signal vector $\boldsymbol{\xi}$ if and only if $\rho_{\text{init}}\geq\rho_{c}(\alpha)$.

III-B Numerical Simulations

The numerical results presented in this section aim to verify the validity of Theorems 2 and 3. First, Figure 6 illustrates the NMSE of PhaseMax as a function of the input cosine similarity given in (3), for two different values of the oversampling ratio $\alpha$. For the simulations, we solve the convex optimization problem (2) using the techniques introduced in and we set the signal dimension as $n=1000$. Note that the asymptotic prediction of Theorem 2 is in excellent agreement with the actual performance of the PhaseMax method in finite dimensions. Of course, the same holds true for the recovery condition of Theorem 3: the theoretical values $\rho_{\text{init}}(\alpha=3)\approx 0.63$ and $\rho_{\text{init}}(\alpha=5)\approx 0.37$ perfectly match with the simulation results Next, Figure 6 plots the NMSE of PhaseMax as a function of the oversampling ratio, for two different values of the input cosine similarity. Again, the figure highlights the sharpness of the results in Theorems 2 and 3.

IV Nonconvex Formulation and New Algorithms

In this section, we propose and study an improved algorithm over PhaseMax, which we call PhaseMax. The natural idea behind PhaseLamp is to solve a sequence of PhaseMax problems. Interestingly, we provide an interpretation of this algorithm as an iterative method for solving the non-convex phase-retrieval problem formulation in (7). This interpretation leads to strong performance guarantees for PhaseLamp in Section IV-B. Finally, in Section IV-C we propose yet one more recovery algorithm, which is also based on optimization over polytopes and which appears to outperform both PhaseMax and PhaseLamp in numerical simulations.

We begin our exposition by arguing in Proposition 2 that, given enough measurements, the solution to the system of quadratic equations in (1) can be found by solving the optimization problem in (7). The proof is in Appendix B-B.

Assume that $\alpha$ satisfies $\alpha>\frac{\pi}{\pi-2}$. Then, the set of optimal solutions of the PhaseLamp problem (7) converges to the set $\{\boldsymbol{\xi},-\boldsymbol{\xi}\}$ in the sense that $\sup_{\widehat{\boldsymbol{x}}\in\mathcal{S}_{\text{lamp}}}(\min\{\mathinner{\!\left\lVert\widehat{\boldsymbol{x}}-\boldsymbol{\xi}\right\rVert}_{2},\mathinner{\!\left\lVert\widehat{\boldsymbol{x}}+\boldsymbol{\xi}\right\rVert}_{2}\})$ converges to zero in probability, where $\mathcal{S}_{\text{lamp}}$ denotes the set of optimal solutions of the PhaseLamp problem given in (7).

According to the proposition, if the number of measurements satisfies $\alpha>\frac{\pi}{\pi-2}\approx 2.752$, then one can hope of solving the phase-retrieval problem by finding the optimal solution of the optimization problem in (7). Unfortunately, (7) is clearly non-convex since it involves maximizing a convex function over a convex set.

In this section, we propose solving (7) by using a standard minorization-maximization (MM) approach as follows. Start by observing that because of convexity the cost function $\mathinner{\!\left\lVert\boldsymbol{x}\right\rVert}_{2}^{2}$ satisfies

Equivalently, the function $m_{k}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{x}\to\boldsymbol{x}_{k}^{T}\boldsymbol{x}_{k}+2\boldsymbol{x}_{k}^{T}\left(\boldsymbol{x}-\boldsymbol{x}_{k}\right)$ is a minorizer of the function $m\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{x}\to\mathinner{\!\left\lVert\boldsymbol{x}\right\rVert}_{2}^{2}$. Moreover, the function $m_{k}$ satisfies $m_{k}(\boldsymbol{x}_{k})=m(\boldsymbol{x}_{k})$. Hence, it is natural to attempt solving (7), via the following iterative scheme

which is of course equivalent to the following:

However, due to the non-convexity of (7), PhaseLamp is not guaranteed in general to converge to the desired global optimal solution of (7). The main theoretical result of this section involves identifying sufficient conditions under which this is indeed the case. Before formalizing those in Section IV-B, it is instructive to consider the performance of PhaseLamp on two different problem instances as shown in Figure 7. Specifically, we present simulation results for the following two cases: (a) $\alpha=3$ and $\rho_{\text{init}}=0.1$ (Figure 7), and (b) $\alpha=4$ and $\rho_{\text{init}}=0.1$ (Figure 7). First, in both instances $\alpha>2.752$; hence Proposition 2 guarantees that the optimal solutions of (7) coincide with the target vectors $\boldsymbol{\xi}$ or $-\boldsymbol{\xi}$. However, as mentioned PhaseLamp is not always guaranteed to find the optimal solutions of (7). For example, it fails to do so in Figure 7, but it succeeds in Figure 7. The sufficient conditions derived in the next section provide rigorous theoretical justifications to these observations.

IV-B Performance Guarantees for PhaseLamp

Clearly, PhaseLamp in the form of (34) can be naturally viewed as an iterative and bootstrapped version of the PhaseMax method (2) where at each iteration $k\geq 1$ , the optimal solution at the previous iteration is used as an (improved) initial guess for a new iteration of PhaseMax. In other words, the cosine similarity $\rho^{k}_{\text{out}}$ between the PhaseMax solution at iteration $k$ and the target signal vector $\boldsymbol{\xi}$ serves as the input cosine similarity $\rho^{k+1}_{\text{init}}$ at iteration $k+1$. One may then imagine leveraging the analysis of PhaseMax in Section III to obtain similar sharp results for PhaseLamp. Unfortunately, more effort and several new arguments are required; the challenge becomes that, after the first iteration of PhaseLamp, the initial guess vector becomes dependent on the sensing vectors $\{\boldsymbol{a}_{i},1\leq i\leq m\}$.

In this section, we overcome these challenges, thus obtaining strong performance guarantees for PhaseLamp. Our arguments are geometric in nature, similar in nature (but somewhat more involved) to the proof of Proposition 1 in Section II-C.

On the one hand, the solution to each iteration of (34) is constrained to live in the feasibility set of PhaseMax. Thus, the same is true for the converging solution (cc. fixed point) of PhaseLamp. Combining this with Theorem 1, which obtains a sharp characterization of the high-dimensional geometry of this feasibility set (cc. its two dimensional projection), we conclude that the fixed points of PhaseLamp belong with probability approaching 1 as $n\rightarrow\infty$ to the set ${\mathcal{D}}_{\rm feas}^{\epsilon}$, for any $\epsilon>0$.

On the other hand, any fixed point $\widehat{\boldsymbol{x}}$ of PhaseLamp satisfies

From this and feasibility of the target vector $\boldsymbol{\xi}$, it holds

Concluding, the fixed points of PhaseLamp belong to the following set

Overall, we have shown that in the limit of high-dimensions, the set of all possible fixed points of PhaseLamp belongs to the intersection of the two sets ${\mathcal{D}}^{\epsilon}_{\rm feas}$ and ${\mathcal{D}}_{\text{opt}}$. The sets ${\mathcal{D}}_{\rm feas}$ and ${\mathcal{D}}_{\text{opt}}$ are illustrated in Figure 8 for $\alpha=7.$ Note that any $\epsilon$-perturbation of the shaded region union the points $(1,0)$ and $(-1,0)$ (ie., the set ${\mathcal{D}}^{\epsilon}_{\rm feas}\cap{\mathcal{D}}_{\text{opt}}$, for any $\epsilon>0$) represents the set of possible fixed points of PhaseLamp. Clearly PhaseLamp is successful when it escapes the shaded region of “bad” stationary points. Hence, the question becomes: for given $\alpha$, what values of initial correlation $\rho_{\text{init}}$ guarantee escaping the bad region? We answer this in Section VI-D; we defer the details to that latter section and only present the final result below.

The theorem below provides an efficient sufficient condition for perfect recovery using PhaseLamp.

where $\theta^{\ast}_{\alpha}$ is the unique solution in the interval $(0,\pi/2)$ of the following equation:

Note that the sufficient condition for PhaseMax and PhaseLamp given in (17) and (37), respectively, are valid for any initial guess vector $\boldsymbol{x}_{\text{init}}$, which can depend on the sensing vectors $\{\boldsymbol{a}_{i},1\leq i\leq m\}$ and the target signal $\boldsymbol{\xi}$. It can be noticed that the PhaseLamp largely improves the performance of the PhaseMax method for dependent and independent initial guess vector $\boldsymbol{x}_{\text{init}}$.

For a better interpretation of the theorem, we have depicted the sufficient recovery condition in Figure 9. In particular, the theorem guarantees that all pairs $(\alpha,\rho_{\text{init}})$ that are above the blue dashed curve lead to perfect recovery performance of PhaseLamp. In the same figure, we also depict in red dashed line the corresponding sufficient condition of PhaseMax from Proposition 1. Clearly, these results indicate that PhaseLamp outperforms PhaseMax in the sense that it achieves perfect recovery for a larger range of input parameters $(\alpha,\rho_{\text{init}})$.

IV-B2 Independent initialization

Similar to Section III, if the initial vector $\boldsymbol{x}_{0}=\boldsymbol{x}_{\text{init}}$ is independent of the sensing vectors and of the target vector, then we can obtain sharper recovery guarantees as shown in the proposition below. For the statement of the proposition it is convenient to first define the following:

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

where $g_{\alpha}(\cdot)$ is the function defined in (26).

The sufficient condition of the theorem is depicted in blue solid line in Figure 9. Observe that it is a ”stronger” condition than that of Theorem 4 when the initialization vector is independent of the sensing vectors and of the true signal. Also, observe by comparison with the red solid line, which represents the result of Theorem 3, that PhaseLamp outperforms Phasemax. In fact, this statement is provable since the condition of Theorem 3 is not only sufficient, but also necessary.

Finally, despite the condition of Theorem 5 being only a sufficient one, the simulation results in Figure 1 suggest that it still provides a reasonably tight bound on the actual performance of PhaseLamp.

IV-C Weighted PhaseLamp

The Weighted PhaseLamp (WPhaseLamp) method is an alternative nonconvex formulation of the phase retrieval problem. Specifically, it consists of formulating the phase retrieval problem as a quadratic program

and $\omega(\cdot)$ is a preprocessing function. Note that the cost function of the optimization problem (41) is a weighted version of the PhaseLamp problem formulated in (7) where the weights depend on the sensing vectors $\{\boldsymbol{a}_{i},0\leq i\leq m\}$ and the target signal vector $\boldsymbol{\xi}$. It can also be noticed that the problem given in (41) is the spectral initialization problem where only the unit norm constraint is replaced by the linear PhaseMax constraints.

In general, the preprocessing function $\omega(\cdot)$ can have negative output values . Hence, the matrix $\boldsymbol{D}_{m}$ is indefinite in general. This means that the cost function of the problem (41) is not convex or concave in general. Write the matrix $\boldsymbol{D}_{m}$ as follows

where $\boldsymbol{D}_{m}^{(1)}$ is constructed using the negative eigenvalues of $\boldsymbol{D}_{m}$ and $\boldsymbol{D}_{m}^{(2)}$ is constructed using the negatives of the positive eigenvalues of $\boldsymbol{D}_{m}$. This means that the cost function of the Weighted PhaseLamp problem (41) can be expressed as a difference of concave functions. Therefore, one can use the convex-concave procedure to efficiently solve the Weighted PhaseLamp problem (41). Specifically, the proposed Weighted PhaseLamp algorithm consists of the following iterative scheme

The analysis of the Weighted PhaseLamp method is left for future work. Next, we provide a simulation example to compare the recovery performance of the Weighted PhaseLamp, the PhaseLamp and the PhaseMax methods. To this end, we set the signal dimension to $n=200$ and we initialize the algorithms randomly. Figure 10 plots the NMSE as a function of the oversampling ratio $\alpha$. It can be noticed that the Weighted PhaseLamp method provides a better recovery performance as compared to the PhaseLamp method and the PhaseMax method for the considered preprocessing functions. Note that the critical oversampling ratio $\alpha_{c}$ needed by the Weighted PhaseLamp method for $\omega\mathrel{\mathop{\mathchar 58\relax}}y\to y^{2}$ is around $\alpha_{c}\approx 3$. Whereas, it is around $\alpha_{c}\approx 4$ for the PhaseLamp method. Additionally note that the optimal preprocessing function introduced in outperforms the preprocessing function $\omega\mathrel{\mathop{\mathchar 58\relax}}y\to y^{2}$.

V Additional Numerical Results

In this section, we present additional simulation results and we compare the performance of polytope-optimization based methods (i.e, PhaseMax, PhaseLamp, WPhaseLamp) to other existing recovery methods in the literature; in particular, Fienup , Wirtinger Flow (WF) , Truncated amplitude flow (TAF) , PhaseLift . All algorithms are initialized using the optimal spectral initialization proposed in and the optimization problems are solved using the PhasePack . In our simulations we consider the following two cases on the measurement vectors: (1) random complex Gaussian measurements, and (2) coded diffraction patterns.

First, we consider sensing vectors that follow a circularly symmetric normal distribution, i.e., $\boldsymbol{a}_{i}\sim\mathcal{CN}(0,{\boldsymbol{I}_{n}}),1\leq i\leq m$. In Figure 11, we plot the NMSE values (average over independent problem realizations) as a function of the oversampling ratio $\alpha$. Observe that WPhaseLamp appears to outperform the rest of the recovery methods. Also, note that PhaseLamp behaves worse than the PhaseMax for small values of $\alpha$ (cf. gets stuck in the bad regime of fixed points discussed in Section IV-B), but it achieves perfect recovery earlier than the latter.

V-B Fourier Measurements

Next, we consider a type of measurements that falls under the category of coded diffraction patterns, where the measurement vectors $\boldsymbol{a}_{i}$’s are the pointwise products between the $k^{th}$ Fourier vector $\boldsymbol{f}_{k}$ and a random modulation pattern with i.i.d. symmetric Bernoulli entries $\boldsymbol{\phi}_{l}$, where $i=(k,l)$ and $1\leq k\leq n$, $1\leq l\leq\alpha$. The simulation results are presented in Figure 12. Note that the PhaseLamp and the WPhaseLamp methods provide similar recovery performance. Moreover, their performance is superior to the rest of the algorithms for $\alpha\geq 3$.

VI Technical Details: Gaussian Min-Max Inequalities

The Gordon’s Gaussian comparison inequality compares the min-max value of two doubly indexed Gaussian processes based on how their autocorrelation functions compare. The inequality is quite general (see ), but for our purposes we only need its application to the following two Gaussian processes:

Put in words: a high-probability lower bound on the AO is a high-probability lower bound on the PO. The premise is that it is often much simpler to lower bound the AO rather than the PO.

VI-A2 Convex Gaussian Min-Max Theorem (CGMT)

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$.

VI-B High-dimensional Analysis

We apply the CGMT and the GMT to characterize the asymptotic NMSE of the PhaseMax optimization in (2) as in (29) and (32) and to prove Theorem 1, respectively. To show Theorem 1, we study the asymptotic behavior of the following optimization problem

where the function $c_{d}$ is defined in (10) and its closed-form expression is given in (62) and where ${\mathcal{S}}_{\rm feas}$ is a random set defined in Section II-B.

To achieve the above goals, we start by writing the optimization problems (2), (7) and (49) in the form of a PO as in (46a), which in turn leads to a corresponding AO optimization problem. Then, we analyze the AO problem. First, define the following general optimization problem

where $p$ and $\mathcal{G}_{\boldsymbol{x}}$ are a general cost function and a general feasibility set, respectively. In this section, we are interested in the analysis of the following three cases:

In this section, we assume that $\alpha>2$. Next, the objective is to precisely analyze the problem (50) in the large system limit when ${\bf C}_{1}$ holds using the CGMT framework. Moreover, the objective is to provide a high-probability lower bound on the problem (50) when ${\bf C}_{2}$ or ${\bf C}_{3}$ holds using the GMT framework. Specifically, we show that the conditions of the CGMT (when ${\bf C}_{1}$ holds) and GMT (when ${\bf C}_{2}$/${\bf C}_{3}$ holds) are satisfied. Then, we formulate, simplify and analyze the corresponding AO.

VI-B2 Formulating the PO

Note that the GMT and CGMT assume that the feasibility sets are compact. We start our theoretical analysis by showing that the compactness assumption is guaranteed.

Assume that $\alpha>2$ and $\mathcal{K}_{n}$ is the feasibility set of the PhaseMax problem formulated in (2). Then, there exists $\tau>0$ such that

where $\tau$ is a finite constant independent of $n$.

The proof of Lemma 3 is deferred to appendix A-E. Based on Lemma 3, the optimization problem given in (50) is equivalent to the following problem with probability going to one as $n$ goes to $\infty$

where $\lambda_{n}$ is deterministic, finite and dependent on $n$. Moreover, assume that $\mathcal{V}_{n}$ denotes the set of optimal solutions of the problem (51) and $\mathcal{V}_{n}(\lambda_{n})$ denotes the set of optimal solutions of the minimization problem in (52).

If ${\bf C}_{2}$ or ${\bf C}_{3}$ holds, we have $V_{n}(\lambda_{n})\leq V_{n}$.

The proof of the above proposition is deferred to Appendix A-B. Proposition 3 shows that under condition ${\bf C}_{1}$, the precise high-dimensional analysis of the optimization problem (50) can be achieved by analyzing the problem (52). Moreover, it shows that under conditions ${\bf C}_{2}$ or ${\bf C}_{3}$, deriving a high-probability lower bound on (52) leads to a high-probability lower bound on (51). Note that the above proposition also guarantees the compactness assumption of the GMT and CGMT.

Based on Proposition 3, we proceed with analyzing the optimization problem (52). Note that the set $\mathcal{S}_{\boldsymbol{x}}$ can be rewritten as follows

Further, note that the constraint sets are convex compact and $\psi$ is convex-concave on $\mathcal{S}_{\boldsymbol{x}}\times\mathcal{S}_{\boldsymbol{u}}(n)$ if ${\bf C}_{1}$ holds, i.e. $p({x}_{1},\widetilde{\boldsymbol{x}})=-\eta_{1}x_{1}-\widetilde{\boldsymbol{\eta}}^{T}\widetilde{\boldsymbol{x}}$, where $\boldsymbol{x}=[x_{1}~{}\widetilde{\boldsymbol{x}}^{T}]^{T}$.

VI-B3 Formulating and simplifying the AO

We are now ready to formulate the corresponding AO problem as follows

Moreover, define the following optimization problem

where $p(s,r)=-s^{2}-r^{2}$ if ${\bf C}_{2}$ holds and $p(s,r)=r^{2}-\alpha c_{d}(s,r)$ if ${\bf C}_{3}$ holds. In addition, assume that ${Z}_{n}$ is the optimal objective and $\mathcal{Z}_{n}$ is the projected set of optimal solutions of the minimization problem in (VI-B3), i.e.

where ${\mathcal{S}}_{n}$ is the set of optimal solutions of the minimization problem in (VI-B3). Also, assume that $\widetilde{Z}_{1,n}$ $\widetilde{Z}_{2,n}$ are the optimal objectives and $\mathcal{\widetilde{Z}}_{1,n}$ and $\mathcal{\widetilde{Z}}_{2,n}$ are the sets of optimal $(s,r)$ of the problems (VI-B3) and (60), respectively.

The proof of the above proposition is deferred to Appendix A-C. It essentially shows that under condition ${\bf C}_{1}$, it suffices to precisely analyze the optimization problem (VI-B3) in the large system limit to determine the properties of the problem (VI-B3). Moreover, it shows that under conditions ${\bf C}_{2}$ or ${\bf C}_{3}$, deriving a high-probability lower bound on (60) leads to a high-probability lower bound on (VI-B3).

Now that we have simplified the AO to a minimization problem as in (VI-B3) and (60), we are ready to study its asymptotic behavior in the regime $m,n\rightarrow\infty,m/n\rightarrow\alpha$.

VI-C CGMT for the PhaseMax Method

In this part, we focus on the PhaseMax problem which means that we assume that the cost function $p$ is given by $p({\boldsymbol{x}})=-\eta_{1}x_{1}-\widetilde{\boldsymbol{\eta}}^{T}\widetilde{\boldsymbol{x}}$, where $\boldsymbol{x}=[x_{1}~{}\widetilde{\boldsymbol{x}}^{T}]^{T}$.

Note that Section VI-B shows that the precise high-dimensional analysis of the problem (VI-B3) can be achieved by precisely analyzing the problem (VI-B3). Hence, the main objective of this part is to study the asymptotic properties of the optimization problem (VI-B3). To this end, define the following deterministic optimization problem

where the function $c_{d}$ is defined in (10) and its closed-form expression is given by

and where the function $\Upsilon$ can be expressed as follows

The following proposition studies the asymptotic properties of the optimization problem (VI-B3) in detail. The proof of the proposition is provided in Appendix A-D.

Assume that the oversampling ratio satisfies $\alpha>2$. Let $\mathcal{V}^{\ast}_{n}$ and $V^{\ast}_{n}$ be the set of optimal solutions and the optimal objective value of the problem (VI-B3) and let $\mathcal{V}^{\ast}$ and $V^{\ast}$ be the set of optimal solutions and the optimal objective value of the deterministic problem formulated in (61). Then, we have

The above proposition shows that the set of optimal solutions and the optimal objective value of the problem (VI-B3) concentrate around the set of optimal solutions and the optimal objective value of the deterministic problem (61).

VI-C2 Solving the scalar performance optimization

In what follows, we focus on simplifying the deterministic problem (61). The following lemma, which is proved in Appendix B-C, simplifies the deterministic optimization problem (61).

The optimization problem (61) admits a unique solution in the variable $z$ which is given by

Additionally, it is equivalent to the two-dimensional problem

We call the deterministic two-dimensional optimization problem in (65) as the scalar performance optimization (SPO).Recall that the SPO in (65) is the converging limit of the problem in (VI-B3). In what follows, we solve the SPO problem for the optimal $s$ and $r$. The following lemma, which is proved in Appendix B-D, further simplifies the optimization problem (65) by showing that it has a unique optimal $r$ for any feasible variable $s$.

Fix $s$ such that $\mathinner{\!\left\lvert s\right\rvert}\leq 1$ and $\alpha>2$. Then, the following optimization problem

admits a unique global optimal solution given by

Based on P.2 in Lemma 1, the set $\mathcal{D}_{\text{feas}}$ is compact. Hence, we can always find a large enough constant $\widetilde{B}>0$ such that $s^{2}+r_{\alpha}(s)^{2}<\widetilde{B}$, for all $s$ such that $\mathinner{\!\left\lvert s\right\rvert}\leq 1$. Therefore, choosing $B$ in (65) such that $B=\widetilde{B}$ guarantees that the optimal value of $r$ in (65) is given by (67). Substituting this value back in (65) and using P.2 in Lemma 1, we can now optimize over $s$ by solving the following:

where $g_{\alpha}(s)=(r_{\alpha}(s))^{2}-\alpha~{}c_{d}(r_{\alpha}(s),s)$. A few algebraic manipulations show that the function $g_{\alpha}$ is as given in (26) and show that (68) is equivalent to (25) in the statement of Theorem 2. To show the equivalence, further note that $\eta_{1}$ and $\widetilde{\boldsymbol{\eta}}$ in (68) 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 (68) is a strictly concave program as shown in the following lemma.

For any fixed $\alpha>2$, the optimization problem formulated in (68) is strictly concave.

The proof of the above lemma is detailed in Appendix B-E. Based on Lemmas 4, 5 and 6, the deterministic optimization problem (61) has a unique global optimal solution. Based Propositions 4 and 5, the optimal objective value and the projected set of optimal solutions of the AO problem (VI-B3) concentrate around the optimal objective value and the set of optimal $(s,r)$ of the deterministic problem (61). Again, given the uniqueness of the solution of the problem (61), based on the proof of Proposition 5 and using the CGMT, the optimal objective value and the projected set of optimal solutions of the PO problem (52) concentrate around the optimal objective value and the set of optimal $(s,r)$ of the deterministic problem (61). Now, using the result stated in Proposition 3, the optimal objective value and the projected set of optimal solutions of PhaseMax (2) concentrate around the optimal objective value and the set of optimal $(s,r)$ of the deterministic problem (61).

Therefore, the optimal objective value of the PhaseMax problem (2) converges in probability to the optimal objective value of the problem (68), i.e,

Moreover, any optimal solution $\widehat{\boldsymbol{x}}$ of the PhaseMax problem (2) satisfies the following

where $s^{\ast}$ is the solution of the problem (68), $r_{\alpha}(s)$ is given in (67), $r(\widehat{\boldsymbol{x}})=\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})\right\rVert}_{2}$, and $\widehat{\boldsymbol{x}}=[s(\widehat{\boldsymbol{x}})~{}~{}\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})^{T}]^{T}$. Note that the above convergence results are valid for $\alpha>2$. This then gives us the statement of Theorem 2.

VI-C3 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 $\widehat{}\boldsymbol{x}$ of PhaseMax is, with high probability, equal to $\boldsymbol{\xi}=\boldsymbol{e}_{1}$. Mapping this to the SPO in (65), we seek conditions under which $s^{\ast}=1$ and $r^{\ast}=0$.

Assume that $\alpha>2$. From the strict concavity result in Lemma 6, perfect recovery happens if and only if the derivative of the cost function of the optimization problem (25) at $s=1$ is nonnegative. By performing a Taylor expansion of the function $s\to\sqrt{(r_{\alpha}(s))^{2}-\alpha~{}c_{d}(r_{\alpha}(s),s)}$ at $s=1$, the derivative of the cost function of the optimization problem (68) at $s=1$ can be expressed as follows

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 the condition given in (32). This then gives us the statement of Theorem 3.

VI-D Sufficient Condition for PhaseLamp

In this subsection, we focus on the PhaseLamp problem. We prove the sufficient conditions for perfect recovery of PhaseLamp stated in Theorems 4 and 5. To this end, fix the oversampling ratio $\alpha$ such that $\alpha>2$.

Note that the fixed points of the PhaseLamp algorithm are elements of the following deterministic set

The following lemma, which is proved in Appendix B-F, analyzes the system given in (73).

The system given in (73) has a unique solution. Moreover, $\text{bd}(\mathcal{D}_{\text{opt}})\not\subset\mathcal{D}_{\text{feas}}$ .

Note that the point $(0,0)$ is in the set $\text{bd}(\mathcal{D}_{\text{feas}})$ and it is also in the set $\text{bd}(\mathcal{D}_{\text{opt}})$. Moreover, observe that the target signal vector $(1,0)$ is in $\text{bd}(\mathcal{D}_{\text{feas}})$ and $\text{bd}(\mathcal{D}_{\text{opt}})$. Based on Lemma 7, the intersection between $\text{bd}({\mathcal{D}}_{\text{feas}})$ and $\text{bd}({\mathcal{D}}_{\text{opt}})$ is $\{(0,0),(\widehat{s},\widehat{r}),(1,0)\}$ where $(\widehat{s},\widehat{r})$ is the unique solution of the system in (73). Given that the solutions of (73) satisfies $r^{2}+s^{2}=s$ and $s\in(0,~{}1)$, we have $\widehat{s}\in(0,~{}1)$ and $\widehat{r}>0$. Also, Lemma 1 shows that the maximum radius of $\mathcal{D}_{\text{feas}}$ is strictly positive for $s=0$. Hence, Lemma 7 essentially shows that all the points $(s,r)$ satisfying $\widehat{s}<s<1$ are not elements of the following set $\mathcal{D}_{\text{feas}}\cap\mathcal{D}_{\text{opt}}$.

Assuming that $s=\sin(\theta)^{2}$ where $\theta\in(0,\pi/2)$, the system in (73) can be rewritten as follows

Based on Lemma 7, equation (VI-D1) has a unique solution in $(0,~{}\pi/2)$. Denote this solution by $\theta_{\alpha}^{\ast}$.

VI-D2 PhaseMax properties

The PhaseLamp method solves a PhaseMax problem as given in (34) at iteration $k+1$. Given that the target signal vectors $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are feasible for the problem (34), the optimal solution ${\boldsymbol{x}}_{k+1}$ at iteration $k+1$ satisfies the following inequality

where we express ${\boldsymbol{x}}_{k}^{T}=[\eta_{1k}~{}~{}{\widetilde{\boldsymbol{\eta}}_{k}}^{T}]$ and ${\boldsymbol{x}}_{k+1}^{T}={[s~{}~{}{{\widetilde{\boldsymbol{x}}}_{k+1}}^{T}]}$. Given that the vectors $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are both valid targets, one can assume without loss of generality that $\eta_{1k}\geq 0$, for all $k\geq 0$. Based on the Cauchy Schwarz inequality, (75) can be rewritten as follows

where $r=\mathinner{\!\left\lVert{{\widetilde{\boldsymbol{x}}}_{k+1}}\right\rVert}_{2}$. Now, define the input cosine similarity $\rho_{\text{init}}^{k}$ at iteration $k+1$ as follows

where $\boldsymbol{x}_{0}=\boldsymbol{x}_{\text{init}}$ denotes the initial guess of PhaseMax and $\rho_{\text{init}}^{0}$ is the input cosine similarity of PhaseMax, i.e. $\rho_{\text{init}}^{0}=\rho_{\text{init}}$. Note that the following equality holds for any $k\geq 0$ (recall: $\boldsymbol{\xi}=\boldsymbol{e}_{1}$.)

Hence, any optimal solution of PhaseLamp at iteration $k+1$ satisfies the following inequality $r\geq\chi_{k}(1-s)$. This implies that any optimal solution of PhaseLamp at iteration $k+1$ belongs to the following set

Now, we provide another property which guarantees that PhaseLamp escapes the bad set of stationary points and converge to the target signal vector. To this end, fix the iteration index $k\geq 0$. Based on P.3 in Lemma 1, the intersection between the boundary of the set ${\mathcal{D}}_{\rm feas}$ and the boundary of the set ${\mathcal{D}}_{\text{fp}}(\rho_{\text{init}}^{k})$ for $s\in(0,1)$ and $r>0$ satisfies

where $\chi_{k}$ is defined in (76) and it satisfies $\chi_{k}\geq 0$. Note that if $\chi_{k}=0$, the boundary of the set ${\mathcal{D}}_{\text{fp}}(\rho_{\text{init}}^{k})$ is the set of $(s,r)$ such that $r=0$. Therefore the system given in (77) has no solutions. The following lemma, which is proved in Appendix B-G, analyzes the system given in (77) in further details.

The system given in (77) has at most one solution. When (77) has a solution, the intersection between the set $\text{bd}({\mathcal{D}}_{\text{fp}}(\rho_{\text{init}}^{k}))$ and the line $s=0$ is not in the set ${\mathcal{D}}_{\rm feas}$.

VI-D3 Sufficient condition for general initialization

Now, define $\widehat{\rho}_{s}(\alpha)$ such that

Given that $\theta_{\alpha}^{\ast}\in(0~{}\pi/2)$, we have $0<\widehat{\rho}_{s}(\alpha)<1$. The following lemma shows that selecting the input cosine similarity of PhaseMax such that it is higher than $\widehat{\rho}_{s}(\alpha)$ guarantees that all the input cosine similarities of the PhaseLamp procedure are higher than $\widehat{\rho}_{s}(\alpha)$.

Select the input cosine similarity of PhaseMax $\rho_{\text{init}}$ such that $\rho_{\text{init}}>\widehat{\rho}_{s}(\alpha)$. Then, the input cosine similarity $\rho_{\text{init}}^{k}$ at iteration $k+1$ of PhaseLamp satisfy the following

The proof of the above lemma is deferred to Appendix B-H. Based on Lemma 9, we obtain $\rho_{\text{init}}^{k}>\widehat{\rho}_{s}(\alpha)$ for any $k\geq 0$. Therefore, we conclude that $0\leq s\leq\widehat{s}$ are not elements of the set $\mathcal{D}_{\text{feas}}\cap\mathcal{D}_{\text{opt}}\cap{\mathcal{D}}_{\text{fp}}(\rho_{\text{init}}^{k})$. Now, based on Lemma 7, all the points $(s,r)$ satisfying $\widehat{s}<s<1$ are not elements of the set $\mathcal{D}_{\text{feas}}\cap\mathcal{D}_{\text{opt}}$. Based on P.2 in Lemma 1, we conclude that selecting the input cosine similarity of PhaseMax $\rho_{\text{init}}$ in this way guarantees that

VI-D4 Sufficient condition for independent initialization

Note that the sufficient condition $\rho_{\text{init}}>\widehat{\rho}_{s}(\alpha)$ is valid for any initial guess vectors $\boldsymbol{x}_{\text{init}}$, which can dependent on the sensing vectors $\{\boldsymbol{a}_{i}\mathrel{\mathop{\mathchar 58\relax}}1\leq i\leq m\}$ and the target signal vector $\boldsymbol{\xi}$. Next, we focus on the case when the initial guess vector $\boldsymbol{x}_{\text{init}}$ is independent of the sensing vectors $\{\boldsymbol{a}_{i}\mathrel{\mathop{\mathchar 58\relax}}1\leq i\leq m\}$ and the target signal vector $\boldsymbol{\xi}$. To improve the above condition, we further exploit the properties of the problem (34) and PhaseMax (2) as given in the following property.

Property: The optimization problem (34) is scale invariant for any $k\geq 0$. Based on Theorem 2, the optimal solution $\widehat{\boldsymbol{x}}$ of PhaseMax satisfies the following

with $s\in[0,~{}1]$, $c_{\alpha}=1/\text{tan}\left(\pi/\alpha\right)$, $\widehat{\boldsymbol{x}}=[s(\widehat{\boldsymbol{x}})~{}~{}\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})^{T}]^{T}$ and $r(\widehat{\boldsymbol{x}})=\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}(\widehat{\boldsymbol{x}})\right\rVert}_{2}$. Based on Section VI-C, we know that $r_{\alpha}(s)$ is the unique solution to the optimization problem (66), for any $s\in[0,~{}1]$. This means that $\forall~{}s\in[0,~{}1]$, $(s,r_{\alpha}(s))\in\mathcal{D}_{\text{feas}}$.

The above property shows that it suffices to select $\rho_{\text{init}}>\rho_{s}(\alpha)$ to guarantee (81), where $\rho_{s}(\alpha)$ is determined such that the optimal solution of the following optimization problem

is $\widehat{s}_{\alpha}$ the unique solution of the following system of equations

where the function $g_{\alpha}$ is defined in (26).

Figure 13 illustrates the above sufficient condition for $\alpha=5$. Note that due to the scale invariance of the optimization problem (34), it is sufficient to select $\rho_{s}(\alpha)$ such that the solution of the PhaseMax problem in the large system limit is determined by the intersection between the equations $r_{\alpha}(s)=\sqrt{c_{\alpha}^{2}+1-s^{2}}-c_{\alpha}$ (cyan curve) and $r={\sqrt{1-\widehat{\rho}_{s}(\alpha)^{2}}/\widehat{\rho}_{s}(\alpha)}~{}s$ (magenta curve).

Note that the unique solution $\widehat{s}_{\alpha}$ of (84) can be expressed as follows

where $a_{\alpha}=\frac{\sqrt{1-\widehat{\rho}_{s}(\alpha)^{2}}}{\widehat{\rho}_{s}(\alpha)}$. Given that $\widehat{\rho}_{s}(\alpha)$ is selected such that (78) is satisfied, we have $\widehat{\rho}_{s}(\alpha)^{2}=\sin(\theta_{\alpha}^{\ast})^{2}$ where $\theta_{\alpha}^{\ast}$ is the unique solution of (VI-D1). This means that $a_{\alpha}=\tan(\theta_{\alpha}^{\ast})^{-1}$ and $\widehat{s}_{\alpha}$ can be rewritten as follows

To ensure that $\widehat{s}_{\alpha}$ is the optimal solution of the optimization problem (83), the first derivative of the cost function of the problem (83) should be zero at $\widehat{s}_{\alpha}$. Note that the first derivative of the cost function of problem (83) can be expressed as

This means that the sufficient input cosine similarity $\rho_{s}(\alpha)$ satisfies the following

VI-D5 Convergence analysis

Now, assume that the input cosine similarity satisfies $\rho_{\text{init}}>\widehat{\rho}_{s}(\alpha)$ for general initial guess and it satisfies $\rho_{\text{init}}>{\rho}_{s}(\alpha)$ for independent initial guess. This means that (81) is satisfied. Based on Lemma 1, an input cosine similarity selected in this way ensures that for any $\epsilon>0$ there exists $\delta>0$ such that

for any $j\geq 0$, where $\lim_{j\to 0}\delta_{\epsilon_{j}}=0$ and $\mathcal{F}_{\text{lamp}}$ denotes the set of fixed points of the PhaseLamp algorithm. This implies that the set $\mathcal{F}_{\text{lamp}}$ converges to the set $\{\boldsymbol{\xi}\}$ in the sense that $\sup_{\widehat{\boldsymbol{x}}\in\mathcal{F}_{\text{lamp}}}\mathinner{\!\left\lVert\widehat{\boldsymbol{x}}-\boldsymbol{\xi}\right\rVert}_{2}$ converges to zero in probability. This then gives us the statement of Theorems 4 and 5.

VII 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 PhaseMax method.

Appendix A Probabilistic Analysis

and where $\rho\mathrel{\mathop{\mathchar 58\relax}}x\to\max(x,0)$. Define the following problem

Next, we study the asymptotic properties of the problem (A-A). Specifically, we study the convergence properties (with growing $n$) of the formulation given in (A-A). Based on the proof of Proposition 5 provided in Appendix A-D, we have the following convergence

Using GMT and based on Section VI-B, it holds that

Now, define the deterministic set ${\mathcal{D}}_{\rm feas}^{\epsilon}$ as follows

Therefore, the convergence result in (98) is equivalent to the following

A-B Proof of Proposition 3

First, we appropriately write the optimization problem in (51) as a min-max program. Start with the following equivalent formulation:

where the function $x\to\mathds{1}\{x\}$ is defined as follows: $\mathds{1}\{x\}=0$ if $x\leq 0$ and $\mathds{1}\{x\}=+\infty$ if $x>0$. Therefore, (101) is equivalent to the following optimization problem

The GMT and the CGMT assumes that the feasibility sets of the optimization variables $\boldsymbol{x}$ and $\boldsymbol{u}$ are compact. Clearly, this assumption is not satisfied by the min-max problem (A-B) since the feasibility set of the variable $\boldsymbol{u}$ is not compact. Case 1: Assume that ${\bf C}_{2}$ or ${\bf C}_{3}$ holds. It can be noticed that the optimal objective of (52) is smaller than the optimal objective of (A-B) with probability one, i.e. $V_{n}(\lambda_{n})\leq V_{n}$. Case 2: Assume that ${\bf C}_{1}$ holds. Define the following optimization problem

Let $\widetilde{\mathcal{K}}_{n}$ be the feasibility set of the optimization problem (104). Clearly, the feasibility set $\widetilde{\mathcal{K}}_{n}$ is a polytope with nonempty extreme point set. Moreover, the cost function of the problem in (104) is lower bounded in the feasibility set $\widetilde{\mathcal{K}}_{n}$. Then, using the result in [32, Corollary 32.3.4], the optimal objective value $C_{n}(\lambda_{n})$ is achieved at one of the vertices of the polytope $\widetilde{\mathcal{K}}_{n}$. Define the set $\widetilde{\mathcal{E}}$ as follows

where the set $\mathcal{E}(\widetilde{\mathcal{K}}_{n})$ denotes the set of all extreme points of the polytope $\widetilde{\mathcal{K}}_{n}$. Since the polytope $\widetilde{\mathcal{K}}_{n}$ has a finite number of extreme points, the set $\widetilde{\mathcal{E}}$ has a finite cardinality.

Assume that $\lambda_{n}\geq\mathinner{\!\left\lVert\boldsymbol{x}_{\text{init}}\right\rVert}_{2}w_{n}\sqrt{n}/\zeta_{n}$ where $\zeta_{n}$ is defined as follows

Case 2.b: Assume that the set $\widetilde{\mathcal{E}}$ is nonempty. Then, for any extreme point $(\boldsymbol{x}^{T},\boldsymbol{z}^{T})$ of the polytope $\widetilde{\mathcal{K}}_{n}$ which belongs to the set $\widetilde{\mathcal{E}}$, we have

Finally, consider the event $E_{n}=\{\boldsymbol{A}\mathrel{\mathop{\mathchar 58\relax}}~{}\mathinner{\!\left\lVert\mathcal{K}_{n}\right\rVert}_{2}\leq\tau\}$, then, we have the following

where the function $\rho\mathrel{\mathop{\mathchar 58\relax}}x\to\max(x,0)$. Denote by $\widehat{V}_{n}$ the cost function of the problem (107) and $\widetilde{V}_{n}$ the cost function of the problem (108). Further, assume that $\widehat{\boldsymbol{x}}_{n}$ is an optimal solution of the problem (107) and $\widetilde{\boldsymbol{x}}_{n}$ is an optimal solution of the problem (108). It is clear that $\widehat{V}_{n}(\widehat{\boldsymbol{x}}_{n})=\widetilde{V}_{n}(\widehat{\boldsymbol{x}}_{n})$ and also

Since the all zero vector is in the polytope $\mathcal{K}_{n}$, $V_{n}$ is finite with probability one. Moreover, since

we obtain the following convergence result

where $\mathcal{V}_{n}$ denotes the set of optimal solutions of the problem (107) and $\mathcal{V}_{n}(\lambda_{n})$ denotes the set of optimal solutions of the problem (108).

The above two cases give us the statement in Proposition 3.

A-C Proof of Proposition 4

It can be noticed that the optimization problem (VI-B3) can be rewritten as follows:

Next, observe that if we fix $\mathinner{\!\left\lvert\boldsymbol{u}\right\rvert}$, then the optimal $\boldsymbol{u}$ satisfies $\text{sign}(\boldsymbol{u})=\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

Define the following optimization problem

Next, we show that the analysis of the optimization problem

can be achieved by analyzing the following problem

where $\boldsymbol{y}_{n}\in\mathcal{X}_{2}^{\ast}(\lambda_{n})$ and $\widetilde{\boldsymbol{y}}_{n}\in\mathcal{X}_{2}^{\ast}(\sqrt{n}\lambda_{n})$. This implies that

Now, for $n_{j}\geq\max(n_{j_{0}},n_{0})$, we have

where this is true for any $\epsilon>0$. Taking $\epsilon\leq\delta$ gives a contradiction. This implies that

In what follows, we analyze the problem (A-C) where the sequence $\lambda_{n}<\infty$ satisfies $\lambda_{n}\overset{n\to\infty}{\longrightarrow}\infty$. In the optimization problem (A-C), one can fix the norm of $\boldsymbol{u}$ and optimize over its direction. This leads to the following optimization problem

where the function $h$ is defined in (58). Therefore, (A-C) is equivalent to the following problem

where the function $\rho\mathrel{\mathop{\mathchar 58\relax}}x\to\max(x,0)$. Now, we distinguish between two cases: Case 1: Assume that ${\bf C}_{1}$ holds, i.e. $p({x}_{1},\widetilde{\boldsymbol{x}})=-\eta_{1}x_{1}-\widetilde{\boldsymbol{\eta}}^{T}\widetilde{\boldsymbol{x}}$. 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}}$. First, fix $\mathinner{\!\left\lVert\widetilde{\boldsymbol{x}}\right\rVert}_{2}=r$ such that $s^{2}+r^{2}\leq B,~{}r\geq 0$ and fix $z=\widetilde{\boldsymbol{\eta}}^{T}\widetilde{\boldsymbol{x}}/\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}$ and solve the following optimization problem

To solve the optimization problem (129), we write $\boldsymbol{h}$ as follows $\boldsymbol{h}=\left(\frac{\boldsymbol{h}^{T}\widetilde{\boldsymbol{\eta}}}{\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}}\right)\frac{\widetilde{\boldsymbol{\eta}}}{\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}}+\left(\boldsymbol{w}^{T}\boldsymbol{h}\right)\boldsymbol{w}$, where $\frac{\widetilde{\boldsymbol{\eta}}}{\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}}$ and $\boldsymbol{w}$ form an orthonormal basis for the two dimensional subspace spanned by $\widetilde{\boldsymbol{\eta}}$ and $\boldsymbol{h}$. Thus, (129) becomes

It is clear that the optimal $\widetilde{\boldsymbol{x}}$ should be in the span of $\widetilde{\boldsymbol{\eta}}$ and $\boldsymbol{h}$ which means that

Therefore, the optimal objective value of the problem (130) can be expressed as follows

where $\mathinner{\!\left\lvert z\right\rvert}\leq r$. Assume that $\widehat{\boldsymbol{\eta}}=\frac{\widetilde{\boldsymbol{\eta}}}{\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}}$, then, the optimization problem reduces to the following problem

where the function $c_{n}$ is defined in (57), $\widetilde{\lambda}_{n}=\sqrt{m(n)}\lambda_{n}$, and where $p(s,r)=-s^{2}-r^{2}$ if ${\bf C}_{2}$ holds and $p(s,r)=r^{2}-\alpha c_{d}(s,r)$ if ${\bf C}_{3}$ holds.

Note that (118) and (126) hold for all cases ${\bf C}_{1}$, ${\bf C}_{2}$ and ${\bf C}_{3}$. This, then, leads us to the statement in Proposition 4.

A-D Proof of Proposition 5

Assume that the oversampling ratio satisfies $\alpha>2$. We show Proposition 5 in three steps. The first two steps study the asymptotic properties of the random function $c_{n}$. Then, these properties are used to prove Proposition 5 in the final step. To this end, consider the random function $f_{n}$ defined as follows

Step 1: We start by showing that the functions $f_{n}$ and $\widetilde{c}_{n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to-c_{n}(s,r)/\sqrt{m(n)}$ have the same pointwise limit. Moreover, the function $f_{n}$ converges pointwise to the function $(s,r)\to\sqrt{c_{d}(s,r)}$, where the function $c_{d}$ is defined in (10) and its closed-form expression is given in (62). To prove the above property, fix $s$ and $r$ such that $(s,r)\in\mathcal{S}$. Using the weak law of large number (WLLN), we have

Now, fix $\epsilon>0$ and consider the probability event $E_{n}=\{(\boldsymbol{q},\boldsymbol{g})\mathrel{\mathop{\mathchar 58\relax}}~{}\mathinner{\!\left\lvert\widetilde{c}_{n}(s,r)-\sqrt{c_{d}(s,r)}\right\rvert}>\epsilon\}$. Then, we have

where $A$ is the event $\{(\boldsymbol{q},\boldsymbol{g})\mathrel{\mathop{\mathchar 58\relax}}~{}\min(\mathinner{\!\left\lvert\boldsymbol{q}\right\rvert}-\mathinner{\!\left\lvert r\boldsymbol{g}+s\boldsymbol{q}\right\rvert})\leq 0\}$. The probability of the event $A^{c}$ is given by

Equation (136) can be rewritten as follows

Given that $\{q_{i}\}_{i=1}^{m(n)}$ are i.i.d. standard Gaussian random variable, the above equation can be rewritten as follows

where $\Phi$ denotes the cumulative distribution function of the standard normal random variable. We know that $\epsilon>0$ and $\mathinner{\!\left\lvert s\right\rvert}<1$, hence, we obtain the following inequality $2\Phi\left(-{\sqrt{m(n)}{\epsilon}}/{(1-\mathinner{\!\left\lvert s\right\rvert})}\right)<1$, for all $m(n)>0$. This leads to the following

Hence, we can conclude that $\widetilde{c}_{n}(s,r)$ converges in probability to $\sqrt{c_{d}(s,r)}$ for any $s$ and $r$ such that $r=0$ and $\mathinner{\!\left\lvert s\right\rvert}<1$. Case 2: if $r=0$ and $\mathinner{\!\left\lvert s\right\rvert}\geq 1$ or $r\neq 0$. In this case, note that

We can conclude that $-c_{n}(s,r)/\sqrt{m(n)}$ converges in probability to $\sqrt{c_{d}(s,r)}$ in this case.

Based on Case 1 and Case 2, the functions $f_{n}$ and $\widetilde{c}_{n}$ have the same pointwise limit which is the function $(s,r)\to\sqrt{c_{d}(s,r)}$. Moreover, the function $c_{d}$ is given by

It can be checked that the function $c_{d}$ is as given in (62).

Step 2: The first step mainly shows that the functions $f_{n}$ and $\widetilde{c}_{n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to-c_{n}(s,r)/\sqrt{m(n)}$ have the same pointwise limit. In the second step, we show that they also converge uniformly in probability to the same function. First, assume that the function $f_{n}$ converges uniformly to some function $f^{\ast}$ and fix $\epsilon>0$. This means that

Consider the following three functions $f_{1,n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to\widetilde{c}_{n}(s,r)-f^{\ast}(s,r)$, $f_{2,n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to\widetilde{c}_{n}(s,r)-f_{n}(s,r)$ and $f_{3,n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to f_{n}(s,r)-f^{\ast}(s,r)$. It is clear that

Consider the following probability events

Consider the following two functions $g_{1,n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to\min(\mathinner{\!\left\lvert\boldsymbol{q}\right\rvert}-\mathinner{\!\left\lvert r\boldsymbol{g}+s\boldsymbol{q}\right\rvert})$ and $g_{2,n}\mathrel{\mathop{\mathchar 58\relax}}(s,r)\to\mathinner{\!\left\lVert(\mathinner{\!\left\lvert\boldsymbol{q}\right\rvert}-\mathinner{\!\left\lvert r\boldsymbol{g}+s\boldsymbol{q}\right\rvert})\wedge\mathbf{0}\right\rVert}_{2}$. The function $f_{2,n}$ can be expressed as follows

Define the set $\mathcal{C}$ as $\mathcal{C}=\{(s,r)\in\mathcal{S}\mathrel{\mathop{\mathchar 58\relax}}~{}g_{1,n}(s,r)>0\}$ and $P$ as the probability of the event $\{(\boldsymbol{q},\boldsymbol{g})\mathrel{\mathop{\mathchar 58\relax}}\sup\limits_{\begin{subarray}{c}(s,r)\in\mathcal{S}\end{subarray}}\mathinner{\!\left\lvert f_{2,n}(s,r)\right\rvert}>\frac{\epsilon}{2}\}$. Then, we have

Given that $\{q_{i}\}_{i=1}^{m(n)}$ are i.i.d. standard Gaussian random variable, we get

where $\Phi$ denotes the cumulative distribution function of the standard normal random variable. Since $\epsilon>0$, we obtain the following inequality $2\Phi\left(-\sqrt{{m(n)}}\frac{\epsilon}{2}\right)<1,~{}\forall m(n)>0$, which means that

which means that the function $\widetilde{c}_{n}$ converges uniformly to the function $f^{\ast}$. Now, if we repeat the above steps with $f_{n}$ replaced by $\widetilde{c}_{n}$, we obtain the second direction.

where $\widehat{\boldsymbol{\eta}}=\widetilde{\boldsymbol{\eta}}/\mathinner{\!\left\lVert\widetilde{\boldsymbol{\eta}}\right\rVert}_{2}$ and define the deterministic function $Q$ on the set $\mathcal{D}$ as follows

Fix $(s,r,z)$ in the set $\mathcal{D}$. Given that $\boldsymbol{h}$ is independent of $\widehat{\boldsymbol{\eta}}$, we have $\boldsymbol{h}^{T}\widehat{\boldsymbol{\eta}}/\sqrt{m(n)}\xrightarrow[]{n\to\infty}0$. Furthermore, using the WLLN, we have $\mathinner{\!\left\lVert\boldsymbol{h}\right\rVert}_{2}^{2}/{m(n)}\xrightarrow[]{n\to\infty}\frac{1}{{\alpha}}$. Therefore, based on the first step, the function $Q_{n}$ converges pointwise to the function $Q$.