Phase retrieval from power spectra of masked signals

Afonso S. Bandeira, Yutong Chen, Dustin G. Mixon

Introduction

In many applications, one wishes to reconstruct a signal from the magnitudes of its Fourier coefficients. This problem is known as phase retrieval, and it has been instrumental to many important scientific advances, including the Nobel Prize–winning work that leveraged X-ray diffraction to establish the double helix structure of DNA . Of course, given only the magnitudes of a signal’s Fourier coefficients, one does not have enough information to recover the signal—while the Fourier transform is injective, the point-wise absolute value is not. As such, one is inclined to use a priori knowledge of the signal, and hope it is then uniquely determined by the Fourier magnitudes. For example, to deduce the structure of DNA, Watson and Crick applied certain chemical assumptions, along with a knowledge of van der Waals interactions between atoms.

In order to image more exotic molecules, such assumptions are difficult to apply, and so there has been quite a bit of work attempting to exploit more general assumptions (e.g., positivity or support constraints). To account for this sort of prior information, the most popular phase retrieval algorithms are modifications of Gerchberg and Saxton’s original approach , which alternates between the time and frequency domains, iteratively correcting the current guess by imposing time-domain assumptions or scaling Fourier coefficients to match the measured data. Marchesini surveys and compares the various modifications, but they all have a tendency to stall in local minima.

To address this issue, Candès, Eldar, Strohmer and Voroninski proposed an alternative methodology whereby nonuniqueness is overcome not by prior information, but by additional illuminations. For each illumination, a different mask (or grating) is used to distort the appearance of the object in question; in mathematical parlance, each mask acts as a multiplication operator on the desired signal before the Fourier transform. Furthermore, Candès et al. constructed three masks such that the corresponding illuminations uniquely determine almost every signal up to a global phase factor; a fourth illumination is necessary to uniquely determine all signals . However, the phase retrieval algorithms they propose are not guaranteed to reconstruct the desired signal from these three illuminations, even if the signal is uniquely determined by the illuminations.

How to polarize Fourier masks

and the existence of the requisite connected component is given by the following classical result in spectral graph theory (e.g., see Lemma 5.2 of ):

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

In summary, in order to apply the polarization trick in the setting of masked DFTs, it suffices to (i) construct a full spark ensemble $\Phi_{V}$ with masked DFTs, and (ii) construct a graph $G$ with a sufficiently large spectral gap and whose corresponding edge vectors $\Phi_{E}$ can also be implemented with masked DFTs. We quickly address (i) with the following result:

The $M\times KM$ matrix whose columns are $D_{k}f_{m}$ has Vandermonde form, as does each of its $M\times M$ submatrices. These submatrices are all invertible precisely when their determinants are nonzero, that is, when $\{\alpha_{k}e^{2\pi im/M}\}_{k=0,}^{K-1}\ \!{}_{m=0}^{M-1}$ are distinct (by the Vandermonde determinant formula). The lemma immediately follows. ∎

The condition above is satisfied with probability $1$ if the $\alpha_{k}$ ’s are drawn uniformly at random from the complex unit circle. For a deterministic alternative, it suffices to take $\alpha_{k}=e^{2\pi ik/KM}$ . Now that we have established how to construct masks $D_{k}$ such that $\Phi_{V}=\{D_{k}f_{m}\}_{k=0,}^{K-1}\ \!{}_{m=0}^{M-1}$ is full spark, we turn to solving (ii). Before describing our construction of $G$ , we first motivate it by illustrating how polarized combinations of our vertex vectors can be expressed using masked DFTs:

Take $\omega=e^{2\pi i/3}$ , let $E=\operatorname{diag}\{e^{2\pi im/M}\}_{m=0}^{M-1}$ denote the modulation operator, and consider $\Phi_{V}=\{D_{k}f_{m}\}_{k=0,}^{K-1}\ \!{}_{m=0}^{M-1}$ , as defined in Lemma 2. Then

Then $D_{k}f_{m}+\omega^{r}D_{k^{\prime}}f_{m^{\prime}}=D_{k}f_{m}+\omega^{r}E^{m^{\prime}-m}D_{k^{\prime}}f_{m}=(D_{k}+\omega^{r}E^{m^{\prime}-m}D_{k^{\prime}})f_{m}$ . ∎

When implementing the modulation trick of Lemma 3 using DFTs, we will have to fix $m^{\prime}-m$ (so as to apply a fixed mask) and let $m$ vary (since we will mask an entire DFT). As such, we intend to use auxiliary masks of the form $\{D_{k}+\omega^{r}E^{a}D_{k^{\prime}}\}_{r=0}^{2}$ , thereby drawing an edge between every $(k,m)$ and $(k^{\prime},m^{\prime})$ such that $m^{\prime}-m=a\bmod M$ . This informs our decision of how to construct the graph $G$ , since it enables a masked-DFT implementation:

Fourier bias is used in additive combinatorics to measure pseudorandomness; here, the main idea is that correlation with any complex sinusoid indicates regularity, which is not typically exhibited by random sets. As indicated earlier, the spectral gap of $G$ is intimately related to the Fourier bias of $A$ :

Consider the graph $G$ defined in Definition 4. The spectral gap of $G$ is

where $\|\cdot\|_{u}$ denotes Fourier bias.

To determine the spectral gap, we first determine the adjacency matrix $W$ of $G$ . For any pair $k,k^{\prime}\in\{0,\ldots,K-1\}$ , the adjacency rule for $(k,m)$ and $(k^{\prime},m^{\prime})$ is whether $m^{\prime}-m\in A$ . As such, the $(k,k^{\prime})$ th block of $W$ is circulant, with each row being a translation of $\mathbf{1}_{A}$ . This gives the expression $W=J\otimes\operatorname{circ}(\mathbf{1}_{A})$ , where $J$ is the $K\times K$ all-ones matrix and $\otimes$ denotes the Kronecker product. One useful property of the Kronecker product is that its eigenvalues are products of eigenvalues:

Here, we note that the only nontrivial eigenvalue of $J$ is $K$ , and the eigenvalues of $\operatorname{circ}(\mathbf{1}_{A})$ are the entries of the Fourier transform $MF^{*}\mathbf{1}_{A}$ . As such, the nontrivial eigenvalues of $W$ are given by

with equality when $m=0$ . This means the largest eigenvalue of $W$ is $\lambda_{1}=K|A|$ , which corresponds to the all-ones eigenvector, and so the spectral gap of $G$ is

Small sets with small Fourier bias

Suppose the entries of $\mathbf{1}_{B}$ are independent, identical Bernoulli random variables with mean $\frac{c\log M}{M}$ . Then the following simultaneously hold with high probability:

$\frac{1}{2}c\log M\leq|B|\leq\frac{3}{2}c\log M$

Next, we have $\operatorname{Pr}(0\in B)=\frac{c\log M}{M}$ . For (iii), we apply a multiplicative form of the Chernoff bound (see (6) and (7) in ): Let $X$ be a sum of independent 0-1 random variables. Then

The result then follows from a union bound. ∎

Note that in the event of Lemma 6, we have $\mathbf{1}_{A}=\mathbf{1}_{B}+\mathbf{1}_{-B}$ , and so

where the last step applied a complex conjugate to $(F^{*}\mathbf{1}_{-B})(m)$ . As such, it suffices to show that random sets $B$ have small Fourier bias:

Pick $c\geq 4$ and suppose the entries of $\mathbf{1}_{B}$ are independent, identical Bernoulli random variables with mean $\frac{c\log M}{M}$ . Then $\|B\|_{u}\leq 3\sqrt{c}\cdot\frac{\log M}{M}$ with high probability.

To prove this lemma, we will apply the following version of the Chernoff bound:

Let $X$ and $Y$ denote the real and imaginary parts of $Z$ . Then a union bound gives

Denote $\sigma_{X}:=\sqrt{\operatorname{Var}(X)}$ and $\sigma_{Y}:=\sqrt{\operatorname{Var}(Y)}$ . Since $\operatorname{Var}(X)+\operatorname{Var}(Y)=\operatorname{Var}(Z)$ , we then have $\sigma_{X}\leq\sigma$ and $\sigma_{Y}\leq\sigma$ , and so

The result then follows from applying Theorem 1.8 of to both terms of the right-hand side. ∎

From here, we apply the binomial theorem to get

Here, another application of the binomial theorem gives

and so by the geometric sum formula, we have

Having calculated the expected value and variance of $Z$ , we are now ready to use Lemma 8 (and the fact that $Z=\frac{M}{2}(F^{*}\mathbf{1}_{B})(m)$ ). A union bound gives

We select $t=(\sqrt{2c(1+\varepsilon)}\log M)/\sigma$ and simplify the exponents in our probability bound:

Since $c\geq 4$ , we therefore have that $\|B\|_{u}\leq 3\sqrt{c}\cdot\frac{\log M}{M}$ with high probability. ∎

We can now combine Lemmas 6 and 7 to produce a graph $G$ of the form in Definition 4 with a large spectral gap:

Pick $c\geq 4$ and suppose the entries of $\mathbf{1}_{B}$ are independent, identical Bernoulli random variables with mean $\frac{c\log M}{M}$ . Take $A:=B\cup(-B)\setminus\{0\}$ and define $G$ according to Definition 4. Then

We will assume that the events of Lemmas 6 and 7 hold simultaneously, as they will with high probability. Then $|A|=2|B|$ , and furthermore by (2), $\|A\|_{u}\leq 2\|B\|_{u}$ . Starting with Lemma 5, we then have

where the second inequality applies Lemmas 6(iii) and 7. ∎

We now apply Lemma 1 to show that $\mathcal{O}(\log M)$ Fourier masks suffice for injectivity:

Take $K=12$ , $c=144$ , suppose the entries of $\mathbf{1}_{B}$ are independent, identical Bernoulli random variables with mean $\frac{c\log M}{M}$ , and take $A:=B\cup(-B)\setminus\{0\}$ . Then with high probability, the $\leq 2\cdot 10^{5}\cdot\log M$ masks

as described in Lemmas 2 and 3, lend injective intensity measurements.

Note that $G$ has $n=KM$ vertices and is $d$ -regular with $d=K|A|$ . We need to be robust to the removal of any $M-1$ vertices, which in turn removes $d(M-1)$ edges, and so it suffices to be robust to the removal of any $dM$ edges. As such, we observe Lemma 1 and take $\varepsilon$ to satisfy

We also want a connected component of size $M$ after the removal of these edges, and so by Lemma 1, it suffices to have

Rearranging (3) and (4) produces the following specification on the spectral gap of $G$ for sufficient connectivity:

where the equality is valid provided $K\geq 2$ . To meet this specification, based on Theorem 9, it suffices to have $1-6/\sqrt{c}\geq 6/K$ , which requires $K>6$ and is equivalent to having $c\geq(1/6-1/K)^{-2}$ . Using Lemma 6, the total number of masks is

Setting $c=(1/6-1/K)^{-2}$ , then $K=12$ minimizes the coefficient of $\log M$ . ∎

At this point, we note that $2\cdot 10^{5}\cdot\log M$ is a rather large number of Fourier masks. Certainly, the $10^{5}$ might be an artifact of our analysis—perhaps it could be decreased by leveraging better approximations. However, as the next result shows, the log factor is necessary for the Fourier masks to lend polarization-based recovery:

Take $G$ as defined in Definition 4. Then $\operatorname{spg}(G)>\varepsilon$ only if

Define $V$ to be the $|A|\times M$ matrix built from taking the rows of $F$ indexed by $A$ and then scaling the columns to have unit norm. Then the inner product between any two columns of $V$ is given by

As such, the worst-case coherence between columns of $V$ can be expressed in terms of the Fourier bias of $A$ (and the spectral gap of $G$ by Lemma 5):

where the last step is by (5). Thus $\delta^{2}\geq 2\varepsilon$ , with which we continue (6):

Here, the last inequality can be verified using the fact that $\varepsilon\leq 1$ . ∎

Numerical simulations

In this section, we compare the polarization method of this paper to a state-of-the-art phase retrieval algorithm known as PhaseLift . The main idea behind PhaseLift is that the intensity measurement $x\mapsto|\langle x,\varphi\rangle|^{2}$ can be viewed as a linear measurement if the signal $x$ is “lifted” to the outer product $xx^{*}$ in the real vector space of self-adjoint $M\times M$ matrices. Indeed, the intensity measurement is a Hilbert-Schmidt inner product in this vector space:

However, this larger vector space has $M^{2}$ dimensions, and so $M^{2}$ inner products are necessary to identify any member of the space—that is, unless more information is available. In this case, we know that the desired self-adjoint matrix $xx^{*}$ is positive semidefinite with rank $1$ , and so one could seek to minimize rank over the positive semidefinite matrices with the given intensity measurements. Rank minimization tends to be a difficult program to solve, so one is inclined to relax it:

First, we describe how we implement phase retrieval with polarization. Starting with an ensemble of vertex measurement vectors $\Phi_{V}$ and edge measurement vectors $\Phi_{E}$ , then for each edge $(i,j)\in E$ , we apply (1) to get

where $\varepsilon_{ijr}$ denotes the noise added to the $(i,j,r)$ th edge measurement. Later, we will normalize $w_{ij}$ in order to estimate the relative phase between the measurements at vertices $i$ and $j$ , but if $w_{ij}$ is small, this normalization will be particularly susceptible to noise. As such, we first remove vertices which are adjacent to small edge weights so as to promote reliability in the edges.

Now that we have isolated a subgraph with reliable edges, we want to find a further subgraph which is sufficiently connected so that its vertices can democratically agree on the relative phases. To find this subgraph, we iteratively remove vertices implicated by spectral clustering until the spectral gap is sufficiently large.

At this point, our graph has reliable edges and is well connected. We now run angular synchronization to reach a consensus on the phases of the vertex measurements, up to a global phase factor.

Having estimated the phases of the vertex measurements corresponding to $V^{\prime\prime}$ , we now multiply the square roots of the vertex measurements $\{\sqrt{|\langle x,\varphi_{i}\rangle|^{2}+\varepsilon_{i}}\}_{i\in V^{\prime\prime}}$ by these phases to estimate the inner products $\{\langle x,\varphi_{i}\rangle\}_{i\in V^{\prime\prime}}$ , and then produce a least-squares estimate for the desired signal $x$ . As established in , this implementation of the polarization method is stable to noise when the vertex measurement vectors are complex Gaussian. In this section, we run numerical simulations to illustrate stability in the Fourier masks setting.

By contrast, PhaseLift offers a lot more flexibility with the number and types of masks used, and this flexibility allows us to run a collection of choice comparisons. We will run experiments at three noise levels, and for each level, we will compare polarization to three different implementations of PhaseLift. In the first implementation, we will only give PhaseLift the measurements corresponding to the original $3$ masks we give to polarization. Next, we will give PhaseLift all $18|A|+3$ of the masks we described in the previous paragraph. For the last comparison, we will give PhaseLift the original $3$ masks along with $18|A|$ additional masks whose diagonal entries are also independent with distribution $\mathcal{N}(0,1)$ . Based on discussions in , this last setup appears to be the intended design of Fourier masks for PhaseLift, so in a sense, this last comparison will allow both algorithms to compete with their own “home-field advantage.”

In the comparison between polarization and PhaseLift, we use two performance metrics: run time and relative error of reconstruction, defined in this setting by

Here, $c$ is playing the role of the global phase factor we lost in the intensity measurement process. For each noise level tested $\sigma^{2}\in\{0,0.1,1\}$ , and for each type of mask ensemble given to PhaseLift (original $3$ , same $18|A|+3$ , random $18|A|+3$ ), we considered multiple signal dimensions $M\in\{2^{5},2^{6},2^{7},2^{8},2^{9}\}$ . In each of these scenarios, we ran $30$ realizations of the following experiment:

Draw each entry of the signal $x$ independently from $\mathcal{N}(0,1)$ .

Add independent $\mathcal{N}(0,\sigma^{2})$ noise to each intensity measurement $|\langle x,\varphi\rangle|^{2}$ .

Run both polarization and PhaseLift to estimate $x$ and record the run time and relative error.

We ran these experiments on an Intel® Core™2 Quad CPU Q9550 at 2.83GHz with 3.6 GB of memory.

Figure 1 presents the results of the noiseless scenarios ( $\sigma^{2}=0$ ). Here, perhaps the most striking distinction between polarization and PhaseLift is the relative error; indeed, polarization produces an estimate with almost no error ( $\sim 10^{-14}$ ), whereas PhaseLift produces noticeable error ( $\sim 10^{-2}$ ). In addition, polarization beats PhaseLift when using the same number of masks. For higher dimension signals, this advantage translates to minutes versus hours of running time. Letting $N$ denote the number of intensity measurements, then PhaseLift takes $\mathcal{O}(N^{3.5})$ operations, whereas the computational bottleneck in polarization is the least-squares step, which takes $\mathcal{O}(N^{3})$ operations, but this does not account for the disparity in run time. Rather, the disparity is indicative of a large constant in front of the $N^{3.5}$ . In comparing the relative error in reconstruction, we note that PhaseLift terminates when the step size is sufficiently small, so this accounts for the significant difference in performance. Overall, in the noiseless case, polarization handily defeats PhaseLift, even when PhaseLift is given the random masks it prefers.

Next, Figure 2 illustrates how polarization and PhaseLift perform in the present of some noise ( $\sigma^{2}=0.1$ ). In this case, PhaseLift continues to be rather slow, and it is clear that polarization enjoys greater stability from the $18|A|$ additional masks when compared to PhaseLift’s performance with just the original $3$ random masks. In , it is claimed that $3$ random masks are sufficient for PhaseLift to reconstruct typical signals (and this is corroborated with our experiments in the noiseless case), but Figure 2 illustrates that PhaseLift requires more masks in order to reconstruct stably. And indeed, when PhaseLift is given more masks, it performs much better. Specifically, when they use the same masks, polarization and PhaseLift produce estimates with statistically comparable relative error, and PhaseLift performs slightly better when given the random masks it prefers. Intuitively, it makes sense that PhaseLift should be at least slightly more stable than polarization since it leverages all of the measurements in a democratic fashion, whereas polarization disenfranchises the vertex and edge measurements that are too small (i.e., “wasting” measurements).

Figure 3 presents the results of our high-noise scenario ( $\sigma^{2}=1$ ), which are very similar to the results corresponding to $\sigma^{2}=0.1$ .

Discussion

In this paper, we showed how to leverage the polarization-based phase retrieval technique to construct $\Theta(\log M)$ Fourier masks that uniquely determine every $M$ -dimensional signal, and then we used numerical simulations to illustrate the stability of polarization with such masks. At this point, we offer two conclusions:

The polarization technique for phase retrieval is flexible enough to provably accommodate certain measurement design criteria (such as Fourier masks).

If you have the ability to use polarization instead of PhaseLift, you will gain considerable speedups in run time at the price of a slight increase in relative error.

It remains to be rigorously proved that stability holds in the Fourier masks setting. We did not exploit the inherent Fourier structure to further speed up reconstruction, though this could very well be possible. The number of masks might be decreased if additional information about the signal were leveraged, as in . Finally, we still seek Fourier masks–based performance guarantees for PhaseLift and its modifications .

Acknowledgments