Outlier-Robust Phase Retrieval in Nearly-Linear Time

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 06:11:07

Abstract: Phase retrieval is a fundamental problem in signal processing, where the goal is to recover a (complex-valued) signal from phaseless intensity measurements. It is well-known that natural nonconvex formulations of phase retrieval do not have spurious local optima. However, the theoretical analyses of such landscape results often rely on strong assumptions, such as the sampled measurements are (complex) Gaussian. In this paper, we propose and study the problem of outlier robust phase retrieval. We focus on the real-valued case, where we seek to recover a vector $x \in \mathbb{R}^d$ from $n$ intensity measurements $y_i = (a_i^\top x )^2$, under the assumption that the $a_i$'s are initially i.i.d. Gaussian but a small fraction of the $(y_i, a_i)$ pairs are adversarially corrupted. Our main result is a near sample-optimal nearly-linear time algorithm that provably recovers the ground-truth vector $x$ in the presence of outliers. We first solve a lightweight convex program to find a vector close to the ground truth. We then run robust gradient descent starting from this initial solution, leveraging recent advances in high-dimensional robust statistics. Our approach is conceptually simple and provides a framework for developing robust algorithms for tractable nonconvex problems.

Main paper: 
 introduction: 
 Phase retrieval is a fundamental problem in signal processing with applications in various fields, including electron microscopy [32] , crystallography [33, 36] , astronomy [11] , and optical imaging [37] .In these applications, one often has access to only the magnitudes of the Fourier transforms of a complex signal. This is because measuring magnitude (e.g., by aggregating energy over time) is much easier than measuring phase (which requires detecting rapid changes). We refer the reader to the survey articles [37, 26] for more details about the theory and applications of phase retrieval.In this paper, we focus on the real-valued generalized phase retrieval problem, where the Fourier transform is replaced by a general linear operator. We first give a formal definition of this problem. Definition 1.1 (Phase Retrieval). Let x ∈ R d be the ground-truth vector. Let a 1 . . . a n ∈ R d be n sampling vectors and let y i = ⟨a i , x⟩ 2 ∈ R be the corresponding intensity measurements. Given (a i , y i ) n i=1 as input, the task is to recover x.Note that it is impossible to distinguish between x and -x, so it is sufficient to recover either one.Under certain assumptions (e.g., when the a i 's are Gaussian distributed), the phase retrieval problem in Definition 1.1 can be solved in polynomial time with provable recovery guarantees. This was first achieved via approaches based on semidefinite programming (SDP) relaxations (see, e.g., Candès et al. [5] ). In practice, the problem is often solved using first-order optimization algorithms such as gradient descent. It is well-established that, although many natural formulations of phase retrieval have nonconvex objectives, all local optima are globally optimal under certain assumptions [34, 3, 40] .An example of such objective function is the following:minimize f (z) = n i=1 (y i -⟨a i , z⟩ 2 ) 2 subject to z ∈ R d .Submitted to 38th Conference on Neural Information Processing Systems (NeurIPS 2024). Do not distribute.However, existing analyses of such landscape results often rely on strong assumptions, such as the sampling vectors a i 's are i.i.d. Gaussian. Our work is motivated by the following questions: Can we relax the assumptions used in proving landscape results in many tractable nonconvex problems? In the context of phase retrieval, what happens if a small fraction of the (a i , y i )'s are changed adversarially?We focus on the following strong contamination model (see, e.g., [13] ).). An algorithm first specifies the number of samples n, and n samples are drawn independently from some unknown distribution D. The adversary is allowed to replace up to ϵn samples with arbitrary points. The modified set of n samples is then given to the algorithm as input. We say that a set of samples is ϵ-corrupted if it is generated by the above process. 1Under the ϵ-corruption model for high-dimensional data, a common goal is to design efficient algorithms that can achieve dimension-independent error guarantees. Early work in robust statistics [42, 23, 25] provided sample-efficient estimators for various tasks, but with runtimes exponential in the dimension. A recent line of work, initiated by [13, 28] , has developed computationally efficient robust algorithms for many fundamental high-dimensional tasks. There has been significant progress in the algorithmic aspects of robust high-dimensional statistics (see, e.g., [12] ).We now formally define the main problem that we pose and study in this paper.Problem 1.3 (Outlier-Robust Phase Retrieval). Let ϵ > 0. Let x ∈ R d be the ground-truth vector with ∥x∥ 2 = 1. First, n sampling vectors (a i ) n i=1 are drawn i.i.d. from N (0, I) ∈ R d . Let y i = ⟨a i , x⟩ 2 be the corresponding intensity measurements. Then, an adversary arbitrarily corrupts an ϵ-fraction of the (a i , y i )'s. Finally, the corrupted (a i , y i )'s are given to the algorithm as input. The task is to find a vector z ∈ R d such that min{∥z -x∥ 2 , ∥z + x∥ 2 } ≤ ∆ for some precision parameter ∆ > 0.Note that we allow corruption in both the sampling vectors a i ∈ R d and the intensity measurements y i ∈ R. We would like to answer the following algorithmic question:Can we design a provably robust and near sample-optimal algorithm for the ϵ-corrupted phase retrieval problem (Problem 1.3) that runs in nearly-linear time?In this paper, we answer the above question affirmatively. We first state the main result of our paper.Theorem 1.4 (Main, Informal). Consider the outlier-robust phase retrieval problem (Problem 1.3).Let ∆ > 0. Given an ϵ-corrupted set of n = Ω(d log 2 (1/∆)) samples, we can compute z ∈ R d in time O(nd) such that min(∥z -x∥ 2 , ∥z + x∥ 2 ) ≤ ∆ with probability at least 0.8.Our algorithm has near-optimal sample complexity, because even without corruption, recovering the ground-truth vector x in general requires Ω(d) samples because there are d degrees of freedom in x.Moreover, our algorithm runs in time nearly-linear in the size of the input, and provably recovers the ground-truth vector x with arbitrary precision ∆. The formal version of Theorem 1.4 is stated as Theorem 3.1 in Section 3.We remark that the success probability of Theorem 1.4 can be boosted to 1 -δ for any δ > 0 by incurring an additional factor of T = O(log(1/δ)) in the sample complexity and runtime. We can randomly partition the input into T equal-sized disjoint sets and run our algorithm on each set to obtain T solutions Z = {z 1 , . . . , z T }. If we output a solution z ⋆ that has the maximum number of points in Z within distance 2∆, we can show that r(z ⋆ ) ≤ 3∆ with probability at least 1 -δ.Our main conceptual contribution is to propose and study the outlier-robust phase retrieval problem, where a small fraction of the input data is adversarially corrupted. Note that we allow arbitrary corruption in both the sampling vectors a i ∈ R d and the intensity measurements y i ∈ R. The nonconvex optimization landscape of phase retrieval is well understood when the a i 's are Gaussian distributed, but the adversarial robustness of such landscape results is largely unexplored.Our main technical contributions include the design and analysis of a near sample-optimal and nearlylinear time algorithm that provably solves the phase retrieval problem in the presence of outliers.Our approach provides a conceptually simple two-step framework for developing outlier-robust algorithms for tractable nonconvex problems that combines the robustness of spectral initialization and the efficiency of the subsequent robust gradient descent.When there are infinite samples and no corruption, the objective function f (z) can be simplified asEQUATIONEven though f (z) is nonconvex, we know that it has no spurious local optima [34, 3, 40] .Our approach follows the general structure of Candès et al. [3] , which uses a two-step procedure.The first step uses spectral techniques to find an initial guess that is close enough to the ground truth.The second step applies gradient descent to converge to the final solution. However, both steps are susceptible to adversarial corruption. We develop nearly-linear time and provably robust algorithms for both steps and combine them to get our main result.Step 1: Robust Spectral Initialization. When there is no adversarial corruption, the empiricalsecond-moment matrix Y = (1/n) n i=1 y i a i a ⊤ i has expectation E [Y ] = I + 2xx ⊤ , so its topeigenvector is close to x. However, the adversary can arbitrarily change the top eigenvector.To circumvent this issue, we assign a (nonnegative) weight w i to each sample, and let Y w denote the weighted intensity-based second-moment matrixY w = n i=1 w i y i a i a ⊤ i .Ideally, if the weights w are uniformly distributed on the remaining clean samples, the top eigenvector of Y w will align with x.We propose a novel optimization problem that can be used to find a weighting w such that Y w must be close to the unknown unbiased expectation I + 2xx ⊤ . Moreover, we show that such a weight w can be computed in nearly-linear time.Step 2: Approximate Gradient Descent. Starting with the initial guess z 1 ∈ R d produced by the robust spectral initialization, we want to apply gradient descent to recover the ground truth x ∈ R d .Without corruption, if the initialization is close enough to x, each iteration will bring z closer to x by a constant factor. This convergence guarantee can be compromised by the corrupted samples.At a high level, approximating the gradient at a specific point amounts to a robust mean estimation problem (for the underlying distribution of the gradients). When the input data is ϵ-corrupted, the gradients of the n samples can be viewed as an ϵ-corrupted set of vectors. We can approximate the true gradient using this ϵ-corrupted set of n gradients using robust mean estimation algorithms.Phase Retrieval. The problem of phase retrieval arises in many areas of science and engineering [11, 33] . Early research on this problem proposes error-reduction algorithms [22, 17, 18] . Convex and nonconvex optimization with various objective functions were later proposed and achieved exact recovery [43, [3] [4] [5] 38] . Follow-up works generalize to robust phase retrieval where the observations are subject to perturbations [45, 27, 7, 6, 31] .Nonconvex Optimization. Even though optimizing a nonconvex function is NP-Hard in general, recent works showed that many nonconvex functions are locally optimizable due to discrete or rotational symmetry. Besides phase retrieval, it is known that all local optima are globally optimal for natural nonconvex formulations of a wide range of machine learning problems, such as matrix completion [21] , matrix sensing [2] , phase synchronization [1] , dictionary learning [39] , and tensor decomposition [20] (see also Chapter 7 of [44] ). Closely related to our work, a recent line of work explored the robustness of these landscape results: [30] studied matrix sensing in the ϵ-corrupted model and [8, 19] studied matrix completion and matrix sensing in semi-random models.High-Dimensional Robust Statistics. Recent works in high-dimensional robust statistics developed nearly-linear time algorithms for the problem of robust mean estimation [9, 16, 29] . Prior works [35, 14] developed meta-algorithms for finding first-order stationary points with dimension-independent accuracy guarantees, which is closely related to the robust gradient descent procedure that we use.We first introduce notations and background in Section 2. Then we give an overview of our approach in Section 3. Next, we focus on how to get an initialization that is close enough to the ground truth x in Section 4. After the initialization, we use robust mean algorithms to estimate gradients to converge to the desired accuracy in Section 5. Finally, we conclude in Section 6 and discuss open problems. 
preliminary and background: 
 Notation. We write [n] for the set of integers {1, . . . , n}. We use {e 1 , . . . , e d } for the standard unit vector basis in R d and I for the identity matrix. For a vector x, we denote its ℓ 1 , ℓ 2 and ℓ ∞ norm as ∥x∥ 1 , ∥x∥ 2 and ∥x∥ ∞ , respectively, and write the i th coordinate in x as x i . For vectors x, y ∈ R d , we denote its inner product as ⟨x, y⟩ = x ⊤ y. For a matrix A, we use ∥A∥ 2 , ∥A∥ * , and ∥A∥ F as its operator norm, nuclear norm, and Frobenius norm, respectively. We write λ k (A) as the kth-largest eigenvalues of A, and λ k (A) as the sum of the k largest eigenvalues. A symmetric n × n matrix A is said to be positive semidefinite (PSD) if for all vectors x ∈ R n , x ⊤ Ax ≥ 0. For two symmetric matrices A and B, we write A ⪯ B when B -A is positive semidefinite.Packing SDP. We will use nearly-linear time solvers for the following packing SDP. 10] ). Given an instance of optimization (*) with semi-positive definite matrices max w ∥w∥ 1 subject to n i=1 w i A i ⪯ I, λ k n i=1 w i B i ≤ k, w i ≥ 0, ∀i. (*) Lemma 2.1 ([A i ∈ R d1×d1 and B i ∈ R d2×d2 with A i = C i C ⊤ i , B i = D i D ⊤ i for all i = 1, 2outputs w ′ with ∥w ′ ∥ 1 ≥ (1 -ϵ 0 )OPT where OPT is optimal value, with probability at least 1 -δ 0 .Computing the Top Eigenvector. We use power method to compute the top eigenvector of a matrix.Lemma 2.2 (Power Method for Top Eigenvector, e.g., [41] ). Let A ∈ R d×d and let λ 1 be its largest eigenvalue. For any δ ∈ (0, 1), there exists an algorithm that takes A and outputs a unit vectorx ∈ R d in time O(t log(d)/δ) such that x T Ax ≥ (1 -δ)λ1 with probability at least 0.99, where t is the time required to compute Av for an arbitrary v ∈ R d .Robust Mean Estimation. Another tool we use is robust mean estimation in the ϵ-corruption model for distributions with bounded covariance. We use robust mean estimation algorithms to approximate the true gradient under adversarial corruption.Lemma 2.3 (Robust Mean Estimation, e.g., [15] ). Let D be a distribution on R d with unknown mean µ and unknown covariance matrix Σ where Σ ⪯ σ 2 I. Let ϵ 0 be a sufficiently small universal constant. Let 0 < ϵ < ϵ 0 and δ > 0. Given an ϵ-corrupted set of n samples drawn from D, we can output a vector µ ∈ R d in time O(nd log(1/δ)) such that, with probability at least 1 -δ -exp(-nϵ),we have ∥ µ -µ∥ 2 = O √ ϵ + d nδ + d(log d+log 1/δ) n σ.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We point out to all the assumptions used in our paper. There are a discussion of the limitations in the conclusion. Our algorithm can only handle a corruption ϵ < ϵ 0 for some universal constant ϵ 0 . Additionally, it is possible that the dependency on ∆ for the sample complexity can be removed.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. Guidelines: 
overview: 
 We first state a formal version of our main result. We first show that the spectral initialization can be done in nearly linear time with high probability, the proof of which can be found in Section 4.Lemma 3.2 (Robust Spectral Initialization). Under the setting of Problem 1.3, for any 0 < ϵ < ϵ ′ for some universal constant ϵ ′ > 0, given an ϵ-corrupted set of n = Ω(d) samples, we can compute a vec-tor z 0 ∈ R d of the ground truth x in time O(nd) such that r(z 1 ) = min{∥z 1 -x∥ 2 , ∥z 1 + x∥ 2 } ≤ 1 8with probability at least 0.95.Then, with such initialization results, we can proceed to show that an approximate gradient descent algorithm can be used to find an arbitrary approximation of the ground truth in Section 5.Consider the setting of Problem 1.3. Let ∆ > 0 be the desired precision. Let 0 < ϵ < ϵ 0 for some sufficiently small universal constant ϵ 0 . Given an ϵ-corrupted set of n = Ω(d log 2 (1/∆)) samples and an initial guess z 1 such that r(z 1 ) = min(∥z 1 -x∥ 2 , ∥z 1 + x∥ 2 ) ≤ 1/8, we can compute a vector z ∈ R d in time O(nd) such that r(z) ≤ ∆ with probability at least 0.95.For technical reasons, we cannot use the same set of samples for both the robust spectral initialization and the approximate gradient descent. Therefore, we partition the ϵ-corrupted set of 2n samples into two equally sized disjoint sets, using one set for each algorithm. • The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
robust spectral initialization: 
 We dedicate this section to proving Lemma 3.2: Given an ϵ-corrupted set of (a i , y i )'s, we can compute an initial guess z 1 ∈ R d that is close to the ground truth x, where min(∥z1 -x∥ 2 , ∥z 1 + x∥ 2 ) ≤ 1/8.To build some intuition, consider the following intensity-based covariance matrix Y =1 n n i=1 y i a i a ⊤ i ,where each a i is drawn independently from N (0, I) andy i = ⟨a i , x⟩ 2 . The expectation of this matrix is E[Y ] = I + 2xx ⊤ .In other words, when there are enough samples and no adversarial corruption, we can obtain a good guess of the ground truth x (or -x) by computing the top eigenvector of Y . However, we cannot rely on this approach in adversarial settings.To tackle this issue, we propose a nearly-linear time preprocessing step (Algorithm 1) that can recover the true expectation of Y under adversarial corruptions. Algorithm 1 assigns a non-negative weight to each sample. For a weight vector w ∈ R n and a set of indices S ⊆ [n], the weighted intensity-based covariance matrix is defined as Y S,w . = i∈S w i y i a i a ⊤ i , and we omit S when S = [n]. The feasible region for the weight vector is:∆ n,ϵ := w ∈ R n : ∥w∥ 1 = 1 and ∀i ∈ [n], 0 ≤ w i ≤ 1 (1-ϵ)n .A weight w defines an empirical distributions over the samples (a i , y i ) n i=1 , where the largest probability assigned to any point is 1 (1-ϵ)n . Ideally, we would like to find a weight vector w * ∈ ∆ n,ϵ that assigns its weight uniformly to all the uncorrupted samples, i.e.,w * i = 1 (1-ϵ)n • 1 i∈G . To find asuitable weighting w, we use the following optimization problem (**) in which λ 2 returns the sum of the top two eigenvalues (commonly known as the Ky Fank norm for k = 2). min w λ 2 n i=1 w i y i a i a ⊤ i subject to 0 ≤ w i ≤ 1 (1-ϵ)n , ∀i ∈ [n], n i=1 w i = 1 . (**)At a high level, our main observation is that y i a i a ⊤ i is always a positive semidefinite matrix as y i ≥ 0.Consequently, the adversary can only add extra directions with large eigenvalues, but will not be able to remove the eigendirection of x. By minimizing the Ky Fan 2 norm, we can remove any directions added by the adversary and make sure that the only remaining large eigendirection is close to x.Let δ ≥ 0 be some constant to be determined. We show that we can obtain a robust spectral initialization by solving the packing SDP problem (*), which can be solved efficiently using Lemma 2.1.In particular, to fit the reweighting problem of (**) into the framework of the generalized packing problem (*), we define the following constraint matrices for all i ∈ [n] :EQUATIONThe matrices (A i ) n i=1 are used to implement the constraint that w ∈ ∆ n,ϵ . The matrices (B i ) n i=1 help make sure the sum of the top two eigenvalues of Y w must be at most roughly 4, because λ 2 (Y w * ) ≈ 4.Input: ϵ-corrupted samples (a i , y i ) i∈[n] Output: The initial guess z ′ ∈ R d 1: function ROBUSTINIT({(a i , y i )} i∈[n] ) 2:{A i , B i } ← Constraint matrices as defined in Equation (2) 3:w ′′ ← Solution to Optimization (*) with constraints {A i , B i }, k = 2, and precision ϵ 0 = 0.9 as in Lemma 2.1 4:w ′ ← w ′′ / ∥w ′′ ∥ 1 5:z 1 ← TOPEIGENVECTOR(Y w ′ ) as in Lemma 2.2 with sufficiently small constant δ. EQUATIONIn order to show Lemma 4.1, we need the following auxiliary Lemma 4.2, the proof of which can be found in Section A. Intuitively, Lemma 4.2 suggests that any weight w in the feasible region ∆ n,2ϵwill not have a huge impact on the properties of uncorrupted measurements.Lemma 4.2. For any δ 0 > 0, and sufficiently small ϵ ≥ 0, given a set of n ϵ-corrupted samples with n > Ω(d), with probability at least 0.98, we have Y G,w -(I + 2xx ⊤ ) 2 ≤ δ 0 for all w ∈ ∆ n,2ϵ .Using Lemma 4.2, we provide a proof sketch for Lemma 4.1, and defer the details to Section A.Proof. We condition on the fact that the event of Lemma 4.2 holds (with probability at least 0.98)for δ 0 = δ. Thus, for the remaining of the proof, we assume that for all w ∈ ∆ n,2ϵ , it holds thatY G,w -(I + 2xx ⊤ ) 2 ≤ δ.Let λ 1 and λ 2 be the top two eigenvalues of Y w ′ , with v 1 and v 2 to be their corresponding eigenvectors.Note that the largest eigenvalue of I + 2xx ⊤ is 3, and the rest of the eigenvalues are all 1. In the proof, we show that the eigenvalues of Y w ′ are also close to the ones of I + 2xx ⊤ . Our proof consists of two parts. We first establish lower bounds for λ 1 and λ 2 , and then find an upper bound for λ 1 + λ 2 .Lower Bound. Since y i a i a ⊤ i ⪰ 0 for any i ∈ [n], for any positive weight vector w ∈ ∆ n,ϵ , we have Y G,w ⪯ Y w . Thus a lower bound on eigenvalues of Y G,w ′ will also be a lower bound on Y w ′ .For the top eigenvalue λ 1 of Y w ′ , it holdsEQUATION) Similarly, for the second largest eigenvalue λ 2 of Y w ′ , we have:EQUATION)Upper Bound. Through the optimization problem (*), a weight w ′′ is calculated such that Y w ′′ are operator-norm upper-bounded by the constraint parameters. Let OPT be the value of the optimal solution of the optimization problem (*). The desired uniform weight vector over the good samples w * ∈ ∆ n,2ϵ is also a feasible solution to this optimization problem because Y w * satisfy the optimization constraints due to Lemma 4.2. Since w ′′ is an ϵ 0 -approximation to the problem, we have  and consequently ,∥w ′′ ∥ 1 ≥ (1 -ϵ 0 )OPT ≥ (1 -ϵ 0 ) ∥w * ∥ 1 = 1 -ϵ 0 By optimization constraints, the Ky Fan 2-norm of i w ′′ i B i = 1 2 (1 -δ)Y w ′′ ≤ 2,EQUATIONBy combining inequalities (4), (5) , and (6), we have shown that the top two eigenvalues of Y w ′ are close to 3 and 1. Since the rest of the eigenvalues of Y w ′ can also be bounded, we can conclude thatY w ′ -(I + 2xx ⊤ ) 2 = O(δ).We can now show the closeness between the top eigenvector of Y w ′ and the ground truth.There exists an universal constant ϵ ′ such that if 0 ≤ ϵ ≤ ϵ ′ , and Algorithm 1 receives in input an ϵ-corrupted set of samples, then it outputs z 1 ∈ R d such that with probability at least 0.95 it holds r(z 1 ) ≤ 1 8 .Proof. We condition on the fact that the event of Lemma 4.1 holds (with probability at least 0.98).Let the eigendecomposition ofY w ′ be Y w ′ = i∈[d] λ i v i v ⊤ i , where λ 1 ≥ . . . ≥ λ d . Under the basis {v 1 , • • • , v d }, the ground truth x can be represented as x = i∈[d] α i v i . Note that ∥x∥ 2 2 = i∈[d] α 2 i = 1. By Lemma 4.1, we have Y w ′ -(I + 2xx ⊤ ) 2 = O(δ). Thus, we havex ⊤ Y w ′ x ≥ 3 -O(δ) and x ⊤ Y w ′ x ≤ λ 1 α 2 1 + λ 2 (1 -α 2 1 ) ≤ (3 + O(δ))α 2 1 + (1 + O(δ))(1 -α 2 1 ) ≤ 1 + 2α 2 + O(δ). This implies α 2 1 ≥ 1 -O(δ). As a result, r 2 (v 1 ) = min{∥v 1 -x∥ 2 2 , ∥v 1 + x∥ 2 2 } = min{(α 1 -1) 2 , (α 1 + 1) 2 } + d i=2 α 2 i = min{2 -2α 1 , 2 + 2α 1 } = O(δ).The last inequality holds as long as δ is sufficiently small. Let z 1 = i∈[d] β i v i be the unit vector approximating v 1 returned by the algorithm. By Lemma 2.2, we have thatz ⊤ 1 Y w ′ z 1 ≥ (1 -δ)λ 1with probability at least 0.99. Thus, we have:z ⊤ 1 Y w ′ z 1 ≥ (1 -δ)λ 1 ≥ (1 -δ)(3 -O(δ)) ≥ 3 -O(δ + δ) and z ⊤ 1 Y w ′ z 1 ≤ λ 1 α 2 1 + λ 2 (1 -α 2 1 ) ≤ 1 + 2β 2 1 + O(δ). Again, this implies that β 2 1 ≥ 1 -O(δ + δ). We can show that min{∥v 1 -z 1 ∥ 2 2 , ∥v 1 + z 1 ∥ 2 2 } = O(δ + δ).By the triangle inequality, we can conclude that r 2 (z 1 ) = O(δ + δ) ≤ 1/64, where the last inequality is obtained by choosing sufficiently small δ and δ. Therefore, there exists an universal constant ϵ ′ ≥ 0 such that for all 0 ≤ ϵ ≤ ϵ ′ , Algorithm 1 takes n = Ω(d) samples and outputs z 1 such that r(z 1 ) ≤ 1/8 with probability at least 0.95. We can directly combine Lemma 4.3 and Lemma 4.4 to finish the proof of Lemma 3.2. 
robust gradient descent: 
 After the robust spectral initialization in Section 4, we have an initial guess z 1 ∈ R d that is close to the ground truth x or -x. Without loss of generality, we can assume that z 1 is closer to x than to -x.In this section, we prove Lemma 3.3: Given an initial guess z 1 with ∥z 1 -x∥ 2 ≤ 1/8, we can use a robust gradient descent algorithm (Algorithm 2) to recover x to any desire precision ∆ > 0. It is well-known that gradient descent can achieve geometric convergence rates in non-adversarial settings.We show that Algorithm 2 achieves a similar convergence rate even when the input is ϵ-corrupted.Consider the natural nonconvex formulation:min z∈R d n i=1 f i (z) where f i (z) . = ⟨a i , z⟩ 2 -y i 2 .Let g i denote the gradient of f i with respect to z. Let D z denote the distribution of g i (z) ∈ R d when there is no adversarial corruption. Formally, g(z) ∼ D z is distributed asEQUATION)To run gradient descent, we want to estimate the expected true gradient µ z . = E g(z) using samples.The challenge is that the input samples {(a i , y i )} i∈[n] are ϵ-corrupted, and consequently the gradients{g i (z)} i∈[n]is an ϵ-corrupted set of vectors drawn from D z . To address this, we use robust mean estimation algorithms (e.g., [16] ) to approximate µ z , the true mean of D z .Algorithm 2 Robust Gradient Descent Input: ϵ > 0, an ϵ-corrupted set of n samples {(a i , y i )} i∈[n] , initial guess z 1 ∈ R d with ∥z 1 -x∥ ≤ 1/8, and desired precision ∆ > 0.Output: z ∈ R d such that ∥z -x∥ 2 ≤ ∆ where x is the ground truth.1: procedure ROBUSTGD(ϵ, {(a i , y i )} i∈[n] , z 1 , ∆) 2: T ← O(log(1/∆)), η ← 1/300 3: {N 1 , N 2 , • • • , N T } ← a random disjoint partition of [n] such that |N t | = n/T for all t ∈ [T ] 4:for t = 1, 2, . . . , T do 5:µ zt ← Robust mean estimation on input {g i (z t )} i∈Nt using Lemma 5. Proof of Lemma 5.2. Since 2ϵ < ϵ 0 , we can view the (2ϵ)-corrupted set of m samples as ϵ 0 -corrupted.We need to replace 2ϵ with ϵ 0 to reduce the failure probability of Lemma 2.3. This weakens the error guarantee of Lemma 2.3, but the resulting µ z is still accurate enough for our algorithm.We apply Lemma 2.3 to the ϵ 0 -corrupted set of m vectors drawn from D z . By Lemma 5.1, the covariance matrix of D z satisfies Σ z ⪯ O(∥z -x∥ O √ ϵ 0 + d mδ + d(log d+log(1/δ)) m ∥z -x∥ 2 ≤ c ∥z -x∥ 2 .The success probability is at least1 -δ -exp(-ϵ 0 m) = 1 -O(δ).Lemma 5.2 shows that even with a (2ϵ)-corrupted set of gradients, the true gradient µ z can be estimated up to an additive error proportional to the distance between z and x. The next lemma shows that such an approximate gradient is sufficient for gradient descent to converge, reducing the distance to the ground truth x by a constant factor in each iteration. . When ∥ µ z -µ z ∥ 2 ≤ c ∥z -x∥ 2 , ⟨ µ z , z -x⟩ ≥ (7.5 -c) ∥z t -x∥ 2 2 and ∥ µ zt ∥ 2 ≤ (29 + c) ∥z t -x∥ 2 .When c < 7.5, we can choose an appropriate step size η such that the distance between z t and x decreases by a constant factor in each iteration.We are now ready to prove Lemma 3.3, which states the performance guarantee, sample complexity, runtime, and success probability of Algorithm 2. We restate Lemma 3.3 before proving it. Lemma 3.3 (Robust Gradient Descent). Consider the setting of Problem 1.3. Let ∆ > 0 be the desired precision. Let 0 < ϵ < ϵ 0 for some sufficiently small universal constant ϵ 0 . Given an ϵ-corrupted set of n = Ω(d log 2 (1/∆)) samples and an initial guess z 1 such that r(z 1 ) = min(∥z 1 -x∥ 2 , ∥z 1 + x∥ 2 ) ≤ 1/8, we can compute a vector z ∈ R d in time O(nd) such that r(z) ≤ ∆ with probability at least 0.95.Proof of Lemma 3.3. First, we prove the success probability of Algorithm 2. Algorithm 2 can fail in two ways: (i) if some N t has more than (2ϵ)-fraction of corrupted samples, or (ii) if Lemma 5.2 fails in some iteration t. The probability of event (i) is at most 0.01 for our choice of n, which follows from a standard application of Hoeffding's inequality and a union bound over T iterations. For event (ii), we choose a sufficiently small δ = O(1/T ) in Lemma 5.2, so each robust gradient estimation fails with probability at most O(δ) = 0.01/T , and overall the probability of event (ii) is at most 0.01.For the rest of the proof, we assume these bad events do not happen.Next, we show the correctness of Algorithm 2. Without loss of generality, we can assume that z 1 is closer to the ground truth x than to -x, which implies ∥z 1 -x∥ 2 ≤ 1/8. By Lemma 5.2, we can obtain an approximation µ z1 of the true gradientµ z1 at z 1 such that ∥ µ z1 -µ z1 ∥ 2 ≤ c ∥z 1 -x∥ 2 .Then by Lemma 5.3, we know that ∥z 2 -x∥ 2 ≤ 0.99 ∥z 1 -x∥ 2 after one step of gradient descent.We can iteratively apply these two lemmas to show that, after T = O(log(1/∆)) iterations, we have ∥z T +1 -x∥ 2 ≤ ∆. One technical issue is that in iteration t, we need to use a fresh subset of samples N t . By the principle of deferred decisions, we can view (a i , y i ) i∈Nt as being generated (and corrupted) after z t is chosen, which forms a (2ϵ)-corrupted set of gradients at z t .Finally, we analyze the sample complexity and runtime of Algorithm 2. Algorithm 2 requires in total n = mT = Ω(d log d log 2 (1/∆)) samples. A random partition can be computed in O(n) time by shuffling the input. In each iteration, the m gradients in N t can be computed using Equation (7)  
conclusions and future directions: 
 In this paper, our main conceptual contribution is to propose and study the outlier-robust phase retrieval problem, where a constant fraction of the input data is corrupted. Notably, we allow corruption in both the sampled frequencies a i ∈ R d and the intensity measurements y i ∈ R. Our main technical contribution is the design and analysis of a near-sample-optimal and nearly-linear-time algorithm that solves this problem with provably guarantees.An immediate technical question is whether our sample complexity can be tightened by removing some log(1/∆) factors. One potential approach is to open robust mean estimation algorithms instead of using them in a black-box manner. One could examine the stability conditions that these algorithms require, and see if these stability conditions can be proved without partitioning the samples and using fresh samples in each iteration. More broadly, we believe our framework can be used to develop outlier-robust algorithms for other tractable nonconvex problems, by first finding an initial solution in a saddle-free region near a global optimum and then converging to this global optimum using robust gradient descent.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [NA]Justification: The paper is a theory paper, and it has no experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [NA]Justification: The paper is a theory paper, and it has no experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [NA]Justification: The paper is a theory paper, and it has no experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA]Justification: The paper provides a theoretical result. We do not believe there are potential societal consequences of our work, aside from advancing the field of Machine Learning.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The paper only contains theoretical results.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA]Justification: The paper only contains theoretical results. Previous work is properly cited.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: Lemma A.1. [Lemma 4.2, Formal]. For any δ 0 > 0, there exists constants ϵ 0 , c > 0 such that when n > c • d log d and we are given a set of n ϵ-corrupted samples, where 0 ≤ ϵ ≤ ϵ 0 , then with probability at least 0.98, it holdsEQUATION)Proof of Lemma A.1. We recall the definition of Y G,w = i∈G w i y i a i a ⊤ i . Let ℓ = ϵ • n and let {(a n+i , y n+i )} ℓ i=1 be the set of samples that were removed by the ϵ-corruption adversary. Let G ′ = G ∪ {n + 1, . . . , n + ℓ}, n ′ = n + ℓ, and ϵ ′ = ϵ/(1 + ϵ). Note that without loss of generality, we can assume that |G| = (1 -ϵ)n and|G ′ | = (1 -ϵ ′ )n ′ = n.We define a mapping σ :EQUATION)In other words, all the weights are the same for the samples with index in the set [n], and are equal to 0 for the samples removed by the adversary. We can verify that σ(w)∈ ∆ n ′ ,3ϵ ′ for all w ∈ ∆ n,2ϵ since σ(w) i ≤ w i ≤ 1/(1 -2ϵ)n = 1/(1 -3ϵ ′ )n ′ for all i ∈ [n ′ ],and ∥σ(w)∥ 1 = ∥w∥ 1 = 1.Furthermore, we have Y G,w = Y G ′ ,σ(w) for all w ∈ ∆ n,2ϵ . We denote with w * ∈ ∆ n ′ ,3ϵ ′ the desired uniform weighting of the samples with index in G ′ , i.e.,w * i = 1 (1-ϵ ′ )n ′ 1 i∈G ′ . It suffices to show both Y G ′ ,w * -(I + 2xx ⊤ ) 2 ≤ δ 0 /2 and Y G ′ ,σ(w)-w * 2 ≤ δ 0 /2.By triangle inequality, for any w ∈ ∆ n,2ϵ , it holdsY G ′ ,σ(w) -(I + 2xx ⊤ ) 2 ≤ Y G ′ ,w * -(I + 2xx ⊤ ) 2 + Y G ′ ,σ(w)-w * 2 , Thus, it suffices to show both Y G ′ ,w * -(I + 2xx ⊤ ) 2 ≤ δ 0 /2 and Y G ′ ,σ(w)-w * 2 ≤ δ 0 /2.We upper bound the first term. By using the definition of w * , note that Y G ′ ,σ(w) -(I + 2xx ⊤ ) 2 = i∈G ′ w * i y i a i a ⊤ i -(I + 2xx ⊤ ) 2 = i∈G ′ 1 |G ′ | y i a i a ⊤ i -(I + 2xx ⊤ ) 2 . Since E(y i a i a ⊤ i ) = I + 2xx ⊤ for any i ∈ G ′ ,EQUATION)It remains to show a high-probability upper bound to the second termY G ′ ,w * -σ(w) 2 ≤ δ 0 /2.The first observation is that the weighting w * and σ(w) for any w ∈ ∆ n,2ϵ cannot too different. In particular, we can show the following upper bound:EQUATIONWe observe that ∆ n ′ ,3ϵ ′ can be seen as the convex combination of all possible uniform weighting over subsets of n ′ (1 -3ϵ ′ ) samples. Thus, the maximum distance will be between two points of the convex hull, and we can upper bound (11) as:EQUATION)For a fixed unit vector z ∈ S d-1 with z = px + qu where u ∈ S d-1 and ⟨u, x⟩ = 0, we have:max w∈∆n,2ϵ z ⊤ Y G ′ ,w * -σ(w) z = max w i∈G ′ (w * i -σ(w) i ) ⟨a i , x⟩ 2 ⟨a i , z⟩ 2 = max w i∈G ′ (w * i -σ(w) i ) ⟨a i , x⟩ 2 (p ⟨a i , x⟩ + q ⟨a i , u⟩) 2 ≤ 2 max w i∈G ′ (w * i -σ(w) i )(⟨a i , x⟩ 4 + ⟨a i , x⟩ 2 ⟨a i , u⟩ 2 ) ≤ 2 max w i∈G ′ |w * i -σ(w) i | ⟨a i , x⟩ 4 + 2 max w i∈G ′ |w * i -σ(w) i | ⟨a i , x⟩ 2 ⟨a i , u⟩ 2 .For ease of notation, letβ i . = |w * i -σ(w) i |. Observe that 0 ≤ β i ≤ 1 (1-2ϵ)n for all i, and i β i ≤ 6ϵ due to (11) . We have:EQUATIONInequality (13) follows by assigning the maximum possible β i to the largest entries of the sum until we hit the budget 6ϵ due to (11) .We bound max L i∈L ⟨a i , x⟩ 4 first. Let X i = ⟨a i , x⟩ ∼ N (0, 1) for i ∈ G ′ , and define the threshold function h r (z) = 0, z ≤ r z, z > r for r = C 2 • ln 2 (1/ϵ) with constant C > 0 to be determined.Note that z ≤ r + h r (z) for all z > 0. Therefore,max L⊆G ′ ,|L|=6ϵn 1 n i∈L X 4 i ≤ max L 1 n i∈L r + max L 1 n i∈L h r (X 4 i ) ≤ 6ϵr + 1 n i∈G ′ h r (X 4 i ).Then, we consider to bound exp i∈G ′ c • h r (X 4 i ) for some c > 0 to be determined. For any i ∈ G ′ and z ≥ 1, with C = 6, for all ϵ > 0, we havePr exp c • h r (X 4 i ) ≥ z ≤ Pr h r (X 4 i ) ≥ 0 = Pr X i ≥ r 1/4 ≤ exp(- √ r/2) ≤ exp(ln ϵ • C/2) ≤ ϵ 3 .At the same time, with c < 1/200, for all z ≥ 1, we havePr exp c • h r (X 4 i ) ≥ z ≤ Pr exp(c • X 4 i ) ≥ z ≤ Pr X 4 i ≥ ln z c ≤ exp - ln z c /2 ≤ z -3 .Therefore, Pr exp c • h r (X 4 i ) ≥ z ≤ min{z -3 , ϵ 3 }, and for all ϵ < 1/2, we haveE exp c • h r (X 4 i ) = ∞ 0 Pr exp c • h r (X 4 i ) ≥ z dz ≤ 1 0 1dz + 1/ϵ 1 ϵ 3 dz + ∞ 1/ϵ z -3 dz = 1 + ϵ 2 -ϵ 3 + 1 2 ϵ 2 ≤ 1 + ϵ 2 ≤ exp(ϵ 2 ).Since {X i } i∈G ′ are independent, we haveE exp i∈G ′ c • h r (X 4 i ) = E i∈G ′ exp c • h r (X 4 i ) ≤ exp(ϵ 2 n).Finally, by Markov's inequality, for any constant δ 1 > 0, as long as ϵ ≤ √ δ 1 c, we havePr i∈G ′ h r (X 4 i ) ≥ 2δ 1 n = Pr exp i∈G ′ c • h r (X 4 i ) ≥ exp(2δ 1 cn) ≤ exp(ϵ 2 n -2δ 1 cn) ≤ exp(-δ 1 cn).As a result, with sufficiently large n ≥ c 2 (δ 1 ) • d log d and sufficiently small ϵ such that 6ϵr ≤ δ 1 ,Pr ϵr + 1 n i∈G ′ h r (X 4 i ) ≥ 3δ 1 ≤ exp(-δ 1 cn) ≤ 0.995 • 9 -d .Also, max L i∈L ⟨a i , x⟩ 2 ⟨a i , u⟩ 2 has a similar tail bound and therefore can be bounded in the same way. We can then bound the operator norm via an epsilon-net argument. Set δ 1 = δ 0 /24. By an1/4-net on S d-1 with |N | ≤ 9 d , we have that Pr max z∈S d-1 z ⊤ Y G ′ ,w * -σ(w) z ≥ δ 0 /2 ≤ 9 d • 0.99 • 9 -d ≤ 0.99.By combining the above inequality with Equation (10), we can conclude that, for any δ 0 > 0, there exists ϵ 0 > 0, such that when n ≥ max{c 1 (δ 0 ), c 2 (δ 0 )} • d log d and 0 ≤ ϵ ≤ ϵ 0 , with probability at least 0.98,∀w ∈ ∆ n,2ϵ , Y G,w -(I + 2xx ⊤ ) 2 = Y G ′ ,σ(w) -(I + 2xx ⊤ ) 2 ≤ δ 0 .Therefore, with probability at least 0.98, we have Y G,w -(I + 2xx ⊤ ) 2 ≤ δ 0 for all w ∈ ∆ n,2ϵ .For the undisturbed samples, we have the following concentration result.Lemma A.2 ([3] Section A.4.2). Let x ∈ R d .For any δ > 0, there exists a constant C(δ) > 0 such that when n > C • d log d and we are given a set of n samples {(a i , y i )} n i=1 with a i ∼ N (0, I) independently and y i = ⟨a i , x⟩ 2 for all i ∈ [n], then with probability at least 0.99, it holds1 n n i=1 y i a i a ⊤ i -(I + 2xx ⊤ ) 2 ≤ δ.B Omitted Proofs in Section 5Lemma 5.1. Let D z be the distribution of gradients at z as defined in Equation (7) .EQUATION2 )I.Proof of Lemma 5.1. Recall that g ∼ D z is distributed asg = ∂ ∂z ⟨a, z⟩ 2 -⟨a, x⟩ 2 2= -4 ⟨a, z⟩ 2 -⟨a, x⟩ 2 ⟨a, z⟩a where a ∼ N (0, 1) .Because Eg∼D z [g] = µ z , the spectral norm of Σ z can be upper bounded as follows:∥Σ z ∥ 2 = E g∼Dz gg ⊤ -µ z µ ⊤ z 2 ≤ E g∼Dz gg ⊤ 2 .Consequently, it suffices to upper bound Eg∼D z gg ⊤ 2 . Let h = z -x.EQUATION= 16 max We can compute the mean of D z using moments of Gaussian:EQUATIONConsider one step of gradient descent in Algorithm 2: z t+1 = z t -η µ zt , where µ zt is an approximate gradient. We have∥z t+1 -x∥ 2 2 = ∥z t -η µ zt -x∥ 2 2 = ∥z t -x∥ 2 2 -2η⟨ µ zt , z t -x⟩ + η 2 ⟨ µ zt , µ zt ⟩To prove convergence, we need to lower bound ⟨ µ zt , z t -x⟩ and upper bound ⟨ µ zt , µ zt ⟩.We write z = z t and h = z -x to simplify notation. We can substitute z = x + h in Equation ( 14): Justification: We provide a formal discussion of our theoretical claim in the first section of the paper, and include a formal statement of our resul in Section 3.µ z = 16⟨x, h⟩ + 12 ∥h∥Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.• The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.• It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.We write G ⊆ [n] for the (remaining) good samples, and B = [n] \ G for the corrupted samples. 


Decision: The paper studies a variant of phase retrieval, in which we observe y_i = |<a_i, x >| and the goal is to recover x. Here, a fraction of the pairs (a_i, y_i) can be adversarially corrupted. The paper proposes a two-phase algorithm which deploys tools of robust statistics in each phase. In the first, initialization, phase, the proposal is to modify the standard spectral initialization by weights which control the influence of adversarially modified pairs (a_i,y_i) on the lead eigenvector. These weights are chosen by solving a certain SDP, which admits efficient special-purpose solvers. In the second, estimation phase, the paper deploys a robust mean estimator to reduce the effect of outlying samples on the gradient. 

Reviewers produced a mixed evaluation of the paper. On the positive side, the paper analyzes robust phase retrieval under a novel setting, in which both a_i and y_i can be modified; while a number of past works utilize robust estimators for nonconvex robust phase retrieval, these works focus on corrupted measurements y_i, rather than pairs (a_i,y_i). The paper is crisply written. Reviewers raised questions about the role of the corruption level \epsilon, the absence of experimental validation and the relationship to the existing literature on phase retrieval. The authors convincingly addressed questions about the level of corruption; it is indeed a standard assumption in robust statistics that the corruption level is bounded by an absolution constant. Other reviewer concerns centered around the significance of this formulation of robust phase retrieval vis-a-vis existing models, which could be better addressed in future versions of this submission. After considering author feedback, reviewers retained a mixed evaluation of the work, putting it below the bar for acceptance.

Reviews: review: paper_summary: paper_summary: The authors study the phase retrieval problem for retrieving a real signal under the influence of arbitrary corruption. The corruption is allowed to be present in labels or features. They propose a two-step solution. First, they ensure that the initialization is robust to the corruption and second, they show that the gradient descent updates can be made resilient. Authors claim that their method can recover the true signal (with possibly sign mismatch) to an arbitrary precision.
strength_weakness: strength_weakness: 1. The authors claim that corruption level up to some universal constant $\epsilon'$ can be handled through their method. Although, to the best of my understanding, this quantity is not characterized in the main paper. What is the maximum value for $\epsilon'$? 
2. There is no discussion on the dependency of $\epsilon'$ on $n$ or $d$.   
3. The claim of signal recovery to an arbitrary precision puzzles me. It is known that for $\epsilon$-corrupted vectors, the robust mean estimation can only be done up to $\Omega(\sqrt{\epsilon})$ error. Despite that, the authors claim signal recovery (with possibly a flipped sign) to an arbitrary precision. Can the authors comment on how this is achieved?   
4. The claim in line 213 says that $y_i$ is always greater than $0$. Why is that true when the adversary can corrupt $y_i$ arbitrarily? As far as I can tell, the algorithm does not discard negative $y_i$s.
questions: questions: Please see above.
limitations: limitations: Please see above.
paper_summary: paper_summary: This paper studies a classical problem called phase retrieval. The goal is to obtain unknown $d$-dimensional vector $x$ from $n$ datapoints $(a_i, \langle a_i, x \rangle^2)$. This work assumes that $a_i$ are iid Gaussian vectors, but also that a small $\varepsilon$ fraction of the data is corrupted. The authors suggest a two-stage process to identify the vector up to a small error. First, a spectral based algorithm is used to have a small constant error. Further, a robust gradient descent is used to approximate the initial guess up to a small error.
strength_weakness: strength_weakness: 1. Results are limited to the Gaussian setting.
2. The method for RME that is used assumes that variance $\sigma$ is known? But in the way it is used here, $\sigma$ depends on the distance between current solution and the true vector. Authors do not comment on this issue.
3. $\tilde O, \tilde \Omega$ notation is not defined.
4.  Intuition in line 98 in my interpretation contradicts more exact version in line 214 (In the end, if I understand correctly, the crucial reason why the spectral initialization algorithm works is that the adversary cannot change the top eigendirection, but can only add new directions).
questions: questions: 1. How applicable are the techniques beyond Gaussian iid setting?
2. Why using $k = 1$ is not sufficient for the optimization problem in Algorithm 1?
limitations: limitations: There are no ethical limitations of this work.
paper_summary: paper_summary: This paper focuses on the problem of outlier robust phase retrieval, whose goal is to recover a vector $x \in \mathbb{R}^d$ from $n$ intensity measurements $y_i = (a_i^\top x)^2$ when a small fraction of the samples are adversarially corrupted. The authors propose and study this problem, providing a nearly sample-optimal and nearly linear-time algorithm to recover the ground-truth vector $x$ in the presence of outliers.
strength_weakness: strength_weakness: 1. The analysis appears to be more incremental in nature compared to those in previous theoretical works related to phase retrieval. More precisely, the key novelty of the spectral initialization step resides in the assignment of a nonnegative weight to each sample. Regarding the gradient descent step, the problem seems to be simplified to the analysis of robust mean estimation algorithms.

2. Important references are missing. For instance, the authors ought to cite the works related to robust compressed sensing (and it would be better to discuss in more detail the disparities between the analysis for the robust gradient descent step in this work and the analysis in these relevant works), such as

- Liu, Liu, Yanyao Shen, Tianyang Li, and Constantine Caramanis. "High dimensional robust sparse regression." In International Conference on Artificial Intelligence and Statistics, pp. 411-421. PMLR, 2020.
- Liu, Liu, Tianyang Li, and Constantine Caramanis. "High Dimensional Robust $ M $-Estimation: Arbitrary Corruption and Heavy Tails." arXiv preprint arXiv:1901.08237 (2019).

3. The authors solely consider the scenario of noiseless intensity measurements and fail to take into account the noisy case.

4. The paper does not include experimental results, which could limit the confidence in the practical effectiveness of the proposed approach.
questions: questions: Why in Theorem 1.4 (or its formal version in Theorem 3.1), there is no dependence on $\epsilon$ in both the sample complexity and the upper bound of $\min (\\|z-x\\|_2,\\|z+x\\|_2)$?
limitations: limitations: No experimental validations.
paper_summary: paper_summary: This paper addresses the challenge of achieving outlier robustness in phase retrieval, specifically focusing on the recovery of real-valued signals from intensity measurements that have been corrupted by adversarial outliers. The contribution of this work is the development of a nearly-linear time algorithm that is nearly sample-optimal and can accurately recover the true vector despite the presence of outliers. This is done through a two-step process that involves robust spectral initialization and robust gradient descent, utilizing recent results in high-dimensional robust statistics.
strength_weakness: strength_weakness: - The paper's theoretical results assume the corruption level $\epsilon$ as a constant in Theorem 3.1, despite it being initially introduced as a variable in Definition 1.2 to denote the extent of sample corruption. This constant treatment affects the robustness of the results by neglecting the influence of $\epsilon$ on the sample complexity. A thorough analysis that explicitly considers the variability of $\epsilon$ would significantly strengthen the theoretical foundation.
- Lack of numerical validations for the theoretical claims
- The idea and design of the two-stage robust phase retrieval algorithm is not novel, by adapting existing reweighting phase retrieval algorithms (throught different design of the reweights) to handle outliers through robust statistics. The contamination model is also from existing works, which has been recent studied in a number of different contexts; e.g., robust linear/nonlinear estimation by e.g., Diakonikolas and coauthors.
- Lack of numerical validation and comparison of the proposed algorithm with respect to existing (robust) Gaussian phase retrieval algorithms. It is difficult to judge if the proposed algorithm is of practical interest (or only of theoretical interest).
- Quite a lot statements in the paper are rather confusing, e.g., "we propose the problem of outlier robust phase retrieval"; robust phase retrieval has been long studied in the literature; as far as I understand, the paper studies a new robust phase retrieval problem by considering also adversarial a_i's on top of existing formualtions. ii) "It is well-known that natural nonconvex formulations of phase retrieval do not have spurious local optima." which is not precise enough and is true under very stringent assumptions. iii) "This is first achieved via approaches based on semidefinite programming (SDP) relaxations (see, e.g., Cand`es et al. (2015c))." I guess the first Gaussian phase retrieval algorithm was the AltMin algorithm (Netrapalli, et al 2013. Phase retrieval using alternating minimization. Advances in Neural Information Processing Systems, 26.) which was interestingly not cited in the submission. iv) "Similar landscape results are known for other natural nonconvex formulations of phase retrieval as well (e.g., min f(z) = P i( √ yi − |⟨ai, z⟩|)2 (Soltanolkotabi, 2019))." The first work to study the Gaussian phase retrieval based on the magnitude-based least-squares nonconvex formulation and achieve provable guarantees is (Wang et al (2017). Solving systems of random quadratic equations via truncated amplitude flow. IEEE Transactions on Information Theory. 64(2):773-94.) Please be careful and fair in stating related results.
questions: questions: i) The main motivation of the paper is to study a new robust phase retrieval (RPR) problem where in addition to the existing well-studied RPR (earlier in e.g., Hand, et al. 2016. Corruption robust phase retrieval via linear programming. arXiv preprint arXiv:1612.03547.), the sampling vectors a_i's can also be adversarially manipuated. Nonetheless, since PR is quite "simple" model, it can be imagined some of the effect of the adversarial a_i's can also be lumped into the adversarial y_i's. I am not sure whether the problem here is indeed more challenging than existing corruption pr considering only adversarial y_i's. If not rigorously theoretically analyzed, numerical comparisons with these existing rpr algorithms will be needed and appreciated. 
ii) The designed algorithm is a two-stage reweighting algorithm that consists of a reweighting spectral initialization and reweighting gradient descent algorithm; nonetheless, this unique design and use of reweighting to improve nonconvex PR performance has been initially documented in [Yuan, et al (2017). Phase retrieval via reweighted Wirtinger flow. Applied Optics, 56(9), 2418-2427; Wang, et al (2018). Phase retrieval via reweighted amplitude flow. IEEE Transactions on Signal Processing. 66(11):2818-33.]; the form of using the weights here is exactly the same as in the two papers]; The only difference is the specific option/design of the values of the weights, certainly for considerations of handling the outliers here. It is very interesting to understand what is the original motivation of using and where does the reweighting technique come from. 
iii) Last but not the least, from a practical perspective, say e.g., the Fourier/modulated phase retrieval problem using random masks, it makes sense that possiblem (adversarial) noise may be present in measuring the intensity measurements; but it is not clear why the sampling vectors that correspond to rows of DFT matrices (or multiplied by random masking +1 or -1 values) could also be manupulated experimentally. Please provide a more practical setup validating the motivation of studying such a RPR problem. 
iv) A lot of typos or inaccurate statements. e.g., "we propose and study..."; the claim "into a saddle-free region" is not theoretically and rigorously established.
limitations: limitations: NA

Comments: title: Response to rebuttals
comment: The TV-distance between $P_1$ and $P_2$ is always $1$ as long as $\theta_1 \neq \theta_2$, as a result the TV-distance between $D_1$ and $D_2$ cannot be less than $1-2\epsilon$ (you gave an upper bound of $1-\epsilon$ but what is more important is the lower bound), so in the 1-D case, we can achieve exact recovery with for example $n \ge 3$ samples when $\epsilon \le 1/3$.

Note that our task is not to distinguish whether a single sample come from $D_1$ or $D_2$. Our task (the property testing version) is, given an $\epsilon$-corrupted set of $O(d)$ samples drawn from either $D_1$ or $D_2$, decide whether they are drawn from $D_1$ or $D_2$ with high constant probability (which is what we claim in Theorem 3.1).
Thank you for your detailed follow-up. We will respond using the setting and notations you provided.

In the 1-D case, let $P_1$ and $P_2$ denote the clean distributions corresponding to the ground-truth parameters $\theta_1$ and $\theta_2$ respectively. Specifically, $(x, y) \sim P_1$ is distributed as first drawing $x \sim \mathcal{N}(0, 1)$, and then setting $y = (\theta_1 x)^2$. The distribution $P_2$ is defined similarly where $y = (\theta_2 x)^2$.

As long as $\theta_1 \neq \pm \theta_2$, the total variation distance between $P_1$ and $P_2$ is always $1$. These two distributions have essentially disjoint support. This means that an adversary cannot make $P_1$ and $P_2$ indistinguishable by corrupting an $\epsilon$-fraction of the samples.

In particular, if the corruption level is less than half, one can compute the (multi-)set $S = \{\sqrt{y_i/x_i^2}\}$, and the most frequent value in $S$ will be the correct $\theta$ with probability $1$.
Thank you for your response. However, I still have reservations about exact recovery in the constant corruption proportion case with a strong adversary model. Intuitively, a strong adversary could make an $\epsilon$ proportion of inputs identical for inputs coming from two different distributions. These distributions would be, at most, $1 - \epsilon$ away from each other in terms of total variation distance. This is problematic when $\epsilon$ is a constant (non-vanishing). Let me elaborate.

Consider $\epsilon \in (0, 1]$ to be the constant corruption proportion. We will use the following phase retrieval model in one dimension:
$$ y = (x \theta)^2 + w_{\theta, x}$$
where $w_{\theta, x}$ denotes the adversarial corruption added by a strong adversary who has access to both $x$ and $\theta$. We draw $x$ from a standard normal distribution. Consider two parameters $\theta_1 > 0$ and $\theta_2 > 0$ with $| \theta_1 - \theta_2 | > \delta$ for some $\delta > 0$. Let $D_1$ and $D_2$ be distributions over $\mathbb{R} \times \mathbb{R}$ corresponding to linear models $y = (x\theta_1)^2 + w_{\theta_1, x}$ and $y = (x\theta_2)^2 + w_{\theta_2, x}$ respectively. Since the adversary can only change $\epsilon$ fraction of inputs, we assume the following marginal distribution for y for $i \in \\{1, 2\\}$:
$$D_i(y | x) = \begin{cases} 1 - \epsilon, \text{ when } y = (x\theta_i)^2 \\\\ \frac{\epsilon}{\sigma}, \text{ when } y \in [ \sigma,  2\sigma] \\\\ 0, \text{ otherwise} \end{cases}$$

We want to be able to differentiate between $D_1$ and $D_2$ based on the inputs drawn from either $D_1$ or $D_2$. By reduction to a hypothesis testing problem and using Neyman-Pearson lemma:
$$ \inf_{\hat{\theta}} \sup_{\theta \in \Theta} P[|\hat{\theta} - \theta|^2 > \delta^2] \geq \frac{1}{2}(1 - TV(D_1, D_2))$$
where $TV(D_1, D_2)$ is the total variation difference between distributions $D_1$ and $D_2$.

In our illustrative example,
$$TV(D_1, D_2) = \frac{1}{2}\int_{\mathbb{R} \times \mathbb{R}} |D_1(x, y) - D_2(x, y)|$$
$$TV(D_1, D_2) = \frac{1}{2}\int_{\mathbb{R} \times \mathbb{R}} D_1(x) |D_1(y|x) - D_2(y|x)|$$

Notice that $D_1(y|x)$ and $D_2(y|x)$ can only differ when $(x\theta_1)^2 \ne (x\theta_2)^2$ and contribute $|D_1(y|x) - D_2(y|x)| \leq 2(1 - \epsilon)$ correspondingly. Overall,

$$TV(D_1, D_2) \leq 1 - \epsilon$$

It follows that,
$$ \inf_{\hat{\theta}} \sup_{\theta \in \Theta} P[|\hat{\theta} - \theta|^2 > \delta^2] \geq \frac{\epsilon}{2}$$

Note that $\epsilon$ is constant, and this minimax rate seems unavoidable to me. Am I missing something here?
Thank you for your detailed reply. I carefully read other reviews together with the authors' responses and remain in favour of acceptance for this work.  
In contrast to other reviewers, I do not find the fact that the error / sample complexity does not depend on $\varepsilon'$ surprising: note that in robust mean estimation, even when $\varepsilon = 0$ (no adversarial data), it is impossible to recover mean *exactly*, whereas in the phase retrieval it is possible, as the authors responded in one of the comments.  
Therefore, it is not counterintuitive that the algorithm proposed by the authors obtains accurate estimates independent from the corruption parameter.
Thank you for the responses. I persist in considering that the theoretical results are somewhat incremental, yet they do not appear significant enough to warrant a purely theoretical submission. Experimental validations of the proposed algorithm (e.g., numerically verifying the counter-intuitive result that there is no dependence on $\epsilon$ in both the sample complexity and the upper bound ) would be advantageous for this submission. Therefore, I am inclined to maintain my current score.
Here is a proof that information-theoretically, $O(d)$ samples are sufficient for exact recovery for phase retrieval, even under $\epsilon$-corruption.

-----

Let $G^*$ denote the original $n$ clean input pairs $(a_i, y_i)$. Each pair specifies a constraint $y_i = \langle a_i, x\rangle^2$. 

We can assume that for any subset $S \subseteq G^*$ with $|S| \ge (1-2\epsilon)n$, only the ground-truth vector $\pm x$ can satisfy all constraints in $S$. Intuitively, this holds because $x$ has $d$ degrees of freedom, and we have $(1-2\epsilon)n > d$ constraints with essentially i.i.d. Gaussian $a_i$'s.

Consider the following exponential-time algorithm: Given an $\epsilon$-corrupted set of $(a_i, y_i)$ pairs,
* Enumerate all subset $T$ of size $(1-\epsilon)n$.
* For each $T$, run a standard phase retrieval algorithm on $T$,
  * If a solution $z$ is returned and $z$ satisfies all constraints in $T$, return $z$.

First observe that this algorithm must return a solution, because $T$ will eventually hit the set of $(1-\epsilon)n$ clean pairs.

Next observe that this algorithm cannot return anything other than $\pm x$, this is because $|T| = (1-\epsilon)n$ and there are at most $\epsilon n$ corrupted pairs, so $T$ contains at least $(1-2\epsilon)n$ clean pairs. As we assumed earlier, only $\pm x$ satisfies any subset of $(1-2\epsilon)n$ clean pairs.

-----

To our knowledge, there is no existing reference for this specific information-theoretic lower bound. (Although the above exponential-time algorithm is forklore in robust statistics.) We are the first to study this problem. Our paper provides another proof that exact recovery is possible with $O(d)$ samples (and we show how to achieve it in nearly-linear time).
Thank you for the response. Could you please provide the reference containing an information-theoretic lower bound for the phase-retrieval problem under corruption?
In our work, we study the problem of robust phase retrieval, which is different from robust mean estimation.

At a high level, robust phase retrieval has an information-theoretic lower bound of error $0$ (i.e., exact recovery), whereas robust mean estimation has an information-theoretic lower bound of $O(\sqrt{\epsilon})$.

For robust phase retrieval, $O(d)$ samples are sufficient for exact recovery, even when $\epsilon$-fraction of the input is corrupted. This is possible due to the special structure of the input: $y_i = \langle a_i, x \rangle^2$. (Our contribution is to show that this can be done in nearly-linear time.)

For robust mean estimation (for bounded covariance distributions in $\mathbb{R}^d$), as the reviewer correctly noted, the information-theoretic optimal error is $O(\sqrt{\epsilon})$, and it takes $O(d/\epsilon)$ samples to achieve this error.
Thank you for your response. If you go through your suggested references, say [1, 2, 4], you will observe two important features in their theorem statements:
1. Their sample complexities depend on the amount of corruption $\epsilon$.   
2. Their convergence results depend on the amount of corruption.
For example, check Theorem 1.3.1 or D.2 of [1], where $n = \mathcal{O}(\frac{d^2}{\epsilon})$ and $\\| \hat{\mu} - \mu_D \\| \leq \sqrt{\epsilon}$. This shows a clear dependency of sample complexity and the convergence rate on the amount of corruption. Your results claim recovery up to arbitrary precision with an unknown $\epsilon$. To justify such a strong statement, in my opinion, it is important to characterize $\epsilon'$ in terms of $n$ and $d$.