Solving (most) of a set of quadratic equalities:
Composite optimization for robust phase retrieval

arXiv:1705.02356v1 [] 5 May 2017

John C. Duchi
Feng Ruan
Stanford University
May 9, 2017

Abstract
We develop procedures, based on minimization of the composition f (x) = h(c(x)) of a
convex function h and smooth function c, for solving random collections of quadratic equalities,
applying our methodology to real-valued phase retrieval problems. We show that the proxlinear algorithm we develop can solve (robust) phase retrieval problems (even with adversarially
faulty measurements) with high probability under appropriate random measurement models as
soon as the number of measurements m is a constant factor larger than the dimension n of
the signal to be recovered. The algorithm requires essentially no tuning—it requires solution
of a sequence of convex problems—and it is implementable without any particular assumptions
on the measurements taken. We provide substantial experiments investigating our methods,
indicating the practical effectiveness of the procedures and showing that they succeed with high
probability as soon as m/n ≥ 2.

1

Introduction

We wish to solve the following problem: we have a set of m vectors ai ∈ Rn and nonnegative scalars
bi ∈ R+ , i = 1, . . . , m, and wish to find a vector x ∈ Rn such that
bi = hai , xi2 for most i ∈ {1, . . . , m}.

(1)

As stated, this is a combinatorial problem that is, in the worst case, NP-hard [4]. Yet it naturally
arises in a number of real-world situations, including phase retrieval [18, 19, 21], in which one
receives measurements of the form
bi = hai , x? i2
for known measurement vectors ai ∈ Rn , while x? ∈ Rn is unknown. The problem in phase
retrieval arises due to limitations of optical sensors, where one illuminates an object x? , which
yields diffraction patter Ax? , but sensors may measure only the amplitudes b = (Ax? )2 , where (·)2
denotes elementwise squaring. See, for example, the paper [31] for a review of the phase retrieval
problem. In the case in which some measurements may be corrupted, the problem is even more
challenging.
A natural objective for the problem (1) is the exact penalty formulation [23], which replaces
the equality constraint bi = hai , xi2 with a non-differentiable cost measuring the error bi − hai , xi2 ,
yielding the formulation
m

minimize f (x) :=
x

1 X
1
| hai , xi2 − bi | =
(Ax)2 − b
m
m
i=1

1

1

.

(2)

This objective is a natural replacement of the equality constrained problem (1), and (as is wellknown in the statistics and optimization literature, where `1 -based losses yield median-based estimators), the `1 -loss handles gross errors in the measurements bi in a relatively benign way. Moreover,
in the case when bi = hai , x? i2 for all i, it is clear that ±x? globally minimize f (x), where we
note that it is only possible to recover x? up to a sign flip. Candès, Strohmer, and Voroninski
[7] and Eldar and Mendelson [17], as well as results we discuss later in the paper, show (roughly)
that f (x) stably identifies x? , in that it grows very quickly as x ± x? grows. The objective is,
unfortunately, non-smooth, non-convex—not even locally convex near x? , as is clear in the special
case f (x) = |x2 − 1|, so that a local analysis based on convexity is impossible—and at least f (x) a
priori seems difficult to minimize.
Nonetheless, the objective (2) enjoys a number of structural properties that, as we explore
below, make solving problem (1) tractable as long as the measurement vectors ai are sufficiently
random. In particular, we can write f as the composition f (x) = h(c(x)) of a convex function h
and smooth function c, a structure known in the optimization literature to be amenable to efficient
algorithms [20, 5, 15]. This compositional structure lends itself nicely to the prox-linear algorithm, a
variant of the Gauss-Newton procedure, which we describe briefly here. The composite optimization
problem, which Fletcher and Watson [20] originally develop and a number of researchers [5, 6, 15, 14]
have studied further, is to minimize
minimize f (x) := h(c(x)) subject to x ∈ X
x

(3)

where the function h : Rm → R is convex, c : Rn → Rm is smooth, and X is a convex set. It is clear
that this general form encompasses our objective (2), which we write as f (x) = h(c(x)) by taking
1
h(z) = m
kzk1 and c(x) = [hai , xi2 −bi ]m
i=1 . Using the common idea of most optimization schemes—
trust region, gradient descent, Newton’s method—to build a simpler to optimize local model of the
objective and repeatedly minimize this model, we can replace h(c(x)) in problem (3) by linearizing
only c. This immediately gives a convex surrogate and leads to the prox-linear algorithm developed
by Burke and Ferris [6, 5], among others [15, 13, 14]. Fixing x ∈ Rn , for any y ∈ Rn we define the
local “linearization” of f at x by

fx (y) := h c(x) + ∇c(x)T (y − x) ,
(4)
where ∇c(x) ∈ Rn×m denotes the Jacobian transpose of c at x. This is evidently convex in y, and
the prox-linear algorithm proceeds iteratively x1 , x2 , . . . by minimizing regularized models


1
2
xk+1 = argmin fxk (x) +
kx − xk k2 ,
(5)
2αk
x∈X
where αk > 0 is a stepsize (which may be fixed in advance or chosen by a line search [6, 15]).
1
If h is L-Lipschitz and ∇c is β-Lipschitz, then choosing any stepsize α ≤ βL
guarantees that the
method (5) is a descent method and finds approximate stationary points of the problem (3) [13, 15].
We briefly summarize our main contribution as follows. We show how to apply prox-linear
method (5) to any measurement matrix A with no tuning parameters except that the stepsize
1
AT A op )−1 . Each iteration requires solving a QP in n variables, which is effisatisfies α ≤ ( m
ciently solvable using standard convex programming approaches. We show that—with extremely
high probability under appropriate random measurement models—our prox-linear method converges quadratically as soon as the number of measurements m/n is greater than some numerical
constant, meaning we must solve only log2 log2 1 such convex problems to find an estimate x
b of
?
?
?
x such that min{kb
x − x k, kb
x + x k} ≤ . In practice, this is 5 convex quadratic programs. Our
procedure applies both in the noiseless setting and when a (constant but random) fraction of the
measurements are even adversarially corrupted.
2

1.1

Related work and approaches to phase retrieval

Our work should be viewed in the context of the recent and successful collection of work on phase
retrieval. A natural strategy for problem (1), when we wish to find x satisfying hai , xi2 = bi for all
i ∈ [m], is to lift the problem into a semidefinite program by setting X = xxT , relaxing the rank
one constraint, and solving
minimize tr(X) subject to X  0, tr(Xai aTi ) = bi .
X

This (and variants thereof) is the approach that a number of convex approaches to phase retrieval
take [11, 9, 7, 8, 38]. This convex optimization problem is computationally challenging for large
n, as it requires storing and manipulating an n × n matrix variable. Moreover, computation times
to achieve -accurate solutions to this problem generally scale as something like n3 /poly(), where
poly() denotes a polynomial in .
These difficulties have led a number of researchers to consider non-convex approaches to the
phase retrieval problem that—as we do—maintain only a vector x ∈ Rn , rather than forming a
full matrix X ∈ Rn×n . We necessarily give an only partial overview, focusing on recent work
giving schemes that are provably convergent. Early work in computational approaches to phase
retrieval is based on (non-convex) alternating projection approaches, notably those by Gerchberg
and Saxton [21] and Fienup [19]. Motivated by the challenges of convex approaches and the success
of alternating minimization [21, 19], Netrapalli et al. [27] develop an algorithm (AltMinPhase) that
alternates between minimizing kAx − Cbk2 in x, where C is a diagonal matrix of phases (signs), and
in C over the diagonal matrices where the diagonals have modulus one. Their algorithm is elegant,
but requires resampling a new measurement matrix A and measurements b in each iteration for
analysis. More recently, Candès et al. [10] develop Wirtinger flow, a gradient-based method that
performs a careful modification of gradient descent on the objective
m

1 X
F (x) :=
(|hai , xi|2 − bi )2 ,
2m
i=1

where x ∈ Cn may be complex. Wang et al. [39] build on this work by attacking a modification of
this objective, showing how to perform a generalized descent method on
m

p 2
1 X
|hai , xi| − bi ,
F (x) :=
2m
i=1

and providing arguments for the convergence of their method. Wang et al. achieve striking empirical
results, achieving better than 50% perfect signal recovery when the measurement ratio m/n = 2,
which is essentially at the threshold for injectivity of the measurements b = (Ax? )2 . Zhang et al.
[40] also study a variant of Wirtinger flow based on median estimates that handles some outliers.
Unfortunately, these procedures rely fairly strongly on Gaussianity assumptions, and their gradient
descent approaches require subsampling schemes (to select “good” terms in the sum), which have
parameters chosen carefully to reflect Gaussianity in the measurement matrices A; it is not always
clear how to extend the procedures to non-Gaussian measurements.

1.2

Our contributions and outline

In this paper, we focus on prox-linear methods, the iterations (5) for the non-smooth non-convex
problem (2). In addition to being (to us at least) aesthetically pleasing, as we minimize the natural
objective (2), our approach yields a number of theoretical and practical benefits.
3

In the literature on signal recovery from phaseless measurements, stability of the reconstruction
of a signal is of paramount importance. To solve the phase retrieval problem at all, one requires
intectivity of the measurements b = (Ax)2 , which for real A ∈ Rm×n in general position necessitates
m ≥ 2n − 1 (cf. [1]). Stability makes this injectivity more robust: Eldar and Mendelson [17] say
that a measurement matrix A ∈ Rm×n is λ ≥ 0 stable if
(Ax)2 − (Ay)2

1

≥ λ kx − yk2 kx + yk2 for all x, y ∈ Rn .

(6)

Such conditions, which hold with high probability for suitable designs A (e.g. Gaussian or isotropic
sub-Gaussian [7, 17]), are also common in semidefinite relaxation approaches to phase retrieval;
cf. Candès et al. [7, Lemma 3.2]. (See also the paper [2].) This condition means that distant signals
x, x0 cannot be confused in the measurement domain {(Ay)2 | y ∈ Rn } ⊂ Rm
+ because A does a good
job of separating them; the more stable a measurement matrix, the “easier” the recovery problem
should be. We provide stability guarantees for a much more general class of random matrices
using variants of Mendelson’s “small ball” techniques [26] (Sec. 3.1). Most literature on non-convex
approaches to phase retrieval requires such a stability condition—and usually more because of the
quadratic objectives often used—to guarantee signal recovery. In contrast, our procedure requires
1
essentially only the stability condition (6), a mild bound on the operator norm m
|||A|||2op of A, and
an initialization within some constant factor of ±x? to guarantee both fast convergence and exact
signal recovery.
With this in mind, in Section 2 we develop purely optimization-based deterministic results,
which build off of a number of researchers’ works on composite optimization, and that rely on the
stability condition (6). By identifying the conditions required for fast convergence and recovery, we
can then spend the remainder of the paper showing how various measurement models guarantee
sufficient conditions for our convergence results. In particular, in Section 3, we show how a number
of sensing matrices A suffice to guarantee convergence and signal recovery in the noiseless setting,
that is, when b = (Ax? )2 . In Section 4, we extend these results to the case when a constant
fraction of the measurements bi may be arbitrarily corrupted, showing that stability and a somewhat
stronger condition on |||A|||op are still sufficient to guarantee signal recovery; again, these results hold
for our basic algorithm with no tuning parameters.
In the final sections of the paper, we provide a substantial empirical evaluation of our proposed
algorithms. While our method in principle requires no tuning—it solves a sequence of explicit convex problems—we confess that there is some art in developing efficient methods for the solution of
the sequence of convex optimization problems we solve. In Section 5, we describe these implementation details. Finally, in Section 6 we provide experimental evidence of the success of our proposed
approach. In reasonably high-dimensional settings (n ≥ 1000), with random Gaussian measurements our method achieves perfect signal recovery in about 80-90% of cases even when m/n = 2,
with improvements for larger n. The method also handles outlying measurements well, substantially improving state-of-the-art performance, and we give applications with measurement matrices
that demonstrably fail all of our conditions, but for which the method is still straightforward to
implement and empirically very successful.
Notation We collect our common notation here. We let k·k and k·k2 denote the usual vector
`2 -norm, and for a matrix A ∈ Rm×n , |||A|||op denotes its `2 -operator norm. Let quantα ({ci }) denote
the α-quantile of a vector c ∈ Rm , that is, if c(1) ≤ c(2) ≤ · · · ≤ c(m) , the αth quantile linearly
interpolates c(bmαc) and c(dmαe) . For a random variable X, quantα (X) denotes its αth quantile.

4

Data: Initializer x0 , stepsize α, tolerance  > 0, convex h : Rm → R, smooth c : Rn → Rm
set k = 0
repeat
1
xk+1 = argminx∈X {fxk (x) + 2α
kx − xk k22 } where fx (y) = h(c(x) + ∇c(x)T (y − x))
k =k+1
until kxk−1 − xk k2 ≥ α;
return xk
Algorithm 1: Prox-linear algorithm for problem (3)

2

Composite Optimization, Algorithm and Convergence Analysis

We begin our development by providing convergence guarantees—under appropriate conditions—
for the prox-linear algorithm (the iteration (5)) applied to the composite optimization problem (3),
which re recall is to
minimize f (x) := h(c(x)) subject to x ∈ X,
x

Rm

where h :
→ R is convex and c : Rn → Rm is smooth. We summarize the algorithm in Alg. 1 for
1
further reference. In our application to quadratic constraints and phase retrieval, h(z) = m
kzk1
and c(x) = (Ax)2 − b, so that the iteration (5) is the solution of a quadratic problem.
A number of researchers have studied convergence and stopping conditions for Algorithm 1,
showing that it converges to stationary points [5], as well as demonstrating that the stopping
condition kxk − xk−1 k2 ≤ α holds after O(−2 ) iterations and guarantees approximate stationarity [13, 15]. Algorithm 1 and the iteration (5) sometimes enjoy fast (local) convergence rates
as well. To describe this phenomenon, we say that h has weak sharp minima if it grows linearly
away from its minima, meaning h(z) ≥ inf z h(z) + λ dist(z, argmin h) for some λ > 0. Under this
condition, Burke and Ferris [6] show that convergence of the prox-linear algorithm near points
in X ? := {x : c(x) ∈ argmin h} is quadratic, because the model (4) of f is quadratically good,
but h(c(x)) grows linearly away from X ? . We build from this elegant development—though our
problems do not actually satisfy the weak sharp minima conditions because of outliers—to show
how the prox-linear algorithm provides an effective, easy-to-implement, and elegant approach to
problems involving solution of quadratic equalities, specifically focusing on phase retrieval.

2.1

Quadratic convergence and the prox-linear method for phase retrieval

We turn now to an analysis of the prox-linear algorithm for phase retrieval problems, providing
conditions on the function f sufficient for quadratic convergence. We introduce two conditions on
the function f (x) and its linearized form fx (y) that suffice for this desidaratum; as we show in the
sequel, these conditions hold with extremely high probability for a number of random measurement
models. These conditions are the keystones of our analysis of the (robust) phase retrieval problem.
As motivation for our first condition, recall Eldar and Mendelson’s stability condition (6). If
the measurement matrix A satisfies condition (6) and the measurements bi are noiseless, so that
b = (Ax? )2 , then evidently we have f (x) − f (x? ) ≥ λ kx − x? k2 kx + x? k2 . When the measurements
have noise or outliers, this may still hold, prompting us to define the following
Condition C1. There exists a parameter λ > 0 such that for all x ∈ Rn we have
f (x) − f (x? ) ≥ λ kx − x? k kx + x? k .
This condition is a close cousin of Burke and Ferris’s sharp minima condition [6], though it does
not require that c(x? ) ∈ argminz h(z); based on their work, it is intuitive that it should prove
5

useful in establishing fast convergence of the prox-linear algorithm. The second condition, which
is essentially automatically satisfied for the linear approximation (4), is a requirement that the
linearized function fx (y) is quadratically close to f (y).
Condition C2. There exists a parameter L < ∞ such that for all x, y ∈ Rn
|f (y) − fx (y)| ≤

L
kx − yk22 .
2

Locally, Condition C2 holds for any composition f (x) = h(c(x)) of a convex h with smooth c,
but the phase retrieval objective (2) satisfies the bound globally. Indeed, for y ∈ Rn we have
hai , yi2 = hai , xi2 + 2 hai , xi hai , y − xi + hai , y − xi2 ,
so the linearization (4) of f around x is
m

fx (y) =

1 X
| hai , xi2 − bi + 2 hai , xi hai , y − xi |.
m

(7)

i=1

]T

Rm×n

Letting A = [a1 · · · am ∈
denote the measurement matrix, then using the preceding
2
expansion of hai , yi , we have immediately by the triangle inequality that


m
1 X
1 T
T
A A (x − y) + fx (y) ≤
| hai , yi2 − bi | = f (y)
−(x − y)
m
m
i=1


1 T
T
≤ fx (y) + (x − y)
A A (x − y).
m
1 T
That is, Condition C2 holds with L = 2||| m
A A|||op :

1
AT A op kx − yk22 .
(8)
m
Given Conditions C1 and C2, we now turn to convergence guarantees for the prox-linear Algorithm 1, which in our case requires solving a sequence of convex quadratic programs. An implementation of Alg. 1 that solves iteration (5) exactly may be computationally challenging. Thus, we
allow inaccuracy in the solutions, assuming there exists a sequence of additive accuracy parameters
k ≥ 0 such that the iterates xk satisfy


L
L
2
2
fxk (xk+1 ) + kxk+1 − xk k2 ≤ inf fxk (x) + kx − xk k2 + k .
(9)
x
2
2
|f (y) − fx (y)| ≤

Defining the symmetrized distance measure dist(x, y) := min{kx + yk , kx − yk}, we have the following theorem, whose proof we provide in Section 2.2.
Theorem 1. Let Conditions C1 and C2 hold. Assume that in each step of Algorithm 1, we
solve the intermediate optimization problem to accuracy k , and define the relative error measures
2k
βk = λkx
. Then
? k2
2

dist(xk , x? )
λ
≤
max
?
kx k2
2L

(

2L dist(x0 , x? )
·
λ
kx? k2

2k


, max

0≤j<k

4L
· βj
λ

If k = 0 in the solution quality inequality (9), then

k
dist(xk , x? )
λ L dist(x0 , x? ) 2
≤
·
.
kx? k2
L λ
kx? k2
6

2k−j−1 )
.

Theorem 1 motivates our approach for the remainder of the paper: we can guarantee exact,
accurate, and fast solutions to the phase retrieval problem under the three conditions
1. Stability (Condition C1),
2. Quadratic approximation (Condition C2), via an upper bound on
inequality (8) and

AT A

op

and application of

3. An initializer x0 of the iterations that is good enough, meaning that it satisfies the constant
relative error guarantee dist(x0 , x? ) ≤ kx? k2 Lλ .
In the coming sections, we show that each of these three conditions holds in both noiseless measurement models (Section 3) and with adversarially perturbed measurements bi (Section 4).
Before continuing with this approach, we provide a few brief remarks on Theorem 1. First, if
k are very small because we solve the intermediate steps (5) to (near) machine precision, then
for all intents and purposes about five iterations suffice for machine precision accurate solutions.
Quadratic convergence is also achievable with errors in inequality (9); if the minimization accuracies
k
decrease quickly enough that k ≤ 2−2 , then we certainly still have quadratic convergence. More
broadly, Theorem 1 shows that the accuracy of solution in iteration j need not be very high to
guarantee high accuracy reconstruction of the signal x? ; only in the last few iterations is moderate
to high accuracy necessary. If it is computationally cheap, it is thus advantageous—as we explore
in our experimental work—to solve early iterations of the prox-linear method inaccurately.

2.2

Proof of Theorem 1

We prove the result in two steps: we first provide a per-iteration progress guarantee, and then we
use this guarantee to show quadratic convergence.
The function x 7→ fxk (x) + L2 kx − xk k22 is L-strongly convex in x. If we define x?k+1 to be the
exact minimizer of fxk (x) + L2 kx − xk k22 , the standard optimality conditions for strongly convex
minimization then imply
fxk (xk+1 ) +

L
L ?
2
− xk 2 + k
kxk+1 − xk k22 ≤ fxk (x?k+1 ) +
x
2
2 k+1
L
L ?
≤ fxk (x? ) + kx? − xk k22 −
x − x?k+1
2
2

2
2

+ k .

Using the approximation Condition C2, so that L2 kxk − xk+1 k22 +fxk (xk+1 ) ≥ f (xk+1 ) and fxk (x? ) ≤
f (x? ) + L2 kxk − x? k22 , we have by substituting in the preceding inequality that
L ?
L ?
2
kx − xk k22 −
x − x?k+1 2 + k
2
2
L ?
2
≤ f (x? ) + L kxk − x? k22 −
x − x?k+1 2 + k .
2

f (xk+1 ) ≤ fxk (x? ) +

Performing a parallel derivation with −xk replacing xk and simplifying, we have
L ?
x − x?k+1
2
L ?
f (xk+1 ) − f (x? ) +
x + x?k+1
2

f (xk+1 ) − f (x? ) +

7

2
2

≤ L kxk − x? k22 + k

2
2

≤ L kxk + x? k22 + k .

By applying the stability Condition C1, we immediately obtain the progress guarantee
L ?
2
x
− x? 2 ≤ L kxk − x? k22 + k
2 k+1
(10)
L ?
?
?
? 2
? 2
λ kxk+1 − x k2 kxk+1 + x k2 +
x
+ x 2 ≤ L kxk + x k2 + k .
2 k+1
We now transform the guarantee (10) into one involving only xk , xk+1 , and x? , rather than x?k+1 ,
by bounding the difference between x?k+1 and xk+1 . The L-strong convexity of fxk (·) + L2 k· − xk k22
implies that
λ kxk+1 − x? k2 kxk+1 + x? k2 +

k + fxk (x?k+1 ) +

whence we obtain
2 xk+1 − x?k+1

2
2

L
2

L ?
x
− xk
2 k+1

x?k+1 − xk+1

+ 2 x?k+1 − x?

2
2

2
,
2

2
2

L
kxk+1 − xk k22
2
L ?
L
2
≥ fxk (x?k+1 ) +
xk+1 − xk 2 +
xk+1 − x?k+1
2
2
≥ fxk (xk+1 ) +

2
,
2

≤ k . Using the standard quadratic inequality kxk+1 − x? k22 ≤
we thus have by expression (10) that

L
kxk+1 − x? k22 ≤ L kxk − x? k22 + 2k
4
(11)
L
?
?
? 2
? 2
λ kxk+1 − x k2 kxk+1 + x k2 + kxk+1 + x k2 ≤ L kxk + x k2 + 2k .
4
Now, we note that by the triangle inequality, we have 2 kx? k2 ≤ kxk − x? k2 + kxk + x? k2 , so
λ kxk+1 − x? k2 kxk+1 + x? k2 +

kxk − x? k2 kxk + x? k2 = dist(xk , x? ) max{kxk − x? k2 , kxk + x? k2 } ≥ kx? k2 dist(xk , x? ).
Substituting into inequality (11), we find that
λ kx? k2 dist(xk+1 , x? ) ≤ L dist(xk , x? )2 + 2k .
Dividing each side by λ kx? k22 yields
dist(xk+1 , x? )
L dist(xk , x? )2
2k
≤
+
.
2
?
?
kx k2
λ
kx k2
λ kx? k22

(12)

Inductively applying inequality (12) when k = 0 yields the second statement of the theorem.
When k > 0, we require a brief technical lemma, which shows the convergence rate of sequences
satisfying inequalities like that in (12):
Lemma 2.1. Let the sequences ak ≥ 0, k ≥ 0 satisfy ak ≤ κa2k−1 + k−1 . Then


2k −1 2k
2k−j−1 −1
2k−j−1
ak ≤ max (2κ)
a0 , max (2κ)
(2j )
.
0≤j<k

Proof The proof is by induction. For a1 , we certainly have a1 ≤ 2κ ∨ 20 because both sequences
are non-negative. For the general case, assume the result holds for ak , where k is arbitrary. Then

ak+1 ≤ κa2k + k ≤ max 2κa2k , 2k
n
o
k+1
k+1
k−j
k−j
≤ 2κ(2κ)2 −2 a02
∨ max 2κ(2κ)2 −2 (2j )2
∨ 2k
0≤j<k


2k+1 −1 2k+1
2k−j −1
2k−j
= max (2κ)
a0 , max (2κ)
(2j )
0≤j<k+1

as desired.
Applying Lemma 2.1 in inequality (12) with κ =
8

L
λ

and ak =

dist(xk ,x? )
kx? k2

yields the theorem.

3

Noiseless Phase Retrieval Problem

We begin our discussion of the phase retrieval problem by considering the noiseless case, that is,
when the observations bi = hai , x? i2 . Based on Theorem 1, to show that the method success it is
sufficient to show the following three results: (i) the objective function (2) satisfies the stability
condition C1; (ii) the objective (2) satisfies the quadratic approximation condition C2 by upper
bounding |||AT A|||op , as in inequality (8); and (iii) that we have a good initializer x0 . In the coming
three sections, we address each of these in turn, providing progressively stronger assumptions that
are sufficient for each condition to hold with high probability as soon as the number of measurements
m/n > c, where c is a numerical constant. In Section 3.4 we provide a summary theorem that
encapsulates our results. For readability, we defer proofs to Sec. A.

3.1

Stability

Our first step in showing the success of the composite optimization approach to phase retrieval, in
the noiseless case, is to provide conditions under which stability holds with high probability. With
this in mind, we make the following assumption.
Assumption A1. There exists a function κst : Rn × Rn → R+ and a constant p0 > 0 such that
for all u, v ∈ Rn , we have
P (| ha, ui ha, vi | ≥ κst (u, v)) ≥ p0 ,
where
κ?st :=

inf

u6=0,v6=0

κst (u, v)
> 0.
kuk2 kvk2

Intuitively, Assumption A1 says that the measurement vectors a ∈ Rn have sufficient support in
all directions u, v ∈ Rn . As an example, isotropic Gaussian vectors satisfy Assumption A1.
Example 1 (Gaussian stability): Let a ∼ N(0, In ). We claim that


1
1
P | ha, ui ha, vi | ≥ q ≥
for all u, v ∈ Sn−1 ,
(13)
2
4
where q ≥ .73 is the median of |X 2 − Y 2 | for X, Y are independent standard normal. That is,
standard Gaussians satisfy Assumption A1 with κst (u, v) = .365 kuk2 kvk2 and p0 = 41 .
Indeed the random variables Zu = ha, ui and Zv = ha, vi are bivariate normal with Var(Z) = 1
and Cov(Zu , Zv ) = hu, vi. Letting (X, Y ) ∼ N(0, I2 ), we have for any set C ⊂ R that
inf

u,v∈Sn−1

P (|Zu Zv | ∈ C) =

P (|(cos(θ)X + sin(θ)Y )(cos(θ)X − sin(θ)Y )| ∈ C)

= inf P |ρX 2 − (1 − ρ)Y 2 | ∈ C .
inf

θ∈[0,π/2]

ρ∈[0,1]

For any c > 0 and any ρ ∈ [0, 1], the triangle inequality implies



P X 2 − Y 2 ≥ 2c ≤ P ρX 2 − (1 − ρ)Y 2 ≥ c + P (1 − ρ)X 2 − ρY 2 ≥ c

= 2P ρX 2 − (1 − ρ)Y 2 ≥ c .
Choosing c = 12 med(X 2 − Y 2 ) yields the claim (13). ♣
Assumption A1, however, does not require any type of light tails: just that the probability of
the magnitude of ha, ui ha, vi being large is a constant. By extending Example 1, we see that a
similar conclusion holds for random vectors a = Z kZkk2 , where Z ∼ N(0, In ) and k ≥ 1 is arbitrary.
9

As we note in the discussion precedindg Condition C1, for the objective (2) it is immediate that
in the noiseless case, we have
1
f (x) − f (±x ) = f (x) =
(Ax)2 − (Ax? )2
m

m

?

1

1 X
=
| hai , x − x? i hai , x + x? i |.
m
i=1

Thus,
Pm if we can show the stability condition (6) of Eldar and Mendelson [17]—equivalently, that
i=1 |hai , uihai , vi| ≥ λm—then the data matrix A guarantees that Condition C1 holds. To that
end, we provide the following guarantee, which we prove in Sec. A.1.
Proposition 1. Let Assumption A1 hold. There exists a numerical constant c < ∞ such that for
any t ≥ 0,
!
r


m
√
1 X
n
| hai , ui hai , vi | ≥ κ?st p0 − c
− 2t for all u, v ∈ Sn−1 ≥ 1 − 2e−mt .
P
m
m
i=1

Proposition 1 immediately yields the following corollary, which shows that the stability condition C1 holds with high probability for m/n & 1.
Corollary 3.1. Let Assumption A1 hold. Then there exists a numerical constant c < ∞ such that
if mp20 ≥ cn, then




1 ?
mp20
?
?
?
n
P f (x) − f (x ) ≥ κst p0 kx − x k2 kx + x k2 for all x ∈ R
.
≥ 1 − 2 exp −
2
32

3.2

Quadratic Approximation

With the stability condition C1 in place, we turn to a discussion of the approximation condition C2.
As implied by the estimate in inequality (8), the quadratic approximation condition is satisfied
1 T
with parameter L = 2||| m
A A|||op . To control this quantity, we require that the rows of the matrix
m×n
A∈R
be sufficiently light-tailed. One standard assumption is sub-Gaussianity of the ai , which
we define in terms of Orlicz norms (following [37, 36, Ch. 2.2]).
Definition 3.1. The random vector a ∈ Rn is σ 2 -sub-Gaussian if for all v ∈ Sn−1
"

#
ha, vi2
E exp
≤ e.
σ2
Assumption A2. The measurement vectors ai are independent and σ 2 -sub-Gaussian.
It is of course possible to bound |||A|||op when the rows have heavier tails using appropriate symmetrization techniques and matrix Khintchine inequalities (cf. [37, Sec. 2.6]), but we choose not to
address such issues here.
2
Assumption A2 is not satisfied by all choices of measurement vectors ai . Using that E[eλZ ] =
−1

2

[1 − 2λ]+ 2 for Z ∼ N(0, 1), we see that it holds for ai ∼ N(0, In ) with σ 2 = e2e
2 −1 ≈ 2.313 and for ai
1
n−1
2
uniform on S
with σ = O(1) · n . In practice, however, it may be useful to apply our algorithm
2 m
to the transformed data {ai / kai k2 }m
i=1 and {bi / kai k2 }i=1 , which (in the noiseless case or case with
infrequent but arbitrary corruptions of bi ) is likely to make the problem better conditioned. There
are two heuristic motivations for this: first, those measurement vectors with larger magnitudes kai k
tend to place a higher weight in the optimization problem (2), and thus normalization can make
10

the observations “comparable” to each other; second, normalization guarantees the measurement
vectors {ai }m
i=1 satisfy Assumption A2, yielding easier verification of Condition C2. (It may be
more challenging to verify Assumption A1, but if the ai are sufficiently isotropic this presents no
special difficulties.)
Standard results guarantee that the random matrices A have well-behaved singular vectors
whenever Assumption A2 holds; we provide one standard such result due to Vershynin [37, Thm. 39,
Eq. (25)] here, with constants that are achievable by tracing his proof.
Lemma 3.1. Let Assumption A2 hold and Σ = E[aaT ]. Then for all t ≥ 0,
(r
)!
1 T
4n
4n
P
A A−Σ
≥ 11σ 2 max
+ t,
+t
≤ exp(−mt).
m
m
m
op
Moreover, |||Σ|||op ≤ σ 2 and E[|ha, vi|k ] ≤ Γ( k2 + 1)eσ k for all k ≥ 0.
Thus, we have the following corollary of Lemma 3.1, which guarantees that Condition C2 holds
with high probability for m/n & 1.
Corollary 3.2. Let Assumption A2 hold. Then there exists a numerical constant c < ∞ such that
whenever m ≥ cn

 m

P |fx (y) − f (y)| ≤ 2σ 2 kx − yk22 holds uniformly for x, y ∈ Rn ≥ 1 − exp −
.
c
p
Proof Assume m > 4n, and choose t small enough in Lemma 3.1 that 11σ 2 4n/m + t ≤ σ 2 .
1 T
1 T
Applying the triangle inequality to ||| m
A A − Σ|||op ≤ ||| m
A A|||op + σ 2 then yields the result by
equation (8).

3.3

Initialization

The last ingredient in achieving strong convergence guarantees for our prox-linear procedure for
phase retrieval is to provide a good initialization. There are a number of initialization strategies in
the literature [10, 39, 40] based on spectral techniques, which work as follows. First, we decompose
the initialization into two ingredients: we (i) find an estimate of the direction direction d? :=
x? / kx? k2 , and (ii) estimateP
the magnitude r? := kx? k2 . The latter is easy: assuming that E[ai aTi ] =
2
In , one simply uses rb2 = m
i=1 bi , which is unbiased and tightly concentrated. The former, the
direction estimate, is somewhat trickier.
Wang, Giannakis, and Eldar [39] provide an empirically excellent initialization whose heuristic
justification is as follows. First, for random vectors ai in high dimensions, we expect ai to usually
be orthogonal to the direction d? . Thus, by extracting the smallest magnitude bi = hai , x? i2 , we
have the vectors ai that are “most” orthogonal to the direction d? ; letting Isel be these small
indices, the eigenvector corresponding
to the smallest eigenvalue (for simplicity, we simply call
P
T
?
this the smallest eigenvector) of
i∈Isel ai ai should be close to the direction d . A variant of
P
m
1
T
this procedure P
is to note that m
i=1 ai ai ≈ In when the ai are isotropic, so that the largest
eigenvector of i6∈Isel ai aTi should also be close to d? . This initalization strategy has the added
benefit that—unlike the original
initialization schemes developed by Candès et al. [10],
P spectral
T that may not concentrate at sub-Gaussian rates (as the
which rely on eigenvectors of m
b
a
a
i=1 i i i
P
sum involves fourth moments of random vectors)—the sums i ai aTi are tightly concentrated.

11

Data: Measurement matrix A ∈ Rm×n and signals b = (Ax? )2
begin
1 Pm
Set rb2 = m
i=1 bi
Select indices Isel := {i ∈ [m] : bi ≤ 21 rb2 }
Construct directional estimates db by
X
X init :=
ai aTi , db = argmin dT X init d
d∈Sn−1

i∈Isel

b
return (b
r, d)
end
Algorithm 2: Initialization procedure for non-noisy data
Unfortunately, we believe Wang et al.’s proof that this initialization works contains a mistake:
letting U ∈ Rn×n−1 be an orthogonal matrix whose columns P
are all orthogonal to d? , in the proof
T
c |−1
of Lemma 2 (Eqs. (68)–(70) in [39]) they assert that |Isel
c U ai ai U → In−1 when the ai
i∈Isel
√ n−1
are uniform on nS
. This is not true (nor does appropriate normalization by n or n − 1 make
it true), as it ignores the subtle effects of conditioning in the construction of Isel . (One can see
this by a quick simulation with d? = e1 , the first standard basis vector.) In spite of this issue, the
initialization they propose works remarkably well, and as we show presently, it provably provides
a good estimate db of d? . We include the initialization procedure in Algorithm 2.
With the previous discussion in mind, we provide a general assumption that is sufficient for
Alg. 2 to return a direction and magnitude estimate sufficiently accurate for phase retrieval.
Assumption A3. For some 0 ∈ [0, 1] and p0 (d? ) > 0, the following hold.
(i) For all  ∈ (0, 0 ], the following continuity and directional likelihood conditions hold:





1− 1+
1−
? 2
? 2
P ha, d i ∈
,
≤ κ and P ha, d i ≤
≥ p0 (d? ) > 0.
2
2
2
(ii) There exist functions φ : [0, 0 ] → R+ and ∆ : [0, 0 ] → Rn×n such that


1−
T
? 2
= In − φ()d? (d? )T + ∆() for  ∈ [0, 0 ].
E aa | ha, d i ≤
2
(iii) E[aaT ] = In .
Assumption A3 on its face seems fairly technically complex. Nonetheless, each of its components
is not too stringent. Part (i) essentially says that ha, d? i2 has no point mass at ha, d? i2 = 21 and
that ha, d? i2 has reasonable probability of being smaller than 12 . Part (iii) simply states that in
expectation, a is isotropic (and a rescaling of a can guarantee this). Part (ii) is the most subtle
and essential for our derivation; it says that a ∈ Rn is reasonably isotropic, even if we condition
on ha, d? i being near zero for some direction d? , so that most mass of aaT is distributed uniformly
in the orthogonal directions In − d? d? T . The error terms φ and ∆ allow non-trivial latitude in
this condition, so that Assumption A3 holds for more than just Gaussian vectors. That said, for
concreteness we provide the following example.
iid
Example 1 (continued): We continue our running example for ai ∼ N(0, In ), showing that such
ai satisfy Assumption A3 for any 0 ∈ (0, 1) with residual error ∆ ≡ 0. Clearly Part (iii) holds.
12

1

For Part (i), note that ha, d? i2 is χ21 -distributed with density f (t) = e−t/2
(2πt)− 2 . Integrating the
p
0 −1
?
0
density using its upper
bound, we may set κ ≤ f ( 1−
2 ) = exp( 4 )/ π(1 − 0 ) and p0 (d ) =
1−0
1√
2
P(χ1 ≤ 2 ) ≥ 2 1 − 0 . Part (ii) of Assumption A3 is all that remains. By the rotational
invariance of a ∼ N(0, In ), we see for any t ∈ R+ that
E[aaT | ha, d? i2 ≤ t] = In − d? d? T + E[ha, d? i2 | ha, d? i2 ≤ t]d? d? T .
We claim the following lemma, whose proof we provide in Appendix C.1.
Lemma 3.2. Let Z be a random variable with a density symmetric about zero, where the density
2
is decreasing on R+ . Then for any c ∈ R we have E[Z 2 | Z 2 ≤ c2 ] ≤ c3 .
Thus, setting φ() = 1 − E[ha, d? i2 | ha, d? i2 ≤ 1−
2 ] ≥ 1−
Assumption A3 holds with φ() ≥ 65 and ∆() ≡ 0. ♣

(1−)2
6

≥

5
6,

we see that part (ii) of

We now state our main proposition of this section.
Proposition 2. Let Assumptions A2 and A3 hold and let 0 and p0 (d? ) be as in Assumption A3.
Define the error measure
"
#
1
4e(1 + log 1∧κ
)σ 2 κ 2κ(1 + )
ν() := 1 +
+
.
p0 (d? )
p0 (d? )2
b be the estimated
There exists a numerical constant c > 0 such that the following holds. Let (b
r, d)
b For any  ∈ [0, 0 ], if
magnitude and direction of Alg. 2 and define x0 = rbd.
c·
then

σ 4 log2 p0 (d? )
1
m
≥
∨
?
2
n
p0 (d )
(κ)2

|||∆()|||op + ν()
dist(x0 , x? ) p
i +
≤ 2(1 + ) h
?
kx k2
φ() − |||∆()||| − ν()
op

(14)

+

with probability at least







mp0 (d? )2
cmκ2
cm2
2 2
1 − exp − 4
− 2 exp −cm κ − exp −
− exp − 4
.
σ
2
σ log2 p0 (d? )
We prove Proposition 2 in two parts. In the first part (Sec. A.2), we define a number of events and
proceed conditionally, showing that if each of the events occurs then the conclusion (14) holds. In
the second part (Sec. A.3) we show that the events occur with high probability.
We provide a few remarks to make the result clearer. Let us make the simplifying assumptions
that the constants in Assumption A3 are absolute constants (which is certainly satisfied for Gaussian
measurement vectors), that is, that that σ 2 = O(1), p0 (d? ) = Ω(1), and κ = Ω(1). Then for
numerical constants c > 0, C < ∞ we have for any  ∈ [0, 0 ] that
ν() ≤ C  log

|||∆()|||op + ν()
m
C
dist(x0 , x? )
1
i +
and
≥ 2 implies
≤ Ch
?

n

kx k2
φ() − |||∆()||| − ν()
op

+

with probability at least 1 − 5 exp(−cm2 ). Here, we see three competing terms. The first two,
the separation φ() and error ∆(), arise from the conditional expectation of Assumption A3(ii),
13

? ? T + ∆(). This is intuitive: the larger the
where we have E[aaT | ha, d? i2 ≤ 1−
2 ] = In − φ()d d
separation φ() from uniformity in the conditional expectation of aaT , the easier it should be for
spectral initialization to succeed; larger error ∆() will hide the directional signal d? . The last term
is the error ν() .  log 1 , which approaches 0 nearly as quickly as  as  → 0.
In the case that the error ∆() = 0 and gap φ() is bounded away from zero, which holds for
elliptical distributions with identity covariance—the Gaussian distribution and uniform distribution
√
−2 we have relative
on nSn−1 being the primary examples—we thus see that as soon as m
n & 
error
dist(x0 , x? )
1
.  log
with probability ≥ 1 − 5 exp(−cm2 ).
(15)
?
kx k2


That is, we can construct an arbitrarily good initialization with large enough sample size. (This
proves that the initialization scheme of Wang et al. [39] also succeeds with high probability.) On the
other hand, when ∆() 6= 0 for all  ∈ [0, 0 ], then Proposition 2 cannot guarantee arbitrarily good
initialization: the error term |||∆()|||op is never zero. However, if it is small enough, we still achieve
initializers that are within constant relative distance of x? , which is good enough for Theorem 1.

3.4

Summary and success guarantees

We have now provided guarantees of stability, quadratic approximation, and good initialization
for appropriate measurement matrices A ∈ Rm×n , at least when the observations b = (Ax? )2 are
noiseless. We provide a summary theorem showing that the composite optimization procedure we
propose works as soon as the sample size is large enough. In stating the theorem, we assume that
each of Assumptions A1, A2, and A3 holds with all of their constants actually numerical constants.
That is, in Assumption A1, we assume that κ?st & 1 and p0 & 1; in Assumption A2 we assume that
σ 2 . 1; and in Assumption A3 we assume that κ . 1 and p0 (d? ) & 1. The one somewhat technical
assumption we require is that relating the sub-Gaussian parameter σ 2 , the stability parameters κ?st
and p0 , and the error ∆() and directional separation constants φ(). In particular, we assume that
for a suitably small numerical constant c > 0, we have
|||∆()|||op ≤ c

κ?st p0
,
σ2

φ() ≥

1
|||∆()|||op , and φ() ≥ φ? > 0
2

for all  ∈ [0, 0 ]. We then have the following theorem, which follows by combining our convergence
Theorem 1 with Proposition 1, Corollary 3.2, and Proposition 2.
Theorem 2. Let the conditions of the preceding paragraph hold. There exists a numerical constant
2
n
c > 0 such that if m
< c, then (i) the initializer x0 returned by Alg. 2 satisfies dist(x0 , x? ) ≤ 8κσ? p0
st
and (ii) assuming no error in the minimization steps of the prox-linear method (Alg. 1),
k

dist(xk , x? ) ≤ kx? k2 2−2
for all k ∈ N with probability at least 1 − 5e−cm .

4

Phase retrieval with outliers

The objective (2) is the analogue of the least-absolute deviation estimator of location—the median
in R—so in analogy with our rich understanding of robustness [24], it is natural to expect it should
be robust to outliers. We turn to this question now, showing that this is the case, and the proxlinear method we have developed is quite effective. We consider the following corruption model: we
14

let {ξi } ⊂ R be an arbitrary (potentially random) sequence, and given the m measurement vectors
ai , we observe
(
hai , x? i2 if i ∈ I in
bi =
ξi
if i ∈ I out ,
where I out ⊂ [m] and I in ⊂ [m] denote the outliers and
respectively. We assume there is
 inliers,

a pre-specified measurement failure probability pfail ∈ 0, 21 , and |I out | = pfail m, and the actual
indices i ∈ I out are chosen randomly. That is, measurement failures are random, though the noise
sequence ξi may depend on ai (even adversarially), as we specify presently. We assume no prior
knowledge of which indices i ∈ [m] actually satisfy i ∈ I out , or even of pfail .
We consider the following two models for errors:
m
Model M1. The measurement vectors {ai }m
i=1 are independent of the all the values {ξi }i=1 .

Model M2. The inlying measurement vectors {ai }i∈I in are independent of the values {ξi }i∈I out of
the corrupted observations.
Model M1 requires full independence between the noise and measurements: the adversary may
only corrupt ξi without observing ai . Model M2 relaxes this requirement, allowing completely
arbitrary dependence between the corrupted data and the measurement vectors ai for i ∈ I out (i.e.
the corrupted indices). This is natural, as the type of corruption may depend on the individual
measurement ai being taken.
The arbitrary corruption causes some technical challenges, but we may still follow the outline
established in our analysis of phase retrieval without noise in Sec. 3. As we show in Section 4.1,
the objective f (x) is still stable (Condition C1) as long as the measurement vectors are light-tailed,
though Gaussianity is unnecessary. The quadratic approximation conditions (Condition C2) are
completely identical to those in Sec. 3.2, so we ignore them. Thus, as long as |||A|||op is not too
large (meaning fx (y) ≈ f (y)) and we can find a good initializer, the prox-linear iterations (5) will
converge quadratically to x? . Finding a good initializer x0 is somewhat trickier, but in Section 4.2
we provide a spectral method, inspired by Wang et al. [39], that works with high probability as
soon as m/n ≥ C for some numerical constant C. We defer our arguments, which are somewhat
long and technical to Appendix B.

4.1

Stability

The outlying indices, even when corruptions are chosen adversarially, have limited effect on the
1
growth and identification behavior of f (x) = m
(Ax)2 − b 1 . In particular, for pfail smaller than
a numerical constant, which we can often specify, the stability condition C1 holds with high probability whenever m/n is large. More precisely, we have the following proposition, which applies to
any σ 2 -sub-Gaussian vector or measurement matrix A with independent rows.
Proposition 3. Let Assumption A2 hold and κst = inf u,v∈Sn−1 E[| ha, ui ha, vi |]. Then under either
of the models M1 or M2, there are numerical constants c > 0 and C < ∞ such that
r


?
?
2 3 n
2
f (x) − f (x ) ≥ κst − 2pfail − Cσ
− Cσ t kx − x? k2 kx + x? k2 for all x ∈ Rn
m
2

with probability at least 1 − 2e−cm − 2e−mt .

15

Data: Measurement matrix A ∈ Rm×n and corrupted signals b
begin
Define indices Isel := {i ∈ [m] : bi ≤ med ({bi }m
i=1 )}.
b
Construct directional and norm estimates d and rb by
m

X

init

1 X
ai aTi 1 {i ∈ Isel } ,
:=
m

db = argmin dT X init d
d∈Sn−1

i=1

rb2 := argmin G(r) :=
r

1
m

m
X

2
bi − r ai , db .

i=1

b
return (b
r, d)
end
Algorithm 3: Initialization procedure with outliers
Let us continue our running example of Gaussian random variables to motivate the proposition.
iid
Example 1 (continued): We claim that for ai ∼ N(0, In ) we have
κ?st :=

inf

u,v∈Sn−1

E[|ha, uiha, vi|] =

2
.
π

Indeed, let Zu = ha, ui and Zv = ha, vi, which are jointly Gaussian with Var(Zu ) = Var(Zv ) = 1,
and let X, Y be independent N(0, 1). Then
inf

u,v∈Sn−1

E[|Zu Zv |] =

E [|(cos(θ)X + sin(θ)Y )(cos(θ)X − sin(θ)Y )|]


= inf E ρX 2 − (1 − ρ)Y 2 .
ρ∈[0,1] |
{z
}
inf

θ∈[0,π/2]

=:f (ρ)

The function f (·) is convex, and it is symmetric around
f (ρ) ≥ f (1/2) = 2/π. ♣

1
2

because f (ρ) = f (1 − ρ). We thus have

In the Gaussian measurement case, whenever pfail < π1 ≈ .318, with high probability we have the
stability f (x) − f (x? ) ≥ λ kx − x? k kx + x? k for a numerical constant λ > 0 as long as m/n is larger
than some numerical constant.

4.2

Initialization

The last ingredient for achieving strong convergence guarantees for the prox-linear algorithm for
phase retrieval is to provide a good initialization x0 ≈ x? . The strategies in the noiseless setting
in Section 3 will fail because of corruptions. With this in mind, we present Algorithm 3, which
provides an initializer in corrupted problems.
Before turning to the analysis, we provide some intuition for the algorithm. We must construct
two estimates: an estimate db of the direction d? = x? / kx? k2 and an estimate rb of the radius, or
magnitude, of the signal r? = kx? k2 . For the former, a variant of the spectral initialization Wang
et al. [39] propose suffices. If we take the Isel ⊂ [m] to be the set of indices Isel corresponding
to the smallest (say) bi , in either model M1 or M2 the
indices i ∈ I out are independent of the
P
measurements ai , so we expect that X init = |Isel |−1 i∈Isel ai aTi = zIn − z 0 d? d? T + ∆, where
z, z 0 are random positive constants—with z 0 large because vectors i with hai , x? i2 small should
16

be nearly orthogonal to d? —and ∆ is an error matrix coming from both randomness in the ai
and the corruptions. As long as the error ∆ is small, the minimum eigenvector of X init should be
approximately d? . Once we have a good initializer db ≈ d? , a natural idea to estimate r? is to pretend
√
that db is the direction of the signal, substitute the variable x = rdb into the objective (2), and
solve for r to get a robust estimate of the signal strength x? . As we show presently, this procedure
succeeds with high probability (and the estimate rb is good even when the data are non-Gaussian).
We first show that our estimate rb of the signal size kx? k2 is robust (See Sec. B.2 for a proof).
Proposition 4. Let Assumption A2 hold and E[aaT ] = In . Let δ ∈ [0, 1] and pfail ∈ [0, 21 ]. There
exist numerical constants 0 < c, C < ∞ such that if db is an estimate of d? for which
δ :=

Cσ 2
b d? ) ≤ 1,
dist(d,
1 − 2pfail

then with probability at least 1−2e−cm(1−2pfail )
defined in Alg. 3, satisfy rb2 ∈ [1 ± δ] kx? k22 .

2 /σ 4

(16)

all minimizers rb2 of G(r) =

1
m

Pm

b2
i=1 |bi −rhai , di |,

Given Proposition 4, finding a good initialization of x? reduces to finding a good estimate
db of the direction d? = x? / kx? k2 . To make this precise, let r? := kx? k and assume that δ =
Cσ 2
b ?
r − r? | ≤ δr? .
1−2pfail dist(d, d ) ≤ 1 as in Proposition 4; assume also the relative error guarantee |b
Using the triangle inequality and Proposition 4, for x0 = rbdb we have
h
i
b d? ) + |r? − rb| ≤ r? (1 + δ) dist(d,
b d? ) + δ ≤ 2δr? =
dist(x0 , x? ) ≤ rb dist(d,

Cσ 2
b d? )r? ,
dist(d,
1 − 2pfail

as claimed. We turn to the directional estimate; to make the analysis cleaner we make the normality
iid

Assumption A4. The measurement vectors ai ∼ N(0, In ).
b we require additional notation for quantiles of Gaussian and χ2 -random
To state our guarantee on d,
variables. Let W ∼ N(0, 1) and define the constant qfail and its associated χ2 -quantile wq2 by
P(W 2 ≤ wq2 ) = qfail :=

1
1 − 2pfail
3 − 2pfail
+
=
< 1.
2(1 − pfail ) 4(1 − pfail )
4(1 − pfail )

(17)

Second, define the constant δq = 1 − E[W 2 | W 2 ≤ wq2 ]. We have the following guarantee.
1 Pm
T
Proposition 5. Let db be the smallest eigenvector of X init := m
i=1 ai ai 1 {i ∈ Isel }, as in Alg. 3.
Let the constant M = 0 if Model M1 holds and M = 1 if Model M2 holds. There are numerical
constants 0 < c, C < ∞ such that for t ≥ 0, with probability at least
!
mδq2
1 − exp(−mt) − exp(−c(1 − 2pfail )m) − exp −c 2
wq
we have
b d? ) ≤ 
dist(d,

C

pn

m +t
 .
p
n
(1 − 2pfail )δq − C m
+ t − M pfail +

For intuition, we provide a few simplifications of Proposition 5 by bounding the quantities wq
and δq defined in Eq. (17). Using the conditional expectation bound in Lemma 3.2 and a more
careful calculation for Gaussian random variables (see Lemma C.1 in the appendices) we have
(
)
 2

wq2
1 2
2
2
2
2
E ai,1 | ai,1 ≤ wq ≤ min
, 1 − wq P(ai,1 ≥ wq ) < 1,
3
2
17

w2

3−2pfail
9
so δq ≥ max{1 − 3q , 21 wq2 P(W 2 ≥ wq2 )}. For pfail ≤ 83 , we have qfail = 4(1−p
≤ 10
and we may
fail )
1
2
take wq ≤ 2.71 and δq > 11 . More generally, a standard Gaussian calculation that Φ−1 (p) ≥
p
1
fail )
| log(1 − p)| as p → 1 shows that wq2 = Φ−1 ( qfail2+1 )2 ≥ log 8(1−p
1−2pfail as pfail → 2 , so δq ≥
8(1−pfail )
1−2pfail
8(1−pfail ) log 1−2pfail as well. Under Model M1, then, as long as the sample is large enough and
b d? ) with probability
pfail < 12 , we can achieve constant accuracy in the directional estimate dist(d,
−cm
of failure e
for c a numerical constant. Under the more adversarial noise model M2, we require

a bit more; more precisely, we must have (1 − 2pfail )δq − pfail > 0 to achieve accurate estimates. A
numerical calculation shows that if pfail < 14 , then this condition holds, so that under Model M2
b d? ) with high probability.
we can achieve constant accuracy in the directional estimate dist(d,

4.3

Summary and success guarantees

With the guarantees of stability, quadratic approximation, and good initialization for suitably
random matrices A ∈ Rm×n , we can provide a theorem showing that the prox-linear approach to the
composite optimization phase retrieval objective succeeds with high probability. Roughly, once m/n
is large than a numerical constant, the prox-linear method with noisy initialization succeeds with
exponentially high probability, even with outliers. Indeed, combining our convergence Theorem 1
with Propositions 3, 4, 5, and Corollary 3.2 of the preceding section, we have the following theorem.
Theorem 3. Let Assumptions A4 hold. There exist numerical constants c > 0 and C < ∞ such
that the following hold for any t ≥ 0. Let the independent outliers Model M1 hold and pfail < π1 or
the adversarial outliers Model M2 hold and pfail < 14 . Let x0 be the initializer returned by Alg. 3,
and assume the iterates xk of Alg. 1 are generated without error. Then
r
n
C
k
?
dist(x0 , x ) ≤
+ t · kx? k2 and dist(xk , x? ) ≤ kx? k2 2−2 for all k ∈ N
2
(1 − 2pfail )
m
2m

with probability at least 1 − 4e−mt − 5e−c(1−2pfail )

.

Thus, we see that the method succeeds with high probability as long as the sample size is large
enough, though there is non-trivial degradation with substantial outliers.

5

Optimization methods

In principle, the sub-problems (5) to minimize fx (y)+ L2 kx − yk22 are convex optimization problems,
but in practice, we require some care to solve the problems at large scale. To that end, in this
section, we describe the three schemes we use to solve the sub-problem (5) (one is simply using the
industrial Mosek solver). We evaluate these more carefully in the experimental section to come.
To fix notation, let φ(·) be the elementwise square operator, that is, for any vector x ∈ Rn
1
we have [φ(x)]j = x2j and φ(x) ∈ Rn+ . Recalling that f (x) = m
kφ(Ax) − bk1 , we perform a few
simplifications to more easily describe the schemes we use to solve the prox-linear sub-problems.
Assuming we begin an iteration at point x0 , defining the diagonal matrix D = 2 diag(Ax0 ) ∈ Rm×m ,
and setting c = b + φ(Ax0 ), we recall inequality (8) and have
f (y) ≤ fx0 (y) +

1
1
(y − x0 )T AT A(y − x0 ) ≤
kDAy − ck1 + m−1 AT A
m
m

18

op

ky − x0 k22 .

Rewriting this with appropriately rescaled diagonal matrices D and vector c, implementing Algorithm 1 becomes equivalent to solving a sequence of optimization problems of the form
minimize kDAx − ck1 +
x

1
kxk22 .
2

(18)

In small scale scenarios, the problem (18) is straightforward to solve via standard interior point
method software; we use Mosek via the Convex.jl toolbox [35]. We do not describe this further.
In larger-scale scenarios, we use more specialized methods, which we now describe.

5.1

Graph splitting methods for the prox-linear sub-problem

When the matrices are large, we use a variant of the Alternating Directions Method of Multipliers
(ADMM) procedure known as the proximal operator graph splitting (POGS) method of Parikh
and Boyd [28, 29], which minimizes objectives of the form f (x) + g(y) subject to a linear constraint
Bx = y. Let the matrix B = DA for shorthand; evidently, problem (18) has precisely this form
and is equivalent to
1
minimize
ky − ck1 + kxk22 subject to Bx = y.
n
m
x∈R ,y∈R
2
The POGS method iterates to solve problem (18) as follows. Introduce dual variables λk ∈ Rn and
νk ∈ Rm associated to x and y, and consider the iterations


1
ρ
2
2
xk+ 1 = argmin
kxk2 + kx − (xk − λk )k2
2
2
2
x
n
o
ρ
yk+ 1 = argmin ky − ck1 + ky − (yk − νk )k22
2
2
y
(19)
#
"



 
−1
In B T xk+ 12 + λk
xk+1
In B T
=
0
0
yk+1
B −Im
yk+ 1 + νk
2

λk+1 = λk + (xk+ 1 − xk+1 ), νk+1 = νk + (yk+ 1 − yk+1 ).
2

2

Each of the steps of the method (19) is trivial except for the matrix inversion, or the “graph
projection” step, which projects the pair (xk+ 1 + λk , yk+ 1 + νk ) to the set {x, y : Bx = y}. The
2
2
first two updates amount to
xk+ 1 =
2

ρ
(xk − λk ) and yk+ 1 = c + sgn(yk − νk − c)
2
1+ρ

[|yk − νk − c| − 1/ρ]+

where denotes elementwise multiplication and each operation is element-wise. The matrix B is
tall, so setting vk = xk+ 1 + λk + B T (yk+ 1 + νk ) ∈ Rn , then the solution of the system
2



In
B

2

   
BT
x
v
= k is xk+1 = (In + B T B)−1 vk and yk+1 = Bxk+1 .
−In y
0n

(20)

In this iteration, it is straightforward to cache the matrix (In +B T B)−1 , or a Cholesky factorization
of the matrix, so that we can repeatedly compute the multiplication (20) in time n2 + nm.
Following Parikh and Boyd [29], we define the primal and dual residuals
"
#


x
−
x
1
k+1
xk − xk+1
k+ 2
pri
dual
rk+1 :=
and rk+1 := ρ
.
yk+1 − yk+ 1
yk − yk+1
2

19

These residuals define a natural stopping criterion [29], where one terminates the iteration (19)
once the residuals satisfy


√
√
rkpri <  n + max{kxk k2 , kyk k2 } and rkdual <  n + max{kλk k2 , kνk k2 }
(21)
2

2

for some  > 0, which must be specified. In our case, the quadratic convergence guarantees of
Theorem 1 suggest a strategy of decreasing  across iterative solutions of the sub-problem (18).
That is, we begin with some  = 0 . We perform iterations of the prox-linear Algorithm 1 using
the POGS iteration (19) (until the residual criterion (21) holds) to solve the inner problem (18).
Periodically, in the outer prox-linear iterations of Alg. 1, we decrease  by some large multiple.
In our experiments with the POGS method for sub-problem solutions, we perform two phases.
In the first phase, we iterate the prox-linear method (Alg. 1) until either k = 25 or kxk − xk+1 k2 ≤
δ/ (1/m)AT A op , where δ = 10−3 , using accuracy parameter  = 10−5 for POGS (usually, the
iterations terminate well-before k = 25). We then decrease this parameter to  = 10−8 , and begin
the iterations again.

5.2

Conjugate gradient methods for sub-problems

In a number of scenarios, it is possible to multiply by the matrix A quickly, though the direct
computation (In +AT D2 A)−1 in expression (20) (recall B = DA) may be difficult. For example, if A
is structured (say, a Fourier or Hadamard transform matrix or or sparse), computing multiplication
by In + AT D2 A quickly is possible. This suggests [33, Part VI] using conjugate gradient methods.
We make this more explicit here to mesh with our experiments to come. Let Hn be an n × n
orthogonal matrix for which computing the multiplication Hn v is efficient (e.g. a Hadamard or
discrete cosine transform matrix). We assume that the measurement matrix A takes the form of
repeated randomized measurements under Hn , that is,

T
A = Hn S1 Hn S2 · · · Hn Sk
(22)
where Sl ∈ Rn×n are (random) diagonal matrices, so that A ∈ Rm×n with m = kn. In this case, the
expensive part of the system (20) is the solution of (I + AT D2 A)x = v. This is a positive definite
system, and it is possible to compute the multiplication (I + AT D2 A)x in time O(kn + kTmult ),
where Tmult is the time to muliply an n-vector by the matrix Hn and HnT . Moreover, when S is a
random sign matrix and Hn is a Hadamard or FFT matrix, the matrix AT D2 A is well-conditioned,
as is known from the analysis of random (subsampled) Hadamard and Fourier transforms (cf. [34]).
Thus, in our experiments with structured matrices we use the conjugate gradient method without
preconditioning [33], iterating until we have relative error (I + AT DA)x − v ≤ cg / kvk where
cg = 10−6 .

6

Experiments

In this section, we describe the experiments we perform to evaluate our methods. We perform a
number of simulations as well as experiments on real images to evaluate our method and compare
with other state-of-the-art (non-convex) methods for phase retrieval problems. To standardize
notation and remind the reader, in each our experiments, we generate data via (variants) of the
following process. We take a measurement matrix A ∈ Rm×n , from one of a few distributions, and
generate a signal x? ∈ Rn either by drawing x? ∼ N(0, In ) or by taking x? ∈ {−1, 1}n uniformly at
random. We receive observations of the form
bi = hai , x? i2 ,
20

and with probability pfail ∈ [0, 21 ], we corrupt the measurements bi .
In our experiments, to more carefully isolate the relative performance of the iterative algorithms,
rather than initialization used, we compare three initializations. The
Gaussianity
P first two rely
Pon
m
−1
? 2
of the measurement matrix A, because in this case we have m−1 m
b
=
m
i=1 i
i=1 hai , x i ≈
x? T E[a1 aT1 ]x? = kx? k22 .
(i) Big: The first, which we call “Big” initialization, is that of Wang et al. [39]. Then we define
o
n
X
I0 := i ∈ [m] : bi / kai k22 ≥ quant5/6 ({bi / kai k22 }) and X init :=
i∈I0

1
T
2 ai ai .
kai k2

1
1 Pm
2 b
We then set the direction db = argmaxkdk2 =1 dT X init d and x0 = ( m
i=1 bi ) d.

(ii) Median: The second initialization, the “median” initialization, is due to Zhang et al. [40].
We set rb2 = quant 1 ({bi })/.455 and define
2

I0 := {i ∈ [m] : |bi | ≤ 9λ0 } and X init :=

X

bi ai aTi .

i∈I0

b
We then set the direction db = argmaxkdk2 =1 dT X init d and x0 = rbd.
(iii) Small: The final initialization is the outlier-aware initialization we describe in Sec. 4.2.
Our convergence results in Section 2 rely on four quantities: the quality of the initialization,
dist(x0 , x? ); the stability parameter of the objective, that is, the λ such that f (x) − f (x? ) ≥
λ kx − x? k2 kx + x? k2 ; the quadratic upper bound guaranteed by our linearized models fx (y), so
that |fx (y) − f (y)| ≤ L2 kx − yk22 ; and the accuracy to which we solve the optimization problems.
The random matrix A governs the middle two quantities—stability λ and closeness L—and we take
2
L= m
|||A|||2op . Thus, in our experiments we directly vary the initialization scheme to generate x0
and the accuracy to which we solve the sub-problems (5), that is,
(
)
m
1 X
L
2
2
xk+1 = argmin
hai , xk i + 2 hai , xk i hai , x − xk i − bi + kx − xk k2 .
m
2
x
i=1

We perform each of our experiments on a server with a 16 core 2.6 GhZ Intel Xeon processor
with 128 GB of RAM using julia as our language, with OpenBLAS as the BLAS library. We restrict
each optimization method to use 4 of the cores (OpenBLAS is multi-threaded).

6.1

Simulations with zero noise and Gaussian matrices

For our first collection of experiments, we evaluate the performance of our method for recovering
random signals x? ∈ Rn using Gaussian measurement matrices A ∈ Rm×n , varying the number of
measurements m over the ten values m ∈ {1.8n, 1.9n, . . . , 2.7n}. (Taking m ≥ 2.7n yielded 100%
exact recovery in all the methods we experiment with.) We assume a noiseless measurement model,
so that for A = [a1 · · · am ]T ∈ Rm×n , we have bi = hai , x? i2 for each i ∈ [m]. We perform experiments with dimensions n ∈ {400, 600, 800, 1000, 1500, 2000, 3000}, and within each experiment we
vary a number of problem parameters, including initialization scheme and the algorithm we use to
solve the sub-problems (5).
To serve as our baseline for comparison, we use Wang, Giannakis, and Eldar’s Truncated Amplitude Flow (TAF) [39], given that it outperforms other non-convex iterative methods for phase
21

retrieval, including
Wirtinger Flow [10], Truncated Wirtinger
Flow [12], and Amplitude Flow [39].
√
1 Pm
Setting ψi = bi , TAF tries to minimize the loss 2m i=1 (ψi − | hai , xi |)2 via a carefully designed
(generalized) gradient method. For convenience, we replicate the method (as described by Wang
et al. [39]) in Alg. 4.
Data: Initializer x0 ∈ Rn , i.i.d. standard Gaussian matrix A ∈ Rm×n , ψi = | hai , x? i |, and
parameters γ = .7 and α = .6
for k = 0, 1, . . . , K − 1 do
1
Set Ik = {i ∈ [m] : | hai , xk i | ≥ 1+γ
ψi }


P
hai ,xk i
α
Update xk+1 = xk − m
ha
,
x
i
−
ψ
i
i
k
i∈Ik
|hai ,xk i| ai
end
return xK

6.1.1

Algorithm 4: Truncated amplitude flow

Low-dimensional experiments with accurate prox-linear steps

We begin by describing our experiments with dimensions n ∈ {400, 600, 800}. For each of these,
we solve the iterative sub-problems to machine precision, using Mosek and the Julia package
Convex.jl [35]. We plot representative results for n = 400 and n = 800 in Figure 1, which
summarize 400 independent experiments, and we generate the true x? by taking x? ∈ {−1, 1}n
uniformly at random. (In separate experiments, we drew x? ∼ N(0, In ), and the results were essentially identical.) In these figures, we use the “big” initialization (i) (the initialization of Wang
et al.), as it yields the best empirical performance. Following Wang et al. [39], we perform 1000
iterations of TAF (Alg. 4), and within each experiment, both TAF and the prox-linear method use
identical initializer and data matrix A. We run the prox-linear method until sequential updates
xk and xk+1 satisfy kxk − xk+1 k2 ≤  = 10−5 , which in every successful experiment we perform
(with these data settings) requires 6 or fewer iterations. We declare an experiment successful if the
output x
b of the algorithm satisfies
dist(b
x, x? ) = min {kb
x − x? k2 , kb
x − x? k2 } ≤ acc kx? k2 where acc = 10−5 .

(23)

In Fig. 1(a), we plot the number of times that one method succeeds (out of the 400 experiments)
while the other does not as a function of the ratio m/n. We see that for these relatively small
dimensional problems, the prox-linear method has somewhat better performance for m/n small
than does truncated amplitude flow. In Fig. 1(b), we plot the fraction of successful runs of the
algorithms, again against m/n for n = 400, and in Fig. 1(c) we plot the fraction of successful runs
for n = 800. We see that even when m/n = 2, the prox-linear method has success rate of around
.6, which meets the state-of-the-art performance of TAF, with mildly better performance for small
problems.
6.1.2

Medium dimensional experients with inaccurate prox-linear steps

We now shift to a description of our experiments in dimensions n ∈ {1000, 1500, 2000, 3000}, again
without noise in the measurements; we perform between 100 and 400 experiments for each dimension
in this regime, and we use the “big” initialization (i) [39]. In this case, we use Parikh and Boyd’s
proximal operator graph splitting (POGS) method [28] for the prox-linear steps, as we describe in
Section 5.1. As described in Sec. 5.1, we perform two phases of the prox-linear method: the first
performing POGS until the residual errors (21) are less than  = 10−5 within the prox-linear steps,
22

1.0

Prox
TAF

0.8

0.8

0.10
0.08
0.06
0.04

0.7

P(success)

0.12

1.0

0.9

P(success)

Fraction of one-sided success

0.14

0.6
0.5
0.4
0.3

0.02

0.6
0.4
0.2

0.2

0.00
1.8

2.0

2.2

m/n

2.4

2.6

2.8

0.1
1.8

2.0

2.2

(a)

m/n

2.4

2.6

2.8

0.0
1.8

2.0

2.2

(b)

m/n

2.4

2.6

2.8

(c)

Figure 1. Zero-noise experiments in dimensions n = 400 and n = 800. (a) Fraction of times (in 400
experiments) that one method succeeded while other failed for n = 400. (b) Fraction of successful
recoveries for n = 400. (c) Fraction of successful recoveries for n = 800.

0.40

0.30
0.25
0.20
0.15
0.10
0.05
0.00
1.8

1.9

2.0

2.1

m/n

2.2

2.3

2.4

2.5

1.0

Prox
TAF

0.35

0.8

0.30
0.25

P(success)

Prox
TAF

Fraction of one-sided success

Fraction of one-sided success

0.35

0.20
0.15
0.10

0.6
0.4
0.2

Prox
TAF

0.05
0.00
1.80

1.85

1.90

1.95

(a)

2.00

m/n

(b)

2.05

2.10

2.15

2.20

0.0
1.80

1.85

1.90

1.95

2.00

m/n

2.05

2.10

2.15

(c)

Figure 2. Zero-noise experiments. (a) Fraction of trials (of 400) in which one method succeeds
while other fails for n = 1000. (b) Fraction of trials (of 100) in which one method succeeds while
other fails for n = 3000. (c) Fraction of successful recoveries for n = 3000.

√
the second to accuracy  = 10−8 . We apply the prox-linear method to the matrix A/ m with
√
data b/m, which is equivalent but numerically more stable because A/ m and In are comparable
(see the recommendations [28]). In these accuracy regimes, the time for solution of the prox-linear
method and that required for 1000 iterations of truncated amplitude flow are comparable; TAF
requires about 1.5 seconds while the two-phase prox-linear method requires around 4 seconds in
dimension n = 1000, and TAF requires about 5 seconds while the the two-phase prox-linear method
requires around 20 seconds in dimension n = 2000, each with m = 2n.
We provide two sets of plots for these results, where again we measure success by the criterion (23), kx ± x? k2 ≤ acc = 10−5 . In the first, Fig. 2, we show performance of prox-linear against
TAF specifically for dimensions n = 1000 and 3000. In Fig. 2(a), we plot the number of trials
in which one method succeeds (out of 400 experiments) while the other does not as a function of
the ratio m/n. In Fig. 2(b) we plot the number of trials in which the prox-linear or TAF method
succeeds, while the other does not, for n = 3000 with x? chosen either N(0, In ) or uniform in {±1}n
(O and + markers, respectively). We ignore ratios m/n ≥ 2.2 as both methods succeed in all of
our trials. Out of 100 trials, there is only one (with m/n = 2) in which TAF succeeds but the
prox-linear method does not. In Fig. 2(c), we plot the fraction of successful runs of the algorithms,
again against m/n for n = 3000. For these larger problems, there is a substantial gap in recovery
23

2.20

0.30
0.25

0.8

P(success)

Fraction of one-sided success

1.0

Prox
TAF

0.35

0.20
0.15

0.6
0.4

0.10
0.2

Prox
TAF

0.05
0.00
1.8

2.0

2.2

m/n

2.4

0.0
1.8

2.6

2.0

2.2

m/n

2.4

2.6

Figure 3: Collected zero-noise experiments in dimensions n ∈ {1000, 1500, 2000}.

probability between prox-linear method and TAF, where with m/n = 2 the prox-linear method
achieves recovery more than 78% (±4), 88% (±6), and 91% (±6) of the time, with 95% confidence
intervals, for n = 1000, 2000, and 3000, respectively. The prox-linear algorithm provides substantial
gains over TAF for these higher dimensional regimes.
In Figure 3, we plot results for all of our experiments with n ∈ {1000, 1500, 2000} simultaneously, where we include both x? ∼ Uni({−1, 1}n ) and x? ∼ N(0, In ), and all experiments with the
initialization. The left plot shows the fraction of the experiments in which one method succeeds
while the other does not. The other plots the probability of success of the methods versus the ratio
m/n. Each line consists of the data of between 100 (n = 2000, 3000) and 400 (n = 1000, 1500)
experiments. The prox-linear method appears to have stronger performance and more frequent
recovery, even with the same data and same initialization.

6.2

Phase retrieval with outlying measurements

One of the claimed advantages of the objective (2) is that it should be robust to outliers. Zhang et al.
[40] develop a method, which they term median-truncated Wirtinger flow, for handling outlying
measurements, which (roughly) sets the index set Ik used for updates in Alg. 4 to be a set of
measurements near the median values of the bi . Their algorithm requires a number of parameters,
which strongly rely on the Gaussianity assumptions of the measurement matrix A, whose entries
are assumed i.i.d. N(0, 1); in contrast, the objective (2) and prox-linear method we investigate are
straightforward to implement without any particular assumptions on A (and require no particular
parameter tuning, relying on an explicit sequence of convex optimizations).
Nonetheless, to make comparisons between the algorithms as fair as possible, we implement
their procedure and perform experiments in dimension n ∈ {100, 200} with i.i.d. standard normal
data matrices A. We perform 100 experiments as follows. Within each experiment, we evaluate
each m ∈ {1.8n, 2n, 3n, 4n, 6n, 8n} and failure probability pfail ∈ {0, .01, .02, . . . , .29, .3}. For fixed
m, n, we draw a data matrix A ∈ Rm×n , then choose x? ∈ Rn either by drawing x? ∼ Uni({−1, 1}n )
or x? ∼ N(0, In ). We then generate bi = hai , x? i2 for i ∈ [m]. For our experiments with n = 100,
we simply set bi = 0, which is more difficult for our initialization strategy, as it corrupts a large
fraction of the vectors ai used to initialize x0 . For our experiments with n = 200, we draw bi from
24

2.8

0.0

0.25

0.5

0.75

1.0

P(success)
1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

1.8
2.0
2.5
3.0
4.0
6.0
8.0

(a)
1.8
2.0
2.5
3.0
4.0
6.0
8.0

(b)

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

1.8
2.0
2.5
3.0
4.0
6.0
8.0

(c)
1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

(d)

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3
pfail

(e)

1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3
pfail

(f)

Figure 4. Probability of success for median-truncated Wirtinger flow (left) and the prox-linear
method (right) for different initializations, where dimension n = 100. Horizontal axis indexes pfail
while vertical axis indexes the measurement ratio m/n. Each pixel represents the fraction of successful
recoveries (defined as dist(b
x, x? )/ kx? k2 ≤ 10−5 ) over 100 experiments. (a) & (b): “Big” initialization,
with rb estimated by Alg. 3. (c) & (d): “Median” initialization. (e) & (f): “Small” initialization.

a Cauchy distribution. Each problem setting (bi Cauchy vs. zeroing and x? discrete or normal)
yields qualitatively similar results.
Figures 4, 5, and 6 summarize our results. In Figures 4 and 5, we display success rates of the
median truncated Wirtinger flow (the left column in each figure) and our composite optimizationbased procedure (right column) for a number of initializations; each plot represents results of
100 experiments. Within each plot, a white square indicates that 100% of trials were successful,
meaning the signal is recovered to accuracy dist(b
x, x? ) ≤ 10−5 kx? k2 , while black squares indicate
0% success rates. Within each row of the figure, we present results for the two methods using
the same initialization scheme. Figure 4 gives results with n = 100, using precise (Mosek-based)
solutions of the prox-linear updates, while Figure 5 gives results with n = 200 using the POGSbased updates (recall step (19)), with identical parameters as in the previous section. It is clear
from the figures that the composite objective yields better recovery.
In Figure 6, we present a different view into the behavior of the prox-linear and MTWF methods. In the left two plots, we show the recovery probability (over 100 trials) for our composite
optimization method (Fig. 6(a)) and MTWF (Fig. 6(c)). The success probability for the composite
method is higher. In Fig. 6(b), we plot the average number of iterations (along with standard
error bars) the prox-linear method performs when the dimension n = 100 and we use accurate
sub-problem solves. We give iteration counts only for those trials that result in successful recovery; the iteration counts on unsuccessful trials are larger (indeed, if the method is not converging
rapidly, this serves as a proxy for failure). In the high measurement regime, m/n ≥ 2.5 or so, we
25

1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

1.8
2.0
2.5
3.0
4.0
6.0
8.0

(a)
1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

(b)

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

(c)

1.8
2.0
2.5
3.0
4.0
6.0
8.0

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

(d)

Figure 5. Probability of success for median-truncated Wirtinger flow (left) and the prox-linear
method (right), using the approximate proximal graph splitting method for sub-problem solves.
Dimension n = 200. (a) & (b): “Big” initialization, with rb estimated by Alg. 3. (c) & (d): “Small”
initialization.

see that if pfail = 0 no more than 7 iterations are required: this is the quadratic convergence of
the method. (Indeed, for m/n = 8, for pfail ≤ .15 each execution of the prox-linear method uses
precisely 5 iterations, never more, and never fewer.) In Fig. 6(d), we show the number of matrix
multiplications by the inverse matrix (In +AT DA)−1 the method uses (recall the update (19) in the
proximal graph operator splitting method). We see that for well-conditioned problems and those
with little noise, the methods require relatively few matrix multiplications, while for more outliers
and when m/n shrinks, there is a non-trivial increase.

6.3

Recovery of real images

Our final collection of experiments investigates recovery of real-world images using more specialized
√
measurement matrices. In this case, we let Hn ∈ {−1, 1}n×n / n denote a normalized Hadamard
matrix, where the multiplication Hn v requires time n log n, and Hn satisfies Hn = HnT and Hn2 = In .
For some k ∈ N, we then take k i.i.d. diagonal sign matrices S1 , . . . , Sk ∈ diag({−1, 1}n ), uniformly
at random, and define A = [Hn S1 Hn S2 · · · Hn Sk ]T ∈ Rkn×n , as in expression (22). We note
that this matrix explicitly does not satisfy the stability conditions that we use for our theoretical
guarantees. Indeed, letting e1 and e2 be the first standard basis vectors, we have Hn S(e1 + e2 ) ⊥
Hn S(e1 − e2 ) no matter the sign matrix S; there are similarly pathological vectors for FFT, discrete
cosine, and other structured matrices. Nonetheless, we perform experiments with this structured
A matrix as follows.1
Given an image X represented as a matrix, we define x? = Vec(X), the vectorized representation
of X. (In the case of colored images, where X ∈ Rn1 ×n2 ×3 because of the 3 RGB channels, we
vectorize the channels as well.) We then set b = φ(Ax), where φ(·) denotes elementwise squaring,
and corrupt a fraction pfail ∈ [0, .2] of the measurements bi by zeroing them. We then follow our
standard experimental protocol, initializing x0 by the “small” initialization scheme (iii), with the
slight twist that now we use the proximal operator graph splitting method (POGS) with conjugate
gradient methods (Sec. 5.2) to solve the graph projection step (20).
We first give results on a collection of 500 images of handwritten digits (using k = 3), available
1

We do not compare to other methods designed for outliers because we could not set the constants the methods
require in such a way as to yield good performance.

26

1.0

14

0.8

12

0.6

10

0.4

8

0.2

6

0.0
0.00

0.05

0.10

0.15

0.20

0.25

0.30

4
0.00

m/n = 1.8
m/n = 2.0
m/n = 2.5
m/n = 3.0
m/n = 4.0
m/n = 6.0
m/n = 8.0

0.05

0.10

(a)

0.15

0.20

0.25

0.30

(b)

1.0

m/n = 1.8
m/n = 2.0
m/n = 2.5
m/n = 3.0
m/n = 4.0
m/n = 6.0
m/n = 8.0

8000

0.8

7000
6000

0.6

5000

0.4

4000
3000

0.2

2000

0.0
0.00

0.05

0.10

0.15

0.20

0.25

0.30

(c)

1000
0.00

0.05

0.10

0.15

0.20

0.25

(d)

Figure 6.
Comparison of median-truncated Wirtinger flow and composite minimization with
“Median” initialization. Horizontal axis of each plot indexes pfail ∈ [0, .3]. (a) Proportion of successful
solves for prox-linear method. (b) Number of iterations of prox-linear step (5) over successful solves
(with error bars). (c) Proportion of successful solves for median-truncated Wirtinger flow method.
(d) Number of matrix multiplications of the form x 7→ (I + AT D2 A)−1 Ax for the prox-linear method
with proximal graph solves.

27

0.30

Figure 7. Example recoveries of digits. Left digit is true digit, middle is initialization, right is
recovered image.
18000

1.0

16000
14000

Matrix multiply count

Success probability

0.9
0.8

12000

0.7

10000

0.6

8000

0.5

6000

0.4

4000

0.3
0.00

0.05

0.10

pfail

0.15

2000
0.20

Figure 8. Performance of composite optimization scheme using conjugate gradient method and
proximal graph splitting to solve sub-problems (Sections 5.1–5.2).

on the website for the book [22]. We provide example results of the execution of our procedure in
Fig. 7, which shows that while there is signal in the initialization, there is substantial work that
the prox-linear method must perform to recover the images. For each of the 500 images, we vary
pfail ∈ {0, .025, .05, . . . , .175, .2}, then execute the prox-linear method. We plot summary results
in Fig. 8, which in the blue curve with square markers gives the probability of successful recovery
of the digit (left axis) versus pfail (horizontal axis). The right axis indexes the number of matrix
multiplications the method executes until completion (black line with circular marks). We see that
in spite of the demonstrated failure of the matrix A to satisfy our stability assumptions, we have
frequent recovery of the images to accuracy 10−3 or better.
Finally, we perform experiments with eight real color images with sizes up to 1024 × 1024
(yielding n = 222 -dimensional problems), where we use k = 3 random Hadamard sensing matrices.
The prox-linear method successfuly recovers each of the 8 images to relative accuracy at least 10−4 ,
and performs an average of 15100 matrix-vector multiplications (i.e. fast Hadamard transforms)
over the eight experiments, with a standard deviation of 2600 multiplications. To give a sense of
the importance of different parts of our procedure, and the relative accuracies to which we solve
sub-problems (5), we display one example in Fig. 9. In this example, we perform phase retrieval on
an image of RNA nanoparticles in cancer cells [30]. In Fig. 9(a), we show the result of initialization,
which displays non-trivial structure, though is clearly noisy. In Fig. 9(b), we show the result of 10
steps of the prox-linear method, solving each step using POGS until the residual errors (21) are
below  = 10−3 . There are clear artifacts in this image. We then perform one refinement step with
higher accuracy ( = 10−7 ), which results in Fig. 9(c). This is indistinguishable, at least to our
eyes, from the original image (Fig. 9(d)).
28

(a)

(b)

(c)

(d)

Figure 9. Reconstruction of RNA image (n = 222 ) from m = 3n measurements. (a) Initialization
x0 of prox-linear method. (b) Result of 10 inaccurate ( = 10−3 ) solutions of prox-linear step (5).
(c) Result of one additional accurate ( = 10−7 ) solution of prox-linear step (5). (d) Original image.

Acknowledgments We thank Dima Drusvyatskiy and Courtney Paquette for a number of inspirational and motivating conversations. JCD and FR were partially supported by NSF-CAREER
Award 1553086 and the Toyota Research Institute. FR was additionally supported by a Stanford
Graduate Fellowship.
29

A

Proofs for noiseless phase retrieval

In this appendix, we collect the proofs of the propositions and results in Section 3.
Because we will use it multiple times in what follows, we state a standard eigenvector perturbation result—a variant of the Davis-Kahane sin-Θ theorem—here.
Lemma A.1 (Stewart and Sun [32], Theorem 3.6). For vectors u, v ∈ Sn−1 , define the angle
θ(u, v) = cos−1 hu, vi. Let X ∈ Rn×n be symmetric, ∆ a symmetric perturbation, and Z = X + ∆,
and define gap(X) = λ1 (X) − λ2 (X) to be the eigengap of X. Let v1 and u1 be the first eigenvectors
of X and Z, respectively. Then
q
1 − hu1 , v1 i2 = | sin θ(u1 , v1 )| ≤ h

A.1

|||∆|||op
gap(A) − |||∆|||op

i .
+

Proof of Proposition 1

Our proof essentially follows Mendelson’s “small-ball” techniques for concentration [26] along with
control over a particular VC-dimension condition. If we define the random function
hi (u, v) := κ?st 1 {| hai , ui hai , vi | ≥ κ?st } ,
then we certainly have

m

m

i=1

i=1

1 X
1 X
| hai , ui hai , vi | ≥
hi (u, v).
m
m
The class of functions F := {a 7→ f (a) = aT uv T a | u, v ∈ Rn } has VC-dimension at most 2(n + 2),
and the associated thresholds 1 {| ha, ui ha, vi | ≥ κ?st } have VC-dimension at most 2(n + 2) (see
Lemma C.2 for a proof of this essentially standard fact). Applying standard VC concentration
inequalities (e.g. [36, Ch. 2.6] or [3, Thm. 6]), we immediately obtain that
!
r


m
1 1 X
n
mt2
P
sup
hi (u, v) − E[hi (u, v)] ≥ c
+ t ≤ 2 exp −
?
m
2
u,v∈Sn−1 κst m
i=1

for all t ≥ 0, where c < ∞ is a numerical constant. Substitute t2 7→ 2t to achieve the result.

A.2

Proof of Proposition 2: deterministic part

Our proof of Proposition 2 eventually reduces to applying eigenvector perturbation laws to the
random matrices X init . To motivate our approach, note that
X init :=

1 X
ai aTi = In − φ() · d? d? T + ∆,
|Isel |

(24)

i∈Isel

where ∆ is a random error term that we will show has small norm. From this quantity—at
least heuristically—we can apply standard eigenvector perturbation arguments to derive that the
directional estimate db = argmind∈Sn−1 dT X init d satisfies db ≈ d? , which will hold so long as |||∆|||op is
small because there is substantial separation in the eigenvalues of In and In − φ()d? d? T . Showing
that rb ≈ kx? k2 is essentially trivial, so that demonstrating equality (24) is enough for the result.

30

With this goal in mind, we define two index sets that (with high probability) surround Isel . Let




1−
1+
? 2
? 2
I− := i ∈ [m] | hai , d i <
and I+ := i ∈ [m] | hai , d i ≤
.
2
2
We now define five events, showing that conditional
Pm onTthese five events, the result of the proposition
holds. These events roughly guarantee that i=1 ai ai is well-behaved, that I+ \ I− is small, and
thus that most of the vectors ai for indices i ∈ Isel are indeed close enough to uniform on the
subspace perpendicular to d? that we have a good directional estimate. Now, let q ∈ [1, ∞] and let
1/p + 1/q = 1, so that p is its conjugate. We define


1
T
E1 :=
A A op ∈ [1 − , 1 + ] , E2 := {|I+ | ≤ |I− | + 2κm} ,
m






X
1
1
1
E3 := |I− | ≥ mp0 (d? ) , E4 :=
ai aTi
≤ 4qσ 2 (κ) p
(25)
 m

2
op
i∈I+ \I−






X
1
T
? ?T
E5 :=
≤ |||∆()|||op +  .
ai ai − In − φ()d d
 |I− |

op
i∈I−

With these five events defined, we prove that the result of the proposition holds whenever each
of the five E1 , . . . , E5 occur. Decompose the matrix X init into the following four parts:
X init = In − φ()d? d? T
|
{z
}
:=Z0





 X

X
1
1
1
1
T
? ? T 
T

ai aTi .
ai ai − (In − φ()d (d ) ) +
−
+
ai ai −
|I− |
|I− |
|I− | |Isel |
i∈Isel
i∈I−
i∈Isel \I−
{z
}
{z
} |
|
{z
} |
X

:=Z3

:=Z2

:=Z1

We bound the operator norms of each of Z1 , Z2 , Z3 in turn. On the event E5 , we immediately have
|||Z1 |||op ≤ |||∆()|||op + .

(26)

We turn to the error matrix Z2 . On the event E1 , we evidently have rb ∈ kx? k2 (1 ± ) by definition
1
of rb2 = m
kAx? k22 , so that
I− ⊂ Isel ⊂ I+ .
Using that the summands ai aTi are all positive semidefinite, we thus obtain the upper bound
|||Z2 |||op =

1
|I− |

X
i∈Isel \I−

ai aTi

≤
op

1
|I− |

X
i∈I+ \I−

(i)

ai aTi

≤
op

1
1
· 4qσ 2 (κ) p ,
?
p0 (d )

(27)

1
P
where in inequality (i) we use that 2|I− | ≥ p0 (d? )m on E3 and ||| i∈I+ \I− ai aTi |||op ≤ 4qmσ 2 (κ) p
on E4 . Lastly, we provide an upper bound on |||Z3 |||op . Again using that I− ⊂ Isel ⊂ I+ on event
E1 , we have
1
1
1
1
|I+ | − |I− |
2κ
−
≤
−
=
≤
,
|I− | |Isel |
|I− | |I+ |
|I− ||I+ |
p0 (d? )2 m

31

where in the last inequality, we use that |I+ | − |I− | ≤ 2κm on E2 and that |I+ | ≥ |I− | ≥
p0 (d? )m/2 on E3 . Thus, by the definition of E1 , we have
|||Z3 |||op ≤

2κ
p0 (d? )2

m

1 X
ai aTi
m
i=1

≤
op

2(1 + )κ
.
p0 (d? )2

(28)

Combining inequalities (26), (27) and (28) on the error matrices Zi , the triangle inequality gives
"
#
1
4qσ 2 (κ) p
2(1 + )κ
.
|||Z1 + Z2 + Z3 |||op ≤ |||∆()|||op +  +
+
p0 (d? )
p0 (d? )2
|
{z
}
:=ν()

This implies equality (24) with error bound |||∆|||op ≤ |||∆()|||op + ν(). Recall the definition of
Z0 = In − φ()d? d? T , which has smallest eigenvector d? and eigengap φ(). Lastly, we simplify ν()
by a specific choice of p and q in the definition (25) of E4 . Without loss of generality, we assume
1
κ < 1 (recall Assumption A3 on κ), and define p = 1 + log1 1 and q = 1 + log κ
. Using that for
κ

any z < 0 we have

exp
1

z
1 − 1/z





1
= exp z +
1 − 1/z



z



= e exp

z
z−1



≤ ez+1 ,

1

1
)(κ) in the error term ν().
we have (κ) p ≤ eκ, allowing us to bound q(κ) p ≤ e(1 + log κ
We now apply the eigenvector perturbation inequality of Lemma A.1. Using that ku − vk22 ≥
2 − 2| hu, vi | ≥ 2 − 2 hu, vi2 for u, v ∈ Sn−1 , a minor rearrangement of Lemma A.1 applied to
X init = Z0 + ∆ for ∆ = Z1 + Z2 + Z3 yields
√
n
o
2(|||∆()|||op + ν())
b d? ) = min kdb − d? k2 , kdb + d? k2 ≤ h
i .
dist(d,
φ() − (|||∆()|||op + ν())
+

Finally, using that x0 = rbdb and defining r? = kx? k, we have by the triangle inequality that
(i)

b d? ) + |r? − rb| ≤
dist(x0 , x? ) ≤ rb dist(d,

h√

i
b d? ) +  r? ,
1 +  dist(d,

where inequality (i) uses E1 , which we recall implies (1 − ) kx? k22 ≤ rb2 ≤ (1 + ) kx? k22 , where
 ∈ [0, 1]. This is the claimed consequence (14) in Proposition 2.

A.3

Proof of Proposition 2: high probability events

It remains to demonstrate that each of the events E1 , . . . , E5 (recall definition (25)) holds with high
probability, to which we dedicate the remainder of this argument in the next series of lemmas, each
of which argues that one of the five events occurs with high probability. Before the statement of
each lemma, we recall the corresponding event whose high probability we wish to demonstrate.
1
|||AT A|||op ∈ [1 ± ]}.
We begin with E1 := { m
Lemma A.2. There exists a numerical constant c > 0 such that P(E1 ) ≥ 1 − exp(−cm2 /σ 4 ) for
m large enough that m/n ≥ σ 4 /(c2 ).

32

2

Proof Set t = c σ 4 in Lemma 3.1, noting that we must have σ 2 ≥ 1 because E[aaT ] = In and
E[aaT ] op ≤ σ 2 (recall the final part of Lemma 3.1). Moreover,  ∈ [0, 1] by assumption. Then
taking c small enough, once we have

n
m

2

≤ c σ 4 we obtain the result.

The event E2 := {|I− | ≥ |I+ | − 2κm} likewise holds with high probability.
Lemma A.3. We have P(E2 ) ≥ 1 − exp(−22 κ2 m).
Proof

We always have that

I− ⊂ I+ and I+ \ I− =


1
1
? 2
i ∈ [m] : (1 − ) ≤ hai , d i ≤ (1 + ) .
2
2

Therefore, the difference in cardinalities of |I− | and |I+ is
m
o
1
1 X n
(|I+ | − |I− |) =
1 1 −  ≤ 2 hai , d? i2 ≤ 1 +  .
m
m
i=1

The right hand side is an average of i.i.d. Bernoulli random variables with means bounded by
κ by Assumption A3(i). Hoeffding’s inequality gives the result that P(|I+ | > |I− | + 2κm) ≤
2 2
e−2κ  m .
Lemma A.4. There exists a numerical constant c > 0 such that the event E3 := {|I− | ≥ 21 mp0 (d? )}
satisfies P(E3 ) ≥ 1 − exp(− 21 mp0 (d? )2 ).
Proof
As in Lemma A.3, this result is immediate from Hoeffding’s inequality. We have |I− | =
Pm
? 2
? 2
?
1{ha
i , d i ≤ (1 − )/2}, an i.i.d sum of Bernoullis with P(hai , d i ≤ (1 − )/2) ≥ p0 (d ) by
i=1
Assumption A3(i). Hoeffding’s inequality gives the result.
Showing that events E4 and E5 in Eq. (25) each happen with high probability requires a little
more work. We begin with E4 , which is defined in terms of a conjugate pair p, q ≥ 1 with 1/p+1/q =
1
P
1 by E4 := {||| i∈I+ \I− ai aTi |||op ≤ 4qσ 2 (κ) p }.
− p2

Lemma A.5. There exists a numerical constant c > 0 such that m/n > c−1 (κ)
2
1 − exp(−cm(κ) p ).

implies P(E4 ) ≥
2

Note that it is no loss of generality to assume that κ ≤ 1 by Assumption A3, so (κ) p ≥ (κ)2 .
2
Proof The {ai }m
i=1 are independent σ -sub-Gaussian random vectors by Assumption A2, and
for any random variable Bi with |Bi | ≤ 1, which may depend on ai , it is clear that the collection
{Bi ai }m
i=1 are mutually independent and still satisfy Definition 3.1. To that end, define the Bernoulli
1+
variables B(a) = 1{ha, d? i2 ∈ [ 1−
2 , 2 ]} (letting Bi = B(ai ) = 1 {i ∈ I+ \ I− } for shorthand).
Then Lemma 3.1 implies for a numerical constant C < ∞ that
r
!
m
1 X
n
n
ai aTi Bi − E[aaT B(a)]
≥ Cσ 2 max
+ t, + t
≤ e−mt .
(29)
P
m
m
m
op
i=1

Now, note by Hölder’s inequality that




1− 1+
E[hv, ai B(a)] ≤ E[hv, ai ] P ha, d i ∈
,
2
2
2

2q

1
q

1

? 2
1

1

≤ (σ 2q Γ(q + 1)e) q (κ) p ≤ qe1/q σ 2 (κ) p
33

 1

p

where we have applied Assumption A3(i) and Lemma 3.1 to bound E[ha, d? i2q ]. Using the triangle
2
1
(κ) p into inequality (29), we find that for any q ∈ (1, ∞) and
inequality and substituting t = 4C
1/p + 1/q = 1, if

n
m

≤

2

1
p
4C (κ)

we have

m

1 X
ai aTi Bi − E[aaT B(a)]
m

P

i=1

2

≥ qσ (κ)

1
p

!

op

2

m(κ) p
≤ exp −
4C

!
.

Applying the triangle inequality and that 1 + e1/q < 4 gives the result.
ThePfinal high probability guarantee is the most complex, and applies to the event E5 :=
1
T
? ?T
{||| |I−
i∈I− ai ai − (In − φ()d d )|||op ≤ |||∆()|||op + }.
|
Lemma A.6. Let E3 = {|I− | ≥
c > 0 such that

1
?
2 mp0 (d )}

c
?

as in Eq. (25). There exists a numerical constant

σ 4 log2 p0 (d? )
m
≥
n
p0 (d? )2
2

0 (d )
implies P(E5 | E3 ) ≥ 1 − exp(−c σ4mp
).
log2 p (d? )
0

Proof
of aaT :

For notational simplicity, define the following shorthand for the conditional expectation


1−
T
? 2
?
= In − φ()d? d? T + ∆(),
Ed := E aa | ha, d i ≤
2

where the equality uses Assumption A3(ii). The main idea of the proof is to show the following
e
crucial fact: define the new sub-Gaussian parameter τ 2 = σ 2 log p0 (d
? ) ≥ 1. Then there exists a
numerical constant 1 ≤ C < ∞ such that for all t ≥ 0,


r

X
n
n
1
≥ Cτ 2 max
ai aTi − Ed?
+ t,
+ t | I−  ≤ exp(−|I− |t). (30)
P
|I− |
|I
|
|I
−
− |
op
i∈I−

Suppose that the bound (30) holds. On the event E3 , we have that |I− | ≥ 21 mp0 (d? ), and so
choosing
2
m
2C 2 τ 4
t=
<
1,
and
letting
≥
4C 2 τ 4
n
p0 (d? )2
yields that

n
|I− |

+ t < 1 on E3 and



P

1
|I− |

X
i∈I−

ai aTi

− Ed?
op



mp0 (d? )2

≥  | I− ≤ exp −
.
8C 2 τ 4

By the definition of Ed? and the triangle inequality we have


X
1
≥ |||∆()|||op +  | I− 
P
ai aTi − (In − φ () d? d? T )
|I− |
op
i∈I−

1 X
≤ P
ai aTi − Ed?
|I− |
i∈I−

34


≥  | I−  .
op

The lemma follows from the fact that E3 is measurable with respect to the indices I− .
Now, we show the key inequality (30). The main idea is to show that, conditioning on the set
I− , the distribution {ai }i∈I− is still conditionally independent and sub-Gaussian. To do so, we
introduce a bit of (more or less standard) notation. For a random variable X, let L(X) denote
the law of distribution of X. Using the independence of the vectors ai , we have the fact that for
any fixed subset I ⊂ [m], the collection {ai }i∈I is independent of {ai }i6∈I . Therefore, using the
definition that I− = {i ∈ [m] : hai , d? i2 ≤ 1−
2 }, we have the key identity




1 X
1 X
1−
dist
L
.
ai aTi | I− = I = L
ai aTi | max hai , d? i2 ≤
i∈I
|I− |
|I|
2
i∈I−

i∈I

This implies that, conditioning on I− = I, the vectors {ai }i∈I are still conditionally independent,
and their conditional distribution is identical to the law L(a | ha, d? i2 ≤ 1−
2 ). The claimed
inequality (30) will thus follow by the matrix concentration inequality in Lemma 3.1, so long as we
can demonstrate appropriate sub-Gaussianity of the conditional law L(a | ha, d? i2 ≤ 1−
2 ).
2 -sub-Gaussian, let J denote all
Indeed, let us temporarily assume that a | ha, d? i2 ≤ 1−
is
τ
2
subsets I ⊂ [m] such that |I| ≥ 12 mp0 (d? ), and define the shorthand Ed? = E[aaT | ha, d? i2 ≤ 1−
2 ].
Then by summing over J , we have on the event E3 that for a numerical constant C < ∞,


r

X
n
1
n
ai aTi − Ed?
≥ Cτ 2 max
+ t,
+ t | E3 
P
|I− |
|I
|
|I
−
− |
op
i∈I−
!
r

X
1 X
n
n
2
T
≥ Cτ max
=
P
ai ai − Ed?
+ t,
+ t | I− = I P(I = I− | E3 )
|I|
|I|
|I|
op
I∈J
i∈I
!
r

X
1 X
n
n
1
−

=
P
≥ Cτ 2 max
+ t,
+ t | max hai , d? i2 ≤
ai aTi − Ed?
P(I = I− | E3 )
i∈I
|I|
|I|
|I|
2
op
I∈J

(i)

≤

X

i∈I

1

? )t

e−|I|t · P(I = I− | E3 ) ≤ e− 2 mp0 (d

,

I∈J

where inequality (i) is an application of Lemma 3.1. This is evidently inequality (30) with appropriate choice of τ 2 .
e
2
2
We thus show that L(a | ha, d? i2 ≤ 1−
2 ) is subgaussian with parameter τ = σ log p0 (d? ) by
bounding the conditional moment generating function. Let λ ∈ [1, ∞] and λ0 be conjugate, so that
1/λ + 1/λ0 = 1. Then by Hölder’s inequality, for any v ∈ Sn−1 we have
 n
oi
h

"
!
#
ha,vi2
2
1−
2
1
ha,
vi
≤
E
exp
2
2
ha, vi
1−
λσ


E exp
| ha, vi2 ≤
=
2
1−
λσ 2
2
P ha, vi ≤
2

1
2
1
2
1− λ0
λ
)])
(E[exp( ha,vi
P
ha,
vi
≤
2
σ2


≤
P ha, vi2 ≤ 1−
2
"
!# 1 
 1

1
λ
λ
ha, vi2
1 −  −λ
e
2
= E exp
P
ha,
vi
≤
≤
,
σ2
2
p0 (d? )


e
where the final inequality uses the σ 2 -sub-Gaussianity of a. Set λ = log p0 (d
? ) to see that condi-

tional on ha, d? i2 ≤

1−
2 ,

e
the vector a is σ 2 log p0 (d
? ) -sub-Gaussian, as desired.

35

B

Proofs for phase retrieval with outliers

In this section, we collect the proofs of the various results in Section 4. Before providing the proofs,
we state one inequality that we use frequently that will be quite useful. Let Wi ∈ {0, 1} satisfy
Wi = 1 if i ∈ I out and Wi = 0 otherwise, so thatPW indexes the outlying measurements. Then
because Wi are independent of the ai vectors and i Wi = pfail m, we have by Lemma 3.1 that for
a numerical constant C < ∞, for all t ≥ 0 we have
!
r
m
n
1 X
n
P
Wi ai aTi − pfail E[ai aTi ]
+ t, + t
≤ e−mt .
(31)
≥ Cσ 2 max
m
m
m
op
i=1

B.1

Proof of Proposition 3

Recalling the set I out of outlying indices, we evidently have
f (x) − f (x? ) =

m

1 X
1 X 
| hai , xi2 − hai , x? i2 | +
| hai , xi2 − ξi | − | hai , xi2 − hai , x? i2 | − f (x? )
m
m
out
i=1

=

(Ax)2

i∈I

(Ax? )2

−
m

1


1 X 
+
| hai , xi2 − ξi | − | hai , xi2 − hai , x? i2 | − | hai , x? i2 − ξi |
m
out
i∈I

(Ax)2 − (Ax? )2
≥
m

1

−

2 X
| hai , xi2 − hai , x? i2 |.
m
out

(32)

i∈I

Now, we note the trivial fact that if we define Wi = 1 for i ∈ I out and Wi = 0 for i ∈ I in , then we
have
X

| hai , xi2 − hai , x? i2 | =

i∈I out

m
X

Wi |(x − x? )T ai aTi (x + x? )| ≤ kdiag(W )A(x − x? )k2 kdiag(W )A(x + x? )k2

i=1

≤ AT diag(W )A op kx − x? k2 kx + x? k2 .
1 Pm
T
Inequality (31) shows that the matrix m
i=1 Wi ai ai is well concentrated.
p
1
T
Now, let E denote the event that m |||A diag(W )A|||op ≤ pfail + Cσ 2 n/m + t, where t is chosen
so that n/m + t ≤ 1. Returning to inequality (32), we have
p
1
f (x) − f (x? ) ≥
(Ax)2 − (Ax? )2 1 − 2(pfail + Cσ 2 n/m + t) kx − x? k2 kx + x? k2
m
for all x ∈ Rn with probability at least 1 − e−mt by inequality (31). We finish with the following
lemma, which is a minor sharpening of Theorem 2.4 of Eldar and Mendelson [17] so that we have
sharp concentration in all dimensions n. We provide a proof for completeness in Section C.2.
Lemma B.1. Let ai be independent σ 2 -sub-Gaussian vectors, and define κst (u, v) := E[| ha, ui ha, vi |]
for u, v ∈ Rn . Then there exist a numerical constants c > 0 and C < ∞ such that
r
m
1 X
2 3 n
| hai , ui hai , vi | ≥ κst (u, v) − Cσ
− σ 2 t for all u, v ∈ Sn−1
m
m
i=1

2

with probability at least 1 − e−cmt − e−cm when m/n ≥ C.
36

Noting that | hai , xi2 − hai , x? i2 | = | hai , x − x? i hai , x + x? i | and substituting the result of the
lemma into the preceding display, we have
r


?
?
2 3 n
2
f (x) − f (x ) ≥ κst − 2pfail − Cσ
− Cσ t kx − x? k2 kx + x? k2
m
2

uniformly in x with probability at least 1 − 2e−cmt − e−cm .

B.2

Proof of Proposition 4

We first state a lemma providing a deterministic bound on the errors of the minimizing radius.
Lemma B.2. Let

6 kx? k22 AT A op
b d? ).
dist(d,
δ=P
P
? i2 −
? i2
ha
,
x
ha
,
x
in
out
i
i
i∈I
i∈I

(33)

If δ ≤ 1, then all minimizers of G(·) belong to the set [1 ± δ] kx? k22 .
Temporarily assuming the conclusions of the lemma, let us show that the random quantities in
the bound (33) are small with high probability. We apply the matrix concentration inequality (31)
n
to see that for a numerical constant C < ∞ and all t ∈ [0, 1 − m
], we have
r
r




1 X
n
1 X
n
2
2
2
2
hai , vi ≥ (1 − pfail ) − Cσ
+ t and
hai , vi ≤ pfail + Cσ
+t
m
m
m
m
out
in
i∈I

i∈I

pn
1
with probability at least 1 − e−mt for all vectors v ∈ Sn−1 , and m
|||AT A|||op ≤ σ 2 (1 + C m
+ t)
n
with the same probability. That is, for t ∈ [0, 1 − m
] with probability at least 1 − 2e−mt we have
that δ in expression (33) satisfies
δ≤
If we assume that
c, we find that

n
m

6σ 2
b d? ).
pn
dist(d,
1 − 2pfail − Cσ 2 m
+t

≤ c(1−2pfail )2 /σ 4 and replace t with c(1−2pfail )2 /σ 4 for small enough constant
Cσ 2
b d? )
dist(d,
1 − 2pfail

δ≤

2

4

with probability at least 1 − 2e−cm(1−2pfail ) /σ , where C is a numerical constant. This is our desired
result.
Proof We define a few pieces of notation for shorthand. Let
1 X
1
1 X
σ
b2 :=
AT A op and σ
bI2 out :=
hai , x? i2 −
hai , x? i2 ,
m
m
m
out
in
i∈I

i∈I

and define the functions gb(δ) = G((1 + δ) kx? k22 ), equivalently
m

1 X
b + δ)|,
gb(δ) :=
|bi − kx? k22 hai , di(1
m
i=1

and a slightly more accurate counterpart
m
1 X
1 X
1 X
g(δ) :=
bi − (1 + δ) hai , x? i2 =
hai , x? i2 − (1 + δ) hai , x? i2 +
bi − (1 + δ) hai , x? i2 .
m
m
m
out
in
i=1

i6∈I

i∈I

37

Note that if δ minimizes gb(δ), then (1 + δ) kx? k22 minimizes G(r). By inspection we find that the
subgradients of g with respect to δ are
∂δ g(δ) =

1 X
1 X
sgn(δ) hai , x? i2 −
sgn((1 + δ) hai , x? i2 − bi ) hai , x? i2 ,
m
m
out
in
i∈I

i∈I

where sgn(t) = 1 if t > 0, \1 if t < 0, and sgn(0) = [\1, 1]. Evidently, for δ > 0 we have g 0 (δ) ≥ σ
bI2 out
and g 0 (−δ) ≤ −b
σI2 out , so that
g(δ) ≥ σ
bI2 out |δ| + g(0).
(34)
Now, we consider the gaps between gb and g: for δ ∈ [\1, 1], we have the gap
m
(1 + δ) kx? k22 X
b 2|
|hai , d? i2 − hai , di
|g(δ) − gb(δ)| ≤
m
i=1

≤

(1 +

δ) kx? k22

AT A

m

op

2
? b
b 2 kd? + dk
b 2 ≤ 4 kx? k2 σ
kd? − dk
2 b dist(d , d),

where we have used the triangle inequality and Cauchy-Schwarz. Thus we obtain
b d? ),
gb(δ) − gb(0) ≥ g(δ) − g(0) + gb(δ) − g(δ) + g(0) − gb(0) ≥ σ
bI2 out |δ| − 6 kx? k22 σ
b2 dist(d,
where we have applied inequality (34). Rearranging, we have that if gb(δ) ≤ gb(0) we must have
|δ| ≤

6 kx? k22 σ
b2
b d? ).
dist(d,
σ
bI2 out

By convexity, any minimizer of gb must thus lie in the above region, which gives the result when we
recall that minimizers δ of gb are equivalent to minimizers (1 + δ) kx? k22 of G.

B.3

Proof of Proposition 5

We introduce a bit of notation before giving the proof proper. Recall that I out ⊂ [m] denotes the
outliers, or failed measurements, and I in = [m] \ I out the set of i such that bi = hai , x? i2 (the
inliers). Recalling the selected set of indices Isel , we define the shorthand
in
out
Isel
:= I in ∩ Isel and Isel
:= I out ∩ Isel .

for the chosen inliers and outliers.
We decompose the matrix X init into four matrices, each of which we control to guarantee
that X init ≈ In − cd? d? T for some constant c, thus guaranteeing db ≈ d? . Let P = d? d? T and
P⊥ = In − d? d? T be the projection operator onto the span of d? and its orthogonal complement.
Then we may decompose the matrix X init into the four parts
m

X

init

m


 
1 X
1 X
in
in
=
P ai aTi P 1 i ∈ Isel
+
P ai aTi P⊥ + P⊥ ai aTi P 1 i ∈ Isel
m
m
{z
} | i=1
{z
}
| i=1
1
+
m
|

:=Z0
n
X

:=Z1

P⊥ ai aTi P⊥ 1



i∈

in
Isel

i=1

{z

:=Z2

38

1
+
m
} |

m
X

ai aTi 1

(35)


out
i ∈ Isel
.

i=1

{z

:=Z3

}

Let us briefly motivate this decomposition. We expect that Z0 should be small because we choose
indices Isel by taking the smallest values of bi , which should be least correlated with d? (recall the
P = d? d? T ). We expect Z1 to be small because of the independence of the vectors P ai and P⊥ ai
out should be not too large. This
for Gaussian measurement vectors, and Z3 to be small because Isel
leaves Z2 , which (by Gaussianity) we expect to be some multiple of In − d? d? T , which will then
allow us to apply eigenvector perturbation guarantees using the eigengap of the matrix X init .
The rotational invariance of the Gaussian means that it is no loss of generality to assume that
d? = e1 , the first standard basis vector, so for the remainder of the argument we assume this without
further comment. This means that we may decompose ai as ai = [ai,1 ai,2 · · · ai,n ]T = [ai,1 aTi,\1 ]T ,
which we will do without further comment for the remainder of the proof.
We now present four lemmas, each controlling one of the terms Zl in the expansion (35). We
defer proofs of each of the lemmas to the end of this argument. We begin by considering the Z0
term, which (because P is rank one) satisfies
 X

m

1
? 2
in
Z0 =
hai , d i 1 i ∈ Isel
d? d? T .
m
i=1
|
{z
}
:=z0

Recalling the definition (17) of the constants δq and wq in the statement of the proposition, we
have the following lemma.

1
1 Pm
? 2
in
in
Lemma B.3. Define the random quantities z0 = m
i=1 hai , d i 1 i ∈ Isel and z2 = m |Isel |.
Then for t ∈ [0, 1],
!


 m
mt2 δq2
1
P (z2 ≥ z0 + (1 − 2pfail )(1 − t)δq ) ≥ 1 − exp − (1 − 2pfail )m − 2 exp −
− exp −
.
10
4
4wq2
See Sec. B.3.1 for a proof of the lemma. We thus see that it is likely that z0 is substantially smaller
than the rough fraction of inlying indices selected.
We now argue that Z1 is likely to be small because it is the sum of products of independent
vectors.
Lemma B.4. For t ≥ 0 we have
r




 m
n
mt2
P |||Z1 |||op ≥ 2
+ t ≤ exp −
+ exp −
.
m
8
2
See Section B.3.2 for a proof. We can also show that Z2 is well-behaved in the sense that it is
approximately a scaled multiple of (I − d? d? T ).
1
in |. There exists a numerical constant C
Lemma B.5. Let z2 be the random quantity z2 := m
|Isel
such that for t ∈ [0, 1] we have

r

n
? ?T
≥C
P
Z2 − z2 (In − d d )
+t
≤ exp(−mt).
m
op

See Section B.3.3 for a proof of the lemma.
Finally, we control the size of the error matrix Z3 in the expansion (35), which corresponds to
the contribution of the outlying measurements ai that are included in the initialization matrix X init .
We provide two slightly different guarantees, depending on the model (strength) of adversarial noise
ξ assumed.
39

Lemma B.6. Define the random quantity z3 :=
C such that the following hold.

1
out ∩I |
sel
m |I

≤ pfail . There is a numerical constant

(i) Under the independent noise model M1, for all t ∈ [0, 1],
r


n
P |||Z3 − z3 In |||op ≥ C
+ t ≤ exp(−mt).
m
(ii) Under the adversarial noise model M2, for all t ∈ [0, 1],
r


n
+ t ≤ exp(−mt).
P |||Z3 |||op ≥ pfail + C
m
See Section B.3.4 for a proof.
With these four lemmas in hand, we can prove the result of the proposition by applying the
eigenvector perturbation
Lemma A.1. The expansion (35), coupled with the four lemmas defining
1 P
1
in |, and z = 1 |I out |, guarantees that
the constants z0 = m i∈I in hai , d? i2 , z2 = m
|Isel
3
m sel
sel

X

init

= Z0 + Z1 + Z2 + Z3 = z0 d? d? T + z2 (In − d? d? T ) + z3 In + ∆
= (z2 + z3 )In − (z2 − z0 )d? d? T + ∆,

where the perturbation ∆ ∈ Rn×n satisfies
? ?T

|||∆|||op ≤ |||Z1 |||op + Z2 − z2 (In − d d

(
|||Z3 − z3 In |||op
)
+
op
|||Z3 |||op

under model M1
under model M2.

On the event that z2 > z0 , the minimal eigenvector of (z2 + z3 )In − (z2 − z0 )d? d? T is d? with
eigengap z2 − z0 . Applying Lemma A.1 gives that db = argmind∈Sn−1 dT X init d satisfies
q
|||∆|||op
?
− 21
b
b d? i2 ≤ h
i .
(36)
2 dist(d, d ) ≤ 1 − hd,
z2 − z0 − |||∆|||op
+

Applying Lemmas B.4 and B.5, we have for some numerical constant C < ∞ that
r
n
? ?T
|||Z1 |||op + Z2 − z2 (In − d d )
≤C
+ t with probability ≥ 1 − e−mt − e−m/2
m
op
for any t ≥ 0. We now p
consider the two noise models in turn. Under Model M1, Lemma B.6 then
n
implies that |||∆|||op ≤ C m
+ t with probability at least 1−e−mt −e−m/2 . Recalling Lemma B.3, we
fail
have for the constants wq2 and δq > 0 defined in the lemma that z2 ≥ z0 + 1−2p
δq with probability
2
2
2
at least 1 − exp(−c(1 − 2pfail )m) − exp(−cmδq /wq ), where c > 0 is a numerical constant. That is,
the perturbatoin inequality (36) implies that under Model M1, we have with probability at least
1 − exp(−mt) − exp(−c(1 − 2pfail )m) − exp(−cmδq2 /wq2 ) that
p
C n/m + t
− 21
?
b
i ,
2 dist(d, d ) ≤ h
p
(1 − 2pfail )δq − C n/m + t
+

where 0 < c, C < ∞ are numerical constants. Under Model M2, we can bound |||Z3 |||op by pfail +
p
C n/m + t with probability 1−e−mt −e−m/2 (recall Lemma B.6), so that with the same probability
as above, we have
p
pfail + C n/m + t
?
− 21
b
i under Model M2.
2 dist(d, d ) ≤ h
p
(1 − 2pfail )δq − C n/m + t − pfail
+

This is the proposition.
40

B.3.1

Proof of Lemma B.3

As noted earlier, it is no loss of generality
to assume that d? = e1 , the first standard basis vector,
P
1
1
? 2
in
so that using our definitions of z0 = m
i∈I in hai , d i and z2 = m |Isel |, we have
sel

m


1 X
in
z2 − z0 =
(1 − a2i,1 )1 i ∈ Isel
.
m

(37)

i=1

Given that we choose in indices Isel by a median, it is helpful to have the following median concentration result.
iid

2
2
Lemma B.7. Let {Wi }m
i=1 ∼ N(0, 1). Fix p ∈ (0, 1) and choose wq ≥ 0 so that q := P(W ≤ wq ) =
2(1 − Φ(wq )) > p. Then




m(q − p)2
2 m
2
P quantp {Wi }i=1 ≥ wq ≤ exp −
.
2(q − p)/3 + 2q(1 − q)
 2
1 Pm
2
Proof Note that quantp ({Wi2 }m
wq2 ≤ p. Using that
i=1 ) ≥ wq if and only if m
i=1 1 Wi ≤P
2
2
2
2
Var(1 Wi ≤ wq ) = q(1−q), Bernstein’s inequality applied to the i.i.d. sum m
i=1 (1{Wi ≤ wq }−q)
gives the result.

We now control the median of the perturbed vector b ∈ Rm . Since we have |I out | ≤ pfail m, we
have deterministic result


? 2
med ({bi }m
)
≤
quant
{ha
,
x
i
}
,
1
in
i
i∈I
i=1
2(1−pfail )

so by upper bounding the right hand quantity we can upper bound med({bi }), which in turn allows
1
us to control Isel . By the definition of wq and qfail in Eq. (17), which satisfies δ = qfail − 2(1−p
=
fail )
1−2pfail
4(1−pfail ) ,

Lemma B.7 with q = qfail and the fact that |I in | = (1 − pfail )m then implies





? 2
2
? 2
2
? 2
P med({bi }m
i=1 ) ≥ wq kx k2 ≤ P quant(2(1−pfail ))−1 ({hai , x i }i∈I in ) ≥ wq kx k2






3(1 − 2pfail )m
3
m(1 − pfail )δ 2
≤ exp −
= exp −
≤ exp − (1 − 2pfail )m . (38)
2δ/3 + 2δ(1 − δ)
4(2 + 6(1 − δ))
32
We now consider the indices i that are inliers for which hai , d? i2 is small; again letting wq ≥ 0
be defined as in the quantile (17), we define
n
o n
o
Iqin := i ∈ I in | hai , d? i2 ≤ wq2 = i ∈ I in | hai , x? i2 ≤ wq2 kx? k22 .
Lemma B.8. Let the set of inliers Iqin be defined as above, and let δq = 1 − E[W 2 | W 2 ≤ wq2 ] for
W ∼ N(0, 1). Then for all t ≥ 0 we have


!
in |t2
X
2|I
1
q
P  in
(1 − a2i,1 ) ≤ δq − t | Iqin  ≤ exp −
.
|Iq |
wq2
in
i∈Iq

We defer the proof of Lemma B.8, continuing on to give the proof of Lemma B.3.

41

We now integrate out the conditioning in Lemma B.8. Recalling the definition (17) of wq in
3−2pfail
terms of pfail , we have that qfail = 4(1−p
> 12 for pfail ∈ [0, 1/2). By Hoeffding’s inequality we
fail )
have P(|Iqin | ≤ m/8) ≤ e−m/4 because |I in | ≥ (1 − pfail )m ≥ m/2, whence we obtain




X
mt2
1
2


(1 − ai,1 ) ≤ δq − t ≤ exp − 2 P(|Iqin | ≥ m/8) + P(|Iqin | ≤ m/8)
P
|Iqin |
4wq
i∈Iqin


 m
mt2
≤ exp − 2 + exp −
.
4wq
4
Using the notation δq in Lemma B.8, for t ∈ [0, 1] we define the event




X
m
1
2
2
? 2
in
(1
−
a
)
≥
(1
−
t)δ
.
E := med({bi }m
)
≤
w
kx
k
,
|I
|
≥
,
q
i,1
i=1
q
q
2


8 |Iqin |
in
i∈Iq

We immediately find that
!

 m
mt2 δq2
1
P(E) ≥ 1 − exp − (1 − 2pfail )m − 2 exp −
− exp −
10
4
4wq2


(39)

by the preceding display and inequality (38). Recalling the set Iqin = {i ∈ I in : hai , d? i2 ≤ wq2 },
in ⊂ I in (because med(b) ≤ w 2 kx? k2 ).
we have on the event E that the selected inliers satisfy Isel
q
q
2
in as the smallest b in the sample, we have that
Because of our selection mechanism with Isel
i
1 X
1 X
2
(1 − a2i,1 ).
(1
−
a
)
≥
i,1
in |
in |
|I
|Isel
q
in
in
i∈Iq |

i∈Isel

in | ≥ |I in | + |I | − m ≥
Moreover, on the event E the rightmost sum is positive, and using that |Isel
sel
1
( 2 − pfail )m, we obtain that on E we have

1 X
1 − 2pfail X
1 − 2pfail X
1 − 2pfail
2
(1
−
a
)
≥
δq .
(1 − a2i,1 ) ≥
(1 − a2i,1 ) ≥
i,1
in
in
m
2|Iq |
2
2|Isel |
in
in
in
i∈Iq

i∈Isel

i∈Isel

Recalling expression (37) thus gives Lemma B.3.
Proof of Lemma B.8
Fix any set of indices I0 ⊂ [m], and note that by Hoeffding’s inequality
for bounded random variables we have




X
1
2|I0 |t2
2
2
2
2
2
2


P
ai,1 − E[ai,1 | ai,1 ≤ wq ] ≥ t | ai,1 ≤ wq for i ∈ I0 ≤ exp −
|I0 |
wq2
i∈I0

for t ≥ 0. Recalling the definition δq := 1 − E[W 2 | W 2 ≤ wq2 ] for W ∼ N(0, 1), we thus find that
for any index set I0 ⊂ [m], we have


X
1
P  in
(1 − a2i,1 ) ≤ δq − t | Iqin = I0 
|Iq |
i∈Iqin




2
X
1
2|I
|t
0
2
2
2
= P
(1 − ai,1 ) ≤ δq − t | ai,1 ≤ wq for i ∈ I0  ≤ exp −
.
|I0 |
wq2
i∈I0

42

Noticing that the random vectors {ai }i∈I0 are independent of {ai }i6∈I0 , we have for any measurable
set C ⊂ R that




X
X
(1 − a2i,1 ) ∈ C | Iqin = I0  = P 
(1 − a2i,1 ) ∈ C | a2i,1 ≤ wq2 for i ∈ I0  .
P
i∈Iqin

i∈I0

Combining the preceding two displays yields Lemma B.8.

B.3.2

Proof of Lemma B.4

By our assumption (w.l.o.g.) that d? = e1 , we have
n

|||Z1 |||op =


1 X
in
ai,1 ai,\1 1 i ∈ Isel
m
i=1

(40)
2

iid

Letting {a0i }m
i=1 ∼ N(0, In ) be an independent collection of vectors, the collections {ai }i∈I in and
in ⊂ I in , that for any
{ξi }i∈I out are independent, as are ai,1 and ai,\1 for each i. Thus, because Isel
n−1
measurable set C ⊂ R
we have




X
X
in 
P
ai,1 ai,\1 ∈ C | {ai,1 }i∈I in , {ξi }i∈I out  = P 
ai,1 a0i,\1 ∈ C | {ai,1 }i∈I in , Isel
. (41)
in
i∈Isel

in
i∈Isel

Now, we use the standard result [25] that if f is an L-Lipschitz function with respect to the `2 -norm,
t2
then for any standard Gaussian vector Z we have P(f (Z) − E[f (Z)] ≥ t) ≤ exp(− 2L
2 ) for all t ≥ 0.
1 P
2
2
b
Thus, defining the random Lipschitz constant L = m i∈I in ai,1 we obtain
"
#
m
2

1 X
1 X 2
n − 1 b2
in
in
E
L ,
ai,1 a0i,\1 1 i ∈ Isel
| Isel
, {ai,1 }i∈I in = 2
ai,1 (n − 1) ≤
m
m
m
2
in
i=1

i∈Isel

and thus for t ≥ 0
m

P


1 X
in
ai,1 a0i,\1 1 i ∈ Isel
m
i=1

2

!
r


n−1
mt2
in
b
≥L
+ t | {ai,1 }i∈I in , Isel ≤ exp −
.
b2
m
2L

b is {ai,1 }i∈I in -measurable, and P(L
b ≥ 2) ≤ exp(−m/2), again by the Lipschitz concenMoreover, L
tration of Gaussian random variables. We thus apply Eq. (41) and obtain


r
X
1
n−1
P
ai,1 ai,\1 ≥ 2
+ t
m
m
2
in
i∈Isel


r
X
1
n−1
b ≤ 2 P(L
b ≤ 2) + P(L
b ≥ 2)
≤ P
ai,1 ai,\1 ≥ 2
+t|L
m
m
2
in
i∈Isel


 m
mt2
≤ exp −
+ exp −
.
8
2
Recalling the equality (41) on |||Z1 |||op shows that the previous display gives the lemma.
43

B.3.3

Proof of Lemma B.5

Our first observation is simply the definition of the projection operator P⊥ , which gives
Zm,2 − z2 (In − d? d? T )

op

1 X
ai,\1 aTi,\1 − z2 In−1
m
in

=

i∈Isel

.
op

iid

As in the proof of Lemma B.4, we let {a0i }m
i=1 ∼ N(0, In ) be an independent copy of the ai . The
independence of {ai }i∈I in and {ξi }i∈I out imply that for any measurable set C we have





X
X 
T
in 
P
.
(ai,\1 aTi,\1 − z2 In−1 ) ∈ C | {ai,1 }i∈I in , {ξi }i∈I out  = P 
a0i,\1 a0i,\1 − In−1 ∈ C | Isel
in
i∈Isel

in
i∈Isel

As {a0i,\1 }m
i=1 are i.i.d. Gaussian vectors, Lemma 3.1 implies that for t ∈ [0, 1]:


r
m 

X

1
n
T
in
in 
≥C
P
a0i,\1 a0i,\1 − In−1 1 i ∈ Isel
+ t | Isel
≤ exp(−mt)
m
m
i=1

op

for some constant C. The claim of the lemma follows by integration of the above inequality.
B.3.4

Proof of Lemma B.6

out = I out ∩ I .
For notational convenience, let us denote the selected outlying indices by Isel
sel
We begin our proof by considering Model M1, in which case the proof is essentially identical to
Lemma B.5 (see Sec. B.3.3). Indeed, we have the identity
m

 
1 X
out
ai aTi − In 1 i ∈ Isel
m

|||Z3 − z3 In |||op =

i=1

,
op
iid

and so following the completely parallel route of introducing the independent sample {a0i }m
i=1 ∼
N(0, In ), we have for any t ≥ 0 that
i


h 
out
P |||Z3 − z3 In |||op ≥ t = E P |||Z3 − z3 In |||op ≥ t | Isel
"
!#
m

1 X 0 0T
out
out
=E P
≥ t | Isel
.
(ai ai − In )1 i ∈ Isel
m
op
i=1

out gives the result.
Applying Lemma 3.1 to the vectors a0i 1 i ∈ Isel
In the case in which we assume Model M2, we have the deterministic upper bound



m

|||Z3 |||op =


1 X
out
ai aTi 1 i ∈ Isel
m
i=1

m

≤
op


1 X
ai aTi 1 i ∈ I out
m
i=1

.
op

As ai 1 i ∈ I out is an O(1)-sub-Gaussian vector, Lemma 3.1 implies


r
m
X
 
1
n
out 
P
ai aTi − In 1 i ∈ I out
≥C
+ t | Isel
≤ exp(−mt).
m
m


i=1

op

Since |I out |/m = pfail , the desired claim follows via the triangle inequality.
44

C
C.1

Technical results
Proof of Lemma 3.2

Rc
Let p(z) denote the density of z and pc (z 0 ) = p(z 0 )/ −c p(z)dz for z 0 ∈ [−c, c]. Let Fc (z) =
Rz
0
0
−c pc (z )dz be the CDF of pc , and define the mapping T : [−c, c] → [−c, c] by
T (u) = Fc−1



u+c
2c


= z for the z s.t. Fc (z) =

u+c
.
2c

Evidently, for Z ∼ Pc we have T (Z) ∼ Uni[−c, c], and by the symmetry and monotonicity properties
of p and pc we have |T (z)| ≥ |z| for z ∈ [−c, c]. In particular, letting U ∼ Uni[−c, c], we have
Var(U ) = E[U 2 ] = E[T (Z)2 | |Z| ≤ c] ≥ E[Z 2 | |Z| ≤ c].
Using that Var(U ) =

C.2

1
2c

Rc

2
−c u du

=

c2
3

gives the result.

Proof of Lemma B.1

The proof is an essentially standard concentration and covering number argument, with a few minor
wrinkles. First, we note that
(i)

Var(| hu, ai i hv, ai i |) ≤ E[| hu, ai i hv, ai i |2 ] ≤ sup E[hu, ai i4 ] ≤ 2eσ 4 ,
u∈Sn−1

where inequality (i) is Cauchy-Schwarz and the final inequality follows Lemma 3.1. Thus, the
lower tail Bernstein’s inequality applied to the positive random variables | hu, ai i hv, ai i | ≥ 0 with
variance bounded by 2eσ 4 implies that
!


m
1 X
mt2
P
.
(42)
| hu, ai i hv, ai i | − κst (u, v) ≤ −t ≤ exp −
m
4eσ 4
i=1

Pm

1
n
Define
have
Pm Zu,v := m i=1 | hu, ai i hv, ai i | for shorthand. Using that for vectors w, x ∈0 R0 wen−1
,
i=1 | hai , wi hai , xi | ≤ kAwk2 kAxk2 by Cauchy-Schwarz, we see that for any u, v, u , v ∈ S
m

|Zu,v − Zu0 ,v0 | ≤
≤

1 X
| hai , ui hai , vi − ai , u0
m
1
m

i=1
m
X

ai , v 0 |


| ai , u − u0 hai , vi | + | hai , ui ai , v − v 0 |

i=1

1
A(u − u0 )
≤
m

1
kAvk2 +
A(v − v 0 )
2
m

kAuk2 ≤
2

|||A|||2op
m

u − u0

+ v − v0
2


2

Now, let N be an -cover of the sphere Sn−1 , which we use to control inf u,v∈Sn−1 {Zu,v − E[Zu,v ]}.
Using the previous display, we obtain by the definition of an -cover argument that
!


2 AT A op
P
inf {Zu,v − κst (u, v)} ≤ −t ≤ P min {Zu,v − κst (u, v)} ≤ −t +
 .
u,v∈N
m
u,v∈Sn−1

45

.

From Lemma 3.1, we know that |||AT A|||op is well-concentrated, and thus by considering the event
that |||AT A|||op & mσ 2 that for some numerical constant C < ∞ we have


P
inf {Zu,v − κst (u, v)} ≤ −t
u,v∈Sn−1




2
≤ P min {Zu,v − κst (u, v)} ≤ −t + 2Cσ  + P AT A op ≥ Cmσ 2
u,v∈N

2 !
cm t − Cσ 2  +
2
+ exp(−m),
(43)
≤ card(N ) exp −
σ4
where the first term comes from Bernstein’s inequality (42) and the second from Lemma 3.1.
Now, let us assume that N is an -cover of Sn−1 with minimal cardinality, which by standard
volume arguments [37, Lemma 2] satisfies N () := card(N ) ≤ (1 + 2/)n for  > 0. Noting that
(1 + 2/)2n ≤ exp(n/) for  ≤ 0 := .21398, we may replace inequality (43) with

 !


2 2
n cm t − Cσ  +
P
inf {Zu,v − κst (u, v)} ≤ −t ≤ exp
+ e−m
−

σ4
u,v∈Sn−1
p
valid for all  ≤ 0 . Now, if we set  = 3 n/m and define t = σ 2 t + (C + 1)σ 2 , then we find that if
m/n ≥ −3
0 , we have
r 



2
2
2 3 n
P
inf {Zu,v − κst (u, v)} ≤ −σ t − (C + 1)σ
≤ exp cmt + e−m .
m
u,v∈Sn−1
This is the desired result.

C.3

Properties of Gaussian Random Variable

Lemma C.1. Let W ∼ N(0, 1). Then for all c ∈ R+


1
E W 2 | W 2 ≤ c ≤ 1 − cP(W 2 > c).
2
Proof By the law of total probability, we have
E[W 2 ] = E[W 2 | W 2 ≤ c]P(W 2 ≤ c) + E[W 2 | W 2 > c]P(W 2 > c).
Using that E[W 2 ] = 1 yields

1 − E[W 2 | W 2 ≤ c] = P(W 2 > c) E[W 2 | W 2 > c] − E[W 2 | W 2 ≤ c] .
For t ∈ [0, c], we have P(W 2 ∈ [c − t, c]) ≤ P(W 2 ∈ [0, t]), so that for such t we have
P(W 2 ≥ c−t | W 2 ≤ c) ≤ P(W 2 ≤ t | W 2 ≤ c), or P(W 2 ≥ t | W 2 ≤ c)+P(W 2 ≥ c−t | W 2 ≤ c) ≤ 1.
Performing the standard change of variables to compute E[W 2 | W 2 ≤ c], we thus obtain
Z c
2
2
E[W | W ≤ c] =
P(W 2 ≥ t | W 2 ≤ c)dt
Z0 c

1
=
P(W 2 ≥ t | W 2 ≤ c) + P(W 2 ≥ c − t | W 2 ≤ c) dt
2
Z0 c
1
c
≤
dt = .
2
0 2
Using that E[W 2 | W 2 > c] ≥ c gives us the lemma.

46

C.4

Combinatorics

Lemma C.2. Then for any c ≥ 0, the VC-dimension of the collection of sets
G := {{x ∈ Rn | | hx, ui hx, vi | ≥ c} | u ∈ Rn , v ∈ Rn }
is upper bounded by 2(n + 2).
Proof

Define the collection of scaled spheres
nn
o
o
Gsphere :=
x ∈ Rn | a kxk2 + hu, xi + b ≥ 0 | a ∈ R, b ∈ R, u ∈ Rn .

Then a standard result of Dudley [16] is that the VC-dimension of Gsphere is bounded by n + 2. Of
course, we have the containment
G ⊂ {G1 ∪ G2 | G1 , G2 ∈ Gsphere } .
Standard stability guarantees of VC-dimension [36, Lemma 2.6.17] imply that the VC-dimension
of G is bounded by twice that of Gsphere .

References
[1] R. Balan, P. Casazza, and D. Edidin. On signal reconstruction without phase. Applied and
Computational Harmonic Analysis, 20:345–356, 2006.
[2] A. Bandeira, J. Cahill, D. Mixon, and A. Nelson. Saving phase: Injectivity and stability for
phase retrieval. Applied and Computational Harmonic Analysis, 37(1):106–125, 2014.
[3] P. L. Bartlett and S. Mendelson. Rademacher and Gaussian complexities: Risk bounds and
structural results. Journal of Machine Learning Research, 3:463–482, 2002.
[4] A. Ben-Tal and A. Nemirovski. Lectures on Modern Convex Optimization. SIAM, 2001.
[5] J. Burke. Descent methods for composite nondifferentiable optimization problems. Mathematical Programming, 33:260–279, 1985.
[6] J. Burke and M. Ferris. A Gauss-Newton method for convex composite optimization. Mathematical Programming, 71:179–194, 1995.
[7] E. Candès, T. Strohmer, and V. Voroninski. PhaseLift: exact and stable signal recovery from
magnitude measurements via convex programming. Communications on Pure and Applied
Mathematics, 66(8), 2013.
[8] E. J. Candès and X. Li. Solving quadratic equations via PhaseLift when there are about as
many equations as unknowns. Foundations of Computational Mathematics, 14(5):1017–1026,
2014.
[9] E. J. Candès, Y. Eldar, T. Strohmer, and V. Voroninski. Phase retrieval via matrix completion.
SIAM Journal on Imaging Sciences, 6(1):199–225, 2013.
[10] E. J. Candès, X. Li, and M. Soltanolkotabi. Phase retrieval via Wirtinger flow: Theory and
algorithms. IEEE Transactions on Information Theory, 61(4):1985–2007, 2015.
47

[11] A. Chai, M. Moscoso, and G. Papanicolaou. Array imaging using intensity-only measurements.
Inverse Problems, 27(1):015005, 2011.
[12] Y. Chen and E. Candès. Solving random quadratic systems of equations is nearly as easy as
solving linear systems. Communications on Pure and Applied Mathematics, To appear, 2015.
[13] D. Drusvyatskiy and A. Lewis. Error bounds, quadratic growth, and linear convergence of
proximal methods. arXiv:1602.06661 [math.OC], 2016.
[14] D. Drusvyatskiy and C. Paquette. Efficiency of minimizing compositions of convex functions
and smooth maps. arXiv:1605.00125 [math.OC], 2016.
[15] D. Drusvyatskiy, A. Ioffe, and A. Lewis. Nonsmooth optimization using Taylor-like models:
error bounds, convergence, and termination criteria. arXiv:1610.03446 [math.OC], 2016. URL
https://arxiv.org/abs/1610.03446.
[16] R. M. Dudley. Balls in Rk do not cut all subsets of k + 2 points. Advances in Mathematics,
31(3):306–308, 1979.
[17] Y. Eldar and S. Mendelson. Phase retrieval: Stability and recovery guarantees. Applied and
Computational Harmonic Analysis, 36(3):473–494, 2014.
[18] J. R. Fienup. Reconstruction of an object from the modulus of its Fourier transform. Optics
Letters, 3(1):27–29, 1978.
[19] J. R. Fienup. Phase retrieval algorithms: a comparison. Applied Optics, 21(15):2758–2769,
1982.
[20] R. Fletcher and G. A. Watson. First and second order conditions for a class of nondifferentiable
optimization problems. Mathematical Programming, 18:291–307, 1980.
[21] R. W. Gerchberg and W. O. Saxton. A practical algorithm for the determination of phase
from image and diffraction. Optik, 35:237–246, 1972.
[22] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning. Springer,
second edition, 2009.
[23] J. Hiriart-Urruty and C. Lemaréchal.
Springer, New York, 1993.

Convex Analysis and Minimization Algorithms I.

[24] P. J. Huber and E. M. Ronchetti. Robust Statistics. John Wiley and Sons, second edition,
2009.
[25] M. Ledoux. The Concentration of Measure Phenomenon. American Mathematical Society,
2001.
[26] S. Mendelson. Learning without concentration. In Proceedings of the Twenty Seventh Annual
Conference on Computational Learning Theory, 2014.
[27] P. Netrapalli, P. Jain, and S. Sanghavi. Phase retrieval using alternating minimization. In
Advances in Neural Information Processing Systems 26, 2013.
[28] N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):
123–231, 2013.
48

[29] N. Parikh and S. Boyd. Block splitting for distributed optimization. Mathematical Programming Computation, 6:77–102, 2014.
[30] F. Pi, H. Zhang, and P. Guo. RNA nanoparticles in cancer cells. National Cancer Institute
Visuals Online, 2016. URL https://visualsonline.cancer.gov/details.cfm?imageid=
11167. Accessed April 1, 2017.
[31] Y. Schechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev. Phase
retrieval with application to optical imaging. IEEE Signal Processing Magazine, pages 87–109,
May 2015.
[32] G. W. Stewart and J.-G. Sun. Matrix Perturbation Theory. Academic Press, 1990.
[33] L. N. Trefethen and D. Bau III. Numerical Linear Algebra. SIAM, 1997.
[34] J. A. Tropp. Improved analysis of the subsampled randomized Hadamard transform. Advances
in Adaptive Data Analysis, 3(1–2):115–126, 2011. Special issue on Sparse Representation of
Data and Images.
[35] M. Udell, K. Mohan, D. Zeng, J. Hong, S. Diamond, and S. Boyd. Convex optimization in
Julia. In First Workshop on High Performance Technical Computing in Dynamic Languages,
pages 18–28. IEEE, 2014.
[36] A. W. van der Vaart and J. A. Wellner. Weak Convergence and Empirical Processes: With
Applications to Statistics. Springer, New York, 1996.
[37] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. In Compressed
Sensing: Theory and Applications, chapter 5, pages 210–268. Cambridge University Press,
2012.
[38] I. Waldspurger, A. d’Aspremont, and S. Mallat. Phase recovery, maxcut and complex semidefinite programming. Mathematical Programming, Series A, 149(1–2):47–81, 2015.
[39] G. Wang, G. B. Giannakis, and Y. C. Eldar. Solving systems of random quadratic equations
via truncated amplitude flow. arXiv:1605.08285 [stat.ML], 2016.
[40] H. Zhang, Y. Chi, and Y. Liang. Provable non-convex phase retrieval with outliers: Median
truncated Wirtinger flow. In Proceedings of the 32nd International Conference on Machine
Learning, 2016.

49

