Fundamental Limits of Weak Recovery
with Applications to Phase Retrieval
Marco Mondelli∗ and

Andrea Montanari†

arXiv:1708.05932v2 [stat.ML] 28 Sep 2017

October 2, 2017

Abstract
In phase retrieval we want to recover an unknown signal x ∈ Cd from n quadratic measurements of
the form yi = |hai , xi|2 + wi where ai ∈ Cd are known sensing vectors and wi is measurement noise.
We ask the following weak recovery question: what is the minimum number of measurements n needed
to produce an estimator x̂(y) that is positively correlated with the signal x? We consider the case of
Gaussian vectors ai . We prove that – in the high-dimensional limit – a sharp phase transition takes place,
and we locate the threshold in the regime of vanishingly small noise. For n ≤ d − o(d) no estimator
can do significantly better than random and achieve a strictly positive correlation. For n ≥ d + o(d) a
simple spectral estimator achieves a positive correlation. Surprisingly, numerical simulations with the
same spectral estimator demonstrate promising performance with realistic sensing matrices. Spectral
methods are used to initialize non-convex optimization algorithms in phase retrieval, and our approach
can boost the performance in this setting as well.
Our impossibility result is based on classical information-theory arguments. The spectral algorithm
computes the leading eigenvector of a weighted empirical covariance matrix. We obtain a sharp characterization of the spectral properties of this random matrix using tools from free probability and generalizing a recent result by Lu and Li. Both the upper and lower bound generalize beyond phase retrieval to
measurements yi produced according to a generalized linear model. As a byproduct of our analysis, we
compare the threshold of the proposed spectral method with that of a message passing algorithm.

1

Introduction

In this work, we consider the problem of recovering a signal x of dimension d, given n generalized linear
measurements. More specifically, the measurements are drawn independently according to the conditional
distribution
yi ∼ p(y | hx, ai i),
i ∈ {1, . . . , n},
(1)
where h·, ·i denotes the inner product, {ai }1≤i≤n is a set of known sensing vector, and p(· | hx, ai i)
is a known probability density function. This model appears in many problems in signal processing and
statistical estimation, e.g., photon-limited imaging [UE88, YLSV12], signal recovery from quantized measurements [RG01], and phase retrieval [Fie82, SEC+ 15]. For the problem of phase retrieval, the model (1)
is specialized to
yi = |hx, ai i|2 + wi ,
i ∈ {1, . . . , n} ,
(2)
∗
†

Department of Electrical Engineering, Stanford University
Department of Electrical Engineering and Department of Statistics, Stanford University

1

where wi is noise. Applications of phase retrieval arise in several areas of science and engineering, including
X-ray crystallography [Mil90, Har93], microscopy [MISE08], astronomy [FD87], optics [Wal63], acoustics
[BCE06], interferometry [DJ17], and quantum mechanics [Cor06].
Popular methods to solve the phase retrieval problem are based on semi-definite programming relaxations [CESV15, CLS15a, CSV13, WdM15]. However, these algorithms rapidly become prohibitive from a
computational point of view when the dimension d of the signal increases, which makes them impractical
in most of the real-world applications. For this reason, several algorithms have been developed in order
to solve directly the non-convex least-squares problem, including the error reduction schemes dating back
to Gerchberg-Saxton and Fienup [Ger72, Fie82], alternating minimization [NJS13], approximate message
passing [SR15], Wirtinger Flow [CLS15b], iterative projections [LGL15], the Kaczmarz method [Wei15],
and a number of other approaches [CC17, ZL16, CLM16, WGE16, WG16, Sol17, DR17, WGSC17]. Furthermore, recently a convex relaxation that operates in the natural domain of the signal was independently
proposed by two groups of authors [GS16, BR17]. All these techniques require an initialization step, whose
goal is to provide a solution x̂ that is positively correlated with the unknown signal x. To do so, spectral
methods are widely employed: the estimate x̂ is given by the principal eigenvector of a suitable matrix
constructed from the data. A similar stategy (initialization step followed by an iterative algorithm) has
proved successful for many other estimation problems, e.g., matrix completion [KMO10, JNS13], blind
deconvolution [LLJB17, LLSW16], sparse coding [AGMM15] and joint alignment from pairwise noisy observations [CC16].
We focus on a regime in which both the number of measurement n and the dimension of the signal d
tend to infinity, but their ratio n/d tends to a positive constant δ. The weak recovery problem requires to
provide an estimate x̂(y) that has a positive correlation with the unknown vector x:


|hx̂(y), xi|
lim inf E
> ,
(3)
n→∞
kx̂(y)k2 kxk2
for some  > 0.
In this paper, we consider either x ∈ Rd or x ∈ Cd and assume that the measurement vectors ai are
standard Gaussian (either real or complex). In the general setting of model (1), we present two types of
results:
1. We develop an information-theoretic lower bound δ` : for δ < δ` , no estimator can output non-trivial
estimates. In other words, the weak recovery problem cannot be solved.
2. We establish an upper bound δu based on a spectral algorithm: for δ > δu , we can achieve weak
recovery (see (3)) by letting x̂ be the principal eigenvector of a matrix suitably constructed from the
data.
The values of the thresholds δ` and δu depend on the conditional distribution p(· | hx, ai i). For the special
case of phase retrieval (see (2)), we evaluate these bounds and we show that they coincide in the limit of
vanishing noise.
Theorem. Let x ∈ Cd be such that kxk22 = d, and assume that {ai }1≤i≤n ∼i.i.d. CN(0d , Id /d). Let
y ∈ Rn be given by (2), with {wi }1≤i≤n ∼ N(0, σ 2 ), and n, d → ∞ with n/d → δ ∈ (0, +∞). Then,
• For δ < 1, no algorithm can provide non-trivial estimates on x;
• For δ > 1, there exists σ0 (δ) > 0 and a spectral algorithm that returns an estimate x̂ satisfying (3),
for any σ ∈ [0, σ0 (δ)].
2

When x ∈ Rd with kxk22 = d and {ai }1≤i≤n ∼i.i.d. N(0d , Id /d), we show that an analogous result
holds and that the threshold moves from 1 to 1/2. This is reminiscent of how the injectivity thresholds are
δ = 4 and δ = 2 in the complex and the real case, respectively [BCE06, BCMN14, CEHV15]. A possible
intuition for this halving phenomenon comes from the fact that the complex problem has twice as many
variables but the same amount of equations of the real problem. Hence, it is reasonable that the complex
case requires twice the amount of data with respect to the real case.
The lower bound is proved by estimating the conditional entropy via the second moment method.
As in earlier work (see Section 1.1), the spectral algorithm computes the eigenvector corresponding to
the largest eigenvalue of a matrix of the form:
n

1X
Dn =
T (yi )ai a∗i ,
n

(4)

i=1

where T : R → R is a pre-processing function. For δ large enough (and a suitable choice of T ), we
expect the resulting eigenvector x̂(y) to be positively correlated with the true signal x. The recent paper
[LL17] computed exactly the threshold value δu , under the assumption that the measurement vectors are real
Gaussian, and T is non-negative.
Here, we generalize the result of [LL17] by removing the assumption that T (y) ≥ 0 and by considering
the complex case. Armed with this result, we compute the optimal1 pre-processing function Tδ∗ (y) for the
general model (1). Our upper bound δu is the phase transition location for this optimal spectral method. In
the case of phase retrieval (as σ → 0), this pre-processing function is given by
Tδ∗ (y) =

y−1
√
,
y+ δ−1

(5)

and achieves weak recovery for any δ > δu = 1. In the limit δ ↓ 1, this converges to the limiting function
T ∗ (y) = 1 − (1/y).
While the expression (5) is remarkably simple, it is somewhat counter-intuitive. Earlier methods [CLS15b,
CC15,LL17] use T (y) ≥ 0 and try to extract information from the large values of yi . The function (5) has a
large negative part for small y, in particular when δ is close to 1. Furthermore, it extracts useful information
from data points with yi small. One possible interpretation is that the points in which the measurement
vector is basically orthogonal to the unknown signal are not informative, hence we penalize them.
Our analysis applies to Gaussian measurement matrices. However, the proposed spectral method works
well also on real images and realistic measurement matrices. To illustrate this fact, in Figure 1 we test our
algorithm on a digital photograph of the painting “The birth of Venus” by Sandro Botticelli. We consider
a type of measurements that falls under the category of coded diffraction patterns (CDP) [CLS15a, CC17]:
the measurement matrix is given by the product of δ copies of a Fourier matrix and a diagonal matrix with
entries i.i.d. and uniform in {1, −1, i, −i}, where i denotes the imaginary unit. We compare our method
with the truncated spectral initialization proposed in [CC17], which consists in discarding the measurements larger than an assigned threshold and leaving the others untouched. The proposed choice of the
pre-processing function allows to recover a good estimate of the original image already when δ = 4, while
the truncated spectral initialization of [CC17] requires δ = 12 to obtain similar results. As mentioned above,
these estimates are to be used as initializations in a non-convex optimization algorithm that improves the
reconstruction error.
The rest of the paper is organized as follows. In Section 2, after introducing the necessary notation, we
define formally the problem. We then state our general information-theoretic lower bound and our spectral
1

Here optimality is understood with respect to the weak recovery threshold.

3

(a) Original image.

(b) proposed – δ = 4.

(c) truncated – δ = 4.

(d) proposed – δ = 6.

(e) truncated – δ = 6.

(f) proposed – δ = 12.

(g) truncated – δ = 12.

Figure 1: Performance comparison between the proposed spectral method and the truncated spectral initialization of [CC17] for the recovery of a digital photograph from coded diffraction patterns.

4

upper bound for the case of complex signal x and complex measurement vectors ai . The main results for
the real case are stated in Section 3. In Sections 4 and 5, we present the proof of the information-theoretic
lower bound and of the spectral upper bound, respectively. In Section 6, we compare the spectral approach
to a message passing algorithm. In particular we show that the latter cannot have a better threshold than
δu and that δu is the threshold for a linearized version of message passing. In Section 7, we present some
numerical simulations that illustrate the behavior of the proposed spectral method for the phase retrieval
problem. The proofs of several results are deferred to the various appendices.

1.1

Related Work

Precise asymptotic information on high-dimensional regression problems has been obtained by several
groups in recent years [DMM11,BM12,OTH13,BLM15,DM16,Kar13,SC16,ZK16]. In particular, informationtheoretically optimal estimation was considered for compressed sensing [DJM13], and random linear estimation [RP16, BMDK17]. Minimax optimal estimation is considered, among others, in [DMM11, SC16,
VJ17].
The performance of the spectral methods for phase retrieval was first considered in [NJS13]. In the
present notation, [NJS13] uses T (y) = y and proves that there exists a constant c1 such that weak recovery
can be achieved for n > c1 · d · log3 d. The same paper also gives an iterative procedure to improve over
the spectral method, but the bottleneck is in the spectral step. The sample complexity of weak recovery
using spectral methods was improved to n > c2 · d · log d in [CSV13] and then to n > c3 · d in [CC17], for
some constants c2 and c3 . Both of these papers also prove guarantees for exact recovery by suitable descent
algorithms. The guarantees on the spectral initialization are proved by matrix concentration inequalities, a
technique that typically does not return exact threshold values. Furthermore, in [GS16], the authors introduce the PhaseMax relaxation and prove an exact recovery result for phase retrieval, which depends on the
correlation between the true signal and the initial estimate given to the algorithm. By using our spectral
method to obtain the initial estimate, it should be possible to improve the existing upper bounds on the
number of samples needed for exact recovery.
As previously mentioned, our analysis of spectral methods builds on the recent work of Lu and Li
[LL17] that compute the exact spectral threshold for a matrix of the form (4) with T (y) ≥ 0. Here we
generalize this result to signed pre-processing functions T (y), and construct a function of this type that
achieves the information-theoretic threshold for phase retrieval. Our proof indeed implies that non-negative
pre-processing functions lead to an unavoidable gap with respect to the ideal threshold.
Finally, while this paper was under completion, two works appeared that address related problems.
In [BKM+ 17], the authors characterize the information-theoretically optimal estimation error for a broad
class of models of the form (1). However, note that this analysis does not prove – in general – the existence
of an efficient estimation algorithm (for instance in the case of phase retrieval). The paper [DL17] studies
the PhaseMax approach [GS16, BR17] to phase retrieval and uses the non-rigorous replica method from
statistical physics to derive exact thresholds for this algorithm.

2
2.1

Main Results: Complex Case
Notation and System Model

We use [n] as a shortcut for {1, . . . , n}. We use upper-case letters (e.g., X, Y, Z, . . .) to denote random
variables when we are taking operators such as expectation, variance or mutual information. We denote by
0n the vector consisting of n 0s. Given a vector x, we denote by kxk2 its `2 norm. Given a matrix A,
5

we denote by kAkF its Frobenius norm, by kAkop its operator norm, by AT its transpose, and by A∗ its
P
conjugate transpose. Given two vectors x, y ∈ Cd , we denote by hx, yi = di=1 xi yi∗ their scalar product.
We take logarithms in the natural basis and we measure entropies in nats. Given c ∈ C, we denote by
P
a.s.
<(c) and =(c) its real and imaginary part, respectively. We use −→ and −→ to denote the convergence in
probability and the almost sure convergence, respectively.
√
Let x ∈ Cd be chosen uniformly at random on the d-dimensional complex sphere with radius d, i.e.,
√
x ∼ Unif( dSd−1
(6)
C ).
Let the sensing vectors {ai }1≤i≤n , with ai ∈ Cd , be independent and identically distributed according
to a circularly-symmetric complex normal distribution with variance 1/d, i.e.,
{ai }1≤i≤n ∼i.i.d. CN(0d , Id /d).

(7)

Given gi = hx, ai i, the vector of measurements y ∈ Rn is obtained by drawing each component
independently according to the following distribution:
yi ∼ p(y | |gi |),

i ∈ [n].

(8)

For the special case of phase retrieval, the measurements are given by the squared scalar product corrupted
by additive Gaussian noise with variance σ 2 , i.e.,


1
(y − |gi |2 )2
.
(9)
pPR (y | |gi |) = √ exp −
2σ 2
σ 2π
Let δn = n/d and assume that, as n → ∞, δn → δ for some δ ∈ (0, ∞).

2.2

Information-Theoretic Lower Bound

The main result of this section establishes the following: there is a critical value δ` such that, for any
δ < δ` , the optimal estimator has the same performance as a trivial estimator that does not have access to
any measurement. The value of δ` depends on the distribution (8) of the measurements and we provide an
expression to compute it.
In order to state formally the result, we need to introduce a few definitions. Consider the function
f : [0, 1] → R, given by
Z
EG1 ,G2 {p(y | |G1 |)p(y | |G2 |)}
f (m) =
dy,
(10)
EG {p(y | |G|)}
R
with


G ∼ CN(0, 1),

(G1 , G2 ) ∼ CN 02 ,



1 c
c∗ 1


,

(11)

and m = |c|2 . Note that the RHS of (10) depends only on m = |c|2 . Indeed, by applying the transformation (G1 , G2 ) → (eiθ1 G1 , eiθ2 G2 ), f (m) does not change, but the correlation coefficient c is mapped into
cei(θ1 −θ2 ) . A more explicit formula for f (m) is provided by Lemma 6 in Appendix A. Furthermore, set
Fδ (m) = δ log f (m) + log(1 − m).

6

(12)

Note that, when m = 0, G1 and G2 are independent. Hence, f (0) = 1, which implies that Fδ (0) = 0 for
any δ > 0. We define the information-theoretic threshold δ` as the largest value of δ such that the maximum
of Fδ (m) is attained at m = 0, i.e.,
δ` = sup{δ | Fδ (m) < 0 for m ∈ (0, 1]}.

(13)

Let us now define the error metric. The setting is the following: we observe the vector of n measurements
y and, given a new sensing vector an+1 , we want to estimate some function φ(|hx, an+1 i|) given by
Z
ϕ(y)p(y | |hx, an+1 i|) dy.
(14)
φ(|hx, an+1 i|) =
R

Then, the minimum mean square error is defined as

2 

MMSE(δn ) = E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|) Y , {Ai }1≤i≤n
,

(15)

where E {φ(|hX, An+1 i|) | Y , {Ai }1≤i≤n } represents the optimal estimator of the quantity φ(|hx, an+1 i|)
and the expectation of the square error is to be intended over all the randomness of the system, i.e., over
X, An+1 , Y , and {Ai }1≤i≤n . Note that this error metric depends on the choice of the function φ. Furthermore, observe that, if we do not have access to the vector of measurements Y , the trivial estimator
E {φ(|hX, An+1 i|)} has a mean square error given by

2 


= Var φ(|hX, An+1 i|) .
(16)
E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|)
At this point we are ready to state our main result, which is proved in Section 4.1.
Theorem 1 (Information-Theoretic Lower Bound for General Complex Sensing Model). Let x, {ai }1≤i≤n+1 ,
and y be distributed according to (6), (7), and (8), respectively. Let n/d → δ and define δ` as in (13). Furthermore, assume that the function ϕ that appears in (14) is bounded. Then, for any δ < δ` , we have
that

lim MMSE(δn ) = Var φ(|hX, An+1 i|) .
(17)
n→∞

Let us point out that the requirement that the function ϕ is bounded can be relaxed when the tails of the
distribution of Y are sufficiently light (e.g., sub-Gaussian). Indeed, this is what happens for the special case
of phase retrieval, which is considered immediately below.
For the special case of phase retrieval, a more explicit error metric is given by the matrix minimum
mean square error, defined as


2

1
∗
∗
MMSEPR (δn ) = 2 E XX − E XX | Y , {Ai }1≤i≤n
.
(18)
d
F
Indeed, the vector x can be recovered only up to a sign change, since we observe a function of the scalar
products |hx, ai i|. Clearly, MMSE(δn ) ∈ [0, 1] and MMSE(δn ) = 1 implies that the optimal estimator
coincides with the trivial estimator that outputs the all-0 vector.
The corollary below provides the exact value of δ` for the case of phase retrieval and it is proved in
Appendix A.
Corollary 1 (Information-Theoretic Lower Bound for Phase Retrieval). Let x, {ai }1≤i≤n , and y be distributed according to (6), (7), and (9), respectively. Let n/d → δ. Then, for any δ < 1, we have that
lim lim MMSEPR (δn ) = 1.

σ→0 n→∞

7

(19)

2.3

Upper Bound via Spectral Method

The main result of this section establishes the following: there is a critical value δu such that, for any δ > δu ,
the principal eigenvector of a suitably constructed matrix, call it Dn , provides an estimate x̂ that satisfies
(3).
The threshold δu is defined as
δu = Z
R

1
,

2
EG p(y | |G|)(|G|2 − 1)
dy
EG {p(y | |G|)}

(20)

with G ∼ CN(0, 1). Given the measurements {yi }1≤i≤n , we construct the matrix Dn as
n

Dn =

1X
T (yi )ai a∗i ,
n

(21)

i=1

where T : R → R is a pre-processing function.
At this point we are ready to state our main result, which is proved in Section 5.
Theorem 2 (Spectral Upper Bound for Complex General Sensing Model). Let x, {ai }1≤i≤n , and y be
distributed according to (6), (7), and (8), respectively. Let n/d → δ and define δu as in (20). Let x̂ be the
principal eigenvector of the matrix Dn defined in (21). For any δ > δu , set the pre-processing function T
to the function Tδ∗ given by
√
δ · T ∗ (y)
∗
√
√u √
Tδ (y) =
,
(22)
δ − ( δ − δu )T ∗ (y)
where
T ∗ (y) = 1 −

EG {p(y | |G|)}
.
EG {p(y | |G|) · |G|2 }

(23)

Then, we have that, almost surely,
|hx̂, xi|
> ,
n→∞ kx̂k kxk
2
2
lim

(24)

for some  > 0.
The corollary below provides the exact value of δu and an explicit expression for Tδ∗ (y) for the case of
phase retrieval. Its proof is contained in Appendix B. Note that, for phase retrieval, δu = δ` = 1, i.e., the
spectral upper bound matches the information-theoretic lower bound.
Corollary 2 (Spectral Upper Bound for Phase Retrieval). Let x, {ai }1≤i≤n , and y be distributed according
to (6), (7), and (9), respectively. Let n/d → δ. Let x̂ be the principal eigenvector of the matrix Dn defined in
(21). For any δ > 1, set the pre-processing function T to the function Tδ∗ given by (with y+ ≡ max(0, y)):
Tδ∗ (y) =

y+ − 1
√
.
y+ + δ − 1

(25)

lim lim

|hx̂, xi|
> ,
kx̂k2 kxk2

(26)

Then, we have that, almost surely,
σ→0 n→∞

for some  > 0.
8

Notice that this statement is stronger than the claim that δu (σ 2 ) → 1 as σ 2 → 0, where δu (σ 2 ) is the
spectral threshold at noise level σ 2 . Indeed it requires proving that the scalar product |hx̂, xi| stays bounded
away from 0, as σ 2 → 0. Furthermore, this is achieved with pre-processing function (25) that does not
require to estimate σ, which can be challenging with real data.

3

Main Results: Real Case

d
Let us now briefly discuss what happens
√ in the real case. Let x ∈ R be chosen uniformly at random on the
d-dimensional real sphere with radius d, i.e.,
√
x ∼ Unif( dSd−1
(27)
R ).

Let the sensing vectors {ai }1≤i≤n , with ai ∈ Rd being independent and identically distributed according to
a normal distribution with zero mean and variance 1/d, i.e.,
{ai }1≤i≤n ∼i.i.d. N(0d , Id /d).

(28)

Given gi = hx, ai i, the vector of measurements y ∈ Rn is obtained by drawing each component independently according to the following distribution:
yi ∼ p(y | gi ),

i ∈ [n].

(29)

We can define the “real” phase retrieval model, whereby the measurements are given by the squared scalar
product corrupted by additive Gaussian noise with variance σ 2 , i.e.,


1
(y − gi2 )2
pPR (y | gi ) = √ exp −
.
(30)
2σ 2
σ 2π
We first present the information-theoretic lower bound. Consider the function f : [−1, 1] → R, given
by

Z
f (m) =
R

EG1 ,G2 {p(y | G1 )p(y | G2 )}
dy,
EG {p(y | G)}

with


G ∼ N(0, 1),

(G1 , G2 ) ∼ N 02 ,



1 m
m 1

(31)


.

(32)

Furthermore, set

1
log(1 − m2 ).
(33)
2
Again, Fδ (0) = 0 for any δ > 0. We define the information-theoretic threshold δ` as the largest value of δ
such that the maximum of Fδ (m) is attained at m = 0, i.e.,
Fδ (m) = δ log f (m) +

δ` = sup{δ | Fδ (m) < 0 for m ∈ [−1, 1] \ {0}}.

(34)

As for the error metric, we observe the vector of n measurements y and, given a new sensing vector
an+1 , we want to estimate some function φ(hx, an+1 i) given by
Z
φ(hx, an+1 i) =
ϕ(y)p(y | hx, an+1 i)dy.
(35)
R

9

Then, the minimum mean square error is defined as

2 

MMSE(δn ) = E φ(hX, An+1 i) − E φ(hX, An+1 i) Y , {Ai }1≤i≤n
.

(36)

Recall that, if we do not have access to the vector of measurements y, the trivial estimator E {φ(hX, An+1 i)}
has a mean square error given by

2 


E φ(hX, An+1 i) − E φ(hX, An+1 i)
= Var φ(hX, An+1 i) .
(37)
At this point we are ready to state the information-theoretic lower bound, which is proved in Section
4.2.
Theorem 3 (Information-Theoretic Lower Bound for Real General Sensing Model). Let x, {ai }1≤i≤n+1 ,
and y be distributed according to (27), (28), and (29), respectively. Let n/d → δ and define δ` as in (34).
Furthermore, assume that the function ϕ that appears in (35) is bounded. Then, for any δ < δ` , we have
that

lim MMSE(δn ) = Var φ(hX, An+1 i) .
(38)
n→∞

Remark 1 (Information-Theoretic Lower Bound for Real Phase Retrieval). For the special case of phase
retrieval, a more explicit error metric is given by the matrix minimum mean square error, defined as


2

1
MMSEPR (δn ) = 2 E XX T − E XX T | Y , {Ai }1≤i≤n
.
(39)
d
F
By calculations similar to those in Lemma 7 contained in Appendix A, one can prove that, if the distribution
p(· | G) appearing in (31) is given by (30), then
lim δ` (σ 2 ) = 1/2.

(40)

σ→0

Consequently, by following a proof analogous to that of Corollary 1 in Appendix A, we conclude that, for
any δ < 1/2,
lim lim MMSEPR (δn ) = 1.
(41)
σ→0 n→∞

Let us now move to the spectral upper bound. The threshold δu is defined as
1

δu = Z

EG p(y | G)(G2 − 1)
EG {p(y | G)}


R

,

2

(42)

dy

with G ∼ N(0, 1). Given the measurements {yi }1≤i≤n , we construct the matrix Dn as
n

Dn =

1X
T (yi )ai aT
i ,
n
i=1

where T : R → R is a pre-processing function.
The proof of the following spectral upper bound is discussed in Remark 6 at the end of Section 5.

10

(43)

Theorem 4 (Spectral Upper Bound for Real General Sensing Model). Let x, {ai }1≤i≤n , and y be distributed according to (27), (28), and (29), respectively. Let n/d → δ and define δu as in (42). Let x̂ be the
principal eigenvector of the matrix Dn defined in (43). For any δ > δu , set the pre-processing function T
to the function Tδ∗ given by
√
δ · T ∗ (y)
∗
√u √
Tδ (y) = √
,
(44)
δ − ( δ − δu )T ∗ (y)
where
T ∗ (y) = 1 −

EG {p(y | G)}
.
EG {p(y | G) · G2 }

(45)

Then, we have that, almost surely,
lim

n→∞

|hx̂, xi|
> ,
kx̂k2 kxk2

(46)

for some  > 0.
Remark 2 (Spectral Upper Bound for Real Phase Retrieval). By calculations similar to those in Lemma 8
contained in Appendix B, one can prove that, if the distribution p(· | G) appearing in (31) is given by (30),
then
lim δu (σ 2 ) = 1/2.
(47)
σ→0

Furthermore, by following a proof analogous to that of Corollary 2 in Appendix B, one can prove the
following result. For any δ > 1/2, set the pre-processing function T to the function Tδ∗ given by (with
y+ = max(y, 0))
y+ − 1
√
Tδ∗ (y) =
.
(48)
y+ + 2δ − 1
Then, we have that, almost surely,
|hx̂, xi|
> ,
σ→0 n→∞ kx̂k2 kxk2
lim lim

(49)

for some  > 0. Note that, for real phase retrieval, the spectral upper bound matches the informationtheoretic lower bound.

4
4.1

Proof of Theorems 1 and 3: Information-Theoretic Lower Bound
Complex Case

The crucial point of the proof consists in the computation of the conditional entropy H(Y | A), which
is contained in Lemma 1. Then, we use this result to compute the mutual information for the considered
model. Finally, we provide the proof of Theorem 1.
√ d−1
Lemma 1 (Conditional Entropy). Let x ∼ Unif( dSC
), A = (a1 , . . . , an ) with {ai }1≤i≤n ∼i.i.d.
CN(0d , Id /d), and y = (y1 , . . . , yn ) with yi ∼ p(· | |gi |) and gi = hx, ai i. Let n/d → δ and define δ` as
in (13). Then, for any δ < δ` , we have that
lim

n→∞

1
H(Y | A) = H(Y1 ).
n

11

(50)

Proof. We divide the proof into two steps. The first step consists in showing that
n
o
Z

2
E
(p(y
|
A))
A
1
1
− 
dy − 1 ≤ H(Y | A) − H(Y1 ) ≤ 0,
n
EA {p(y | A)}
n
Rd

(51)

which holds for all n ∈ N and for all δ > 0. The proof of (51) does not require any assumption on the
distribution of x and on the distribution of {ai }1≤i≤n (as long as the vectors {ai }1≤i≤n are independent).
The second step consists in showing that
n
o
Z

2
E
(p(y
|
A))
A
1
dy − 1 = 0.
(52)
lim
n→+∞ n
EA {p(y | A)}
Rd
It is clear that (51) and (52) imply the thesis.
First step. By definition of conditional entropy, we have that
Z
1
1
H(Y | A) =
EA {−p(y | A) log p(y | A)} dy.
n
n

(53)

By using the definition of yi and the fact that they are independent, we can rewrite EA {p(y | A)} as follows:
EA {p(y | A)} = EA,X {p(y | A, X)}
( n
)
Y
= EA,X
p(yi | |hX, Ai i|)

(54)

i=1

=

n
Y

EGi {p(yi | |Gi |)} ,

i=1

where we set Gi = hX, Ai i.
Let us now give an upper bound on the RHS of (53):
1
n

Z

(a)

1
EA {−p(y | A) log p(y | A)} dy ≤
n
Rd

Z

1
=
n

Z

(b)

Rd

−EA {p(y | A)} log EA {p(y | A)} dy
−

Rd
n

1X
=
n

n
Y

EGi {p(yi | |Gi |)}

i=1

n
X

log EGj {p(yj | Gj )} dy

j=1

Z
−EGi {p(yi | |Gi |)} log EGi {p(yi | |Gi |)} dyi
R

i=1

= H(Y1 ),
where in (a) we apply Jensen’s inequality as the function g(x) = −x log x is concave, and in (b) we use
(54). This immediately implies that
1
H(Y | A) − H(Y1 ) ≤ 0.
n
12

(55)

Note that the upper bound (55) is based on the inequality
EA {−p(y | A) log p(y | A)} − (−EA {p(y | A)} log EA {p(y | A)}) ≤ 0.
Let us now find a lower bound to this quantity:
EA {−p(y | A) log p(y | A)} − (−EA {p(y | A)} log EA {p(y | A)})

= EA −p(y | A) log

p(y | A)
EA {p(y | A)}



(a)

= EA {p(y | A)} EZ {−Z log Z}

(b)

= EA {p(y | A)} EZ {−Z log Z + Z − 1}

≥ −EA {p(y | A)} EZ (Z − 1)2


(d)
= −EA {p(y | A)} EZ Z 2 − 1
n
o


EA (p(y | A))2
= −
− EA {p(y | A)} ,
EA {p(y | A)}
(c)

(56)

where in (a) we set Z = p(y | A)/EA {p(y | A)}, in (b) we use that EZ {Z} = 1, in (c) we use that
−z log z + z − 1 ≥ −(z − 1)2 for any z ≥ 0, and in (d) we use again that EZ {Z} = 1. Therefore,
Z
1
1
H(Y | A) − H(Y1 ) =
(EA {−p(y | A) log p(y | A)} − (−EA {p(y | A)} log EA {p(y | A)})) dy
n
n
o

Z  n
EA (p(y | A))2
(a) 1

≥−
− EA {p(y | A)} dy,
n Rd
EA {p(y | A)}
n
o
Z

2
E
(p(y
|
A))
A
(b) 1
dy − 1 ,
=− 
n
E
{p(y
|
A)}
A
d
R
where in (a) we use (56) and in (b) we use that the integral of p(y | A) is 1. This concludes the proof of
(51).
√
Second step. As X ∼ Unif( dSd−1
C ) and Ai ∼ CN(0d , Id /d), we have that
{Gi }1≤i≤n ∼i.i.d. CN (0, 1) .
n
o
Let us rewrite the quantity EA (p(y | A))2 as follows:

( n
)!2 


n
o
Y
2
EA (p(y | A)) = EA
E
p(yi | |hX, Ai i|)
 X

i=1
(
( n
))
Y
(a)
= EA EX1 ,X2
p(yi | |hX1 , Ai i|) · p(yi | |hX2 , Ai i|)
i=1

(b)

= EC

( n
Y

)
EGi,1 ,Gi,2 {p(yi | |Gi,1 |) · p(yi | |Gi,2 |)} ,

i=1

13

(57)

where in (a) X1 and X2 are independent, and in (b) we set Gi,1 = hX1 , Ai i, Gi,2 = hX2 , Ai i, and
C=

hX1 , X2 i
.
kX1 k2 kX2 k2

√ d−1
Then, given C = c, as X1 , X2 ∼i.i.d. Unif( dSC
) and Ai ∼ CN(0d , Id /d), we have that
 

1 c
{(Gi,1 , Gi,2 )}1≤i≤n ∼i.i.d. CN 02 ,
.
c∗ 1
Hence,
1
n

Z
Rd

n
o
EA (p(y | A))2
EA {p(y | A)}

( n
)
Y EGi,1 ,Gi,2 {p(y | |Gi,1 |)p(y | |Gi,2 |)}
dy
EC
EGi {p(y | |Gi |)}
Rd
i=1
( n Z
)
Y
EGi,1 ,Gi,2 {p(y | |Gi,1 |)p(y | |Gi,2 |)}
1
= EC
dyi
n
EGi {p(y | |Gi |)}
R

1
dy =
n
(a)

Z

i=1

1
= EM {(f (M ))n }
n
Z 1
(c) d − 1
=
(f (m))n (1 − m)d−2 dm,
n
0

(b)

where in (a) we use (54) and (57), in (b) we use the fact that f depends only on m = |c|2 , which is
clear from the explicit expression provided by Lemma 6 contained in Appendix A, and in (c) we use that
M ∼ Beta(1, d − 1) by Lemma 9 contained in Appendix C.
Set d0 = d − 2 and δn0 = n/d0 . Thus,
Z 1
Z 1

n
d−2
(f (m)) (1 − m)
dm =
exp n · Fδn0 (m) dm,
(58)
0

0

where Fδn0 (m) is given by (12). Define
F̃δ (m) = δ max(log f (m), 0) + log(1 − m).

(59)

As δ < δ` and n/d0 → δ, there exists δ∗ ∈ (δ, δ` ) such that δn0 < δ∗ for n sufficiently large. As Fδn0 (m) ≤
F̃δn0 (m) and F̃δ (m) is non-decreasing in δ, we have that
Z 1
Z 1



exp n · Fδn0 (m) dm ≤
exp n · F̃δ∗ (m) dm.
(60)
0

0

Note that F̃δ∗ (m) < 0 if and only if Fδ∗ (m) < 0. Thus, by definition of δ` , we have that F̃δ∗ (m) < 0 for
m ∈ (0, 1] when n is sufficiently large. Furthermore, F̃δ∗ (0) = 0 and F̃δ∗ is a continuous function. As a
result, by Lemma 11, the integral in (60) tends to 0 as n → ∞ and the claim immediately follows.
Remark 3 (Mutual Information). An immediate consequence of Lemma 1 is that one can compute the
mutual information I(X; Y , A) for any δ < δ` :
lim

n→+∞

1
I(X; Y , A) = H (EG {p(· | |G|)}) − EG {H(p(· | |G|))} ,
n

where G ∼ CN(0, 1).
14

(61)

Proof of Theorem 1. Define y1:n = (y1 , . . . , yn ) and a1:n = (a1 , . . . , an ). We divide the proof into two
steps. The first step consists in showing that the mutual information between the next observation yn+1 and
the previous observations y1:n tends to 0. More formally, we will prove that
I(Yn+1 ; Y1:n , A1:n | An+1 ) = on (1).

(62)

The second step consists in showing that the estimate obtained on φ(|hx, an+1 i|) given the observations
y1:n is similar to the estimate on φ(|hx, an+1 i|) when no observation is available. This means that the
observations y1:n do not provide any help. More formally, we will prove that

2 


= on (1),
(63)
EY1:n ,A1:n+1 E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|) Y1:n , A1:n
where φ is defined in (14).
Furthermore, we have that

2 

E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|) Y1:n , A1:n
 
2

− E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|) Y1:n , A1:n

2 

= E φ(|hX, An+1 i|) − E φ(|hX, An+1 i|)

= Var φ(|hX, An+1 i|) .

(64)

By applying (63) and (64), the proof of Theorem 1 follows.
First step. By using the chain rule of entropy and that yi is independent from ai+1:n+1 , we obtain that
n+1

1
1 X
H(Y1:n+1 | A1:n+1 ) =
H(Yi | Y1:i−1 , A1:n+1 )
n+1
n+1
=

1
n+1

i=1
n+1
X

H(Yi | Y1:i−1 , A1:i ).

i=1

The sequence sn = H(Yn | Y1:n−1 , A1:n ) is decreasing, as conditioning reduces entropy. Hence sn has a
limit, and this limit must be equal to H(Y1 ) by Lemma 1. Since the Yi are i.i.d., we obtain that
H(Yn+1 | Y1:n , A1:n+1 ) = H(Yn+1 ) + on (1).
By using again that conditioning reduces entropy, we also obtain that
H(Yn+1 | An+1 ) = H(Yn+1 ) + on (1).
By putting these last two equations together, we deduce that (62) holds.

15

Second step. Given two probability distributions p and q, let DKL (p||q) and kp − qkTV denote their
Kullback-Leibler divergence and their total variation distance, respectively. Then,
I(Yn+1 ;Y1:n , A1:n | An+1 )
= EY1:n ,A1:n+1 {DKL (p(yn+1 | Y1:n , A1:n+1 )||p(yn+1 | An+1 ))}
n
o
(a) 1
≥ · EY1:n ,A1:n+1 (kp(yn+1 | Y1:n , A1:n+1 ) − p(yn+1 | An+1 ))kTV )2
2
(Z
(b) 1
· EY1:n ,A1:n+1
p(yn+1 | Y1:n , A1:n+1 )ϕ(yn+1 ) dyn+1
≥
2K 2
R
2 )
Z
−
p(yn+1 | An+1 )ϕ(yn+1 ) dyn+1
R

(c)

=

1
· EY1:n ,A1:n+1
2K 2

(Z

Z

Cd

Z
−
=

R

2 )

Z
p(yn+1 | x, An+1 )ϕ(yn+1 ) dyn+1 dx

p(x)
Cd

(d)

p(yn+1 | x, Y1:n , A1:n+1 )ϕ(yn+1 ) dyn+1 dx

p(x | Y1:n , A1:n )

R

n
o
1
2
·
E
(E
{φ(|hX,
A
i|)}
−
E
{φ(|hX,
A
i|)
|
Y
,
A
})
.
n+1
n+1
1:n
1:n
Y
,A
1:n
1:n+1
2K 2
(65)

where in (a) we use Pinsker’s inequality, in (b) we use that ϕ is bounded and we set kϕk∞ = K, in (c) we
use that X and An+1 are independent, and in (d) we use the definition (14). By combining (62) and (65),
(63) immediately follows.

4.2

Real Case

The proof is very similar to the one provided in Section 4.1 for the complex case. In particular, the crucial
point consists in showing that
1
lim H(Y | A) = H(Y1 ),
(66)
n→∞ n
√
where x ∼ Unif( dSd−1
R ), a = (a1 , . . . , an ) with {ai }1≤i≤n ∼i.i.d. N(0d , Id /d), and y = (y1 , . . . , yn )
with yi ∼ p(· | gi ) and gi = hx, ai i. Then, the proof of Theorem 3 follows similar passages as the proof of
Theorem 1.
In order to prove (66), we show that (51) and (52) hold. The proof of (51) follows the same passages
as the first step of the proof of Lemma 1, hence it is omitted. The proof of (52) is slightly different and we
detail what changes in the remaining part of this section.
Similarly to (54), we have that
EA {p(y | A)} =

n
Y

EGi {p(yi | Gi )} ,

i=1

where Gi = hX, Ai i ∼ N(0, 1). Furthermore, similarly to (57), we also have that
( n
)
n
o
Y
EA (p(y | A))2 = EM
EGi,1 ,Gi,2 {p(yi | Gi,1 ) · p(yi | Gi,2 )} ,
i=1

16

where Gi,1 = hX1 , Ai i, Gi,2 = hX2 , Ai i, and we define
M=

hX1 , X2 i
.
kX1 k2 kX2 k2

√ d−1
) and Ai ∼ N(0d , Id /d), we have that
Then, given M = m, as X1 , X2 ∼i.i.d. Unif( dSR
 

1 m
{(Gi,1 , Gi,2 )}1≤i≤n ∼i.i.d. N 02 ,
.
m 1
Hence,
1
n

Z
Rd

n
o
EA (p(y | A))2
EA {p(y | A)}

1
EM {(f (M ))n }
n
Z 1
Γ( d2 )
d−3
(b) 1
= √
(f (m))n (1 − m2 ) 2 dm,
n πΓ( d−1
2 ) −1
(a)

dy =

(67)

where in (a) we use the definition (33) of f and in (b) we plug in the distribution of M obtained from Lemma
10 contained in Appendix C. Note that
lim

d→∞ d
2

Γ( d2 )
· Γ( d−1
2 )

= 1.

Therefore, by showing that the integral in the RHS of (67) tends to 0, the claim immediately follows.
Set d0 = d − 3 and δn0 = n/d0 . Thus,
Z

1
n

2

(f (m)) (1 − m )

d−3
2

Z

1

dm =
−1

−1


exp n · Fδn0 (m) dm,

(68)

where Fδn0 (m) is defined in (33). Define
F̃δ (m) = δ max(log f (m), 0) +

1
log(1 − m2 ).
2

(69)

As δ < δ` and n/d0 → δ, there exists δ∗ ∈ (δ, δ` ) such that δn0 < δ∗ for n sufficiently large. As Fδn0 (m) ≤
F̃δn0 (m) and F̃δ (m) is non-decreasing in δ, we have that
Z
0

1


exp n · Fδn0 (m) dm ≤

Z

1



exp n · F̃δ∗ (m) dm.

(70)

0

Note that F̃δ∗ (m) < 0 if and only if Fδ∗ (m) < 0. Thus, by definition of δ` , we have that F̃δ∗ (m) < 0 for
m 6= 0 when n is sufficiently large. Furthermore, F̃δ∗ (0) = 0 and F̃δ∗ is a continuous function. As a result,
by Lemma 11, the integral in (70) tends to 0 as n → ∞ and the claim immediately follows.

17

5

Proof of Theorems 2 and 4: Spectral Upper Bound

We will consider the complex case. The proof for the real case is essentially the same and it is briefly
discussed in Remark 6 at the end of this section.
A crucial ingredient of the proof consists in Lemma 2, which is a generalization of Theorem 1 of [LL17].
Before stating this result, we need some definitions. Let G ∼ CN(0, 1), Y ∼ p(· | |G|), and Z = T (Y ).
Assume that Z has bounded support and let τ be the supremum of this support, i.e.,
τ = inf{z : P(Z ≤ z) = 1}.

(71)

For λ ∈ (τ, ∞) and δ ∈ (0, ∞), define

φ(λ) = λ · E
and


ψδ (λ) = λ

Z · |G|2
λ−Z

1
+E
δ





Z
λ−Z

,

(72)


.

(73)

Note that φ(λ) is a monotone non-increasing function and that ψδ (λ) is a convex function. Let λ̄δ be the
point at which ψδ attains its minimum, i.e.,
λ̄δ = arg min ψδ (λ).

(74)

ζδ (λ) = ψδ (max(λ, λ̄δ )).

(75)

λ≥τ

For λ ∈ (τ, ∞), define also
Lemma 2 (Generalization of Theorem 1 of [LL17]). Let x ∼ Unif(Sd−1
C ), {ai }1≤i≤n ∼i.i.d. CN(0d , Id ),
and y be distributed according to (8). Let n/d → δ, G ∼ CN(0, 1) and define Z = T (Y ) for Y ∼ p( · | |G|).
Assume that Z satisfies P(Z = 0) < 1 and that it has bounded support. Let τ be defined in (71). Assume
further that, as λ approaches τ from the right, we have




Z
Z · |G|2
lim E
= lim E
= ∞.
(76)
(λ − Z)2
λ−Z
λ→τ +
λ→τ +
Let x̂ be the principal eigenvector of the matrix Dn , defined as in (21). Then, the following results hold:
(1) The equation
ζδ (λ) = φ(λ)

(77)

admits a unique solution, call it λ∗δ , for λ > τ .
(2) As n → ∞,
|hx̂, xi|2
kx̂k22

kxk22

a.s.

−→



0,






if ψδ0 (λ∗δ ) ≤ 0,

ψδ0 (λ∗δ )
,
ψδ0 (λ∗δ ) − φ0 (λ∗δ )

where ψδ0 and φ0 denote the derivatives of these two functions.

18

if ψδ0 (λ∗δ ) > 0,

(78)

n
n
denote the two largest eigenvalues of Dn . Then, as n → ∞,
≥ λD
(3) Let λD
2
1

a.s.

n
−→ ζδ (λ∗δ ),
λD
1

a.s.

n
−→ ζδ (λ̄δ ).
λD
2

(79)

Before proceeding with the proof, we discuss these results in more detail and we describe in what sense
Lemma 2 provides a generalization of Theorem 1 of [LL17].
Remark 4 (Two different regimes). The results of Lemma 2 imply that, according to the value of δ, we can
distinguish between two possible regimes.
On the one hand, suppose that φ(λ̄δ ) > ψδ (λ̄δ ). Recall that φ(λ) is non-increasing and that λ̄δ is
the point in which ψδ (λ) attains its minimum. Thus, λ̄δ < λ∗δ , which implies that ψδ0 (λ∗δ ) > 0 and that
ζδ (λ∗δ ) > ζδ (λ̄δ ). This means that the scalar product |hx̂, xi| is bounded away from zero and that there is
a strictly positive gap between the two largest eigenvalues of Dn . In this regime, the spectral method that
outputs x̂ solves the weak recovery problem and (24) holds for some  > 0.
On the other hand, suppose that φ(λ̄δ ) ≤ ψδ (λ̄δ ). Thus, λ̄δ ≥ λ∗δ , which implies that ψδ0 (λ∗δ ) ≤ 0 and
that ζδ (λ∗δ ) = ζδ (λ̄δ ). In words, this means that the scalar product |hx̂, xi| converges to zero and that there
is no strictly positive gap between the two largest eigenvalues of Dn . In this regime, the spectral method
that outputs x̂ does not solve the weak recovery problem.
Remark 5 (Lemma 2 and Theorem 1 of [LL17]). Lemma 2 generalizes Theorem 1 of [LL17] in the following
two regards:
• x and {ai }1≤i≤n are complex vectors, while Theorem 1 of [LL17] considers the real case;
• Z can also be negative, while Theorem 1 of [LL17] assumes that Z ≥ 0.
The first generalization does not require additional work as the whole argument of [LL17] generalizes
in the natural way to the complex case: Gaussian random variables become circularly-symmetric complex Gaussian random variables, transposes of vectors and matrices become conjugate transposes, squares
become modulus squares, and so on.
On the contrary, the second generalization is more challenging, as it requires the result of Lemma 3,
which is stated below and proved in Appendix D.

As a final observation, let us point out that Theorem 1 of [LL17] assumes also that E Z · |G|2 >
E {Z}. A careful check shows that this hypothesis is never used in the proof of that theorem, but it is
required only in the proof of some additional results of [LL17].
Lemma 3 (Generalization of [BY12] to non-PSD matrices). Consider the random matrix
Sn =

1
U Mn U ∗ ,
n

(80)

n
where the entries of U ∈ C(d−1)×n are ∼i.i.d. CN(0, 1), and Mn ∈ Cn×n is independent of U . Let λM
1
denote the largest eigenvalue of Mn . Assume that the empirical spectral measure of the eigenvalues of Mn
almost surely converges weakly to the probability distribution H, where H is the law of the random variable
Z. Let ΓH be the support of H and let τ be the supremum of ΓH . Assume also that, as n → ∞,

a.s.

n
λM
−→ α∗ 6∈ ΓH .
1

19

(81)

n
Let n/d → δ, denote by λS
1 the largest eigenvalue of the matrix (80), and define ψδ as in (73). Then, as
n → ∞,

a.s.

if ψδ0 (α∗ ) > 0,

n
λS
1 −→ ψδ (α∗ ),

(82)

a.s.

0
n
λS
1 −→ min ψδ (λ), if ψδ (α∗ ) ≤ 0.
λ>τ

Proof of Lemma 2. In this proof, we follow closely the approach detailed in Section III of [LL17]. First of
all, let us write the matrix Dn defined in (21) as
Dn =

1
AZA∗ ,
n

(83)

where A = [a1 , . . . , an ], Z is a diagonal matrix with entries zi = T (yi ) for i ∈ [n], the random variables
yi are independent and distributed according to p(· | |gi |), and {gi }1≤i≤n ∼i.i.d. CN(0, 1). As the sensing
vectors {ai }1≤i≤n are drawn from the circularly-symmetric complex normal distribution, we can assume
without loss of generality that x = e1 , where e1 is the first element of the canonical basis of Cd .
Consider a matrix U ∈ C(d−1)×n independent of {gi }1≤i≤n and Z. Let the elements of U be ∼i.i.d.
CN(0, 1). Define
1
Pn = U ZU ∗ ,
(84)
n
and
1
qn = U v,
(85)
n
where v = [z1 g1 , . . . , zn gn ]∗ . Then, (83) can be rewritten as


an qn∗
Dn =
,
(86)
qn Pn
Pn
2
2
where an =
i=1 zi |gi | /n is a scalar that converges almost surely to E(Z · |G| ) as n → ∞, with
G ∼ CN(0, 1).
Next, consider a parametric family of matrices {Pn + µqn qn∗ } and let Ln (µ) denote their largest eigenvalues, i.e.,
Ln (µ) = λ1 (Pn + µqn qn∗ )
n
The idea is to compute the largest eigenvalue of Dn , call it λD
1 , and the scalar product between X̂ and e1
via a fixed-point equation involving Ln (µ).
To do so, we first need an intermediate result holding for any matrix D that can be written in the form


a q∗
D=
,
q P

where a ∈ R, P ∈ C(d−1)×(d−1) is a Hermitian matrix and q ∈ Cd−1 is such that kqk =
6 0. Note that the
matrix Dn defined in (21) fulfills such requirements, since the matrix Pn defined in (84) is Hermitian and
qn defined in (85) is such that kqn k =
6 0 with high probability, as P(Z = 0) < 1.
P ≥ · · · ≥ λP
Let λP
≥
λ
be
the
set of eigenvalues of P , and let w1 , w2 , . . . , wd−1 be a correspond1
2
d−1
ing set of eigenvectors. For λ ∈ (max{λP
6 0}, ∞), define
i : hq, wi i =
∗

−1

R(λ) = q (P − λI)

q=

d−1
X
|hq, wi i|2
i=1

20

λP
i −λ

.

(87)

Note that R(λ) increases monotonically from −∞ to 0. Hence, it admits an inverse, call it R−1 (x), for
x < 0. Then, the maximum eigenvalue L(µ) = λ1 (P + µqq ∗ ) is given by
L(µ) = max(R−1 (−1/µ), λP
1 ).

(88)

The proof of (88) is standard, cf. e.g. Lemma 1 in [LL17]. Note that L(µ) is a non-decreasing function such
that limµ→∞ L(µ) = ∞. Indeed, by construction, R−1 (−1/µ) is strictly increasing and
lim R−1 (−1/µ) = ∞.

µ→∞

Furthermore, L(µ) is convex since it is the maximum of a set of linear functions, as
L(µ) = λ1 (P + µqq ∗ ) = max x∗ (P + µqq ∗ )x.
x:kxk=1

Let µ∗ > 0 be the solution to the fixed-point equation
µ = (L(µ) − a)−1 .

(89)

This solution is unique, since L(µ) is a non-decreasing function with limµ→∞ L(µ) = ∞. Then,

and

∗
λD
1 = L(µ ),

(90)


∂+ L(µ∗ )
∂− L(µ∗ )
,
,
|hx̂, e1 i| ∈
∂− L(µ∗ ) + (1/µ∗ )2 ∂+ L(µ∗ ) + (1/µ∗ )2

(91)

2



where ∂− L(µ∗ ) and ∂+ L(µ∗ ) denote the left and right derivative of L(µ), respectively. In particular, if L(µ)
is differentiable at µ∗ , then
L0 (µ∗ )
|hx̂, e1 i|2 = 0
.
(92)
L (µ) + (1/µ∗ )2
The proof of (90), (91), and (92) uses the characterization (88) and it is analogous to the proof of Proposition
2 in [LL17].
At this point, we need to compute Ln (µ) for the matrix Dn defined in (86). The eigenvalues of a low
rank perturbation of a random matrix are studied in [BGN11]. However, we cannot apply those results, as
Pn and qn are dependent. Hence, we write
Pn + µqn qn∗ =

1
U Mn U ∗ ,
n

where Mn is independent of U with
Mn = Z +

µ ∗
vv .
n

n
n
n be the set of eigenvalues of M
We start by studying the spectrum of Mn . Let λM
≥ λM
≥ · · · ≥ λM
n
n
1
2
and let
n
1 X
f Mn =
δλMn
i
n−1

i=2

be the empirical spectral measure of the last n − 1 eigenvalues.

21

Then, standard interlacing theorems (see [HJ12, Section 4.3]) yield that f Mn almost surely converges
weakly to the probability law of Z. Furthermore, by using the characterization (88), we can show that
a.s.

λ1Mn −→ λµ = Q−1 (1/µ),

(93)

where Q−1 is the inverse of the function

Q(λ) = E

Z 2 · |G|2
λ−Z


.

The proof of these results is the same as the proof of Proposition 3 in [LL17].
Note that Q(λ) is defined for λ ∈ (τ, ∞), it is continuous and strictly decreasing with Q(∞) = 0.
Furthermore, by hypothesis (76), we have that limλ→τ + Q(λ) = ∞. Thus, Q(λ) admits an inverse and
Q−1 (1/µ) is well-defined for all µ > 0.
Let us now consider the matrix n1 U Mn U ∗ . First, if Z ≥ 0, then Mn is positive semi-definite (PSD)
and we can apply results from [BY12] to compute the limit of Ln (µ). If Mn is not necessarily PSD, we use
Lemma 3 with α∗ = λµ to conclude that
a.s.

Ln (µ) −→ ψδ (λµ ),

if ψδ0 (λµ ) > 0,

a.s.

Ln (µ) −→ min ψδ (λ), if ψδ0 (λµ ) ≤ 0.

(94)

λ>τ

The remaining part of the proof follows the argument of Section III-D in [LL17]. For the sake of
readability, we reproduce it below.
We start by proving the first claim of the lemma. For n ≥ 1, let µn be the unique solution to the
fixed-point equation (89). Then,
Ln (µn ) − 1/µn = an .
Now, fix any µ > 0. Then, by using the definition (75) and the fact that λµ = Q−1 (1/µ), (94) immediately
implies that, as n → ∞,
a.s.
Ln (µ) − 1/µ −→ ζδ (Q−1 (1/µ)) − 1/µ.
(95)
a.s.

Note that, as n → ∞, an −→ E(Z · |G|2 ). Furthermore, as Ln (µ) and ζδ (µ) are non-decreasing, the two
functions on both sides of (95) are strictly increasing. Consequently, by Lemma 3 in Appendix E of [LL17],
we conclude that
a.s.
µn −→ µ∗ ,
(96)
where µ∗ is the unique fixed point such that
ζδ (Q−1 (1/µ∗ )) = E(Z · |G|2 ) + 1/µ∗ .

(97)

λ∗ = Q−1 (1/µ∗ ).

(98)

ζδ (λ∗ ) = E(Z · |G|2 ) + Q(λ∗ ) = φ(λ∗ ),

(99)

Define
Then, (97) can be rewritten as

where φ is defined in (72). By construction, ζδ (λ) is a non-decreasing continuous function on (τ, ∞)
and φ(λ) is a strictly decreasing continuous function. Furthermore, by hypothesis (76), we have that
22

limλ→τ + φ(λ) = ∞. Hence, the existence and the uniqueness of λ∗ satisfying (99) is guaranteed. This
suffices to prove the first claim of the lemma.
Let us now move on to the proof of the second claim of the lemma. Suppose that ζδ (Q−1 (1/µ)) is
differentiable at µ = µ∗ . Then, as Ln (µ) is convex for any n ≥ 1, by Lemma 4 in Appendix E of [LL17],
we have that
−1
−ζδ0 (Q−1 (1/µ∗ ))
a.s. dζδ (Q (1/µ))
∂− Ln (µn ) −→
= 0 −1
.
dµ
Q (Q (1/µ∗ )) · (µ∗ )2
µ=µ∗
Similarly,
a.s.

∂+ Ln (µn ) −→

−ζδ0 (Q−1 (1/µ∗ ))
.
Q0 (Q−1 (1/µ∗ )) · (µ∗ )2

By using (91), we obtain that
a.s.

|hx̂, e1 i|2 −→

ζδ0 (Q−1 (1/µ∗ ))
ζδ0 (λ∗ )
=
,
ζδ0 (Q−1 (1/µ∗ )) − Q0 (Q−1 (1/µ∗ ))
ζδ0 (λ∗ ) − φ0 (λ∗ )

where the equality follows from the definition (98) of λ∗ and from the fact that Q0 (λ) = φ0 (λ). In order
to prove the second claim of the lemma, it suffices to note that, by its definition in (75), ζδ0 (λ) = ψδ0 (λ) if
ψδ0 (λ) > 0, and ζδ0 (λ) = 0 if ψδ0 (λ) < 0.
n
Finally, let us prove the third claim of the lemma. By using (90), we immediately obtain that λD
=
1
Ln (µn ). By applying (96) and Lemma 3 in Appendix E of [LL17], we conclude that
a.s.

n
λD
−→ ζδ (λ∗ ).
1

As Pn is obtained by deleting the first row and column of Dn , by applying Cauchy interlacing theorem (see,
e.g., [HJ12, Theorem 4.3.17]), we also have that
Dn
n
n
λP
≤ λP
2 ≤ λ2
1 .

Furthermore, the upper edge of the support of the limiting spectral distribution of Pn is given by [SC95,
Section 4] and [BY12, Lemma 3.1]
min ψδ (λ) = ζδ (λ̄δ ),
λ>τ

where λ̄δ is defined in (74). Therefore,

a.s.

n
λD
−→ ζδ (λ̄δ ),
2

which concludes the proof.
At this point, we are ready to prove our spectral upper bound.
Proof of Theorem 2. Note that the normalization of x and {ai }1≤i≤n required in Lemma 2 is different from
the normalization required in Theorem 2. However, the scalar product hx, ai i is the same and the data
matrix Dn changes by a factor d. Hence, the principal eigenvector x̂ is not affected by this change in the
normalization.
Let G ∼ CN(0, 1), Y ∼ p(· | |G|) and Z = T (Y ), where p is defined in (8) and T is some preprocessing function that we will choose later on. We will that assume the supremum τ of the support of Z
is strictly positive and that conditions (76) are satified, and will verify later that our choice of the function
T satisfies these requirements. Recall that the function ψδ (λ) defined in (73) is convex and that it attains
23

its minimum at the point λ̄δ . Since by condition (76) ψδ (λ) ↑ ∞ as λ ↓ 0, we have λ̄δ ∈ (τ, ∞). Hence,
ψδ0 (λ̄δ ) = 0. By calculating the derivative of ψδ (λ) and setting it to 0, we have

E

Z2
(λ̄δ − Z)2



1
= .
δ

(100)

Furthermore, as pointed out in Remark 4, (24) holds for some  > 0 if and only if
φ(λ̄δ ) > ψδ (λ̄δ ).

(101)

As τ > 0, we also have that λ̄δ > 0. Consider now the matrix Dn0 = Dn /α for some α > 0. Then, the
principal eigenvector of Dn0 is equal to the principal eigenvector of Dn . Hence, we can assume without loss
of generality that λ̄δ = 1. Consequently, the conditions (100) and (101) can be respectively rewritten as


1
Z2
= ,
(102)
E
2
(1 − Z)
δ


1
Z(|G|2 − 1)
> .
(103)
E
1−Z
δ
Furthermore, as Z = T (Y ), we also obtain that
 Z 
2
Z2
T (y)
E
=
EG {p(y | |G|)} dy,
(1 − Z)2
R 1 − T (y)

 Z

Z(|G|2 − 1)
T (y)
E
=
EG p(y | |G|) · (|G|2 − 1) dy.
1−Z
R 1 − T (y)


(104)

Let T ∗ (y) be defined in (23). Note that, if we substitute T (y) = T ∗ (y) into the RHS of (104), then




Z2
Z(|G|2 − 1)
1
E
=E
= ,
2
(1 − Z)
1−Z
δu
where δu is defined in (20). Let Tδ∗ (y) be defined in (22). Then,
r
Tδ∗ (y)
δu T ∗ (y)
=
,
1 − Tδ∗ (y)
δ 1 − T ∗ (y)
which immediately implies that

(Tδ∗ (Y ))2
1
E
= ,
∗
2
(1 − Tδ (Y ))
δ
 ∗

Tδ (Y )(|G|2 − 1)
1
1
> .
E
=√
1 − Tδ∗ (Y )
δ
δ · δu


(105)
(106)

As a result, in order to complete the proof, we need to show that the function Tδ∗ (y) fulfills the following
requirements:
(1) Tδ∗ (y) is bounded;
(2) P(Tδ∗ (Y ) = 0) < 1;
24

(3) the supremum τ of the support of Tδ∗ (Y ) is strictly positive;
(4) the condition (76) holds.
Note that Tδ∗ (y) is bounded, as T ∗ (y) ≤ 1. Furthermore, if

EG {p(y | |G|)} = EG p(y | |G|)|G|2 ,

(107)

identically, then δu = ∞ and the claim of Theorem 2 trivially holds. Hence, we can assume that (107)
does not hold, which implies that the function T ∗ is not equal to the constant value 0. Consequently,
P(Tδ∗ (Y ) = 0) < 1.
By definition (23) of T ∗ , we have that

 Z


1
EY
=
EG p(y | |G|) · |G|2 dy = EG |G|2 = 1.
(108)
∗
1 − T (Y )
R
Hence, P(T ∗ (Y ) > 0) > 0, which implies that P(Tδ∗ (Y ) > 0) > 0. Consequently, the supremum τ of the
support of Tδ∗ (Y ) is strictly positive.
If P(Tδ∗ (Y ) = τ ) > 0, then the condition (76) is satisfied. Suppose now that P(Tδ∗ (Y ) = τ ) = 0. Then,
for any 1 > 0, there exists ∆1 (1 ) such that

0 < P Tδ∗ (Y ) ∈ (τ − ∆1 (1 ), τ ) ≤ 1 .
(109)
Define
Tδ∗ (y, 1 )

=

 ∗
 Tδ (y),


if Tδ∗ (y) ≤ τ − ∆1 (1 ),
(110)

τ − ∆1 (1 ),

otherwise.

Clearly, the random variable Tδ∗ (Y, 1 ) has a point mass, hence the condition (76) is satisfied.
As a final step, we show that we can take 1 = 0. Define
n

Dn (1 ) =

1X ∗
Tδ (yi , 1 )ai a∗i .
n
i=1

Define also

n

Dn =

1X ∗
Tδ (yi )ai a∗i .
n
i=1

Let x̂(1 ) and x̂ be the principal eigenvectors of Dn (1 ) and of Dn , respectively. Then,
kDn (1 ) − Dn kop ≤ C1 · ∆1 (1 ),

(111)

where the constant C1 depends only on n/d. By Lemma 2, there is a strictly positive gap, call it θ, between
the first and the second eigenvalue of Dn (1 ). Consequently, by the Davis-Kahan theorem [DK70], we
conclude that
kx̂(1 ) − x̂k2 ≤ C2 · ∆1 (1 ),
(112)
where the constant C2 depends only on n/d and on θ. In words, for any n, as 1 tends to 0, the principal
eigenvector of Dn (1 ) tends to the principal eigenvector of Dn . This means that we can set T = Tδ∗ and
the proof is complete.
25

Remark 6 (Proof of Spectral Upper Bound for the Real Case). First, we need to prove a result analogous
to that of Lemma 2, where x ∼ Unif(Sd−1
R ), {ai }1≤i≤n ∼i.i.d. N(0d , Id ), y is distributed according to
(29), and G ∼ N(0, 1). To do so, one can follow the proof of Theorem 1 of [LL17]. The technical difficulty
consists in the fact that the matrix Mn is not necessarily PSD. In order to solve this issue, we apply the
version of Lemma 3 for the real case discussed in Remark 7 at the end of Appendix D. At this point, the
proof of Theorem 4 follows from the same argument as the proof of Theorem 2.

6

Comparison with Message Passing Algorithms

6.1

Motivation and Background

Message passing algorithms have proved successful in a broad range of statistical estimation problems, including high-dimensional regression [BM12], robust regression [DM16], low-rank matrix estimation [DM13,
KMZ13, MR16, KKM+ 16], and network structure estimation [DKMZ11, MNS14, MX16]. A bold conjecture from statistical physics suggests that – for these and other problems – message passing approaches
achieve optimal statistical performances among polynomial-time algorithms. In view of this conjecture, it
is interesting to compare our spectral approach to message passing algorithms. We will present two types of
results (with δu the spectral threshold defined in (42)):
1. We prove that, for δ < δu (i.e. in the regime in which the spectral approach fails), message passing
converges to an un-informative fixed point, even if initialized in a state that is correlated with the true
signal x.
2. Vice versa, for δ > δu (when the spectral algorithm achieves weak recovery), we consider a linearized
message passing algorithm, and prove that the un-informative fixed point is unstable. The proof of
this fact builds on the analysis contained in the previous pages.
For the sake of simplicity, we will assume that the signal x and the measurement matrix A are real.
Of particular interest for the present setting is approximate message passing (AMP) [DMM09, BM11]:
this is a broad class of iterative methods that operates with dense random matrices (as the sensing matrix
A in the present case). In particular, in [Ran11] it was proposed a “generalized approximate message
passing” (GAMP) scheme, which is an AMP algorithm for Bayesian estimation in non-linear regression
models. This approach was further developed in the context of phase retrieval in [SR15]. We will follow the
same Bayesian formulation here, by considering an AMP algorithm that is equivalent to GAMP although
somewhat simpler.
In order to minimize technical overhead, we assume throughout this section that the conditional density
p(y | g) is bounded and two times differentiable with respect to g. Denote by ∂g p(y | g) and ∂g2 p(y | g) the
first and the second derivative of p(y | g), respectively. Let G ∼ N(0, 1) and define the function
√
EG {∂g p(y | q̄ x + q̄G)}
√
F(x, y; q̄) =
.
(113)
EG {p(y | q̄ x + q̄G)}
We further define the following “state evolution” recursion:
µt+1 = δ · h(qt ) ,
µt
qt =
,
1 + µt
26

(114)

where

(

Z
h(q) =

EG0
R

2 )
√
√
EG1 {∂g p(y | qG0 + 1 − qG1 )}
√
dy ,
√
EG1 {p(y | qG0 + 1 − qG1 )}

(115)

with G0 , G1 ∼i.i.d. N(0, 1).
Given the sensing matrix A = (a1 , . . . , an )T ∈ Rn×d , and the vector of measurements y = (y1 , . . . , yn ) ∈
n
R , the message
passing algorithm updates iteratively the estimate z t ∈ Rd of the signal x ∈ Rd , with
√
kxk2 = d, according to the iteration
z t+1 = AT ft (ẑ t ; y) − bt z t ,

(116)

ẑ t = Az t − ft−1 (ẑ t−1 ; y) .

Here, the function ft (ẑ; y) = (ft (ẑ1 ; y1 ), . . . , ft (ẑn ; yn )) is understood to be applied component-wise to its
arguments and bt it is defined as
ft (ẑ; y) = F(ẑ, y; 1 − qt ) ,
(117)
and the “Onsager coefficient” bt is defined as
bt = δ · E{ft0 (µt G0 +

√

µt G1 ; Y )} ,

(118)

where ft0 (ẑ; y) denotes the derivative of ft (ẑ; y) with respect to ẑ, and the expectation is with respect to
G0 , G1 ∼i.i.d. N(0, 1) and Y ∼ p( · |G0 ). The recursion (116) is initialized with z 0 ∈ Rd and it is understood that f−1 (· ; ·) = 0n .
State evolution precisely tracks the asymptotics of AMP. The next statement is a consequence of [BM11,
JM13]. We refer to Appendix E for its proof.
Lemma √
4 (State Evolution for AMP Iteration (116)). Let x ∈ Rd denote the unknown signal such that
kxk2 = d, A = (a1 , . . . , an )T ∈ Rn×d with {ai }1≤i≤n ∼i.i.d. N(0d , Id /d), and y = (y1 , . . . , yn ) with
yi ∼ p(· | hx, ai i). Consider the AMP iterates z t , ẑ t defined in (116), where ft (ẑ; y) and bt are given by
(117) and (118), respectively. Assume that the initialization z 0 is independent of A and that, almost surely,
1
hx, z 0 i = µ0 ,
n→∞ d
lim

1 0 2
kz k = µ20 + µ0 .
n→∞ d
lim

(119)

Let the state evolution recursion qt , µt be defined as in (114) with initialization µ0 . Then, for any t, and for
any function ψ : R2 → R such that |ψ(u) − ψ(v)| ≤ L(1 + kuk2 + kvk2 )ku − vk2 for some L ∈ R, we
have that, almost surely,
n

1X
√
lim
ψ(xi , zit ) = E {ψ(X0 , µt X0 + µt G)} ,
n→∞ n

(120)

i=1

where the expectation is taken with respect to X0 , G ∼i.i.d. N(0, 1).
Informally, this lemma states that z t is a noisy version of the signal x, namely z t ≈ µt x +
g ∼ N(0d , Id ), and that this approximation holds for empirical averages.

27

√

µt g, with

6.2

Results

In order to obtain a non-vanishing weak recovery threshold, we assume that the observation model satisfies
the condition
EG {∂g p(y | G)} = 0 ,
(121)
where the expectation is with respect to G ∼ N(0, 1). Notice that this implies h(0) = 0, therefore µt =
qt = 0 is a fixed point of state evolution. Furthermore, F(0, y; 1) = 0, therefore qt = 0, z t = 0d is a fixed
point of the message passing algorithm. We will refer to this as to the “un-informative fixed point”. Note
that the condition (121) holds – among others – for the phase retrieval problem.
Vice versa, if EG {∂g p(y | G)} =
6 0, then µ1 > 0 even if µ0 = 0, for any δ > 0. Thanks to Lemma 4,
this implies that weak recovery is possible for all δ > 0. Hence, we will assume that the condition (121)
holds.
The first result of this section establishes the following: for δ < δu , the message passing algorithm fails
even if the initial condition has a positive correlation with the unknown signal. We refer to Appendix E for
its proof.
Theorem
5 (Message Passing Fails for δ < δu ). Let x ∈ Rd denote the unknown signal such that kxk2 =
√
d. Let A = (a1 , . . . , an )T ∈ Rn×d with {ai }1≤i≤n ∼i.i.d. N(0d , Id /d), and y = (y1 , . . . , yn ) with
yi ∼ p(· | hx, ai i). Let n/d → δ and define δu as in (42). Let G ∼ N(0, 1) and assume that the
condition (121) holds for any y ∈ R.
Consider the AMP algorithm defined in (116), and assume that the initial condition z 0 is such that
hz 0 , xi
= .
n→∞ kz 0 k kxk
2
2
lim

(122)

Then, for any δ < δu , there exists 0 (δ) such that for any  ∈ (0, 0 (δ)), almost surely,
1 t
z
t→∞ n→∞ d
lim lim

2

(123)

= 0d .

Next, we consider the case δ > δu and we linearize the iteration (116) around the non-informative fixed
point.
Lemma 5 (Linearized AMP Equations). Consider the one-step map defined in (116) and assume that condition (121) holds. Define Rt = kz t k2 + kẑ t−1 k2 . Then, as Rt → 0 and qt → 0, we obtain
 t+1 


zt
z
+ o(Rt ) + Rt oqt (1) ,
(124)
= Ln
ẑ t−1
ẑ t
where Ln ∈ R(n+d)×(n+d) is defined as

Ln =

AT J A −AT J 2
A
−J


,

(125)

and J ∈ Rn×n is a diagonal matrix with entries ji = F0 (0, yi ; 1) for i ∈ [n], with F0 denoting the derivative
of F with respect to the first argument.
The second result of this section establishes the following: for δ > δu , the un-informative fixed point is
unstable for the iteration (116), i.e., the matrix Ln has an eigenvalue that is larger than 1 in modulus. We
refer to Appendix F for its proof.
28

Theorem 6 (Message Passing Escapes from
for δ > δu ). Let x ∈ Rd denote
√ Un-informative Fixed Point
the unknown signal such that kxk2 = d. Let A = (a1 , . . . , an )T ∈ Rn×d with {ai }1≤i≤n ∼i.i.d.
N(0d , Id /d), and y = (y1 , . . . , yn ) with yi ∼ p(· | hx, ai i). Let n/d → δ and define δu as in (42).
Furthermore, assume that (121) holds for any y ∈ R. Define Ln ∈ R(n+d)×(n+d) as in (125), where
J ∈ Rn×n is a diagonal matrix with entries ji = G(0, yi ; 1) for i ∈ [n]. Then, the eigenvalues of Ln are
n
real and the largest of them, call it λL
1 , is such that, for any δ > δu ,
n
lim λL
1 > 1.

n→∞

7

(126)

Numerical Experiments

We focus on the phase retrieval problem and present some numerical results to illustrate the performance
achieved by the proposed spectral method. First, we consider the case in which the unknown vector is
chosen uniformly at random and the sensing vectors are Gaussian. Then, we consider the more practical
scenario in which the unknown vector is an image and the sensing vectors come from a coded diffraction
model.

7.1

Gaussian Sensing Vectors for Synthetic Data

Let us consider the complex case. In our experiments,
the vector x is chosen uniformly at random on the d√
dimensional complex sphere with radius d, the sensing vectors {ai }1≤i≤n are i.i.d. circularly-symmetric
normal with variance 1/d and, for i ∈ [n], the measurement yi is equal to |hx, ai i|2 . We take d = 4096
and the numerical simulations are averaged over nsample = 40 independent trials. The results are plotted in
Figure 2a.
The red curve corresponds to the proposed pre-processing function given by
T (y) =

y−1
p
.
y + δ̃ − 1

(127)

We pick δ̃ = 1.001 and, as shown by the figure, weak recovery is possible for values of δ very close to 1.
The black curve corresponds to the pre-processing function given by
(
1,
for y > t,
T (y) =
(128)
0,
otherwise.
This choice was proposed in [WGE16] and it is also considered in [LL17], where the authors refer to it as
the “subset algorithm”. For each value of t, we can compute the smallest value of δ, call it δ ∗ (t), that yields
a strictly positive scalar product according to the result of Lemma 2. Hence, we pick t = 2 that corresponds
to the smallest value of δ ∗ (t) over t ∈ {0.25, 0.5, 0.75, . . . , 10}.
The blue curve corresponds to the pre-processing function given by
(
y,
for y ≤ t,
T (y) =
(129)
0,
otherwise.
This choice corresponds to the truncated spectral initialization proposed in [CC17] and it is also considered
in [LL17], where the authors refer to it as the “trimming algorithm”. For each value of t, we can compute the
29

smallest value of δ, call it δ ∗ (t), that yields a strictly positive scalar product according to the result of Lemma
2. Hence, we pick t = 5.25 that corresponds to the smallest value of δ ∗ (t) over t ∈ {0.25, 0.5, 0.75, . . . , 10}.
Note that the numerical simulations follow closely the theoretical prediction given by (78). Furthermore,
the choice of the pre-processing function (127) yields a remarkable performance gain with respect to both
the subset algorithm and the trimming algorithm.
Similar considerations apply to the √
real case. Here, the vector x is chosen uniformly at random on the
d-dimensional real sphere with radius d and the sensing vectors {ai }1≤i≤n are i.i.d. normal with zero
mean and variance 1/d. We pick d = 4096 and nsample = 40. The results are plotted in Figure 2b. Again,
the numerical simulations follow closely the theoretical prediction. The red curve corresponds to the preprocessing function given by (127), where we pick δ̃ = 1.001. Note that weak recovery is possible for
values of δ very close to 1/2. The blue curve corresponds to the pre-processing function given by (128),
where we pick t = 2 which yields the smallest value of δ ∗ (t) over t ∈ {0.25, 0.5, 0.75, . . . , 10}. The black
curve corresponds to the pre-processing function given by (129), where we pick t = 7 which yields the
smallest value of δ ∗ (t) over t ∈ {0.25, 0.5, 0.75, . . . , 10}.

7.2

Coded Diffraction Model for Natural Images

We consider a model of coded diffraction patterns in which the sensing vectors {ar }1≤r≤n are obtained as
follows. For t1 ∈ [d1 ] and t2 ∈ [d2 ], denote by ar (t1 , t2 ) the (t1 , t2 )-th component of the vector ar ∈ Cd ,
with d = d1 · d2 . Then,
ar (t1 , t2 ) = d` (t1 , t2 ) · ei2πk1 t1 /d1 · ei2πk2 t2 /d2 ,
(130)
where i denotes the imaginary unit. The index r ∈ [n] is associated to a pair (`, k), with ` ∈ [L]; the index
k ∈ [d] is associated to a pair (k1 , k2 ) with k1 ∈ [d1 ] and k2 ∈ [d2 ]. As usual, the measurement yr of
an unknown d-dimensional vector x is equal to |hx, ar i|2 . As an immediate consequence, the number of
measurements n is equal to L · d, therefore δ = L ∈ N. In words, for a fixed `, we collect the magnitude
of the diffraction pattern of x modulated by d` . By varying ` and changing the modulation pattern d` , we
generate L distinct views. The vectors {d` }1≤`≤L are i.i.d. and their entries are also i.i.d. drawn uniformly
from the set {1, −1, i, −i}.
We test the spectral method on the digital photograph represented in Figure 1a. Each color image can
be viewed as a d1 × d2 × 3 array. We run the spectral algorithm separately on the vectors xj ∈ Rd , where
j ∈ {1, 2, 3}. In our example, d1 = 820 and d2 = 1280. Let x̂j be the estimate of xj provided by the
spectral method. Then, we employ as a performance metric the average squared normalized scalar product
3

1 X |hx̂j , xj i|2
.
3
kx̂j k2 kxj k2

(131)

j=1

Note that the scalar product between the input and the measurement vectors can be interpreted as a 2dimensional Fourier transform, hence it can be computed with an FFT algorithm. In order to evaluate the
principal eigenvector of the data matrix, we use the power method with a random initialization, as described
in Appendix B of [CLS15b]. As a stopping criterion, we require that one of the following two conditions is
fulfilled: either the number of iterations reaches the maximum value of 10000, or the modulus of the scalar
product between the estimate at the current iteration T and at the iteration T − 10 is larger than 1 − 10−7 .
The results are summarized in Figure 3. The red curve corresponds to the proposed pre-processing function. In this case, the eigenvalues of the data matrix can be negative. Recall that the power method outputs
the eigenvector associated to the largest eigenvalue in modulus, while we are interested in the eigenvector
30

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

1

2

3

4

5

6

7

8

3

3.5

4

(a) Complex case.

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

0.5

1

1.5

2

2.5

(b) Real case.

Figure 2: Performance of the spectral method for the phase retrieval problem where the unknown vector is
uniformly random on the sphere and the sensing vectors are Gaussian. On the x-axis, we have the ratio δ
between the number of samples and the dimension of the signal; on the y-axis, we have the square of the
normalized scalar product between the unknown signal x and the estimate x̂. Note that the proposed choice
of the pre-processing function (red curve) provides a significant performance improvement with respect to
the subset algorithm considered in [WGE16, LL17] (black curve) and the truncated spectral initialization
considered in [CC17, LL17] (blue curve).

31

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

2

4

6

8

10

12

14

16

18

20

Figure 3: Performance of the spectral method for the phase retrieval problem where the unknown vector is
a digital photograph and the sensing vectors are obtained from a coded diffraction model. On the x-axis,
we have the ratio δ between the number of samples and the dimension of the signal; on the y-axis, we have
the square of the normalized scalar product between the unknown signal x and the estimate x̂ (averaged on
the three RGB components of the image). Note that the proposed choice of the pre-processing function (red
curve) provides a significant performance improvement with respect to the truncated spectral initialization
considered in [CC17] (blue curve).
associated to the largest eigenvalue. To address this issue, we add to the data matrix a multiple α of the
identity matrix. However, as α grows, the convergence of the power method becomes slower and slower. In
order to reduce the negative tail of the distribution of eigenvalues and, consequently, the value of α, we pick
the pre-processing function given by
T1 (y) = max(T (y), −M ),

(132)

where T (y) is defined in (127), δ̃ = 1.001, and M = 40. In this way, by taking α = 100, the largest
eigenvalue in modulus has positive sign.
The blue curve corresponds to the truncated spectral initialization in [CC17], i.e., the pre-processing
function is given by (129) with t = 9.
The numerical simulations for the optimal pre-processing function follow closely the theoretical predictions (78) obtained for a Gaussian measurement matrix, with the exception of the point δ = 2. On the
contrary, the numerical simulations for the truncated spectral initialization show a different behavior with
respect to the Gaussian model. Our algorithm provides weak recovery of the original image for δ ≥ 3,
while the truncated spectral initialization requires δ ≥ 6. Furthermore, for any value of δ, the proposed
choice of the pre-processing function yields a better performance than the choice in [CC17]. For a visual
representation of these results, see Figure 1.

32

Acknowledgement
M. M. was supported by an Early Postdoc.Mobility fellowship from the Swiss National Science Foundation.
A. M. was partially supported by grants NSF DMS-1613091 and NSF CCF-1714305.

A

Proof of Corollary 1

We start by providing in Lemma 6 a less compact, but more explicit form of the expression (10). This
more explicit expression is employed to prove Lemma 7, which yields the value of δ` for the case of phase
retrieval. Finally, we provide the proof of Corollary 1.
Lemma 6 (Explicit Formula for f (m) - Complex Case). Consider the function f : [0, 1] → R defined in
(10). Then, f (m) is given by the following expression:

 2

√ 
Z +∞ Z +∞
2r1 r2 m
1
r + r22
· I0
dr1 dr2
4r1 r2 · p(y | r1 )p(y | r2 ) · exp − 1
1−m 0
1−m
1−m
0
dy,
f (m) =
Z +∞

2
2r · p(y | r) · exp −r dr
R

Z

0

(133)

where I0 denotes the modified Bessel function of the first kind, given by
Z
1 π
I0 (x) =
exp (x cos θ) dθ.
π 0

(134)

Proof. Let us rewrite G as
(R)

G=G

(I)

(R)

+ jG ,

with (G




1
, G ) ∼ N 0d , I2 ,
2
(I)

i.e., G(R) and G(I) are i.i.d. Gaussian random variables with mean 0 and variance 1/2. Set
q
2
2
R=
G(R) + G(I) .
√
Then, R follows a Rayleigh distribution with scale parameter 1/ 2, hence
Z +∞

EG {p(y | |G|)} = ER {p(y | R)} =
2r · p(y | r) · exp −r2 dr.
0

Let us rewrite (G1 , G2 ) as
(R)

(G1 , G2 ) = (G1
with



(I)

(R)

+ jG1 , G2


1


1
(R)
(R)
(I)
(I)
 <(c)
(G1 , G2 , G1 , G2 ) ∼ N 
 0d , 2  0
−=(c)

33

(I)

+ jG2 ),

<(c)
1
=(c)
0

0
=(c)
1
<(c)


−=(c)

0
 ,
<(c) 
1

(135)

and consider the following change of variables:
 (R)
G1 = R1 cos θ1



 (R)
G2 = R2 cos θ2
.
(I)

G
=
R
sin
θ

1
1
1

 (I)
G2 = R2 sin θ2
Then, after some algebra, we have that
Z +∞ Z +∞ Z 2π Z 2π
1
r1 r2 · p(y | r1 )p(y | r2 )·
EG1 ,G2 {p(y | |G1 |) · p(y | |G2 |)} = 2
π (1 − |c|2 ) 0
0
0
0
 2

r1 + r22 − 2r1 r2 (<(c) cos(θ2 − θ1 ) − =(c) sin(θ2 − θ1 ))
exp −
dr1 dr2 dθ1 dθ2 .
1 − |c|2
(136)
By writing (<(c), =(c)) = (|c| cos θc , |c| sin θc ) and by using the definition (134), we can further simplify
the RHS of (136) as
 2



Z +∞ Z +∞
1
r1 + r22
2r1 r2 |c|
4r1 r2 · p(y | r1 )p(y | r2 ) · exp −
· I0
dr1 dr2 .
(137)
1 − |c|2 0
1 − |c|2
1 − |c|2
0
From (135) and (137), the claim easily follows.
Lemma 7 (Computation of δ` for Phase Retrieval). Let δ` (σ 2 ) be defined as in (13) and assume that the
distribution p(· | |G|) appearing in (10) is given by (9). Then,
lim δ` (σ 2 ) = 1.

σ→0

(138)

Proof. For the special case of phase retrieval, it is possible to compute explicitly the function f (m) defined
in (10) and simplified in Lemma 6. Indeed,
Z +∞
Z +∞

√
(a)
2r · pPR (y | r) · exp −r2 dr =
pPR (y | z) · exp (−z) dz
0
0


Z +∞
1
(y − z)2
(b)
√ exp −
=
· exp (−z) · H(z) dz (139)
2σ 2
−∞ σ 2π
(c)

= EZ {exp(−Z)H(Z)} ,

where in (a) we do the change of variables z = r2 ; in (b) we use the definition (9) and we define H(x) = 1
if x ≥ 0 and H(x) = 0 otherwise; and in (c) we define Z ∼ N(y, σ 2 ). In the limit σ 2 → 0, we have that
EZ {exp(−Z)H(Z)} → exp(−y) · H(y),

34

(140)

by Lebesgue’s dominated convergence theorem. Similarly,

 2

√ 
Z +∞ Z +∞
2r1 r2 m
r1 + r22
4r1 r2 · pPR (y | r1 )pPR (y | r2 ) · exp −
· I0
dr1 dr2
1−m
1−m
0
0
√ 
 √


Z +∞ Z +∞
2 z1 z2 m
√
√
z1 + z2
(a)
pPR (y | z1 )pPR (y | z2 ) · exp −
=
· I0
dz1 dz2
1−m
1−m
0
0

2

Z +∞ Z +∞ 
1
(y − z1 )2 + (y − z2 )2
(b)
√
=
exp −
2σ 2
σ 2π
−∞
−∞
√ 
 √


2 z1 z2 m
z1 + z2
· I0
· H(z1 )H(z2 )dz1 dz2
· exp −
1−m
1−m


 √


√ 
Z1 + Z2
2 Z1 Z2 m
(c)
= EZ1 ,Z2 exp −
· I0
· H(Z1 )H(Z2 ) ,
1−m
1−m
where in (a) we do the change of variables z1 = r12 and z2 = r22 ; in (b) we use the definition (9) and we
define H(x) = 1 if x ≥ 0 and H(x) = 0 otherwise; and in (c) we define (Z1 , Z2 ) ∼i.i.d. N(y, σ 2 ). In the
limit σ 2 → 0, we have that



 √

√ 
Z1 + Z2
2 Z1 Z2 m
EZ1 ,Z2 exp −
· I0
· H(Z1 )H(Z2 )
1−m
1−m


 √ 
2y
2y m
→ exp −
· I0
· H(y),
1−m
1−m
by Lebesgue’s dominated convergence theorem. As a result, by using (134), we obtain that


 √ 
Z +∞
2y
2y m
1
σ 2 →0
exp −
· I0
· exp(y) dy
f (m) −→
1−m 0
1−m
1−m
 

√
Z π Z +∞
1
1 + m − 2 m cos θ
=
exp −y
dy dθ
π(1 − m) 0 0
1−m
Z
1 π
1
1
√
=
.
dθ =
π 0 1 + m − 2 m cos θ
1−m
Consequently,
σ 2 →0

Fδ (m) −→ (1 − δ) log(1 − m),
which implies the desired result.
Proof of Corollary 1. We follow the proof of Theorem 1 presented in Section 4. The first step is exactly the
same, i.e., by applying Lemma 1, we show that (62) holds. On the contrary, the second step requires some
modifications, since the definition of the error metric is different. In particular, we will prove that
n
o
1
2
∗
∗
E
kE
{XX
}
−
E
{XX
|
Y
,
A
}k
= on (1).
1:n
1:n
Y
,A
1:n
1:n
F
d2

35

(141)

Furthermore, we have that
n
o
kE {XX ∗ } − E {XX ∗ | Y1:n , A1:n }k2F + E kXX ∗ − E {XX ∗ | Y1:n , A1:n }k2F
o
(a) n
≥ E kE {XX ∗ } − XX ∗ k2F
n
o
(b)
= E kId − XX ∗ k2F

(142)

(c)

= E {trace (Id − 2XX ∗ + XX ∗ XX ∗ )}

(d)

= d − 2d + d2 = d2 − d,

∗
where in (a) we use the triangle inequality,
in (c) we
p in (b) we use that E {XX } = Id by Lemma 12,
∗
∗
∗
usen that, for
o any matrix A, kAkF = trace(AA ), and in (d) we use that E {trace (XX XX )} =

E kXk4 = d2 . By applying (141) and (142), the proof of Corollary 1 is complete.
Let us now give the proof of (141). Similarly to (65), we have that

I(Yn+1 ; Y1:n , A1:n |An+1 )
(Z
Z
1
·
E
p(x
|
Y
,
A
)
pPR (yn+1 | x, Y1:n , A1:n+1 )ϕPR (yn+1 )dyn+1 dx
≥
1:n
1:n
Y1:n ,A1:n+1
2K 2
Cd
R
)
Z
Z
2
−
p(x) pPR (yn+1 | x, An+1 )ϕPR (yn+1 ) dyn+1 dx
,
Cd

R

(143)
where we define ϕPR (x) = x for |x| ≤ M , and ϕPR (x) = M · sign(x) otherwise. Then,
Z
Z
p(x | Y1:n , A1:n ) pPR (yn+1 | x, Y1:n , A1:n+1 ) · ϕPR (yn+1 ) dyn+1 dx
Cd
RZ
Z
(a)
=
p(x | Y1:n , A1:n ) pPR (yn+1 | x, Y1:n , A1:n+1 ) · yn+1 dyn+1 dx + E1
d
R
ZC
(b)
=
p(x | Y1:n , A1:n ) · |hAn+1 , xi|2 dx + E1
Cd
Z

(c)
= hAn+1 ,
p(x | Y1:n , A1:n ) · xx∗ dx An+1 i + E1
Cd

= hAn+1 , E {XX ∗ | Y1:n , A1:n } An+1 i + E1 ,
(144)
where in (a) we set
Z
Z
E1 =
p(x | Y1:n , A1:n ) pPR (yn+1 | x, Y1:n , A1:n+1 ) · (ϕPR (yn+1 ) − yn+1 ) dyn+1 dx,
Cd

R

in (b) we use the definition (9), and in (c) we use that |hAn+1 , xi|2 = hAn+1 , xx∗ An+1 i. Similarly, we
have that
Z
Z
p(x) pPR (yn+1 | x, An+1 )ϕPR (yn+1 ) dyn+1 dx = hAn+1 , E {XX ∗ } An+1 i + E2 ,
(145)
Cd

R

36

with

Z

Z

pPR (yn+1 | x, An+1 ) · (ϕPR (yn+1 ) − yn+1 ) dyn+1 dx.

p(x)

E2 =
Cd

R

By applying (144) and (145), we can rewrite the RHS of (143) as
n
o
1
2
∗
∗
·
E
E
|hA
,
(E
{XX
|
Y
,
A
}
−
E
{XX
})
A
i
+
E
−
E
|
n+1
1:n
1:n
n+1
1
2
Y1:n ,A1:n An+1
2K 2

n
o
 (146)


1
2
2
2
≥
· EY1:n ,A1:n EAn+1 |hAn+1 , M An+1 i| − EAn+1 |E1 | − EAn+1 |E2 |
,
2K 2

where we define M = E {XX ∗ | Y1:n , A1:n } − E {XX ∗ }. As K goes large, EAn+1 |Ei |2 tends to 0,
for i ∈ {1, 2}. Furthermore, we have that
d
n
o (a) X
1
∗
EAn+1 |hAn+1 , M An+1 i|2 =
Mij Mkl
· 2 (δij · δkl + δil · δjk )
d
i,j,k,l=1

1 
= 2 |trace(M )|2 + kM k2F
d
(b) 1
= 2 kM k2F ,
d

where in (a) we use the following definition of the Kronecker delta:

1, if a = b,
δab =
0, otherwise,

(147)

(148)

and in (b) we use that
trace(M ) =

d
X




E |Xi |2 | Y1:n , A1:n − E |Xi |2

i=1

=E

( d
X

)
|Xi |2 | Y1:n , A1:n

i=1

−E

( d
X

(149)

)
|Xi |2

= 0.

i=1

As a result, we conclude that (141) holds.

B

Proof of Corollary 2

First, we evaluate the RHS of (20). Then, we give the proof of Corollary 2.
Lemma 8 (Computation of δu for Phase Retrieval). Let δu (σ 2 ) be defined as in (20) and T ∗ (y) as in (23).
Assume that the distribution p(· | |g|) is given by (9). Then,
lim δu (σ 2 ) = 1.

σ→0

(150)

Proof. The proof boils down to computing expected values and integrals. By using (135), (139) and (140),
we immediately obtain that
Z +∞
 σ2 →0
EG {pPR (y | |G|)} =
2r · pPR (y | r) · exp −r2 dr −→ exp(−y) · H(y).
0

37

Similarly, we have that

EG pPR (y | |G|)|G|2 =

Z

+∞

 σ2 →0
2r3 · pPR (y | r) · exp −r2 dr −→ y · exp(−y) · H(y).

0

Thus,

Z
R


2
Z +∞
EG pPR (y | |G|)(|G|2 − 1)
(y − 1)2 exp(−y)dy = 1,
dy =
EG {p(yPR | |G|)}
0

which implies the desired result.
Proof of Corollary 2. Pick σ sufficiently small. Let G ∼ CN(0, 1), Y ∼ pPR (· | |G|) and Z = T (Y ),
where pPR is defined in (9) and T is a pre-processing function (possibly dependent on σ) that we will
choose later on. Assume that
(1) T (y) is upper and lower bounded by constants independent of σ;
(2) P(Z = 0) ≤ c1 < 1 and c1 is independent of σ;
(3) the condition (76) holds.
Then, by Lemma 2, we have that, as n → ∞,


0,


|hx̂, xi|2 a.s.
−→ ρ =
ψδ0 (λ∗δ )

kx̂k22 kxk2

,
 0 ∗
ψδ (λδ ) − φ0 (λ∗δ )

if ψδ0 (λ∗δ ) ≤ 0,
if ψδ0 (λ∗δ ) > 0,

(151)

where λ∗δ is the unique solution of the equation ζδ (λ) = φ(λ), and φ, ψδ and ζδ are defined in (72), (73) and
(75), respectively.
Let τ be the supremum of the support of Z. Assume also that, for λ̄δ < λ∗δ ,
(4) τ ≥ c2 > 0 and c2 is independent of σ;
(5) φ0 (λ∗δ ) is lower bounded by a constant independent of σ;
(6)

min

λ∈(min(λ̄δ ,λ∗δ ))

ψδ00 (λ) is lower bounded by a strictly positive constant independent of σ.

Let λ̄δ be the point at which ψδ attains its minimum. Then,
(a)

φ(λ̄δ ) − ψδ (λ̄δ ) = φ(λ̄δ ) − φ(λ∗δ ) + ζδ (λ∗δ ) − ψδ (λ̄δ )
(b)

= φ(λ̄δ ) − φ(λ∗δ ) + ζδ (λ∗δ ) − ζδ (λ̄δ )

(c) 0
= ζδ (x1 ) − φ0 (x1 ) · (λ∗δ − λ̄δ )

0
0

(d) ζδ (x1 ) − φ (x1 )
=
· ψδ0 (λ∗δ ) − ψδ0 (λ̄δ )
00
ψδ (x2 )

0
0
(e) ζδ (x1 ) − φ (x1 )
=
· ψδ0 (λ∗δ )
ψδ00 (x2 )
(f)

≤ c3 · ψδ0 (λ∗δ ),
38

(152)

where in (a) we use that ζδ (λ∗δ ) = φ(λ∗δ ), in (b) we use that ζδ (λ̄δ ) = ψδ (λ̄δ ), (c) holds for some x1 ∈
(λ̄δ , λ∗δ ) by the mean value theorem, (d) holds for some x2 ∈ (λ̄δ , λ∗δ ) by the mean value theorem, and in (e)
we use that ψδ0 (λ̄δ ) = 0. Note that (f) holds for some constant c3 independent of σ, as ζδ0 (x1 ) ≥ 0, ψδ00 (x2 )
is bounded, and φ0 (x2 ) < 0 since P(Z = 0) < 1.
As φ0 (λ∗δ ) is bounded, from (151) and (152) we deduce that

ρ ≥ c4 · φ(λ̄δ ) − ψδ (λ̄δ ) ,
(153)
for some constant c4 independent of σ. Notice that, if λ∗δ ≤ λ̄δ , then the right hand side is non-positive and
hence the lower bound still holds.
As τ > 0, we also have that λ̄δ > 0. Consider now the matrix Dn0 = Dn /α for some α > 0. Then, the
principal eigenvector of Dn0 is equal to the principal eigenvector of Dn . Hence, we can assume without loss
of generality that λ̄δ = 1. This condition can be rewritten as


Z2
1
E
= ,
(154)
2
(1 − Z)
δ
and (153) can be rewritten as


 
1
Z(|G|2 − 1)
−
.
ρ ≥ c4 · E
1−Z
δ
We set
T (y) = Tδ∗ (y, σ) ,

y+ − 1
p
,
y+ + δ c(σ) − 1

(155)

(156)

where y+ = max(y, 0) and c(σ) is a function of σ to be set as to satisfy Eq. (154). By substituting (156)
into (154), we get



Z2
1
E
=
E
(Y+ − 1)2 .
(157)
(1 − Z)2
δc(σ)
Hence, Eq. (154) is satisfied by
n

2 o
c(σ) = E (Y+ − 1)2 = E (|G|2 + σW )+ − 1
,

(158)

where W ∼ N(0, 1) is independent of G. Therefore, c(σ) is always well defined and, by dominated convergence, c(σ) → c(0) = 1, as σ → 0. Furthermore,



1
Z(|G|2 − 1)
=p
E (Y+ − 1)(|G|2 − 1) .
(159)
E
1−Z
δc(σ)
By applying again dominated convergence, we get


1 
Z(|G|2 − 1)
1
lim E
= √ E (|G|2 − 1)2 = √ .
σ→0
1−Z
δ
δ

(160)

Hence, by using (153), we get that, for δ > 1 and σ ≤ σ1 (δ),
lim inf
n→∞

 1
|hx̂σ , xi|2
1
√
≥
c
−
> 0,
5
kx̂σ k22 kxk22
δ δ
39

(161)

where x̂σ denotes the spectral estimator corresponding to the pre-processing function (156). Let us now
verify that, by setting T = Tδ∗ , the requirements stated above are fulfilled. As δ > 1, the function T is
bounded by constants independent of σ. It is also clear that the conditions (2) and (4) hold. Furthermore,
the conditions (5) and (6) follow by showing that φ(λ), ψδ (λ) have well defined uniform limits as σ → 0
that satisfy those conditions: this can be proved by one more application of dominated convergence.
In order to show that the condition (3) holds, we follow the argument presented at the end of the proof of
Lemma 2. First, we add a point mass with associated probability at most 1 , which immediately implies that
(76) is satisfied. Then, by applying the Davis-Kahan theorem [DK70], we show that we can take 1 = 0.
This proves the claim of the corollary for the pre-processing function Tδ∗ (y, σ), defined in (156). Let us
now prove that the same conclusion holds for Tδ∗ (y) defined in (25). Let
x+1
.
x+a

(162)



1
|a − 1|
≤
max
1,
.
(x + a)2
x2

(163)

fa (x) =
Then, for any x, a ∈ R≥0 ,
|fa0 (x)| =

p
Therefore, since Tδ∗ (y, σ) = 1 − fy+ ( δc(σ) − 1), we have that
sup Tδ∗ (y, σ) − Tδ∗ (y) ≤
y∈R

p
√
1
c(σ) − 1 .
p
√
2 δ ·
min( δc(σ) − 1, δ − 1, 1)

(164)

Denote by Dn (σ) and Dn the matrices constructed with the pre-processing functions Tδ∗ (y, σ) and Tδ∗ (y),
respectively. It follows that, for any δ > 1, there exists a function ∆(σ) with ∆(σ) → 0 as σ → 0 such that
kDn (σ) − Dn kop ≤ ∆(σ) .

(165)

Hence, by applying again the Davis-Kahan theorem, we conclude that, for all δ > 1 and σ ≤ σ2 (δ),
lim inf
n→∞

 1
|hx̂, xi|2
1
√
≥
c
−
> 0,
5
kx̂k22 kxk22
δ δ

(166)

where x̂ is the estimator corresponding to the pre-processing function Tδ∗ (y).

C

Auxiliary Lemmas

Lemma 9 (Distribution of Scalar Product of Two Unit Complex Vectors). Let x1 , x2 ∼i.i.d. Unif(Sd−1
C )
and define M = |hx1 , x2 i|2 . Then,
M ∼ Beta(1, d − 1).
(167)
Proof. Without loss of generality, we can pick x2 to be the first element of the canonical base of Cd . Thus,
M is equal to the squared modulus of the first component of x1 . Furthermore, we can think to x1 as being
chosen uniformly at random on the 2d-dimensional real sphere with radius 1. Note that, by taking a vector of
i.i.d. standard Gaussian random variables and dividing it by its norm, we obtain a vector uniformly random
on the sphere of radius 1. Hence,
U2 + U2
M = P12d 2 ,
2
i=1 Ui

with {Ui }1≤i≤2d ∼i.i.d. N(0, 1).
40

P
2
Set A = U12 + U22 and B = 2d
i=3 Ui . Then, A and B are independent, A follows a Gamma distribution
with shape 1 and scale 2, i.e., A ∼ Γ(1, 2), and B follows a Gamma distribution with shape d − 1 and scale
2, i.e., B ∼ Γ(d − 1, 2). Thus, we conclude that
M=

A
∼ Beta(1, d − 1),
A+B

which proves the claim.
Lemma 10 (Distribution of Scalar Product of Two Unit Real Vectors). Let x1 , x2 ∼i.i.d. Unif(Sd−1
R ) and
define M = hx1 , x2 i. Then, the distribution of M is given by
p(m) = √

Γ( d2 )
πΓ( d−1
2 )

(1 − m2 )

d−3
2

,

m ∈ [−1, 1].

(168)

Proof. Without loss of generality, we can pick x2 to be the first element of the canonical base of Rd . Thus,
M is equal to the first component of x1 . Note that, by taking a vector of i.i.d. standard Gaussian random
variables and dividing it by its norm, we obtain a vector uniformly random on the sphere of radius 1. Hence,
U2
M 2 = Pd 1

2
i=1 Ui

with {Ui }1≤i≤d ∼i.i.d. N(0, 1).

,

P
Set A = U12 and B = di=2 Ui2 . Then, A and B are independent, A follows a Gamma distribution with
shape 1/2 and scale 2, i.e., A ∼ Γ(1/2, 2), and B follows a Gamma distribution with shape (d − 1)/2 and
scale 2, i.e., B ∼ Γ((d − 1)/2, 2). Thus, we obtain that
M2 =

A
∼ Beta(1/2, (d − 1)/2).
A+B

A change of variable and the observation that the distribution of M is symmetric around 0 immediately let
us conclude that, for m ∈ [−1, 1],
d−3
p(m) = c · (1 − m2 ) 2 ,
(169)
where the normalization constant c is given by
Z

1
2

(1 − m )

c=
−1

d−3
2

−1
Γ( d2 )
dm
=√
.
πΓ( d−1
2 )

Lemma 11 (Laplace’s Method). Let F : [0, 1] → R be such that
• F is continuous;
• F (x) < 0 for x ∈ (0, 1];
• F (0) = 0.
Then,
Z
lim

n→+∞ 0

1

exp (n · F (x)) dx = 0.
41

(170)

Proof. Pick  > 0 and separate the integral into two parts:
Z
Z 
Z 1
exp (n · F (x)) dx +
exp (n · F (x)) dx =
0

1

exp (n · F (x)) dx.



0

Now, the first integral is at most  since F (x) ≤ 0 for any x ∈ [0, 1], and the second integral tends to 0 as
n → +∞ since F (x) < 0 for x ∈ (0, 1]. Thus, the claim immediately follows.
√ d−1
Lemma 12 (Second Moment of Uniform Vector on Complex Sphere). Let x ∼ Unif( dSC
). Then,
E {XX ∗ } = Id .

(171)

Proof. Let z ∼ CN(0d , Id ) and note that, by taking a vector of i.i.d. standard complex normal random
variables and dividing
it by its norm, we obtain a vector uniformly random on the complex sphere of radius
√
1. Then, x = dz/ kzk.
For i ∈ [d], denote by xi and by zi the i-th component of x and z, respectively. Then, for i 6= j,


Zi Zj∗

E Xi Xj∗ = d · E
= 0,
kZk2
where the last equality holds by symmetry. Furthermore,



|Zi |2
E |Xi |2 = d · E
= 1,
kZk2
as |Zi |2 / kZk2 ∼ Beta(1, d − 1) by the argument of Lemma 9. As a result, the thesis is readily proved.

D

Proof of Lemma 3

Before presenting the proof of the lemma, let us introduce some basic definitions and well-known results.
Let H be a probability measure on [0, +∞). Denote by ΓH the support of H and by τ the supremum of ΓH .
Let sH (g) denote the Stieltjes transform of H, which is defined as
Z
1
sH (g) =
dH(t),
(172)
t−g
and let gH (s) denote its inverse.
Consider a matrix
Sn =

1
U Mn U ∗ ,
d

(173)

and assume that
(1) Mn is PSD for all n ∈ N;
(2) U ∈ Cd×n is a random
whose entries {ui,j }1≤i≤d,1≤j≤n are i.i.d. such that E {Ui,j } = 0,
 matrix
2
4
E |Ui,j | = 1, and E |Ui,j | < ∞ (this includes the cases in which the entries are ∼i.i.d. CN(0, 1)
or are ∼i.i.d. N(0, 1));
(3) The sequence of empirical spectral distributions of Mn ∈ Cn×n converges weakly to a probability
distribution H, as n → +∞;
42

(4) n/d → δ ∈ (0, +∞), as n → ∞;
(5) The sequence of spectral norms of Mn is bounded.
Note that the normalization of (173) differs from the normalization of (80) by a factor of δ. However, since
the form (173) is more common in the literature, we will stick to it for the rest of this section. In order to
obtain the desired result for the matrix (80), it suffices to incorporate a factor 1/δ in the definition of the
function ψδ .
Let Fδ,H be the probability measure on [0, +∞) such that the inverse gFδ,H of its Stieltjes transform
sFδ,H is given by
Z
1
t
gFδ,H (s) = − + δ
dH(t), s ∈ {z ∈ C : =(z) > 0}.
(174)
s
1 + ts
Then, the sequence of empirical spectral distributions of Sn converges weakly to Fδ,H [MP67], [SB10,
Chapter 4].
For α 6∈ ΓH and α 6= 0, let us also define


1
.
(175)
ψFδ,H (α) = gFδ,H −
α
The function ψFδ,H links the support of Fδ,H with the support of the generating measure H (see [SC95,
Section 4] and [BY12, Lemma 3.1]). In particular, if λ 6∈ ΓFδ,H , then sFδ,H (λ) 6= 0 and α = −1/sFδ,H (λ)
satisfies
(1) α 6∈ ΓH and α 6= 0 (so that ψFδ,H (α) is well-defined);
(2) ψF0 δ,H (α) > 0.
Conversely, if α satisfies (1) and (2), then λ = ψFδ,H (α) 6∈ ΓFδ,H .
n
denote the largest eigenvalue of Mn and assume that, as n → ∞,
Let λM
1
a.s.

n
λM
−→ α∗ 6∈ ΓH .
1

(176)

Denote by λ1Sn the largest eigenvalue of Sn . Then, the results in [BY12] prove that
a.s.

if ψF0 δ,H (α∗ ) > 0,

a.s.

if ψF0 δ,H (α∗ ) ≤ 0.

n
λS
1 −→ λ∗ = ψFδ,H (α∗ ),
n
λS
1 −→ min ψFδ,H (α),

α>τ

(177)

n
Informally, the eigenvalue λM
is mapped into the point ψFδ,H (α∗ ), where α∗ = −1/sFδ,H (λ∗ ). This point
1
emerges from the support of Fδ,H if and only if ψF0 δ,H (α∗ ) > 0.
In what follows, we relax the first hypothesis, i.e., we consider the case in which the matrix Mn is not
PSD. We will show that (177) still holds, which implies the claim of Lemma 3.

Proof of Lemma 3. As U is drawn from a rotationally invariant distribution, we can assume without loss of
generality that Mn is diagonal. Then, we have that

 ∗ 
Mn+
0k
U+
Sn = (U+ , U− )
0n−k −Mn−
U−∗
(178)
1
1
+ ∗
− ∗
= U+ Mn U+ − U− Mn U− ,
d
d
43

where Mn+ ∈ Rk×k is the diagonal matrix containing the positive eigenvalues of Mn , Mn− ∈ R(n−k)×(n−k)
is the diagonal matrix containing the negative eigenvalues of Mn with the sign changed, U+ contains the
first k columns of U , and U− contains the remaining n − k columns of U .
Note that U+ and U− are independent. Furthermore, if H is a unitary matrix, then U− and HU− have
the same distribution. Hence, we can rewrite the matrix Sn as
Sn =

1
1
U1 Mn+ U1∗ − HU2 Mn− U2∗ H ∗ ,
d
d

(179)

where U1 and U2 are independent with entries ∼i.i.d. CN(0, 1), and H is a random unitary matrix distributed
according to the Haar measure.
Recall that, by hypothesis, the sequence of empirical spectral distributions of Mn converges weakly
to the probability distribution H, where H is the law of the random variable Z. Then, the sequence of
empirical spectral distributions of Mn+ converges weakly to the probability distribution H + , where H + is
the law of Z + = max(Z, 0). Let Fδ,H + be the probability measure on [0, +∞) such that the inverse gFδ,H +
of its Stieltjes transform sFδ,H + is given by
1
gFδ,H + (s) = − + δ
s

Z

t
dH + (t).
1 + ts

(180)

Define Sn+ = d1 U1 Mn+ U1∗ . Then, as Mn+ is PSD, the sequence of empirical spectral distributions of Sn+
converges weakly to Fδ,H + [MP67], [SB10, Chapter 4].
Similarly, the sequence of empirical spectral distributions of Mn− converges weakly to the probability
distribution H − , where H − is the law of Z − = − min(Z, 0). Let Fδ,H − be the probability measure on
[0, +∞) such that the inverse gFδ,H − of its Stieltjes transform sFδ,H − is given by
1
gFδ,H − (s) = − + δ
s

Z

t
dH − (t).
1 + ts

(181)

Define Sn− = d1 U2 Mn− U2∗ . Then, as Mn− is PSD, the sequence of empirical spectral distributions of Sn−
converges weakly to Fδ,H − [MP67], [SB10, Chapter 4]. Furthermore, the sequence of empirical spectral
distributions of −Sn− converges weakly to the probability measure Fδ,H − such that
inv

gF

δ,H −
inv

where gF

δ,H −
inv

(s) = −gFδ,H − (−s),

denotes the inverse of the Stieltjes transform sF

Define

δ,H −
inv

Fδ,H = Fδ,H +  Fδ,H − ,
inv

(182)
of Fδ,H − .
inv

(183)

where  denotes the free additive convolution. Recall the decomposition (179). Then, the sequence of
empirical spectral distributions of Sn converges weakly to Fδ,H [Voi91, Spe93]. Consequently, the inverse

44

gFδ,H of the Stieltjes transform sFδ,H of Fδ,H can be computed as
(a)

gFδ,H (s) = gFδ,H + F

δ,H −
inv

(s)

1
= gFδ,H + (s) + gF − (s) +
δ,H
s Z
inv
Z
t
t
(c) 1
+
dH (t) − δ
dH − (t)
=− +δ
s
1 + ts
1 − ts
Z
Z
t
t
(d) 1
+
=− +δ
dH (t) + δ
dH − (−t)
s
1 + ts
1 + ts
Z
t
(e) 1
=− +δ
dH(t),
s
1 + ts

(b)

(184)

where in (a) we use (183), in (b) we use that the R-transform of the free convolution is the sum of the
R-transforms of the addends, in (c) we use (180), (181), and (182), in (d) we perform the change of variable
t → −t in the second integral; and in (e) we use the fact that H + (t) is the law of max(Z, 0), H − (−t) is the
law of min(Z, 0), and that t/(1 + ts) = 0 for t = 0.
n a.s.
By hypothesis, λM
−→ α∗ 6∈ ΓH . First, we establish under what condition the largest eigenvalue
1
+

of Sn+ , call it λ1Sn , converges to a point outside the support of Fδ,H + . To do so, define ψFδ,H + (α) =
+

+

a.s.

n
gFδ,H + (−1/α). Then, λ1Sn −→ ψFδ,H + (α∗ ), if ψF0 + (α∗ ) > 0; and λS
converges almost surely to a
1
δ,H
point inside the support of Fδ,H + , otherwise [BY12].
For the moment, assume that ψF0 + (α∗ ) > 0. We now establish under what condition the largest
δ,H

n
eigenvalue of Sn , call it λS
1 , converges to a point outside the support of Fδ,H . To do so, let ω1 and ω2
denote the subordination functions corresponding to the free convolution Fδ,H +  Fδ,H − . These functions
inv
satisfy the following analytic subordination property:

sFδ,H + F

δ,H −
inv

(z) = sFδ,H + (ω1 (z)) = sF

δ,H −
inv

(ω2 (z)).

(185)

Then, by Theorem 2.1 of [BBCF15], we have that the spike ψFδ,H + (α∗ ) is mapped into ω1−1 (ψFδ,H + (α∗ )).
The Stieltjes transform at this point is given by
(a)

sFδ,H + F

δ,H −
inv

(ω1−1 (ψFδ,H + (α∗ ))) = sFδ,H + (ψFδ,H + (α∗ ))
(b)

= sFδ,H + (gFδ,H + (−1/α∗ ))

(c)

= −1/α∗ ,

where in (a) we use (185), in (b) we use the definition of ψFδ,H + , and in (c) we use that gFδ,H + is the
functional inverse of the Stieltjes transform sFδ,H + . As a result, by [SC95, Section 4], we conclude that
ω1−1 (ψFδ,H + (α∗ )) 6∈ ΓFδ,H if and only if ψF0 δ,H (α∗ ) > 0. Furthermore, the condition ψF0 δ,H (α∗ ) > 0 is
more restrictive than the condition ψF0 + (α∗ ) > 0 since
δ,H

ψF0

δ,H

Z 
(α∗ ) = 1 − δ
+

t
α∗ − t

2

dH + ≥ 1 − δ

Z 

t
α∗ − t

2

dH = ψF0 δ,H (α∗ ).

0
n
Hence, λS
1 converges to a point outside the support of Fδ,H if and only if ψFδ,H (α∗ ) > 0 and the proof is
complete.

45

Remark 7 (Lemma 3 for the Real Case). Consider the random matrix n1 U Mn U T , where U ∈ R(d−1)×n is
a random matrix whose entries are ∼i.i.d. N(0, 1) and Mn ∈ Rn×n . Then, the claim of Lemma 3 still holds.
Let us briefly explain why this is the case.
If Mn is PSD, then the results of [BY12] allow us to conclude. If Mn is not PSD, we can write an
expression analogous to (179):
1
1
1
U Mn U T = U1 Mn+ U1T − HU2 Mn− U2T H ∗ ,
d
d
d

(186)

where Mn+ is the diagonal matrix containing the positive eigenvalues of Mn , Mn− is the diagonal matrix
containing the negative eigenvalues of Mn with the sign changed, U1 and U2 are independent with entries
∼i.i.d. N(0, 1), H is a random unitary matrix distributed according to the Haar measure, and we have used
the fact that the eigenvalues of U2 Mn− U2T are the same as the eigenvalues of HU2 Mn− U2T H ∗ since H is
unitary. Hence, the proof follows from the same argument of Lemma 3.

E

Proof of Lemma 4 and Theorem 5

We start by proving a result similar to Lemma 4 for a general AMP iteration, where the function ft (ẑ; y) is
generic.
Lemma √
13 (State Evolution for General AMP Iteration). Let x ∈ Rd denote the unknown signal such that
kxk2 = d, A = (a1 , . . . , an )T ∈ Rn×d with {ai }1≤i≤n ∼i.i.d. N(0d , Id /d), and y = (y1 , . . . , yn ) with
yi ∼ p(· | hx, ai i). Consider the AMP iterates z t , ẑ t defined in (116) for some function ft (ẑ; y), with bt
given by
bt = δ · E{ft0 (µt G0 + τt G1 ; Y )} ,
(187)
where the expectation is with respect to G0 , G1 ∼i.i.d. N(0, 1) and Y ∼ p( · | G0 ). Assume that the
initialization z 0 is independent of A and that, almost surely,
lim

n→∞

1
hx, z 0 i = µ0 ,
d

lim

n→∞

1 0 2
kz k = µ20 + τ02 .
d

Let the state evolution recursion τt , µt be defined as
Z

E ∂g p(y | X0 )ft (µt X0 + τt G; y) dy ,
µt+1 = δ
R
n
2 o
2
τt+1 = δ · E ft (µt X0 + τt G; Y )
,

(188)

(189)

with initialization µ0 and τ0 , where the expectation is taken with respect to X0 , G ∼i.i.d. N(0, 1). Then, for
any t, and for any function ψ : R2 → R such that |ψ(u) − ψ(v)| ≤ L(1 + kuk2 + kvk2 )ku − vk2 for some
L ∈ R, we have that, almost surely,
n

1X
ψ(xi , zit ) = E {ψ(X0 , µt X0 + τt G)} .
n→∞ n
lim

i=1

Proof. For g ∈ R, let H( · ; g) : [0, 1] → R ∪ {+∞, −∞} be the generalized inverse of
Z y
F(y | g) ≡
p(y 0 | g) dy 0 ,
−∞

46

(190)

namely,

H(w; g) ≡ inf y ∈ R : F(y | g) ≥ w .

(191)

With this definition, the model yi ∼ p( · | hai , xi) is equivalent to yi = H(wi ; hai , xi) for {wi }1≤i≤n ∼i.i.d.
Unif([0, 1]) independent of A and x. Let w = (w1 , . . . , wn ) ∈ Rn and denote by [v1 | · · · | vk ] ∈ Rm×k
the matrix obtained by stacking column vectors v1 , . . . , vk ∈ Rm .
For t ≥ 0, define r t = 0d , r̂ t = Ax, and introduce the extended state variables st ∈ Rd×2 and
t
ŝ ∈ Rn×2 , defined as
st = [z t | r t ],
ŝt = [ẑ t | r̂ t ].

(192)

We further define the functions ht = [ht,1 | ht,2 ] : R2 × R → R2 and ĥt = [ĥt,1 | ĥt,2 ] : R2 × R → R2 by
setting
ht (s1 , s2 ; x) ≡ [s1 | x ] ,
ĥt (ŝ1 , ŝ2 ; w) ≡ [ft (ŝ1 ; H(w; ŝ2 )) | 0 ] .

(193)

With these notations, the iteration (116) is equivalent to
st+1 = AT ĥt (ŝt ; w) − ht (st ; x)B̂t ,
ŝt = Aht (st ; x) − ĥt−1 (ŝt−1 ; w)Bt−1 ,

(194)

where the functions ht (st ; x) and ĥt (ŝt ; w) are understood to be applied component-wise to their arguments
and Bt , B̂t ∈ R2×2 are defined by
)
(
∂ ĥt,k
(µt X0 + τt G, X0 ; W ) ,
(B̂t )j,k = δ · E
∂ŝj
(195)


∂ht,k
(Bt )j,k = δ · E
(µt X0 + τt G, 0; X0 ) .
∂sj
The iteration (194) satisfies the assumptions of [JM13][Proposition 5]. By applying that result, the claim
follows.
At this point, first we present the proof of Lemma 4 and then of Theorem 5.
Proof of Lemma 4. Consider the state evolution recursion defined in (114) with initialization µ0 . Let ft be
defined as in (117) with F given by (113). Suppose that, for any t, (189) holds with µt = τt2 . Then, by
Lemma 13, the claim immediately follows.
The remaining part of the proof is devoted to show that (189) holds with µt = τt2 , for t ≥ 0. First,
we prove by induction that µt = τt2 , for t ≥ 0. The basis of the induction, i.e., µ0 = τ02 , is true by the
2 . Set
hypothesis of the Lemma. Now, we assume that µt = τt2 and we show that µt+1 = τt+1
Z = µt X0 + τt G,
and note that Z ∼ N(0, µ2t + τt2 ). Then, we can re-write X0 as
e
X0 = aZ + bG,
47

(196)

e ∼ N(0, 1) and independent from Z. In order to compute the coefficients a and
for some a, b ∈ R, where G
2
b, we evaluate E{X0 } and E{X0 · Z}, thus obtaining the equations
a2 (µ2t + τt2 ) + b2 = 1,
a(µ2t + τt2 ) = µt ,
which can be simplified as
a=

µ2t

µt
,
+ τt2
τt

.
b= p
µ2t + τt2
Furthermore, by using the inductive hypothesis µt = τt2 and that qt = µt /(1 + µt ), we obtain that
p
e
X0 = (1 − qt ) Z + 1 − qt G.

(197)

Hence, the following chain of equalities holds:
Z n
2 o
2 (a)
τt+1 = δ
E p(y | X0 ) · ft (µt X0 + τt G; y)
dy
ZR n
p

 o
(b)
e · ft (Z; y) 2 dy
=δ
E p y | (1 − qt ) Z + 1 − qt G
ZR n
o
p

2 
(c)
e Z dy
=δ
E ft (Z; y) · E p y | (1 − qt ) Z + 1 − qt G
ZR (
)
√
e | Z} 2 
p

E{∂g p(y | (1 − qt ) Z + 1 − qt G)
(d)
e Z
=δ
E
dy
· E p y | (1 − qt ) Z + 1 − qt G
√
e | Z}
E{p(y
|
(1
−
q
)
Z
+
1
−
q
G)
t
t
R
Z (
)
√
e | Z} 
p

E{∂g p(y | (1 − qt ) Z + 1 − qt G)
e Z
=δ
E
· E ∂g p y | (1 − qt ) Z + 1 − qt G
dy
√
e
E{p(y
|
(1
−
q
)
Z
+
1
−
q
G)
|
Z}
t
t
R
Z n
o
p


(e)
e Z dy
=δ
E ft (Z; y) · E ∂g p y | (1 − qt ) Z + 1 − qt G
ZR n
o
(f)
=δ
E ft (µt X0 + τt G; y) · ∂g p(y | X0 ) dy = µt+1 ,
R

(198)
where in (a) we use that Y ∼ p(· | X0 ), in (b) we use (196) and (197), in (c) we condition with respect to
Z, in (d) we use the definition (117) of ft , in (e) we use again the definition (117) of ft , and in (f) we use
again (196) and (197).
Finally, we prove that µt+1 satisfies (189). Indeed, the following chain of equalities holds:
Z (
 )
√
e | Z} 2
E{∂g p(y | (1 − qt ) Z + 1 − qt G)
(a)
µt+1 = δ
E
dy
√
e | Z}
E{p(y
|
(1
−
q
)
Z
+
1
−
q
G)
t
t
ZR (
p
2 )
√
E{∂g p(y | (1 − qt ) µ2t + τt2 G0 + 1 − qt G1 ) | G0 }
(b)
p
(199)
dy
=δ
E
√
E{p(y | (1 − qt ) µ2t + τt2 G0 + 1 − qt G1 ) | G0 }
R
Z (
2 )
√
√
E{∂g p(y | qt G0 + 1 − qt G1 ) | G0 }
(c)
√
E
dy = δ · h(qt ),
=δ
√
E{p(y | qt G0 + 1 − qt G1 ) | G0 }
R
48

p
e and G0 = Z/ µ2 + τ 2 , and in (c) we use that µt = τ 2
where in (a) we use (198), in (b) we set G1 = G
t
t
t
and that qt = µt /(1 + µt ).
Proof of Theorem 5. In view of Lemma 4, it is sufficient to show that (q, µ) = (0, 0) is an attractive fixed
point of the recursion (114).
First of all, let us check that (q, µ) = (0, 0) is a fixed point. This happens if and only if
2
Z
EG1 {∂g p(y | G1 )}
h(0) =
dy = 0,
(200)
EG1 {p(y | G1 )}
R
which holds because of the condition (121).
Let us now prove that this fixed point is stable. We start by re-writing the function h(q) defined in (115)
as
(
2 )
√
Z
hnum ( q, y)
h(q) =
EG0
dy ,
(201)
√
hden ( q, y)
R
where
n
o
p
hnum (x, y) = EG1 ∂g p(y | x · G0 + 1 − x2 G1 ) ,
o
n
p
hden (x, y) = EG1 p(y | x · G0 + 1 − x2 G1 ) .

(202)

Note that hnum (0, y) = 0 by assumption (121). Then,
√ ∂hnum (x, y)
√
hnum ( q, y) = q
∂x

+
x=0

q ∂ 2 hnum (x, y)
2
∂2x

√ ∂hden (x, y)
hden (x, y) = hden (0, y) + q
∂x

,
x=x1

(203)

,
x=x2

√
for some x1 , x2 ∈ [0, q]. Furthermore, by applying Stein’s lemma, we have that
n
o
p
1
hnum (x, y) = √
EG1 G1 · p(y | x · G0 + 1 − x2 G1 ) .
1 − x2

(204)

By using (204), we can re-write (203) as
1
q
√
√
hnum ( q, y) = q G0 · EG1 {G1 · ∂g p(y | G1 )} +
EG {fnum (G0 , G1 , x1 )} ,
2 (1 − x21 )5/2 1
√
hden (x, y) = EG1 {p(y | G1 )} + q EG1 {fden (G0 , G1 , x2 )} ,
where
fnum (G0 , G1 , x1 ) = G1 (1 +

2x21 ) p(y

q
| x1 · G0 + 1 − x21 G1 )

q
q

2
2
− 2G0 ·
− 1) + G1 1 − x1 (1 + 2x1 ) ∂g p(y | x1 · G0 + 1 − x21 G1 )
!
q
q
(205)
+ (x21 − 1)(G20 (x21 − 1) − G21 x21 + 2G0 G1 x1 1 − x21 )∂g2 p(y | x1 · G0 + 1 − x21 G1 ) ,
x1 (x21

fden (G0 , G1 , x2 ) =

x2

G0 − p
G1
1 − x22

!
· ∂g p(y | x2 · G0 +

49

q

1 − x22 G1 ).

By applying again Stein’s lemma and by using that the conditional density p(y | g) is bounded, we note that
EG1 {fnum (G0 , G1 , x1 )} and EG1 {fden (G0 , G1 , x2 )} are bounded. Hence, by dominated convergence, we
obtain that
2
Z
EG1 {∂g2 p(y | G1 )}
h(q) = q ·
dy + o(q).
EG1 {p(y | G1 )}
R
Therefore, in a neighborhood of the fixed point we have
qt = µt + o(µt ) ,
2
Z
EG1 {∂g2 p(y | G1 )}
µt+1 = δ · qt ·
dy + o(qt ) .
EG1 {p(y | G1 )}
R

(206)

Furthermore, by applying twice Stein’s lemma, we also have that
EG1 {∂g2 p(y | G1 )} = EG1 {p(y | G1 )(G21 − 1)}.

(207)

By using (206), (207) and by recalling the definition (42) of δu , we conclude that
qt = µt + o(µt ) ,
δ
µt+1 =
qt + o(qt ).
δu

(208)

As δ < δu , the fixed point is stable.

F

Proof of Lemma 5 and Theorem 6

For the proofs in this section, it is convenient to introduce the function
√


√
EG {∂g2 p(y | q̄ x + q̄G)}
EG {∂g p(y | q̄ x + q̄G)} 2
√
√
G(x, y; q̄) =
−
.
EG {p(y | q̄ x + q̄G)}
EG {p(y | q̄ x + q̄G)}

(209)

First, we present the proof of Lemma 5 and then of Theorem 6.
Proof of Lemma 5. The condition (121) implies that
F(0, y; 1) = 0.

(210)

Furthermore, we have that
(

(a)

EY {G(0, Y ; 1)} = EY
(b)

EG {∂g2 p(Y | G)}
EG {p(Y | G)}

Z

)

EG {∂g2 p(y | G)} dy
 Z

2
= EG ∂g
p(y | G) dy = 0,
=

(211)

R

R

where in (a) we use (210) and the definition (209) of G(0, y; 1), and in (b) we use the fact that y has density
EG {p(y | G)}.
50

Denote by F0 (x, y; q̄) the derivative of F with respect to its first argument. Then, we have
F0 (x, y; q̄) = q̄G(x, y; q̄) .

(212)

Hence,
bt = δ · (1 − qt ) · E{G(µt G0 +

√

µt G1 , Y ; 1 − qt )} = δ · E{G(0; Y ; 1)} + oqt (1) = oqt (1) .

(213)

By using (212) and (213), we linearize the recursion (116) around the fixed point z t = 0d and ẑ t−1 = 0n
as
z t+1 = AT J ẑ t + oqt (1)(kz t k2 + kẑ t k2 ) + o(kẑ t k2 ) ,
ẑ t = Az t − J ẑ t−1 + oqt (1) kẑ t−1 k2 + o(kẑ t−1 k2 ) ,

(214)
(215)

where J ∈ Rn×n is a diagonal matrix with entries ji = F0 (0, yi ; 1) for i ∈ [n]. By substituting the
expression (215) for ẑ t into the RHS of (214), the result follows.
Proof of Theorem 6. By definition, α is an eigenvalue of Ln if and only if
det(Ln − αIn+d ) = 0.
Recall that, when D is invertible,


A B
= det(D) · det(A − BD −1 C).
det
C D

(216)

(217)

Then, after some calculations, we obtain that (216) is equivalent to
αd · det(−J − αIn ) · det(Id − AT (In + αJ −1 )−1 A) = 0.

(218)

From (218), we immediately deduce that the eigenvalues of Ln are real if and only if all the solutions to
det(Id − AT (In + αJ −1 )−1 A) = 0

(219)

are real. We will prove that in fact this equation does not have any solution for α ∈ C \ R.
Let U ΣV T be the SVD of A. Then, (219) is equivalent to
det(U Σ−2 − (In + αJ −1 )−1 U ) = 0.
Using the fact that det(Σ) 6= 0, and det(In + αJ −1 ) 6= 0 for α ∈ C \ R, Eq. (219) is equivalent to
det((In + αJ −1 )U − U Σ2 ) = 0 ,
or equivalently
det(In + αJ −1 − AAT ) = det(In + αJ −1 − U Σ2 U T ) = 0 .
Given that the solutions of this equations are generalized eigenvalues for the pairs of symmetric matrices
AAT − In and J −1 , they must be real. We conclude that the eigenvalues of Ln are real.

51

Note that
EG {∂g2 p(y | G)}
EG {p(y | G)}
2
(b) EG {p(y | G)(G − 1)}
=
EG {p(y | G)}
∗
(c) T (y)
,
=
1 − T ∗ (y)
(a)

G(0, y; 1) =

where in (a) we use that F(0, 1; y) = 0 as (121) holds, in (b) we apply twice Stein’s lemma, and in (c) we
use the definition (45) of T ∗ . Then, (219) can be re-written as
!
n
X
T ∗ (yi )
det Id −
ai aT
= 0.
(220)
i
T ∗ (yi ) + α(1 − T ∗ (yi ))
i=1

D∗

Let λ1 n (α) be the largest eigenvalue of the matrix Dn∗ (α) defined as
Dn∗ (α) =

n
X
i=1

T ∗ (yi )
ai aT
i .
T ∗ (yi ) + α(1 − T ∗ (yi ))

(221)

Note that, as α → +∞, the entries of Dn∗ (α) tend to 0 with high probability. Since the eigenvalues of a
matrix are continuous functions of the elements of the matrix, we also obtain that
D∗

lim λ1 n (α) = 0.

α→+∞
D∗

D∗

Hence, if there exists ᾱ > 1 such that λ1 n (ᾱ) > 1, then there exists also ᾱ0 > ᾱ > 1 such that λ1 n (ᾱ0 ) =
1. Consequently, there exists α > 1 that satisfies (220),pwhich implies the result of the theorem.
The rest
p of the proof consists in showing that ᾱ = δ/δu satisfies the desired requirements. First of all,
note that δ/δu > 1, as δ > δu . Furthermore, we have that
Dn∗ (ᾱ) =

n
X

Tδ∗ (yi )ai aT
i ,

(222)

i=1

√
where Tδ∗ is defined in (44). Recall that, by hypothesis, x is such that kxk2 = d and {ai }1≤i≤n ∼i.i.d.
√
√
en
N(0d , Id /d). Let x̃ = x/ d and ãi = d · ai . Then, hx, ai i = hx̃, ãi i. Let λD
1 be the largest eigenvalue
e n defined as
of the matrix D
n
X
en = 1
D
Tδ∗ (yi )ãi ãT
(223)
i .
n
i=1

e n = D ∗ (ᾱ)/δ, it remains to prove that
Since D
n
1
e n a.s.
λD
−→ λ̃ > .
1
δ

(224)

To do so, we apply a result analogous to that of Lemma 2 for the real case with T = Tδ∗ . For the
moment, assume that Tδ∗ fulfills the hypotheses of Lemma 2 (we will prove later that this is the case). Then,
n
λD
converges almost surely to ζδ (λ∗δ ).
1

e

52

Recall that
ζδ (λ) = ψδ (max(λ, λ̄δ )),
where λ̄δ is the point of minimum of the convex function ψδ (λ) defined as



Tδ∗ (Y )
1
.
ψδ (λ) = λ
+E
δ
λ − Tδ∗ (Y )
Notice also that this minimum is the unique local minimizer since ψδ is convex and analytic.
Furthermore, λ∗δ is the unique solution to the equation ζδ (λ∗δ ) = φ(λ∗δ ), where φ(λ) is defined as
 ∗

Tδ (Y ) · G2
φ(λ) = λ · E
.
λ − Tδ∗ (Y )
By setting the derivative of ψδ (λ) to 0, we have that


(Tδ∗ (Y ))2
1
= .
E
∗
2
δ
(λ̄δ − Tδ (Y ))
By using the definition (44) of Tδ∗ and the definition (45) of T ∗ , we verify that
r
Tδ∗ (Y )
δu EG {p(y | G)(G2 − 1)}
=
.
1 − Tδ∗ (Y )
δ
EG {p(y | G)}

(225)

Hence, by using the definition (42) of δu , we obtain that

E

(Tδ∗ (Y ))2
(1 − Tδ∗ (Y ))2



δu
=
δ

Z
R

2
EG {p(y | G)(G2 − 1)}
1
dy = ,
EG {p(y | G)}
δ

which immediately implies that
λ̄δ = 1.

(226)

By using (225), one also obtains that
r
 r Z

Tδ∗ (Y )
δu
δu
2
=
EG {p(y | G)(G − 1)} dy =
EG {G2 − 1} = 0,
E
1 − Tδ∗ (Y )
δ R
δ
which implies that
1
ψδ (1) = .
δ

(227)

Furthermore, we have that
2
 ∗
 r Z
EG {p(y | G)(G2 − 1)}
Tδ (Y )(G2 − 1)
δu
1
1
=
E
dy = √
> ,
1 − Tδ∗ (Y )
δ R
EG {p(y | G)}
δ
δ · δu
which implies that
φ(1) = √

1
1
> ,
δ
δ · δu

(228)

as δ > δu . By putting (226), (227), and (228) together, we obtain that
φ(λ̄δ ) > ζδ (λ̄δ ).
53

(229)

Recall that ζδ (λ) is monotone non-decreasing and φ(λ) is monotone non-increasing. Consequently, (229)
implies that λ∗δ > λ̄δ . Thus, we conclude that
n
= ζδ (λ∗δ ) = ψδ (λ∗δ )
lim λD
1

e

n→∞

> ψδ (λ̄δ ) = ψδ (1) =

1
.
δ

(230)

Now, we show that Tδ∗ fulfills the hypotheses of Lemma 2 by using arguments similar to those at the
end of the proof of Theorem 2. First of all, since T ∗ (y) ≤ 1, we have that Tδ∗ (y) is bounded. Furthermore,
if Tδ∗ (y) is equal to the constant value 0, then δu = ∞ and the claim of Theorem 6 trivially holds. Hence,
we can assume that P(Tδ∗ (Y ) = 0) < 1. Let τ be the supremum of the support of Tδ∗ (Y ). If P(Tδ∗ (Y ) =
τ ) > 0, then the condition (76) is satisfied and the proof is complete. Otherwise, for any 1 > 0, there exists
∆1 (1 ) such that Eq. (109) holds. Define Tδ∗ (y, 1 ) as in (110). Clearly, the random variable Tδ∗ (Y, 1 ) has
a point mass at δ, hence the condition (76) is satisfied. As a final step, we show that we can take 1 ↓ 0.
Define
n
1X ∗
e
Dn (1 ) =
Tδ (yi , 1 )ãi ã∗i .
n
i=1

Then,
e n (1 ) − D
en
D

op

≤ C1 · ∆1 (1 ),

(231)

where the constant C1 depends only on n/d. Consequently, by using (231) and Weyl’s inequality, we
conclude that
e ( )
en
D
(232)
|λ1 n 1 − λD
1 | ≤ C1 · ∆1 (1 ).
e n (1 ) tends to the largest eigenvalue of D
e n,
Hence, for any n, as 1 tends to 0, the largest eigenvalue of D
which concludes the proof.

References
[AGMM15] Sanjeev Arora, Rong Ge, Tengyu Ma, and Ankur Moitra, Simple, efficient, and neural algorithms for sparse coding, Conference on Learning Theory (COLT) (Paris, France), July 2015,
pp. 113–149.
[BBCF15]

Serban T. Belinschi, Hari Bercovici, Mireille Capitaine, and Maxime Février, Outliers in the
spectrum of large deformed unitarily invariant models, arXiv:1412.4916, 2015.

[BCE06]

Radu Balan, Pete Casazza, and Dan Edidin, On signal reconstruction without phase, Applied
and Computational Harmonic Analysis 20 (2006), no. 3, 345–356.

[BCMN14] Afonso S. Bandeira, Jameson Cahill, Dustin G. Mixon, and Aaron A. Nelson, Saving phase:
Injectivity and stability for phase retrieval, Applied and Computational Harmonic Analysis 37
(2014), no. 1, 106–125.
[BGN11]

Florent Benaych-Georges and Raj R. Nadakuditi, The eigenvalues and eigenvectors of finite,
low rank perturbations of large random matrices, Advances in Mathematics 227 (2011), no. 1,
494–521.
54

[BKM+ 17] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová, Phase
transitions, optimal errors and optimality of message-passing in generalized linear models,
arXiv:1708.03395, 2017.
[BLM15]

Mohsen Bayati, Marc Lelarge, and Andrea Montanari, Universality in polytope phase transitions and message passing algorithms, Annals of Applied Probability 25 (2015), no. 2, 753–
822.

[BM11]

Mohsen Bayati and Andrea Montanari, The dynamics of message passing on dense graphs,
with applications to compressed sensing, IEEE Trans. Inform. Theory 57 (2011), 764–785.

[BM12]

Mohsen Bayati and Andrea Montanari, The LASSO risk for Gaussian matrices, IEEE Trans.
Inform. Theory 58 (2012), no. 4, 1997–2017.

[BMDK17] Jean Barbier, Nicolas Macris, Mohamad Dia, and Florent Krzakala, Mutual information and
optimality of approximate message-passing in random linear estimation, arXiv:1311.2445,
2017.
[BR17]

Sohail Bahmani and Justin Romberg, Phase retrieval meets statistical learning theory: A flexible convex relaxation, Proc. of the 20th International Conference on Artificial Intelligence and
Statistics (AISTATS) (Fort Lauderdale, FL), June 2017, pp. 252–260.

[BY12]

Zhidong Bai and Jianfeng Yao, On sample eigenvalues in a generalized spiked population
model, Journal of Multivariate Analysis 106 (2012), 167–177.

[CC15]

Yuxin Chen and Emmanuel J. Candès, Solving random quadratic systems of equations is nearly
as easy as solving linear systems, Advances in Neural Information Processing Systems, 2015,
pp. 739–747.

[CC16]

, The projected power method: An efficient algorithm for joint alignment from pairwise
differences, arXiv:1609.05820, 2016.

[CC17]

, Solving random quadratic systems of equations is nearly as easy as solving linear
systems, Communications on Pure and Applied Mathematics 70 (2017), 0822–0883.

[CEHV15] Aldo Conca, Dan Edidin, Milena Hering, and Cynthia Vinzant, An algebraic characterization
of injectivity in phase retrieval, Applied and Computational Harmonic Analysis 38 (2015),
no. 2, 346–356.
[CESV15]

Emmanuel J. Candès, Yonina C. Eldar, Thomas Strohmer, and Vladislav Voroninski, Phase
retrieval via matrix completion, SIAM Review 57 (2015), no. 2, 225–251.

[CLM16]

T. Tony Cai, Xiaodong Li, and Zongming Ma, Optimal rates of convergence for noisy sparse
phase retrieval via thresholded Wirtinger flow, The Annals of Statistics 44 (2016), no. 5, 2221–
2251.

[CLS15a]

Emmanuel J. Candès, Xiaodong Li, and Mahdi Soltanolkotabi, Phase retrieval from coded
diffraction patterns, Applied and Computational Harmonic Analysis 39 (2015), no. 2, 277–
299.

55

[CLS15b]

, Phase retrieval via Wirtinger flow: Theory and algorithms, IEEE Trans. Inform.
Theory 61 (2015), no. 4, 1985–2007.

[Cor06]

John V. Corbett, The Pauli problem, state reconstruction and quantum-real numbers, Reports
on Mathematical Physics 57 (2006), no. 1, 53–68.

[CSV13]

Emmanuel J. Candès, Thomas Strohmer, and Vladislav Voroninski, Phaselift: Exact and stable
signal recovery from magnitude measurements via convex programming, Communications on
Pure and Applied Mathematics 66 (2013), no. 8, 1241–1274.

[DJ17]

Laurent Demanet and Vincent Jugnon, Convex recovery from interferometric measurements,
IEEE Trans. Computational Imaging 3 (2017), no. 2, 282–295.

[DJM13]

David L. Donoho, Adel Javanmard, and Andrea Montanari, Information-theoretically optimal compressed sensing via spatial coupling and approximate message passing, IEEE Trans.
Inform. Theory 59 (2013), no. 11, 7434–7464.

[DK70]

Chandler Davis and William M. Kahan, The rotation of eigenvectors by a perturbation. III,
SIAM Journal on Numerical Analysis 7 (1970), no. 1, 1–46.

[DKMZ11] Aurelien Decelle, Florent Krzakala, Cristopher Moore, and Lenka Zdeborová, Asymptotic
analysis of the stochastic block model for modular networks and its algorithmic applications,
Physical Review E 84 (2011), no. 6, 066106.
[DL17]
[DM13]

Oussama Dhifallah and Yue M. Lu, Fundamental limits of PhaseMax for phase retrieval: A
replica analysis, arXiv:1708.03355, 2017.
p
Yash Deshpande and Andrea Montanari, Finding hidden cliques of size N/e in nearly linear
time, Foundations of Computational Mathematics (2013), 1–60.

[DM16]

David L. Donoho and Andrea Montanari, High dimensional robust M-estimation: Asymptotic
variance via approximate message passing, Probability Theory and Related Fields 166 (2016),
no. 3–4, 935–969.

[DMM09]

David L. Donoho, Arian Maleki, and Andrea Montanari, Message Passing Algorithms for
Compressed Sensing, Proceedings of the National Academy of Sciences 106 (2009), 18914–
18919.

[DMM11]

David L. Donoho, Arian Maleki, and Andrea Montanari, The noise-sensitivity phase transition
in compressed sensing, IEEE Trans. Inform. Theory 57 (2011), no. 10, 6920–6941.

[DR17]

John C. Duchi and Feng Ruan, Solving (most) of a set of quadratic equalities: Composite
optimization for robust phase retrieval, arXiv:1705.02356, 2017.

[FD87]

J. R. Fienup and J. C. Dainty, Phase retrieval and image reconstruction for astronomy, Image
Recovery: Theory and Application (1987), 231–275.

[Fie82]

J. R. Fienup, Phase retrieval algorithms: A comparison, Applied Optics 21 (1982), no. 15,
2758–2769.

56

[Ger72]

Ralph W. Gerchberg, A practical algorithm for the determination of the phase from image and
diffraction plane pictures, Optik 35 (1972), 237–246.

[GS16]

Tom Goldstein and Christoph Studer, Phasemax: Convex phase retrieval via basis pursuit,
arXiv:1610.07531, 2016.

[Har93]

Robert W. Harrison, Phase problem in crystallography, J. Optical Soc. America A 10 (1993),
no. 5, 1046–1055.

[HJ12]

Roger A. Horn and Charles R. Johnson, Matrix analysis, Cambridge University Press, 2012.

[JM13]

Adel Javanmard and Andrea Montanari, State evolution for general approximate message passing algorithms, with applications to spatial coupling, Information and Inference (2013), 115–
144.

[JNS13]

Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi, Low-rank matrix completion using alternating minimization, Proc. of the 45th Ann. ACM Symp. on Theory of Computing (STOC)
(Palo Alto, CA), ACM, June 2013, pp. 665–674.

[Kar13]

Noureddine El Karoui, Asymptotic behavior of unregularized and ridge-regularized highdimensional robust regression estimators: Rigorous results, arXiv:1311.2445, 2013.

[KKM+ 16] Yoshiyuki Kabashima, Florent Krzakala, Marc Mézard, Ayaka Sakata, and Lenka Zdeborová,
Phase transitions and sample complexity in bayes-optimal matrix factorization, IEEE Trans.
Inform. Theory 62 (2016), no. 7, 4228–4265.
[KMO10]

Raghunandan H. Keshavan, Andrea Montanari, and Sewoong Oh, Matrix completion from a
few entries, IEEE Trans. Inform. Theory 56 (2010), no. 6, 2980–2998.

[KMZ13]

Florent Krzakala, Marc Mézard, and Lenka Zdeborová, Phase diagram and approximate message passing for blind calibration and dictionary learning, Proc. of the IEEE Int. Symposium
on Inform. Theory (ISIT) (Istanbul, Turkey), July 2013, pp. 659–663.

[LGL15]

Gen Li, Yuantao Gu, and Yue M. Lu, Phase retrieval using iterative projections: Dynamics
in the large systems limit, Proc. of the 53rd Annual Allerton Conf. on Commun., Control, and
Computing (Allerton) (Monticello, IL), Oct. 2015, pp. 1114–1118.

[LL17]

Yue M. Lu and Gen Li, Phase transitions of spectral initialization for high-dimensional nonconvex estimation, arXiv:1702.06435, 2017.

[LLJB17]

Kiryung Lee, Yanjun Li, Marius Junge, and Yoram Bresler, Blind recovery of sparse signals
from subsampled convolution, IEEE Trans. Inform. Theory 63 (2017), no. 2, 802–821.

[LLSW16] Xiaodong Li, Shuyang Ling, Thomas Strohmer, and Ke Wei, Rapid, robust, and reliable blind
deconvolution via nonconvex optimization, arXiv:1606.04933, 2016.
[Mil90]

Rick P. Millane, Phase retrieval in crystallography and optics, J. Optical Soc. America A 7
(1990), no. 3, 394–411.

[MISE08]

Jianwei Miao, Tetsuya Ishikawa, Qun Shen, and Thomas Earnest, Extending X-ray crystallography to allow the imaging of noncrystalline materials, cells, and single protein complexes,
Annu. Rev. Phys. Chem. 59 (2008), 387–410.
57

[MNS14]

Elchanan Mossel, Joe Neeman, and Allan Sly, Belief propagation, robust reconstruction and
optimal recovery of block models, Conference on Learning Theory (COLT) (Barcelona, Spain),
June 2014, pp. 356–370.

[MP67]

Vladimir A. Marčenko and Leonid A. Pastur, Distribution of eigenvalues in certain sets of
random matrices, Mat. Sb. (N.S.) 72 (1967), 457–483 (in Russian).

[MR16]

Andrea Montanari and Emile Richard, Non-negative principal component analysis: Message
passing algorithms and sharp asymptotics, IEEE Trans. Inform. Theory 62 (2016), no. 3,
1458–1484.

[MX16]

Elchanan Mossel and Jiaming Xu, Density evolution in the degree-correlated stochastic block
model, Conference on Learning Theory (COLT) (New York, NY), June 2016, pp. 1319–1356.

[NJS13]

Praneeth Netrapalli, Prateek Jain, and Sujay Sanghavi, Phase retrieval using alternating minimization, Advances in Neural Information Processing Systems, 2013, pp. 2796–2804.

[OTH13]

Samet Oymak, Christos Thrampoulidis, and Babak Hassibi, The squared-error of generalized
LASSO: A precise analysis, Proc. of the 51st Annual Allerton Conf. on Commun., Control, and
Computing (Allerton) (Monticello, IL), Oct. 2013, pp. 1002–1009.

[Ran11]

Sundeep Rangan, Generalized Approximate Message Passing for Estimation with Random Linear Mixing, Proc. of the IEEE Int. Symposium on Inform. Theory (ISIT) (St. Petersburg), Aug.
2011, pp. 2168–2172.

[RG01]

Sundeep Rangan and Vivek K. Goyal, Recursive consistent estimation with bounded noise,
IEEE Trans. Inform. Theory 47 (2001), no. 1, 457–464.

[RP16]

Galen Reeves and Henry D. Pfister, The replica-symmetric prediction for compressed sensing
with Gaussian matrices is exact, Proc. of the IEEE Int. Symposium on Inform. Theory (ISIT)
(Barcelona, Spain), July 2016, pp. 665–669.

[SB10]

Jack W. Silverstein and Zhidong Bai, Spectral Analysis of Large Dimensional Random Matrices (2nd edition), Springer, 2010.

[SC95]

Jack W. Silverstein and Sang-Il Choi, Analysis of the limiting spectral distribution of largedimensional random matrices, Journal of Multivariate Analysis 54 (1995), no. 2, 295–309.

[SC16]

Weijie Su and Emmanuel J. Candès, Slope is adaptive to unknown sparsity and asymptotically
minimax, Annals of Statistics 44 (2016), no. 3, 1038–1068.

[SEC+ 15]

Yoav Shechtman, Yonina C. Eldar, Oren Cohen, Henry N. Chapman, Jianwei Miao, and
Mordechai Segev, Phase retrieval with application to optical imaging: a contemporary
overview, IEEE Signal Processing Magazine 32 (2015), no. 3, 87–109.

[Sol17]

Mahdi Soltanolkotabi, Structured signal recovery from quadratic measurements: Breaking
sample complexity barriers via nonconvex optimization, arXiv:1702.06175, 2017.

[Spe93]

Roland Speicher, Free convolution and the random sum of matrices, Publ. Res. Inst. Math. Sci.
29 (1993), 731–744.

58

[SR15]

Philip Schniter and Sundeep Rangan, Compressive phase retrieval via generalized approximate
message passing, IEEE Transactions on Signal Processing 63 (2015), no. 4, 1043–1055.

[UE88]

Michael Unser and Murray Eden, Maximum likelihood estimation of linear signal parameters
for Poisson processes, IEEE Trans. Acoust., Speech, and Signal Process. 36 (1988), no. 6,
942–945.

[VJ17]

Ramji Venkataramanan and Oliver Johnson, Strong converse bounds for high-dimensional estimation, arXiv:1706.04410, 2017.

[Voi91]

Dan Voiculescu, Limit laws for random matrices and free products, Inventiones Mathematicae
104 (1991), 201–220.

[Wal63]

Adriaan Walther, The question of phase retrieval in optics, Journal of Modern Optics 10 (1963),
no. 1, 41–49.

[WdM15]

Irène Waldspurger, Alexandre d’Aspremont, and Stéphane Mallat, Phase recovery, maxcut and
complex semidefinite programming, Mathematical Programming 149 (2015), no. 1-2, 47–81.

[Wei15]

Ke Wei, Solving systems of phaseless equations via Kaczmarz methods: A proof of concept
study, Inverse Problems 31 (2015), no. 12.

[WG16]

Gang Wang and Georgios B. Giannakis, Solving random systems of quadratic equations via
truncated generalized gradient flow, Advances in Neural Information Processing Systems,
2016, pp. 568–576.

[WGE16]

Gang Wang, Georgios B. Giannakis, and Yonina C. Eldar, Solving systems of random quadratic
equations via truncated amplitude flow, arXiv:1605.08285, 2016.

[WGSC17] Gang Wang, Georgios B. Giannakis, Yousef Saad, and Jie Chen, Solving almost all systems of
random quadratic equations, arXiv:1705.10407, 2017.
[YLSV12]

Feng Yang, Yue M. Lu, Luciano Sbaiz, and Martin Vetterli, Bits from photons: Oversampled
image acquisition using binary Poisson statistics, IEEE Trans. Image Process. 21 (2012), no. 4,
1421–1436.

[ZK16]

Lenka Zdeborová and Florent Krzakala, Statistical physics of inference: Thresholds and algorithms, Advances in Physics 65 (2016), no. 5, 453–552.

[ZL16]

Huishuai Zhang and Yingbin Liang, Reshaped Wirtinger Flow for solving quadratic system of
equations, Advances in Neural Information Processing Systems, 2016, pp. 2622–2630.

59

