# Towards Sample-Optimal Compressive Phase Retrieval with Sparse and
Generative Priors

Zhaoqiang Liu, Subhroshekhar Ghosh, Jonathan Scarlett  Z. Liu is with the
Department of Computer Science, National University of Singapore (email:
dcslizha@nus.edu.sg). S. Ghosh is with the Department of Mathematics, National
University of Singapore (email: subhrowork@gmail.com). J. Scarlett is with the
Department of Computer Science and the Department of Mathematics, National
University of Singapore (email: scarlett@comp.nus.edu.sg). This work was
supported by the Singapore National Research Foundation (NRF) under grant
R-252-000-A74-281, and the Singapore Ministry of Education (MoE) under grants
R-146-000-250-133 and R-146-000-312-114.

###### Abstract

Compressive phase retrieval is a popular variant of the standard compressive
sensing problem, in which the measurements only contain magnitude information.
In this paper, motivated by recent advances in deep generative models, we
provide recovery guarantees with order-optimal sample complexity bounds for
phase retrieval with generative priors. We first show that when using i.i.d.
Gaussian measurements and an $L$-Lipschitz continuous generative model with
bounded $k$-dimensional inputs, roughly $O(k\log L)$ samples suffice to
guarantee that the signal is close to any vector that minimizes an amplitude-
based empirical loss function. Attaining this sample complexity with a
practical algorithm remains a difficult challenge, and a popular spectral
initialization method has been observed to pose a major bottleneck. To
partially address this, we further show that roughly $O(k\log L)$ samples
ensure sufficient closeness between the signal and any globally optimal
solution to an optimization problem designed for spectral initialization
(though finding such a solution may still be challenging). We adapt this
result to sparse phase retrieval, and show that $O(s\log n)$ samples are
sufficient for a similar guarantee when the underlying signal is $s$-sparse
and $n$-dimensional, matching an information-theoretic lower bound. While our
guarantees do not directly correspond to a practical algorithm, we propose a
practical spectral initialization method motivated by our findings, and
experimentally observe significant performance gains over various existing
spectral initialization methods of sparse phase retrieval.

##  I Introduction

In this paper, we consider the (real-valued) phase retrieval problem, which
aims to recover a signal $\mathbf{x}\in\mathbb{R}^{n}$ from noisy magnitude-
only measurements:

$y_{i}=|\langle\mathbf{a}_{i},\mathbf{x}\rangle|+\eta_{i},\quad
i=1,2,\ldots,m,$ (1)

where $\mathbf{a}_{i}\in\mathbb{R}^{n}$ is the $i$-th sensing vector, and
$\eta_{i}$ is the additive noise term. This problem arises naturally in areas
such as diffraction imaging, X-ray crystallography, microscopy, optics and
astronomy [1], where it is often difficult or even impossible to observe the
linear measurements directly, and one can only record the magnitudes or
intensities (squared magnitudes).

In many real applications, to reduce the required number of measurements, it
is of interest to exploit structure in the signal being estimated. In
particular, for applications related to signal processing and imaging, it is
well-known that the underlying signal typically admits a sparse or
approximately sparse representation in some known basis [2]. Motivated by
this, and considering the popularity of the standard compressive sensing (CS)
problem [3], the sparse phase retrieval problem has attracted significant
research interest.

Moreover, inspired by the successful applications of deep generative models in
many fields [4], recently, a new perspective of CS has emerged, for which the
sparsity assumption is replaced by a generative model assumption. That is,
rather than being sparse, the signal is instead assumed to be close to the
range of a generative model [5]. In addition to the theoretical developments
in [5], the authors also provide impressive numerical results showing that for
some imaging applications, using generative priors can significantly reduce
the required number of measurements (e.g., by a factor of 5 to 10) for
recovering the signal up to a given accuracy. Follow-up works of [5] include
[6, 7, 8, 9, 10, 11], just to name a few.

In this paper, we focus on providing recovery guarantees with order-optimal
sample complexity bounds for the compressive phase retrieval problem, for
which the underlying signal is assumed to be sparse or lie in (or close to)
the range of a generative model.

###  I-A Related Work

In this subsection, we provide a summary of some relevant works, which can
roughly be divided into (i) the general phase retrieval problem with no
structural assumptions on the signal, (ii) sparse phase retrieval, and (iii)
phase retrieval with generative models.

General phase retrieval: A wide range of approaches have been designed to
solve the phase retrieval problem. Error-reduction approaches [12, 13] work
well in practice, but they lack provable guarantees. Convex methods such as
PhaseLift [14] and PhaseCut [15] lift the phase retrieval problem to a higher
dimension and suffer from high computation cost. Convex relaxations that
operate in the natural domain of the signal are proposed in [16, 17]. However,
these convex relaxation based methods are not empirically competitive against
widely-used non-convex methods.

Following the seminal work of Netrapalli et al. [18], several works have
studied theoretically-guaranteed non-convex optimization algorithms for phase
retrieval [1, 19, 20, 21]. These algorithms start from a spectral
initialization method, and then use an iterative algorithm (e.g., alternating
minimization or gradient descent) to further decrease the approximation error.
For the general phase retrieval problem, the optimal sample complexity $O(n)$
can be achieved for some non-convex methods [19, 20, 21]. Interestingly,
optimal (up to a logarithmic factor) sample complexity guarantees have also
been provided for the case of using a random initialization [22, 23]. It is
worth mentioning that some variants of the classic Wirtinger Flow (WF)
algorithm [1, 19] and trust-region methods [22] consider squared measurements:

$y_{i}=|\langle\mathbf{a}_{i},\mathbf{x}\rangle|^{2}+\bm{\epsilon}_{i},\quad
i=1,2,\ldots,m,$ (2)

and minimize the intensity based empirical loss function:

$f_{I}(\mathbf{w}):=\frac{1}{2}\left\|\mathbf{y}-|\mathbf{A}\mathbf{w}|^{2}\right\|_{2}^{2},$
(3)

where $\mathbf{A}\in\mathbb{R}^{m\times n}$ is the sensing matrix with its
$i$-th row being $\mathbf{a}_{i}^{T}$, and $|\cdot|$ is applied element-wise.
Nonetheless, numerical results suggest that algorithms minimizing the
amplitude based empirical risk (with the measurement model (1))

$f(\mathbf{w}):=\|\mathbf{y}-|\mathbf{A}\mathbf{w}|\|_{2}^{2}$ (4)

are usually more efficient in computation [20, 21, 24]. Based on this, we
focus on the measurement model (1) and the associated loss function (4).

Sparse phase retrieval: Various algorithms have also been devised for sparse
phase retrieval [25, 26, 27]. In particular, there exist various non-convex
optimization based algorithms, including thresholded/projected WF [28, 24],
sparse truncated amplitude flow [29], compressive phase retrieval with
alternating minimization [30], and sparse phase retrieval by hard iterative
pursuit [31]. All of these approaches are analyzed under the assumption of
using a sensing matrix with i.i.d. real Gaussian entries. In addition, similar
to non-convex methods for the general phase retrieval problem, all of these
approaches start with a spectral initialization step, and then use an
iterative algorithm to refine the initial vector. More specifically, in the
initialization step, these approaches try to accurately estimate the support
of the underlying signal.

Unlike the general phase retrieval setting in which optimal sample complexity
guarantees are provided for some practical algorithms, for non-convex methods
of sparse phase retrieval, the typical sample complexity requirement is
$O(s^{2}\log n)$, where $s$ is the sparsity level. This does not match the
information theoretic lower bound [32, 33]; the bottleneck comes from the
accurate estimation of the support set [29, 24, 30]. Improved guarantees are
also available in the noiseless setting under more restrictive assumptions on
the magnitudes of the non-zero entries [34]. In addition, for sparse phase
retrieval, recovery guarantees with respect to globally optimal solutions of
both the intensity and amplitude based loss functions have been provided in
[35, 36, 37, 38, 39], achieving the optimal sample complexity $O(s\log n)$.
However, designing a computationally efficient algorithm achieving the optimal
sample complexity for $s$-sparse vectors remains an important open problem.

Phase retrieval with generative models: Phase retrieval using generative
models has been studied in [40, 41, 42, 43, 44]. In particular, the authors of
[40] model natural signals as being in the range of a fully-connected ReLU
neural network with random weights and no offsets, and study the global
geometry for gradient methods. Extensive numerical experiments for phase
retrieval with generative models, as well as both Gaussian and coded
diffraction pattern measurements, have been presented in [43]. Algorithms
proposed in both works [40, 43] minimize the objective function directly over
the latent variable $\mathbf{z}\in\mathbb{R}^{k}$ using gradient descent,
where $k$ is the latent dimension. The corresponding objective function is
highly non-convex, and performing gradient descent directly over $\mathbf{z}$
limits the explorable solution space, which may lead to getting stuck in local
minima.

Recovery guarantees for non-convex optimization algorithms for phase retrieval
with pre-trained ReLU neural network or un-trained neural network generative
models have been provided in [41, 42], with assumptions of noiseless
measurements and the existence of a good initial vector. Recovery guarantees
for phase retrieval with generative models have been provided in [44], within
a more general setting of using differentiable but unknown link functions. The
weaker assumption of an unknown link function comes as a price. In particular,
for an unknown nonlinearity, it is an open problem to handle signals with
representation error [45]. Moreover, the assumption of differentiability makes
it is not directly applicable to the measurement model (1), and the loss
function considered in [44] is different from both (3) and (4), which are
widely adopted for phase retrieval.

###  I-B Contributions

  * •

We show that for compressive phase retrieval using i.i.d. Gaussian
measurements, when the signal is close to the range of an $L$-Lipschitz
continuous generative model with bounded $k$-dimensional inputs, roughly
$O(k\log L)$ samples suffice for ensuring the closeness between the signal and
any vector minimizing the amplitude-based empirical loss function (4).

  * •

To address the sample complexity barrier posed by spectral initialization, we
propose a relevant optimization problem, and show that roughly $O(k\log L)$
samples are sufficient to guarantee that any globally optimal solution of the
optimization problem is close to the true signal. This suggests the
plausibility of practical spectral initialization algorithms that are able to
find such a global optimum, though we do not attempt to provide one that can
provably do so.

  * •

We adapt our analysis to the case of sparse phase retrieval, and show that
roughly $O(s\log n)$ samples suffice for a similar recovery guarantee, where
$s$ is the sparsity level and $n$ is the ambient dimension. This matches the
information-theoretic lower bound.

  * •

Motivated by our theoretical findings for sparse phase retrieval, we propose a
spectral initialization method based on sparse principal component analysis.
We verify on synthetic experiments that our initialization method
significantly outperforms several popular spectral initialization methods of
sparse phase retrieval (in the sense of relative error). This further
corroborates the optimal sample complexity guarantee in our theory.

###  I-C Notation

We use upper and lower case boldface letters to denote matrices and vectors
respectively. We write $[N]=\\{1,2,\cdots,N\\}$ for a positive integer $N$. We
define the $\ell_{2}$-ball
$B_{2}^{k}(r):=\\{\mathbf{z}\in\mathbb{R}^{k}:\|\mathbf{z}\|_{2}\leq r\\}$,
and the unit sphere
$\mathcal{S}^{n-1}:=\\{\mathbf{x}\in\mathbb{R}^{n}:\|\mathbf{x}\|_{2}=1\\}$.
We use $G$ to denote an $L$-Lipschitz continuous generative model from
$B_{2}^{k}(r)$ to $\mathbb{R}^{n}$. For a set $B\subseteq B_{2}^{k}(r)$, we
write $G(B)=\\{G(\mathbf{z})\,:\,\mathbf{z}\in B\\}$. The sensing matrix
$\mathbf{A}\in\mathbb{R}^{m\times n}$ is assumed to have i.i.d. standard
normal entries, i.e., $a_{ij}\overset{i.i.d.}{\sim}\mathcal{N}(0,1)$. For any
$\mathbf{s}\in\mathbb{R}^{n}$, we use
$\bar{\mathbf{s}}=\frac{\mathbf{s}}{\|\mathbf{s}\|_{2}}$ to denote the
corresponding normalized vector. The support (set) of a vector is the index
set of its non-zero entries. For any $\mathbf{X}\in\mathbb{R}^{m\times n}$ and
any index set $I\subseteq[m]$, $\mathbf{X}_{I}$ denotes the $|I|\times n$ sub-
matrix of $\mathbf{X}$ that only keeps the rows indexed by $I$. We use
$\|\mathbf{X}\|_{2\to 2}$ to denote the spectral norm of a matrix
$\mathbf{X}$. For any positive integer $N$, we use $\mathbf{I}_{N}$ to denote
the identity matrix in $\mathbb{R}^{N\times N}$. We use the generic notations
$C$ and $C^{\prime}$ to denote large positive constants, and we use $c$ and
$c^{\prime}$ to denote small positive constants. Their values may differ from
line to line.

##  II Optimal Solutions of Amplitude-Based Loss Function for Generative
Priors

In this section, we provide recovery guarantees with respect to optimal
solutions of the amplitude based loss function (4) for phase retrieval with
generative models. We let the generative model
$G\,:\,B_{2}^{k}(r)\to\mathbb{R}^{n}$ be any $L$-Lipschitz continuous
function, and we suppose that the underlying signal $\mathbf{x}$ is close to
(but does not need to lie in) the range of $G$. We are interested in the case
$k\ll n$ (i.e., relatively small latent dimension). The sensing matrix
$\mathbf{A}$ is assumed to have i.i.d. standard normal entries. Let
$\mathbf{q}\in\mathrm{Range}(G)$ minimize (4) to within additive $\tau>0$ of
the optimum, i.e.,

$\|\mathbf{y}-|\mathbf{A}\mathbf{q}|\|_{2}^{2}\leq\min_{\mathbf{w}\in\mathrm{Range}(G)}\|\mathbf{y}-|\mathbf{A}\mathbf{w}|\|_{2}^{2}+\tau,$
(5)

and define $\mathbf{p}\in\mathrm{Range}(G)$ to be the point that is closest to
the signal $\mathbf{x}$, i.e.,

$\mathbf{p}=\arg\min_{\mathbf{w}\in\mathrm{Range}(G)}\|\mathbf{x}-\mathbf{w}\|_{2}.$
(6)

With the above settings in place, we present the following theorem, proved in
Appendix B.

######  Theorem 1.

Consider the observation model (1) with i.i.d. Gaussian measurements, and
$\mathbf{p},\mathbf{q}$ satisfying (5)–(6) with parameter $\tau$. Then, for
any $\delta>0$, we have that if $m=\Omega(k\log\frac{Lr}{\delta})$,111Here and
subsequently, all $m=\Omega(\cdot)$ statements are assumed to have a large
enough implied constant. with probability $1-e^{-\Omega(m)}$, it holds that

$\min\\{\|\mathbf{q}-\mathbf{x}\|_{2},\|\mathbf{q}+\mathbf{x}\|_{2}\\}\leq
O\left(\|\mathbf{p}-\mathbf{x}\|_{2}+\frac{\|\bm{\eta}\|_{2}+\sqrt{\tau}}{\sqrt{m}}+\delta\right).$
(7)

Since we assume that $\mathbf{x}$ is close to the range of $G$, the
representation error $\|\mathbf{p}-\mathbf{x}\|_{2}$ is small (compared to
$\|\mathbf{x}\|_{2}$). The optimization error $\tau$ is also ideally small. In
addition, it is common to assume that $\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c\|\mathbf{x}\|_{2}$ (see, e.g., [21]), where $c>0$ is also small. A typical
$d$-layer neural network generative model has $\mathrm{poly}(n)$-bounded
weights in each layer, and thus its Lipschitz constant is $L=n^{O(d)}$ [5],
where $d$ is the depth. Then, the values of $r$ and $\delta$ can be as extreme
as $r=n^{O(d)}$ and $\delta=\frac{1}{n^{O(d)}}$ without affecting the final
scaling. In this case, Theorem 1 reveals that roughly $O(k\log L)$ samples
suffice to guarantee that any optimal solution of the amplitude-based loss
function (4) is close to the true signal (up to an unavoidable possible global
sign flip). By the analysis of information-theoretic lower bounds in [46, 47],
this scaling cannot be improved without further assumptions.

##  III Spectral Initializations with Generative Priors

In Theorem 1, we provided recovery guarantees for phase retrieval with
generative models with respect to optimal solutions of the amplitude based
loss function. However, it is not clear whether such a guarantee is achievable
by practical algorithms.

A notable recent work is [41], showing that for a depth-$d$, width-$w$ ReLU
neural network, roughly $O(kd\log w)$ samples222This matches the $O(k\log L)$
scaling upon substituting $L=O\big{(}w^{d}\big{)}$ for ReLU networks [5].
suffice for their alternating phase projected gradient descent (APPGD)
approach to obtain an accurate estimate of the signal, as long as there exists
a proper initial vector $\mathbf{x}^{(0)}$, in the sense that
$\min\\{\|\mathbf{x}-\mathbf{x}^{(0)}\|_{2},\|\mathbf{x}+\mathbf{x}^{(0)}\|_{2}\\}\leq
c\|\mathbf{x}\|_{2}$ for some small positive constant $c$. However, as we know
for sparse phase retrieval, designing a practical algorithm to find a proper
initial vector with optimal sample complexity requirements is more difficult
than designing the subsequent iterative algorithm that refines the initial
guess, and it still remains as an open problem.

Existing initialization methods for sparse phase retrieval typically first
estimate the support of the sparse signal, and then perform a power method on
a sub-matrix corresponding to the estimated support set to calculate an
initial vector. Without extra assumptions on the structure of the sparse
signal, these methods have a sub-optimal sample complexity of $O(s^{2}\log
n)$. For phase retrieval with generative models, the situation is even worse –
there is no theoretically-guaranteed initialization method at all. To address
this gap, we provide recovery guarantees for a spectral initialization method
of phase retrieval with generative models.

We assume that $\mathbf{x}$ is in the range of
$G\,:\,B_{2}^{k}(r)\to\mathbb{R}^{n}$ (again assumed to be $L$-Lipschitz
continuous), and we assume that the noise terms
$\eta_{1},\eta_{2},\ldots,\eta_{m}$ are bounded as follows:

$\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq c_{0}\|\mathbf{x}\|_{2},$ (8)

where $c_{0}>0$ is a sufficiently small constant. This assumption states that
the signal to noise ratio (SNR) is sufficiently large, and has also been made
in relevant existing works (see, e.g., [21, Theorem 3]). Moreover, we assume
that

$\|\bm{\eta}\|_{\infty}\leq c_{1}\|\mathbf{x}\|_{2},$ (9)

where $c_{1}>0$ is also a sufficiently small constant. This assumption is
satisfied, for example, when
$\|\bm{\eta}\|_{\infty}=O\big{(}\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\big{)}$
(assuming (8)), i.e., when none of the noise entries are unusually large.
Again, similar assumptions have been made in existing works such as [19,
Theorem 2], and [48, Theorem 2].333In these papers, the authors consider an
intensity based measurement model, and accordingly, the corresponding
assumption is for $\frac{\|\bm{\eta}\|_{\infty}}{\|\mathbf{x}\|_{2}^{2}}$
instead of for $\frac{\|\bm{\eta}\|_{\infty}}{\|\mathbf{x}\|_{2}}$.

Let $\lambda$ be defined as

$\lambda:=\sqrt{\frac{\pi}{2}}\cdot\frac{1}{m}\sum_{i=1}^{m}y_{i}.$ (10)

Using sub-Gaussian concentration [49, Proposition 5.10], $\lambda$ is close to
$\|\mathbf{x}\|_{2}$ with high probability given enough samples. In the
following, we focus on estimating the normalized signal vector
$\bar{\mathbf{x}}:=\frac{\mathbf{x}}{\|\mathbf{x}\|_{2}}$. For this purpose,
similar to the idea in [50], we consider a normalized generative model
$\tilde{G}\,:\,\mathcal{D}\to\mathcal{S}^{n-1}$, where
$\mathcal{D}:=\\{\mathbf{z}\in
B_{2}^{k}(r)\,:\,\|G(\mathbf{z})\|_{2}>R_{\min}\\}$ for some
$R_{\min}>0$,444As discussed in [50], the dependence on $R_{\min}$ in the
sample complexity is very mild. Under the typical scaling $L=n^{O(d)}$ for a
$d$-layer neural network, the scaling laws remain unchanged even with
$R_{\min}=\frac{1}{n^{O(d)}}$. $\mathcal{S}^{n-1}$ denotes the unit sphere in
$\mathbb{R}^{n}$, and
$\tilde{G}(\mathbf{z})=\frac{G(\mathbf{z})}{\|G(\mathbf{z})\|_{2}}$. Then,
$\tilde{G}$ is $\tilde{L}$-Lipschitz continuous with
$\tilde{L}=\frac{L}{R_{\min}}$.

Inspired by the initialization method used in [21] for the general phase
retrieval problem, we consider the following matrix defined with a truncation
operation (see Remark 2 below for a discussion comparing to the non-truncated
counterpart):

$\mathbf{V}:=\frac{1}{m}\sum_{i=1}^{m}y_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l\lambda<y_{i}<u\lambda\\}},$
(11)

where $l,u$ are positive constants satisfying $1<l<u$, and $\lambda$ is given
in (10). Note that for the noiseless case with
$y_{i}=|\langle\mathbf{a}_{i},\mathbf{x}\rangle|$, in accordance with the
above-mentioned closeness between $\lambda$ and $\|\mathbf{x}\|_{2}$, it is
useful to consider the following expectation:

$\mathbf{J}:=\mathbb{E}\left[\frac{1}{m}\sum_{i=1}^{m}y_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l\|\mathbf{x}\|_{2}<y_{i}<u\|\mathbf{x}\|_{2}\\}}\right]=\|\mathbf{x}\|_{2}(\gamma_{0}\mathbf{I}_{n}+\beta_{0}\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}),$
(12)

where for $g\sim\mathcal{N}(0,1)$, we define
$\gamma_{0}:=\mathbb{E}[|g|\mathbf{1}_{\\{l<|g|<u\\}}]$ and
$\beta_{0}:=\mathbb{E}\big{[}|g|^{3}\mathbf{1}_{\\{l<|g|<u\\}}\big{]}-\gamma_{0}$
(cf. Lemma 8 in Appendix C-B). Based on the observation that
$\bar{\mathbf{x}}$ is a leading eigenvector of $\mathbf{J}$ and the assumption
$\mathbf{x}\in\mathrm{Range}(G)$, we consider using the following optimization
problem to find an $\hat{\mathbf{x}}$ that approximates $\bar{\mathbf{x}}$:

$\hat{\mathbf{x}}:=\arg\max_{\mathbf{w}\in\tilde{G}(\mathcal{D})}\quad\mathbf{w}^{T}\mathbf{V}\mathbf{w}.$
(13)

In addition, in order to take the norm of $\mathbf{x}$ into account, we set
the the initial vector $\mathbf{x}^{(0)}$ as

$\mathbf{x}^{(0)}=\lambda\hat{\mathbf{x}}.$ (14)

We show that roughly $O(k\log L)$ samples suffice for ensuring the closeness
between $\mathbf{x}^{(0)}$ and $\mathbf{x}$ (up to a global signal flip).
Formally, we have the following.

######  Theorem 2.

Consider the model (1) with i.i.d. Gaussian measurements and noise satisfying
(8)–(9), and let $\hat{\mathbf{x}}$ be as defined in (13). Suppose that555This
is weaker than the assumption that $\|\mathbf{x}\|_{2}>R_{\min}$ made in [50].
Moreover, even if $\|\mathbf{x}\|_{2}\leq R_{\min}$, the condition
$\bar{\mathbf{x}}\in\tilde{G}(\mathcal{D})$ may be satisfied (e.g., when $G$
is a ReLU neural network generative model with no offsets, then the range of
$G$ becomes a cone, and if $0<\|\mathbf{x}\|_{2}\leq R_{\min}$, we can choose
a sufficiently large $C$ such that $C\|\mathbf{x}\|_{2}>R_{\min}$, with
$C\mathbf{x}\in{\rm Range}(G)$). $\bar{\mathbf{x}}\in\tilde{G}(\mathcal{D})$
and $u>l>1+c_{1}$, where
$c_{1}\geq\frac{\|\bm{\eta}\|_{\infty}}{\|\mathbf{x}\|_{2}}$ appears in (9).
Given a sufficiently small positive constant $c$, when
$m=\Omega\big{(}k\log(\tilde{L}nr)\big{)}$, we have with probability
$1-e^{-\Omega(m)}$ that

$\min\\{\|\hat{\mathbf{x}}-\bar{\mathbf{x}}\|_{2},\|\hat{\mathbf{x}}+\bar{\mathbf{x}}\|_{2}\\}<0.9c.$
(15)

In addition, we have with probability $1-e^{-\Omega(m)}$ that
$\mathbf{x}^{(0)}=\lambda\hat{\mathbf{x}}$ satisfies

$\min\\{\|\mathbf{x}^{(0)}-\mathbf{x}\|_{2},\|\mathbf{x}^{(0)}+\mathbf{x}\|_{2}\\}<c\|\mathbf{x}\|_{2}.$
(16)

The proof is given in Appendix C. Using well-established chaining arguments as
those in [5, 50, 51], the sample complexity $\Omega(k\log(\tilde{L}nr))$ can
be reduced to $\Omega(k\log(\tilde{L}r))$. However, since for a typical
$d$-layer neural network, we have the Lipschitz constant $L=n^{O(d)}$ (and
thus $\tilde{L}$ and $r$ can also be of order $n^{O(d)}$) [5], such a
reduction is typically of minor importance.

##  IV Sparse Phase Retrieval

The proof of Theorem 2 relies on the fact that $\bar{\mathbf{x}}$ lies in a
subset of the unit sphere with the $\delta$-covering number being upper
bounded by $\big{(}\frac{\tilde{L}r}{\delta}\big{)}^{k}$. Letting
$\Sigma_{s}^{n}$ is the set of all $s$-sparse unit vectors in
$\mathbb{R}^{n}$, we know that the $\delta$-covering number of
$\Sigma_{s}^{n}$ is upper bounded by
$\binom{n}{s}\big{(}\frac{1}{\delta}\big{)}^{s}\leq\big{(}\frac{en}{\delta
s}\big{)}^{s}$ [52]. Based on this observation, we can readily adapt the
analysis in Theorem 2 and obtain the following theorem.

######  Theorem 3.

Consider the model (1) with i.i.d. Gaussian measurements and noise satisfying
(8)–(9), and suppose that $\mathbf{x}\in\mathbb{R}^{n}$ is $s$-sparse. Let
$\lambda$ be defined as in (10). Suppose that $u>l>1+c_{1}$, where
$c_{1}\geq\frac{\|\bm{\eta}\|_{\infty}}{\|\mathbf{x}\|_{2}}$ appears in (9),
and $\mathbf{V}$ is defined in (11). Let $\hat{\mathbf{x}}$ be defined as

$\hat{\mathbf{x}}:=\arg\max_{\mathbf{w}\in\Sigma_{s}^{n}}\mathbf{w}^{T}\mathbf{V}\mathbf{w}.$
(17)

Then, when $m=\Omega\big{(}s\log n\big{)}$, we have with probability
$1-e^{-\Omega(m)}$ that

$\min\\{\|\hat{\mathbf{x}}-\bar{\mathbf{x}}\|_{2},\|\hat{\mathbf{x}}+\bar{\mathbf{x}}\|_{2}\\}<0.9c,$
(18)

where $c>0$ is a sufficiently small constant. In addition, we have with
probability $1-e^{-\Omega(m)}$ that $\mathbf{x}^{(0)}=\lambda\hat{\mathbf{x}}$
satisfies

$\min\\{\|\mathbf{x}^{(0)}-\mathbf{x}\|_{2},\|\mathbf{x}^{(0)}+\mathbf{x}\|_{2}\\}<c\|\mathbf{x}\|_{2}.$
(19)

Implications regarding sparse PCA and support recovery: We note that (17) is
essentially a problem of sparse principal component analysis (SPCA) [53, 54].
Existing spectral initialization methods for sparse phase retrieval typically
start from estimating the support of $\mathbf{x}$, and an accurate estimate of
the support requires a sub-optimal $O(s^{2}\log n)$ sample complexity [28, 29,
30], which is the major bottleneck of the sample complexity upper bound of the
whole procedure (i.e., the initialization step and the subsequent iterative
step). It is unknown whether there is a practical and sample-optimal approach
for estimating the support. In fact, even the following question has not been
answered previously:

Key question: Does there exist any spectral-based method (computationally
efficient or otherwise) attaining an accurate initial vector with order-
optimal sample complexity?

We provide a positive answer to this question, showing that $O(s\log n)$
samples suffice to ensure any optimal solution of the SPCA problem (17) is
close to the underlying (normalized) signal.

Motivated by the strong sample complexity guarantee provided in Theorem 3, we
propose a spectral initialization method for sparse phase retrieval as
described in Algorithm 1. In particular, after calculating $\mathbf{V}$ from
$\mathbf{A}$ and $\mathbf{y}$, we use algorithms for SPCA to solve (17). We
observe via numerical experiments in Section V that our initialization method
significantly outperforms several popular spectral initialization methods used
for sparse phase retrieval.

######  Remark 1.

In various works on SPCA, it has been shown that there exists a statistical-
computational gap for achieving consistency and optimal convergence rates
uniformly over a parameter space ( e.g., see [55, 56, 28]). However, the
models therein are all different from ours (e.g., the spiked covariance model
[57]), so it is unclear to what extent we should expect similar computational
difficulties.

######  Remark 2.

The truncation function $\bm{1}_{\\{l\lambda<y_{i}<u\lambda\\}}$ in (11) plays
an important role in reducing the required sample complexity in our analysis.
Without the truncation function, if we use the same proof technique to derive
a sample complexity bound, we need to make use of concentration inequalities
for random variables that have heavier tail than sub-exponential random
variables, e.g., sub-Weibull random variables of order $\alpha<1$ (see, e.g.,
[58, Theorem 3.1]). This leads to a larger sample complexity bound. In Section
V, we present numerical results for applying SPCA to the non-truncated
weighted empirical matrix

$\tilde{\mathbf{V}}=\frac{1}{m}\sum_{i=1}^{m}y_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{T}.$
(20)

We will see in our experiments that the method corresponding to the truncated
empirical matrix $\mathbf{V}$ significantly outperforms the method
corresponding to the non-truncated empirical matrix $\tilde{\mathbf{V}}$,
although the expectations of these two matrices have a very similar structure.
Thus, we expect that truncation has the potential of overcoming the
computational challenges discussed above.

We note that [59, 60] consider solving sparse phase retrieval via practical
SPCA methods, in a more general setting considering an unknown link function.
They consider a weighted empirical covariance matrix with no truncation, and
the derived sample complexity has a quadratic dependence on $s$.

##  V Numerical Experiments

Algorithm 1 A spectral initialization method for sparse phase retrieval based
on SPCA (PRI-SPCA)

Input: $\mathbf{A}$, $\mathbf{y}$, $l$, $u$

Output: An initial vector $\mathbf{x}^{(0)}$

1) Calculate $\lambda$ from (10) and then calculate $\mathbf{V}$ from (11).

2) Use an SPCA algorithm to compute $\hat{\mathbf{x}}$ according to (17).

3) Let $\mathbf{x}^{(0)}=\lambda\hat{\mathbf{x}}$.

In this section, we perform synthetic experiments to complement our
theoretical results. We focus only on the sparse setting, since numerous
existing algorithms are known for sparse priors but not for generative priors.
Specifically, we compare our method (PRI-SPCA; see Algorithm 1) with spectral
initialization methods used in several popular algorithms for sparse phase
retrieval: Thresholded Wirtinger flow (ThWF) [28], sparse truncated amplitude
flow (SPARTA) [29], and compressive phase retrieval with alternating
minimization (CoPRAM) [30].666We use the MATLAB package shared by the authors
of [30] at https://github.com/GauriJagatap/model-copram. The package not only
contains the codes for CoPRAM, but also for ThWF and SPARTA. We also compare
with the approach of solving an SPCA problem with respect to
$\tilde{\mathbf{V}}$ (cf. (20)), where $\tilde{\mathbf{V}}$ is a non-truncated
version of $\mathbf{V}$. The corresponding method is termed PRI-SPCA-NT, with
‘NT’ standing for no truncation. For PRI-SPCA, we set $l=1$ and $u=5$, as
previously suggested in [21].

Various algorithms have been developed to tackle the SPCA problem, see, e.g.,
[61, 62, 63]. In our experiments, we present the numerical results using the
generalized Rayleigh quotient iteration (GRQI) method [64], but we observed
that other popular SPCA methods such as the truncated power method (TPower)
[65] give similar averaged relative errors.777Although using TPower also
results in small relative error for our PRI-SPCA method, in our experiments,
as we want a fair comparison and do not carefully design the starting point of
TPower, it is not as stable as GRQI. For GRQI, we set the deflation parameter
to be $0.2$, and the total number of iterations to be $100$. Since the
sparsity level $s$ is typically assumed to be known a priori for SPARTA and
CoPRAM (ThWF assumes the knowledge of a parameter that plays the similar role
as $s$), for GRQI, we set the parameter about the maximum number of non-zero
indices as $s$. Note that the time complexity of each iteration of GRQI is
$O(ns+s^{3})$, which scales well with the subsequent $O(s^{2}n\log n)$ time
complexity of the iterative methods used in the three approaches (for ThWF,
the corresponding time complexity is $O(n^{2}\log n)$. See [30, Table I]). The
number of power method steps used in all the methods is fixed to be $100$.

The authors of [64] suggest using the column of maximal norm as a starting
point of GRQI. Because this policy has a time complexity of $O(n^{2})$, we
instead choose the column corresponding to the largest diagonal entry as the
starting point. Considering the structure of $\mathbf{V}$ (cf. (12)), this
modified strategy is in fact identical to the suggested policy, but with a
smaller time complexity requirement. For a fair comparison, all power methods
used in ThWF, SPARTA, and CoPRAM are similarly initialized. For the
initialization method used in ThWF, we set the tuning parameter $\alpha=0.1$,
as suggested in [28]. Similarly, we set
$|\bar{\mathcal{I}}|=\lceil\frac{1}{6}m\rceil$ for the initialization method
of SPARTA.

The support of the signal vector $\mathbf{x}$ is uniformly random, and the
non-zero entries are generated from i.i.d. standard normal distributions. The
noise vector $\bm{\eta}$ is generated from the distribution
$\mathcal{N}(\bm{0},\sigma^{2}\|\mathbf{x}\|_{2}^{2}\mathbf{I}_{m})$, where
$\sigma>0$ is a parameter indicating the noise level. The initial relative
error is defined as

$\frac{\min\\{\|\mathbf{x}-\mathbf{x}^{(0)}\|_{2},\|\mathbf{x}+\mathbf{x}^{(0)}\|_{2}\\}}{\|\mathbf{x}\|_{2}},$
(21)

where $\mathbf{x}^{(0)}$ is an initial vector produced by each of the above
mentioned initialization methods. The reported relative error is averaged over
$50$ random trials. All numerical experiments are conducted using MATLAB
R2014b on a computer with an Intel Core i5 CPU at 1.8 GHz and 8GB RAM.

  1. 1.

Sample size effect: In our first experiment, similarly to [29], we consider
noiseless measurements, and we fix the signal dimension as $n=1000$. We
consider $s$ equaling both $10$ and $20$. The number of measurements $m$
increases in $\\{100,200,\dotsc,3000\\}$. From Figure 3, we observe that our
PRI-SPCA method almost always achieves the smallest initial relative error,
thus outperforming all other methods, including the PRI-SPCA-NT method, which
uses the non-truncated empirical matrix $\tilde{\mathbf{V}}$.

  2. 2.

Sparsity effect: We consider noiseless measurements and fix $n=1000$, and set
$m=1000$ or $2000$. We vary $s$ in $\\{5,10,\dotsc,50\\}$. We observe from
Figure 3 that for this experiment, our PRI-SPCA method also gives the best
numerical performance.

  3. 3.

Noise effect: To demonstrate that PRI-SPCA also performs well in the noisy
case, we vary $\sigma$ in $\\{0.1,0.2,\dotsc,1\\}$. In addition, similarly to
[29], we fix $n=1000$, $m=3000$, and consider $s$ equaling both $10$ and $20$.
For this experiment, similarly to that in [30], we do not compare with ThWF,
as it uses quadratic measurements. From Figure 3, we observe that PRI-SPCA
outperforms SPARTA and CoPRAM, and it is better than PRI-SPCA-NT at low noise
levels.

In Appendix E, we additionally compare with other initialization methods in
terms of the effectiveness when combining with the subsequent iterative
algorithm of CoPRAM. In particular, we compare the relative error and
empirical success rate for the noiseless case, and for the noisy case, we
compare the relative error when using approximately the same running time. We
observe that for most cases, our PRI-SPCA method is superior to all other
methods we compare with. This again suggests that PRI-SPCA can serve as a good
initialization method in practice.

(a)  (b)

Figure 1: The relation between the average relative error and the number of
measurements $m$.

(a)  (b)

Figure 2: The relation between the average relative error and the sparsity
level $s$.

(a)  (b)

Figure 3: The relation between the average relative error and the noise level
$\sigma$.

##  VI Conclusion and Future Work

We have provided an order-optimal recovery guarantee for amplitude-based loss
minimization for phase retrieval with generative priors. In addition,
motivated by the bottleneck of the spectral initialization for compressive
phase retrieval, we have provided order-optimal recovery guarantees with
respect to the optimal solutions of optimization problems designed for
spectral initialization.

An immediate future research direction is to design provably sample-optimal
and computationally efficient spectral initialization methods that build on
our theoretical results. In addition, extensions to complex-valued phase
retrieval would be of significant interest.

##  Appendix A General Auxiliary Results

In this section, we provide some useful auxiliary results that will be used
throughout our analysis. First, we have the following basic concentration
inequality.

######  Lemma 1.

([66, Lemma 1.3]) Fix fixed $\mathbf{x}\in\mathbb{R}^{n}$, we have for any
$\epsilon\in(0,1)$ that

$\displaystyle\mathbb{P}\left((1-\epsilon)\|\mathbf{x}\|_{2}^{2}\leq\Big{\|}\frac{1}{\sqrt{m}}\mathbf{A}\mathbf{x}\Big{\|}_{2}^{2}\leq(1+\epsilon)\|\mathbf{x}\|_{2}^{2}\right)\geq
1-2e^{-\epsilon^{2}(1-\epsilon)m/4}.$ (22)

Next, we state the following standard definition.

######  Definition 1.

A random variable $X$ is said to be sub-Gaussian if there exists a positive
constant $C$ such that $\left(\mathbb{E}\left[|X|^{p}\right]\right)^{1/p}\leq
C\sqrt{p}$ for all $p\geq 1$. The sub-Gaussian norm of a sub-Gaussian random
variable $X$ is defined as $\|X\|_{\psi_{2}}:=\sup_{p\geq
1}p^{-1/2}\left(\mathbb{E}\left[|X|^{p}\right]\right)^{1/p}$.

According to [49, Proposition 5.10], we have the following concentration
inequality for sub-Gaussian random variables.

######  Lemma 2.

(Hoeffding-type inequality [49, Proposition 5.10]) Let $X_{1},\ldots,X_{N}$ be
independent zero-mean sub-Gaussian random variables, and let
$K=\max_{i}\|X_{i}\|_{\psi_{2}}$. Then, for any
$\bm{\alpha}=[\alpha_{1},\alpha_{2},\ldots,\alpha_{N}]^{T}\in\mathbb{R}^{N}$
and any $t\geq 0$, it holds that

$\mathbb{P}\left(\Big{|}\sum_{i=1}^{N}\alpha_{i}X_{i}\Big{|}\geq
t\right)\leq\exp\left(1-\frac{ct^{2}}{K^{2}\|\bm{\alpha}\|_{2}^{2}}\right),$
(23)

where $c>0$ is a constant.

Alongside the sub-Gaussian notion in Definition 1, we use the following
definition of a sub-exponential random variable and sub-exponential norm.

######  Definition 2.

A random variable $X$ is said to be sub-exponential if there exists a positive
constant $C$ such that
$\left(\mathbb{E}\left[|X|^{p}\right]\right)^{\frac{1}{p}}\leq Cp$ for all
$p\geq 1$. The sub-exponential norm of $X$ is defined as

$\|X\|_{\psi_{1}}=\sup_{p\geq
1}p^{-1}\left(\mathbb{E}\left[|X|^{p}\right]\right)^{\frac{1}{p}}.$ (24)

We have the following concentration inequality for sums of independent sub-
exponential random variables.

######  Lemma 3.

(Bernstein-type inequality [49, Proposition 5.16]) Let $X_{1},\ldots,X_{N}$ be
independent zero-mean sub-exponential random variables, and
$K=\max_{i}\|X_{i}\|_{\psi_{1}}$. Then for every
$\bm{\alpha}=[\alpha_{1},\ldots,\alpha_{N}]^{T}\in\mathbb{R}^{N}$ and
$\epsilon\geq 0$, it holds that

$\mathbb{P}\bigg{(}\Big{|}\sum_{i=1}^{N}\alpha_{i}X_{i}\Big{|}\geq\epsilon\bigg{)}\leq
2\exp\left(-c\cdot\mathrm{min}\Big{(}\frac{\epsilon^{2}}{K^{2}\|\bm{\alpha}\|_{2}^{2}},\frac{\epsilon}{K\|\bm{\alpha}\|_{\infty}}\Big{)}\right).$
(25)

##  Appendix B Proof of Theorem 1 (Recovery Guarantee for Amplitude-Based Loss
Minimization)

Before proving the theorem, we present some additional auxiliary results.

###  B-A Useful Lemmas

First, based on a basic two-sided concentration bound for standard Gaussian
matrices (cf., Lemma 1 in Appendix A), and the well-established chaining
arguments used in [5, 51], we have the following lemma, which essentially
gives a two-sided Set-Restricted Eigenvalue Condition (S-REC) [5].

######  Lemma 4.

[5, Lemma 4.1], [51, Lemma 2] For $\alpha<1$ and $\delta>0$, if
$m=\Omega\big{(}\frac{k}{\alpha^{2}}\log\frac{Lr}{\delta}\big{)}$, then with
probability $1-e^{-\Omega(\alpha^{2}m)}$, we have for all
$\mathbf{x}_{1},\mathbf{x}_{2}\in G(B_{2}^{k}(r))$ that

$(1-\alpha)\|\mathbf{x}_{1}-\mathbf{x}_{2}\|_{2}-\delta\leq\frac{1}{\sqrt{m}}\|\mathbf{A}(\mathbf{x}_{1}-\mathbf{x}_{2})\|_{2}\leq(1+\alpha)\|\mathbf{x}_{1}-\mathbf{x}_{2}\|_{2}+\delta.$
(26)

In addition, we have the following lemma, which is similar to Lemma 1 in
Appendix A.

######  Lemma 5.

[67, Lemma 4.3] Suppose that $X_{1},X_{2},\ldots,X_{N}$ are i.i.d. standard
normal random variables. Then, for $0\leq\epsilon\leq\mu$, we have with
probability $1-e^{-\Omega(m\epsilon^{2})}$ that

$\sqrt{m}(\mu-\epsilon)<\sqrt{\sum_{i=1}^{\lceil
m/2\rceil}|X|_{(i)}^{2}}<\sqrt{m}(\mu+\epsilon),$ (27)

where $\mu\geq\frac{1}{18}\sqrt{\frac{\pi}{2}}$ is a positive constant, and
$|X|_{(1)}\leq|X|_{(2)}\leq\ldots\leq|X|_{(N)}$, i.e., $|X|_{(k)}$ is the
$k$-th smallest entry in $|X_{1}|,|X_{2}|,\ldots,|X_{N}|$.

Based on Lemma 5 and using a chaining argument as in [5, 50], we arrive at the
following lemma, which can be viewed as another (one-sided) variant of the
S-REC in Lemma 4. The proof of Lemma 6 follows easily from [5, 50]; for
completeness, we provide an outline in Appendix D.

######  Lemma 6.

For $\alpha<\mu$ and $\delta>0$, if
$m=\Omega\big{(}\frac{k}{\alpha^{2}}\log\frac{Lr}{\delta}\big{)}$, then with
probability $1-e^{-\Omega(\alpha^{2}m)}$, we have for all
$\mathbf{x}_{1},\mathbf{x}_{2}\in G(B_{2}^{k}(r))$ that

$(\mu-\alpha)\|\mathbf{x}_{1}-\mathbf{x}_{2}\|_{2}-\delta\leq\min_{I\subseteq[m],|I|\geq\frac{m}{2}}\frac{1}{\sqrt{m}}\|\mathbf{A}_{I}(\mathbf{x}_{1}-\mathbf{x}_{2})\|_{2},$
(28)

where for any index set $I\subseteq[m]$, $\mathbf{A}_{I}$ denotes the
$|I|\times n$ sub-matrix of $\mathbf{A}$ that only keeps the rows indexed by
$I$.

###  B-B Proof of Theorem 1

Because $\mathbf{p}\in\mathrm{Range}(G)$, we have

$\displaystyle\|\mathbf{y}-|\mathbf{A}\mathbf{p}|\|_{2}^{2}+\tau$
$\displaystyle\geq\min_{\mathbf{w}\in\mathrm{Range}(G)}\|\mathbf{y}-|\mathbf{A}\mathbf{w}|\|_{2}^{2}+\tau$
(29) $\displaystyle\geq\|\mathbf{y}-|\mathbf{A}\mathbf{q}|\|_{2}^{2}$ (30)
$\displaystyle=\|(\mathbf{y}-|\mathbf{A}\mathbf{p}|)+(|\mathbf{A}\mathbf{p}|-|\mathbf{A}\mathbf{q}|)\|_{2}^{2},$
(31)

where (30) follows from (5). Let $S,T$ be the following index sets:

$\displaystyle
S:=\\{i\in[m]\,:\,\mathrm{sign}(\mathbf{a}_{i}^{T}\mathbf{p})=\mathrm{sign}(\mathbf{a}_{i}^{T}\mathbf{q})\\},$
(32) $\displaystyle
T:=\\{i\in[m]\,:\,\mathrm{sign}(\mathbf{a}_{i}^{T}\mathbf{p})\neq\mathrm{sign}(\mathbf{a}_{i}^{T}\mathbf{q})\\}.$
(33)

Without loss of generality, we assume that $|S|\geq\frac{m}{2}$ (otherwise, we
have $|T|\geq\frac{m}{2}$ and we can derive an upper bound for
$\|\mathbf{q}+\mathbf{x}\|_{2}$ instead of for
$\|\mathbf{q}-\mathbf{x}\|_{2}$). Expanding the squares in (31), we obtain

$\displaystyle\||\mathbf{A}\mathbf{p}|-|\mathbf{A}\mathbf{q}|\|_{2}^{2}$
$\displaystyle\leq
2\langle\mathbf{y}-|\mathbf{A}\mathbf{p}|,|\mathbf{A}\mathbf{q}|-|\mathbf{A}\mathbf{p}|\rangle+\tau$
(34) $\displaystyle\leq
2\|\mathbf{y}-|\mathbf{A}\mathbf{p}|\|_{2}\||\mathbf{A}\mathbf{q}|-|\mathbf{A}\mathbf{p}|\|_{2}+\tau$
(35) $\displaystyle\leq
2\|\mathbf{y}-|\mathbf{A}\mathbf{p}|\|_{2}\|\mathbf{A}(\mathbf{q}-\mathbf{p})\|_{2}+\tau$
(36) $\displaystyle\leq
2(\|\bm{\eta}\|_{2}+\|\mathbf{A}(\mathbf{p}-\mathbf{x})\|_{2})\|\mathbf{A}(\mathbf{q}-\mathbf{p})\|_{2}+\tau,$
(37)

where (36) is from the inequality $||a|-|b||\leq\min\\{|a-b|,|a+b|\\}$ for
$a,b\in\mathbb{R}$, and (37) is from (1) and the triangle inequality. Note
that $\mathbf{p}-\mathbf{x}$ is a fixed vector. From a basic concentration
bound for standard Gaussian matrices (cf., Lemma 1 in Appendix A), we have
with probability $1-e^{-\Omega(m)}$ that

$\|\mathbf{A}(\mathbf{p}-\mathbf{x})\|_{2}\leq
2\sqrt{m}\|\mathbf{p}-\mathbf{x}\|_{2}.$ (38)

In addition, setting $\alpha=\frac{1}{2}$ in Lemma 4, we obtain that if
$m=\Omega\big{(}k\log\frac{Lr}{\delta}\big{)}$, with probability
$1-e^{-\Omega(m)}$, it holds that

$\|\mathbf{A}(\mathbf{q}-\mathbf{p})\|_{2}\leq\sqrt{m}\left(\frac{3}{2}\|\mathbf{q}-\mathbf{p}\|_{2}+\delta\right).$
(39)

Moreover, setting $\alpha=\frac{\mu}{2}$ in Lemma 6, we have that if
$m=\Omega\big{(}k\log\frac{Lr}{\delta}\big{)}$, with probability
$1-e^{-\Omega(m)}$, it holds that

$\displaystyle\||\mathbf{A}\mathbf{p}|-|\mathbf{A}\mathbf{q}|\|_{2}^{2}$
$\displaystyle=\|\mathbf{A}_{S}(\mathbf{p}-\mathbf{q})\|_{2}^{2}+\|\mathbf{A}_{T}(\mathbf{p}+\mathbf{q})\|_{2}^{2}$
(40) $\displaystyle\geq\|\mathbf{A}_{S}(\mathbf{p}-\mathbf{q})\|_{2}^{2}$ (41)
$\displaystyle\geq
m\left(\frac{\mu}{2}\|\mathbf{p}-\mathbf{q}\|_{2}-\delta\right)^{2}$ (42)
$\displaystyle\geq
m\left(\frac{\mu^{2}}{4}\|\mathbf{p}-\mathbf{q}\|_{2}^{2}-\mu\delta\|\mathbf{p}-\mathbf{q}\|_{2}\right),$
(43)

where (40) is from the definitions of $S$ and $T$, and (42) is from Lemma 6
and our assumption $|S|\geq\frac{m}{2}$. Combining (37), (38), (39), and (43),
we obtain

$\frac{\mu^{2}}{4}\|\mathbf{p}-\mathbf{q}\|_{2}^{2}-\mu\delta\|\mathbf{p}-\mathbf{q}\|_{2}\leq
2\left(\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}+2\|\mathbf{p}-\mathbf{x}\|_{2}\right)\left(\frac{3}{2}\|\mathbf{p}-\mathbf{q}\|_{2}+\delta\right)+\frac{\tau}{m}.$
(44)

Simplifying terms in (44), we obtain the desired result.

##  Appendix C Proof of Theorem 2 (Spectral Initializations with Generative
Priors)

Before proving Theorem 2, we provide some intuition and some useful auxiliary
results. Recall also that throughout the analysis, $c$ and $c^{\prime}$
represent small constants whose values may differ from line to line.

###  C-A Ideas Behind the Proof

Since the full proof of Theorem 2 is rather lengthy, we first provide some
intuition behind it. To illustrate our ideas, we consider the simple case that
$\bm{\eta}=\bm{0}$ (i.e, the noiseless case) and $\lambda$ exactly equals its
expectation $\|\mathbf{x}\|_{2}$. Then, we have
$y_{i}=|\langle\mathbf{a}_{i},\mathbf{x}\rangle|$ and $\mathbf{V}$ is reduced
to

$\mathbf{V}_{0}=\frac{1}{m}\sum_{i=1}^{m}y_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{T}\bm{1}_{\\{l<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u\\}},$
(45)

and we obtain (cf. (12))

$\mathbf{J}=\mathbb{E}[\mathbf{V}_{0}]=\|\mathbf{x}\|_{2}(\gamma_{0}\mathbf{I}_{n}+\beta_{0}\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}).$
(46)

In addition, in this simplified setting,
$\hat{\mathbf{x}}\in\mathcal{S}^{n-1}$ is the vector that satisfies

$\hat{\mathbf{x}}=\arg\max_{\mathbf{w}\in\tilde{G}(\mathcal{D})}\mathbf{w}^{T}\mathbf{V}_{0}\mathbf{w}.$
(47)

We aim to derive lower and upper bounds for
$\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$ and then combine these
bounds to derive our desired result.

On the one hand, from $\bar{\mathbf{x}}\in\tilde{G}(\mathcal{D})$, we have

$\displaystyle\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$
$\displaystyle\geq\bar{\mathbf{x}}^{T}\mathbf{V}_{0}\bar{\mathbf{x}}$ (48)
$\displaystyle=\bar{\mathbf{x}}^{T}(\mathbf{V}_{0}-\mathbf{J})\bar{\mathbf{x}}+\bar{\mathbf{x}}^{T}\mathbf{J}\bar{\mathbf{x}}$
(49)
$\displaystyle=\bar{\mathbf{x}}^{T}(\mathbf{V}_{0}-\mathbf{J})\bar{\mathbf{x}}+\|\mathbf{x}\|_{2}(\gamma_{0}+\beta_{0}).$
(50)

With $y_{i}=|\mathbf{a}_{i}^{T}\mathbf{x}|$, we have that
$y_{i}(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\bm{1}_{\\{l<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u\\}}$
is upper bounded by $\|\mathbf{x}\|_{2}u^{3}$. From (50) we can use the
concentration inequality for the sum of sub-Gaussian random variables (cf.
Lemma 2) to derive a lower bound of
$\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$.

On the other hand, we have

$\displaystyle\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$
$\displaystyle=\hat{\mathbf{x}}^{T}(\mathbf{V}_{0}-\mathbf{J})\hat{\mathbf{x}}+\hat{\mathbf{x}}^{T}\mathbf{J}\hat{\mathbf{x}}$
(51)
$\displaystyle=\hat{\mathbf{x}}^{T}(\mathbf{V}_{0}-\mathbf{J})\hat{\mathbf{x}}+\|\mathbf{x}\|_{2}\left(\gamma_{0}+\beta_{0}(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\right).$
(52)

Note that $\hat{\mathbf{x}}$ is dependent on $\mathbf{V}_{0}$. To derive an
upper bound of (52), we construct a $\delta$-net, and write $\hat{\mathbf{x}}$
as

$\hat{\mathbf{x}}=(\hat{\mathbf{x}}-\tilde{\mathbf{x}})+\tilde{\mathbf{x}},$
(53)

where $\tilde{\mathbf{x}}$ is in the $\delta$-net satisfying
$\|\tilde{\mathbf{x}}-\hat{\mathbf{x}}\|_{2}<\delta$. For any
$\mathbf{s}\in\mathbb{R}^{n}$, we have that
$y_{i}(\mathbf{a}_{i}^{T}\mathbf{s})^{2}\bm{1}_{\\{l<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u\\}}$
is sub-exponential with the sub-exponential norm being upper bounded by
$\|\mathbf{x}\|_{2}uC\|\mathbf{s}\|_{2}^{2}$, where $C$ is an absolute
constant. Using a concentration inequality for the sum of sub-exponential
random variables (cf. Lemma 3) and taking a union bound over the $\delta$-net,
we obtain an upper bound for
$\tilde{\mathbf{x}}^{T}(\mathbf{V}_{0}-\mathbf{J})\tilde{\mathbf{x}}$. Using a
high-probability upper bound on $\|\mathbf{a}_{i}\|_{2}$ and the fact that
$\|\tilde{\mathbf{x}}-\hat{\mathbf{x}}\|_{2}<\delta$, we can control the terms
$(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{V}_{0}-\mathbf{J})\hat{\mathbf{x}}$
and
$(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{V}_{0}-\mathbf{J})(\hat{\mathbf{x}}-\tilde{\mathbf{x}})$.
Then, from (52), we obtain an upper bound on
$\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$. Combining the upper and
lower bounds on $\hat{\mathbf{x}}^{T}\mathbf{V}_{0}\hat{\mathbf{x}}$, we can
derive the desired result.

In the full analysis in Section C-C, we additionally carefully deal with the
noise terms and the case that $\lambda$ approximately equals
$\|\mathbf{x}\|_{2}$.

###  C-B Additional Lemmas

Based on Lemma 2, we present the following simple lemma showing that $\lambda$
approximates $\|\mathbf{x}\|_{2}$.

######  Lemma 7.

For any fixed $c\in(0,1)$, with probability $1-e^{-\Omega(m)}$, we have that
$\lambda$ defined in (10) satisfies

$1-c<\frac{\lambda}{\|\mathbf{x}\|_{2}}<1+c.$ (54)

###### Proof.

Since $|\langle\mathbf{a}_{i},\mathbf{x}\rangle|$ is sub-Gaussian with the
sub-Gaussian norm upper bounded by $C\|\mathbf{x}\|_{2}$ and
$\mathbb{E}[|\langle\mathbf{a}_{i},\mathbf{x}\rangle|]=\sqrt{\frac{2}{\pi}}\|\mathbf{x}\|_{2}$,
from Lemma 2, we have with probability $1-e^{-\Omega(m)}$ that

$\left|\frac{1}{m}\sum_{i=1}^{m}|\langle\mathbf{a}_{i},\mathbf{x}\rangle|-\sqrt{\frac{2}{\pi}}\|\mathbf{x}\|_{2}\right|\leq
c\|\mathbf{x}\|_{2}.$ (55)

In addition, by the Cauchy–Schwarz inequality and the upper bound for the
noise term provided in (8), we have

$\left|\frac{1}{m}\sum_{i=1}^{m}\eta_{i}\right|\leq\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c_{0}\|\mathbf{x}\|_{2}.$ (56)

Then, using (55), (56), and a triangle inequality, we obtain

$\left(\sqrt{\frac{2}{\pi}}-c^{\prime}\right)\|\mathbf{x}\|_{2}\leq\frac{1}{m}\sum_{i=1}^{m}y_{i}=\frac{1}{m}\sum_{i=1}^{m}(|\langle\mathbf{a}_{i},\mathbf{x}\rangle|+\eta_{i})\leq\left(\sqrt{\frac{2}{\pi}}+c^{\prime}\right)\|\mathbf{x}\|_{2}.$
(57)

Therefore, we obtain

$(1-c)\|\mathbf{x}\|_{2}<\lambda=\sqrt{\frac{\pi}{2}}\left(\frac{1}{m}\sum_{i=1}^{m}y_{i}\right)<(1+c)\|\mathbf{x}\|_{2}$
(58)

for a choice of $c$ possibly different from that above, but still arbitrarily
small. ∎

For convenience, we present the following lemma concerning the expectation of
the weighted empirical matrix.

######  Lemma 8.

Suppose that $\mathbf{a}\sim\mathcal{N}(\bm{0},\mathbf{I}_{n})$. For any
$\mathbf{s}\in\mathbb{R}^{n}$ and $\eta\in\mathbb{R}$, choosing
$\alpha_{1},\alpha_{2}\in\mathbb{R}$ such that
$1+\frac{\eta}{\|\mathbf{s}\|_{2}}<\alpha_{1}<\alpha_{2}$ and setting
$y=|\langle\mathbf{a},\mathbf{s}\rangle|+\eta$, we have

$\displaystyle\mathbb{E}\left[y\mathbf{a}\mathbf{a}^{T}\mathbf{1}_{\\{\alpha_{1}\|\mathbf{s}\|_{2}<y<\alpha_{2}\|\mathbf{s}\|_{2}\\}}\right]=\|\mathbf{s}\|_{2}\left(\gamma_{1}\mathbf{I}_{n}+\beta_{1}\bar{\mathbf{s}}\bar{\mathbf{s}}^{T}\right)+\eta(\check{\gamma}_{1}\mathbf{I}_{n}+\check{\beta}_{1}\bar{\mathbf{s}}\bar{\mathbf{s}}^{T}),$
(59)

where $\bar{\mathbf{s}}:=\frac{\mathbf{s}}{\|\mathbf{s}\|_{2}}$, and for
$g\sim\mathcal{N}(0,1)$ and
$\phi_{1}(x)=\mathbf{1}_{\\{\alpha_{1}-\frac{\eta}{\|\mathbf{s}\|_{2}}<|x|<\alpha_{2}-\frac{\eta}{\|\mathbf{s}\|_{2}}\\}}$,
we define $\gamma_{1}=\mathbb{E}[|g|\phi_{1}(g)]$,
$\beta_{1}=\mathbb{E}[|g|^{3}\phi_{1}(g)]-\gamma_{1}$,
$\check{\gamma}_{1}=\mathbb{E}[\phi_{1}(g)]$ and
$\check{\beta}_{1}=\mathbb{E}[g^{2}\phi_{1}(g)]-\check{\gamma}_{1}$. (Note
that from the assumption that
$1+\frac{\eta}{\|\mathbf{s}\|_{2}}<\alpha_{1}<\alpha_{2}$, we have
$\beta_{1}>0$ and $\check{\beta}_{1}>0$.)

###### Proof.

Let $g=\langle\mathbf{a},\bar{\mathbf{s}}\rangle\sim\mathcal{N}(0,1)$. We have

$\displaystyle
y\mathbf{a}\mathbf{a}^{T}\mathbf{1}_{\\{\alpha_{1}\|\mathbf{s}\|_{2}<y<\alpha_{2}\|\mathbf{s}\|_{2}\\}}$
$\displaystyle=(\|\mathbf{s}\|_{2}|\mathbf{a}^{T}\bar{\mathbf{s}}|+\eta)\mathbf{a}\mathbf{a}^{T}\bm{1}_{\\{\alpha_{1}-\frac{\eta}{\|\mathbf{s}\|_{2}}<|\langle\mathbf{a},\bar{\mathbf{s}}\rangle|<\alpha_{2}-\frac{\eta}{\|\mathbf{s}\|_{2}}\\}}$
(60)
$\displaystyle=(\|\mathbf{s}\|_{2}|g|+\eta)\phi_{1}(g)\mathbf{a}\mathbf{a}^{T}.$
(61)

For any $i\in[n]$, we have $\mathrm{Cov}[a_{i},g]=\bar{s}_{i}$, and thus
$a_{i}$ can be written as

$a_{i}=\bar{s}_{i}g+t_{i},$ (62)

where $t_{i}\sim\mathcal{N}(0,1-\bar{s}_{i}^{2})$ is independent of $g$. Then,
we have

$\displaystyle\mathbb{E}\left[|g|\phi_{1}(g)a_{i}^{2}\right]$
$\displaystyle=\mathbb{E}\left[|g|\phi_{1}(g)(\bar{s}_{i}g+t_{i})^{2}\right]$
(63)
$\displaystyle=\bar{s}_{i}^{2}\mathbb{E}\left[|g|^{3}\phi_{1}(g)\right]+\left(1-\bar{s}_{i}^{2}\right)\mathbb{E}[|g|\phi_{1}(g)]$
(64) $\displaystyle=\gamma_{1}+\beta_{1}\bar{s}_{i}^{2},$ (65)

and

$\displaystyle\mathbb{E}\left[\phi_{1}(g)a_{i}^{2}\right]$
$\displaystyle=\mathbb{E}\left[\phi_{1}(g)(\bar{s}_{i}g+t_{i})^{2}\right]$
(66)
$\displaystyle=\bar{s}_{i}^{2}\mathbb{E}\left[g^{2}\phi_{1}(g)\right]+\left(1-\bar{s}_{i}^{2}\right)\mathbb{E}[\phi_{1}(g)]$
(67) $\displaystyle=\check{\gamma}_{1}+\check{\beta}_{1}\bar{s}_{i}^{2}.$ (68)

Moreover, for $1\leq i\neq j\leq n$, we have
$0=\mathbb{E}[a_{i}a_{j}]=\mathbb{E}[(\bar{s}_{i}g+t_{i})(\bar{s}_{j}g+t_{j})]=\bar{s}_{i}\bar{s}_{j}+\mathbb{E}[t_{i}t_{j}]$,
which gives $\mathbb{E}[t_{i}t_{j}]=-\bar{s}_{i}\bar{s}_{j}$. Then, similarly
to (65) and (68), we have

$\mathbb{E}\left[|g|\phi_{1}(g)a_{i}a_{j}\right]=\mathbb{E}\left[|g|\phi_{1}(g)(\bar{s}_{i}g+t_{i})(\bar{s}_{j}g+t_{j})\right]=\beta_{1}\bar{s}_{i}\bar{s}_{j},$
(69)

and

$\mathbb{E}\left[\phi_{1}(g)a_{i}a_{j}\right]=\mathbb{E}\left[\phi_{1}(g)(\bar{s}_{i}g+t_{i})(\bar{s}_{j}g+t_{j})\right]=\check{\beta}_{1}\bar{s}_{i}\bar{s}_{j}.$
(70)

Combining (61), (65), (68), (69), (70), we obtain the desired result. ∎

###  C-C Proof of Theorem 2

Let $\mathcal{E}$ be the event that

$1-c\leq\frac{\lambda}{\|\mathbf{x}\|_{2}}\leq 1+c.$ (71)

From Lemma 7, we know that $\mathcal{E}$ occurs with probability
$1-e^{-\Omega(m)}$.

Lower bounding $\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}$: Since
$\bar{\mathbf{x}}\in\tilde{G}(\mathcal{D})$, and $\hat{\mathbf{x}}$ is a
solution of (13), we have under event $\mathcal{E}$ that

$\displaystyle\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}\geq\bar{\mathbf{x}}^{T}\mathbf{V}\bar{\mathbf{x}}$
(72)
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}y_{i}\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l\lambda<y_{i}<u\lambda\\}}$
(73)
$\displaystyle=\frac{\|\mathbf{x}\|_{2}}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|^{3}\mathbf{1}_{\\{l\lambda<y_{i}<u\lambda\\}}+\frac{1}{m}\sum_{i=1}^{m}\eta_{i}\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l\lambda<y_{i}<u\lambda\\}}$
(74)
$\displaystyle\geq\frac{\|\mathbf{x}\|_{2}}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|^{3}\mathbf{1}_{\\{l(1+c^{\prime})\|\mathbf{x}\|_{2}<y_{i}<u(1-c^{\prime})\|\mathbf{x}\|_{2}\\}}$
$\displaystyle\quad-\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l(1-c^{\prime})\|\mathbf{x}\|_{2}<y_{i}<u(1+c^{\prime})\|\mathbf{x}\|_{2}\\}}$
(75)
$\displaystyle=\frac{\|\mathbf{x}\|_{2}}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|^{3}\mathbf{1}_{\\{l(1+c^{\prime})-\frac{\eta_{i}}{\|\mathbf{x}\|_{2}}<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c^{\prime})-\frac{\eta_{i}}{\|\mathbf{x}\|_{2}}\\}}$
$\displaystyle\quad-\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l(1-c^{\prime})-\frac{\eta_{i}}{\|\mathbf{x}\|_{2}}<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c^{\prime})-\frac{\eta_{i}}{\|\mathbf{x}\|_{2}}\\}}$
(76)
$\displaystyle\geq\frac{\|\mathbf{x}\|_{2}}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|^{3}\mathbf{1}_{\\{l(1+c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c)\\}}-\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}},$
(77)

where we use (71) in (75), we use
$y_{i}=|\mathbf{a}_{i}^{T}\mathbf{x}|+\eta_{i}$ in (74) and (76), and we use
(9) as well as $u>l>1+c_{1}$ to derive (77). We aim to bound the two terms in
(77). Let

$\displaystyle\mathbf{U}=\frac{1}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\mathbf{x}|\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l(1+c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c)\\}},$
(78)
$\displaystyle\mathbf{W}=\frac{1}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\mathbf{x}|\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}},$
(79)
$\displaystyle\check{\mathbf{U}}=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l(1+c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c)\\}},$
(80)
$\displaystyle\check{\mathbf{W}}=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\mathbf{a}_{i}\mathbf{a}_{i}^{T}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}.$
(81)

According to the proof of Lemma 8, we have

$\displaystyle\mathbb{E}[\mathbf{U}]=\|\mathbf{x}\|_{2}\left(\gamma\mathbf{I}_{n}+\beta\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}\right),$
(82)
$\displaystyle\mathbb{E}[\mathbf{W}]=\|\mathbf{x}\|_{2}\left(\gamma^{\prime}\mathbf{I}_{n}+\beta^{\prime}\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}\right),$
(83)
$\displaystyle\mathbb{E}[\check{\mathbf{U}}]=\bar{\eta}\left(\check{\gamma}\mathbf{I}_{n}+\check{\beta}\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}\right),$
(84)
$\displaystyle\mathbb{E}[\check{\mathbf{W}}]=\bar{\eta}\left(\check{\gamma}^{\prime}\mathbf{I}_{n}+\check{\beta}^{\prime}\bar{\mathbf{x}}\bar{\mathbf{x}}^{T}\right).$
(85)

where for $g\sim\mathcal{N}(0,1)$ and

$\displaystyle\phi(x):=\bm{1}_{\\{l(1+c)<|x|<u(1-c)\\}},$ (86)
$\displaystyle\psi(x):=\bm{1}_{\\{l(1-c)<|x|<u(1+c)\\}},$ (87)

we set $\gamma=\mathbb{E}[|g|\phi(g)]$,
$\beta=\mathbb{E}[|g|^{3}\phi(g)]-\gamma$,
$\gamma^{\prime}=\mathbb{E}[|g|\psi(g)]$,
$\beta^{\prime}=\mathbb{E}[|g|^{3}\psi(g)]-\gamma^{\prime}$,
$\bar{\eta}=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|$,
$\check{\gamma}=\mathbb{E}[\phi(g)]$,
$\check{\beta}=\mathbb{E}[g^{2}\phi(g)]-\check{\gamma}$,
$\check{\gamma}^{\prime}=\mathbb{E}[\psi(g)]$, and
$\check{\beta}^{\prime}=\mathbb{E}[g^{2}\psi(g)]-\check{\gamma}^{\prime}$.
Then, in (77), we have

$\displaystyle\frac{\|\mathbf{x}\|_{2}}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|^{3}\mathbf{1}_{\\{l(1+c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c)\\}}$
$\displaystyle=\bar{\mathbf{x}}^{T}\mathbf{U}\bar{\mathbf{x}}$ (88)
$\displaystyle=\bar{\mathbf{x}}^{T}(\mathbf{U}-\mathbb{E}[\mathbf{U}])\bar{\mathbf{x}}+\bar{\mathbf{x}}^{T}\mathbb{E}[\mathbf{U}]\bar{\mathbf{x}}$
(89)
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\left(|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\phi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)-\mathbb{E}\left[|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\phi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)\right]\right)+\|\mathbf{x}\|_{2}(\gamma+\beta),$
(90)

where we use (82) and (86) to obtain (90). Note that
$|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\phi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)=|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\bm{1}_{\\{l(1+c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1-c)\\}}\leq
u^{3}(1-c)^{3}\|\mathbf{x}\|_{2}=O(\|\mathbf{x}\|_{2})$. Then, from Lemma 2,
we have with probability $1-e^{-\Omega(m)}$ that

$\left|\frac{1}{m}\sum_{i=1}^{m}\left(|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\phi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)-\mathbb{E}\left[|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\bar{\mathbf{x}})^{2}\phi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)\right]\right)\right|\leq
c\|\mathbf{x}\|_{2}.$ (91)

In addition, we have

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
$\displaystyle\quad=\bar{\mathbf{x}}^{T}\check{\mathbf{W}}\bar{\mathbf{x}}$
(92)
$\displaystyle\quad=\bar{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\bar{\mathbf{x}}+\bar{\mathbf{x}}^{T}\mathbb{E}[\check{\mathbf{W}}]\bar{\mathbf{x}}$
(93)
$\displaystyle\quad=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)-\mathbb{E}\left[\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)\right]\right)+\bar{\eta}(\check{\gamma}^{\prime}+\check{\beta}^{\prime})$
(94)
$\displaystyle\quad\leq\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)-\mathbb{E}\left[\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)\right]\right)+c_{0}\|\mathbf{x}\|_{2}(\check{\gamma}^{\prime}+\check{\beta}^{\prime}),$
(95)

where we use
$\bar{\eta}=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\leq\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c_{0}\|\mathbf{x}\|_{2}$ (see (8)) in (95). From the definition of $\psi$ in
(87), we have
$\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}^{2}\psi\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}\leq
u^{2}(1+c)^{2}=O(1)$. Setting $\alpha_{i}=\frac{|\eta_{i}|}{m}$ and
$X_{i}=\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}^{2}\psi\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}-\mathbb{E}\big{[}\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}^{2}\psi\big{(}\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\big{)}\big{]}$
in Lemma 2, we obtain with probability $1-e^{-\Omega(m)}$ that

$\left|\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)-\mathbb{E}\left[\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)^{2}\psi\left(\mathbf{a}_{i}^{T}\bar{\mathbf{x}}\right)\right]\right)\right|\leq\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c_{0}\|\mathbf{x}\|_{2}.$ (96)

Combining (77), (90), (91), (95) and (96), we have with probability
$1-e^{-\Omega(m)}$ that

$\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}\geq(\gamma+\beta-c)\|\mathbf{x}\|_{2}.$
(97)

Upper bounding $\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}$: We have

$\displaystyle\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}=\frac{1}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\hat{\mathbf{x}})^{2}\bm{1}_{\\{l\lambda<y_{i}<u\lambda\\}}+\frac{1}{m}\sum_{i=1}^{m}\eta_{i}(\mathbf{a}_{i}^{T}\hat{\mathbf{x}})^{2}\bm{1}_{\\{l\lambda<y_{i}<u\lambda\\}}$
(98)
$\displaystyle\leq\frac{1}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\hat{\mathbf{x}})^{2}\bm{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}+\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|(\mathbf{a}_{i}^{T}\hat{\mathbf{x}})^{2}\bm{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
(99)
$\displaystyle=\hat{\mathbf{x}}^{T}\mathbf{W}\hat{\mathbf{x}}+\hat{\mathbf{x}}^{T}\check{\mathbf{W}}\hat{\mathbf{x}}$
(100)
$\displaystyle=\hat{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\hat{\mathbf{x}}+\hat{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\hat{\mathbf{x}}+\|\mathbf{x}\|_{2}\left(\gamma^{\prime}+\beta^{\prime}\left(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}}\right)^{2}\right)+\bar{\eta}(\check{\gamma}^{\prime}+\check{\beta}^{\prime}\left(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\right)$
(101)
$\displaystyle\leq\hat{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\hat{\mathbf{x}}+\hat{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\hat{\mathbf{x}}+\|\mathbf{x}\|_{2}\left(\gamma^{\prime}+\beta^{\prime}\left(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}}\right)^{2}\right)+c\|\mathbf{x}\|_{2},$
(102)

where (99) is derived similarly to (77), (100) comes from the definitions of
$\mathbf{W}$ and $\check{\mathbf{W}}$ in (79) and (81), (101) is from (83) and
(85), and (102) is deduced from $\bar{\eta}\leq c_{0}\|\mathbf{x}\|_{2}$.

We aim to provide an upper bound for (102). For $\delta>0$, let $M$ be a
$(\delta/\tilde{L})$-net of $\mathcal{D}\subseteq B_{2}^{k}(r)$. There exists
such a net with [49, Lemma 5.2]

$\log|M|\leq k\log\frac{4\tilde{L}r}{\delta}.$ (103)

By the $\tilde{L}$-Lipschitz continuity of $\tilde{G}$, we have that
$\tilde{G}(M)$ is a $\delta$-net of $\tilde{G}(\mathcal{D})$. We write

$\hat{\mathbf{x}}=(\hat{\mathbf{x}}-\tilde{\mathbf{x}})+\tilde{\mathbf{x}},$
(104)

where $\tilde{\mathbf{x}}\in\tilde{G}(M)$ with
$\|\hat{\mathbf{x}}-\tilde{\mathbf{x}}\|_{2}\leq\delta$. Then, we have

$\hat{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\hat{\mathbf{x}}=\tilde{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\tilde{\mathbf{x}}+2(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\tilde{\mathbf{x}}+(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])(\hat{\mathbf{x}}-\tilde{\mathbf{x}}).$
(105)

We bound the three terms in (105) separately.

  * •

The first term: For any $\mathbf{s}\in\tilde{G}(M)$, we have from (79) that

$\displaystyle\mathbf{s}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\mathbf{s}=\frac{1}{m}\sum_{i=1}^{m}\left(|\mathbf{a}_{i}^{T}\mathbf{x}|\langle\mathbf{a}_{i},\mathbf{s}\rangle^{2}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}-\mathbf{s}^{T}\mathbb{E}[\mathbf{W}]\mathbf{s}\right).$
(106)

Note that for $\mathbf{s}\in\tilde{G}(M)$, $\|\mathbf{s}\|_{2}=1$. Then,
$|\mathbf{a}_{i}^{T}\mathbf{x}|\langle\mathbf{a}_{i},\mathbf{s}\rangle^{2}\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
is sub-exponential with the sub-exponential norm upper bounded by
$Cu(1+c)\|\mathbf{x}\|_{2}$. From Lemma 3, we have with probability at least
$1-e^{-\Omega(m\epsilon^{2})}$ that

$\left|\mathbf{s}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\mathbf{s}\right|\leq\epsilon\|\mathbf{x}\|_{2}.$
(107)

Taking a union bound over all $\tilde{G}(M)$, we obtain that when
$m=\Omega\big{(}\frac{k}{\epsilon^{2}}\log\frac{\tilde{L}r}{\delta}\big{)}$,
with probability $1-e^{-\Omega(m\epsilon^{2})}$, for all
$\mathbf{s}\in\tilde{G}(M)$,

$\left|\mathbf{s}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\mathbf{s}\right|\leq\epsilon\|\mathbf{x}\|_{2}.$
(108)

Hence, since $\tilde{\mathbf{x}}\in\tilde{G}(M)$, setting $\epsilon$ to be a
sufficiently small absolute constant, we obtain that when
$m=\Omega\big{(}k\log\frac{\tilde{L}r}{\delta}\big{)}$, with probability
$1-e^{-\Omega(m)}$,

$\left|\tilde{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\tilde{\mathbf{x}}\right|\leq
c\|\mathbf{x}\|_{2}.$ (109)

  * •

The second term: From Lemma 2 and a union bound over $[m]$, we have with
probability $1-me^{-\Omega(n)}$ that

$\max_{i\in[m]}\|\mathbf{a}_{i}\|_{2}\leq\sqrt{2n}.$ (110)

Conditioned on the event in (110), we have

$\displaystyle(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\tilde{\mathbf{x}}$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\left(|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\tilde{\mathbf{x}})(\mathbf{a}_{i}^{T}(\hat{\mathbf{x}}-\tilde{\mathbf{x}}))\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}-(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}\mathbb{E}\left[\mathbf{W}\right]\tilde{\mathbf{x}}\right)$
(111)
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}|\mathbf{a}_{i}^{T}\mathbf{x}|(\mathbf{a}_{i}^{T}\tilde{\mathbf{x}})(\mathbf{a}_{i}^{T}(\hat{\mathbf{x}}-\tilde{\mathbf{x}}))\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
$\displaystyle\qquad\qquad\qquad-\|\mathbf{x}\|_{2}\left(\gamma^{\prime}(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}\tilde{\mathbf{x}}+\beta^{\prime}((\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}\bar{\mathbf{x}})(\tilde{\mathbf{x}}^{T}\bar{\mathbf{x}})\right)$
(112)
$\displaystyle\leq(2n(1+c)u+\gamma^{\prime}+\beta^{\prime})\|\mathbf{x}\|_{2}\delta.$
(113)

Therefore, we obtain with probability $1-e^{-\Omega(n)}$ that

$\left|2(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\tilde{\mathbf{x}}\right|\leq
2(2n(1+c)u+\gamma^{\prime}+\beta^{\prime})\|\mathbf{x}\|_{2}\delta.$ (114)

  * •

The third term: By an argument similar to (114), we have with probability
$1-e^{-\Omega(n)}$ that

$\left|(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])(\hat{\mathbf{x}}-\tilde{\mathbf{x}})\right|\leq(2n(1+c)u+\gamma^{\prime}+\beta^{\prime})\|\mathbf{x}\|_{2}\delta^{2}.$
(115)

Combining (105), (109), (114), and (115), and setting $\delta=\frac{c}{n}$, we
obtain that when $m=\Omega(k\log(\tilde{L}nr))$, with probability
$1-e^{-\Omega(m)}$, it holds that

$\left|\hat{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\hat{\mathbf{x}}\right|\leq
c\|\mathbf{x}\|_{2}.$ (116)

Next, we provide an upper bound for
$\hat{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\hat{\mathbf{x}}$.
The strategy is mostly similar to that for upper bounding
$\hat{\mathbf{x}}^{T}(\mathbf{W}-\mathbb{E}[\mathbf{W}])\hat{\mathbf{x}}$ in
(102), but we provide the details below for completeness. Similarly to (105),
we have

$\hat{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\hat{\mathbf{x}}=\tilde{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\tilde{\mathbf{x}}+2(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\tilde{\mathbf{x}}+(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])(\hat{\mathbf{x}}-\tilde{\mathbf{x}}).$
(117)

We bound the three terms in (117) separately.

  * •

The first term: From the definition of $\check{\mathbf{W}}$ in (81), we have
for any $\mathbf{s}\in\tilde{G}(M)$ that

$\displaystyle\mathbf{s}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\mathbf{s}=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|(\mathbf{a}_{i}^{T}\mathbf{s})^{2}\bm{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}-\mathbf{s}^{T}\mathbb{E}[\check{\mathbf{W}}]\mathbf{s}.$
(118)

Note that $\|\mathbf{s}\|_{2}=1$. For any $i\in[m]$,
$(\mathbf{a}_{i}^{T}\mathbf{s})^{2}\bm{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
is sub-exponential with the sub-exponential norm upper bounded by an absolute
constant $C$. Then, from Lemma 3 (or see [58, Theorem 3.1]), we obtain that
for $t>2$, with probability $1-e^{-t}$, it holds that

$\left|\mathbf{s}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\mathbf{s}\right|\leq
C\left(\frac{\|\bm{\eta}\|_{2}}{m}\sqrt{t}+\frac{\|\bm{\eta}\|_{\infty}}{m}t\right).$
(119)

Setting $t=m$ in (119), and from $\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c_{0}\|\mathbf{x}\|_{2}$ (cf. (8)) and $\|\bm{\eta}\|_{\infty}\leq
c_{1}\|\mathbf{x}\|_{2}$ (cf. (9)), we obtain with probability $1-e^{-m}$ that

$\left|\mathbf{s}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\mathbf{s}\right|\leq
c\|\mathbf{x}\|_{2}.$ (120)

Taking a union bound over all $\tilde{G}(M)$, we obtain that when
$m=\Omega\big{(}k\log\frac{\tilde{L}r}{\delta}\big{)}$, with probability
$1-e^{-\Omega(m)}$, for all $\mathbf{s}\in\tilde{G}(M)$,

$\left|\mathbf{s}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\mathbf{s}\right|\leq
c\|\mathbf{x}\|_{2}.$ (121)

Hence, since $\tilde{\mathbf{x}}\in\tilde{G}(M)$, we have

$\left|\tilde{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\tilde{\mathbf{x}}\right|\leq
c\|\mathbf{x}\|_{2}.$ (122)

  * •

The second term: Conditioned on the event in (110), we have

$\displaystyle(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\tilde{\mathbf{x}}$
$\displaystyle\quad=\frac{1}{m}\sum_{i=1}^{m}|\eta_{i}|\left(\mathbf{a}_{i}^{T}(\hat{\mathbf{x}}-\tilde{\mathbf{x}})\right)(\mathbf{a}_{i}^{T}\tilde{\mathbf{x}})\mathbf{1}_{\\{l(1-c)<|\mathbf{a}_{i}^{T}\bar{\mathbf{x}}|<u(1+c)\\}}$
$\displaystyle\qquad\qquad\qquad-\bar{\eta}(\check{\gamma}^{\prime}(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}\tilde{\mathbf{x}}+\check{\beta}^{\prime}((\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}\bar{\mathbf{x}})(\tilde{\mathbf{x}}^{T}\bar{\mathbf{x}}))$
(123)
$\displaystyle\quad\leq\bar{\eta}\delta(2n+\check{\gamma}^{\prime}+\check{\beta}^{\prime}).$
(124)

Therefore, we obtain with probability $1-e^{-\Omega(n)}$ that

$\left|2(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\tilde{\mathbf{x}}\right|\leq
2\bar{\eta}\delta(2n+\check{\gamma}^{\prime}+\check{\beta}^{\prime}).$ (125)

  * •

The third term: By an argument similar to (125), we have with probability
$1-e^{-\Omega(n)}$ that

$\left|(\hat{\mathbf{x}}-\tilde{\mathbf{x}})^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])(\hat{\mathbf{x}}-\tilde{\mathbf{x}})\right|\leq(2n(1+c)u+\gamma^{\prime}+\beta^{\prime})\bar{\eta}\delta^{2}.$
(126)

Note that $\bar{\eta}\leq\frac{\|\bm{\eta}\|_{2}}{\sqrt{m}}\leq
c_{0}\|\mathbf{x}\|_{2}$ (cf. (8)). Combining (117), (122), (125), and (126),
and setting $\delta=\frac{c}{n}$, we obtain that when
$m=\Omega(k\log(\tilde{L}nr))$, with probability $1-e^{-\Omega(m)}$, it holds
that

$\left|\hat{\mathbf{x}}^{T}(\check{\mathbf{W}}-\mathbb{E}[\check{\mathbf{W}}])\hat{\mathbf{x}}\right|\leq
c\|\mathbf{x}\|_{2}.$ (127)

Combining (102), (116), and (127), we obtain that when
$m=\Omega(k\log(\tilde{L}nr))$, with probability $1-e^{-\Omega(m)}$, it holds
that

$\hat{\mathbf{x}}^{T}\mathbf{V}\hat{\mathbf{x}}\leq\left(c+\gamma^{\prime}+\beta^{\prime}(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\right)\|\mathbf{x}\|_{2}.$
(128)

Combining and simplifying: Finally, combining (97) and (128), we obtain that
when $m=\Omega(k\log(\tilde{L}nr))$, with probability $1-e^{-\Omega(m)}$,

$(\gamma+\beta-c)\|\mathbf{x}\|_{2}\leq\left(c+\gamma^{\prime}+\beta^{\prime}(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\right)\|\mathbf{x}\|_{2}.$
(129)

Simplifying (129), we obtain

$\beta^{\prime}(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\geq\beta+(\gamma-\gamma^{\prime})-2c.$
(130)

Then, for $g\sim\mathcal{N}(0,1)$, we have

$\displaystyle\gamma^{\prime}-\gamma$
$\displaystyle=\mathbb{E}[|g|\mathbf{1}_{\\{l(1-c)\leq|g|\leq
l(1+c)\\}}]+\mathbb{E}[|g|\mathbf{1}_{\\{u(1-c)\leq|g|\leq u(1+c)\\}}]$ (131)
$\displaystyle\leq 2\mathbb{E}[|g|\mathbf{1}_{\\{u(1-c)\leq|g|\leq
u(1+c)\\}}]$ (132) $\displaystyle\leq 4u(1+c)\cdot\mathbb{P}(u(1-c)<g<u(1+c))$
(133) $\displaystyle\leq 5u\cdot(\Phi(u(1+c))-\Phi(u(1-c)))$ (134)
$\displaystyle\leq 10C_{0}u^{2}c=O(c),$ (135)

where
$\Phi(x)=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{t^{2}}{2}}\mathrm{d}t$
is the standard normal distribution function and
$C_{0}=\max_{x\in\mathbb{R}}\Phi^{\prime}(x)=\frac{1}{\sqrt{2\pi}}$.
Similarly, we have $\beta\geq\beta^{\prime}-c$. Recall from (71) that the
event $\mathcal{E}$ occurs with probability $1-e^{-\Omega(m)}$. Therefore,
from (130), we obtain that when $m=\Omega(k\log(\tilde{L}nr))$, with
probability $1-e^{-\Omega(m)}$, it holds that

$(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\geq 1-c^{\prime}.$ (136)

Without loss of generality, we assume that
$\hat{\mathbf{x}}^{T}\bar{\mathbf{x}}>0$ (otherwise, we derive an upper bound
for $\|\hat{\mathbf{x}}+\bar{\mathbf{x}}\|_{2}$). We have

$\|\hat{\mathbf{x}}-\bar{\mathbf{x}}\|_{2}^{2}=2(1-\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})\leq
2\left(1-(\hat{\mathbf{x}}^{T}\bar{\mathbf{x}})^{2}\right)\leq 2c^{\prime}.$
(137)

Thus, we obtain (15) as desired. Using Lemma 7, we also obtain (16).

##  Appendix D Proof of Lemma 6 (Variant of the S-REC)

For fixed $\delta>0$ and a positive integer $\ell$, let $M=M_{0}\subseteq
M_{1}\subseteq\ldots\subseteq M_{\ell}$ be a chain of nets of $B_{2}^{k}(r)$
such that $M_{i}$ is a $\frac{\delta_{i}}{L}$-net with
$\delta_{i}=\frac{\delta}{2^{i}}$. There exists such a chain of nets with [49,
Lemma 5.2]

$\log|M_{i}|\leq k\log\frac{4Lr}{\delta_{i}}.$ (138)

By the $L$-Lipschitz assumption on $G$, we have for any $i\in[\ell]$ that
$G(M_{i})$ is a $\delta_{i}$-net of $G(B_{2}^{k}(r))$. We write

$\displaystyle\mathbf{x}_{1}=(\mathbf{x}_{1}-\mathbf{s}_{\ell})+\sum_{i=1}^{\ell}(\mathbf{s}_{i}-\mathbf{s}_{i-1})+\mathbf{s}_{0},\quad\mathbf{x}_{2}=(\mathbf{x}_{2}-\mathbf{t}_{\ell})+\sum_{i=1}^{\ell}(\mathbf{t}_{i}-\mathbf{t}_{i-1})+\mathbf{t}_{0},$
(139)

where $\mathbf{s}_{i},\mathbf{t}_{i}\in G(M_{i})$ for all $i\in[\ell]$, and
$\|\mathbf{x}_{1}-\mathbf{s}_{\ell}\|\leq\frac{\delta}{2^{\ell}}$,
$\|\mathbf{x}_{2}-\mathbf{t}_{\ell}\|\leq\frac{\delta}{2^{\ell}}$,
$\|\mathbf{s}_{i}-\mathbf{s}_{i-1}\|_{2}\leq\frac{\delta}{2^{i-1}}$,
$\|\mathbf{t}_{i}-\mathbf{t}_{i-1}\|_{2}\leq\frac{\delta}{2^{i-1}}$ for all
$i\in[\ell]$. Therefore, the triangle inequality gives

$\|\mathbf{x}_{1}-\mathbf{s}_{0}\|_{2}<2\delta,\quad\|\mathbf{x}_{2}-\mathbf{t}_{0}\|_{2}<2\delta.$
(140)

For any index $I\subseteq[m]$ with $|I|\geq\frac{m}{2}$, from the triangle
inequality, we have

$\displaystyle\|\mathbf{A}_{I}(\mathbf{x}_{1}-\mathbf{x}_{2})\|_{2}$
$\displaystyle\geq\|\mathbf{A}_{I}(\mathbf{s}_{0}-\mathbf{t}_{0})\|_{2}-\sum_{i=1}^{\ell}\left(\|\mathbf{A}_{I}(\mathbf{s}_{i}-\mathbf{s}_{i-1})\|_{2}+\|\mathbf{A}_{I}(\mathbf{t}_{i}-\mathbf{t}_{i-1})\|_{2}\right)$
$\displaystyle\quad-\|\mathbf{A}_{I}(\mathbf{s}_{\ell}-\mathbf{x}_{1})\|_{2}-\|\mathbf{A}_{I}(\mathbf{t}_{\ell}-\mathbf{x}_{2})\|_{2}$
(141)
$\displaystyle\geq\|\mathbf{A}_{I}(\mathbf{s}_{0}-\mathbf{t}_{0})\|_{2}-\sum_{i=1}^{\ell}\left(\|\mathbf{A}(\mathbf{s}_{i}-\mathbf{s}_{i-1})\|_{2}+\|\mathbf{A}(\mathbf{t}_{i}-\mathbf{t}_{i-1})\|_{2}\right)$
$\displaystyle\quad-\|\mathbf{A}(\mathbf{s}_{\ell}-\mathbf{x}_{1})\|_{2}-\|\mathbf{A}(\mathbf{t}_{\ell}-\mathbf{x}_{2})\|_{2}.$
(142)

For any $\mathbf{r}_{1},\mathbf{r}_{2}\in G(M)$ and $0<\alpha<\mu$, from Lemma
5, we have with probability $1-e^{-\Omega(\alpha^{2}m)}$ that

$\|\mathbf{A}_{I}(\mathbf{r}_{1}-\mathbf{r}_{2})\|_{2}\geq\sqrt{m}(\mu-\alpha)\|\mathbf{r}_{1}-\mathbf{r}_{2}\|_{2}.$
(143)

Taking a union bound over $G(M)\times G(M)$, we have that when
$m=\Omega\big{(}\frac{k}{\alpha^{2}}\log\frac{Lr}{\delta}\big{)}$, with
probability $1-e^{-\Omega(\alpha^{2}m)}$, for all
$\mathbf{r}_{1},\mathbf{r}_{2}\in G(M)$,

$\|\mathbf{A}_{I}(\mathbf{r}_{1}-\mathbf{r}_{2})\|_{2}\geq\sqrt{m}(\mu-\alpha)\|\mathbf{r}_{1}-\mathbf{r}_{2}\|_{2}.$
(144)

Therefore, we have that when
$m=\Omega\big{(}\frac{k}{\alpha^{2}}\log\frac{Lr}{\delta}\big{)}$, with
probability $1-e^{-\Omega(\alpha^{2}m)}$,

$\frac{1}{\sqrt{m}}\|\mathbf{A}_{I}(\mathbf{s}_{0}-\mathbf{t}_{0})\|_{2}\geq(\mu-\alpha)\|\mathbf{s}_{0}-\mathbf{t}_{0}\|_{2}\geq(\mu-\alpha)(\|\mathbf{x}_{1}-\mathbf{x}_{2}\|_{2}-4\delta).$
(145)

In addition, using the results in [5, 50], for the standard Gaussian matrix
$\mathbf{A}$, we have
$m=\Omega\big{(}\frac{k}{\alpha^{2}}\log\frac{Lr}{\delta}\big{)}$, with
probability $1-e^{-\Omega(\alpha^{2}m)}$,

$\sum_{i=1}^{\ell}\left(\|\mathbf{A}(\mathbf{s}_{i}-\mathbf{s}_{i-1})\|_{2}+\|\mathbf{A}(\mathbf{t}_{i}-\mathbf{t}_{i-1})\|_{2}\right)\leq
O(\sqrt{m}\alpha\delta)=O(\sqrt{m}\delta).$ (146)

Moreover, by [49, Corollary 5.35], we have
$\big{\|}\frac{1}{\sqrt{m}}\mathbf{A}\big{\|}_{2\to 2}\leq
2+\sqrt{\frac{n}{m}}$ with probability at least $1-e^{-m/2}$. Hence, choosing
$\ell=\lceil\log_{2}n\rceil$, we have with probability at least $1-e^{-m/2}$
that

$\Big{\|}\frac{1}{\sqrt{m}}\mathbf{A}(\mathbf{t}_{\ell}-\mathbf{x}_{2})\Big{\|}_{2}+\Big{\|}\frac{1}{\sqrt{m}}\mathbf{A}(\mathbf{s}_{\ell}-\mathbf{x}_{1})\Big{\|}_{2}\leq
2\left(2+\sqrt{\frac{n}{m}}\right)\frac{\delta}{2^{l}}=O(\delta),$ (147)

where we used the fact that
$\|\mathbf{t}_{\ell}-\mathbf{x}_{2}\|\leq\frac{\delta}{2^{\ell}}$ and
$\|\mathbf{s}_{\ell}-\mathbf{x}_{1}\|\leq\frac{\delta}{2^{\ell}}$. Combining
(142), (145), (146) and (147), we obtain the desired result.

##  Appendix E Further Experiments: Comparing the Iterative Algorithm in
CoPRAM Combined with Various Initialization Methods

In this section, we provide additional numerical results comparing with other
initialization methods, in that case that the initialization is combined with
the subsequent iterative algorithm used in CoPRAM. In addition to PRI-SPCA-NT
and the spectral initialization methods used in ThWF, SPARTA, and CoPRAM, we
also compare with random initialization with the initial vector being
$\lambda\frac{\mathbf{g}}{\|\mathbf{g}\|_{2}}$, where
$\mathbf{g}\sim\mathcal{N}(\bm{0},\mathbf{I}_{n})$. The worst performance of
random initialization reveals the importance of using spectral initialization
methods for sparse phase retrieval.

We first consider the noiseless case and compare the relative error and
empirical success rate of different methods. More specifically, after
obtaining initial vectors from different initialization methods using the
procedure discussed in Section V, we run the iterative algorithm of CoPRAM for
$T=100$ iterations to further refine the estimated vectors and obtain
$\mathbf{x}^{(T)}$. We observe that for our experiments, $100$ iterations are
usually sufficient for convergence. A success is declared for a trial if the
relative error

$\frac{\min\\{\|\mathbf{x}^{(T)}-\mathbf{x}\|_{2},\|\mathbf{x}^{(T)}+\mathbf{x}\|_{2}\\}}{\mathbf{x}}$
(148)

is less than $0.01$. We fix $n=1000$, and consider $s=10$ or $s=20$, while
increasing $m$ from $100$ to $1000$ by $50$. For each of the experiments in
this section, we repeat $50$ trials for $10$ times, and calculate the standard
deviation over these $10$ times to produce the error bars.

We report the relative error in Figure 4, and we report the empirical success
rate of different methods in Figure 5. From Figures 4 and 5, we observe that
for most cases, our PRI-SPCA method leads to the smallest relative error and
the largest empirical success rate.

(a)  (b)

Figure 4: Relative error versus the number of measurements $m$.

(a)  (b)

Figure 5: Empirical success rate versus the number of measurements $m$.

Next, we consider the noisy case, and compare the relative error for using
approximately the same time cost. The time cost is calculated as the sum of
the running time of each initialization method and the running time of the
subsequent iterative algorithm. As we have mentioned in Section V, the time
complexity of each iteration in GRQI is $O(s^{3}+sn)$, while for the iterative
algorithm of CoPRAM, the time complexity of each iteration is $O(s^{2}n\log
n)$, which dominates the total time complexity. In this experiment, we fix
$n=1000$, $m=500$, $s=20$, and set $\sigma=0.1$ or $0.2$. Since the running
times of the initialization methods are typically less than $0.1$ second, we
vary the time cost (in seconds) $t$ in $\\{0.1,0.12,\ldots,0.48,0.5\\}$. For
each $t$ and each of the methods we consider, we find the number of the
iteration whose time cost is closest to $t$, and record the corresponding
relative error. Note that as mentioned in Section V, for noisy measurements,
we do not compare with ThWF because it corresponds to quadratic measurements.
The results are reported in Figure 6, from which we observe that when combined
with the iterative algorithm of CoPRAM and using approximately the same time
cost, our PRI-SPCA method gives smallest relative error for most cases.

(a)  (b)

Figure 6: Relative error for different initialization combined with the
iterative algorithm of CoPRAM when using approximately the same time cost.

## References

  * [1] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: Theory and algorithms,” _IEEE Trans. Inf. Theory_ , vol. 61, no. 4, pp. 1985–2007, 2015.
  * [2] S. Mallat, _A wavelet tour of signal processing_. Elsevier, 1999.
  * [3] S. Foucart and H. Rauhut, _A Mathematical Introduction to Compressive Sensing_. Springer New York, 2013.
  * [4] D. Foster, _Generative Deep Learning : Teaching Machines to Paint, Write, Compose and Play_. O’Reilly Media, Inc, USA, 2019.
  * [5] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” in _Int. Conf. Mach. Learn. (ICML)_ , 2017, pp. 537–546.
  * [6] D. Van Veen, A. Jalal, M. Soltanolkotabi, E. Price, S. Vishwanath, and A. G. Dimakis, “Compressed sensing with deep image prior and learned regularization,” _https://arxiv.org/1806.06438_ , 2018.
  * [7] R. Heckel and P. Hand, “Deep decoder: Concise image representations from untrained non-convolutional networks,” in _Int. Conf. Learn. Repr. (ICLR)_ , 2019.
  * [8] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, and R. Willett, “Deep learning techniques for inverse problems in imaging,” _IEEE J. Sel. Areas Inf. Theory_ , vol. 1, no. 1, pp. 39–56, 2020.
  * [9] J. Whang, Q. Lei, and A. G. Dimakis, “Compressed sensing with invertible generative models and dependent noise,” _https://arxiv.org/2003.08089_ , 2020\.
  * [10] J. Cocola, P. Hand, and V. Voroninski, “Nonasymptotic guarantees for spiked matrix recovery with generative priors,” in _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 33, 2020.
  * [11] A. Jalal, L. Liu, A. G. Dimakis, and C. Caramanis, “Robust compressed sensing using generative models,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , 2020\.
  * [12] R. W. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” _Optik_ , vol. 35, pp. 237–246, 1972\.
  * [13] J. R. Fienup, “Phase retrieval algorithms: A comparison,” _Appl. Opt._ , vol. 21, no. 15, pp. 2758–2769, 1982.
  * [14] E. J. Candès, T. Strohmer, and V. Voroninski, “PhaseLift: Exact and stable signal recovery from magnitude measurements via convex programming,” _Comm. Pure Appl. Math._ , vol. 66, no. 8, pp. 1241–1274, 2013.
  * [15] I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite programming,” _Math. Program._ , vol. 149, no. 1, pp. 47–81, 2015.
  * [16] T. Goldstein and C. Studer, “Phasemax: Convex phase retrieval via basis pursuit,” _IEEE Trans. Inf. Theory_ , vol. 64, no. 4, pp. 2675–2689, 2018\.
  * [17] S. Bahmani and J. Romberg, “Phase retrieval meets statistical learning theory: A flexible convex relaxation,” in _Int. Conf. Artif. Intell. Stat. (AISTATS)_. PMLR, 2017, pp. 252–260.
  * [18] P. Netrapalli, P. Jain, and S. Sanghavi, “Phase retrieval using alternating minimization,” _IEEE Trans. Sig. Proc._ , vol. 63, no. 18, pp. 4814–4826, 2015.
  * [19] Y. Chen and E. J. Candès, “Solving random quadratic systems of equations is nearly as easy as solving linear systems,” _Comm. Pure Appl. Math._ , vol. 70, no. 5, pp. 822–883, 2017.
  * [20] G. Wang, G. B. Giannakis, and Y. C. Eldar, “Solving systems of random quadratic equations via truncated amplitude flow,” _IEEE Trans. Inf. Theory_ , vol. 64, no. 2, pp. 773–794, 2017.
  * [21] H. Zhang, Y. Zhou, Y. Liang, and Y. Chi, “A nonconvex approach for phase retrieval: Reshaped Wirtinger flow and incremental algorithms,” _J. Mach. Learn. Res._ , vol. 18, 2017.
  * [22] J. Sun, Q. Qu, and J. Wright, “A geometric analysis of phase retrieval,” _Found. Comput. Math._ , vol. 18, no. 5, pp. 1131–1198, 2018.
  * [23] Y. Chen, Y. Chi, J. Fan, and C. Ma, “Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval,” _Math. Program._ , vol. 176, no. 1, pp. 5–37, 2019.
  * [24] M. Soltanolkotabi, “Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization,” _IEEE Trans. Inf. Theory_ , vol. 65, no. 4, pp. 2374–2400, 2019.
  * [25] H. Ohlsson, A. Yang, R. Dong, and S. Sastry, “CPRL–An extension of compressive sensing to the phase retrieval problem,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 25, pp. 1367–1375, 2012.
  * [26] Y. Shechtman, A. Beck, and Y. C. Eldar, “GESPAR: Efficient phase retrieval of sparse signals,” _IEEE Trans. Sig. Proc._ , vol. 62, no. 4, pp. 928–938, 2014.
  * [27] P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” _IEEE Trans. Sig. Proc._ , vol. 63, no. 4, pp. 1043–1055, 2014.
  * [28] T. T. Cai, X. Li, Z. Ma _et al._ , “Optimal rates of convergence for noisy sparse phase retrieval via thresholded Wirtinger flow,” _Ann. Stat._ , vol. 44, no. 5, pp. 2221–2251, 2016.
  * [29] G. Wang, L. Zhang, G. B. Giannakis, M. Akçakaya, and J. Chen, “Sparse phase retrieval via truncated amplitude flow,” _IEEE Trans. Sig. Proc._ , vol. 66, no. 2, pp. 479–491, 2017.
  * [30] G. Jagatap and C. Hegde, “Sample-efficient algorithms for recovering structured signals from magnitude-only measurements,” _IEEE Trans. Inf. Theory_ , vol. 65, no. 7, pp. 4434–4456, 2019.
  * [31] J.-F. Cai, J. Li, X. Lu, and J. You, “Sparse signal recovery from phaseless measurements via hard thresholding pursuit,” _https://arxiv.org/abs/2005.08777_ , 2020.
  * [32] J. Scarlett and V. Cevher, “An introductory guide to Fano’s inequality with applications in statistical estimation,” _https://arxiv.org/1901.00555_ , 2019.
  * [33] L. V. Truong and J. Scarlett, “Support recovery in the phase retrieval model: Information-theoretic fundamental limits,” _IEEE Trans. Inf. Theory_ , vol. 66, no. 12, pp. 7887–7910, 2020.
  * [34] F. Wu and P. Rebeschini, “Hadamard Wirtinger flow for sparse phase retrieval,” in _Int. Conf. Artif. Intell. Stat. (AISTATS)_. PMLR, 2021, pp. 982–990.
  * [35] Y. C. Eldar and S. Mendelson, “Phase retrieval: Stability and recovery guarantees,” _Appl. Comput. Harmon. Anal._ , vol. 36, no. 3, pp. 473–494, 2014.
  * [36] M. Iwen, A. Viswanathan, and Y. Wang, “Robust sparse phase retrieval made easy,” _Appl. Comput. Harmon. Anal._ , vol. 1, no. 42, pp. 135–142, 2017\.
  * [37] M. Huang and Z. Xu, “Performance bound of the intensity-based model for noisy phase retrieval,” _https://arxiv.org/abs/2004.08764_ , 2020.
  * [38] ——, “The estimation performance of nonlinear least squares for phase retrieval,” _IEEE Trans. Inf. Theory_ , vol. 66, no. 12, pp. 7967–7977, 2020\.
  * [39] Y. Xia and Z. Xu, “The recovery of complex sparse signals from few phaseless measurements,” _Appl. Comput. Harmon. Anal._ , vol. 50, pp. 1–15, 2021.
  * [40] P. Hand, O. Leong, and V. Voroninski, “Phase retrieval under a generative prior,” in _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , 2018, pp. 9154–9164.
  * [41] R. Hyder, V. Shah, C. Hegde, and M. S. Asif, “Alternating phase projected gradient descent with generative priors for solving compressive phase retrieval,” in _Int. Conf. Acoust. Sp. Sig. Proc. (ICASSP)_. IEEE, 2019, pp. 7705–7709.
  * [42] G. Jagatap and C. Hegde, “Algorithmic guarantees for inverse imaging with untrained network priors,” in _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 32, 2019.
  * [43] F. Shamshad and A. Ahmed, “Compressed sensing-based robust phase retrieval via deep generative priors,” _IEEE Sens. J._ , vol. 21, no. 2, pp. 2286–2298, 2020.
  * [44] X. Wei, Z. Yang, and Z. Wang, “On the statistical rate of nonlinear recovery in generative models with heavy-tailed data,” in _Int. Conf. Mach. Learn. (ICML)_ , 2019, pp. 6697–6706.
  * [45] Y. Plan and R. Vershynin, “The generalized Lasso with non-linear observations,” _IEEE Trans. Inf. Theory_ , vol. 62, no. 3, pp. 1528–1537, 2016.
  * [46] A. Kamath, S. Karmalkar, and E. Price, “Lower bounds for compressed sensing with generative models,” _https://arxiv.org/abs/1912.02938_ , 2019.
  * [47] Z. Liu and J. Scarlett, “Information-theoretic lower bounds for compressive sensing with generative models,” _IEEE J. Sel. Areas Inf. Theory_ , vol. 1, no. 1, pp. 292–303, 2020.
  * [48] H. Zhang, Y. Chi, and Y. Liang, “Provable non-convex phase retrieval with outliers: Median truncated Wirtinger flow,” in _Int. Conf. Mach. Learn. (ICML)_. PMLR, 2016, pp. 1022–1031.
  * [49] R. Vershynin, “Introduction to the non-asymptotic analysis of random matrices,” _https://arxiv.org/abs/1011.3027_ , 2010.
  * [50] Z. Liu, S. Gomes, A. Tiwari, and J. Scarlett, “Sample complexity bounds for 1-bit compressive sensing and binary stable embeddings with generative priors,” in _Int. Conf. Mach. Learn. (ICML)_ , 2020.
  * [51] Z. Liu and J. Scarlett, “The generalized Lasso with nonlinear observations and generative priors,” in _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 33, 2020.
  * [52] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “A simple proof of the restricted isometry property for random matrices,” _Constr. Approx._ , vol. 28, no. 3, pp. 253–263, 2008.
  * [53] H. Zou, T. Hastie, and R. Tibshirani, “Sparse principal component analysis,” _J. Comput. Graph. Stat._ , vol. 15, no. 2, pp. 265–286, 2006.
  * [54] T. Ma and A. Wigderson, “Sum-of-squares lower bounds for sparse PCA,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 28, pp. 1612–1620, 2015.
  * [55] Q. Berthet and P. Rigollet, “Complexity theoretic lower bounds for sparse principal component detection,” in _Conf. Learn. Theory (COLT)_. PMLR, 2013, pp. 1046–1066.
  * [56] T. Wang, Q. Berthet, R. J. Samworth _et al._ , “Statistical and computational trade-offs in estimation of sparse principal components,” _Ann. Stat._ , vol. 44, no. 5, pp. 1896–1930, 2016.
  * [57] C. Gao, Z. Ma, H. H. Zhou _et al._ , “Sparse CCA: Adaptive estimation and computational barriers,” _Ann. Stat._ , vol. 45, no. 5, pp. 2074–2101, 2017.
  * [58] B. Hao, Y. Yadkori, Z. Wen, and G. Cheng, “Bootstrapping upper confidence bound,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , 2019.
  * [59] M. Neykov, Z. Wang, and H. Liu, “Agnostic estimation for misspecified phase retrieval models,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , 2016.
  * [60] Y. S. Tan, “Sparse phase retrieval via sparse PCA despite model misspecification: A simplified and extended analysis,” _https://arxiv.org/1712.04106_ , 2017.
  * [61] B. Moghaddam, Y. Weiss, and S. Avidan, “Spectral bounds for sparse PCA: Exact and greedy algorithms,” _Conf. Neur. Inf. Proc. Sys. (NeurIPS)_ , vol. 18, p. 915, 2006.
  * [62] M. Journée, Y. Nesterov, P. Richtárik, and R. Sepulchre, “Generalized power method for sparse principal component analysis.” _J. Mach. Learn. Res._ , vol. 11, no. 2, 2010.
  * [63] Z. Ma _et al._ , “Sparse principal component analysis and iterative thresholding,” _Ann. Stat._ , vol. 41, no. 2, pp. 772–801, 2013.
  * [64] V. Kuleshov, “Fast algorithms for sparse principal component analysis based on Rayleigh quotient iteration,” in _Int. Conf. Mach. Learn. (ICML)_. PMLR, 2013, pp. 1418–1425.
  * [65] X.-T. Yuan and T. Zhang, “Truncated power method for sparse eigenvalue problems,” _J. Mach. Learn. Res._ , vol. 14, no. 4, 2013.
  * [66] S. S. Vempala, _The random projection method_. American Mathematical Soc., 2005, vol. 65.
  * [67] V. Voroninski and Z. Xu, “A strong restricted isometry property, with an application to phaseless compressed sensing,” _Appl. Comput. Harmon. Anal._ , vol. 40, no. 2, pp. 386–395, 2016.