Irène Waldspurger

Phase retrieval with random Gaussian sensing vectors by alternating projections

We consider a phase retrieval problem, where we want to reconstruct a n-dimensional vector from its phaseless scalar products with m sensing vectors, independently sampled from complex normal distributions. We show that, with a suitable initialization procedure, the classical algorithm of alternating projections (Gerchberg-Saxton) succeeds with high probability when m ≥ Cn, for some C > 0. We conjecture that this result is still true when no special initialization procedure is used, and present numerical experiments that support this conjecture.

I. INTRODUCTION

A phase retrieval problem consists in recovering an unknown vector x 0 ∈ C n from m phaseless linear measurements, of the form

b k = | a k , x 0 |, k = 1, . . . , m.
Such problems naturally appear in various applications, notably in optical imaging [START_REF] Schechtman | Phase retrieval with application to optical imaging: a contemporary overview[END_REF]. A lot of efforts have thus been made to design efficient algorithms that could numerically solve these problems.

The oldest reconstruction algorithms [START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF], [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF] were iterative: they started from a random initial guess of x 0 , and tried to iteratively refine it by various heuristics. These methods sometimes succeed, but can also fail to converge towards the correct solution: they may get stuck in stagnation points, whose existence is due to the non-convexity of the problem. When they converge and when they do not is not clearly understood.

To overcome convergence problems, convexification methods have been introduced [START_REF] Chai | Array imaging using intensity-only measurements[END_REF], [START_REF] Candès | Phaselift: exact and stable signal recovery from magnitude measurements via convex programming[END_REF]. Their principle is to lift the non-convex problem to a matricial space, where it can be approximated by a convex problem. These methods provably reconstruct the unknown vector x 0 with high probability if the sensing vectors a k are "random enough" [START_REF] Candès | Solving quadratic equations via phaselift when there are about as many equations as unknowns[END_REF]- [START_REF] Gross | A partial derandomization of PhaseLift using spherical designs[END_REF]. Numerical experiments show that they also perform well on more structured, non-random phase retrieval problems [START_REF] Sun | Estimating a signal from a magnitude spectrogram via convex optimization[END_REF], [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF].

Unfortunately, this good precision comes at a high computational cost: optimizing an n × n matrix is much slower than directly reconstructing a n-dimensional vector. More recently, a new family of algorithms has thus been developed, which enjoy similar theoretical guarantees as convexified methods, but have a much smaller computational complexity. The algorithms of this family rely on the following two-step scheme:

(1) an initialization step, that returns a point close to the solution;

CEREMADE, Université Paris Dauphine, and Équipe Mokaplan, Inria Paris; e-mail address: waldspurger@ceremade.dauphine.fr.

(2) a gradient descent (with possible additional refinements) over a well-chosen non-convex cost function. The intuitive reason why this scheme works is that the cost function, although globally non-convex, enjoys some good geometrical property in a neighborhood of the solution (like convexity or a weak form of it [START_REF] White | The local convexity of solving systems of quadratic equations[END_REF]). When the initial point belongs to this neighborhood, the gradient descent then converges to the correct solution.

A preliminary form of this scheme appeared in [START_REF] Netrapalli | Phase retrieval using alternating minimization[END_REF], with an alternating minimization in step (2) instead of a gradient descent. [START_REF] Candès | Phase retrieval via wirtinger flow: theory and algorithms[END_REF] then proved, for a specific cost function, that this two-step scheme was able to exactly reconstruct the unknown x 0 with high probability, in the regime m = O(n log n), if the sensing vectors were independent and Gaussian. In [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], [START_REF] Kolte | Phase retrieval via incremental truncated Wirtinger flow[END_REF], the same result was shown in the regime m = O(n) for a slightly different cost function, with additional truncation steps. In [START_REF] Zhang | Reshaped Wirtinger flow for solving quadratic systems of equations[END_REF], [START_REF] Wang | Solving random systems of quadratic equations via truncated generalized gradient flow[END_REF], it was extended to a different, non-smooth, cost function.

These new methods enjoy much stronger theoretical guarantees than "traditional" algorithms. However, it is not clear whether they really perform better in applications, or whether they actually behave similarly, and are simply easier to theoretically study. Traditional algorithms are well-known, simple and very easy to implement; understanding how they compare to more modern methods is of much value for applications.

In this article, we take a first step towards this goal, by considering the very classical alternating projections algorithm, introduced by Gerchberg and Saxton in [START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF], arguably the most simple and widely used method for phase retrieval. We show that, in the setting where sensing vectors are independent and Gaussian, it performs as well as gradient descent over a suitable cost function: it converges linearly to the true solution with high probability, provided that it is correctly initialized.

Theorem (See Corollary III.8). There exist absolute constants

C 1 , C 2 , M > 0, δ ∈]0; 1[ such that, if m ≥ M n, if
the sensing vectors are independently chosen according to complex normal distributions, and if the alternating projections are carefully initialized, for example with the method described in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], then the sequence of iterates (z t ) t∈N they produce satisfies

∀t ∈ N * , inf φ∈R ||e iφ x 0 -z t || ≤ δ t ||x 0 ||,
with probability at least

1 -C 1 exp(-C 2 m).
Several authors have already tried to establish properties of this kind, but, compared to ours, their results were significantly suboptimal in various respects. Using transversality arguments, [START_REF] Noll | On local convergence of the method of alternating projections[END_REF], [START_REF] Chen | Phase retrieval with one or two diffraction patterns by alternating projection with null initialization[END_REF] have shown the local convergence of alternating projections for relatively general families of sensing vectors. Unfortunately, transversality arguments give no control on the convergence radius of the algorithm, which can be extremely small. Lower bounding the radius requires using the statistical properties of the sensing vectors. This was first attempted in [START_REF] Netrapalli | Phase retrieval using alternating minimization[END_REF]. The result obtained by these authors was very similar to ours, but required resampling the sensing vectors at each step of the algorithm, an operation that is almost never done in practice. For a non resampled version, a preliminary result was given in [START_REF] Soltanolkotabi | Algorithms and theory for clustering and nonconvex quadratic programming[END_REF], but this result did not capture the correct scaling of the convergence radius as a function of m and n. As a consequence, it only established global convergence to the correct solution for a suboptimal number of measurements (m = O(n log 2 n)), and with a complex initialization procedure.

To have theory and practice exactly coincide, the role of the initialization procedure should also be examined: in applications, alternating projections are often used with a random initialization, and not with a carefully chosen one. Our numerical experiments indicate that removing the careful initialization does not significantly alter the convergence of the algorithm (still in the setting of Gaussian independent sensing vectors). This fact is related to the observations of [START_REF] Sun | A geometric analysis of phase retrieval[END_REF], who proved that, at least in the regime m ≥ O(n log 3 n) and for a specific cost function, the initialization part of the two-step scheme is not necessary in order for the algorithm to converge. In the context of alternating projections, we were however not able to prove a similar result 1 , and leave it as a conjecture.

Conjecture (See Conjecture IV.1). Let any > 0 be fixed. When m ≥ Cn, for C > 0 large enough, alternating projections, starting from a random initialization (chosen according to a rotationally invariant distribution), converge to the true solution with probability at least 1 -.

The article is organized as follows. Section II precisely defines phase retrieval problems and the alternating projections algorithm. Section III states and proves the main result: the global convergence of alternating projections, with proper initialization, for m = O(n) independent Gaussian measurements. Section IV contains numerical experiments, and presents our conjecture about the non-necessity of the initialization procedure.

This article only considers the most simple setting, where sensing vectors are independent and Gaussian, and measurements are not noisy. We made this choice in order to keep the technical content simple, but we hope that our results extend to more realistic settings, and future work should examine this issue.

A. Notations

For any z ∈ C, |z| is the modulus of z. We extend this notation to vectors: if z ∈ C k for some k ∈ N * , then |z| is the vector of (R + ) k such that

|z| i = |z i |, ∀i = 1, . . . , k.
For any z ∈ C, we set E phase (z) to be the following subset of C:

E phase (z) = z |z| if z ∈ C -{0}; = {e iφ , φ ∈ R} if z = 0.
We extend this definition to vectors z ∈ C k :

E phase (z) = k i=1 E phase (z i ).
For any z ∈ C, we define phase(z) by

phase(z) = z |z| if z ∈ C -{0}; = 1 if z = 0,
and extend this definition to vectors z ∈ C k , as for the modulus.

We denote by the pointwise product of vectors: for all a, b ∈ C k , (a b) is the vector of C k such that

(a b) i = a i b i , ∀i = 1, . . . , k.
We define the operator norm of any matrix A ∈ C n1×n2 by

|||A||| = sup v∈C n 2 ,||v||=1
||Av||.

We denote by A † its Moore-Penrose pseudo-inverse. We note that AA † is the orthogonal projection onto Range(A).

II. PROBLEM SETUP A. Phase retrieval problem

Les n, m be positive integers. The goal of a phase retrieval problem is to reconstruct an unknown vector x 0 ∈ C n from m measurements with a specific form.

We assume a 1 , . . . , a m ∈ C n are given; they are called the sensing vectors. We define a matrix A ∈ C m×n by

A =    a * 1 . . . a * m    .
This matrix is called the measurement matrix. The associated phase retrieval problem is:

reconstruct x 0 from b def = |Ax 0 |. (1)
As the modulus is invariant to multiplication by unitary complex numbers, we can never hope to reconstruct x 0 better than up to multiplication by a global phase. So, instead of exactly reconstructing x 0 , we want to reconstruct x 1 such that

x 1 = e iφ x 0 , for some φ ∈ R.
In all this article, we assume the sensing vectors to be independent realizations of centered Gaussian variables with identity covariance:

(a i ) j ∼ N 0, 1 2 +N 0, 1 2 i, ∀1 ≤ i ≤ m, 1 ≤ j ≤ n.
(

) 2
The measurement matrix is in particular independent from x 0 .

[23] and [START_REF] Conca | Algebraic characterization of injectivity in phase retrieval[END_REF] have proved that when m ≥ 4n -4, Problem (1) always has a unique solution, up to a global phase, for all measurement matrices A belonging to some open dense subset of C m×n . So when A is chosen according to a probability density that is absolutely continuous with respect to Lebesgue measure, Problem (1) has a unique solution with probability 1 (because the complement of an open dense subset has zero Lebesgue measure). In particular, with our measurement model (2), the reconstruction is guaranteed to be unique, with probability 1, when m ≥ 4n -4.

B. Alternating projections

The alternating projections method has been introduced for phase retrieval problems by [START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF]. It focuses on the reconstruction of Ax 0 ; if A is injective, this then allows to recover x 0 .

To reconstruct Ax 0 , it is enough to find z ∈ C m in the intersection of the following two sets.

(1) z ∈ {z ∈ C m , |z | = b};

(2) z ∈ Range(A). Indeed, when the solution to Problem (1) is unique, Ax 0 is the only element of C m that simultaneously satisfies these two conditions (up to a global phase).

A natural heuristic to find such a z is to pick any initial guess z 0 , then to alternatively project it on the two constraint sets. In this context, we call projection on a closed set E ⊂ C m a function P :

C m → E such that, for any x ∈ C m , ||x -P (x)|| = inf e∈E ||x -e||.
The two sets defining constraints (1) and ( 2) admit projections with simple analytical expressions, which leads to the following formulas:

y k = b phase(y k ); (Projection onto set (1)) (3a) y k+1 = (AA † )y k . (Projection onto set (2)) (3b)
If, for each k, we define z k as the vector such that y k = Az k2 , an equivalent form of these equations is:

z k+1 = A † (b phase(Az k )).
The hope is that the sequence (y k ) k∈N converges towards Ax 0 . Unfortunately, it can get stuck in stagnation points. The following proposition (proven in Appendix A) characterizes these stagnation points.

Proposition II.1. For any y 0 , the sequence (y k ) k∈N is bounded. Any accumulation point y ∞ of (y k ) k∈N satisfies the following property:

∃u ∈ E phase (y ∞ ), (AA † )(b u) = y ∞ .
In particular, if y ∞ has no zero entry,

(AA † )(b phase(y ∞ )) = y ∞ .
Despite the relative simplicity of this characteristic property, it is extremely difficult to exactly compute the stagnation points, determine their attraction basin or avoid them when the algorithm happens to run into them.

III. ALTERNATING PROJECTIONS WITH GOOD

INITIALIZATION

In this section, we prove our main result: for m = O(n) Gaussian independent sensing vectors, the method of alternating projections converges to the correct solution with high probability, if it is carefully initialized.

A. Local convergence of alternating projections

We begin with a key result, that we will need to establish our statement. This result is a local contraction property of the alternating projections operator x → A † (b phase(Ax)).

Theorem III.1. There exist , C 1 , C 2 , M > 0, and δ ∈]0; 1[ such that, if m ≥ M n, then, with probability at least

1 -C 1 exp(-C 2 m),
the realizations of A and b = |Ax 0 | are such that the following property holds: for any x ∈ C n such that

inf φ∈R ||e iφ x 0 -x|| ≤ ||x 0 ||, we have inf φ∈R ||e iφ x 0 -A † (b phase(Ax))|| ≤ δ inf φ∈R ||e iφ x 0 -x||. (4)
Proof. For any x ∈ C n , we can write Ax as

Ax = λ x (Ax 0 ) + µ x v x , (5)
where λ x ∈ C, µ x ∈ R + , and v x ∈ Range(A) is a unitary vector orthogonal to Ax 0 .

Outline of the proof: As we will see, the norm of |||A † ||| can be upper bounded by a number arbitrarily close to 1, so to prove Inequality (4), it is enough to focus on the following quantity:

inf φ∈R ||e iφ Ax 0 -b phase(Ax)|| = inf φ∈R ||e iφ |Ax 0 | phase(Ax 0 ) -|Ax 0 | phase(Ax)||.
As Ax is close to Ax 0 (up to a global phase), we can use a kind of mean value inequality, formalized by Lemma III.2. This will yield the upper bound

2 |Ax 0 | 1 6|µx/λx||v x |≥|Ax0| + 6 5 Im µ x λ x v x phase(Ax 0 ) . ( 6
)
The imaginary part comes from the derivative of the phase. The term involving an indicator function can be thought of as a second order term.

For the first term of Equation ( 6), we will see that, when is small enough, the norm of 6|µ x /λ x |v x is very small compared to the norm of Ax 0 , so the vector 1 6|µx/λx||v x |≥|Ax0| is very sparse. Using this sparsity, we will be able to upper bound the first term in Equation ( 6) by

12η |µ x | |λ x | , (7)
where η > 0 can be arbitrarily close to zero if is small enough.

For the second term, if we imagined that v x was a random vector, chosen uniformly in the sphere independently from Ax 0 (which is of course not true), elementary computations would show that

Im µ x λ x v x phase(Ax 0 ) ≈ 1 √ 2 |µ x | |λ x | .
Making this intuition rigorous, we will show that the second term in ( 6) can be upper bounded by [START_REF] Conca | Algebraic characterization of injectivity in phase retrieval[END_REF] 25

|µ x | |λ x | . ( 8
)
Combining Equations ( 7) and ( 8), we will then upper bound (6) by

δ |µ x | |λ x | for some δ < 1.
From the definition of λ x and µ x , we will see that

|λ x | ≈ 1 and |µ x | inf φ∈R ||e iφ Ax 0 -Ax||, which leads us to inf φ∈R ||e iφ Ax 0 -b phase(Ax)|| δ inf φ∈R ||e iφ Ax 0 -Ax||.
Applying A † on both sides concludes.

Let us now give the details.

Proof of Equation (6): The following lemma is proven in Paragraph B-A.

Lemma III.2. For any z 0 , z ∈ C,

| phase(z 0 + z) -phase(z 0 )| ≤ 2 • 1 |z|≥|z0|/6 + 6 5 Im z z 0 .
From the lemma, for any i = 1, . . . , m,

| phase(λ x )(Ax 0 ) i -(b phase(Ax)) i | = |phase(λ x )(Ax 0 ) i -|Ax 0 | i phase(λ x (Ax 0 ) i + µ x (v x ) i )| = |Ax 0 | i phase(Ax 0 ) i -phase (Ax 0 ) i + µ x λ x (v x ) i ≤ 2.|Ax 0 | i 1 |µx/λx||v x i |≥|Ax0|i/6 + 6 5 |Ax 0 | i Im µx λx v x i (Ax 0 ) i = 2.|Ax 0 | i 1 |µx/λx||v x i |≥|Ax0|i/6 + 6 5 Im µx λx v x i phase((Ax 0 ) i ) .
As a consequence,

|| phase(λ x )(Ax 0 ) -b phase(Ax)|| ≤ 2.|Ax 0 | 1 |µx/λx||v x |≥|Ax0|/6 + 6 5 Im µ x λ x v x phase(Ax 0 ) ≤ 2 |Ax 0 | 1 6|µx/λx||v x |≥|Ax0| + 6 5 Im µ x λ x v x phase(Ax 0 ) . ( 9
)
Upper bound for |µx| |λx| : We use a classical result, that allows us to control the norms of A and A † .

Proposition III. 3 ( [25], Thm II.13). If A is chosen according to Equation (2), then, for any t, with probability at least

1 -2 exp -mt 2 , we have, for any x ∈ C n , √ m 1 - n m -t ||x|| ≤ ||Ax||; ||Ax|| ≤ √ m 1 + n m + t ||x||. Let us set x = inf φ∈R ||e iφ x0-x|| ||x0||
≤ . We have

x |||A||| ||x 0 || ≥ inf φ∈R ||e iφ Ax 0 -Ax|| = (1 -|λ x |) 2 ||Ax 0 || 2 + |µ x | 2 .
For any two numbers α and β, α 2 + β 2 ≥ max(α, β). Hence, the previous inequality implies

|λ x | ≥ 1 - x |||A||| ||x 0 || ||Ax 0 || and |µ x | ≤ x |||A||| ||x 0 ||.
From there, we get, upper bounding |||A||| and lower bounding ||Ax 0 || according to Proposition III.3,

|µ x | |λ x | ≤ x |||A||| ||x 0 || 1 - x |||A||| ||x0|| ||Ax0|| ≤ x √ m 1 + n/m + t ||x 0 || 1 - x √ m 1+ √ n/m+t √ m 1- √ n/m-t ≤ x √ m ||x 0 || 1 -n/m -t -x 1 + n/m + t .
So when m ≥ 20000n and ≤ 10 -3 , setting t = 10 -3 , we have

|µ x | |λ x | ≤ 1.01 √ m x ||x 0 ||, (10)
with probability at least 1 -2 exp(-C 2 m), for C 2 = t 2 = 10 -6 .

Upper bounds for the two terms of Equation (9) (corresponds to Equations We define γ > 0 as in Lemma III.4. Using Lemmas III.4 (applied to v = 6(µ x /λ x )v x ) and III.5, when m ≥ M n for M large enough, we can upper bound Equation ( 9) by

|| phase(λ x )(Ax 0 ) -b phase(Ax)|| ≤ 12η + 24 25
µ x λ x ≤ 0.98 × 1.01 √ m x ||x 0 || ≤ 0.99 √ m x ||x 0 ||.
This holds uniformly over all x such that inf φ∈R ||e iφ x 0 -x|| ≤ ||x 0 ||, with probability at least

1 -C 1 exp(-C 2 m). (12)
(From Equation ( 10) and Proposition III.3, still assuming m ≥ 20000n, ≤ 10 -3 and t = 10 -3 , the condition ||v|| < γ||Ax 0 || in Lemma III.4 is satisfied when the event described in Proposition III.3 is realized.) Conclusion: From the last inequality, we have

inf φ∈R ||e iφ x 0 -A † (b phase(Ax))|| ≤ || phase(λ x )x 0 -A † (b phase(Ax))|| ≤ |||A † ||| || phase(λ x )(Ax 0 ) -b phase(Ax)|| ≤ 0.99|||A † ||| √ m x ||x 0 || = 0.99|||A † ||| √ m inf φ∈R ||e iφ x 0 -x||.
When the event of Proposition III.3 is realized, still for t = 10 -3 ,

|||A † ||| ≤ 1 √ m(0.999 -n/m) ≤ 0.999 0.99 √ m .
So, with the probability of Equation ( 12

B. Global convergence

In the last paragraph, we have seen that the alternating projections operator is contractive, with high probability, in an ||x 0 ||-neighborhood of the solution x 0 . This implies that, if the starting point of alternating projections is at distance at most ||x 0 || from x 0 , alternating projections converge to x 0 . So if we have a way to find such an initial point, we obtain a globally convergent algorithm.

Several initialization methods have been proposed that achieve the precision we need with an optimal number of measurements, that is m = O(n). Let us mention the truncated spectral initialization by [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] (improving upon the slightly suboptimal spectral initializations introduced by [START_REF] Netrapalli | Phase retrieval using alternating minimization[END_REF] and [START_REF] Candès | Phase retrieval via wirtinger flow: theory and algorithms[END_REF]), the null initialization by [START_REF] Chen | Phase retrieval with one or two diffraction patterns by alternating projection with null initialization[END_REF] and the method described by [START_REF] Gao | Gauss-Newton method for phase retrieval[END_REF]. All these methods consist in computing the largest or smallest eigenvector of

m i=1 α i a i a * i ,
where the α 1 , . . . , α m are carefully chosen coefficients, that depend only on b.

The method of [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], for example, has the following guarantees.

Theorem III.7 (Proposition 3 of [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]). Let > 0 be fixed.

There exist C 1 , C 2 , M > 0 such that, when m ≥ M n, with probability at least

1 -C 1 exp(-C 2 m),
the realizations of A and b = |Ax 0 | are such that the main eigenvector of

1 m m i=1 b 2 i 1 b 2 i ≤ 9 m m j=1 b 2 i a i a * i , (13)
that we denote by z, obeys

inf φ∈R,λ∈R * + ||e iφ x 0 -λz|| ≤ ||x 0 ||.
Combining this initialization procedure with alternating projections, we get Algorithm 1. As shown by the following corollary, it converges towards the correct solution, at a linear rate, with high probability, for m = O(n).

Input : A ∈ C m×n , b = |Ax 0 | ∈ R m , T ∈ N * .
Initialization: set z 0 to be the main eigenvector of the matrix in Equation [START_REF] Candès | Phase retrieval via wirtinger flow: theory and algorithms[END_REF].

∀t ≤ T, inf φ∈R ||e iφ x 0 -z t || ≤ δ t ||x 0 ||, (14)
provided that m ≥ M n.

Proof. Let us fix , δ ∈]0; 1[ as in Theorem III.1. Let us assume that the properties described in Theorems III.1 and III.7 hold; it happens on an event of probability at least

1 -C 1 exp(-C 2 m),
provided that m ≥ M n, for some constants C 1 , C 2 , M > 0.

Let us prove that, on this event, Equation ( 14) also holds. We proceed by recursion. From Theorem III.7, there exist

φ ∈ R, λ ∈ R * + such that ||e iφ x 0 -λz 0 || ≤ ||x 0 ||.
So, from Theorem III.1, applied to x = λz 0 ,

inf φ∈R ||e iφ x 0 -z 1 || = inf φ∈R ||e iφ x 0 -A † (b phase(z 0 ))|| = inf φ∈R ||e iφ x 0 -A † (b phase(λz 0 ))|| ≤ δ inf φ∈R ||e iφ x 0 -λz 0 || ≤ δ||x 0 ||.
This proves Equation ( 14) for t = 1. Let us stress that it is critical here that the property of Theorem III.1 holds uniformly over all x that are ||x 0 ||-close to x 0 : if this property was only true for vectors x independent from A and b, we could not apply it to λz 0 , since λz 0 depends on A and b through the initialization procedure.

The same reasoning can be reapplied to also prove the equation for t = 2, 3, . . . . Remark III.9. From our computations based on the proof in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], explicit values for the constants in Theorem III.7 can be chosen as follows: These theoretical values are absurdly high. They do not match our numerical experiments (Section IV), which suggest that Algorithm 1 succeeds with probability close to 1 as soon as m ≥ 4n.

C 1 = 6;
For Theorem III.7, much tighter constants could probably be deduced from the recent work [START_REF] Lu | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF]. However, for Corollary III.8, it does not seem possible to make the constants close to tight without significantly improving the proof technique.

One particular source of non-tightness in the current technique is that several inequalities that are required hold "uniformly over some large set" (the property in Theorem III.1, for example, holds "for all x ∈ C n that are -close to x 0 "). This requirement is a priori too strong: it is actually sufficient that the inequalities hold for "all the points in the large set that are visited by the algorithm". However, since the points visited by the algorithm do not have a simple expression, it is unclear how to prove only the weaker version.

C. Complexity

Let η > 0 be the relative precision that we want to achieve:

inf φ∈R ||e iφ x 0 -z T || ≤ η||x 0 ||.
Let us compute the number of operations that Algorithm 1 requires to reach this precision.

The main eigenvector of the matrix defined in Equation ( 13) can be computed -up to precision η -in approximately O(log(1/η) + log(n)) power iterations. Each power iteration is essentially a matrix-vector multiplication, and thus requires O(mn) operations. 3 As a consequence, the complexity of the initialization is

O(mn (log(1/η) + log(n))).
Then, at each step of the for loop, the most costly operation is the multiplication by A † . When performed with the conjugate gradient method, it requires O(mn log(1/η)) operations. To reach a precision equal to η, we need to perform O(log(1/η)) iterations of the loop. So the total complexity of Algorithm 1 is O(mn log 2 (1/η) + log(n) ).

As a comparison, Truncated Wirtinger flow, which is currently the most efficient known method for phase retrieval from Gaussian measurements, has an identical complexity, up to a log(1/η) factor (see Figure 1).

Let us note that, when A has some special structure (it encodes Fourier masks [START_REF] Candès | Phase retrieval from coded diffraction patterns[END_REF], for example), the complexity of the algorithm can be further reduced, and exactly matches the one of Truncated Wirtinger flow.

IV. NUMERICAL EXPERIMENTS

In this section, we numerically validate Corollary III.8. We formulate a conjecture about the convergence of alternating projections with random initialization, in the regime m = O(n).

The code used to generate Figures 2, 3 and 4 is available at http://www-math.mit.edu/ ∼ waldspur/code/alternating projections code.zip. 3 These matrix-vector multiplications can be computed without forming the whole matrix (which would require O(mn 2 ) operations), because this matrix factorizes as

1 m A * Diag(b 2 I)A,
where

I ∈ R m is such that ∀i ≤ m, I i = 1 b 2 i ≤ 9

A. Alternating projections with careful initialization

Corollary III.8 states that alternating projections succeed with high probability, when they start from a good initial point, in the regime where the number of measurements is linear in the problem dimension (m = O(n)).

We use the initialization method described in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], as presented in Algorithm 1. We run the algorithm many times, for various choices of n and m. This allows us to compute an empirical probability of success, for each value of (n, m).

The results are presented in Figure 2. They confirm that, when m = Cn, for a sufficiently large constant C > 0, the success probability can be arbitrarily close to 1.

B. Alternating projections without careful initialization

In a second experiment, we measure the probability that alternating projections succeed, when started from a random initial point (sampled from the unit sphere with uniform probability).

As previously, we numerically compute this probability for various pairs (m, n). The results, presented in Figure 3 seem to indicate that alternating projections behave similarly with and without a careful initialization procedure: they converge towards the correct solution as soon as m ≥ Cn, for a suitable constant C > 0. Only the value of C changes. Figure 4 displays similar computations, under a different form: it shows how many measurements m are needed so that alternating projections, with careful or with random initialization, succeeds with probability 1/2. This number apparently grows linearly with n whatever the initialization method.

As a consequence, we have the following conjecture.

Conjecture IV.1. Let any > 0 be fixed. When m ≥ Cn, for C > 0 large enough, alternating projections with a random rotationally invariant initialization succeed with probability at least 1 -.

This conjecture is in the same spirit as several recent works, which show that non-convex iterative algorithms, starting from random initial points, succeed with high probability for a number of problems. These works cover matrix sensing [START_REF] Bhojanapalli | Global optimality of local search for low rank matrix recovery[END_REF], matrix completion [START_REF] Ge | Matrix completion has no spurious local minimum[END_REF], robust PCA [START_REF] Ge | No spurious local minima in nonconvex low rank problems: a unified geometric analysis[END_REF], and synchronization [START_REF]Nonconvex phase synchronization[END_REF]. The closest one to our conjecture is [START_REF] Sun | A geometric analysis of phase retrieval[END_REF], that also considers phase retrieval with Gaussian sensing vectors: in the almost optimal regime m = O(n log 3 n), it shows that gradient descent over a (specific) non-convex function succeeds with high probability, even when randomly initialized.

Unfortunately, the proof method of these works does not seem to adapt our case. Indeed, it consists in showing that the iterative routine has no attractive stagnation point. But in our case, complementary experiments (not shown here), suggest that, in the regime m = O(n), there are attractive stagnation points. These points seem to have small attraction basins, so they are not necessarily an obstacle for the success of the algorithm with high probability, but rigorously establishing this property seems difficult.

APPENDIX A PROOF OF PROPOSITION II.1 Proposition (Proposition II.1). For any y 0 , the sequence (y k ) k∈N is bounded. Any accumulation point y ∞ of (y k ) k∈N satisfies the following property:

∃u ∈ E phase (y ∞ ), (AA † )(b u) = y ∞ .
In particular, if y ∞ has no zero entry,

(AA † )(b phase(y ∞ )) = y ∞ .
Proof of Proposition II.1. The boundedness of (y k ) k∈N is a consequence of the fact that ||y k || = ||b|| for all k, so

||y k+1 || ≤ |||AA † ||| ||b||.
Let us show the second part of the statement. Let y ∞ be an accumulation point of (y k ) k∈N , and let φ : N → N be an extraction such that

y φ(n) → y ∞ when n → +∞.
By compacity, as (y φ(n) ) n∈N and (y φ(n)+1 ) n∈N are bounded sequences, we can assume, even if we have to consider replace φ by a subextraction, that they also converge. We denote by y ∞ and y +1

∞ their limits:

y φ(n) → y ∞ and y φ(n)+1 → y +1 ∞ when n → +∞.
Dimension n Number of measurements m Figure 4: As a function of n, the smallest value of m for which alternating projections succeeds with probability 50%. The solid line is when the algorithm is initialized with the method in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]; the dashed line is with random Gaussian initialization.

Let us define

E b = {y ∈ C m , |y | = b}.
We observe that, for any k,

d(y k-1 , Range(A)) ≥ d(y k , E b ) ≥ d(y k , Range(A)).
Indeed, because the operators y → b phase(y) and y → (AA † )y are projections,

d(y k-1 , Range(A)) = d(y k-1 , y k ) ≥ d(y k , E b ); d(y k , E b ) = d(y k , y k ) ≥ d(y k , Range(A)).
So the sequences (d(y k , E b )) k∈N and (d(y k , Range(A))) k∈N converge to the same non-negative limit, that we denote by δ.

In particular,

d(y ∞ , E b ) = δ = d(y ∞ , Range(A)).
If we pass to the limit the equalities d(y

φ(n) , E b ) = ||y φ(n) - y φ(n) || and d(y φ(n) , Range(A)) = ||y φ(n) -y φ(n)+1 ||, we get ||y ∞ -y ∞ || = ||y ∞ -y +1 ∞ || = δ = d(y ∞ , Range(A)
). As Range(A) is convex, the projection of y ∞ onto it is uniquely defined. This implies

y ∞ = y +1
∞ , and, because ∀n, y φ(n)+1 = (AA † )y φ(n) ,

y ∞ = y +1 ∞ = (AA † )y ∞ .
To conclude, we now have to show that y ∞ = b u for some u ∈ E phase (y ∞ ). We use the fact that, for all n, y φ(n) = b phase(y φ(n) ).

For any i ∈ {1, . . . , m}, if (y ∞ ) i = 0, phase is continuous around (y ∞ ) i , so (y ∞ ) i = b i phase((y ∞ ) i ). We then set u i = phase((y ∞ ) i ), and we have (y

∞ ) i = b i u i . If (y ∞ ) i = 0, we set u i = phase((y ∞ ) i ) ∈ E phase (0) = E phase ((y ∞ ) i ). We then have y ∞ = |y ∞ |u i = b i u i .
With this definition of u, we have, as claimed, y ∞ = b u and u ∈ E phase (y ∞ ).

APPENDIX B TECHNICAL LEMMAS FOR SECTION III

A. Proof of Lemma III.2 Lemma (Lemma III.2). For any z 0 , z ∈ C,

| phase(z 0 + z) -phase(z 0 )| ≤ 2 • 1 |z|≥|z0|/6 + 6 5 Im z z 0 .
Proof. The inequality holds if z 0 = 0, so we can assume z 0 = 0. We remark that, in this case,

| phase(z 0 + z) -phase(z 0 )| = | phase(1 + z/z 0 ) -1|.
It is thus enough to prove the lemma for z 0 = 1, so we make this assumption. When |z| ≥ 1/6, the inequality is valid. Let us now assume that |z| < 1/6. Let θ ∈ -π 2 ; π 2 be such that e iθ = phase(1 + z).

Then | phase(1 + z) -1| = |e iθ -1| = 2| sin(θ/2)| ≤ | tan θ| = |Im (1 + z)| |Re (1 + z)| ≤ |Im (z)| 1 -|z| ≤ 6 5 |Im (z)|.
So the inequality is also valid. Proof. For any S ⊂ {1, . . . , m}, we denote by 1 S the vector of C m such that

(1 S ) j = 1 if j ∈ S = 0 if j / ∈ S.
We use the following two lemmas, proven in Paragraphs B-B1 and B-B2.

Let β > 0 be such that 10 β log(1/β) ≤ η (we can for instance choose β = min 0.01, η 2 300 log(1/η) ). Let M be larger than 3 β , as in Lemma B.2. We set γ = β 3/2 e -1/2 .

We assume that Equations ( 15) and ( 16) hold; from the lemmas, this occurs with probability at least

1 -6 exp(-β 3 m/e),
provided that m ≥ M n. (We have used the fact that, for β ≤ 0.01, β 3 /e ≤ 3β log(1/β).) On this event, for any v ∈ Range(A) such that ||v|| < γ||Ax 0 ||, if we set

S v = {i s.t. |v i | ≥ |Ax 0 | i }, we have that Card S v < βm.
Indeed, if it was not the case, we would have, by Equation [START_REF] Kolte | Phase retrieval via incremental truncated Wirtinger flow[END_REF],

||v|| ≥ ||v 1 Sv || ≥ || |Ax 0 | 1 Sv || ≥ β 3/2 e -1/2 ||Ax 0 || = γ||Ax 0 ||,
which is in contradiction with the way we have chosen v.

So we can apply Equation ( 16), and we get

|| |Ax 0 | 1 |v|≥|Ax0| || ≤ ||v 1 S || ≤ 10 β log(1/β)||v|| ≤ η||v||.
Remark B.3. We can choose the values of C 1 , C 2 , M, γ as follows:

C 1 = 6; C 2 = 10 -6 e min 1, η 6 27 log 3 (1/η) ; M = 300 min 1, η 2 3 log(1/η) ; γ = 10 -3 e -1/2 min 1, η 3
3 3/2 log 3/2 (1/η) .
1) Proof of Lemma B.1:

Proof of Lemma B.1. We first assume S fixed, with cardinality Card S ≥ βm. We use the following lemma.

Lemma B.4 ( [32], Lemma 2.2). Let k 1 < k 2 be natural numbers. Let X ∈ C k2 be a random vector whose coordinates are independent, Gaussian, of variance 1. Let Y be the projection of X onto its k 1 first coordinates. Then, for any t > 0,

Proba ||Y || ||X|| ≤ tk 1 k 2 ≤ exp (k 1 (1 -t + log t)) if t < 1; Proba ||Y || ||X|| ≥ tk 1 k 2 ≤ exp (k 1 (1 -t + log t)) if t > 1.
From this lemma, for any t ∈]0; 1[, because Ax 0 has independent Gaussian coordinates,

P || |Ax 0 | 1 S || ||Ax 0 || ≤ tβ ≤ exp (-βm(t -1 -ln t)) .
In particular, for t = β 2 e ,

P || |Ax 0 | 1 S || ||Ax 0 || ≤ β 3/2 e -1/2 ≤ exp -βm β 2 e -2 ln β . (17)
As soon as m ≥ 1 2β , the number of subsets S of {1, . . . , m} with cardinality βm satisfies

m βm ≤ em βm βm ≤ exp 2mβ log 1 β . (18)
(The first inequality is a classical result regarding binomial coefficients.)

We combine Equations ( 17) and ( 18): Property ( 15) is satisfied for any S of cardinality βm with probability at least

1 -exp - β 3 e m ≥ 1 -2 exp - β 3 e m , provided that m ≥ 1 2β . When m < 1 2β , since β < 1/2, 1 -2 exp(-β 3 m/e) < 1 -2 exp(-β 2 /(2e)) < 1 -2 exp(-1/(8e))
< 0, so the property also holds with probability larger than 1 -2 exp -β 3 m/e . If it is satisfied for any S of cardinality βm , then it is satisfied for any S of cardinality larger than βm, which implies the result.

2) Proof of Lemma B.2:

Proof of Lemma B.2. We first assume S to be fixed, of cardinality exactly βm .

Any vector y ∈ Range(A) is of the form y = Av, for some v ∈ C n . Inequality ( 16) can then be rewritten as:

||A S v|| = ||Diag(1 S )Av|| ≤ 10 β log(1/β)||Av||, (19)
where A S , by definition, is the submatrix obtained from A by extracting the rows whose indexes are in S.

We apply Proposition III.3 to A and A S , respectively for t = 1 2 and t = 3 log(1/β). It guarantees that the following properties hold:

inf v∈C n ||Av|| ||v|| ≥ √ m 1 2 - n m ; sup v∈C n ||A S v|| ||v|| ≤ √ Card S 1 + n Card S + 3 log(1/β) ,
with respective probabilities at least

1 -2 exp - m 4 ; and 1 -2 exp (-9(Card S) log(1/β)) ≥ 1 -2 exp (-9β log(1/β)m) .
Assuming m ≥ M n for some M > 0, we deduce from these inequalities that, for all v ∈ C n ,

||A S v|| ≤ Card S m 1 + n Card S + 3 log(1/β) 1 2 -n m ||Av|| ≤ β + 1 m   1 + 1 βM + 3 log(1/β) 1 2 - 1 M   ||Av||, (20)
with probability at least

1 -2 exp (-9β log(1/β)m) -2 exp - m 4 .
If we choose M ≥ 3 β , remember that β ≤ 0.01 and m ≥ M n ≥ M , we can deduce from Equation (20) that

||A S v|| ≤ β 1 + 1 3   1 + 1 3 + 3 log(1/β) 1 2 - 1 300   ||Av|| ≤ 10 β log(1/β)||Av||,
which is Equation [START_REF] Chen | Phase retrieval with one or two diffraction patterns by alternating projection with null initialization[END_REF].

As in the proof of Lemma B. This ends the proof. Indeed, if Equation ( 19) holds for any set of cardinality βm , it also holds for any set of cardinality Card S < βm, because ||A S v|| ≤ ||A S v|| whenever S ⊂ S. This implies Equation [START_REF] Zhang | Reshaped Wirtinger flow for solving quadratic systems of equations[END_REF]. We take t = m n-1 (0.04) 2 (which is larger than 1 when m ≥ M n with M > 625), and it implies that provided that m ≥ 2500n. (We have used the fact that, when M ≥ 2500, then t ≥ 4, and when t ≥ 4, then t -1 -log t ≥ t/4.) Second, as A 2:n is a random matrix of size m × (n -1), whose entries are independent and distributed according to the law N (0, 1/2)+N (0, 1/2)i, we deduce from Proposition III.3 applied with t = 0.01 that, with probability at least

provided that m ≥ 20000n.

When Equations ( 24) and ( 25) are simultaneously valid, any w = A 2:n w belonging to Range(A 2:n ) satisfies:

We now conclude. Equations ( 23), ( 24) and ( 25) hold simultaneously with probability at least 1 -4 exp -5.10 -5 m , provided that m ≥ 20000n. Let us show that, on this event, Equation ( 21) also holds. Any v ∈ Range(A)∩{Ax 0 } ⊥ , from Equality [START_REF] Waldspurger | Phase retrieval with random gaussian sensing vectors by alternating projections[END_REF], can be written as

v = w - w, A 1 ||A 1 || 2 A 1 ,
for some w ∈ Range(A 2:n ). Using Equation ( 23), then Equation ( 26), we get:

5 ||v|| ( 11 )

 511 [START_REF] Candès | Phase retrieval from coded diffraction patterns[END_REF] and[START_REF] Gross | A partial derandomization of PhaseLift using spherical designs[END_REF] in the outline): Each term is handled in a separate lemma. The first lemma is proved in Paragraph B-B, the second one in Paragraph B-C. Lemma III.4. For any η > 0, there exists C 1 , C 2 , M, γ > 0 such that the inequality || |Ax 0 | 1 |v|≥|Ax0| || ≤ η||v|| holds for any v ∈ Range(A) such that ||v|| < γ||Ax 0 ||, with probability at least 1 -C 1 exp(-C 2 m), when m ≥ M n. Lemma III.5. When m ≥ 20000n, the property ||Im (v phase(Ax 0 ))|| ≤ 4 holds for any v ∈ Range(A) ∩ {Ax 0 } ⊥ , with probability at least 1 -4 exp -5.10 -5 m . Upper bound for Equation (9): Let us set η = 1

  ), the following property holds uniformly over all x satisfying inf φ∈R ||e iφ x 0 -x|| ≤ ||x 0 ||:inf φ∈R ||e iφ x 0 -A † (b phase(Ax))|| ≤ 0.999 inf φ∈R ||e iφ x 0 -x||.Remark III.6. From the proof (and Remark B.3 in appendix, evaluated at η = 1 600 ), we can deduce explicit values for the constants , C 1 , C 2 , M in Theorem III.1: δ = 0.999; = 10 -3 ; M = 3.10 9 ; C 1 = 12; C 2 = 10 -21 .

  for t = 1 to T do Set z t ← A † (b phase(Az t-1 )). end Output: z T . Algorithm 1: Alternating projections with truncated spectral initialization Corollary III.8. There exist C 1 , C 2 , M > 0, δ ∈]0; 1[ such that, with probability at least 1 -C 1 exp(-C 2 m), the realizations of A and b = |Ax 0 | are such that, for any T ∈ N * , Algorithm 1 satisfies

C 2 = 5 . 10 -12 4 ;M = 4 . 10 12 - 4 .

 251044104 Combined with Remark III.6, this guarantees that Corollary III.8 holds for δ = 0.999;C 1 = 18;C 2 = 5.10 -24 ; M = 4.10 24 .

Figure 1 :Figure 2 :

 12 Figure 1: Complexity of alternating projections with initialization, and truncated Wirtinger flow.

Figure 3 :

 3 Figure 3: Probability of success for alternating projections with a random Gaussian initialization, as a function of n and m. Black points indicate a probability equal to 0, and white points a probability equal to 1. The line serves as a reference: it is the line m = 4n.

B.

  Proof of Lemma III.4 Lemma (Lemma III.4). For any η > 0, there exists C 1 , C 2 , M, γ > 0 such that the inequality || |Ax 0 | 1 |v|≥|Ax0| || ≤ η||v|| holds for any v ∈ Range(A) such that ||v|| < γ||Ax 0 ||, with probability at least 1 -C 1 exp(-C 2 m), when m ≥ M n.

Lemma B. 1 . 15 )Lemma B. 2 .

 1152 Let β ∈]0; 1/2[ be fixed. With probability at least 1 -2 exp(-β 3 m/e), the following property holds: for any S ⊂ {1, . . . , m} such that Card(S) ≥ βm, || |Ax 0 | 1 S || ≥ β 3/2 e -1/2 ||Ax 0 ||. (Let β ∈ 0; 1 100 be fixed. For any M ≥ 3 β , if m ≥ M n, then, with probability at least 1 -4 exp -3β log 1 β m , the following property holds: for any S ⊂ {1, . . . , m} such that Card(S) < βm and for any y ∈ Range(A), ||y 1 S || ≤ 10 β log(1/β)||y||.

C 5 ||v|| ( 21 ). 1 -

 5211 . Proof of Lemma III.5 Lemma (Lemma III.5). When m ≥ 20000n, the property ||Im (v phase(Ax 0 ))|| ≤ 4 holds for any v ∈ Range(A) ∩ {Ax 0 } ⊥ , with probability at least 1 -4 exp -5.10 -5 m .Proof. If we multiply x 0 by a positive real number, we can assume ||x 0 || = 1. Moreover, as the law of A is invariant under right multiplication by a unitary matrix, we can assume that Then, if we write A 1 the first column of A, and A 2:n the submatrix of A obtained by removing this first column,Range(A) ∩ {Ax 0 } ⊥ = w -w, A 1 ||A 1 || 2 A 1 , w ∈ Range(A 2:n ) . (22)We first observe thatsup w∈Range(A2:n)-{0} | w, A 1 | ||w||is the norm of the orthogonal projection of A 1 onto Range(A 2:n ). The (n-1)-dimensional subspace Range(A 2:n ) has a rotationally invariant distribution in C m , and is independent of A 1 . Thus, from Lemma B.4 coming from[START_REF] Dasgupta | An elementary proof of a theorem of Johnson and Lindenstrauss[END_REF], for any t exp (-(n -1)(t -1 -log t)) .

1 -

 1 exp(-(n -1)(t -1 -log t)) ≥ 1 -exp(-(n -1)t/4) = 1 -exp(-4.10 -4 m).

1 - 2 1 -

 21 exp -10 -4 m , we have, for any x ∈ C n-1 ,||A 2:n x|| ≥ √ m that m ≥ 10000n. We now set C = Diag(phase(A 1 ))A 2:n .The matrix Im C Re C has size m × (2(n -1)); its entries are independent and distributed according to the law N (0, 1/2). So by[START_REF] Davidson | Local operator theory, random matrices and Banach spaces[END_REF] Thm II.13] (applied with t = 0.01), with probability at least 1 -exp(-5.10 -5 m), we have, for any x ∈ R2(n-1) ,Im C Re C x ≤ m 2 1 + 2(n -1) m + 0.01 ||x|| ≤ 1.02 m 2 ||x||,

≤ 1

 1 Im (w phase(Ax0 )) = ||Im (Cw )|| = Im C Re C Re w Im w

+ w, A 1 ||A 1 || 2 A 1 ≤ 2 ||A 1 || 2 ≥ ( 1 -

 11212121 Im (v phase(Ax 0 )) ≤ Im (w phase(Ax 0 )) Im (w phase(Ax 0 ))+ 0.04||w|| ≤ 0.79||w||.But then, by Equation (23) again,||v|| 2 = ||w|| 2 -w, A 1 (0.04) 2 )||w|| 2 .SoIm (v phase(Ax 0 )) ≤ 0

The proof technique used by[START_REF] Sun | A geometric analysis of phase retrieval[END_REF] collapses when gradient descent is replaced with alternating projections, as detailed in a companion study of this article, available on arXiv[START_REF] Waldspurger | Phase retrieval with random gaussian sensing vectors by alternating projections[END_REF].

which exists and is unique, because y k belongs to Range(A) and A is injective with probability 1