Solving Systems of Random Quadratic Equations via Truncated Amplitude Flow

I Introduction

Consider a system of $m$ quadratic equations

In many fields of physical sciences and engineering, the problem of recovering the phase from intensity/magnitude-only measurements is commonly referred to as phase retrieval . Relevant application domains include X-ray crystallography , optics , array and high-power coherent diffractive imaging , astronomy , and microscopy . In these settings, due to physical limitations, optical sensors and detectors such as charge-coupled device (CCD) cameras, photosensitive films, and human eyes can record only the (squared) modulus of the Fresnel or Fraunhofer diffraction pattern, while losing the phase of the incident light striking the object. It has been shown that reconstructing a discrete, finite-duration signal from its Fourier transform magnitudes is generally NP-complete . Even checking quadratic feasibility (i.e., whether a solution to a given quadratic system exists or not) is itself an NP-hard problem [18, Theorem 2.6]. Thus, despite its simple form and practical relevance across various fields, tackling the quadratic system in (1) is challenging and NP-hard in general.

Adopting the least-squares criterion, the task of recovering $\bm{x}$ from data $y_{i}$ observed in additive white Gaussian noise (AWGN) can be recast as that of minimizing the intensity-based empirical loss

An alternative is to consider an amplitude-based loss, in which $\psi_{i}$ is observed instead of $y_{i}$ in AWGN

Unfortunately, the presence of quadratic terms in (2) or the modulus in (3) renders the corresponding objective function nonconvex. Minimizing nonconvex objectives, which may exhibit many stationary points, is in general NP-hard . In fact, even checking whether a given point is a local minimum or establishing convergence to a local minimum turns out to be NP-complete .

In the classical discretized one-dimensional ($1$D) phase retrieval, the amplitude vector $\bm{\psi}$ corresponds to the $n$-point Fourier transform of the $n$-dimensional signal $\bm{x}$ . It has been shown based on spectral factorization that in general there is no unique solution to $1$D phase retrieval, even if we disregard trivial ambiguities . To overcome this ill-posedness, several approaches have been suggested. One possibility is to assume additional constraints on the unknown signal such as sparsity . Other approaches rely on introducing redundancy into the measurements using for example, the short-time Fourier transform, or masks . Finally, recent works assume random measurements (e.g., Gaussian $\{\bm{a}_{i}\}$ designs) . Henceforth, this paper focuses on random measurements $\{\psi_{i}\}$ obtained from independently and identically distributed (i.i.d.) Gaussian $\{\bm{a}_{i}\}$ designs.

Existing approaches to solving (2) (or related ones using the Poisson likelihood; see, e.g., ) or (3) fall under two categories: nonconvex and convex ones. Popular nonconvex solvers include alternating projection such as Gerchberg-Saxton and Fineup , AltMinPhase , (Truncated) Wirtinger flow (WF/TWF) , and Karzmarz variants as well as trust-region methods . Inspired by WF, other relevant judiciously initialized counterparts have also been developed for faster semidefinite optimization , blind deconvolution , and matrix completion . Convex counterparts on the other hand rely on the so-called matrix-lifting technique or Shor’s semidefinite relaxation to obtain the solvers abbreviated as PhaseLift , PhaseCut , and CoRK . Further approaches dealing with noisy or sparse phase retrieval are discussed in .

In terms of sample complexity, it has been proven thatThe notation $\phi(n)=\mathcal{O}(g(n))$ means that there is a constant $c>0$ such that $|\phi(n)|\leq c|g(n)|$. $\mathcal{O}(n)$ noise-free random measurements suffice for uniquely determining a general signal . It is also self-evident that recovering a general $n$-dimensional $\bm{x}$ requires at least $\mathcal{O}(n)$ measurements. Convex approaches enable exact recovery from the optimal bound $\mathcal{O}(n)$ of noiseless Gaussian measurements ; they are based on solving a semidefinite program with a matrix variable of size $n\times n$, thus incurring worst-case computational complexity on the order of $\mathcal{O}(n^{4.5})$ that does not scale well with the dimension $n$. Upon exploiting the underlying problem structure, $\mathcal{O}(n^{4.5})$ can be reduced to $\mathcal{O}(n^{3})$ . Solving for vector variables, nonconvex approaches achieve significantly improved computational performance. Using formulation (3) and adopting a spectral initialization commonly employed in matrix completion , AltMinPhase establishes exact recovery with sample complexity $\mathcal{O}(n\log^{3}n)$ under i.i.d. Gaussian $\{\bm{a}_{i}\}$ designs with resampling .

Concerning formulation (2), WF iteratively refines the spectral initial estimate by means of a gradient-like update, which can be approximately interpreted as a stochastic gradient descent variant , . The follow-up TWF improves upon WF through a truncation procedure to separate gradient components of excessively extreme (large or small) sizes. Likewise, due to the heavy tails present in the initialization stage, data $\left\{y_{i}\right\}_{i=1}^{m}$ are pre-screened to yield improved initial estimates in the so-termed truncated spectral initialization method . WF allows exact recovery from $\mathcal{O}(n\log n)$ measurements in $\mathcal{O}(mn^{2}\log(1/\epsilon))$ time/flops to yield an $\epsilon$-accurate solution for any given $\epsilon>0$ , while TWF advances these to $\mathcal{O}(n)$ measurements and $\mathcal{O}(mn\log(1/\epsilon))$ time . Interestingly, the truncation procedure in the gradient stage turns out to be useful in avoiding spurious stationary points in the context of nonconvex optimization, as will be justified in Section IV by the numerical comparison between our amplitude flow (AF) algorithms with or without the judiciously designed truncation rule. It is also worth mentioning that when $m\geq Cn\log^{3}n$ for some sufficiently large positive constant $C$, the objective function in (3) is shown to admit benign geometric structure that allows certain iterative algorithms (e.g., trust-region methods) to efficiently find a global minimizer with random initializations . Hence, the challenge of solving systems of random quadratic equations lies in the case where a near-optimal number of equations are involved, e.g., $m=2n-1$ in the real-valued setting.

Although achieving a linear (in the number of unknowns $n$) sample and computational complexity, the state-of-the-art TWF approach still requires at least $4n\sim 5n$ equations to yield stable empirical success rate (e.g., $\geq 99\%$) under the noiseless real-valued Gaussian model [6, Section 3], which is more than twice the known information-limit of $m=2n-1$ . Similar though less obvious results hold in the complex-valued scenario. While the truncated spectral initialization in improves upon the “plain-vanilla” spectral initialization, its performance still suffers when the number of measurements is relatively small and its advantage (over the untruncated one) diminishes as the number of measurements grows; see more details in Fig. 4 and Section II. Furthermore, extensive numerical and experimental validation confirms that the amplitude-based cost function performs significantly better than the intensity-based one ; that is, formulation (3) is superior to (2). Hence, besides enhancing initialization, markedly improved performance in the gradient stage can be expected by re-examining the amplitude-based cost function and incorporating judiciously designed gradient regularization rules.

I-B This paper

Along the lines of suitably initialized nonconvex schemes and inspired by , the present paper develops a linear-time (i.e., the computational time linearly in both dimensions $m$ and $n$) algorithm to minimize the amplitude-based cost function, referred to as truncated amplitude flow (TAF). Our approach provably recovers an $n$-dimensional unknown signal $\bm{x}$ exactly from a near-optimal number of noiseless random measurements, while also featuring near-perfect statistical performance in the noisy setting. TAF operates in two stages: In the first stage, we introduce an orthogonality-promoting initialization that is computable using a few power iterations. Stage two refines the initial estimate by successive updates of truncated generalized gradient iterations.

Our initialization is built upon the hidden orthogonality characteristics of high-dimensional random vectors , which is in contrast to spectral alternatives originating from the strong law of large numbers (SLLN) . Furthermore, one challenge of phase retrieval lies in reconstructing the signs/phases of $\langle\bm{a}_{i},\bm{x}\rangle$ in the real-/complex-valued settings. Our TAF’s refinement stage leverages a simple yet effective regularization rule to eliminate the erroneously estimated phases in the generalized gradient components with high probability. Simulated tests corroborate that the proposed initialization returns more accurate and robust initial estimates than its spectral counterparts in the noiseless and noisy settings. In addition, our TAF (with gradient truncation) markedly improves upon its “plain-vanilla” version AF. Empirical results demonstrate the advantage of TAF over its competing alternatives.

Focusing on the same amplitude-based cost function, an independent work develops the so-termed reshaped Wirtinger flow (RWF) algorithm , which coincides with amplitude flow (AF). A slightly modified variant of spectral initialization is used to obtain an initial guess, followed by a sequence of non-truncated generalized gradient iterations . Numerical comparisons show that the proposed TAF method performs better than RWF especially when the number of equations approaches the information-theoretic limit ($2n-1$ in the real case).

The remainder of this paper is organized as follows. The amplitude-based cost function, as well as the two algorithmic stages is described and analyzed in Section II. Section III summarizes the TAF algorithm and establishes its theoretical performance. Extensive simulated tests comparing TAF with Wirtinger-based approaches are presented in Section IV. Finally, main proofs are given in Section V, while technical details are deferred to the Appendix.

II Truncated Amplitude Flow

To start, let us define the Euclidean distance of any estimate $\bm{z}$ to the solution set: ${\rm dist}(\bm{z},\,\bm{x}):=\min\,\{\left\|\bm{z}+\bm{x}\right\|,\,\left\|\bm{z}-\bm{x}\right\|\}$ for real signals, and ${\rm dist}(\bm{z},\,\bm{x}):={\rm minimize}_{\phi\in[0,2\pi)}\|\bm{z}-\bm{x}{\rm e}^{i\phi}\|$ for complex ones , where $\|\!\cdot\!\|$ denotes the Euclidean norm. Define also the indistinguishable global phase constant in the real-valued setting as

Henceforth, fixing $\bm{x}$ to be any solution of the given quadratic system (1), we always assume that $\phi\left({\bm{z}}\right)=0$; otherwise, ${\bm{z}}$ is replaced by ${\rm e}^{-j\phi\left({\bm{z}}\right)}{\bm{z}}$, but for simplicity of presentation, the constant phase adaptation term ${\rm e}^{-j\phi\left({\bm{z}}\right)}$ will be dropped whenever it is clear from the context.

For brevity, collect all vectors $\{\bm{a}_{i}\}_{i=1}^{m}$ in the $m\times n$ matrix $\bm{A}:=\left[\bm{a}_{1}~{}\cdots~{}\bm{a}_{m}\right]^{\mathcal{T}}$, and all amplitudes $\left\{\psi_{i}\right\}_{i=1}^{m}$ to form the vector $\bm{\psi}:=\left[\psi_{1}~{}\cdots~{}\psi_{m}\right]^{\mathcal{T}}$. One can rewrite the amplitude-based cost function in matrix-vector representation as

[56, Definition 1.1] The generalized gradient of a function $h$ at $\bm{z}$, denoted by $\partial h$, is the convex hull of the set of limits of the form $\lim\nabla h(\bm{z}_{k})$, where $\bm{z}_{k}\to\bm{z}$ as $k\to+\infty$, i.e.,

Having introduced the notion of a generalized gradient, and with $t$ denoting the iteration count, our approach to solving (5) amounts to iteratively refining the initial guess $\bm{z}_{0}$ (returned by the orthogonality-promoting initialization method to be detailed shortly) by means of the ensuing truncated generalized gradient iterations

for some index set $\mathcal{I}_{t+1}\subseteq[m]:=\left\{1,2,\ldots,m\right\}$ to be designed next. The convention $\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}|}:=0$ is adopted, if $\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}=0$. It is easy to verify that the update in (6) with a full generalized gradient in (7) monotonically decreases the objective function value in (5).

for entry-wise product $\odot$, which may have many solutions. Clearly, if $\bm{z}^{\ast}$ is a solution, then so is $-\bm{z}^{\ast}$. Furthermore, both solutions/global minimizers $\bm{x}$ and $-\bm{x}$ satisfy (8) due to the fact that $\bm{A}\bm{x}-\bm{\psi}\odot\frac{\bm{A}\bm{x}}{\left|\bm{A}\bm{x}\right|}=\bm{0}$. Considering any stationary point $\bm{z}^{\ast}\neq\pm\bm{x}$ that has been adapted such that $\phi(\bm{z}^{\ast})=0$, one can write

Thus, a necessary condition for $\bm{z}^{\ast}\neq\bm{x}$ in (9) is $\frac{\bm{A}\bm{z}^{\ast}}{\left|\bm{A}\bm{z}^{\ast}\right|}\neq\frac{\bm{A}\bm{x}}{\left|\bm{A}\bm{x}\right|}$. Expressed differently, there must be sign differences between $\bm{A}\bm{z}^{\ast}$ and $\bm{A}\bm{x}$ whenever one gets stuck with an undesirable stationary point $\bm{z}^{\ast}$. Inspired by this observation, it is reasonable to devise algorithms that can detect and separate out the generalized gradient components corresponding to mistakenly estimated signs $\left\{\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}|}\right\}$ along the iterates $\{\bm{z}_{t}\}$.

Precisely, if $\bm{z}_{t}$ and $\bm{x}$ lie at different sides of the hyperplane $\bm{a}_{i}^{\mathcal{T}}\bm{z}=0$, then the sign of $\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}$ will be different than that of $\bm{a}_{i}^{\mathcal{T}}\bm{x}$; that is, $\frac{\bm{a}_{i}^{\mathcal{T}}\bm{x}}{|\bm{a}_{i}^{\mathcal{T}}\bm{x}|}\neq\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}|}$. Specifically, one can re-write the $i$-th generalized gradient component as

where $\bm{h}:=\bm{z}-\bm{x}$. Intuitively, the SLLN asserts that averaging the first term $\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}\bm{h}$ over $m$ instances approaches $\bm{h}$, which qualifies it as a desirable search direction. However, certain generalized gradient entries involve erroneously estimated signs of $\bm{a}_{i}^{\mathcal{T}}\bm{x}$; hence, nonzero $\bm{r}_{i}$ terms exert a negative influence on the search direction $\bm{h}$ by dragging the iterate away from $\bm{x}$, and they typically have sizable magnitudes as will be further elaborated in Remark 2 shortly.

Figure 1 demonstrates this from a geometric perspective, where the black dot denotes the origin, and the red dot the solution $\bm{x}$; here, $-\bm{x}$ is omitted for ease of exposition. Assume without loss of generality that the $i$-th missing sign is positive, i.e., $\bm{a}_{i}^{\mathcal{T}}\bm{x}=\psi_{i}$. As will be demonstrated in Theorem 1, with high probability, the initial estimate returned by our orthogonality-promoting method obeys $\|\bm{h}\|\leq\rho\|\bm{x}\|$ for some sufficiently small constant $\rho>0$. Therefore, all points lying on or within the circle (or sphere in high-dimensional spaces) in Fig. 1 satisfy $\|\bm{h}\|\leq\rho\|\bm{x}\|$. If $\bm{a}_{i}^{\mathcal{T}}\bm{z}=0$ does not intersect with the circle, then all points within the circle satisfy $\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}|}=\frac{\bm{a}_{i}^{\mathcal{T}}\bm{x}}{|\bm{a}_{i}^{\mathcal{T}}\bm{x}|}$ qualifying the $i$-th generalized gradient as a desirable search (descent) direction in (10). If, on the other hand, $\bm{a}_{i}^{\mathcal{T}}\bm{z}=0$ intersects the circle, then points lying on the same side of $\bm{a}_{i}^{\mathcal{T}}\bm{z}=0$ with $\bm{x}$ in Fig. 1 admit correctly estimated signs, while points lying on different sides of $\bm{a}_{i}^{\mathcal{T}}\bm{z}=0$ with $\bm{x}$ would have $\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}|}\neq\frac{\bm{a}_{i}^{\mathcal{T}}\bm{x}}{|\bm{a}_{i}^{\mathcal{T}}\bm{x}|}$. This gives rise to a corrupted search direction in (10), implying that the corresponding generalized gradient component should be eliminated.

Nevertheless, it is difficult or even impossible to check whether the sign of $\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}$ equals that of $\bm{a}_{i}^{\mathcal{T}}\bm{x}$. Fortunately, as demonstrated in Fig. 1, most spurious generalized gradient components (those corrupted by nonzero $\bm{r}_{i}$ terms) hover around the watershed hyperplane $\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}=0$. For this reason, TAF includes only those components having $\bm{z}_{t}$ sufficiently away from its watershed, i.e.,

Regarding our gradient regularization rule in (11), two observations are in order.

Our truncation rule deviates from the intuition behind TWF, which throws away gradient components corresponding to large-size $\{|\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}|/|\bm{a}_{i}^{\mathcal{T}}\bm{x}|\}$ in (11). As demonstrated by our analysis in Appendix A-E, it rarely happens that a gradient component having large $|\bm{a}_{i}^{\mathcal{T}}\bm{z}_{t}|/|\bm{a}_{i}^{\mathcal{T}}\bm{x}|$ yields an incorrect sign of $\bm{a}_{i}^{\mathcal{T}}\bm{x}$ under a sufficiently accurate initialization. Moreover, discarding too many samples (those for which $i\notin\mathcal{T}_{t+1}$ in TWF [6, Section 2.1]) introduces large bias into $(1/m)\sum_{i\in\mathcal{T}_{t+1}}^{m}\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}\bm{h}$, so that TWF does not work well when $m/n$ is close to the information-limit of $m/n\approx 2$. In sharp contrast, the motivation and objective of our truncation rule in (11) is to directly sense and eliminate gradient components that involve mistakenly estimated signs with high probability.

To demonstrate the power of TAF, numerical tests comparing all stages of (T)AF and (T)WF will be presented throughout our analysis. The basic test settings used in this paper are described next. For fairness, all pertinent algorithmic parameters involved in all compared schemes were set to their default values. Simulated estimates are averaged over $100$ independent Monte Carlo (MC) realizations without mentioning this explicitly each time. Performance of different schemes is evaluated in terms of the relative root mean-square error, i.e.,

and the success rate among $100$ trials, where a success is claimed for a trial if the returned estimate incurs a relative error less than $10^{-5}$ . Simulated tests under both noiseless and noisy Gaussian models are performed, corresponding to \psi_{i}=\big{|}\bm{a}_{i}^{\mathcal{H}}\bm{x}+\eta_{i}\big{|} with $\eta_{i}=0$ and $\eta_{i}\sim\mathcal{N}(0,\sigma^{2})$, respectively, with i.i.d. $\bm{a}_{i}\sim\mathcal{N}(\bm{0},\bm{I}_{n})$ or $\bm{a}_{i}\sim\mathcal{CN}(\bm{0},\bm{I}_{n})$.

Numerical comparison depicted in Fig. 2 using the noiseless real-valued Gaussian model suggests that even when starting with the same truncated spectral initialization, TAF’s refinement outperforms those of TWF and WF, demonstrating the merits of our gradient update rule over TWF/WF. Furthermore, comparing TAF (gradient iterations in (6)-(7) with truncation in (11) initialized by the truncated spectral estimate) and AF (gradient iterations in (6)-(7) initialized by the truncated spectral estimate) corroborates the power of the truncation rule in (11).

II-B Orthogonality-promoting initialization stage

Leveraging the SLLN, spectral initialization methods estimate $\bm{x}$ as the (appropriately scaled) leading eigenvector of $\bm{Y}:=\frac{1}{m}\sum_{i\in\mathcal{T}_{0}}y_{i}\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}$, where $\mathcal{T}_{0}$ is an index set accounting for possible data truncation. As asserted in , each summand $(\bm{a}_{i}^{\mathcal{T}}\bm{x})^{2}\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}$ follows a heavy-tail probability density function lacking a moment generating function. This causes major performance degradation especially when the number of measurements is small. Instead of spectral initializations, we shall take another route to bypass this hurdle. To gain intuition into our initialization, a motivating example is presented first that reveals fundamental characteristics of high-dimensional random vectors.

where $\theta_{i}$ is the angle between vectors $\bm{a}_{i}$ and $\bm{x}$. Consider ordering all $\{\cos^{2}\theta_{i}\}$ in an ascending fashion, and collectively denote them as $\bm{\xi}:=[\cos^{2}\theta_{[m]}~{}\cdots~{}\cos^{2}\theta_{}]^{\mathcal{T}}$ with $\cos^{2}\theta_{}\geq\cdots\geq\cos^{2}\theta_{[m]}$. Figure 3 plots the ordered entries in $\bm{\xi}$ for $m/n$ varying by $2$ from $2$ to $10$ with $n=1,000$. Observe that almost all $\left\{\bm{a}_{i}\right\}$ vectors have a squared normalized inner-product with $\bm{x}$ smaller than $10^{-2}$, while half of the inner-products are less than $10^{-3}$, which implies that $\bm{x}$ is nearly orthogonal to a large number of $\bm{a}_{i}$’s.

This example corroborates the folklore that random vectors in high-dimensional spaces are almost always nearly orthogonal to each other . This inspired us to pursue an orthogonality-promoting initialization method. Our key idea is to approximate $\bm{x}$ by a vector that is most orthogonal to a subset of vectors $\{\bm{a}_{i}\}_{i\in\mathcal{I}_{0}}$, where $\mathcal{I}_{0}$ is an index set with cardinality $|\mathcal{I}_{0}|<m$ that includes indices of the smallest squared normalized inner-products $\left\{\cos^{2}\theta_{i}\right\}$. Since $\left\|\bm{x}\right\|$ appears in all inner-products, its exact value does not influence their ordering. Henceforth, we assume with no loss of generality that $\|\bm{x}\|=1$.

Using data $\left\{(\bm{a}_{i};\,\psi_{i})\right\}$, evaluate $\cos^{2}\theta_{i}$ according to (14) for each pair $\bm{x}$ and $\bm{a}_{i}$. Instrumental for the ensuing derivations is noticing from the inherent near-orthogonal property of high-dimensional random vectors that the summation of $\cos^{2}\theta_{i}$ over all indices $i\in\mathcal{I}_{0}$ should be very small; rigorous justification is deferred to Section V. Therefore, the sum $\sum_{i\in\mathcal{I}_{0}}\cos^{2}\theta_{i}$ is also small, or according to (14), equivalently,

is small. Therefore, a meaningful approximation of $\bm{x}$ can be obtained by minimizing the former with $\bm{x}$ replaced by the optimization variable $\bm{z}$, namely

This amounts to finding the smallest eigenvalue and the associated eigenvector of $\bm{Y}_{0}:=\frac{1}{|\mathcal{I}_{0}|}\sum_{i\in\mathcal{I}_{0}}\frac{\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}}{\|\bm{a}_{i}\|^{2}}\succeq\bm{0}$ (the symbol $\succeq$ means positive semidefinite). Finding the smallest eigenvalue calls for eigen-decomposition or matrix inversion, each typically requiring computational complexity on the order of $\mathcal{O}(n^{3})$. Such a computational burden may be intractable when $n$ grows large. Applying a standard concentration result, we show how the computation can be significantly reduced.

which can be efficiently solved via simple power iterations.

It is clear from (21) that the first term on the right hand side of (22) approximates $\|\bm{x}\|^{2}$. The second term approaches $1$ because the denominator $(1/m)\cdot\sum_{i=1}^{m}\|\bm{a}_{i}\|^{2}\approx n$ appealing to the SLLN again and the fact that $\bm{a}_{i}\sim\mathcal{N}(\bm{0},\bm{I}_{n})$. For simplicity, we choose to work with the first norm estimate

It is worth highlighting that, compared to the matrix $\bm{Y}:=\frac{1}{m}\sum_{i\in\mathcal{T}_{0}}y_{i}\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}$ used in spectral methods, our constructed matrix $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{Y}}_{0}$ in (18) does not depend on the observed data $\{y_{i}\}$ explicitly; the dependence is only through the choice of the index set $\mathcal{I}_{0}$. The novel orthogonality-promoting initialization thus enjoys two advantages over its spectral alternatives: a1) it does not suffer from heavy-tails of the fourth-order moments of Gaussian $\{\bm{a}_{i}\}$ vectors common in spectral initialization schemes; and, a2) it is less sensitive to noisy data.

III Main Results

The TAF algorithm is summarized in Algorithm 1. Default values are set for pertinent algorithmic parameters. Assuming independent data samples $\{(\bm{a}_{i};\psi_{i})\}$ drawn from the noiseless real-valued Gaussian model, the following result establishes the theoretical performance of TAF.

with $\rho={1}/{10}$ (or any sufficiently small positive constant), provided that $m\geq c_{1}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\geq c_{2}n$ for some numerical constants $c_{1},\,c_{2}>0$, and sufficiently large $n$. Furthermore, choosing a constant step size $\mu\leq\mu_{0}$ along with a truncation level $\gamma\geq 1/2$, and starting from any initial guess $\bm{z}_{0}$ satisfying (24), successive estimates of the TAF solver (tabulated in Algorithm 1) obey

for some $0<\nu<1$, which holds with probability exceeding $1-(m+5){\rm e}^{-n/2}-8{\rm e}^{-c_{0}m}-1/n^{2}$.

Typical parameter values for TAF in Algorithm 1 are $\mu=0.6$, and $\gamma=0.7$. The proof of Theorem 1 is relegated to Section V. Theorem 1 asserts that: i) TAF reconstructs the solution $\bm{x}$ exactly as soon as the number of equations is about the number of unknowns, which is theoretically order optimal. Our numerical tests demonstrate that for the real-valued Gaussian model, TAF achieves a success rate of $100\%$ when $m/n$ is as small as $3$, which is slightly larger than the information limit of $m/n=2$ (Recall that $m\geq 2n-1$ is necessary for the uniqueness.) This is a significant reduction in the sample complexity ratio, which is $5$ for TWF and $7$ for WF. Surprisingly, TAF also enjoys a success rate of over $50\%$ when $m/n$ is the information limit $2$, which has not yet been presented for any existing algorithms. See further discussion in Section IV; and, ii) TAF converges exponentially fast with convergence rate independent of the dimension $n$. Specifically, TAF requires at most $\mathcal{O}(\log(1/\epsilon))$ iterations to achieve any given solution accuracy $\epsilon>0$ (a.k.a., ${\rm dist}(\bm{z}_{t},\bm{x})\leq\epsilon\left\|\bm{x}\right\|$), with iteration cost $\mathcal{O}(mn)$. Since the truncation takes time on the order of $\mathcal{O}(m)$, the computational burden of TAF per iteration is dominated by the evaluation of the gradient components. The latter involves two matrix-vector multiplications that are computable in $\mathcal{O}(mn)$ flops, namely, $\bm{A}\bm{z}_{t}$ yields $\bm{u}_{t}$, and $\bm{A}^{\mathcal{T}}\bm{v}_{t}$ the gradient, where $\bm{v}_{t}:=\bm{u}_{t}-\bm{\psi}\odot\tfrac{\bm{u}_{t}}{|\bm{u}_{t}|}$. Hence, the total running time of TAF is $\mathcal{O}(mn\log(1/\epsilon))$, which is proportional to the time taken to read the data $\mathcal{O}(mn)$.

In the noisy setting, TAF is stable under additive noise. To be more specific, consider the amplitude-based data model $\psi_{i}=|\bm{a}_{i}^{\mathcal{T}}\bm{x}|+\eta_{i}$. It can be shown that the truncated amplitude flow estimates in Algorithm 1 satisfy

IV Simulated Tests

In this section, we provide additional numerical tests evaluating performance of the proposed scheme relative to (T)WF Matlab codes directly downloaded from the authors’ websites: http://statweb.stanford.edu/~candes/TWF/algorithm.html; http://www-bcf.usc.edu/~soltanol/WFcode.html. and AF. The initial estimate was found based on $50$ power iterations, and was subsequently refined by $T=1,000$ gradient-type iterations in each scheme. The Matlab implementations of TAF are available at https://gangumn.github.io/TAF/ for reproducibility.

Left panel in Fig. 5 presents the average relative error of three initialization methods on a series of noiseless/noisy real-valued Gaussian problems with $m/n=6$ fixed, and $n$ varying from $500$ to $10^{4}$, while those for the corresponding complex-valued Gaussian instances are shown in the right panel. Clearly, the proposed initialization method returns more accurate and robust estimates than the spectral ones. Under the same condition for the real-valued Gaussian model, Fig. 6 compares the initialization implemented in Algorithm 1 obtained by solving the maximum eigenvalue problem in (19) with the one obtained by tackling the minimum eigenvalue problem in (16) via the Lanczos method . When the number of equations is relatively small (less than about $3n$), the former performs better than the latter. Interestingly though, the latter works remarkably well and almost halves the error incurred by the implemented initialization of Algorithm 1 as soon as the number of equations becomes larger than $4$.

To demonstrate the power of TAF, Fig. 7 plots the relative error of recovering a real-valued signal in logarithmic scale versus the iteration count under the information-limit of $m=2n-1$ noiseless i.i.d. Gaussian measurements . In this case, since the returned initial estimate is relatively far from the optimal solution (see Fig. 4), TAF converges slowly for the first $200$ iterations or so due to elimination of a significant amount of ‘bad’ generalized gradient components (corrupted by mistakenly estimated signs). As the iterate gets more accurate and lands within a small-size neighborhood of $\bm{x}$, TAF converges exponentially fast to the globally optimal solution. It is worth emphasizing that no existing method succeeds in this case. Figure 8 compares the empirical success rate of three schemes under both real-valued and complex-valued Gaussian models with $n=10^{3}$ and $m/n$ varying by $0.1$ from $1$ to $7$, where a success is claimed if the estimate has a relative error less than $10^{-5}$. For real-valued vectors, TAF achieves a success rate of over $50\%$ when $m/n=2$, and guarantees perfect recovery from about $3n$ measurements; while for complex-valued ones, TAF enjoys a success rate of $95\%$ when $m/n=3.4$, and ensures perfect recovery from about $4.5n$ measurements.

To demonstrate the stability of TAF, the relative mean-squared error (MSE)

as a function of the signal-to-noise ratio (SNR) is plotted for different $m/n$ values. We consider the noisy model $\psi_{i}=|\langle\bm{a}_{i},\bm{x}\rangle|+\eta_{i}$ with $\bm{x}\sim\mathcal{N}(\bm{0},\bm{I}_{1,000})$ and real-valued independent Gaussian sensing vectors $\bm{a}_{i}\sim\mathcal{N}(\bm{0},\bm{I}_{1,000})$, in which $m/n$ takes values $\{6,\,8,\,10\}$, and the SNR in dB, given by

is varied from $10$ dB to $50$ dB. Averaging over $100$ independent trials, Fig. 9 demonstrates that the relative MSE for all $m/n$ values scales inversely proportional to SNR, hence justifying the stability of TAF under bounded additive noise.

The next experiment evaluates the efficacy of the proposed initialization method, simulating all schemes initialized by the truncated spectral initial estimate and the orthogonality-promoting initial estimate. Apparently, all algorithms except WF admit a significant performance improvement when initialized by the proposed orthogonality-promoting initialization relative to the truncated spectral initialization. Nevertheless, TAF with our developed orthogonality-promoting initialization enjoys superior performance over all simulated approaches.

where $\bm{F}$ denotes the $n\times n$ discrete Fourier transform matrix, and $\bm{D}^{(k)}$ is a diagonal matrix holding entries sampled uniformly at random from $\{1,\,-1,\,j,\,-j\}$ (phase delays) on its diagonal, with $j$ denoting the imaginary unit. Each $\bm{D}^{(k)}$ represents a random mask placed after the object . With $K=6$ masks implemented in our experiment, the total number of quadratic measurements is $m=6n$. Every algorithm was run independently on each of the three bands. A number $100$ of power iterations were used to obtain an initialization, which was refined by $100$ gradient-type iterations. The relative errors after our orthogonality-promoting initialization and after $100$ TAF iterations are $0.6807$ and $9.8631\times 10^{-5}$, respectively, and the recovered images are displayed in Fig. 11. In sharp contrast, TWF returns images of corresponding relative errors $1.3801$ and $1.3409$, which are far away from the ground truth.

Regarding running times in all performed experiments, TAF converges slightly faster than TWF, while both are markedly faster than WF. All experiments were implemented using MATLAB on an Intel CPU @ $3.4$ GHz ($32$ GB RAM) computer.

V Proofs

This section presents the main ideas behind the proof of Theorem 1, and establishes a few necessary lemmas. Technical details are deferred to the Appendix. Relative to WF and TWF, our objective function involves nonsmoothness and nonconvexity, rendering the proof of exact recovery of TAF nontrivial. In addition, our initialization method starts from a rather different perspective than spectral alternatives, so that the tools involved in proving performance of our initialization deviate from those of spectral methods . Part of our proof is adapted from and .

The proof of Theorem 1 consists of two parts: Section V-A justifies the performance of the proposed orthogonality-promoting initialization, which essentially achieves any given constant relative error as soon as the number of equations is on the order of the number of unknowns, namely, $m\asymp n$.The notations $\phi(n)=\mathcal{O}(g(n))$ or $\phi(n)\gtrsim g(n)$ (respectively, $\phi(n)\lesssim g(n)$) means there exists a numerical constant $c>0$ such that $\phi(n)\leq cg(n)$, while $\phi(n)\asymp g(n)$ means $\phi(n)$ and $g(n)$ are orderwise equivalent. Section V-B demonstrates theoretical convergence of TAF to the solution of the quadratic system in (1) at a geometric rate provided that the initial estimate has a sufficiently small constant relative error as in (24). The two stages of TAF can be performed independently, meaning that better initialization methods, if available, could be adopted to initialize our truncated generalized gradient iterations; likewise, our initialization may be applied to other iterative optimization algorithms.

This section concentrates on proving guaranteed performance of the proposed orthogonality-promoting initialization method, as asserted in the following proposition. An alternative approach may be found in .

for $\rho=1/10$ or any positive constant, with the proviso that $m\geq c_{1}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\geq c_{2}n$ for some numerical constants $c_{1},\,c_{2}>0$ and sufficiently large $n$.

Due to homogeneity in (28), it suffices to consider the case $\|\bm{x}\|=1$. Assume for the moment that $\left\|\bm{x}\right\|=1$ is known and $\bm{z}_{0}$ has been scaled such that $\left\|\bm{z}_{0}\right\|=1$ in (23). The error between the employed $\bm{x}$’s norm estimate $\sqrt{\frac{1}{m}\sum_{i=1}^{m}y_{i}}$ and the unknown norm $\left\|\bm{x}\right\|=1$ will be accounted for at the end of this section. Instrumental in proving Proposition 1 is the following result, whose proof is provided in Appendix A-A.

We now turn to prove Proposition 1. The first step consists in upper-bounding the term on the right-hand-side of (29). Specifically, its numerator is upper bounded, and the denominator lower bounded, as summarized in Lemma 2 and Lemma 3 next; their proofs are provided in Appendix A-B and Appendix A-C, respectively.

In the setup of Lemma 1, if $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\geq c_{1}^{\prime}n$, then

holds with probability at least $1-2{\rm e}^{-c_{K}n}$, where $c_{2}^{\prime}$ and $c_{K}$ are some universal constants.

In the setup of Lemma 1, the following holds with probability at least $1-(m+1){\rm e}^{-n/2}-{\rm e}^{-c_{0}m}-1/n^{2}$,

provided that $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\geq c_{1}^{\prime}n$, $m\geq c_{2}^{\prime}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$, and $m\geq c_{3}^{\prime}n$ for some absolute constants $c_{1}^{\prime},\,c_{2}^{\prime},\,c_{3}^{\prime}>0$, and sufficiently large $n$.

Leveraging the upper and lower bounds in (30) and (31), one arrives at

which holds with probability at least $1-(m+3){\rm e}^{-n/2}-{\rm e}^{-c_{0}m}-1/n^{2}$, assuming that $m\geq c_{1}^{\prime}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$, and $m\geq c_{2}^{\prime}n$, $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\geq c_{3}^{\prime}n$ for some absolute constants $c_{1}^{\prime},\,c_{2}^{\prime},\,c_{3}^{\prime}>0$, and sufficiently large $n$.

Summarizing the two inequalities, we conclude that

V-B Exact recovery from noiseless data

We now prove that with accurate enough initial estimates, TAF converges at a geometric rate to $\bm{x}$ with high probability (i.e., the second part of Theorem 1). To be specific, with initialization obeying (28) in Proposition 1, TAF reconstructs the solution $\bm{x}$ exactly in linear time. To start, it suffices to demonstrate that the TAF’s update rule (i.e., Step 4 in Algorithm 1) is locally contractive within a sufficiently small neighborhood of $\bm{x}$, as asserted in the following proposition.

Consider the noise-free measurements $\psi_{i}=\left|\bm{a}_{i}^{\mathcal{T}}\bm{x}\right|$ with i.i.d. Gaussian design vectors $\bm{a}_{i}\sim\mathcal{N}(\bm{0},\,\bm{I}_{n})$, $1\leq i\leq m$, and fix any $\gamma\geq 1/2$. There exist universal constants $c_{0},\,c_{1}>0$ and $0<\nu<1$ such that with probability at least $1-7{\rm e}^{-c_{0}m}$, the following holds

Proposition 2 demonstrates that the distance of TAF’s successive iterates to $\bm{x}$ is monotonically decreasing once the algorithm enters a small-size neighborhood around $\bm{x}$. This neighborhood is commonly referred to as the basin of attraction; see further discussions in . In other words, as soon as one lands within the basin of attraction, TAF’s iterates remain in this region and will be attracted to $\bm{x}$ exponentially fast. To substantiate Proposition 2, recall the local regularity condition, which was first developed in and plays a fundamental role in establishing linear convergence to global optimum of nonconvex optimization approaches such as WF/TWF .

for all $\bm{z}$ obeying $\left\|\bm{h}\right\|\leq\epsilon\left\|\bm{x}\right\|$. Evidently, if the ${\rm LRC}(\mu,\lambda,\epsilon)$ is proved for TAF, then (37) follows upon letting $\nu:=\lambda\mu$.

On the other hand, standard matrix concentration results confirm that the largest singular value of $\bm{A}=\left[\bm{a}_{1}~{}\cdots~{}\bm{a}_{m}\right]^{\mathcal{T}}$ with i.i.d. Gaussian $\{\bm{a}_{i}\}$ satisfies $\sigma_{1}:=\|\bm{A}\|\leq(1+\delta^{\prime\prime})\sqrt{m}$ for some $\delta^{\prime\prime}>0$ with probability exceeding $1-2{\rm e}^{-c_{0}m}$ as soon as $m\geq c_{1}n$ for sufficiently large $c_{1}>0$, where $c_{1}>0$ is a universal constant depending on $\delta^{\prime\prime}$ [59, Remark 5.25]. Combining (42), (43), and (44) yields

which occupies the remaining of this section. Formally, this can be stated as follows.

Consider the noiseless measurements $\psi_{i}=|\bm{a}_{i}^{\mathcal{T}}\bm{x}|$, and fix any sufficiently small constant $\epsilon>0$. There exist universal constants $c_{0},\,c_{1}>0$ such that if $m>c_{1}n$, then the following holds with probability exceeding $1-4{\rm e}^{-c_{0}m}$:

Before justifying Proposition 3, we introduce the following events.

Fix any $\gamma>0$. For each $i\in[m]$, define

From Fig. 1, it is clear that if $\bm{z}\in\xi_{i}^{2}$, then the sign of $\bm{a}_{i}^{\mathcal{T}}\bm{z}$ will be different than that of $\bm{a}_{i}^{\mathcal{T}}\bm{x}$. The region $\xi_{i}^{2}$ can be readily specified by the conditions that

Under our initialization condition $\left\|\bm{h}\right\|/\left\|\bm{x}\right\|\leq\rho$, it is self-evident that $\mathcal{D}_{i}$ describes two symmetric spherical caps over $\bm{a}_{i}^{\mathcal{T}}\bm{x}=\psi_{i}$ with one being $\xi_{i}^{2}$. Hence, it holds that $\mathcal{E}_{i}\cap\mathcal{K}_{i}=\xi_{i}^{2}\subseteq\mathcal{D}_{i}\cap\mathcal{K}_{i}$. ∎

To prove (47), consider rewriting the truncated gradient in terms of the events defined in Lemma 4:

Using the definitions and properties in Lemma 4, one further arrives at

where the last inequality arises from the property ${\left|\bm{a}_{i}^{\mathcal{T}}\bm{x}\right|}\leq\frac{1+\gamma}{2+\gamma}{\left|\bm{a}_{i}^{\mathcal{T}}\bm{h}\right|}$ by the definition of $\mathcal{D}_{i}$.

Establishing the regularity condition or Proposition 3, boils down to lower bounding the right-hand side of (V-B1), namely, to lower bounding the first term and to upper bounding the second one. By the SLLN, the first term in (V-B1) approximately gives $\left\|\bm{h}\right\|^{2}$ as long as our truncation procedure does not eliminate too many generalized gradient components (i.e., summands in the first term). Regarding the second, one would expect its contribution to be small under our initialization condition in (28) and as the relative error $\left\|\bm{h}\right\|/\left\|\bm{x}\right\|$ decreases. Specifically, under our initialization, $\mathcal{D}_{i}$ is provably a rare event, thus eliminating the possibility of the second term exerting a noticeable influence on the first term. Rigorous analyses concerning the two terms are elaborated in Lemma 5 and Lemma 6, whose proofs are provided in Appendix A-D and Appendix A-E, respectively.

Fix $\gamma\geq 1/2$ and $\rho\leq 1/10$, and let $\mathcal{E}_{i}$ be defined in (48). For independent random variables $W\sim\mathcal{N}(0,\,1)$ and $Z\sim\mathcal{N}(0,\,1)$, set

Then for any $\epsilon>0$ and any vector $\bm{h}$ obeying $\left\|\bm{h}\right\|/\left\|\bm{x}\right\|\leq\rho$, the following holds with probability exceeding $1-2{\rm e}^{-c_{5}\epsilon^{2}m}$:

provided that $m>(c_{6}\cdot\epsilon^{-2}\log\epsilon^{-1})n$ for some universal constants $c_{5},\,c_{6}>0$.

To have a sense of how large the quantities involved in Lemma 5 are, when $\gamma=0.7$ and $\rho=1/10$, it holds that

hence leading to $\zeta_{1}\approx 0.08$.

Having derived a lower bound for the first term in the right-hand side of (V-B1), it remains to deal with the second one.

Fix $\gamma>0$ and $\rho\leq 1/10$, and let $\mathcal{D}_{i}$, $\mathcal{K}_{i}$ be defined in (49), (50), respectively. For any constant $\epsilon>0$, there exists some universal constants $c_{5},\,c_{6}>0$ such that

holds with probability at least $1-2{\rm e}^{-c_{5}\epsilon^{2}m}$ provided that $m/n>c_{6}\cdot\epsilon^{-2}\log\epsilon^{-1}$ for some universal constants $c_{5},\,c_{6}>0$, where $\zeta_{2}^{\prime}=0.9748\sqrt{\rho\tau/(0.99\tau^{2}-\rho^{2})}$ with $\tau=(2+\gamma)/(1+\gamma)$.

With our TAF default parameters $\rho=1/10$ and $\gamma=0.7$, we have $\zeta_{2}^{\prime}\approx 0.2463$. Using (V-B1), (55), and (56), choosing $m/n$ exceeding some sufficiently large constant such that $c_{0}\leq c_{5}\epsilon^{2}$, and denoting $\zeta_{2}:=2\zeta_{2}^{\prime}(1+\gamma)/(2+\gamma)$, the following holds with probability exceeding $1-4{\rm e}^{-c_{0}m}$

for all $\bm{x}$ and $\bm{z}$ such that $\left\|\bm{h}\right\|/\left\|\bm{x}\right\|\leq\rho$ for $0<\rho\leq 1/10$ and any fixed $\gamma\geq 1/2$. This combined with (39) and (41) proves Proposition 2 for appropriately chosen $\mu>0$ and $\lambda>0$.

To conclude this section, an estimate for the working step size is provided next. Plugging the results of (45) and (47) into (40) suggests that

and also that $\lambda=2\left(1-\zeta_{1}-\zeta_{2}-2\epsilon\right)-\mu(1+\delta)^{4}\buildrel\triangle\over{=}\lambda_{0}$ in the local regularity condition in (39). Clearly, it holds that $0<\lambda<2(1-\zeta_{1}-\zeta_{2})$. Taking $\epsilon$ and $\delta$ to be sufficiently small, one obtains the feasible range of the step size for TAF

In particular, under default parameters in Algorithm 1, $\mu_{0}=0.8388$ and $\lambda_{0}=1.22$, thus concluding the proof of Theorem 1.

VI Conclusion

This paper developed a linear-time algorithm termed TAF for solving generally unstructured systems of random quadratic equations. Our TAF algorithm builds on three key ingredients: an orthogonality-promoting initialization, along with a simple yet effective gradient truncation rule, as well as scalable gradient-like iterations. Numerical tests using synthetic data and real images corroborate the superior performance of TAF over state-of-the-art solvers of the same type.

A few timely and pertinent future research directions are worth pointing out. First, in parallel with spectral initialization methods, the proposed orthogonality-promoting initialization can be applied for semidefinite optimization , matrix completion , as well as blind deconvolution . It is also interesting to investigate suitable gradient regularization rules in more general nonconvex optimization settings. Extending the theory to the more challenging case where $\bm{a}_{i}$’s are generated from the coded diffraction pattern model constitutes another meaningful direction.

Appendix A Proofs for Section V

By homogeneity of (28), it suffices to work with the case where $\|\bm{x}\|=1$. It is easy to check that

We next construct the following expression:

Upon letting $\bm{u}=\bm{x}^{\perp}$, the last inequality taken together with (61) concludes the proof of (29).

A-B Proof of Lemma 2

By the argument above, assume without loss of generality that $\bm{x}=\bm{e}_{1}$. Consider now the truncated vector $\bm{s}_{\backslash 1}|(\bm{s}^{\mathcal{T}}\bm{x})^{2}>\tau$, or equivalently, $\bm{s}_{\backslash 1}|s_{1}^{2}>\tau$. It is then clear that $\bm{s}_{\backslash 1}|s_{1}^{2}>\tau$ is bounded, and thus subgaussian; furthermore, the next hold

where (76b) is obtained as a submatrix of the first term in (A-B) since the second term $C_{2}\bm{e}_{1}\bm{e}_{1}^{\mathcal{T}}$ is removed.

The rows of $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0,\backslash 1}$ may therefore be viewed as independent realizations of the conditional random vector $\bm{s}_{\backslash 1}^{\mathcal{T}}|s_{1}^{2}>\tau$, with the threshold $\tau$ being the $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$-largest value in $\{y_{i}/\|\bm{a}_{i}\|^{2}\}_{i=1}^{m}$. Standard concentration inequalities on the sum of random positive semi-definite matrices composed of independent non-isotropic subgaussian rows [59, Remark 5.40] confirm that

for $n\geq 3$. Therefore, one readily concludes that

holds with probability at least $1-2{\rm e}^{-c_{K}n}$, provided that |\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\big{/}n exceeds some constant. Note that $c_{K}$ depends on the maximum subgaussian norm of rows of $\bm{S}$, and we assume without loss of generality $c_{K}\geq 1/2$. Hence, $\|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}\bm{u}\|^{2}$ in (29) is upper bounded simply by letting $\bm{u}=\bm{x}^{\perp}$ in (80).

A-C Proof of Lemma 3

We next pursue a meaningful lower bound for $\|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}\bm{x}\|^{2}$ in (31). When $\bm{x}=\bm{e}_{1}$, one has $\|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}\bm{x}\|^{2}=\|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}\bm{e}_{1}\|^{2}=\sum_{i=1}^{|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|}\bar{s}_{i,1}^{2}$, where $\{\bar{s}_{i,1}\}_{i=1}^{|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|}$ are entries of the first column of $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}$. It is further worth mentioning that all squared entries of any spherical random vector obey the Beta distribution with parameters $\alpha=\frac{1}{2}$, and $\beta=\frac{n-1}{2}$, i.e., $\bar{s}^{2}_{i,j}\sim{\rm Beta}\!\left(\frac{1}{2},\,\frac{n-1}{2}\right)$ for all $i,\,j$, [62, Lemma 2]. Although they have closed-form probability density functions (pdfs) that may facilitate deriving a lower bound, we take another route detailed as follows. A simple yet useful inequality is established first.

Given $m$ fractions obeying $1>\frac{p_{1}}{q_{1}}\geq\frac{p_{2}}{q_{2}}\geq\cdots\geq\frac{p_{m}}{q_{m}}>0$, in which $p_{i},\,q_{i}>0$, $\forall i\in[m]$, the following holds for all $1\leq k\leq m$

where $p_{[i]}$ denotes the $i$-th largest one among $\{p_{i}\}_{i=1}^{m}$, and hence, $q_{}$ is the maximum in $\{q_{i}\}_{i=1}^{m}$.

For any $k\in[m]$, according to the definition of $q_{[i]}$, it holds that $p_{}\geq p_{}\geq\cdots\geq p_{[k]}$, so $\frac{p_{}}{q_{}}\geq\frac{p_{}}{q_{}}\geq\cdots\geq\frac{p_{[k]}}{q_{}}$. Considering $q_{}\geq q_{i}$, $\forall i\in[m]$, and letting $j_{i}\in[m]$ be the index such that $p_{j_{i}}=p_{[i]}$, then $\frac{p_{j_{i}}}{q_{j_{i}}}=\frac{p_{[i]}}{q_{j_{i}}}\geq\frac{p_{[i]}}{q_{}}$ holds for any $i\in[k]$. Therefore, $\sum_{i=1}^{k}\frac{p_{j_{i}}}{q_{j_{i}}}=\sum_{i=1}^{k}\frac{p_{[i]}}{q_{j_{i}}}\geq\sum_{i=1}^{k}\frac{p_{[i]}}{q_{}}$. Note that $\left\{\frac{p_{[i]}}{q_{j_{i}}}\right\}_{i=1}^{k}$ comprise a subset of terms in $\left\{\frac{p_{i}}{q_{i}}\right\}_{i=1}^{m}$. On the other hand, according to our assumption, $\sum_{i=1}^{k}\frac{p_{i}}{q_{i}}$ is the largest among all sums of $k$ summands; hence, $\sum_{i=1}^{k}\frac{p_{i}}{q_{i}}\geq\sum_{i=1}^{k}\frac{p_{[i]}}{q_{j_{i}}}$ yields $\sum_{i=1}^{k}\frac{p_{i}}{q_{i}}\geq\sum_{i=1}^{k}\frac{p_{[i]}}{q_{}}$ concluding the proof. ∎

Without loss of generality and for simplicity of exposition, let us assume that indices of $\bm{a}_{i}$’s have been re-ordered such that

where $a_{i,1}$ denotes the first element of $\bm{a}_{i}$. Therefore, writing $\|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\bm{S}}_{0}\bm{e}_{1}\|^{2}=\sum_{i=1}^{|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|}a_{i,1}^{2}/\|\bm{a}_{i}\|^{2}$, the next task amounts to finding the sum of the $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$ largest out of all $m$ entities in (82). Applying the result (81) in Lemma 7 gives

in which $a_{[i],1}^{2}$ stands for the $i$-th largest entity in $\left\{a^{2}_{i,1}\right\}_{i=1}^{m}$.

To obtain an alternative bound, let us examine first the typical size of the maximum in $\left\{a_{i,1}^{2}\right\}_{i=1}^{m}$. Observe obviously that the modulus $\left|a_{i,1}\right|$ follows the half-normal distribution having the pdf $p(r)=\sqrt{{2}/{\pi}}\cdot{\rm e}^{-{r^{2}}/{2}}$, $r>0$, and it is easy to verify that

Choosing now $\xi:=\sqrt{2\log n}$ leads to

which holds with the proviso that $m/n$ is large enough, and the symbol $o(1)$ represents a small constant probability. Thus, provided that $m/n$ exceeds some large constant, the event $\max_{i\in[m]}a_{i,1}^{2}\geq 2\log n$ occurs with high probability. Hence, one may expect a tighter lower bound than $(1-\epsilon_{0})|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$, which is on the same order of $m$ under the assumption that $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|/m$ is about a constant.

Although $a_{i,1}^{2}$ obeys the Chi-square distribution with $k=1$ degrees of freedom, its cdf is rather complicated and does not admit a nice closed-form expression. A small trick is hence taken in the sequel. Assume without loss of generality that both $m$ and $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|$ are even. Grouping two consecutive $a_{[i],1}^{2}$’s together, introduce a new variable $\vartheta{[i]}:=a_{[2k-1],1}^{2}+a_{[2k],1}^{2}$, $\forall k\in[{m}/{2}]$, hence yielding a sequence of ordered numbers, i.e., $\vartheta_{}\geq\vartheta_{}\geq\cdots\geq\vartheta_{[m/2]}>0$. Then, one can equivalently write the wanted sum as

On the other hand, for i.i.d. standard normal random variables $\left\{a_{i,1}\right\}_{i=1}^{m}$, let us consider grouping randomly two of them and denote the corresponding sum of their squares by $\chi_{k}:=a_{k_{i},1}^{2}+a_{k_{j},1}^{2}$, where $k_{i}\neq k_{j}\in[m]$, and $k\in[m/2]$. It is self-evident that the $\chi_{k}$’s are identically distributed obeying the Chi-square distribution with $k=2$ degrees of freedom, having the pdf

and the following complementary cdf (ccdf)

Ordering all $\chi_{k}$’s, summing the $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|/2$ largest ones, and comparing the resultant sum with the one in (89) confirms that

Furthermore, applying the Hoeffding-type inequality [59, Proposition 5.10] and leveraging the convexity of the ccdf in (91), one readily establishes that

Taking without loss of generality \xi:=0.05\hat{\chi}_{|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|/2}=0.1\log\big{(}m\big{/}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|\big{)} gives

for some universal constants $c_{0},\,c_{\chi}>0$, and sufficiently large $n$ such that ${|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|}/{m}\gtrsim c_{\chi}>0$. The remaining part in this section assumes that this event occurs.

Choosing $\xi:=4\log n$ and substituting this into the ccdf in (91) leads to

Thus, $\left\{\vartheta_{i}\right\}_{i=1}^{m/2}$ are i.i.d. centered and bounded random variables following from the mean-subtraction and the upper-bound truncation. Further, according to the ccdf (91) and the definition of sub-exponential random variables [59, Definition 5.13], the terms $\left\{\vartheta_{i}\right\}_{i=1}^{m/2}$ are sub-exponential. Then, the following

holds with probability at least $1-2{\rm e}^{-c_{s}\min\left({\tau}/{K_{s}},{\tau^{2}}/{K_{s}^{2}}\right)}$, in which $c_{s}>0$ is a universal constant, and $K_{s}:=\max_{i\in[m/2]}\|\vartheta_{i}\|_{\psi_{1}}$ represents the maximum subexponential norm of the $\vartheta_{i}$’s.

Indeed, $K_{s}$ can be found as follows [59, Definition 5.13]:

Choosing $\tau:=8|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}_{0}|/(c_{s}m)\cdot\log^{2}n$ in (100) yields

for some small constant $\epsilon_{s}>0$, which holds with probability at least $1-m{\rm e}^{-n/2}-{\rm e}^{-c_{0}m}-1/n^{2}$ as long as $m/n$ exceeds some numerical constant and $n$ is sufficiently large. Therefore, combining (85), (92), and (A-C), one concludes that the following holds with high probability

Taking $\epsilon_{s}:=0.01$ without loss of generality concludes the proof of Lemma 3.

A-D Proof of Lemma 5

Let us first prove the argument for a fixed pair $\bm{h}$ and $\bm{x}$, such that $\bm{h}$ and $\bm{z}$ are independent of $\{\bm{a}_{i}\}_{i=1}^{m}$, and then apply a covering argument. To start, introduce a Lipschitz-continuous counterpart for the discontinuous indicator function [6, A.2]

By homogeneity and rotational invariance of normal distributions, it suffices to prove the case where $\bm{x}=\bm{e}_{1}$ and $\|\bm{h}\|/\|\bm{x}\|=\|\bm{h}\|\leq\rho$. According to (105), lower bounding the first term in (V-B1) can be achieved by lower bounding $\sum_{i=1}^{m}(\bm{a}_{i}^{\mathcal{T}}\bm{h})^{2}\chi_{E}\left(\left|1+\frac{\bm{a}_{i}^{\mathcal{T}}\bm{h}}{\bm{a}_{i}^{\mathcal{T}}\bm{x}}\right|\right)$ instead. To that end, let us find the mean of $\left(\bm{a}_{i}^{\mathcal{T}}\bm{h}\right)^{2}\chi_{E}\left(\left|1+\frac{\bm{a}_{i}^{\mathcal{T}}\bm{h}}{\bm{a}_{i}^{\mathcal{T}}\bm{x}}\right|\right)$. Note that $\left(\bm{a}_{i}^{\mathcal{T}}\bm{h}\right)^{2}$ and $\chi_{E}\left(\left|1+\frac{\bm{a}_{i}^{\mathcal{T}}\bm{h}}{\bm{a}_{i}^{\mathcal{T}}\bm{x}}\right|\right)$ are dependent. Introduce an orthonormal matrix $\bm{U}_{\bm{h}}$ that contains $\bm{h}^{\mathcal{T}}/\|\bm{h}\|$ as its first row, i.e.,

Direct application of the Berstein-type inequality [59, Proposition 5.16] confirms that for any $\epsilon>0$, the following

holds with probability at least $1-{\rm e}^{-c_{5}m\epsilon^{2}}$ for some numerical constant $c_{5}>0$ provided that $\epsilon\leq\|\varrho\|_{\psi_{1}}$ by assumption.

To obtain uniform control over all vectors $\bm{z}$ and $\bm{x}$ such that $\|\bm{z}-\bm{x}\|\leq\rho$, the net covering argument is applied [59, Definition 5.1]. Let $\mathcal{S}_{\epsilon}$ be an $\epsilon$-net of the unit sphere, $\mathcal{L}_{\epsilon}$ be an $\epsilon$-net of $[0,\,\rho]$, and define

Since the cardinality $\left|\mathcal{S}_{\epsilon}\right|\leq\left(1+2/\epsilon\right)^{n}$ [59, Lemma 5.2], then

due to the fact that $\rho/\epsilon<2/\epsilon<1+2/\epsilon$ for $0<\rho<1$.

Consider now any $\left(\bm{z},\,\bm{h},\,t\right)$ obeying $\|\bm{h}\|=t\leq\rho$. There exists a pair $\left(\bm{z}_{0},\,\bm{h}_{0},\,t_{0}\right)\in\mathcal{N}_{\epsilon}$ such that $\left\|\bm{z}-\bm{z}_{0}\right\|$, $\|\bm{h}-\bm{h}_{0}\|$, and $|t-t_{0}|$ are each at most $\epsilon$. Taking the union bound yields

with probability at least $1-\left(1+2/\epsilon\right)^{2n+1}{\rm e}^{-c_{5}\epsilon^{2}m}\geq 1-{\rm e}^{-c_{0}m}$, which follows by choosing $m$ such that $m\geq\left(c_{6}\cdot\epsilon^{-2}\log\epsilon^{-1}\right)n$ for some constant $c_{6}>0$.