Solving Random Quadratic Systems of Equations Is Nearly as Easy as Solving Linear Systems

Introduction

Imagine we are given a set of $m$ quadratic equations taking the form

This problem is combinatorial in nature as one can alternatively pose it as recovering the missing signs of $\langle\boldsymbol{a}_{i},\boldsymbol{x}\rangle$ from magnitude-only observations. As is well known, many classical combinatorial problems with Boolean variables may be cast as special instances of (1). As an example, consider the NP-hard stone problem in which we have $n$ stones each of weight $w_{i}>0$ ($1\leq i\leq n$), which we would like to divide into two groups of equal sum weight. Letting $x_{i}\in\{-1,1\}$ indicate which of the two groups the $i$th stone belongs to, one can formulate this problem as solving the following quadratic system

However simple this formulation may seem, even checking whether a solution to (2) exists or not is known to be NP-hard.

Continuing this motivating line of thought, in any real-world application recorded intensities are always corrupted by at least a small amount of noise so that observed data are only about $|\langle\boldsymbol{a}_{i},\boldsymbol{x}\rangle|^{2}$; i.e.

Although we present results for arbitrary noise distributions—even for non-stochastic noise—we shall pay particular attention to the Poisson data model, which assumes

The reason why this statistical model is of special interest is that it naturally describes the variation in the number of photons detected by an optical sensor in various imaging applications.

2 Nonconvex optimization

Under a stochastic noise model with independent samples, a first impulse for solving (3) is to seek the maximum likelihood estimate (MLE), namely,

modulo some constant offset. Unfortunately, the log-likelihood is usually not concave, thus making the problem of computing the MLE NP-hard in general.

To alleviate this computational intractability, several convex surrogates have been proposed that work particularly well when the design vectors $\{\boldsymbol{a}_{i}\}$ are chosen at random . The basic idea is to introduce a rank-one matrix $\boldsymbol{X}=\boldsymbol{x}\boldsymbol{x}^{*}$ to linearize the quadratic constraints, and then relax the rank-one constraint. Suppose we have Poisson data, then this strategy converts the problem into a convenient convex program:

Note that the log-likelihood function is augmented by the trace functional $\mathsf{Tr}(\cdot)$ whose role is to promote low-rank solutions. While such convex relaxation schemes enjoy intriguing performance guarantees in many aspects (e.g. they achieve minimal sample complexity and near-optimal error bounds for certain noise models), the computational cost typically far exceeds the order of $n^{3}$. This limits applicability to large-dimensional data.

This paper follows another route: rather than lifting the problem into higher dimensions by introducing matrix variables, this paradigm maintains its iterates within the vector domain and optimize the nonconvex objective directly (e.g. ). One promising approach along this line is the recently proposed two-stage algorithm called Wirtinger Flow (WF) . Simply put, WF starts by computing a suitable initial guess $\boldsymbol{z}^{(0)}$ using a spectral method, and then successively refines the estimate via an update rule that bears a strong resemblance to a gradient descent scheme, namely,

WF achieves exact recovery from $m=O\left(n\log n\right)$ quadratic equations when there is no noise; The standard notation $f(n)=O\left(g(n)\right)$ or $f(n)\lesssim g(n)$ (resp. $f(n)=\Omega\left(g(n)\right)$ or $f(n)\gtrsim g(n)$) means that there exists a constant $c>0$ such that $\left|f(n)\right|\leq c|g(n)|$ (resp. $|f(n)|\geq c\left|g(n)\right|$). $f(n)\asymp g(n)$ means that there exist constants $c_{1},c_{2}>0$ such that $c_{1}|g(n)|\leq|f(n)|\leq c_{2}|g(n)|$.

WF attains $\epsilon$-accuracy—in a relative sense—within $O(mn^{2}\log({1}/{\epsilon}))$ time (or flops);

In the presence of Gaussian noise, WF is stable and converges to the MLE as shown in .

While these results formalize the advantages of WF, the computational complexity of WF is still much larger than the best one can hope for. Moreover, the statistical guarantee in terms of the sample complexity is weaker than that achievable by convex relaxations. M. Soltanolkotabi recently informed us that the sample complexity of WF may be improved if one employs a better initialization procedure.

3 This paper: Truncated Wirtinger Flow

This paper develops an efficient linear-time algorithm, which also enjoys near-optimal statistical guarantees. Following the spirit of WF, we propose a novel procedure called Truncated Wirtinger Flow (TWF) adopting a more adaptive gradient flow. Informally, TWF proceeds in two stages:

Initialization: compute an initial guess $\boldsymbol{z}^{(0)}$ by means of a spectral method applied to a subset $\mathcal{T}_{0}$ of the observations $\{y_{i}\}$;

for some index subset $\mathcal{T}_{t+1}\subseteq\{1,\cdots,m\}$ determined by $\boldsymbol{z}^{(t)}$.

Firstly, we regularize both the initialization and the gradient flow in a data-dependent fashion by operating only upon some iteration-varying index subsets $\mathcal{T}_{t}$. This is a distinguishing feature of TWF in comparison to WF and other gradient descent variants. In words, $\mathcal{T}_{t}$ corresponds to those data $\{y_{i}\}$ whose resulting spectral or gradient components are in some sense not excessively large; see Sections 2 and 3 for details. As we shall see later, the main point is that this careful data trimming procedure gives us a tighter initial guess and more stable search directions.

Secondly, we recommend that the step size $\mu_{t}$ is either taken as some appropriate constant or determined by a backtracking line search. For instance, under appropriate conditions, we can take $\mu_{t}=0.2$ for all $t$.

Hence, $\boldsymbol{A}\boldsymbol{z}$ gives $\boldsymbol{v}$ and $\boldsymbol{A}^{\top}\boldsymbol{v}$ the desired regularized gradient.

A detailed specification of the algorithm is deferred to Section 2.

4 Numerical surprises

To give the readers a sense of the practical power of TWF, we present here three illustrative numerical examples. Since it is impossible to recover the global sign—i.e. we cannot distinguish $\boldsymbol{x}$ from $-\boldsymbol{x}$—we will evaluate our solutions to the quadratic equations through the distance measure put forth in representing the Euclidean distance modulo a global sign: for complex-valued signals,

The first numerical example concerns the following two problems under noiseless real-valued data:

Apparently, (a) involves solving a linear system of equations (or a linear least squares problem), while (b) is tantamount to solving a quadratic system. Arguably the most popular method for solving large-scale least squares problems is the conjugate gradient (CG) method applied to the normal equations. We are going to compare the computational efficiency between CG (for solving least squares) and TWF with a step size $\mu_{t}\equiv 0.2$ (for solving a quadratic system). Set $m=8n$ and generate $\boldsymbol{x}\sim\mathcal{N}({\bf 0},\boldsymbol{I})$ and $\boldsymbol{a}_{i}\sim\mathcal{N}({\bf 0},\boldsymbol{I})$, $1\leq i\leq m$, independently. This gives a matrix $\boldsymbol{A}^{\top}\boldsymbol{A}$ with a low condition number equal to about $(1+\sqrt{1/8})^{2}/(1-\sqrt{1/8})^{2}\approx 4.38$ by the Marchenko-Pastur law. Therefore, this is an ideal setting for CG as it converges extremely rapidly [32, Theorem 38.5]. Fig. 1 shows the relative estimation error of each method as a function of the iteration count, where TWF is seeded through 10 power iterations. For ease of comparison, we illustrate the iteration counts in different scales so that 4 TWF iterations are equivalent to 1 CG iteration.

Recognizing that each iteration of CG and TWF involves two matrix vector products $\boldsymbol{A}\boldsymbol{z}$ and $\boldsymbol{A}^{\top}\boldsymbol{v}$, for such a design we reach a suprising observation:

Even when all phase information is missing, TWF is capable of solving a quadratic system of equations only about 4 times slower than solving a least squares problem of the same size!

To illustrate the applicability of TWF on real images, we turn to testing our algorithm on a digital photograph of Stanford main quad containing $320\times 1280$ pixels. We consider a type of measurements that falls under the category of coded diffraction patterns (CDP) and set

Here, $\boldsymbol{F}$ stands for a discrete Fourier transform (DFT) matrix, and $\boldsymbol{D}^{(l)}$ is a diagonal matrix whose diagonal entries are independently and uniformly drawn from $\{1,-1,j,-j\}$ (phase delays). In phase retrieval, each $\boldsymbol{D}^{(l)}$ represents a random mask placed after the object so as to modulate the illumination patterns. When $L$ masks are employed, the total number of quadratic measurements is $m=nL$. In this example, $L=12$ random coded patterns are generated to measure each color band (i.e. red, green, or blue) separately. The experiment is carried out on a MacBook Pro equipped with a 3 GHz Intel Core i7 and 16GB of memory. We run 50 iterations of the truncated power method for initialization, and 50 regularized gradient iterations, which in total costs 43.9 seconds or 2400 FFTs for each color band. The relative error after regularized spectral initialization and after 50 TWF iterations are 0.4773 and $2.16\times 10^{-5}$, respectively, with the recovered images displayed in Fig. 2. In comparison, the spectral initialization using 50 untruncated power iterations returns an image of relative error 1.409, which is almost like a random guess and extremely far from the truth.

While the above experiments concern noiseless data, the numerical surprise extends to the noisy realm. Suppose the data are drawn according to the Poisson noise model (4), with $\boldsymbol{a}_{i}\sim\mathcal{N}\left({\bf 0},\boldsymbol{I}\right)$ independently generated. Fig. 3 displays the empirical relative mean-square error (MSE) of TWF as a function of the signal-to-noise ratio (SNR), where the relative MSE for an estimate $\hat{\boldsymbol{x}}$ and the SNR are defined asTo justify the definition of SNR, note that the signals and noise are captured by $\mu_{i}=(\boldsymbol{a}_{i}^{\top}\boldsymbol{x})^{2}$ and $y_{i}-\mu_{i}$, $1\leq i\leq m$, respectively. The ratio of the signal power to the noise power is therefore $\frac{\sum_{i=1}^{m}\mu_{i}^{2}}{\sum_{i=1}^{m}{\bf Var}[y_{i}]}=\frac{\sum_{i=1}^{m}|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|^{4}}{\sum_{i=1}^{m}|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|^{2}}\approx\frac{3m\|\boldsymbol{x}\|^{4}}{m\|\boldsymbol{x}\|^{2}}=3\|\boldsymbol{x}\|^{2}.$

Fig. 3 illustrates the empirical performance for this ideal problem. The plots demonstrate that even when all phases are erased, TWF yields a solution of nearly the best possible quality, since it only incurs an extra $1.5$ dB loss compared to ideal MLE computed with all true phases revealed. This phenomenon arises regardless of the SNR!

5 Main results

The preceding numerical discoveries unveil promising features of TWF in three aspects: (1) exponentially fast convergence; (2) exact recovery from noiseless data with sample complexity $O\left(n\right)$; (3) nearly minimal mean-square loss in the presence of noise. This paper offers a formal substantiation of all these findings. To this end, we assume a tractable model in which the design vectors $\boldsymbol{a}_{i}$’s are independent Gaussian:

For concreteness, our results are concerned with TWF designed based on the Poisson log-likelihood function

As explained below, we can often take $\mu_{0}\approx 0.3$.

As will be made precise in Section 5 (and in particular Proposition 1), one can take

for some small quantities $\zeta_{1},\zeta_{2}$ and some predetermined threshold $\alpha_{h}$ that is usually taken to be $\alpha_{h}\geq 5$. Under appropriate conditions, one can treat $\mu_{0}$ as $\mu_{0}\approx 0.3$.

We emphasize that enhanced performance vis-à-vis WF is not the result of a sharper analysis, but rather, the result of key algorithmic changes. In both the initialization and iterative refinement stages, TWF proceeds in a more prudent manner by means of proper regularization, which effectively trims away those components that are too influential on either the initial guess or search directions, thus reducing the volatility of each movement. With a tighter initialization and better-controlled search directions in place, we take the step size in a far more liberal fashion—which is some constant bounded away from 0—compared to a step size which is $O(1/n)$ as explained in . In fact, what enables the movement to be more aggressive is exactly the cautious choice of $\mathcal{T}_{t}$, which precludes adverse effects from high-leverage samples.

To be broadly applicable, the proposed algorithm must guarantee reasonably faithful estimates in the presence of noise. Suppose that

where $\eta_{i}$ represents an error term. We claim that TWF is stable against additive noise, as demonstrated in the theorem below.

Consider the noisy case (14). Suppose that the step size $\mu_{t}$ is either taken to be a positive constant $\mu_{t}\equiv\mu$ or chosen via a backtracking line search. If

then with probability at least $1-c_{2}\exp\left(-c_{3}m\right)$, the truncated Wirtinger Flow estimates (Algorithm 1 with parameters specified in Table 1) satisfy

Under the Poisson noise model (4), one has

with probability approaching one, provided that $\|\boldsymbol{x}\|\geq\log^{1.5}m$.

establishes stability estimates using the WF approach under Gaussian noise. There, the sample and computational complexities are still on the order of $n\log n$ and $mn^{2}$ respectively whereas the computational complexity in Theorem 2 is linear, i.e. on the order of $mn$.

Theorem 2 essentially reveals that the estimation error of TWF rapidly shrinks to $O\left(\frac{\|\boldsymbol{\eta}\|/\sqrt{m}}{\|\boldsymbol{x}\|}\right)$ within logarithmic iterations. Put another way, since the SNR for the model (14) is captured by

we immediately arrive at an alternative form of the performance guarantee:

revealing the stability of TWF as a function of SNR. We emphasize that this estimate holds for any error term $\boldsymbol{\eta}$—i.e. any noise structure, even deterministic. This being said, specializing this estimate to the Poisson noise model (4) with $\left\|\boldsymbol{x}\right\|\gtrsim\log^{1.5}m$ gives an estimation error that will eventually approach a numerical constant, independent of $n$ and $m$.

Encouragingly, this is already the best statistical guarantee any algorithm can achieve. We formalize this claim by deriving a fundamental lower bound on the minimax estimation error.

Suppose that $\boldsymbol{a}_{i}\sim\mathcal{N}({\bf 0},\boldsymbol{I})$, $m=\kappa n$ for some fixed $\kappa$ independent of $n$, and $n$ is sufficiently large. For any $K\geq\log^{1.5}m$, defineHere, 0.1 can be replaced by any positive constant within (0, 1/2).

With probability approaching one, the minimax risk under the Poisson model (4) obeys

where the infimum is over all estimator $\hat{\boldsymbol{x}}$. Here, $\varepsilon_{1}>0$ is a numerical constant independent of $n$ and $m$.

When the number $m$ of measurements is proportional to $n$ and the energy of the planted solution exceeds $\log^{3}m$, Theorem 3 asserts that there exists absolutely no estimator that can achieve an estimation error that vanishes as $n$ increases. This lower limit matches the estimation error of TWF, which corroborates the optimality of TWF under noisy data.

Recall that in many optical imaging applications, the output data we collect are the intensities of the diffractive waves scattered by the sample or specimen under study. The Poisson noise model employs the input $\boldsymbol{x}$ and output $\boldsymbol{y}$ to describe the numbers of photons diffracted by the specimen and detected by the optical sensor, respectively. Each specimen needs to be sufficiently illuminated in order for the receiver to sense the diffracted light. In such settings, the low-intensity regime $\|\boldsymbol{x}\|\leq\log^{1.5}m$ is of little practical interest as it corresponds to an illumination with just very few photons. We forego the details.

It is worth noting that apart from WF, various other nonconvex procedures have been proposed as well for phase retrieval, including the error reduction schemes dating back to Gerchberg-Saxton and Fienup , iterated projections , alternating minimization , generalized approximate message passing , Kaczmarz method , and greedy methods that exploit additional sparsity constraint , to name just a few. While these paradigms enjoy favorable empirical behavior, most of them fall short of theoretical support, except for a version of alternating minimization (called AltMinPhase) that requires fresh samples for each iteration. In comparison, AltMinPhase attains $\epsilon$-accuracy when the sample complexity exceeds the order of $n\log^{3}n+n\log^{2}n\log({1}/{\epsilon})$, which is at least a factor of $\log^{3}n$ from optimal and is empirically largely outperformed by the variant that reuses all samples. In contrast, our algorithm uses the same set of samples all the time and is therefore practically appealing. Furthermore, none of these algorithms come with provable stability guarantees, which are particularly important in most realistic scenarios. Numerically, each iteration of Fienup’s algorithm (or alternating minimization) involves solving a least squares problem, and the algorithm converges in tens or hundreds of iterations. This is computationally more expensive than TWF, whose computational complexity is merely about 4 times that of solving a least squares problem. Interesting readers are referred to for a comparison of several non-convex schemes, and for a discussion of other alternative approaches (e.g. ) and performance lower bounds (e.g. ).

Algorithm: Truncated Wirtinger Flow

For independent samples, the gradient of the real-valued Poisson log-likelihood obeys

where $\nu_{i}$ represents the weight assigned to each $\boldsymbol{a}_{i}$. This forms the descent direction of WF updates.

Hence, to remedy the aforementioned stability issue, it would be natural to separate the small fraction of abnormal gradient components by regularizing the weights $\nu_{i}$, possibly via data-dependent trimming rules. This gives rise to the update rule of TWF:

for some trimming criteria specified by $\mathcal{E}_{1}^{i}\left(\cdot\right)$ and $\mathcal{E}_{2}^{i}\left(\cdot\right)$. In our algorithm, we take $\mathcal{E}_{1}^{i}\left(\boldsymbol{z}\right)$ and $\mathcal{E}_{2}^{i}\left(\boldsymbol{z}\right)$ to be two collections of events given by

where $\alpha_{z}^{\text{lb}}$, $\alpha_{z}^{\text{ub}}$, $\alpha_{z}$ are predetermined thresholds. To keep notation light, we shall use $\mathcal{E}_{1}^{i}$ and $\mathcal{E}_{2}^{i}$ rather than $\mathcal{E}_{1}^{i}\left(\boldsymbol{z}\right)$ and $\mathcal{E}_{2}^{i}\left(\boldsymbol{z}\right)$ whenever it is clear from context.

We emphasize that the above trimming procedure simply throws away those components whose weights $\nu_{i}$’s fall outside some confidence range, so as to remove the influence of outlier components. To achieve this, we regularize both the numerator and denominator of $\nu_{i}$ by enforcing separate trimming rules. Recognize that for any fixed $\boldsymbol{z}$, the denominator obeys

leading up to the rule (24). Regarding the numerator, by the law of large numbers one would expect

and hence it is natural to regularize the numerator by ensuring

As a remark, we include an extra term ${\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\right|}/{\left\|\boldsymbol{z}\right\|}$ in (25) to sharpen the theory, but all our results continue to hold (up to some modification of constants) if we drop this term in (25). Detailed procedures are summarized in Algorithm 1 Careful readers might note that we include some extra factor $\frac{\sqrt{n}}{\|\boldsymbol{a}_{i}\|}$ (which is approximately 1 in the Gaussian model) in Algorithm 1. This occurs since we present Algorithm 1 in a more general fashion that applies beyond the model $\boldsymbol{a}_{i}\sim\mathcal{N}({\bf 0},\boldsymbol{I})$, but all results / proofs continue to hold in the presence of this extra term. .

The proposed paradigm could be counter-intuitive at first glance, since one might expect the larger terms to be better aligned with the desired search direction. The issue, however, is that the large terms are extremely volatile and could have too high of a leverage on the descent directions. In contrast, TWF discards these high-leverage data, which slightly increases the bias but remarkably reduces the variance of the descent direction. We expect such gradient regularization and variance reduction schemes to be beneficial for solving a broad family of nonconvex problems.

2 Truncated spectral initialization

In order for the gradient stage to converge rapidly, we need to seed it with a suitable initialization. One natural alternative is the spectral method adopted in , which amounts to computing the leading eigenvector of $\widetilde{\boldsymbol{Y}}:=\frac{1}{m}\sum_{i=1}^{m}y_{i}\boldsymbol{a}_{i}\boldsymbol{a}_{i}^{\top}$. This arises from the observation that when $\boldsymbol{a}_{i}\sim\mathcal{N}\left({\bf 0},\boldsymbol{I}\right)$ and $\left\|\boldsymbol{x}\right\|=1$,

whose leading eigenvector is exactly $\boldsymbol{x}$ with an eigenvalue of 3.

Unfortunately, this spectral technique converges to a good initial point only when $m\gtrsim n\log n$, due to the fact that $(\boldsymbol{a}_{i}^{\top}\boldsymbol{x})^{2}\boldsymbol{a}_{i}\boldsymbol{a}_{i}^{\top}$ is heavy-tailed, a random quantity which does not have a moment generating function. To be more precise, consider the noiseless case $y_{i}=|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|^{2}$ and recall that $\max_{i}y_{i}\approx 2\log m$. Letting $k=\arg\max_{i}y_{i}$, one can calculate

which is much larger than $\boldsymbol{x}^{\top}\widetilde{\boldsymbol{Y}}\boldsymbol{x}=3$ unless $m/n$ is very large. This tells us that in the regime where $m\asymp n$, there exists some unit vector $\boldsymbol{a}_{k}/\|\boldsymbol{a}_{k}\|$ that is closer to the leading eigenvector of $\widetilde{\boldsymbol{Y}}$ than $\boldsymbol{x}$. This phenomenon happens because the summands of $\widetilde{\boldsymbol{Y}}$ have huge tails so that even one large term could end up dominating the empirical sum, thus preventing the spectral method from returning a meaningful initial guess.

Notably, the aforementioned drawback of the spectral method is not merely a theoretical concern but rather a substantial practical issue. We have seen this in Fig. 2 (main quad example) showing the enormous advantage of truncated spectral initialization. This is also further illustrated in Fig. 6, which compares the empirical efficiency of both methods with $\alpha_{y}=3$ set to be the truncation threshold. For both Gaussian designs and CDP models, the empirical loss incurred by the original spectral method increases as $n$ grows, which is in stark constrast to the truncated spectral method that achieves almost identical accuracy over the same range of $n$.

3 Choice of algorithmic parameters

One implementation detail to specify is the step size $\mu_{t}$ at each iteration $t$. There are two alternatives that work well in both theory and practice:

Backtracking line search with truncated objective. This strategy performs a line search along the descent direction

and determines an appropriate step size that guarantees a sufficient improvement. In contrast to the conventional search strategy that determines the sufficient progress with respect to the true objective function, we propose to evaluate instead a regularized version of the objective function. Specifically, put

Then the backtracking line search proceeds as

where $\beta\in(0,1)$ is some pre-determined constant;

When a backtracking line search is adopted, an extra parameter $\alpha_{p}$ is needed, which we take to be $\alpha_{p}\geq 5$. In all theory presented herein, we assume that the parameters fall within the range singled out in Table 1.

Why TWF works?

Before proceeding, it is best to develop an intuitive understanding of the TWF iterations. We start with a notation representing the (unrecoverable) global phase for real-valued data

despite the global phase uncertainty. For simplicity of presentation, we shall drop the phase term by letting $\boldsymbol{z}$ be $e^{-j\phi\left(\boldsymbol{z}\right)}\boldsymbol{z}$ and setting $\boldsymbol{h}=\boldsymbol{z}-\boldsymbol{x}$, whenever it is clear from context.

The first object to consider is the descent direction. To this end, we find it convenient to work with a fixed $\boldsymbol{z}$ independent of the design vectors $\boldsymbol{a}_{i}$, which is of course heuristic but helpful in developing some intuition. Rewrite

where (i) follows from the identity $a^{2}-b^{2}=(a+b)(a-b)$. The first component of (36), which on average gives $-4\boldsymbol{h}$, makes a good search direction when averaged over all the observations $i=1,\ldots,m$. The issue is that the other term $\boldsymbol{r}_{i}$—which is in general non-integrable—could be devastating. The reason is that $\boldsymbol{a}_{i}^{\top}\boldsymbol{z}$ could be arbitrarily small, thus resulting in an unbounded $\boldsymbol{r}_{i}$. As a consequence, a non-negligible portion of the $\boldsymbol{r}_{i}$’s may exert a very strong influence on the descent direction in an undesired manner.

Such an issue can be prevented if one can detect and separate those gradient components bearing abnormal $\boldsymbol{r}_{i}$’s. Since we cannot observe the individual components of the decomposition (36), we cannot reject indices with large values of $\boldsymbol{r}_{i}$ directly. Instead, we examine each gradient component as a whole and discard it if its size is not absolutely controlled. Fortunately, such a strategy is sufficient to ensure that most of the contribution from the regularized gradient comes from the first component of (36), namely, $-4(\boldsymbol{a}_{i}^{\top}\boldsymbol{h})\boldsymbol{a}_{i}$. As will be made precise in Proposition 2 and Lemma 7, the regularized gradient obeys

Here, one has $(4-\epsilon)\|\boldsymbol{h}\|^{2}$ in (37) instead of $4\|\boldsymbol{h}\|^{2}$ to account for the bias introduced by adaptive trimming, where $\epsilon$ is small as long as we only throw away a small fraction of data. Looking at (37) and (38) we see that the search direction is sufficiently aligned with the deviation $-\boldsymbol{h}=\boldsymbol{x}-\boldsymbol{z}$ of the current iterate; i.e. they form a reasonably good angle that is bounded away from $90^{\circ}$. Consequently, $\boldsymbol{z}$ is expected to be dragged towards $\boldsymbol{x}$ provided that the step size is appropriately chosen.

holds for all $\boldsymbol{z}$ obeying $\left\|\boldsymbol{z}-\boldsymbol{x}\right\|\leq\epsilon\|\boldsymbol{x}\|$, where $0<\epsilon<1$ is some constant. Such an $\epsilon$-ball around $\boldsymbol{x}$ forms a basin of attraction. Formally, under $\mathsf{RC}\left(\mu,\lambda,\epsilon\right)$, a little algebra gives

for any $\boldsymbol{z}$ with $\left\|\boldsymbol{z}-\boldsymbol{x}\right\|\leq\epsilon$. In words, the TWF update rule is locally contractive around the planted solution, provided that $\mathsf{RC}\left(\mu,\lambda,\epsilon\right)$ holds for some nonzero $\mu$ and $\lambda$. Apparently, Conditions (37) and (38) already imply the validity of $\mathsf{RC}$ for some constants $\mu,\lambda\asymp 1$ when $\|\boldsymbol{h}\|/\|\boldsymbol{z}\|$ is reasonably small, which in turn allows us to take a constant step size $\mu$ and enables a constant contraction rate $1-\mu\lambda$.

Numerical experiments

In this section, we report additional numerical results to verify the practical applicability of TWF. In all numerical experiments conducted in the current paper, we set

This is a concrete combination of parameters satisfying our condition (30). Unless otherwise noted, we employ 50 power iterations for initialization, adopt a fixed step size $\mu_{t}\equiv 0.2$ when updating TWF iterates, and set the maximum number of iterations to be $T=1000$ for the iterative refinement stage.

To see how special the real-valued Gaussian designs are to our theoretical finding, we perform experiments on two other types of measurement models. In the first, TWF is applied to complex-valued data by generating $\boldsymbol{a}_{i}\sim\mathcal{N}\left({\bf 0},\frac{1}{2}\boldsymbol{I}\right)+j\mathcal{N}\left({\bf 0},\frac{1}{2}\boldsymbol{I}\right)$. The other is the model of coded diffraction patterns described in (9). Fig. 9 depicts the average success rate for both types of measurements over 100 Monte Carlo trials, indicating that $m>4.5n$ and $m\geq 6n$ are often sufficient under complex-valued Gaussian and CDP models, respectively.

For the sake of comparison, we also report the empirical performance of WF in all the above settings, where the step size is set to be the default choice of , that is, $\mu_{t}=\min\{1-e^{-t/330},0.2\}$. As can be seen, the empirical success rates of TWF outperform WF when $T=1000$ under Gaussian models, suggesting that TWF either converges faster or exhibits better phase transition behavior.

Another series of experiments has been carried out to demonstrate the stability of TWF when the number $m$ of quadratic equations varies. We consider the case where $n=1000$, and vary the SNR (cf. (10)) from 15 dB to 55dB. The design vectors are real-valued independent Gaussian $\boldsymbol{a}_{i}\sim\mathcal{N}\left({\bf 0},\boldsymbol{I}\right)$, while the measurements $y_{i}$ are generated according to the Poisson noise model (4). Fig. 10 shows the relative mean square error—in the dB scale—as a function of SNR, when averaged over 100 independent runs. For all choices of $m$, the numerical experiments demonstrate that the relative MSE scales inversely proportional to SNR, which matches our stability guarantees in Theorem 2 (since we observe that on the dB scale, the slope is about -1 as predicted by the theory (19)).

Exact recovery from noiseless data

This section proves the theoretical guarantees of TWF in the absence of noise (i.e. Theorem 1). We separate the noiseless case mainly out of pedagogical reasons, as most of the steps carry over to the noisy case with slight modification.

with high probability, provided that $m/n$ exceeds some numerical constant. With this in place, it suffices to demonstrate that the TWF update rule is locally contractive, as stated in the following proposition.

Consider the noiseless case (1). Under the condition (30), there exist some universal constants $0<\rho_{0}<1$ and $c_{0},c_{1},c_{2}>0$ such that with probability exceeding $1-c_{1}\exp\left(-c_{2}m\right)$,

provided that $m\geq c_{0}n$ and that $\mu$ is some constant obeying

Proposition 1 reveals the monotonicity of the estimation error: once entering a neighborhood around $\boldsymbol{x}$ of a reasonably small size, the iterative updates will remain within this neighborhood all the time and be attracted towards $\boldsymbol{x}$ at a geometric rate.

As shown in Section 3, under the hypothesis $\mathsf{RC}\left(\mu,\lambda,\epsilon\right)$ one can conclude

Thus, everything now boils down to showing $\mathsf{RC}\left(\mu,\lambda,\epsilon\right)$ for some constants $\mu,\lambda,\epsilon>0$. This occupies the rest of this section.

Before proceeding, we gather a few properties of the events $\mathcal{E}_{1}^{i}$ and $\mathcal{E}_{2}^{i}$, which will prove crucial in establishing $\mathsf{RC}\left(\mu,\lambda,\epsilon\right)$. To begin with, recall that the truncation level given in $\mathcal{E}_{2}^{i}$ depends on $\frac{1}{m}\left\|\mathcal{A}\left(\boldsymbol{x}\boldsymbol{x}^{\top}-\boldsymbol{z}\boldsymbol{z}^{\top}\right)\right\|_{1}$. Instead of working with this random variable directly, we use deterministic quantities that are more amenable to analysis. Specifically, we claim that $\frac{1}{m}\left\|\mathcal{A}\left(\boldsymbol{x}\boldsymbol{x}^{\top}-\boldsymbol{z}\boldsymbol{z}^{\top}\right)\right\|_{1}$ offers a uniform and orderwise tight estimate on $\left\|\boldsymbol{h}\right\|\left\|\boldsymbol{z}\right\|$, which can be seen from the following two facts.

Fix $\zeta\in(0,1)$. If $m>c_{0}n\zeta^{-2}\log\frac{1}{\zeta}$, then with probability at least $1-C\exp(-c_{1}\zeta^{2}m)$,

Since [6, Lemma 3.1] already establishes the upper bound, it suffices to prove the lower tail bound. Consider all symmetric rank-2 matrices $\boldsymbol{M}$ with eigenvalues $1$ and $-t$ for some $-1\leq t\leq 1$. When $t\in$, it has been shown in the proof of [6, Lemma 3.2] that with high probability,

for all such rank-2 matrices $\boldsymbol{M}$, where $f\left(t\right):=\frac{2}{\pi}\left\{2\sqrt{t}+\left(1-t\right)\left(\pi/2-2\text{arc}\tan(\sqrt{t})\right)\right\}$. The lower bound in this case can then be justified by recognizing that $f\left(t\right)/\sqrt{1+t^{2}}\geq 0.9$ for all $t\in$, as illustrated in Fig. 11. The case where $t\in$ is an immediate consequence from [6, Lemma 3.1]. ∎

Take $\boldsymbol{h}=\boldsymbol{z}-\boldsymbol{x}$ and write

When $\|\boldsymbol{h}\|<\frac{1}{2}\|\boldsymbol{z}\|$, the Cauchy-Schwarz inequality gives

Taken together the above two facts demonstrate that with probability $1-\exp\left(-\Omega\left(m\right)\right)$,

holds simultaneously for all $\boldsymbol{z}$ and $\boldsymbol{x}$ satisfying $\left\|\boldsymbol{h}\right\|\leq\frac{1}{11}\left\|\boldsymbol{z}\right\|$. Conditional on (50), the inclusion

holds with respect to the following events

The point of introducing these new events is that the $\mathcal{E}_{3}^{i}$’s (resp. $\mathcal{E}_{4}^{i}$’s) are statistically independent for any fixed $\boldsymbol{x}$ and $\boldsymbol{z}$ and are, therefore, easier to work with.

Note that each $\mathcal{E}_{3}^{i}$ (resp. $\mathcal{E}_{4}^{i}$) is specified by a quadratic inequality. A closer inspection reveals that in order to satisfy these quadratic inequalities, the quantity $\boldsymbol{a}_{i}^{\top}\boldsymbol{h}$ must fall within two intervals centered around and $2\boldsymbol{a}_{i}^{\top}\boldsymbol{z}$, respectively. One can thus facilitate analysis by decoupling each quadratic inequality of interest into two simple linear inequalities, as stated in the following lemma.

2 Proof of the regularity condition

implying $|v_{i}|\lesssim\|\boldsymbol{h}\|$ and hence $\|\boldsymbol{v}\|\lesssim\sqrt{m}\|\boldsymbol{h}\|$. The Marchenko–Pastur law gives $\|\boldsymbol{A}\|\lesssim\sqrt{m}$, whence

A more refined estimate will be provided in Lemma 7.

The above argument essentially tells us that to establish $\mathsf{RC}$, it suffices to verify a uniform lower bound of the form

as formally derived in the following proposition.

Consider the noise-free measurements $y_{i}=|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|^{2}$ and any fixed constant $\epsilon>0$. Under the condition (30), if $m>c_{1}n$, then with probability exceeding $1-C\exp\left(-c_{0}m\right)$,

Here, $c_{0},c_{1},C>0$ are some universal constants, and $\zeta_{1}$ and $\zeta_{2}$ are defined in (30).

The basic starting point is the observation that $(\boldsymbol{a}_{i}^{\top}\boldsymbol{z})-(\boldsymbol{a}_{i}^{\top}\boldsymbol{x})^{2}=(\boldsymbol{a}_{i}^{\top}\boldsymbol{h})(2\boldsymbol{a}_{i}^{\top}\boldsymbol{z}-\boldsymbol{a}_{i}^{\top}\boldsymbol{h})$ and hence

One would expect the contribution of the second term of (62) (which is a second-order quantity) to be small as $\left\|\boldsymbol{h}\right\|/\left\|\boldsymbol{z}\right\|$ decreases.

To facilitate analysis, we rewrite (62) in terms of the more convenient events $\mathcal{D}_{\gamma}^{i,1}$ and $\mathcal{D}_{\gamma}^{i,2}$. Specifically, the inclusion property (51) together with Lemma 3 reveals that

where the parameters $\gamma_{3},\gamma_{4}$ are given by

This taken collectively with the identity (62) leads to a lower estimate

leaving us with three quantities in the right-hand side to deal with. We pause here to explain and compare the influences of these three terms.

To begin with, as long as the trimming step does not discard too many data, the first term should be close to $\frac{2}{m}\sum_{i}|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}|^{2}$, which approximately gives $2\|\boldsymbol{h}\|^{2}$ from the law of large numbers. This term turns out to be dominant in the right-hand side of (65) as long as $\|\boldsymbol{h}\|/\|\boldsymbol{z}\|$ is reasonably small. To see this, please recognize that the second term in the right-hand side is $O(\|\boldsymbol{h}\|^{3}/\|\boldsymbol{z}\|)$, simply because both $\boldsymbol{a}_{i}^{\top}\boldsymbol{h}$ and $\boldsymbol{a}_{i}^{\top}\boldsymbol{z}$ are absolutely controlled on $\mathcal{D}_{\gamma_{4}}^{i,1}\cap\mathcal{E}_{1}^{i}$. However, $\mathcal{D}_{\gamma_{4}}^{i,2}$ does not share such a desired feature. By the very definition of $\mathcal{D}_{\gamma_{4}}^{i,2}$, each nonzero summand of the last term of (65) must obey $\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right|\approx 2\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\right|$ and, therefore, $\frac{\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right|^{3}}{\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\right|}{\bf 1}_{\mathcal{E}_{1}^{i}\cap\mathcal{D}_{\gamma_{4}}^{i,2}}$ is roughly of the order of $\|\boldsymbol{z}\|^{2}$; this could be much larger than our target level $\left\|\boldsymbol{h}\right\|^{2}$. Fortunately, $\mathcal{D}_{\gamma_{4}}^{i,2}$ is a rare event, thus precluding a noticable influence upon the descent direction. All of this is made rigorous in Lemma 4 (first term), Lemma 5 (second term) and Lemma 6 (third term) together with subsequent analysis.

Fix $\gamma>0$, and let $\mathcal{E}_{1}^{i}$ and $\mathcal{D}_{\gamma}^{i,1}$ be defined in (24) and (55), respectively. Set

where $\xi\sim\mathcal{N}\left(0,1\right)$. For any $\epsilon>0$, if $m>c_{1}n\epsilon^{-2}\log\epsilon^{-1}$, then with probability at least $1-C\exp(-c_{0}\epsilon^{2}m)$,

We now move on to the second term in the right-hand side of (65). For any fixed $\gamma>0$, the definition of $\mathcal{E}_{1}^{i}$ gives rise to an upper estimate

For any constant $\gamma>0$, if $m/n\geq c_{0}\cdot\epsilon^{-2}\log\epsilon^{-1}$, then

with probability at least $1-C\exp(-c_{1}\epsilon^{2}m)$ for some universal constants $c_{0},c_{1},C>0$.

It remains to control the last term of (65). As mentioned above, the influence of this term is small since the set of $\boldsymbol{a}_{i}$’s satisfying $\mathcal{D}_{\gamma}^{i,2}$ accounts for a small fraction of measurements. Put formally, the number of equations satisfying $\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right|\geq\gamma\left\|\boldsymbol{h}\right\|$ decays rapidly for large $\gamma$ (at least at a quadratic rate), as stated below.

For any $0<\epsilon<1$, there exist some universal constants $c_{0},c_{1},C>0$ such that

with probability at least $1-C\exp\left(-c_{0}\epsilon^{2}m\right)$. This holds with the proviso $m/n\geq c_{1}\cdot\epsilon^{-2}\log\epsilon^{-1}$.

The constraint $\left|\frac{\boldsymbol{a}_{i}^{\top}\boldsymbol{h}}{\left\|\boldsymbol{h}\right\|}-\frac{2\boldsymbol{a}_{i}^{\top}\boldsymbol{z}}{\left\|\boldsymbol{h}\right\|}\right|\leq\gamma$ of $\mathcal{D}_{\gamma}^{i,2}$ necessarily requires

where the last inequality comes from our assumption on $\gamma$. With Lemma 6 in place, (72) immediately gives

In addition, on $\mathcal{E}_{1}^{i}\cap\mathcal{D}_{\gamma}^{i,2}$, the amplitude of each summand can be bounded in such a way that

where both inequalities are immediate consequences from the definitions of $\mathcal{D}_{\gamma}^{i,2}$ and $\mathcal{E}_{1}^{i}$ (see (56) and (24)). Taking this together with the cardinality bound (74) and picking $\epsilon$ appropriately, we get

one can simplify (77) by observing that $\vartheta_{1}\leq\frac{1}{100}$, which results in

Putting all preceding results in this subsection together reveals that with probability exceeding $1-\exp\left(-\Omega\left(m\right)\right)$,

holds simultaneously over all $\boldsymbol{x}$ and $\boldsymbol{z}$ satisfying

To conclude this section, we provide a tighter estimate about the norm of the regularized gradient.

Fix $\delta>0$, and assume that $y_{i}=(\boldsymbol{a}_{i}^{\top}\boldsymbol{x})^{2}$. Suppose that $m\geq c_{0}n$ for some large constant $c_{0}>0$. There exist some universal constants $c,C>0$ such that with probability at least $1-C\exp\left(-cm\right)$,

Lemma 7 complements the preceding arguments by allowing us to identify a concrete plausible range for the step size. Specifically, putting Lemma 7 and Proposition 2 together suggests that

Taking $\epsilon$ and $\delta$ to be sufficiently small we arrive at a feasible range (cf. Definition (39))

This establishes Proposition 1 and in turn Theorem 1 when $\mu_{t}$ is taken to be a fixed constant.

To justify the contraction under backtracking line search, it suffices to prove that the resulting step size falls within this range (83), which we defer to Appendix D.

Stability

This section goes in the direction of establishing stability guarantees of TWF. We concentrate on the iterative gradient stage, and defer the analysis for the initialization stage to Appendix C.

Before continuing, we collect two bounds that we shall use several times. The first is the observation that

where the last inequality follows from Cauchy-Schwarz. Setting

as usual, this inequality together with the trimming rules $\mathcal{E}_{1}^{i}$ and $\mathcal{E}_{2}^{i}$ gives

where (i) arises from [40, Corollary 5.35].

Unfortunately, (86) does not hold for all $\boldsymbol{z}$ within the neighborhood of $\boldsymbol{x}$ due to the existence of noise. Instead we establish the following:

The condition (86) holds for all $\boldsymbol{h}$ obeying

for some constants $c_{3},c_{4}>0$ (we shall call it Regime 1); this will be proved later. In this regime, the reasoning in Section 3 gives

for some appropriate constants $\mu,\rho>0$ and, hence, error contraction occurs as in the noiseless setting.

However, once the iterate enters Regime 2 where

for some constant $c_{5}>0$. Moreover, as long as $\|\boldsymbol{\eta}\|_{\infty}/\|\boldsymbol{x}\|^{2}$ is sufficiently small, one can guarantee that $c_{5}\frac{\|\boldsymbol{\eta}\|}{\sqrt{m}\|\boldsymbol{x}\|}\leq c_{5}\frac{\|\boldsymbol{\eta}\|_{\infty}}{\|\boldsymbol{x}\|}\leq c_{4}\|\boldsymbol{x}\|$. In other words, if the iterate jumps out of Regime 2, it will still fall within Regime 1.

Below we justify the condition (86) for Regime 1, for which we start by gathering additional properties of the trimming rules. By Cauchy-Schwarz, $\frac{1}{m}\left\|\boldsymbol{\eta}\right\|_{1}\leq\frac{1}{\sqrt{m}}\left\|\boldsymbol{\eta}\right\|\leq\frac{1}{c_{3}}\left\|\boldsymbol{h}\right\|\left\|\boldsymbol{z}\right\|$. When $c_{3}$ is sufficiently large, applying Lemmas 1 and 2 gives

We are now ready to analyze the regularized gradient, which we separate into several components as follows

For each index $i\in\mathcal{G}$, the inclusion property (51) (i.e. $\mathcal{E}_{3}^{i}\subseteq\mathcal{E}_{2}^{i}\subseteq\mathcal{E}_{4}^{i}$) holds. To see this, observe that

Next, letting $w_{i}=\frac{2\eta_{i}}{\boldsymbol{a}_{i}^{\top}\boldsymbol{z}}{\bf 1}_{\mathcal{E}_{1}^{i}\cap\mathcal{E}_{2}^{i}}{\bf 1}_{\{i\in\mathcal{G}\}}$, we see that for any constant $\delta>0$, the noise component obeys

provided that $m/n$ is sufficiently large. Here, (ii) arises from [40, Corollary 5.35], and the last inequality is a consequence of the upper estimate

Since $\|\boldsymbol{h}\|\geq c_{3}\frac{\|\boldsymbol{\eta}\|}{\sqrt{m}\|\boldsymbol{z}\|}$ for some large constant $c_{3}>0$, setting $\epsilon$ to be small one obtains

which finishes the proof of Theorem 2 for general $\boldsymbol{\eta}$.

Up until now, we have established the theorem for general $\boldsymbol{\eta}$, and it remains to specialize it to the Poisson model. Standard concentration results, which we omit, give

with high probability. Substitution into (16) completes the proof.

Minimax lower bound

The basic idea is to adopt the general reduction scheme discussed in [41, Section 2.2], which amounts to finding a finite collection of hypotheses that are minimally separated. Below we gather one result useful for constructing and analyzing such hypotheses.

for all $\boldsymbol{w}^{(l)},\boldsymbol{w}^{(j)}\in\mathcal{M}$,

In words, Lemma 8 constructs a set $\mathcal{M}$ of exponentially many vectors/hypotheses scattered around $\boldsymbol{x}$ and yet well separated. From (ii) we see that each pair of hypotheses in $\mathcal{M}$ is separated by a distance roughly on the order of $1$, and all hypotheses reside within a spherical ball centered at $\boldsymbol{x}$ of radius $3/2+o(1)$. When $\|\boldsymbol{x}\|\geq\log^{1.5}m$, every hypothesis $\boldsymbol{w}\in\mathcal{M}$ satisfies $\|\boldsymbol{w}\|\approx\|\boldsymbol{x}\|\gg 1$. In addition, (iii) says that the quantities ${|\boldsymbol{a}_{i}^{\top}\left(\boldsymbol{w}-\boldsymbol{x}\right)|}/{|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|}$ are all very well controlled (modulo some logarithmic factor). In particular, when $\|\boldsymbol{x}\|\geq\log^{1.5}m$, one must have

In the Poisson model, such a quantity turns out to be crucial in controlling the information divergence between two hypotheses, as demonstrated in the following lemma.

Lemma 9 and (106) taken collectively suggest that on the event $\mathcal{B}\cap\mathcal{C}$ ($\mathcal{B}$ is in Lemma 8 and $\mathcal{C}:=\{\|\boldsymbol{A}\|\leq\sqrt{2m}\}$), the conditional KL divergence (we condition on the $\boldsymbol{a}_{i}$’s) obeys

here, the inequality holds for some constant $c_{3}>0$ provided that $\left\|\boldsymbol{x}\right\|\geq\log^{1.5}m$, and the last inequality is a result of $\mathcal{C}$ (which occurs with high probability). We now use hypotheses as in Lemma 8 but rescaled in such a way that

for some $0<\delta<1$. This is achieved via the substitution $\boldsymbol{w}\longleftarrow\boldsymbol{x}+\delta(\boldsymbol{w}-\boldsymbol{x})$; with a slight abuse of notation, $\mathcal{M}$ denotes the new set.

The hardness of a minimax estimation problem is known to be dictated by information divergence inequalities such as (108). Indeed, suppose that

holds, then the Fano-type minimax lower bound [41, Theorem 2.7] asserts that

Since $M=\exp(n/30)$, (110) would follow from

Hence, we just need to select $\delta$ to be a small multiple of $\sqrt{n/m}$. This in turn gives

Discussion

To keep our treatment concise, this paper does not strive to explore all possible generalizations of the theory. There are nevertheless a few extensions worth pointing out.

More general objective functions. For concreteness, we restrict our analysis to the Poisson log-likelihood function, but the analysis framework we laid out easily carries over to a broad class of (nonconvex) objective functions. For instance, all results continue to hold if we replace the Poisson log-likelihood by the Gaussian log-likelihood; that is, the polynomial function $-\sum_{i=1}^{m}(y_{i}-|\boldsymbol{a}_{i}^{\top}\boldsymbol{z}|^{2})^{2}$ studied in . A general guideline is to first check whether the expected regularity condition

Sub-Gaussian measurements. The theory extends to the situation where the $\boldsymbol{a}_{i}$’s are i.i.d. sub-Gaussian random vectors, although the truncation threshold might need to be tweaked based on the sub-Gaussian norm of $\boldsymbol{a}_{i}$. A more challenging scenario, however, is the case where the $\boldsymbol{a}_{i}$’s are generated according to the CDP model, since there is much less randomness to exploit in the mathematical analysis. We leave this to future research.

Having demonstrated the power of TWF in recovering a rank-one matrix $\boldsymbol{x}\boldsymbol{x}^{*}$ from quadratic equations, we remark on the potential of TWF towards recovering low-rank matrices from rank-one measurements. Imagine that we wish to estimate a rank-$r$ matrix $\boldsymbol{X}\succeq{\bf 0}$ and that all we know about $\boldsymbol{X}$ is

It is known that this problem can be efficiently solved by using more computational-intensive semidefinite programs . With the hope of developing a linear-time algorithm, one might consider a modified TWF scheme, which would maintain a rank-$r$ matrix variable and operate as follows: perform truncated spectral initialization, and then successively update the current guess via a regularized gradient descent rule applied to a presumed log-likelihood function.

Moving away from i.i.d. sub-Gaussian measurements, there is a proliferation of problems that involve completion of a low-rank matrix $\boldsymbol{X}$ from partial entries, where the rank is known a priori. It is self-evident that such entry-wise observations can also be cast as rank-one measurements of $\boldsymbol{X}$. Therefore, the preceding modified TWF may add to recent literature in applying non-convex schemes for low-rank matrix completion , robust PCA , or even a broader family of latent-variable models (e.g. dictionary learning , sparse coding , and mixture problems ). A concrete application of this flavor is a simple form of the fundamental alignment/matching problem . Imagine a collection of $n$ instances, each representing an image of the same physical object but with different shift $r_{i}\in\{0,\cdots,M-1\}$. The goal is to align all these instances from observations on the relative shift between pairs of them. Denoting by $\boldsymbol{X}_{i}$ the cyclic shift by an amount $r_{i}$ of $\boldsymbol{I}_{M}$, one sees that the collection matrix $\boldsymbol{X}:=[\boldsymbol{X}_{i}^{\top}\boldsymbol{X}_{j}]_{1\leq i,j\leq k}$ is a rank-$M$ matrix, and the relative shift observations can be treated as rank-one measurements of $\boldsymbol{X}$. Running TWF over this problem instance might result in a statistically and computationally efficient solution. This would be of great practical interest.

Appendix A Proofs for Section 5

First, we make the observation that $(\boldsymbol{a}_{i}^{\top}\boldsymbol{z})^{2}-(\boldsymbol{a}_{i}^{\top}\boldsymbol{x})^{2}=\left(2\boldsymbol{a}_{i}^{\top}\boldsymbol{z}-\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right)\boldsymbol{a}_{i}^{\top}\boldsymbol{h}$ is a quadratic function in $\boldsymbol{a}_{i}^{\top}\boldsymbol{h}$. If we assume $\gamma\leq\frac{\alpha_{z}^{\text{lb}}\|\boldsymbol{z}\|}{\|\boldsymbol{h}\|}$, then on the event $\mathcal{E}_{1}^{i}$ one has

Solving the quadratic inequality that specifies $\mathcal{D}_{\gamma}^{i}$ gives

Suppose for the moment that $\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\geq 0$, then the preceding two intervals are respectively equivalent to

Assuming (115) and making use of the observations

Setting $\gamma_{1}:=\frac{\gamma}{1+\sqrt{2}}$ gives

Proceeding with the same argument, we can derive exactly the same inner and outer bounds in the regime where $\boldsymbol{a}_{i}^{\top}\boldsymbol{z}<0$, concluding the proof.

A.2 Proof of Lemma 4

By homogeneity, it suffices to establish the claim for the case where both $\boldsymbol{h}$ and $\boldsymbol{z}$ are unit vectors.

Suppose for the moment that $\boldsymbol{h}$ and $\boldsymbol{z}$ are statistically independent from $\{\boldsymbol{a}_{i}\}$. We introduce two auxiliary Lipschitz functions approximating indicator functions:

Since $\boldsymbol{h}$ and $\boldsymbol{z}$ are assumed to be unit vectors, these two functions obey

We proceed to lower bound $\frac{1}{m}\sum_{i=1}^{m}\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right)^{2}\chi_{z}\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\right)\chi_{h}\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right)$.

Firstly, to compute the mean of $(\boldsymbol{a}_{i}^{\top}\boldsymbol{h})^{2}\chi_{z}(\boldsymbol{a}_{i}^{\top}\boldsymbol{z})\chi_{h}(\boldsymbol{a}_{i}^{\top}\boldsymbol{h})$, we introduce an auxiliary orthonormal matrix

whose first row is along the direction of $\boldsymbol{z}$, and set

where the identity (122) arises from (66) and (67). Since $\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right)^{2}\chi_{z}\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{z}\right)\chi_{h}\left(\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right)$ is bounded in magnitude by $\gamma^{2}\left\|\boldsymbol{h}\right\|^{2}$, it is a sub-Gaussian random variable with sub-Gaussian norm $O(\gamma^{2}\left\|\boldsymbol{h}\right\|^{2})$. Apply the Hoeffding-type inequality [40, Proposition 5.10] to deduce that for any $\epsilon>0$,

with probability at least $1-\exp(-\Omega(\epsilon^{2}m))$.

The next step is to obtain uniform control over all unit vectors, for which we adopt a basic version of an $\epsilon$-net argument. Specifically, we construct an $\epsilon$-net $\mathcal{N}_{\epsilon}$ with cardinality $\left|\mathcal{N}_{\epsilon}\right|\leq\left(1+2/\epsilon\right)^{2n}$ (cf. ) such that for any $\left(\boldsymbol{h},\boldsymbol{z}\right)$ with $\left\|\boldsymbol{h}\right\|=\left\|\boldsymbol{z}\right\|=1$, there exists a pair $\boldsymbol{h}_{0},\boldsymbol{z}_{0}\in\mathcal{N}_{\epsilon}$ satisfying $\left\|\boldsymbol{h}-\boldsymbol{h}_{0}\right\|\leq\epsilon$ and $\left\|\boldsymbol{z}-\boldsymbol{z}_{0}\right\|\leq\epsilon$. Now that we have discretized the unit spheres using a finite set, taking the union bound gives

with probability at least $1-(1+2/\epsilon)^{2n}\exp(-\Omega(\epsilon^{2}m))$.

Define $f_{1}(\cdot)$ and $f_{2}(\cdot)$ such that $f_{1}(\tau):=\tau\chi_{h}(\sqrt{\tau})$ and $f_{2}(\tau):=\chi_{z}(\sqrt{\tau})$, which are both bounded functions with Lipschitz constant $O(1)$. This guarantees that for each unit vector pair $\boldsymbol{h}$ and $\boldsymbol{z}$,

Consequently, there exists some universal constant $c_{3}>0$ such that

where (i) results from Lemma 1, and (ii) arises from Lemma 2 whenever $\epsilon<1/2$.

With the assertion (125) in place, we see that with high probability,

for all unit vectors $\boldsymbol{h}$ and $\boldsymbol{z}$. Since $\epsilon$ can be arbitrary, putting this and (119) together completes the proof.

A.3 Proof of Lemma 5

The proof makes use of standard concentration of measure and covering arguments, and it suffices to restrict our attention to unit vectors $\boldsymbol{h}$. We find it convenient to work with an auxiliary function

Apparently, $\chi_{2}\left(\tau\right)$ is a Lipschitz function of $\tau$ with Lipschitz norm $O\left(\gamma\right)$. Recalling the definition of $\mathcal{D}_{\gamma}^{i,1}$, we see that each summand is bounded above by

For each fixed $\boldsymbol{h}$ and $\epsilon>0$, applying the Bernstein inequality [40, Proposition 5.16] gives

with probability exceeding $1-\exp\left(-\Omega\left(\epsilon^{2}m\right)\right)$.

From [40, Lemma 5.2], there exists an $\epsilon$-net $\mathcal{N}_{\epsilon}$ of the unit sphere with cardinality $\left|\mathcal{N}_{\epsilon}\right|\leq\left(1+\frac{2}{\epsilon}\right)^{n}$. For each $\boldsymbol{h}$, suppose that $\left\|\boldsymbol{h}_{0}-\boldsymbol{h}\right\|\leq\epsilon$ for some $\boldsymbol{h}_{0}\in\mathcal{N}_{\epsilon}$. The Lipschitz property of $\chi_{2}$ implies

where (i) arises by combining Lemmas 1 and 2. This demonstrates that with high probability,

for all unit vectors $\boldsymbol{h}$, as claimed.

A.4 Proof of Lemma 6

Without loss of generality, the proof focuses on the case where $\left\|\boldsymbol{h}\right\|=1$. Fix an arbitrary small constant $\delta>0$. One can eliminate the difficulty of handling the discontinuous indicator functions by working with the following auxiliary function

For any fixed unit vector $\boldsymbol{h}$, the above argument leads to an upper tail estimate: for any $0<t\leq 1$,

holds with probability exceeding $1-\exp\left(-\Omega(\epsilon^{2}m)\right)$.

We now proceed to obtain uniform control over all $\boldsymbol{h}$ and $2\leq\gamma\leq 2^{n}$. To begin with, we consider all $2\leq\gamma\leq m$ and construct an $\epsilon$-net $\mathcal{N}_{\epsilon}$ over the unit sphere such that: (i) $\left|\mathcal{N}_{\epsilon}\right|\leq\left(1+\frac{2}{\epsilon}\right)^{n}$; (ii) for any $\boldsymbol{h}$ with $\left\|\boldsymbol{h}\right\|=1$, there exists a unit vector $\boldsymbol{h}_{0}\in\mathcal{N}_{\epsilon}$ obeying $\left\|\boldsymbol{h}-\boldsymbol{h}_{0}\right\|\leq\epsilon$. Taking the union bound gives the following: with probability at least $1-\frac{\log m}{\log\left(1+\delta\right)}\left(1+\frac{2}{\epsilon}\right)^{n}\exp(-\Omega(\epsilon^{2}m))$,

holds simultaneously for all $\boldsymbol{h}_{0}\in\mathcal{N}_{\epsilon}$ and $\gamma_{0}\in\left\{\left(1+\delta\right)^{k}\mid 1\leq k\leq\frac{\log m}{\log\left(1+\delta\right)}\right\}$.

Putting the above results together gives that for all $2\leq\gamma\leq\left(1+\delta\right)^{\frac{\log m}{\log\left(1+\delta\right)}}=m$,

with probability exceeding $1-\frac{\log m}{\log\left(1+\delta\right)}\left(1+\frac{2}{\epsilon}\right)^{n}\exp\left(-c\epsilon^{2}m\right)$. This establishes (71) for all $2\leq\gamma\leq m$.

It remains to deal with the case where $\gamma>m$. To this end, we rely on the following observation:

where (i) comes from [6, Lemmas 3.1]. This basically tells us that with high probability, none of the indicator variables can be equal to 1. Consequently, $\frac{1}{m}\sum_{i=1}^{m}{\bf 1}_{\left\{\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right|\geq m\right\}}=0$, which proves the claim.

A.5 Proof of Lemma 7

Fix $\delta>0$. Recalling the notation v_{i}:=2\Big{\{}2\boldsymbol{a}_{i}^{\top}\boldsymbol{h}-\frac{\left|\boldsymbol{a}_{i}^{\top}\boldsymbol{h}\right|^{2}}{\boldsymbol{a}_{i}^{\top}\boldsymbol{z}}\Big{\}}{\bf 1}_{\mathcal{E}_{1}^{i}\cap\mathcal{E}_{2}^{i}}, we see from the expansion (62) that

as soon as $m\geq c_{1}n$ for some sufficiently large $c_{1}>0$. Here, the norm estimate $\left\|\boldsymbol{A}\right\|\leq\sqrt{m}\left(1+\delta\right)$ arises from standard random matrix results [40, Corollary 5.35].

Everything then comes down to controlling $\|\boldsymbol{v}\|$. To this end, making use of the inclusion (63) yields

The first term is controlled by [6, Lemma 3.1] in such a way that with probability $1-\exp(-\Omega(m))$,

Turning to the remaining terms, we see from the definition of $\mathcal{D}_{\gamma}^{i,1}$ and $\mathcal{D}_{\gamma}^{i,2}$ that

where the last inequality follows from (69) and (78).

Recall that $\gamma_{4}=3\alpha_{h}$. Taken together all these bounds lead to the upper bound

Appendix B Proofs for Section 7

Firstly, we collect a few results on the magnitudes of $\boldsymbol{a}_{i}^{\top}\boldsymbol{x}$ ($1\leq i\leq m$) that will be useful in constructing the hypotheses. Observe that for any given $\boldsymbol{x}$ and any sufficiently large $m$,

To summarize, with probability $1-o(1)$, one has

In the sequel, we will first produce a set $\mathcal{M}_{1}$ of exponentially many vectors surrounding $\boldsymbol{x}$ in such a way that every pair is separated by about the same distance, and then verify that a non-trivial fraction of $\mathcal{M}_{1}$ obeys (105). Without loss of generality, we assume that $\boldsymbol{x}$ takes the form $\boldsymbol{x}=\left[b,0,\cdots,0\right]^{\top}$ for some $b>0$.

The construction of $\mathcal{M}_{1}$ follows a standard random packing argument. Let $\boldsymbol{w}=\left[w_{1},\cdots,w_{n}\right]^{\top}$ be a random vector with

where $z_{i}\overset{\text{ind.}}{\sim}\mathcal{N}\left(0,1\right)$. The collection $\mathcal{M}_{1}$ is then obtained by generating $M_{1}=\exp\left(\frac{n}{20}\right)$ independent copies $\boldsymbol{w}^{(l)}$ ($1\leq l<M_{1}$) of $\boldsymbol{w}$. For any $\boldsymbol{w}^{(l)},\boldsymbol{w}^{(j)}\in\mathcal{M}_{1}$, the concentration inequality [40, Corollary 5.35] gives

Taking the union bound over all ${M_{1}\choose 2}$ pairs we obtain

with probability exceeding $1-2M_{1}^{2}\exp\left(-\frac{n}{8}\right)\geq 1-2\exp\left(-\frac{n}{40}\right)$.

The next step is to show that many vectors in $\mathcal{M}_{1}$ satisfy (105). For any given $\boldsymbol{w}$ with $\boldsymbol{r}:=\boldsymbol{w}-\boldsymbol{x}$, by letting $\boldsymbol{a}_{i,\perp}:=\left[a_{i,2},\cdots,a_{i,n}\right]^{\top}$, $r_{\|}:=r_{1}$, and $\boldsymbol{r}_{\perp}:=\left[r_{2},\cdots,r_{n}\right]^{\top}$, we derive

It then boils down to developing an upper bound on $\frac{\left|\boldsymbol{a}_{i,\perp}^{\top}\boldsymbol{r}_{\perp}\right|^{2}}{\left|a_{i,1}\right|^{2}}$. This ratio is convenient to work with since the numerator and denominator are stochastically independent. To simplify presentation, we reorder $\{\boldsymbol{a}_{i}\}$ in a way that

this will not affect our subsequent analysis concerning $\boldsymbol{a}_{i,\perp}^{\top}\boldsymbol{r}_{\perp}$ since it is independent of $\boldsymbol{a}_{i}^{\top}\boldsymbol{x}$.

To proceed, we let $\boldsymbol{r}_{\perp}^{(l)}$ consist of all but the first entry of $\boldsymbol{w}^{(l)}-\boldsymbol{x}$, and introduce the indicator variables

where $k=\frac{m}{4\log m}$ as before. In words, we divide $\boldsymbol{a}_{i,\perp}^{\top}\boldsymbol{r}_{\perp}^{(l)}$, $1\leq i\leq m$ into two groups, with the first group enforcing far more stringent control than the second group. These indicator variables are useful since any $\boldsymbol{w}^{(l)}$ obeying $\prod_{i=1}^{m}\xi_{i}^{l}=1$ will satisfy (105) when $n$ is sufficiently large. To see this, note that for the first group of indices, $\xi_{i}^{l}=1$ requires

where the second inequality follows from (134). This taken collectively with (130) and (135) yields

Regarding the second group of indices, $\xi_{i}^{l}=1$ gives

where the last inequality again follows from (134). Plugging (138) and (131) into (135) gives

Consequently, (105) is satified for all $1\leq i\leq m$. It then suffices to guarantee the existence of exponentially many vectors obeying $\prod_{i=1}^{m}\xi_{i}^{l}=1$.

Note that the first group of indicator variables are quite stringent, namely, for each $i$ only a fraction $O(1/m)$ of the equations could satisfy $\xi_{i}^{l}=1$. Fortunately, $M_{1}$ is exponentially large, and hence even ${M_{1}}/{m^{k}}$ is exponentially large. Put formally, we claim that the first group satisfies

with probability exceeding $1-\exp\left(-\Omega\left(k\right)\right)-\exp(-\widetilde{M}_{1}/4)$. With this claim in place (which will be proved later), one has

when $n$ and $m/n$ are sufficiently large. In light of this, we will let $\mathcal{M}_{2}$ be a collection comprising all $\boldsymbol{w}^{(l)}$ obeying $\prod_{i=1}^{k}\xi_{i}^{l}=1$, which has size $M_{2}\geq\frac{1}{2}\exp\left(\frac{1}{25}n\right)$ based on the preceding argument. For notational simplicity, it will be assumed that the vectors in $\mathcal{M}_{2}$ are exactly $\boldsymbol{w}^{(j)}$ ($1\leq j\leq M_{2}$).

We now move on to the second group by examining how many vectors $\boldsymbol{w}^{(j)}$ in $\mathcal{M}_{2}$ further satisfy $\prod_{i=k+1}^{m}\xi_{i}^{j}=1$. Notably, the above construction of $\mathcal{M}_{2}$ relies only on $\{\boldsymbol{a}_{i}\}_{1\leq i\leq k}$ and is independent of the remaining vectors $\{\boldsymbol{a}_{i}\}_{i>k}$. In what follows the argument proceeds conditional on $\mathcal{M}_{2}$ and $\{\boldsymbol{a}_{i}\}_{1\leq i\leq k}$. Applying the union bound gives

This combined with Markov’s inequality gives

with probability $1-o(1)$. Putting the above inequalities together suggests that with probability $1-o(1)$, there exist at least

vectors in $\mathcal{M}_{2}$ satisfying $\prod_{l=k+1}^{m}\xi_{i}^{l}=1$. We then choose $\mathcal{M}$ to be the set consisting of all these vectors, which forms a valid collection satisfying the properties of Lemma 8.

Finally, the only remaining step is to establish the claim (139). To start with, consider an $n\times k$ matrix $\boldsymbol{B}:=\left[\boldsymbol{b}_{1},\cdots,\boldsymbol{b}_{k}\right]$ of i.i.d. standard normal entries, and let $\boldsymbol{u}\sim\mathcal{N}\left({\bf 0},\frac{1}{n}\boldsymbol{I}_{n}\right)$. Conditional on the $\{\boldsymbol{b}_{i}$’s,

For sufficiently large $m$, one has $k=\frac{m}{4\log m}\leq\frac{1}{4}n$. Using [40, Corollary 5.35] we get

with probability $1-\exp\left(-\Omega(k)\right)$. Thus, for any constant $0<\epsilon<\frac{1}{2}$, conditional on $\{\boldsymbol{b}_{i}\}$ and (140) we obtain

When it comes to our quantity of interest, the above lower bound (143) indicates that on an event (defined via $\{\boldsymbol{a}_{i}\}$) of probability approaching 1, we have

Since conditional on $\{\boldsymbol{a}_{i}\}$, $\prod\nolimits_{i=1}^{k}\xi_{i}^{l}$ are independent across $l$, applying the Chernoff-type bound [55, Theorem 4.5] gives

with probability exceeding 1-\exp\Big{(}-\frac{1}{8}\frac{M_{1}}{\left(2\pi\right)^{k/2}(1+4\sqrt{k/n})^{k/2}}\left(\frac{1}{\sqrt{2\pi}m}\right)^{k}\Big{)}. This concludes the proof.

B.2 Proof of Lemma 9

Before proceeding, we introduce the $\chi^{2}$-divergence between two probability measures $P$ and $Q$ as

It is well known (e.g. [41, Lemma 2.7]) that

and hence it suffices to develop an upper bound on the $\chi^{2}$ divergence.

The preceding identity (147) arises from the following computation: by definition of $\chi^{2}(\cdot\|\cdot)$,

Set $\boldsymbol{r}:=\boldsymbol{w}_{1}-\boldsymbol{w}_{0}$. To summarize,

Appendix C Initialization via truncated spectral Method

This section demonstrates that the truncated spectral method works when $m\asymp n$, as stated in the proposition below.

for some sufficiently small constant $\varepsilon>0$. With probability exceeding $1-\exp\left(-\Omega\left(m\right)\right)$, the solution $\boldsymbol{z}^{(0)}$ returned by the truncated spectral method obeys

provided that $m>c_{0}n$ for some constant $c_{0}>0$.

By homogeneity, it suffices to consider the case where $\|\boldsymbol{x}\|=1$. Recall from [6, Lemma 3.1] that $\frac{1}{m}\sum_{i=1}^{m}|\boldsymbol{a}_{i}^{\top}\boldsymbol{x}|^{2}\in[1\pm\varepsilon]\|\boldsymbol{x}\|^{2}$ holds with probability $1-\exp(-\Omega(m))$. Under the hypothesis (150), one has