Phase retrieval with polarization

Introduction

These deficiencies have prompted two important lines of research in phase retrieval:

For which measurement designs $\Phi$ is $[x]\mapsto|\Phi^{*}x|^{2}$ injective?

For which injective designs can $[x]$ be reconstructed stably and efficiently?

This leads one to attempt provably stable and efficient reconstruction from measurements of the form (1) with particular ensembles $\Phi$. Until recently, this was only known to be possible in cases where $N=\Omega(M^{2})$ . By contrast, the state of the art comes from Candès, Strohmer and Voroninski , who use semidefinite programming to stably reconstruct from $N=\mathcal{O}(M\log M)$ Gaussian-random measurements. There is other work along this vein which also uses semidefinite programming and provides related guarantees. Typically, semidefinite programs are solved via interior point methods. Since these methods are computationally expensive, in practice, one is inclined to instead use faster numerical methods, but these lack performance guarantees.

Returning to measurements of the form (1), this paper combines ideas from both state-of-the-art theory and state-of-the-art practice by proposing an exchange of sorts: If you already have $\mathcal{O}(M\log M)$ Gaussian-random measurements vectors (as prescribed in ), then we offer a faster reconstruction method with a stable performance guarantee, but at the price of $\mathcal{O}(M\log M)$ additional (non-adaptive) measurements. These new measurement vectors are interferometry-inspired combinations of the originals, and the computational speedups gained in reconstruction come from our use of different spectral methods. While the ideas in this paper can be applied for phase retrieval of 2-D images, we focus on the 1-D case for simplicity. Also, note that the sequel leverages the techniques of this paper to construct masked Fourier measurements, thereby mimicking the illumination methodology of ; we suspect that these ideas can be similarly leveraged to tackle a wide variety of practical instances of the phase retrieval problem. To help motivate our measurement design and phase retrieval procedure, we start in the next section by considering the simpler, noiseless case. In this case, the success of our method follows from a neat trick involving the polarization identity along with some well-known results in the theory of expander graphs. In Section 3, we modify the method to obtain provable stability in the noisy case; here, we exploit some recent developments in spectral graph theory. Our results are then corroborated by simulations in Section 4. We give concluding remarks in Section 5, and provide the more technical proofs in the appendix.

The noiseless case

Having established the utility of the relative phase between coefficients, we now seek some method of extracting this information. To this end, we turn to a special version of the polarization identity:

We start by expanding the right-hand side of (4):

Finally, we apply the following easy-to-verify identities:

Thus, if in addition to $\Phi_{V}$ we measure with $\{\varphi_{i}+\zeta^{k}\varphi_{j}\}_{k=0}^{2}$, we can use (5) to determine $\overline{\langle x,\varphi_{i}\rangle}\langle x,\varphi_{j}\rangle$ and then normalize to get the relative phase:

In pursuit of $\mathcal{O}(M)$ measurements, take some simple graph $G=(V,E)$, arbitrarily assign a direction to each edge, and only take measurements with $\Phi_{V}$ and $\Phi_{E}:=\bigcup_{(i,j)\in E}\{\varphi_{i}+\zeta^{k}\varphi_{j}\}_{k=0}^{2}$. To recover $[x]$, we again arbitrarily assign some nonzero vertex measurement to have positive phase, and then we propagate relative phase information along the edges by multiplication to determine the phase of the other vertex measurements relative to the original vertex measurement:

However, if $x$ is orthogonal to a given vertex vector, then that measurement is zero, and so relative phase information cannot propagate through the corresponding vertex; indeed, such orthogonality has the effect of removing the vertex from the graph, and for some graphs, this will prevent recovery. For example, if $G$ is a star, then $x$ could be orthogonal to the vector corresponding to the internal vertex, whose removal would render the remaining graph edgeless. That said, we should select $\Phi_{V}$ and $G$ so as to minimize the impact of orthogonality with vertex vectors.

First, we can take $\Phi_{V}$ to be full spark, that is, $\Phi_{V}$ has the property that every subcollection of $M$ vectors spans. Full spark frames appear in a wide variety of applications. Explicit deterministic constructions of them are given in . For example, we can select the first $M$ rows of the $|V|\times|V|$ discrete Fourier transform matrix, and take $\Phi_{V}$ to be the columns of the resulting $M\times|V|$ matrix; in this case, the fact that $\Phi_{V}$ is full spark follows from the Vandermonde determinant formula. In our application, $\Phi_{V}$ being full spark will be useful for two reasons. First, this implies that $x\neq 0$ is orthogonal to at most $M-1$ members of $\Phi_{V}$, thereby limiting the extent of $x$’s damage to our graph. Additionally, $\Phi_{V}$ being full spark frees us from requiring the graph to be connected after the removal of vertices; indeed, any remaining component of size $M$ or more will correspond to a subcollection of $\Phi_{V}$ that spans, meaning it has a dual frame to reconstruct with. It remains to find a graph of $\mathcal{O}(M)$ vertices and edges that maintains a size-$M$ component after the removal of any $M-1$ vertices.

To this end, we consider a well-studied family of sparse graphs known as expander graphs. We choose these graphs for their notably strong connectivity properties. There is a combinatorial definition of expander graphs, but we will focus on the spectral definition. Given a $d$-regular graph $G$ of $n$ vertices, consider its adjacency matrix $A$, and define the Laplacian to be $L:=I-\frac{1}{d}A$; if $G$ were not regular, we would consider the diagonal matrix $D$ of vertex degrees and define the Laplacian to be $L:=I-D^{-1/2}AD^{-1/2}$. This is often called the normalized Laplacian in the literature, but we make no distinction here. We are particularly interested in the eigenvalues of the Laplacian: $0=\lambda_{1}\leq\cdots\leq\lambda_{n}$. The second eigenvalue $\lambda_{2}$ of the Laplacian is called the spectral gap of the graph, and as we shall see, this value is particularly useful in evaluating the graph’s connectivity. We say $G$ has expansion $\lambda$ if $\{\lambda_{2},\ldots,\lambda_{n}\}\subseteq[1-\lambda,1+\lambda]$; note that since $1-\lambda\leq\lambda_{2}$, small expansion implies large spectral gap. Furthermore, a family of $d$-regular graphs $\{G_{i}\}_{i=1}^{\infty}$ is a spectral expander family if there exists $c<1$ such that every $G_{i}$ has expansion $\lambda(G_{i})\leq c$. Since $d$ is constant over an expander family, expanders with many vertices have particularly few edges. There are many results which describe the connectivity of expanders, but the following is particularly relevant to our application:

Consider a $d$-regular graph $G$ of $n$ vertices with spectral gap $\lambda_{2}$. For all $\varepsilon\leq\frac{\lambda_{2}}{6}$, removing any $\varepsilon dn$ edges from $G$ results in a connected component of size $\geq(1-\frac{2\varepsilon}{\lambda_{2}})n$.

Note that removing $\varepsilon n$ vertices from a $d$-regular graph necessarily removes $\leq\varepsilon dn$ edges, and so this lemma directly applies. For our application, we want to guarantee that the removal of any $M-1$ vertices maintains a size-$M$ component. To do this, we will ensure both (i) $M-1\leq\varepsilon n$ and (ii) $M-1<(1-\frac{2\varepsilon}{\lambda_{2}})n$, and then invoke the above lemma. Note that since $n\geq M\geq 2$,

where the last inequality is a rearrangement of $\varepsilon\leq\frac{\lambda_{2}}{6}$. Thus $\varepsilon n<(1-\frac{2\varepsilon}{\lambda_{2}})n$, meaning (i) implies (ii), and so it suffices to have $M\leq\varepsilon n+1$. Overall, we use the following criteria to pick our expander graph: Given the signal dimension $M$, use a $d$-regular graph $G=(V,E)$ of $n$ vertices with spectral gap $\lambda_{2}$ such that $M\leq(\frac{\lambda_{2}}{6})n+1$. Then by the previous discussion, the total number of measurements is $N=|V|+3|E|=(\frac{3}{2}d+1)n$. If we think of the degree $d$ as being fixed, then the number of vertices $n$ in the graph is proportional to the total number of measurements $N$ (this is the key distinction from the previous complete-graph case).

Recall that we seek $N=\mathcal{O}(M)$ measurements. To minimize the redundancy $\frac{N}{M}$ for a fixed degree $d$, we would like a maximal spectral gap $\lambda_{2}$, and it suffices to seek minimal spectral expansion $\lambda$. Spectral graph families known as Ramanujan graphs are asymptotically optimal in this sense; taking $\mathcal{G}_{n}^{d}$ to be the set of connected $d$-regular graphs with $\geq n$ vertices, Alon and Boppana (see ) showed that for any fixed $d$,

while Ramanujan graphs are defined to have spectral expansion $\leq\frac{2\sqrt{d-1}}{d}$. To date, Ramanujan graphs have only been constructed for certain values of $d$. One important construction was given by Lubotzky, Phillips, and Sarnak , which produces a Ramanujan family whenever $d-1\equiv 1\bmod 4$ is prime. Among these graphs, we get the smallest redundancy $\frac{N}{M}$ when $M=\lfloor(1-\frac{2\sqrt{d-1}}{d})\frac{n}{6}+1\rfloor$ and $d=6$:

Thus, in such cases, our techniques allow for phase retrieval with only $N\leq 236M$ measurements. However, the number of vertices in each Ramanujan graph from is of the form $q(q^{2}-1)$ or $\frac{q(q^{2}-1)}{2}$, where $q\equiv 1\bmod 4$ is prime, and so any bound on redundancy $\frac{N}{M}$ using these graphs will only be valid for particular values of $M$.

In order to get $N=\mathcal{O}(M)$ in general, we use the fact that random graphs are nearly Ramanujan with high probability. In particular, for every $\varepsilon>0$ and even $d$, a random $d$-regular graph has spectral expansion $\lambda\leq\frac{2\sqrt{d-1}+\varepsilon}{d}$ with high probability as $n\rightarrow\infty$ . Thus, picking $\varepsilon$ and $d$ to satisfy $\frac{2\sqrt{d-1}+\varepsilon}{d}<1$, we may take $M=\lfloor(1-\frac{2\sqrt{d-1}+\varepsilon}{d})\frac{n}{6}+1\rfloor$ to get

and this choice will satisfy $M\leq(\frac{\lambda_{2}}{6})n+1$ with high probability. To see how small this redundancy is, note that taking $\varepsilon=0.1$ and $d=8$ gives $N\leq 240M$. While the desired expansion properties of a random graph are only present with high probability, estimating the spectral gap is inexpensive, and so it is computationally feasible to verify whether a randomly drawn graph is good enough. Moreover, $n$ can be any sufficiently large integer, and so the above bound is valid for all sufficiently large $M$, i.e., our procedure can perform phase retrieval with $N=\mathcal{O}(M)$ measurements in general.

Combining this with the above discussion, we have the following measurement design and phase retrieval procedure:

Fix $d>2$ even and $\varepsilon\in(0,d-2\sqrt{d-1})$.

Given $M$, pick some $d$-regular graph $G=(V,E)$ with spectral gap $\lambda_{2}\geq\lambda^{\prime}:=1-\frac{2\sqrt{d-1}+\varepsilon}{d}$ and $|V|=\lceil\frac{6}{\lambda^{\prime}}(M-1)\rceil$, and arbitrarily direct the edges.

Phase Retrieval Procedure A (noiseless case)

Given $\{|\langle x,\varphi\rangle|^{2}\}_{\varphi\in\Phi}$, delete the vertices $i\in V$ with $|\langle x,\varphi_{i}\rangle|^{2}=0$.

In the remaining induced subgraph, find a connected component of $\geq M$ vertices $V^{\prime}$.

Pick a vertex in $V^{\prime}$ to have positive phase and propagate/multiply relative phases (7), which are calculated by normalizing (5), see (6).

Having $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime}}$ up to a global phase factor, find the least-squares estimate of $[x]$ by applying the Moore-Penrose pseudoinverse of $\{\varphi_{i}\}_{i\in V^{\prime}}$, see (3).

Note that this phase retrieval procedure is particularly fast. Indeed, if we use $E\subseteq V^{2}$ to store $G$, then we can delete vertices $i\in V$ with $|\langle x,\varphi_{i}\rangle|^{2}=0$ by deleting the edges for which (5) is zero, which takes $\mathcal{O}(|E|)$ time. Next, if the members of $E$ are ordered lexicographically, the remaining subgraph can be easily partitioned into connected components in $\mathcal{O}(|E|)$ time by collecting edges with common vertices, and then propagating relative phase in the largest component is performed in $\mathcal{O}(|E|)$ time using a depth- or breadth-first search. Overall, we only use $\mathcal{O}(M)$ time before the final least-squares step of the phase retrieval procedure, which happens to be the bottleneck, depending on the subcollection $\Phi_{V^{\prime}}$. In general, we can find the least-squares estimate in $\mathcal{O}(M^{3})$ time using Gaussian elimination, but if $\Phi_{V^{\prime}}$ has special structure (e.g., it is a submatrix of the discrete Fourier transform matrix), then one might exploit that structure to gain speedups (e.g., use the fast Fourier transform in conjunction with an iterative method). Regardless, our procedure reduces the nonlinear phase retrieval problem to the much simpler problem of solving an overdetermined linear system.

While this measurement design and phase retrieval procedure is particularly efficient, it certainly lacks stability. Perhaps most notably, we have not imposed anything on $\Phi_{V}$ that guarantees stability with inverting $\Phi_{V^{\prime}}$; indeed, we have merely enforced linear independence between vectors, while stability will require well-conditioning. Another noteworthy source of instability is our method of phase propagation, which naturally accumulates error; it would be better if the relative phases were combined using a more democratic process that encourages noise cancellation. In the next section, we will address these concerns (and others) and modify our procedure accordingly; the revised procedure will be stable, but at the price of a log factor in the number of measurements: $N=\mathcal{O}(M\log M)$. As we mention in the concluding remarks, we do not think this log factor is necessary, but we leave this pursuit for future work.

The noisy case

In this section, we consider a noise-robust version of the measurement design and phase retrieval procedure of the previous section. In the end, the measurement design will be nearly identical: vertex measurements will be independent complex Gaussian vectors (thereby being full spark with probability 1), and the edge measurements will be the same sort of linear combinations of vertex measurements. Our use of randomness in this version will enable the vertex measurements to simultaneously satisfy two important conditions with high probability: projective uniformity with noise and numerical erasure robustness. Before defining these conditions, we motivate them by considering a noisy version of our phase retrieval procedure.

Recall that our noiseless procedure starts by removing the vertices $i\in V$ for which $|\langle x,\varphi_{i}\rangle|^{2}=0$. Indeed, since we plan to propagate relative phase information along edges, these -vertices are of no use, as relative phase with these vertices is not well defined. Since we calculate relative phase by normalizing (5), we see that relative phase is sensitive to perturbations when (5) is small, meaning either $\langle x,\varphi_{i}\rangle$ or $\langle x,\varphi_{j}\rangle$ is small. As such, while -vertices provide no relative phase information in the noiseless case, small vertices provide unreliable information in the noisy case, and so we wish to remove them accordingly (alternatively, one might use weights according to one’s confidence in the information, but we decided to use hard thresholds to simplify the analysis). However, we also want to ensure that there are only a few small vertices. In the noiseless case, we limit the number of -vertices by using a full spark frame; in the noisy case, we make use of a new concept we call projective uniformity:

We now explain why only reliable pieces of relative phase information will remain after running the above algorithm, provided $\Phi_{V}$ has sufficient projective uniformity. The main idea is captured in the following:

and so $\cos\theta>0$, i.e., $\theta\in[0,\frac{\pi}{2})$. Finally, by the concavity of $\sin(\cdot)$ and then (8), we conclude that

After applying Algorithm 1, our graph will have slightly fewer vertices, but the remaining edges will correspond to reliable pieces of relative phase information. Recall that we plan to use this information on the edges to determine phases for the vertices, and we want to do this in a stable way. To understand when this is even possible, we first consider a few simple scenarios. Suppose that after removing vertices with Algorithm 1, the graph has a vertex of degree 0. Then we have no information about the phase of this vertex, and it should be removed accordingly. For a less extreme scenario, suppose the vertex has degree 1. Then any noise in the corresponding edge measurement would be passed directly to the vertex, which inherently lacks stability compared to the noise cancellation that would come with more edges. More generally, if the graph has a cut vertex (e.g., the neighbor of a degree-1 vertex), then we would need to rely on the correctness of this lone vertex to ensure consistency between the parts of the graph it connects—this scenario is also rather unstable. After considering these examples, it makes intuitive sense that stability necessitates a high level of connectivity in the graph, regardless of the algorithm used to extrapolate the vertex phases.

As such, we seek to remove a small proportion of vertices so that the remaining graph is very connected, i.e., has large spectral gap. To do this, we will iteratively remove sets of vertices that are poorly connected to the rest of the graph. These sets will be identified using spectral clustering (Algorithm 2), a process which is strongly motivated by an inequality in Riemannian geometry by Cheeger and which has performance guarantees originating with Alon . The main idea of spectral clustering follows the intuition that a random walk on a graph tends to be trapped in sections of the graph which have few connections to the rest of the vertices (this intuition is made more explicit in ). Moreover, the second eigenvector of the corresponding stochastic matrix tends to identify these sections.

In the appendix, we show that for a particular choice of threshold $\tau$, Algorithm 3 recovers a level of connectivity that may have been lost when pruning for reliability in Algorithm 1, and it does so by removing only a small proportion of the vertices.

At this point, we have pruned our graph so that the measured relative phases are reliable and the vertex phases can be stably reconstructed. Now we seek an efficient method to reconstruct these vertex phases from the measured relative phases. Before devising such a method, we first organize the information we have into a matrix. Given the graph output $G^{\prime}=(V^{\prime},E^{\prime})$ of Algorithm 3, we take $A_{1}$ to be a $|V^{\prime}|\times|V^{\prime}|$ weighted adjacency matrix of $G^{\prime}$. Specifically, for each $\{i,j\}\in E^{\prime}$, let $\varepsilon_{ij}$ denote the effective noise in the estimate of $\overline{\langle x,\varphi_{i}\rangle}\langle x,\varphi_{j}\rangle$ using (5), and normalize this noisy estimate to get

Otherwise when $\{i,j\}\not\in E^{\prime}$, take $A_{1}[i,j]=0$. Unlike the noiseless case, here, we account for both directions $(i,j)$ and $(j,i)$ whenever $\{i,j\}\in E^{\prime}$, with the understanding that $\varepsilon_{ji}=\overline{\varepsilon_{ij}}$; this will simplify our analysis since this makes $A_{1}$ self-adjoint. Considering $A_{1}[i,j]$ is an approximation of the relative phase $\omega_{i}^{-1}\omega_{j}$, it seems reasonable to extrapolate the vertex phases $\omega:=\{\omega_{i}\}_{i\in V^{\prime}}$ from $A_{1}$ by minimizing the following quantity:

where $D$ is the diagonal matrix of vertex degrees. Dividing by $\operatorname{vol}(G^{\prime})=\omega^{*}D\omega$, which does not vary with $\omega$, this is equivalent to minimizing

To be clear, the right-hand side above is the first eigenvalue of $L_{1}:=I-D^{-1/2}A_{1}D^{-1/2}$, which we call the connection Laplacian; note that this bears some resemblance to the Laplacian defined in the previous section. In minimizing the above quantity, it makes sense to consider the eigenvector $u$ corresponding to the smallest eigenvalue of $L_{1}$, but we require each coordinate of $D^{-1/2}u$ to have unit modulus. Provided $u$ has no entries which are zero, we can normalize the entries to form an estimate of $\omega$, and as we show in the appendix (using results from ), this estimate is stable provided the spectral gap of $G^{\prime}$ is sufficiently large. This spectral method is known in the literature as angular synchronization , and we summarize the procedure in Algorithm 4

To reiterate, Algorithm 4 will produce estimates for the phases of the inner products $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime}}$. Also, we can take square roots of the vertex measurements $\{|\langle x,\varphi_{i}\rangle|^{2}+\nu_{i}\}_{i\in V^{\prime}}$ to estimate $\{|\langle x,\varphi_{i}\rangle|\}_{i\in V^{\prime}}$. Then we can combine these to estimate $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime}}$. However, note that the largest of these inner products will be most susceptible to noise in the corresponding phase estimate. As such, we remove a small fraction of these largest vertices so that the final collection of vertices $V^{\prime\prime}$ has size $\kappa|V|$, where $V$ was the original vertex set, and $\kappa$ is sufficiently close to $1$.

Now that we have estimated the phases of $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime\prime}}$, we wish to reconstruct $x$ by applying the Moore-Penrose pseudoinverse of $\{\varphi_{i}\}_{i\in V^{\prime\prime}}$. However, since $V^{\prime\prime}$ is likely a strict subset of $V$, it can be difficult in general to predict how stable the pseudoinverse will be. Fortunately, a recent theory of numerically erasure-robust frames (NERFs) makes this prediction possible: If the members of $\Phi_{V}$ are independent Gaussian vectors, then with high probability, every submatrix of columns $\Phi_{V^{\prime\prime}}$ with $\kappa=|V^{\prime\prime}|/|V|$ sufficiently large has a stable pseudoinverse . This concludes the phase retrieval procedure, briefly outlined below together with the measurement design.

Fix $d>2$ even and $\varepsilon\in(0,d-2\sqrt{d-1})$.

Given $M$, pick some $d$-regular graph $G=(V,E)$ with spectral gap $\lambda_{2}\geq\lambda^{\prime}:=1-\frac{2\sqrt{d-1}+\varepsilon}{d}$ and $|V|=cM\log M$ for $c$ sufficiently large, and arbitrarily direct the edges.

Prune the remaining induced subgraph for connectivity, producing the vertex set $V^{\prime}$ (Algorithm 3).

Estimate the phases of the vertex measurements using angular synchronization (Algorithm 4).

Remove the vertices with the largest measurements, keeping only $|V^{\prime\prime}|=\kappa|V|$.

Having estimates for $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime\prime}}$ up to a global phase factor, find the least-squares estimate of $[x]$ by applying the Moore-Penrose pseudoinverse of $\{\varphi_{i}\}_{i\in V^{\prime\prime}}$, see (3).

Having established our measurement design and phase retrieval procedure for the noisy case, we now present the following guarantee of stable performance:

Numerical results

In the previous sections, we described measurement designs and phase retrieval procedures for both the noiseless and noisy cases. This section presents results from numerical simulations to illustrate how well our phase retrieval procedures perform in practice. In particular, we will consider the noiseless and noisy cases separately.

For this case, we consider a slightly different measurement design. Rather than drawing a $d$-regular graph of $n$ vertices at random, we instead draw an Erdős-Rényi random graph. That is, for a fixed $n$ and $c\leq n$, we take $n$ vertices, and place an edge between each pair of vertices independently with probability $p:=c/n$. Note that for this model, the mean degree of each vertex is $p(n-1)=c(n-1)/n\approx c$. As we will see, this slight change to the graph model will not adversely affect the quality of our methods.

For each $(r,d)$ pair, we performed 30 trials of this process, and we populated the corresponding cell in Figure 1 with a shade of gray according to the proportion of good estimates (black indicates that none of the estimates were good, while white indicates that all of the estimates were good).

Having established this phase transition, we can use it to minimize the number of measurements. In particular, the total number of edges in the graph tends to be around $\frac{nc}{2}$, and so the total number of measurements is $N\approx(1+\frac{3}{2}c)n$. As before, we seek to minimize redundancy:

The right-hand side above is minimized when $r\approx 1.28$, in which case we get a redundancy of $\frac{N}{M}\approx 5.02\ll 236$. The reason for this disparity is simple: In the previous expander-graph-based analysis, we were chiefly concerned with ensuring that our measurement vectors $\Phi$ lend injective intensity measurements, so that we could reconstruct any given signal. On the other hand, the above analysis demonstrates that we can get away with far fewer measurement vectors if we only need to be able to reconstruct almost every signal. In this sense, these numerical simulations fail to capture the most challenging feature of measurement design for phase retrieval: injectivity (versus unique representation of almost every signal).

2 The noisy case

Interestingly, Phase Retrieval Procedure B produces relative errors similar to those gotten by doing least-squares estimation from $\Phi^{*}_{V}x+\nu_{V}$. This suggests that the portion of Phase Retrieval Procedure B which reconstructs (most of) the phases in $\Phi^{*}_{V}x+\nu_{V}$ performs rather well. In fact, this shows that the polarization trick of using edges to estimate vertex phases is particularly successful. Notice that performing least-squares estimation from $\Phi^{*}x+\nu$ produces a much smaller relative error, as expected. As far as runtime is concerned, Phase Retrieval Procedure B is slower than the phase oracles because it takes some time to estimate vertex phases with angular synchronization; however, this is not a substantial difference in runtime, as our procedure still produces an estimate in less than one second.

We would like to point out that in pruning for connectivity, instead of directly applying Algorithm 3, it sufficed to find the largest surviving component, since in our trials, this component always had spectral gap larger than $\tau=0.1$. We suspect that this is an artifact of the random graph, as this will certainly not happen in general.

The most striking thing about this simulation is that alternating projections consistently produces a slightly better estimate than Phase Retrieval Procedure B. For comparison, we also ran alternating projections using only the $3M$ vertex measurement vectors $\Phi_{V}$, and the relative errors were consistently on the order of $1$, i.e., alternating projections consistently stalled in this case. This suggests that there is some fundamental quality about the polarized measurement vectors $\Phi$ which makes them particularly well-suited for alternating projections, and we intend to study this in the future. Regardless, alternating projections took a lot longer to terminate; to be clear, we terminated the loop once applying both projections moved the estimate by less than $10^{-3}$, or by the $100$th iteration, whichever occurred first.

Concluding remarks

This paper provides a new way to perform phase retrieval, and our main result (Theorem 5) shows that our method is stable. In comparing with the stability result of , we note that neither result is completely satisfying when viewed from the perspective of application: In the real world, you are given a noise level and an acceptable level of estimate error, and you are asked to meet these specifications with signal processing techniques. For phase retrieval, the available guarantees fail to prescribe a measurement design that overcomes a given noise level—rather, they merely establish that with sufficiently many measurements, there exists some level of stability, i.e., $K$ in Theorem 5 or $C_{0}$ in Theorem 1.2 of . This reveals a gap in what is known about stability in phase retrieval, and we leave this for future work.

Admittedly, there are several gaps remaining between modern theory and application of phase retrieval. For example, thoughout this paper, it is assumed that the user has complete knowledge of the measurement design, but this is not always possible in practice. This can be resolved in part with new stability results which account for “noise” in the measurement design (this is sometimes called mismatch error).

One might feel that our phase retrieval algorithms are slightly unsatisfying because we perform hard thresholds to remove vertices according to how small or large the corresponding measurements are. Alternatively, there could very well be a way to more smoothly weight these measurements according to our confidence in them, and such weightings are already accounted for in the theory of angular synchronization . However, we decided to use hard thresholds because they greatly simplify the analysis of projective uniformity (though the analysis is still rather technical).

While the worst-case analysis we provide here is useful in many applications (and enables a comparison with the worst-case stability results of ), stochastic noise is a more appropriate model in other applications. We believe that the phase retrieval procedure of this paper will perform substantially better in the average case, but we leave this analysis for future work. Also, a notable distinction between our measurement designs in the noiseless and noisy cases is the presence of a log factor in the number of measurements used. However, we believe this factor is an artifact of our current analysis, and we intend to remove it in the future.

Appendix

This section proves the following guarantee:

Take proportions $p\geq q\geq\frac{2}{3}$, and consider a regular graph $G=(V,E)$ with spectral gap $\lambda_{2}>g(p,q):=1-2(q(1-q)-(1-p))$. After Algorithm 1 removes at most $(1-p)|V|$ vertices from $G$, then setting $\tau=\frac{1}{8}(\lambda_{2}-g(p,q))^{2}$, Algorithm 3 outputs a subgraph with at least $q|V|$ vertices.

To prove this theorem, we will apply a graph version of the Cheeger inequality, which provides a guarantee for Algorithm 2:

Consider a graph $G=(V,E)$ with spectral gap $\lambda_{2}$. Then Algorithm 2 outputs a set of vertices $S$ such that $h(S)\leq\sqrt{2\lambda_{2}}$.

First, Algorithm 1 removes a set of vertices, which we denote by $S_{0}$. In applying Algorithm 3, the $i$th step of the while loop removes another set of vertices $S_{i}$. We claim this while loop will end with $|\bigcup_{i\geq 0}S_{i}|<(1-q)|V|$. Supposing to the contrary, consider the first $k$ for which $S:=\bigcup_{i=0}^{k}S_{i}$ has at least $(1-q)|V|$ vertices. Then since each iteration of the while loop removes at most half of the remaining vertices, we have $(1-q)|V|\leq|S|\leq(1-\frac{q}{2})|V|$.

Next, Theorem 7 bounds each $h(S_{i})$ in terms of the spectral gap the remaining graph, which is necessarily less than $\tau$ by the condition of the while loop. Thus, we continue:

For the lower bound, we use the expander mixing lemma, which says

which, as substitution reveals, contradicts our choice for $\tau$. ∎

2 Angular synchronization

This section proves the following guarantee:

where $P:=\min_{\{i,j\}\in E}|\overline{\langle x,\varphi_{i}\rangle}\langle x,\varphi_{j}\rangle+\varepsilon_{ij}|$ and $C$ is a universal constant.

Note that one way to reconstruct $\omega^{*}$ is to minimize this quantity. Indeed, $\eta(\omega)$ is small when the angular differences $\omega_{j}-\omega_{i}$ are close to the measured differences $\rho_{ij}$, and in the noiseless case, $\eta(\omega)=0$ precisely when $\omega=\omega^{*}$, provided the graph $G$ is connected. In terms of this objective function, the following guarantee ensures that the output of Algorithm 4 is no worse than a constant multiple of optimal:

where $C^{\prime}$ is a universal constant.

By abuse of notation, we identify $a$ and $b$ with their coset representatives in $[-\pi,\pi)$. Since

then rearranging gives the following inequality:

where the last equality requires $a-b\in[-\pi,\pi]$. Otherwise, we note that

The reverse triangle inequality gives $|\tfrac{b}{|b|}-b|=|1-|b||=||a|-|b||\leq|a-b|$, and so the triangle inequality gives $|a-\tfrac{b}{|b|}|\leq|a-b|+|b-\tfrac{b}{|b|}|\leq 2|a-b|$. ∎

This relationship will allow us to apply Theorem 9. To this end, for notational convenience, we define $\gamma_{i}:=\arg(u_{i})-\arg(\langle x,\varphi_{i}\rangle)$ for every $i\in V$. The right-hand inequality of (13) and Lemma 10 together give

Denoting $\hat{\omega}_{i}:=\arg(u_{i})$ and $\omega^{*}_{i}:=\arg(\langle x,\varphi_{i}\rangle)$, then the definition of $\eta$ gives

With this, we continue (14) by applying the left-hand inequality of (13):

where the second inequality follows from Theorem 9. Furthermore, the right-hand inequality of (13) gives

i.e., $D^{1/2}w$ is orthogonal to $D^{1/2}1$. Also, since $(D-A)1=0$, we have

and so $D^{1/2}w$ is orthogonal to a first eigenvector of $L$, considering $L$ is positive semidefinite. Thus,

Continuing, we apply the definitions of $D$ and $A$ to get

From here, we proceed in two cases. First, when $\alpha\neq 0$, we may take $\theta:=\arg(\alpha)$. Then the left-hand inequality of (13) and Lemma 11 give

where the last inequality uses the fact that $\deg(i)\geq 1$ for every $i\in V$, i.e., the graph $G$ has no isolated vertex since $G$ is connected, which follows from the fact that $\tau>0$. In the case where $\alpha=0$, we may arbitrarily take $\theta=0$. Then similar analysis yields

The last inequality follows from Lemma 4, taking $z=\overline{\langle x,\varphi_{i}\rangle}\langle x,\varphi_{j}\rangle+\varepsilon_{ij}$ and $\varepsilon=-\varepsilon_{ij}$. ∎

3 Projective uniformity

This section is motivated by Theorem 8 of the previous section, which exhibits significant dependence on the size of $P$. Here, we show how Algorithm 1 ensures that $P$ will not too small, and our guarantee will be in terms of the following noise-robust version of projective uniformity:

for every noise vector $\varepsilon=\{\varepsilon_{ij}\}_{\{i,j\}\in E}$.

To prove this theorem, we apply the following lemma, which follows from concentration-of-measure arguments that we provide later:

then by Lemma 14, we have $\|\varepsilon\|_{2}\geq\frac{C^{\prime}}{\sqrt{M}}$ with overwhelming probability for some constant $C^{\prime}>0$. It remains to consider the case where (19) does not hold. To this end, define

Recalling the definition of $\mathcal{J}(\alpha,x,\varepsilon)$ in Definition 12, we claim that $\mathcal{J}(\alpha,x,\varepsilon)\subseteq\mathcal{J}_{1}(\alpha,x,\varepsilon)$. To see this, label each edge $\{i,j\}\in E$ with $|\overline{\langle x,\varphi_{i}\rangle}\langle x,\varphi_{j}\rangle+\varepsilon_{ij}|$. Then by definition, $V\setminus\mathcal{J}_{1}(\alpha,x,\varepsilon)$ neighbors more of the smallest edges than any other collection of $n-\lceil\alpha n\rceil=\lfloor(1-\alpha)n\rfloor$ vertices. By comparison, Algorithm 1 effectively deletes the smallest edge $\{i,j\}$ by deleting both of its incident vertices $i$ and $j$, one of which must be in $\mathcal{J}_{1}(\alpha,x,\varepsilon)$. The smallest edge in the remaining graph is guaranteed to not touch either $i$ or $j$, but rather some $k\in\mathcal{J}_{1}(\alpha,x,\varepsilon)$, provided $\mathcal{J}_{1}(\alpha,x,\varepsilon)$ has not yet been completely removed from the graph. After $\lfloor(1-\alpha)n\rfloor$ iterations, then by the pigeonhole principle, all vertices in $V\setminus\mathcal{J}_{1}(\alpha,x,\varepsilon)$ will be removed, meaning $\mathcal{J}(\alpha,x,\varepsilon)\subseteq\mathcal{J}_{1}(\alpha,x,\varepsilon)$, as claimed. This implies that

and so minimizing over all unit vectors $x$ gives

Continuing, we use the fact that $\mathcal{J}_{3}(\alpha,x,\varepsilon)\subseteq\mathcal{J}_{2}(\alpha,x,\varepsilon)$ to get

We conclude this section with the rather technical proof of Lemma 14:

thereby implying $|\mathcal{N}_{\delta}|\leq(\frac{2}{\delta}+1)^{2M}$.

As such, we first consider the success probability of these Bernoulli random variables:

where the last step is by the union bound. Since $a_{1}$’s density function is $\leq\sqrt{\frac{M}{\pi}}$, it follows that

For the other term in (24), note that $2M\|\varphi_{i}\|^{2}$ is a sum of independent standard Gaussian random variables. Applying Lemma 1 of then gives that for every $t>0$,

Substituting (25) and (25) into (24) then gives

Now, to bound (23), we will apply Hoeffding’s inequality , which says that the tail probability of a sum of independent Bernoulli random variables $X_{i}$, each with success probability $p$, has the following bound:

Also, note that replacing $p$ in the left-hand side above with some $p^{\prime}\geq p$ will not increase the probability. As such, taking $p^{\prime}$ to be the right-hand side of (27) and $t=1-\alpha-p^{\prime}$, we have

for some constants $c_{0},c_{1}>0$. Considering the above analysis and taking logarithms, it suffices to have

with $\delta^{2}=\frac{C^{\prime}}{M}$. Since $\alpha<1-\frac{1}{2C}$ by assumption, picking $C^{\prime}<\frac{\pi}{72}(1-\frac{1}{2C}-\alpha)^{2}$ will make $2n(1-\alpha)$ the dominant term in the above inequality, thereby proving the result. ∎

4 Removing large vertices

In this section, we prove how well we can remove the vertices with the largest noisy intensity measurements. We start with a lemma:

Taking $\delta=\frac{1}{2}\sqrt{\frac{C}{M}}$, we then have that $|\langle x,\varphi_{i}\rangle|^{2}>\frac{C}{M}$ implies $|\langle v_{x},\varphi_{i}\rangle|^{2}>\frac{C}{4M}$. Recall that we wish to bound the probability of the event $\mathcal{E}$ that (28) is violated for some unit vector $x$. To this end, the above implication allows us to focus on a finite set of points:

As discussed in the proof of Lemma 14, we may take $|\mathcal{N}_{\delta}|\leq(\frac{2}{\delta}+1)^{2M}$, and so the union bound and the symmetric distribution of $\Phi$ both give

where the last step is by the union bound. We continue, using the fact that $a_{1}$ and $b_{1}$ are both distributed as $\mathcal{N}(0,\frac{1}{2M})$:

where $z$ is a standard Gaussian random variable. With this, we now apply Hoeffding’s inequality to (29):

In counting large $z_{i}$’s, we identify which come from large or small inner products $|\langle x,\varphi_{i}\rangle|^{2}$. First,

is of size $<\beta n$ with overwhelming probability by Lemma 15. The rest of the large $z_{i}$’s have indices in

To count these, note that $\nu_{i}=z_{i}-|\langle x,\varphi_{i}\rangle|^{2}\geq\frac{C}{M}\|x\|^{2}$, and so

Rearranging then reveals that $|\mathcal{K}|=\mathcal{O}(M)$, meaning $\mathcal{K}$ has fewer than $\beta n\geq 3M\log M$ members when $M$ is sufficiently large. ∎

5 Main result

This section proves the main result of the paper, which we restate here:

We will prove the result by considering the steps of our phase retrieval process in reverse order. In the last step, we have the following estimates of $\langle x,\varphi_{i}\rangle$ for every vertex $i\in V^{\prime\prime}\subseteq V$ which survives our graph-pruning and large-vertex-removing processes:

Here, $\theta$ is a global phase which is calculated in the proof of Theorem 8. From these estimates, we reconstruct by finding the least-squares estimate of $x$:

As such, we have the following bound on the reconstruction error:

Next, we wish to bound $\|\delta\|$. By definition, we have

Note that the above square root operates under the assumption that $z_{i}\geq 0$ for each $i\in V^{\prime\prime}$, which is ensured when we prune for reliability. Denote $\xi_{i}:=\sqrt{z_{i}}-|\langle x,\varphi_{i}\rangle|$. Then by the triangle inequality, we have

For any $a,b\geq 0$, then since $0\leq(a-b)^{2}=a^{2}-2ab+b^{2}$, we have $(a+b)^{2}=a^{2}+2ab+b^{2}\leq 2(a^{2}+b^{2})$. Applying this inequality to the right-hand side above then gives

Similarly, for any $a,b\geq 0$, then since $ab\geq\min\{a^{2},b^{2}\}$, we have $(a-b)^{2}=a^{2}-2ab+b^{2}\leq|a^{2}-b^{2}|$. Applying this to $\xi_{i}^{2}=(\sqrt{z_{i}}-|\langle x,\varphi_{i}\rangle|)^{2}$ then gives $\xi_{i}^{2}\leq|\nu_{i}|$. Also by Theorem 16, our large-vertex-removing process ensures that $z_{i}\leq\frac{2c_{1}}{M}\|x\|^{2}$ for every $i\in V^{\prime\prime}$. Combined, these facts imply

Applying Theorem 8, Definition 12 and Theorem 13 further gives

Recalling the definition of $\varepsilon$, we have

and so, letting $E^{\prime}$ denote the edges which remained after pruning for connectivity, the Cauchy-Schwarz inequality gives

Finally, we combine (32), (33), (34) and (35), along with $\|\nu_{V}\|_{1}\leq\sqrt{n}\|\nu_{V}\|\leq\sqrt{2CM\log M}\|\nu\|$:

Acknowledgments

The authors thank the anonymous referees for providing thoughtful suggestions that led to a more complete discussion of the context of our results. The authors also thank Prof. Amit Singer for insightful discussions. B. Alexeev was supported by the NSF Graduate Research Fellowship under Grant No. DGE-0646086, A.S. Bandeira was supported by NSF Grant No. DMS-0914892, M. Fickus was supported by NSF Grant No. DMS-1042701 and AFOSR Grant Nos. F1ATA01103J001 and F1ATA00183G003, and D.G. Mixon was supported by the A.B. Krongard Fellowship. The views expressed in this article are those of the authors and do not reflect the official policy or position of the United States Air Force, Department of Defense, or the U.S. Government.

References