The local convexity of solving systems of quadratic equations

Introduction

However, due to the large dimensionality, storing all of the incoming vectors $x_{i}$ might be prohibitive. Instead, we randomly draw a set of sensing vectors $\{a_{k}\}$ which are efficient to store (e.g., they are sparse) and for each incoming data point compute $y_{ki}=(a_{k}^{T}x_{i})^{2}$. We are now only storing $(y_{ki},a_{k})_{k,i}$ which is a sparse data set. Note that if we define

The question posed above is: can we compute the covariance structure of the $\{x_{i}\}$ given only this data?

The example above describes covariance sketching of high-dimensional data streams [DSBN12, CCG13], but there are many other scenarios that fall under our problem setting, e.g., phaseless measurements in physics and optics [RBM94, TLOB12, Ger72, Fie82]. Because this data is invariant under the transformation

In the rank-1 setting in particular, several alternative reconstruction algorithms have been proposed with global phase recovery guarantees which operate directly on the lower-dimensional problem, and thus are more computationally efficient. Notably, [NJS13] considers the nonconvex optimization problem

and proves that after a judiciously chosen initialization, with high probability alternating minimization will converge to the underlying vector $x$ up to phase, assuming random Gaussian measurements. Subsequently [CLS14] used the same initialization to show convergence when followed by gradient descent without requiring resampling. Both of these algorithms provably recover the underlying vector $x$ up to global phase, from a number of measurements $m$ which is optimal up to additional logarithmic factors in $n$. Very recently, the paper [CC15] provides a modified gradient method which removes the additional logarithmic factors of $n$ in the number of measurements.

In a similar vein, many recent works have demonstrated global convergence guarantees for gradient descent on other nonconvex matrix factorization problems. Specifically, in [ZB15] the authors consider gradient descent on the Grassmannian and prove global convergence for a class of SVD problems. In [DSOR14] a stochastic gradient algorithm was shown to converge globally for a low-rank matrix least squares problem. In [SQW15] the authors consider the recovery of a full-rank matrix from sparse linear measurements via manifold optimization over the sphere. In all of these works including ours, the underlying idea is that the lack of convexity can be fixed by operating on an appropriate matrix manifold.

In this paper, we consider the more general version of problem (1) in which the underlying matrix is of rank $r$:

As noted in [CSV13], it seems unlikely that a deterministic RIP condition holds in this setting. In any case, our local convexity results are novel and might shed light on other nonconvex problems unrelated to matrix recovery.

for general $r$, we demonstrate that after $m\geq Cnr(\log n)^{2}$ Gaussian samples, in a quantifiable region the function (2) is strongly convex in directions perpendicular to the manifold of solutions

The size of this region is independent of both the ambient dimension and the rank

with an additional factor of $r^{5}$ samples, a simple spectral initialization will land within this region with high probability and thus standard gradient descent on (2) will linearly converge to a global minimizer

In the real-valued rank one setting, the strong convexity result we present actually holds in much more generality than the initialization result – for sub-gaussian measurements – and we believe this should be of independent interest; in particular, our results hold for Bernoulli measurements and Sparse Gaussian measurements. We note that in the rank-1 setting, recovery results from general sub-gaussian measurements were also provided in [KL15] using convex optimization for reconstruction, and a similar incoherence condition on the underlying $x$ was also required there.

While preparing this manuscript, we became aware of [[Sol14], p.250] which also certifies local convexity of the function (2), for the special case of Gaussian measurements in the rank one setting.

Our results can be viewed as exact recovery guarantees for a special case of a manifold-constrained least squares problem where the manifold is the set of rank $r$ positive semidefinite matrices. This algorithm was studied empirically in [FM15]. Many nonconvex problems of interest can be reformulated as a manifold-constrained least squares problem, and we believe that the exact recovery guarantees presented here should be extendable to a broader class of problems.

Main results

by solving the nonconvex optimization problem

Because the function appearing in (4) is invariant under right multiplication by an orthogonal matrix, there is an entire manifold of solutions given by $\left\{XO:O\in\mathcal{O}(r)\right\}$ where $\mathcal{O}(r)$ is the set of $r\times r$ orthogonal matrices. Our strategy is to establish that a spectral initialization will land (with high probability) in a region of strong convexityStrong convexity here and throughout always refers to convexity in directions orthogonal to the manifold of solutions. around the manifold of global minimizers. An overview of our approach is given in Algorithm 1.

There are two main ingredients to proving performance guarantees for Algorithm 1, namely, the strong convexity of the function $f$ in a region around the manifold of global minimizers at finite sample complexity, and a guarantee that spectral initialization will land within this region. The finite sample convexity result holds in more generality when $r=1$, while for general $r$ we always assume Gaussian measurements.

for i.i.d. standard Gaussian vectors $\{a_{i}\}_{i=1}^{m}$. Now that we have an entire manifold of solutions given by $\left\{XO:O\in\mathcal{O}(r)\right\}$ we will need to consider the quantity

which is well-defined by compactness of the orthogonal group. We note that the minimizer may not be unique, but this is not important for our purposes. We will also need to consider

The main finite sample convexity result is as follows:

Then with probability at least $1-4e^{-rn}-7/m^{2}$, it holds that

We now show how this local strong convexity results in linear convergence to the true $X$ we seek to recover. We have the following theorem which concisely establishes the initialization and performance guarantees of Algorithm 1.

Suppose we take $m\geq C\|X\|_{F}^{8}\lambda_{r}^{-4}nr^{2}(\log n)^{2}$ samples of the form (5), where $\lambda_{1}$ and $\lambda_{r}$ are as in (7). Define the matrix

where $\sigma_{1}\geq\sigma_{2}...\geq\sigma_{r+1}>0$ are the eigenvalues of $M$ and $u_{i}$ are the corresponding normalized eigenvectors. If we iteratively update $U_{k}$ via gradient descent

then with probability at least $1-3e^{-rn}-7/m^{2}$,

For a proof of Theorem 2.2, see Section 3.3.

Note that the quantity $\|X\|_{F}^{8}\lambda_{r}^{-4}$ is scale invariant; however, we have the bounds

One consequence of our result is that the sampling complexity is entirely independent of the desired solution tolerance. That is, the fixed set of $m\geq Cnr^{6}(\log n)^{2}$ samples suffices to produce a global solution up to arbitrary accuracy.

Our numerical results in §4 suggest that in general the sampling complexity only linearly depends on the ambient dimension $n$. Consequently a more refined analysis and initialization procedure such as that found in the recent work of [CC15] for the case of rank-1 recovery is most likely possible also in the general rank-r recovery setting.

This method of analysis should find use in providing recovery guarantees by gradient descent for a broader class of nonconvex problems arising in machine learning applications such as matrix completion, nonnegative matrix factorization, clustering, etc. More generally, such an analysis could possibly be useful towards achieving provable guarantees for machine learning problems which have many unstable saddle points, such as neural networks [DPG+14].

2 Rank-One Matrix Recovery

where $\Sigma$ is the covariance matrix, which we assume is invertible. With this setup, we then consider minimization of the random function

If $m\geq C\|\Sigma\|_{op}^{2}n(\log n)^{3},$ then with probability greater than $1-4/n^{2},$

Above, $C>0$ is a constant which depends only on the sub-gaussian norm of $a_{i}$.

The finite sample convexity result holds for general sub-gaussian measurements satisfying (22), while our initialization results require more restrictive conditions, namely that the fourth moment of the measurements is close to that of Gaussian measurements; for simplicity we have only included the result for Gaussians which follows from Lemma 3.11 in the next section.

The rest of the paper is organized as follows: in §3.1 we prove the main finite sample convexity result Theorem 2.1, which relies on classifying tangent and normal directions to the manifold of solutions $\left\{XO:O\in\mathcal{O}(r)\right\}$ and an explicit formula for the expected Hessian. In §3.2 we prove convexity results for the rank one case under more general randomness assumptions. In §3.3 we prove that with high probability the initialization step produces a matrix in a convex region around the manifold of solutions and establish the convergence of gradient descent. Briefly in §3.4 we describe how our results generalize to the complex setting. Finally, in §4 we conclude with some numerical experiments demonstrating the performance and robustness of the results presented here.

Convexity

Here we present lemmas that are used in the proof of Theorem 2.1, as well as a summary of the proof. For the full proof, we refer the reader to Section 5.1.3.

The main lemma we rely on is the following simple characterization of the normal directions to the manifold of solutions:

Assume $X$ has full column rank and let $\displaystyle O^{*}=\arg\min_{O\in\mathcal{O}(r)}\|XO-U\|_{F}^{2}$, which is not necessarily unique. Then we can write

This basically follows from the solution to the Orthogonal Procrustes Problem [Sch66]. If we write $X^{T}U=ZDV^{T}$ for the singular value decomposition of $X^{T}U$, then we can expand the objective as follows:

is a symmetric positive semidefinite matrix. As $X^{T}A=0$ is equivalent to $P_{\perp}A=A$, we arrive at the stated claim. ∎

This lemma says that if we consider the direction $W=U-XO^{*}$ between $U$ and its closest solution matrix $XO^{*}$ we have that

which is a symmetric matrix. Why symmetry is important will become apparent after the next lemma, which establishes formulas for the expectation of the Hessian of (4):

The gradient of $f(U)=\frac{1}{4m}\sum_{i=1}^{m}(y_{i}-a_{i}^{T}UU^{T}a_{i})^{2}$ is given by

where the $nr\times nr$ block matrices $A$ and $D$ satisfy

For details, see §5.1.1 in the Appendix. We will also need a standard concentration result:

Suppose we collect $m\geq C\delta^{-2}\beta nr\log(nr)$ samples of the form $y_{i}:=a_{i}^{T}XX^{T}a_{i}$, where $\delta$ and $\beta$ are given constants and $r=$ rank$(X)$; then we have that with probability greater than $1-2e^{-\beta rn}-6/m^{2}$

The sampling complexity can be improved, but we state Lemma 3.3 as a general proof-of-concept. For details see §5.1.2.

To complete the proof sketch, observe that

can be written as a convex quadratic polynomial in $\|U-XO^{*}\|_{F}$, where the constant term is given by

and consequently we can bound its smallest positive root using the remarks above (see §3.2.2 for the rank one setting, where this observation is more straightforward). We apply the concentration from above along with the following one-sided martingale bound from [Ben03] (as stated in [CLS14]) to establish the stated non-asymptotic bound. For details see §5.1.3.

where one can take $c_{0}=25$ and $\Phi(\cdot)$ is the CDF for the standard normal.

2 Rank One

where $\Sigma$ is the covariance matrix, which we assume is invertible. Consider the eigenvalue decomposition of the covariance matrix, $\Sigma=\sum_{k=1}^{n}v_{k}v_{k}^{T}$. An important quantity in our analysis will be

a coherence parameter for $\Sigma^{-1/2}x,$ and

We consider convexity of the function $f(u)$ defined in (12) (equivalently, positive semi-definiteness of the Hessian matrix $\nabla^{2}f(u)$) in the neighborhood of $u=x$, first in expectation with respect to the draw of $a_{i}$, or in the limit of infinitely many samples $m$. These results are necessary for the proof of Lemma 3.9.

For details, see §5.3.1 in the Appendix. We then have the following asymptotic convexity result:

where $\displaystyle\tau(x)$ is the coherence of $\Sigma^{1/2}x$ as in (23), and above $[u]_{-}=\min\{u,0\}$ and $[u]_{+}=\max\{u,0\}$.

for some $\delta\leq 1$. In fact we find that a loose bound is given by

This result alone provides enough information to prove performance guarantees for stochastic gradient descent after an initialization procedure. Via a union bound and covering argument, this result along with matrix concentration will also guarantee uniform convexity in this region at finite sample size $m\geq n^{2}$. However, to ensure uniform convexity at finite sample size $m\geq Cn\log(n)$, we will need a more refined analysis based on the structure of the Hessian matrix, as presented in the next section.

2.2 Non-Asymptotic Convexity

Here we present the sketch of the proof of Theorem 2.3. For the full proof, we refer the reader to Section 5.3.5.

As before, we will use the standard concentration result:

This result can be proved by first truncating the norms of the measurements vectors and then applying Matrix Bernstein’s Inequality (e.g., [Tro12]). The sampling complexity can be improved, but we state Lemma 3.8 as a general proof-of-concept. For details see §5.3.4. For $\epsilon$ sufficiently small, this result indicates that we can control the eigenvalues of $\nabla^{2}f(u)$ for $u$ sufficiently close to $x$. In particular, if $\nabla^{2}f(u)$ is positive definite in a region around $x$, then $f(u)$ is strongly convex and $x$ is the unique minimum in this region. It is not immediately clear how to extend such control to a quantifiable region around $x$. However, Theorem 2.3 requires only that we have a lower bound on the eigenvalues.

Assuming that $\Sigma=Id$, the same technique from §3.1 can be applied: first write $u=x+t\hat{w}$ for a unit vector $\|\hat{w}\|_{2}=1$ and observe that

and consequently using Lemma 3.8 and Lemma 3.6 we can control this term. As before, applying Lemma 3.4 to the positive term

where $\displaystyle\tau(x)=\max_{1\leq k\leq n}\frac{|x_{k}|^{2}}{\|x\|_{2}^{2}}$ is the coherence of $x$.

The proof of this uses the fact that the smallest eigenvalue is a concave function of $\mu_{4}$; the proof can be found in §5.3.2. We can now quantify the lower bound appearing in (15) for a large class of sub-gaussian measurements:

Bernoulli: For standard Bernoulli measurement vectors, where $a_{ik}$ are i.i.d. $\pm 1$ with equal probability, $\mu_{4}=1$ and we have a quantifiable strong convexity guarantee so long as $x$ is incoherent, i.e., $\tau(x)<1/2$. This is sharp in the sense that for $x=\begin{bmatrix}1/\sqrt{2}&1/\sqrt{2}\end{bmatrix}$ the expected Hessian has a 0 eigenvalue.

Gaussian: For vectors $a_{i}$ with i.i.d. standard Gaussian entries, $\mu_{4}=3$ and Lemma 3.9 provides the uniform lower bound

for all $\|u-x\|_{2}\leq\frac{1}{15}\|x\|_{2}$.

Sparse Gaussian: Note that (28) holds anytime $\mu_{4}\geq 3$ by Lemma 3.9. This includes sparse Gaussian vectors, whose coordinates are i.i.d. standard normal with probability $p$ and 0 with probability $1-p$. In this case

3 Initialization and Gradient Descent

We have shown that the function is strongly convex in a quantifiable region around the global minimizers. To guarantee results for gradient descent, we will also need the following lemma which bounds the Lipschitz constant of the gradient of our function.

Consider the function $f(U)=\frac{1}{4m}\sum_{i=1}^{m}(y_{i}-\|a_{i}^{T}U\|_{2}^{2})^{2}$. Suppose $m\geq C$. For a universal constant $C>0$, it holds with probability exceeding $1-2m^{-3}$ that for any $U$ within the region of convexity given by (8),

with $B=Cn^{2}\log(m)^{2}\lambda_{r}$. Here, $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{r}$ are the eigenvalues of $XX^{T}$ and $\|U-XO^{*}\|^{2}_{F}=\min_{O\in\mathcal{O}(r)}\|XO-U\|^{2}_{F}.$

By the sub-gaussian assumption, the following holds with probability exceeding $1-2m^{-3}$:

Conditioning on this event, recalling that $\|X\|_{F}^{2}=\lambda_{1}+\dots+\lambda_{r}$, and recalling the formula for the gradient $\nabla f$ in (18), observe the bound

Thus, $\|\nabla f(U)-\nabla f(X)\|_{F}\leq C\lambda_{r}n^{2}\log(m)^{2}\|U-XO^{*}\|_{F}.$ ∎

It remains to certify a point in this region to initialize gradient descent.

Suppose we take $m\geq C\beta\lambda_{r}^{-4}\|X\|_{F}^{8}nr^{2}(\log n)^{2}$ samples of the form (5), where $\lambda_{1}$ and $\lambda_{r}$ are as in (7). Define the matrix

where $\sigma_{1}\geq\sigma_{2}...\geq\sigma_{r+1}>0$ are the eigenvalues of $M$ and $u_{i}$ are the corresponding normalized eigenvectors. Then with probability at least $1-3e^{-\beta rn}-7/m^{2}$ we have that

For the proof, see §5.2. Initializing from a matrix satisfying (30) guarantees we are close enough so that gradient descent will converge. We can now prove the main theorem, Theorem 2.2:

Given the number of samples $m\geq C\beta\lambda_{r}^{-4}\|X\|_{F}^{8}nr(\log n),$ the following events simultaneously occur with the stated probability:

Lemma 3.11 holds, and thus $d(U_{0}):=\min_{O\in\mathcal{O}(r)}\|XO-U_{0}\|^{2}_{F}<\frac{9}{100\|X\|_{F}^{2}}\lambda_{r}^{2}$.

Theorem 2.1 holds, and so considering the Taylor expansion of $f$ around $U$, the following holds for all $U$ satisfying $d(U)<\frac{9}{100\|X\|_{F}^{2}}\lambda_{r}^{2}:$

Lemma 3.10 holds with Lipschitz constant $B=Cn^{2}\log(nr^{5})^{2}\lambda_{r}$.

Let $U^{+}:=U_{0}-\gamma\nabla f(U_{0})$ and $O^{*}=\text{arg}\min_{O\in\mathcal{O}(r)}\|XO-U_{0}\|_{F}^{2}$. Then

4 The Complex Case

where we note $\begin{bmatrix}a_{j}\\ b_{j}\end{bmatrix}\sim\mathcal{N}(0,1/2Id_{2n\times 2n})$. Moreover, the columns of $Z$ are orthogonal, and the map

gives us an isomorphism onto the orientation preserving component of the orthogonal group $\mathcal{O}(2)$. Thus this problem is equivalent to recovering an unknown real-valued rank 2 matrix, and the results in the previous sections reproduce known optimality guarantees for gradient descent in the phase retrieval model as found in e.g., [CLS14]. Specifically, we have shown the following:

Given $m\geq Cn(\log n)^{2}$ noiseless samples of the form

and let $\hat{Z}$ be the output of Algorithm 1 applied to the data $\left\{(y_{i},\mathbf{a}_{i})\right\}_{i=1}^{m}$ with constant step size

Then with probability at least $1-3e^{-\beta n}-7/m^{2}$ we have that

where $k$ is the number of iterations of gradient descent and $Z$ is the matrix (33).

Examples and Experiments

First, we consider the performance of the algorithm (1) in the rank-1 real-valued setting, where the measurements are $y_{i}=(a_{i}^{T}x)^{2}$. Our numerical studies strongly suggest that the algorithm (1) is stable to noise, that is, given measurements of the form $y_{i}=(a_{i}^{T}x)^{2}+\eta_{i},$ the algorithm successfully returns an matrix $\widehat{x}$ up to the noise level $\frac{\|\widehat{x}-x\|_{2}}{\|x\|_{2}}\leq\|\eta\|_{2}$. We consider three different measurement ensembles:

Bernoulli: $a_{i}$ are i.i.d. Bernoulli random vectors

Standard Gaussian: $a_{i}$ are i.i.d. drawn from ${\cal N}(0,Id_{n\times n})$.

Gaussian with covariance: $a_{i}$ are i.i.d drawn from ${\cal N}(0,\Sigma)$ with covariance matrix

In a first experiment, we fix an $n$-dimensional vector $x$ of unit norm with randomly-generated coefficients, and consider noiseless measurements $y_{i}=(a_{i}^{T}x)^{2}$. We implement the meta-algorithm 1, calling Matlab’s built-in function fminunc to find a stationary point starting from the initialization. In the local optimization procedure, we do not provide any information to fminunc other than the function itself; by default Matlab uses a quasi-Newton method for local minimization. We run this experiment using the three different measurement ensembles above, at problem size $n=100$ and at a number of measurements $m=2n,3n,\dots,8n$. If the solution $\widehat{x}$ recovered by the algorithm is within the tolerance $\min\{\|\widehat{x}-x\|_{2},\|\widehat{x}+x\|_{2}\}\leq.001$, we say the algorithm has succeeded in finding the global solution. In Figure 1, the results of this experiment are displayed, averaged over 100 trials.

Next, we analyze numerically the stability of the algorithm to additive measurement noise. For these experiments, we consider noisy measurements of the form

where $\eta_{i}$ are i.i.d. mean-zero uniformly distributed, and normalized such that $\|\eta\|_{2}=\mu\|\sum_{i}(a_{i}^{T}x)^{2}\|_{2}$ for $\mu=.5$ (low signal to noise ratio) and $\mu=2$ (high signal to noise ratio). We observe that the meta-algorithm is robust to such additive noise, with relative reconstruction error $\min\{\|\widehat{x}-x\|_{2},\|\widehat{x}+x\|_{2}\}$ averaging below the signal to noise threshold. We leave a theoretical analysis of this observed noise stability to future work.

where $Z\Sigma V^{T}=X^{T}\widehat{X}$ is the singular value decomposition.n In Figure 1, the results of the experiment are displayed, averaged over 100 trials.

R. Ward and C. White were funded in part by an NSF CAREER Grant and an AFOSR Young Investigator Award. We would like to thank Ju Sun for pointing out a mistake in the original rank one proof, and thank Mahdi Soltanolkotabi and Laurent Jacques for additional helpful comments and corrections. C. White would like to thank Aaron Royer and Ravi Srinivasan for helpful conversations.

References

Appendix

which can be seen by writing out the entries of the matrix individually. This implies

1.2 Proof of Lemma 3.3

We begin with a more general concentration result.

First Term. Note that we have a product of independent subexponential random variables, and so if we condition on the bounds

both of which happen with probability at least $1-1/m^{2}$, we find via Bernstein that

which we can make smaller than $e^{-\alpha rn}$ so long as $m\geq C\delta^{-2}\alpha nr\log(n)$ and we conclude via an $\epsilon$-net argument that

We can make this bound smaller than $e^{-\alpha rn}$ so long as $m\geq C\delta^{-2}\alpha rn$. As (35) happens with probability greater than $1-1/m^{2}$ we find

Second Term. For this term we further decompose $w$ into its $x$-component and its $x_{\perp}$-component. We can apply the same analysis for the first and third terms to the $x_{\perp}$ terms, and have only to deal with

which we can make smaller than $e^{-\alpha r}$ so long as $m\geq C\alpha\delta^{-2}r(\log r)^{2}$. Moreover, note that we can also make (37) smaller than $1/m^{2}$ for $m\geq C\delta^{-2}$. We conclude that

Combining (34), (38), and (36) yields the stated result. ∎

Suppose we collect $m\geq C\delta^{-2}\beta nr\log(n)^{2}$ samples of the form $y_{i}:=a_{i}^{T}XX^{T}a_{i}$, where $\delta$ and $\beta$ are given constants and $r=$ rank$(X)$; then we have that with probability greater than $1-2e^{-\beta rn}-6/m^{2}$

Note that $y_{i}=\sum_{k=1}^{r}(a_{i}^{T}x_{k})^{2}$, and so by Theorem 5.1 we find that with probability greater than $1-2e^{-\beta rn}-6/m^{2}$

Suppose we collect $m\geq C\delta^{-2}\beta nr\log(n)^{2}$ samples of the form $y_{i}:=a_{i}^{T}XX^{T}a_{i}$, where $\delta$ and $\beta$ are given constants and $r=$ rank$(X)$; then we have that with probability greater than $1-2e^{-\beta rn}-6/m^{2}$

1.3 Proof of the Convexity Theorem 2.1

We will rely on the following Lemma from [Ben03], as stated in [CLS14]:

where one can take $c_{0}=25$ and $\Phi(\cdot)$ is the CDF for the standard normal.

Let $\hat{W}:=W/\|W\|_{F}$ be the normalized direction from $U$ to $XO^{*}$ and let $t\geq 0$ be a positive scalar. Moreover, WLOG we will be assuming that $\|X\|_{F}=1$.

Finally, because $f(U)$ is invariant under the action of $\mathcal{O}(r),$ it suffices to consider the case where $O^{*}=Id$.

Consider the single-variable function $f(X+t\hat{W})$ which can be written

which is a convex polynomial in $t$; observe that $f^{\prime\prime}(0)>0$. If the linear term is positive, then clearly (40) is positive for all $t$ and we have nothing to show (the smallest eigenvalue is bounded below by $f^{\prime\prime}(0)$ in the direction $\hat{W}$). Define the following quantities:

Observe that $A_{i}$ is a chi-squared random variable with 1 degree of freedom by the normalization $\|\hat{W}\|_{F}=1$. Thus we have

By the definition of $f^{\prime\prime}(t)$, we have

Next we consider the variance of $Z_{i}(t)$:

we find that if $m\geq 288\alpha\lambda_{r}^{-2}C(t)^{2}nr$ then with probability at least $1-e^{-\alpha nr}$ we have

where we used (44) to lower bound $\mu(t)$ by

Moreover, an $\epsilon$-net argument over all directions $\hat{W}$ shows that (46) holds for an arbitrary $\hat{W}$ with probability at least $1-e^{-\beta nr}$.

Further observe that our condition on $m$ guarantees

with probability at least $1-2e^{-\beta rn}-6/m^{2}$ by Corollary 5.3. This implies that

with probability at least $1-3e^{-\beta rn}-6/m^{2}$ for any direction $\hat{W}$. Thus by a tangent line bound we find that the smallest positive root of $f^{\prime\prime}(t)$ is bounded below by

Thus $f^{\prime\prime}(t)\geq\frac{\lambda_{r}^{2}}{18}$ for all $t\in[0,\frac{3}{10}\lambda_{r}]$ which is the advertised lower bound.

For the upper bound, observe that by Cauchy-Schwarz

and thus we find an upper bound for (40) is given by

Now, $\omega_{i}:=\left(a_{i}^{T}\hat{W}\hat{W}^{T}a_{i}\right)$ is a chi-squared random variable with one degree of freedom; consequently we find

with probability greater than $1-e^{-\beta nr}$.

for all $t\in[0,\frac{3}{10}\lambda_{r}]$.

To finish the proof, observe that we can write

and by Corollary 5.3 (where, given the number of measurements $m$, we may take $\delta\geq C\lambda_{r}/\lambda_{1}$) we have

From everything above, we conclude that for $\frac{\|U\|_{F}}{\|X\|_{F}}\leq\frac{3}{10}\lambda_{r}\left(\frac{XX^{T}}{\|X\|_{F}^{2}}\right)$,

we conclude that for general $X$, it holds for $\|U\|_{F}\leq\frac{3}{10\|X\|_{F}}\lambda_{r}$ that

2 Proofs for 3.3, Initialization and Convergence

It suffices to prove the case $\|X\|_{F}^{2}=1$. By Corollary 5.2 we have that with probability greater than $1-2e^{-\beta rn}-6/m^{2}$

so long as $m\geq C\delta^{-2}\beta nr\log(n)^{2}$. This implies

Let $U_{0}:=U\Sigma^{1/2}$, and $Q:=O_{1}O_{2}^{T}$ where $(X^{T}X)^{-1/2}X^{T}U=O_{1}DO_{2}^{T}$ is the singular value decomposition and observe

where (53c) follows from Theorem 2 in [YWS15] and (53d) follows from

The first equality holds because $X$ has orthogonal columns and thus $X^{T}X$ is a diagonal matrix.

3 Proofs for Rank-one Matrix Recovery, §3.2

If $\Sigma\neq Id$, then observe that if we define $b:=\Sigma^{-1/2}a$,

then the inner term satisfies the assumptions needs for (56), and so we find

where $\Sigma^{1/2}=\begin{bmatrix}v_{1}&v_{2}&...&v_{n}\end{bmatrix}_{n\times n}$. ∎

3.2 Proof of Lemma 3.9

Suppose first that $u_{n}^{2}>0$. Then we can use the determinant formula

If any $u_{j}^{2}=u_{j+1}^{2}$ then $\lambda=-u_{j}^{2}$ is an eigenvalue.

If all of the squared coordinates are distinct then each eigenvalue $\lambda$ satisfies

and because there will be a vertical asymptote at each $-u_{k}^{2}$ we see the eigenvalues of $Z$ interlace the squared coordinates, the smallest occurring somewhere between $(-u_{1}^{2},-u_{2}^{2})$ and the largest somewhere after $-u^{2}_{n}$.

In general, if some of the coordinates are , we see that with the assumed ordering $Z$ will be a block matrix and we can apply (57) to the reduced space where $Z$ acts nontrivially.

The bound $\lambda_{min}\geq-1/2$ follows from

Note that by the Gershgorin Circle Theorem, the largest eigenvalue $\lambda_{max}$ is no larger than $1-u_{n}^{2}$. ∎

Lemma 5.5 above shows $g(1)\geq 1-\min\{2\tau(x),1\}$ and it is clear that $g(3)=1$. By concavity of $\lambda_{min}(\cdot)$ we then find

If $\mu_{4}>3$, then $2xx^{T}+(\mu_{4}-3)\sum_{k=1}^{n}x_{k}^{2}e_{k}e_{k}^{T}$ is a positive semi-definite matrix and thus $g(\mu)\geq 1$. Observing that

3.3 Proof of Lemma 3.6

For (60), we used Lemma 3.9. Lastly, note that

Consequently we can define the polynomials

and by convexity we can bound the smallest positive root by the intercept of the tangent line; the bounds (60) and (61) thus yield the stated conclusion for $\Sigma=Id$.

For general covariance matrices, note that we have just shown that

whenever $\Sigma^{1/2}u$ and $\Sigma^{1/2}x$ are close enough, which implies

3.4 Proof of Lemma 3.8

Begin by assuming $\Sigma=Id$ and $\|x\|_{2}=1$. Note that because of the sub-gaussian assumption we have that for $m\geq C$

where the constants depend on the sub-gaussian norm of $a_{i}$. Consequently

Bernstein’s inequality (Theorem 4.1 in [Tro12]) tells us that

where $C$ is a constant which depends on the moments of $a_{i}$.

Lastly observe that if $\hat{a}_{i}:=a_{i}\chi_{(a_{i}^{T}x)^{2}\|a_{i}\|_{2}^{2}\geq c^{2}n(\log m)^{2}}$ then we can write

where we used Jensen’s inequality for the first line and have assumed $m\geq Cn$. Consequently we find that for $m\geq Cn$

where $\delta:=\epsilon-3/m^{2}$ from (63). Now, all we have left is to show that the exponential can be made less than a power of $m$. Using (62) we find that we need $m$ to satisfy

for which it suffices to require $m\geq C\epsilon^{-2}n(\log n)^{3}$ for some constant $C$ which only depends on the moments of $a_{i}$.

For the more general statement note that our previous work shows

whenever $m\geq C\epsilon^{-2}\|\Sigma\|_{op}^{2}n(\log n)^{3}$. Consequently,

3.5 Proof of Theorem 2.3

Assume without loss that $\|x\|_{2}=1$ and that $\hat{w}$ is the normalized direction from $u$ to $x$. Moreover begin by assuming $\Sigma=Id$. Closely following the proof of Theorem 2.1 we first note that

Note that by the computations done in Lemma 3.5 we have

where $C(t)$ depends only on the subgaussian norm of the $a_{i}$.

Now, for a given $\epsilon\in(0,1)$ define

we find that if $m\geq 288\alpha\lambda^{-2}C(t)^{2}n$ then with probability at least $1-e^{-\alpha n}$ we have

where we used the concentration guaranteed by Lemma 3.8 above with $\epsilon<\lambda/8$ and the fact that

Consequently, using a tangent line bound for the smallest positive root we find that for all

Moreover, an $\epsilon$-net argument over all directions $\hat{w}$ shows that (67) holds for an arbitrary $\hat{w}$ with probability at least $1-e^{-\beta n}$.

For general covariance matrices, apply the previous argument to $\Sigma^{1/2}u$ and $\Sigma^{1/2}x$ with measurements $b_{i}=\Sigma^{-1/2}a_{i}$ as usual.