Phase Retrieval Meets Statistical Learning Theory: A Flexible Convex Relaxation

Introduction

for some absolute constant $\delta\in\left(0,1\right)$. The above geometric intuition is explained in more detail in Section 3.1. While our main result simply assumes that the anchor vector is given by an oracle, which may use the existing measurements, we discuss in Section 1.1 some realistic scenarios where a valid anchor vector exists or can be computed.

With these assumptions in place, we propose the solution to the convex programIn the real case, (3) reduces to a linear program.

as a computationally efficient estimator for $\bm{x}_{\star}$. Of course, the points equal to $\bm{x}$ up to a global phase, namely,

In Lemma 2, below in Section 3, we establish a geometric condition that is sufficient to guarantee accurate estimation of $\bm{x}_{\star}$ via the convex program (3).

The sufficient condition given by Lemma 2 can be interpreted in terms of (non-)existence of a particularly constrained halfspace that includes all of the points $\bm{a}_{i}\bm{a}_{i}^{*}\bm{x}_{\star}$. For random measurement vectors $\bm{a}_{i}$, this interpretation resembles the model and theory of linear classifiers studied in statistical learning theory, albeit in an unusual regime. Borrowing classic results from this area (summarized in Appendix A), we show that with high probability (3) produces an accurate estimate of $\bm{x}_{\star}$.

Specifically, in our main result, Theorem 1 in Section 3, we show that drawing

i.i.d. random measurements, with the hidden constant factor on the right-hand side depending on $\delta$, would suffice for the conditions of Lemma 2 to hold with probability $\geq 1-\varepsilon$. Consequently, solution $\widehat{\bm{x}}$ of (3) would obey

Our approach critically depends on the choice of the anchor vector $\bm{a}_{0}$ that obeys an inequality of the form (2). Below we discuss two interesting scenarios where such a vector would be accessible.

Perhaps the simplest scenario is when the target signal $\bm{x}_{\star}$ is known to be real and non-negative. In usual imaging modalities these model assumptions are realistic as natural images are typically represented by pixel intensities. For these types of signals we can choose $\bm{a}_{0}=\frac{1}{\sqrt{N}}\bm{1}$ for which we obtain $\left|\bm{a}_{0}^{*}\bm{x}_{\star}\right|=\left\lVert\bm{x}_{\star}\right\rVert_{1}/\sqrt{N}$. Then, for (2) to hold it suffices that $\left\lVert\bm{x}_{\star}\right\rVert_{1}\geq\delta\sqrt{N}\left\lVert\bm{x}_{\star}\right\rVert_{2}$ for some absolute constant $\delta\in(0,1)$. In particular, we need $\bm{x}_{\star}$ to have at least $\delta^{2}N$ non-zero entries.

Random measurements:

A more interesting scenario is when we can construct the vector $\bm{a}_{0}$ from the (random) measurements. An effective strategy is to set $\bm{a}_{0}$ to be the principal eigenvector of the matrix $\bm{\varSigma}=\frac{1}{M}\sum_{i=1}^{M}b_{i}\bm{a}_{i}\bm{a}_{i}^{*}$. The principal eigenvector of $\bm{\varSigma}$ and its “truncated” variants have been used previously for initialization of the Wirtinger Flow algorithm (Candès et al., 2015b) and its refined versions (Chen and Candès, 2015; Zhang et al., 2016). For example, the following result is shown in Candès et al. (2015b, Section VII.H).

then (2) holds with probability $\geq 1-O(N^{-2})$.

While Lemma 1 can be refined or extended in various ways, we do not pursue these paths in this paper.

2 Related work

There is a large body of research on phase retrieval addressing various aspect of the problem (see (Jaganathan et al., 2015) and references therein). However, we focus only on the relevant results mostly developed in recent years. Perhaps, among the most important developments are PhaseLift and similar methods that cast the phase retrieval problem as a particular semidefinite program (Candès et al., 2013; Candès and Li, 2014; Waldspurger et al., 2015). The main idea used by Candès et al. (2013) and Candès and Li (2014) is that by lifting the unknown signal using the transformation $\bm{x}\bm{x}^{*}\mapsto\bm{X}$, the (noisy) phaseless measurements (1) that are quadratic in $\bm{x}_{\star}$ can be converted to linear measurements of the rank-one positive semidefinite matrix $\bm{X}_{\star}=\bm{x}_{\star}\bm{x}_{\star}^{*}$. With this observation, these SDP-based methods aim to solve the corresponding linear equations using the trace-norm to induce the rank-one structure in the solution. Inspired by the well-known convex relaxation of Max-Cut problem, PhaseCut method (Waldspurger et al., 2015) considers the measurement phases as the unknown variables and applies a similar lifting transform to formulate a different semidefinite relaxation for phase retrieval. While these SDP-based methods are shown to produce accurate estimates of $\bm{X}_{\star}$ at optimal sample complexity for certain random measurement models, they become computationally prohibitive in medium- to large-scale problems where SDP is practically inefficient.

More recently, there has been a growing interest in non-convex iterative methods for phase retrieval (see e.g., Netrapalli et al., 2013; Candès et al., 2015b; Schniter and Rangan, 2015; Chen and Candès, 2015; Zhang et al., 2016; Wang and Giannakis, 2016; Sun et al., 2016). These methods typically operate in the natural space of the signal and thus do not suffer the drawbacks of the SDP-based methods. With a specific initialization Netrapalli et al. (2013) establish some accuracy guarantees for a variant of the classic methods by Gerchberg and Saxton (1972); Fienup (1982) that iteratively update the estimate assuming the measurements’ phase match that of the previous iterate. The established sample complexity is (nearly) optimal in the dimension of the target signal, but it does not vary gracefully with the prescribed precision. Phase retrieval via the Wirtinger Flow (WF), a non-convex gradient descent method at core, is proposed by Candès et al. (2015b). It is shown that for random measurements that have Normal distribution or certain coded diffraction patterns, with an appropriate initialization the WF iterates exhibit the linear rate of convergence to the target signal. More recent work on the WF method introduce better initialization by excluding the outlier measurements and achieve the optimal sample complexity (Chen and Candès, 2015; Zhang et al., 2016). The WF class of algorithms and our proposed method both achieve optimal sample complexity (up to the constant factor) and have low computational cost. However, the WF methods need careful tuning of a step size parameter and their convergence analysis often relies on Gaussian measurements. This is partly because establishing robustness of non-convex methods generally requires stronger conditions. Our method provably works for a broader set of measurement distributions, has no tuning parameters, and can be implemented in various convex optimization software.

Shortly after a draft of this manuscript was first posted online, a few independent papers proposed and analyzed the same method and its variants. Goldstein and Studer (2016), who dubbed (3) PhaseMax, obtained sharper constants in the sample complexity by assuming a stronger condition that the anchor is independent of the measurements in their analysis. Alternative proofs and variations that rely on matrix concentration inequalities appeared later in (Hand and Voroninski, 2016c, b, a). Another distinctive feature of our analysis compared to the mentioned results is that it is less sensitive to measurement distribution as it relies on VC–type bounds.

3 Variations and Extensions

In this section we discuss several different ways to extend the proposed method that we leave for future research. While the core geometric idea still applies, some modifications of our theoretical arguments would be necessary to analyze these extensions.

The gross noise model considered in this paper can be pessimistic in scenarios where we have random noise or deterministic noise with a different type of bound. In these scenarios, augmenting the estimator by a noise regularization term could result in accuracy bounds that gracefully vary with the considered noise.

Another interesting extension to the proposed method, is to adapt the current theory to the case of blockwise independent measurements as in coded diffraction imaging. Our numerical experiments in Section (2) suggest that the proposed method still performs well with these structured measurements. Nevertheless, to extend the analysis we may need to revise the current simple arguments based on Vapnik-Chervonenkis theory using more sophisticated tools from the theory of empirical processes.

Numerical experiments

Theoretical Analysis

To understand the geometry of (3) it is worthwhile to first consider the noiseless scenario. The feasible set is the intersection of the sets

corresponding to the pairs $\left(\bm{a}_{i},b_{i}\right)$ for $i=1,2,\dotsm,M$. The sets $\mathcal{S}_{i}$ are effectively symmetric “complex slabs”. Denote their intersection by

cannot hold simultaneously. The following lemma provides the desired sufficient condition.

and $\epsilon\geq 0$ be some constant. If every vector $\bm{h}\in\mathcal{R}_{\delta}$ with $\left\lVert\bm{h}\right\rVert_{2}>\epsilon$ violates at least one of the inequalities

then any solution $\widehat{\bm{x}}$ to (3) obeys

It suffices to show that $\bm{h}=\widehat{\bm{x}}-\bm{x}_{\star}$ obeys (5) and it belongs to $\mathcal{R}_{\delta}$. Given that

and $\xi_{i}\leq\eta^{-1}$, we have $\langle\bm{a}_{i}\bm{a}_{i}^{*}\bm{x}_{\star},\bm{h}\rangle\leq\frac{1}{2}\eta^{-1}$. Feasibility of $\bm{x}_{\star}$ also guarantees that $\langle\bm{a}_{0},\bm{h}\rangle\geq 0$. Therefore, we have shown that $\bm{h}$ satisfies (5).

where $\bm{h}_{\perp}=\left(\bm{I}-\bm{x}_{\star}\bm{x}_{\star}^{*}\right)\bm{h}=\bm{h}-\left(\bm{x}_{\star}^{*}\bm{h}\right)\bm{x}_{\star}$. The above inequality completes the proof as it is equivalent to $\bm{h}\in\mathcal{R}_{\delta}$. ∎

2 Guarantees for random measurements

In this section we will show that if the vectors $\bm{a}_{i}$ for $1\leq i\leq M$ are drawn from a random distribution and (2) holds for a sufficiently large constant $\delta$, then with high probability (3) produces an accurate estimate of $\bm{x}_{\star}$. Our strategy is to show that for a sufficiently large $M$ the sufficient condition provided in Lemma 2 holds with high probability.

For $\delta\in(0,1)$ let $\mathcal{C}_{\delta}$ be the convex cone given by

where $\bm{x}_{\star}^{*}\bm{y}$ is implicitly assumed to be a real number. The polar cone of a set $\mathcal{C}$ is defined as

It is easy to verify that the polar cone of $\mathcal{C}_{\delta}$ is

Since $\bm{a}_{0}\in\mathcal{C}_{\delta}$ by assumption, it follows that for every $\bm{h}\in\mathcal{C}_{\delta}^{{}^{\circ}}$ we have $\langle\bm{a}_{0},\bm{h}\rangle\leq 0$. Therefore, the inequality $\langle\bm{a}_{0},\bm{h}\rangle\geq 0$ can hold only for vectors $\boldsymbol{z}$ in the closure of the complement of $\mathcal{C}_{\delta}^{{}^{\circ}}$ which we denote by

A typical positioning of $\bm{a}_{0}$ and $\bm{a}_{i}\bm{a}_{i}^{*}\bm{x}_{\star}$ needed to guarantee unique recovery is illustrated in Figure 4.

for a constant $p_{\min}(\delta,t)\in\left(0,1\right)$ depending only on $\delta$ and $t$.Clearly, the best $p_{\min}\left(\delta,t\right)$ decreases as $t$ increases. For any $\varepsilon>0$, if we have

with the hidden constant factor inversely related to $p_{\min}(\delta,t)$, then with probability $\geq 1-\varepsilon$ the estimate $\widehat{\bm{x}}$ obtained through (3) obeys

which is an approximation of the true probability of the event denoted by

Now, because $p(\bm{h})\geq p_{\min}(\delta,t)$ for all $\bm{h}\in\mathcal{C}^{\prime}_{\delta}\cap\mathcal{R}_{\delta}$ with $\left\lVert\bm{h}\right\rVert_{2}>\left(t\eta\right)^{-1}$, the above inequality implies that

If $M=\frac{8}{p_{\min}^{2}(\delta,t)}\left(c\cdot 2N+2\log\frac{8}{\varepsilon}\right)$, then we have

where we used the inequality $\log u-\log 2=\log\frac{u}{2}\leq\frac{u}{2}-1$ in the second line. Setting $c=2\log\frac{8e}{p_{\min}^{2}\left(\delta,t\right)}$, it follows that

This immediately implies that for $M\overset{\delta,t}{\gtrsim}N+\log\frac{1}{\varepsilon}$ we have

We can consider the case of measurements with normal distribution as a concrete example. To apply the Theorem 1, it suffices to quantify the constant $p_{\min}(\delta,t)$ which can be achieved through Lemma 3 below.

Since $\bm{h}\in\mathcal{R}_{\delta}$ and $\left\lVert\bm{h}\right\rVert_{2}>\left(t\eta\right)^{-1}$ we have

We consider two cases depending on $\left\lVert\bm{h}_{\perp}\right\rVert_{2}=0$ or not. If $\left\lVert\bm{h}_{\perp}\right\rVert_{2}=0$, then $\left|\langle\bm{x}_{\star},\bm{h}\rangle\right|>\left(t\eta\right)^{-1}$. The fact that $\bm{h}\in\mathcal{C}^{\prime}_{\delta}$ as well, implies that $\langle\bm{x}_{\star},\bm{h}\rangle$ is non-negative and thereby $\langle\bm{x}_{\star},\bm{h}\rangle>\left(t\eta\right)^{-1}$. Consequently, we have

If $\left\lVert\bm{h}_{\perp}\right\rVert>0$, then we can invoke Lemma 5 in the Appendix with $\alpha=\frac{\langle\bm{x}_{\star},\bm{h}\rangle}{\left\lVert\bm{h}_{\perp}\right\rVert_{2}}$ to show that

Then by rewriting (7) as $\frac{\left(t\eta\right)^{-2}}{\left\lVert\bm{h}_{\perp}\right\rVert_{2}^{2}}\leq 1+\delta^{-2}+\alpha^{2}$ we have

The fact that $\bm{h}\in\mathcal{C}^{\prime}_{\delta}$, guarantees that $\alpha\geq-\sqrt{\delta^{-2}-1}$. Since $\frac{\alpha}{\sqrt{\alpha^{2}+1}}$ and $-\frac{\sqrt{1+\delta^{-2}+\alpha^{2}}}{\sqrt{\alpha^{2}+1}+\alpha}$ are both increasing in $\alpha$, we obtain

The above lower bound is the smaller one of the two considered cases and thus the proof is complete. ∎

Appendix A Tools from statistical learning theory

For reference, here we provide some of the classic results in statistical learning theory that we employed in our analysis. We mostly follow the exposition of the subject presented by Devroye et al. (2013, chapters 13 and 14).

The $n$-th shatter coefficient (or growth function) of a class $\mathcal{F}$ of binary functions $f:\mathcal{X}\to\left\{0,1\right\}$ is defined as

Intuitively, the shatter coefficient $s(\mathcal{F},n)$ is the largest number of binary patterns that the functions in $\mathcal{F}$ can induce on $n$ points.

The Vapnik–Chervonenkis (VC) dimension of a class $\mathcal{F}$ of binary functions is the largest number $n$ such that $s(\mathcal{F},n)=2^{n}$, namely,

If $\mathcal{F}$ can induce all binary patterns on $n$ points, $\mathcal{F}$ is said to “shatter” $n$ points. Therefore, the VC–dimension of $\mathcal{F}$ is the largest number of points that $\mathcal{F}$ can shatter.

The following theorem is originally due to Vapnik and Chervonenkis (1971). We restate the theorem as presented in Devroye et al. (2013).

Let $\mathcal{F}$ be a class of binary functions and $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{n}$ be i.i.d. copies of an arbitrary random variable $\bm{x}$. Then for every $t>0$ we have

Appendix B Auxiliary Lemma

where the second line follows from the change of variable $v=\frac{\beta}{\sqrt{\alpha^{2}+1}}u^{-1}$. It is straightforward to show that $G(0)=\frac{1}{2\sqrt{\alpha^{2}+1}}$. A simple integration then yields

for $\beta\geq 0$, and since $G(\beta)$ is even for all $\beta$ we have

where the last line follows from the fact that $\frac{g}{\sqrt{v^{2}+g^{2}}}$ has a uniform distribution over $$. Replacing$F(0)$in the expression of$F(\beta)$ and straightforward simplifications yield the desired result. ∎

References