Global Convergence of Stochastic Gradient Descent for Some Non-convex Matrix Problems

Introduction

We analyze an algorithm to solve the stochastic optimization problem

In practice, many people use stochastic gradient descent (SGD) to solve (2). Efficient SGD implementations can scale to very large datasets . However, standard stochastic gradient descent on (2) does not converge globally, in the sense that there will always be some initial values for which the norm of the iterate will diverge (see Appendix A).

People have attempted to compensate for this with sophisticated methods like geodesic step rules and manifold projections ; however, even these methods cannot guarantee global convergence. Motivated by this, we describe Alecton, an algorithm for solving (2), and analyze its convergence. Alecton is an SGD-like algorithm that has a simple update rule with a step size that is a simple function of the norm of the iterate $Y_{k}$. We show that Alecton converges globally. We make the following contributions:

We establish the convergence rate to a global optimum of Alecton using a random initialization; in contrast, prior analyses have required more expensive initialization methods, such as the singular value decomposition of an empirical average of the data.

In contrast to previous work that uses bounds on the magnitude of the noise , our analysis depends only on the variance of the samples. As a result, we are able to be robust to different noise models, and we apply our technique to these problems, which did not previously have global convergence rates:

matrix completion, in which we observe entries of $A$ one at a time (Section 4.1),

subspace tracking, in which $A$ is a projection matrix and we observe random entries of a random vector in its column space (Section 4.4).

Our result is also robust to different noise models.

We describe a martingale-based analysis technique that is novel in the space of non-convex optimization. We are able to generalize this technique to some simple regularized problems, and we are optimistic that it has more applications.

Much related work exists in the space of solving low-rank factorized optimization problems. Foundational work in this space was done by Burer and Monteiro , who analyzed the low-rank factorization of general semidefinite programs. Their results focus on the classification of the local minima of such problems, and on conditions under which no non-global minima exist. They do not analyze the convergence rate of SGD.

Another general analysis in Journée et al. exhibits a second-order algorithm that converges to a local solution. Their results use manifold optimization techniques to optimize over the manifold of low-rank matrices. These approaches have attempted to correct for falling off the manifold using Riemannian retractions , geodesic steps , or projections back onto the manifold. General non-convex manifold optimization techniques tell us that first-order methods, such as SGD, will converge to a fixed point, but they provide no convergence rate to the global optimum. Our algorithm only involves a simple rescaling, and we are able to provide global convergence results.

Our work follows others who have studied individual problems that we consider. Jain et al. study matrix completion and provides a convergence rate for an exact recovery algorithm, alternating minimization. Candès et al. provide a similar result for phase retrieval. In contrast to these results, which require expensive SVD-like operations to initialize, our results allow random initialization. Our provided convergence rates apply to additional problems and SGD algorithms that are used in practice (but are not covered by previous analysis). However, our convergence rates are slower in their respective settings. This is likely unavoidable in our setting, as we show that our convergence rate is optimal in this more general setting (see Appendix E).

A related class of algorithms that are similar to Alecton is stochastic power iteration . These algorithms reconsider (1) as an eigenvalue problem, and uses the familiar power iteration algorithm, adapted to a stochastic setting. Stochastic power iteration has been applied to a wide variety of problems . Oja show convergence of this algorithm, but provides no rate. Arora et al. analyze this problem, and state that “obtaining a theoretical understanding of the stochastic power method, or of how the step size should be set, has proved elusive.” Our paper addresses this by providing a method for selecting the step size, although our analysis shows convergence for any sufficiently small step size.

Shamir provide exponential-rate local convergence results for a stochastic power iteration algorithm for PCA. As they note, it can be used in practice to improve the accuracy of an estimate returned by another, globally-convergent algorithm such as Alecton.

Also recently, Balsubramani et al. and Hardt and Price provide a global convergence rate for the stochastic power iteration algorithm. Our result only depends on the variance of the samples, while both their results require absolute bounds on the magnitude of the noise. This allows us to analyze a different class of noise models, which enables us to do matrix completion, phase retrieval, and subspace tracking in the same model.

Algorithmic Derivation

The low-rank factorization introduces symmetry into the problem. If we let

where $G_{x}$ is the matrix such that for all $u$ and $v$,

where the right side of this equation denotes the Riemannian metric of the manifold at $x$. For (2), the manifold in question is

For Alecton, we are free to pick any Riemannian metric and step size. Inspired by (3), we pick a new step size parameter $\eta$, and let $\alpha_{k}=\frac{1}{4}\eta$ and set

For $p=1$, choosing a Riemannian metric to use with SGD results in the same algorithm as choosing an SGD step size that depends on the iterate $Y_{k}$. The same update rule would result if we substituted

into the standard SGD update formula. We can think of this as the manifold results giving us intuition on how to set our step size.

The reason why selecting this particular step size/metric is useful in practice is that we can run the simpler update rule

We can use (4) to compute the column space (or “angular component”) of the solution, before then recovering the rest of the solution (the “radial component”) using averaging. Doing this corresponds to Algorithm 1, Alecton. Notice that, unlike most iterative algorithms for matrix recovery, Alecton does not require any special initialization phase and can be initialized randomly.

Analyzing this algorithm is challenging, as the low-rank decomposition also introduces symmetrical families of fixed points. Not all these points are globally optimal: in fact, a fixed point will occur whenever

One consequence of the non-optimal fixed points is that the standard proof of SGD’s convergence, in which we choose a Lyapunov function and show that this function’s expectation decreases with time, cannot work. This is because, if such a Lyapunov function were to exist, it would show that no matter where we initialize the iteration, convergence to a global optimum will still occur rapidly; this cannot be possible due to the presence of the non-optimal fixed points. Thus, a standard statement of global convergence, that convergence occurs uniformly regardless of initial condition, cannot hold.

We therefore use martingale-based methods to show convergence. Specifically, our attack involves defining a process $x_{k}$ with respect to the natural filtration $\mathcal{F}_{k}$ of the iteration, such that $x_{k}$ is a supermartingale, that is $\mathbf{E}\left[x_{k+1}\middle|\mathcal{F}_{k}\right]\leq x_{k}$. We then use the optional stopping theorem to bound both the probability and rate of convergence of $x_{k}$, from which we derive convergence of the original algorithm. We describe this analysis in the next section.

Convergence Analysis

First, we need a way to define convergence for the angular phase. For most problems, we want $C(Y_{k})$ to be as close as possible to the span of $u_{1},u_{2},\ldots,u_{p}$. However, for some cases, this is not what we want. For example, consider the case where $p=1$ but $\lambda_{1}=\lambda_{2}$. In this case, the algorithm could not recover $u_{1}$, since it is indistinguishable from $u_{2}$. Instead, it is reasonable to expect $C(Y_{k})$ to converge to the span of $u_{1}$ and $u_{2}$.

To handle this case, we instead want to measure convergence to the subspace spanned by some number, $q\geq p$, of the algebraically largest eigenvectors (in most cases, $q=p$). For a particular $q$, let $U$ be the projection matrix onto the subspace spanned by $u_{1},u_{2},\ldots,u_{q}$, and define $\Delta$, the eigengap, as $\Delta=\lambda_{q}-\lambda_{q+1}$. We now let $\epsilon>0$ be an arbitrary error term, and define an angular success condition for Alecton.

This condition requires that all members of the column space of $Y_{k}$ are close to the desired subspace. We say that success has occurred by time $t$ if success has occurred for some timestep $k<t$. Otherwise, we say the algorithm has failed, and we let $F_{t}$ denote this failure event.

To prove convergence, we need to put some restrictions on the problem. Our theorem requires the following three conditions.

In Section 4, we show several models that satisfy AVC.

Most of the noise models we analyze have rank-1 samples, and so satisfy the rank condition.

This represents a constant step size parameter that is independent of problem scaling. An instance of Alecton satisfies the Alecton Step Size Condition if and only if $\gamma\leq 1$.

Note that the step size condition is only an upper bound on the step size. This means that, even if we do not know the problem parameters exactly, we can still choose a feasible step size as long as we can bound them. (However, smaller step sizes imply slower convergence, so it is a good idea to choose $\eta$ as large as possible.)

We will now define a useful function, then state our main theorem that bounds the probability of failure.

Assume that we run an instance of Alecton that satisfies the variance, rank, and step size conditions. Then for any $t$, the probability that the angular phase will have failed up to time $t$ is

Also, in the radial phase, for any constant $\psi$ it holds that

In particular, if $\sigma_{a}\Delta^{-1}$ does not vary with $n$, this theorem implies convergence of the angular phase with constant probability after $O(\epsilon^{-1}np^{3}\log n)$ iterations and in the same amount of time. Note that since we do not reuse samples in Alecton, our rates do not differentiate between sampling and computational complexity, unlike many other algorithms (see Appendix B). We also do not consider numerical error or overflow: periodically re-normalizing the iterate may be necessary to prevent these in an implementation of Alecton.

Since the upper bound expression uses $Z_{p}$, which is obscure, we plot it here (Figure 1). We also can make a more precise statement about the failure rate for $p=1$.

A proof for Theorem 1 and full formal definitions will appear in Appendix C of this document, but since the method is nonstandard for non-convex optimization (although it has been used in Shamir to show convergence for convex problems), we will outline it here. First, we define a failure event $f_{k}$ at each timestep, that occurs if the iterate gets “too close” to the unstable fixed points. Next, we define a sequence $\tau_{k}$, where

(where $\left|X\right|$ denotes the determinant of $X$); the intuition here is that $\tau_{k}$ is close to $1$ if and only if success occurs, and close to when failure occurs. We show that, if neither success nor failure occurs at time $k$,

for some constant $R$; here, $\mathcal{F}_{k}$ denotes the filtration at time $k$, which contains all the events that have occurred up to time $k$ . If we let $T$ denote the first time at which either success or failure occurs, then this implies that $\tau_{k}$ is a submartingale for $k<T$. We use the optional stopping Theorem (here we state a discrete-time version).

A random variable $T$ is a stopping time with respect to a filtration $\mathcal{F}_{k}$ if and only if $\left\{T\leq k\right\}\in\mathcal{F}_{k}$ for all $k$. That is, we can tell whether $T\leq k$ using only events that have occurred up to time $k$.

If $x_{k}$ is a martingale (or submartingale) with respect to a filtration $\mathcal{F}_{k}$, and $T$ is a stopping time with respect to the same filtration, then $x_{k\wedge T}$ is also a martingale (resp. submartingale) with respect to the same filtration, where $k\wedge T$ denotes the minimum of $k$ and $T$. In particular, for bounded submartingales, this implies that $\mathbf{E}\left[x_{0}\right]\leq\mathbf{E}\left[x_{T}\right]$.

Here, $T$ is a stopping time since it depends only on events occurring before timestep $T$. Applying this to the submartingale $\tau_{k}$ results in

This isolates the probability of the failure event occurring. Next, subtracting $1$ from both sides of (6) and taking the logarithm results in

So, if we let $W_{k}=\log(1-\tau_{k})+R\delta k$, then $W_{k}$ is a supermartingale. We again apply the optional stopping theorem to produce

This isolates the expected value of the stopping time. Finally, we notice that success occurs before time $t$ if $T\leq t$ and $f_{T}$ does not occur. By the union bound, this implies that

Substituting the isolated values for $\underset{}{\mathbf{P}}\left(f_{T}\right)$ and $\mathbf{E}\left[T\right]$ produces the expression above in (5).

Application Examples

One sampling distribution that arises in many applications (most importantly, matrix completion ) is entrywise sampling. This occurs when the samples are independently chosen from the entries of $A$. Specifically,

where $i$ and $j$ are each independently drawn from ${1,\ldots,n}$. It is standard for these types of problems to introduce a matrix coherence bound .

For problems in which the matrix $A$ is of constant rank, and its eigenvalues do not vary with $n$, neither $\left\|A\right\|_{F}$ nor $\mathbf{tr}\left(A\right)$ will vary with $n$. In this case, $\sigma_{a}^{2}$, $\sigma_{r}^{2}$, and $\Delta$ will be constants, and the $O(\epsilon^{-1}n\log n)$ bound on convergence time will hold.

2 Rectangular Entrywise Sampling

Entrywise sampling also commonly appear in rectangular matrix recovery problems. In these cases, we are trying to solve something like

To solve this problem using Alecton, we first convert it into a symmetric matrix problem by constructing the block matrix

it is known that recovering the dominant eigenvectors of $A$ is equivalent to recovering the dominant singular vectors of $M$.

In the case where we can bound the entries of $M$ (this is natural for recommender systems), we can prove the following.

for all $i$ and $j$, then the rectangular entrywise sampling distribution on $M$ satisfies the Alecton variance condition with parameters

As above, for problems in which the singular values of $M$ do not vary with problem size, our big-$O$ convergence time bound will still hold.

3 Trace Sampling

Another common sampling distribution arises from the matrix sensing problem . In this problem, we are given the value of $v^{T}Aw$ for unit vectors $v$ and $w$ selected uniformly at random. (This problem has been handled for the more general complex case in using Wirtinger flow.) Using a trace sample, we can construct an unbiased sample

This lets us bound the variance as follows.

As above, for problems in which the eigenvalues of $A$ do not vary with problem size, our big-$O$ convergence time bound will still hold.

In some cases of the trace sampling problem, instead of being given samples of the form $u^{T}Av$, we know $u^{T}Au$. In this case, we need to use two independent samples $u_{1}^{T}Au_{1}$ and $u_{2}^{T}Au_{2}$, and let $u\propto u_{1}+u_{2}$ and $v\propto u_{1}-u_{2}$ be two unit vectors which we will use in the above sampling scheme. Notice that since $u_{1}$ and $u_{2}$ are independent and uniformly distributed, $u$ and $v$ will also be independent and uniformly distributed (by the spherical symmetry of the underlying distribution). Furthermore, we can compute

This allows us to use our above trace sampling scheme even with samples of the form $u^{T}Au$.

4 Subspace Sampling

Our analysis can handle more complicated sampling schemes. Consider the following distribution, which arises in subspace tracking . Our matrix $A$ is a rank-$r$ projection matrix, and each sample consists of some randomly-selected entries from a randomly-selected vector in its column space. Specifically, we are given $Qv$ and $Rv$, where $v$ is some vector selected uniformly at random from $C(A)$, and $Q$ and $R$ are independent random diagonal projection matrices with expected value $mn^{-1}I$. Using this, we can construct the distribution

This distribution is unbiased since $\mathbf{E}\left[qvv^{T}\right]=A$. When bounding its second moment, we run into the same coherence problem as we did in the entrywise case, which motivates us to introduce a coherence constraint for subspaces.

Using this, we can prove the following facts about the second moment of this distribution.

The subspace sampling distribution, when sampled from a subspace that is incoherent with parameter $\mu$, satisfies the Alecton variance condition with parameters

In many cases of subspace sampling, we are given just some entries of $v$ at each timestep (as opposed to two separate random sets of entries associated with $Q$ and $R$). That is, we are given a random diagonal projection matrix $S$, and the product $Sv$. We can use this to construct a sample of the above form by randomly splitting the given entries among $Q$ and $R$ in such a way that $Q=QS$ and $R=RS$, and $Q$ and $R$ are independent. We can then construct an unbiased sample as

which uses only the entries of $v$ that we are given.

5 Noisy Sampling

Since our analysis depends only on a variance bound, it is straightforward to handle the case in which the values of our samples themselves are noisy. Using the additive property of the variance for independent random variables, we can show that additive noise only increases the variance of the sampling distribution by a constant amount proportional to the variance of the noise. Similarly, using the multiplicative property of the variance for independent random variables, multiplicative noise only multiplies the variance of the sampling distribution by a constant factor proportional to the variance of the noise. In either case, we can show that the noisy sampling distribution satisfies AVC.

6 Extension to Higher Ranks

It is possible to use multiple iterations of the rank-$1$ version of Alecton to recover additional eigenvalue/eigenvector pairs of the data matrix $A$ one-at-a-time. This is a standard technique for using power iteration algorithms to recover multiple eigenvalues. Sometimes, this may be preferable to using a single higher-rank invocation of Alecton (for example, we may not know a priori how many eigenvectors we want). We outline this technique as Algorithm 2.

This strategy allows us to recover the largest $p$ eigenvectors of $A$ using $p$ executions of Alecton. If the eigenvalues of the matrix are independent of $n$ and $p$, we will be able to accomplish this in $O(\epsilon^{-1}pn\log n)$ total steps.

Experiments

We experimentally verify our main claim, that Alecton does converge quickly for practical datasets.

All experiments were run on a machine with a single twelve-core socket (Intel Xeon E5-2697, 2.70GHz), and 256 GB of shared memory. All were written in C++, excepting the Netflix Prize problem experiment, which was written in Julia. No data was collected for the radial phase of Alecton, since the performance of averaging is already well understood.

Figure 2(b) illustrates the performance of Alecton ($p=q=1$ again) on a larger dataset with $n=10^{6}$ as the step size parameter $\eta$ is varied. As we would expect, a smaller value of $\eta$ yields slower, but more accurate convergence. Also notice that the smaller the value of $\eta$, the more the initial value seems to affect convergence time.

Figure 3 demonstrates convergence results on real data from the Netflix Prize problem. This problem involves recovering a matrix with 480,189 columns and 17,770 rows from a training dataset containing 110,198,805 revealed entries. We used the rectangular entrywise distribution described above, then ran Alecton with $\eta=10^{-12}$ and $p=q=1$ for ten million iterations to recover the most significant singular vector. Next, we used Algorithm 2 to recover additional singular vectors of the matrix, up to a maximum of $p=12$. The absolute runtime and RMS errors after the recovery of each subsequent eigenvector are plotted in Figure 3. This plot illustrates that the runtime of the one-at-a-time algorithm does not increase disastrously as the number of recovered eigenvectors expands.

The Hogwild! algorithm is a parallel, lock-free version of stochastic gradient descent that has been shown to perform similarly to sequential SGD on convex problems, while allowing for a good parallel speedup. It is an open question whether a Hogwild! version of Alecton for non-convex problems converges with a good rate, but we are optimistic that it will.

Conclusion

This paper exhibited Alecton, a stochastic gradient descent algorithm applied to a non-convex low-rank factorized problem; it is similar to the algorithms used in practice to solve a wide variety of problems. We prove that Alecton converges globally, and provide a rate of convergence. We do not require any special initialization step but rather initialize randomly. Furthermore, our result depends only on the variance of the samples, and therefore holds under broad sampling conditions that include both matrix completion and matrix sensing, and is also able to take noisy samples into account. We show these results using a martingale-based technique that is novel in the space of non-convex optimization, and we are optimistic that this technique can be applied to other problems in the future.

References

Appendix A Negative Results

Here, we observe what happens when we choose a constant step size for stochastic gradient descent for quartic objective functions. Consider the simple optimization problem of minimizing

This function will have gradient descent update rule

We now prove that, for any reasonable step size rule chosen independently of $x_{k}$, there is some initial condition such that this iteration diverges to infinity.

Assume that we iterate using the above rule, for some choice of $\alpha_{k}$ that is not super-exponentially decreasing; that is, for some $C>1$ and some $\alpha>0$, $\alpha_{k}\geq\alpha C^{-2k}$ for all $k$. Then, if $x_{0}^{2}\geq\alpha^{-1}(C+1)$, for all $k$

We will prove this by induction. The base case follows directly from the assumption, while under the inductive case, if the proposition is true for $k$, then

This proof shows that, for some choice of $x_{0}$, $x_{k}$ will diverge to infinity exponentially quickly. Furthermore, no reasonable choice of $\alpha_{k}$ will be able to halt this increase for all initial conditions. We can see the effect of this in stochastic gradient descent as well, where there is always some probability that, due to an unfortunate series of gradient steps, we will enter the zone in which divergence occurs. On the other hand, if we chose step size $\alpha_{k}=\gamma_{k}x_{k}^{-2}$, for some $0<\gamma_{k}<2$, then

which converges for all starting values of $x_{k}$. This simple example is what motivates us to take $\left\|Y_{k}\right\|$ into account when choosing the step size for Alecton.

Now, we know that $e_{1}$ is the most significant eigenvector of $A$, and that $y=2e_{1}$ is the global solution to the problem. However,

. This implies that if $e_{1}^{T}y_{0}=0$, then $e_{1}^{T}y_{k}=0$ for all $k$, which means that convergence to the global optimum cannot occur. This illustrates that global convergence does not occur for all manifold optimization problems using a low-rank factorization and for all starting points.

We might think that our results can be generalized to give $O(n\log n)$ convergence of low-rank factorized problems with arbitrary constraints. Here, we show that this will not work for all problems by encoding an NP-complete problem as a constrained low-rank optimization problem.

For any graph with node set $N$ and edge set $E$, the MAXCUT problem on the graph requires us to solve

Equivalently, if we let $A$ denote the edge-matrix of the graph, we can represent this as a matrix problem

Since the diagonal of $A$ is zero, if we fix all but one of the entries of $y$, the objective function will have an affine dependence on that entry. In particular, this means that a global minimum of the problem must occur on the boundary where $y_{i}\in\{-1,1\}$, which implies that this problem has the same global solution as the original MAXCUT problem. Furthermore, for sufficiently large values of $\sigma$, the problem

will also have the same solution. But, this problem is in the same form as a low-rank factorization of

where $X=yy^{T}$. Since MAXCUT is NP-complete, it can’t possibly be the case that SGD applied to this low-rank factorized problem converges quickly to the global optimum, because that would imply an efficient solution to this NP-complete problem. This suggests that care will be needed when analyzing problems with constraints, in order to exclude these sorts of cases.

Appendix B Comparison with Other Methods

There are several other algorithms that solve similar matrix recover problems in the literature. In Table B, we list some other algorithms, and their convergence rates, in terms of both number of samples required (sampling complexity) and number of iterations performed (computational complexity). For this table, the data is assumed to be of dimension $n$, and the rank (where applicable) is assumed to be $p$. (In order to save space, factors of $\log\log\epsilon^{-1}$ have been omitted from some formulas.)

Appendix C Proofs of Main Results

In this appendix, we provide rigorous definitions and detail the proof outlined in Section 3.1.

Fleming and Harrington provide the following definitions of filtration and martingale. We state the definitions adapted to the discrete-time case.

Given a measurable probability space $(\Omega,\mathcal{F})$, a filtration is a sequence of sub-$\sigma$-algebras $\left\{\mathcal{F}_{t}\right\}$ for $t\geq 0$, such that for all $s\leq t$,

That is, if an event $A$ is in $\mathcal{F}_{s}$, and $t\geq s$, then $A$ is also in $\mathcal{F}_{t}$. This definition encodes the monotonic increase in available information over time.

Let $\left\{X_{t}\right\}$ be a stochastic process and $\left\{\mathcal{F}_{t}\right\}$ be a filtration over the same probability space. Then $X$ is called a martingale with respect to the filtration if for every $t$, $X_{t}$ is $\mathcal{F}_{t}$-measurable, and

We call $X$ a submartingale if the same conditions hold, except (7) is replaced with

We call $X$ a supermartingale if the same conditions hold, except (7) is replaced with

C.2 Preliminaries

In addition to the quantities used in the statement of Theorem 1, we let

and define sequences $\tau_{k}$ and $\phi_{k}$ as

This agrees with the definition of $\tau_{k}$ stated in the body of the paper. Using this sequence, we define the failure event $f_{k}$ as the event that occurs when

Finally, we define $T$, the stopping time, to be the first time at which either the success event or the failure event occurs.

Now, we state some lemmas we will need in the following proofs. We defer proofs of the lemmas themselves to Appendix D. First, we state a lemma about quadratic rational functions that we will need in the next section.

Next, a lemma about the expected initial value of $\tau$:

If we initialize $Y_{0}$ uniformly as in the Alecton algorithm, then

Next, a lemmas that bounds a determinant expression.

Next, a lemma that bounds $\tau$ in the case that the success condition does not occur.

If we run Alecton, and at timestep $k$, the success condition does not hold, then

Finally, a lemma that relates $\phi$ and $\tau$.

Using the definitions above, for all $k$,

C.3 Main Proofs

We now proceed to prove Theorem 1 in six steps, as outlined in Section 3.1.

First, we prove Lemma 11, the dominant mass bound lemma, which bounds $\mathbf{E}\left[\tau_{k+1}\middle|\mathcal{F}_{k}\right]$ from below by a quadratic function of the step size $\eta$.

We use this to prove Lemma 12, which establishes the result stated in (6).

We use the optional stopping theorem to prove Lemma 13, which bounds the probability of a failure event occurring before success.

We use the optional stopping theorem again to prove Lemma 14, which bounds the expected time until either a failure or success event occurs.

We use Markov’s inequality and the union bound to bound the angular failure probability of Theorem 1.

Finally, we prove the radial phase result stated in Theorem 1.

If we run Alecton under the conditions of Theorem 1, then for any $k$,

From the definition of $\tau$, at the next timestep we will have

Now, since $U$ commutes with $A$, we will have that

Since our instance of Alecton satisfies the variance condition, and $W$ commutes with $A$,

and therefore, since $\mathbf{tr}\left(W\right)\leq q+1$,

Finally, since for our chosen value of $\gamma$,

If we run Alecton under the conditions of Theorem 1, then for any time $k$ at which neither the success event nor the failure event occur,

for sequence $S_{k}$. Now, it can be easily verified that we chose $\gamma$ such that

Substituting this in to our original expression produces

If we run Alecton under the conditions of Theorem 1, then the probability that the failure event will occur before the success event is

To prove this, we use the stopping time $T$, which we defined as the first time at which either the success event or failure event occurs. First, if $k<T$, it follows that neither success nor failure have occurred yet, so we can apply Lemma 12, which results in

Therefore $\tau_{k}$ is a supermartingale for $k<T$. So, we can apply the optional stopping theorem, which produces

where $f_{T}$ is the failure event at time $T$. Applying the definition of the failure event from (8),

Therefore, solving for $\underset{}{\mathbf{P}}\left(f_{T}\right)$,

If we run Alecton under the conditions of Theorem 1, then the expected value of the stopping time $T$ will be

First, as above if $k<T$, we can apply Lemma 12, which results in

Now, if $k<T$, then since failure hasn’t occurred yet, $\tau_{k}>\frac{1}{2}$. So,

Now, since the logarithm function is concave, by Jensen’s inequality we have

Now, we define a new process $\psi_{k}$ as

so $\psi_{k}$ is a supermartingale for $k<T$. We can therefore apply the optional stopping theorem, which states that

Since $1-\tau_{0}<1$, it follows that $\log(1-\tau_{0})<0$. Therefore,

Solving for the expected value of the stopping time,

Finally, substituting $\eta$ in terms of $\gamma$ results in

First, we notice that the total failure event up to time $t$ can be written as

That is, total failure up to time $t$ occurs if either failure happens before success (event $f_{T}$), or neither success nor failure happen before $t$. By the union bound,

Finally, applying Lemmas 13 and 14 produces

Recall that in Alecton, $\bar{R}$ is defined as

Now, computing the expected distance to the mean,

Applying the Alecton variance condition, and recalling that $\mathbf{tr}\left(\hat{Y}\hat{Y}^{T}\right)=p$, results in

We can now apply Markov’s inequality to this expression. This results in, for any constant $\psi>0$,

Appendix D Proofs of Lemmas

First, we prove the lemmas used above to demonstrate the general result.

Dividing both sides by $1+bx+cx^{2}$ (which we can do since this is assumed to be positive) reconstructs the desired identity. ∎

We first note that, by the symmetry of the multivariate Gaussian distribution, initializing $Y_{0}$ uniformly at random such that $Y_{0}^{T}Y_{0}=I$ is equivalent to initializing the entries of $Y_{0}$ as independent standard normal random variables, for the purposes of computing $\tau_{0}$. Under this initialization strategy, $\mathbf{E}\left[\tau_{0}\right]$ is

Since $X$ and $Z$ are selected orthogonally from a Gaussian random matrix, they must be independent, so we can take their expected values independently. Taking the expected value first with respect to $Z$, we notice that $\left|V\right|^{-1}$ is a convex function in $V$, and so by Jensen’s inequality,

We will prove this separately for each case. First, if $m=1$, then $Y$ is a vector, and the desired expression simplifies to

Straightforward evaluation indicates that this expression holds in this case.

Next, we consider the case where $B$ is rank-1. In this case, we can rewrite it as $B=uv^{T}$ for vectors $u$ and $v$, such that $u^{T}Zu=1$. Then,

Applying the matrix determinant lemma, and recalling that

Rewriting this in terms of the matrix $B=uv^{T}$,

Substitution produces the desired result. ∎

First, for the lower bound, we notice that

since the interior of the left expression is a projection matrix. This lets us conclude that

Appling this to the result of Lemma 15 produces the desired lower bound.

For the upper bound, recall that, by the Cauchy-Schwarz inequality, for any rank-1 matrix $A$,

Appling this to the result of Lemma 15 produces the desired upper bound. ∎

For any symmetric matrix $0\preceq X\preceq I$,

If $x_{1},x_{2},\ldots,x_{p}$ are the eigenvalues of $x$, then this statement is equivalent to

If we let $f(X)$ denote this expression, then

It follows that the minimum of $f$ is attained at $X=I$. However, when $X=I$, $f(X)=0$, and so $f>0$, which proves the lemma. ∎

From the definition of $\phi_{k}$, if we let $Z^{2}=\left(Y_{k}^{T}WY_{k}\right)^{-1}$ for $Z$ positive semidefinite, then

Since $0\preceq ZY_{k}^{T}U^{T}UY_{k}Z\preceq I$, we can apply Lemma 16, which produces

and define $\hat{z}$ as the unit vector such that

It follows that $\hat{Y}_{k}^{T}U\hat{Y}_{k}$ has an eigenvalues less than $1-\epsilon$.

Since this is a matrix that has eigenvalues between and $1$, it follows that its determinant is less than each of its eigenvalues. From the analysis above, we can bound one of the eigenvalues of this matrix. Doing this results in

By the definition of expected value, since $x$ is normally distributed,

If we let $\mathcal{F}$ denote the fourier transform, then

Furthermore, since the Gaussian functions are eigenfunctions of the Fourier transform, we know that

Letting $u=\frac{\omega+a}{\sqrt{2}}$ and $d\omega=\sqrt{2}du$, so

This is the desired expression. Furthermore, since for all $x$,

we can also produce the desired upper bound on $Z_{1}$,

Next, we prove the Alecton Variance Conditions lemmas for the distributions mentioned in the body of the paper.

To analyze the entrywise sampling case, we need some lemmas that makes the incoherence condition more accessible.

If matrix $A$ is symmetric and incoherent with parameter $\mu$, and $B$ is a symmetric matrix that commutes with $A$, then $B$ is incoherent with parameter $\mu$.

Since $A$ and $B$ commute, they must have the same eigenvectors. Therefore, the set of eigenvectors that shows that $A$ is incoherent with parameter $\mu$ will also show that $B$ has the same property. ∎

If matrix $A$ is symmetric and incoherent with parameter $\mu$, and $e_{i}$ is a standard basis element, then

Let $u_{1},u_{2},\ldots,u_{n}$ be the eigenvectors guaranteed by the incoherence of $A$, and let $\lambda_{1},\ldots,\lambda_{n}$ be the corresponding eigenvalues. Then,

We recall that the entrywise samples are of the form

where $u$ and $v$ are independently, uniformly chosen standard basis elements. We further recall that $\mathbf{E}\left[uu^{T}\right]=\mathbf{E}\left[vv^{T}\right]=n^{-1}I$. Now, evaluating the desired quantity,

Since $W$ commutes with $A$, by Lemmas 18 and 19, $u^{T}Wu\leq\mu^{2}n^{-1}\mathbf{tr}\left(W\right)$. Therefore,

Since $A^{2}$ commutes with $A$, the same logic shows that $v^{T}A^{2}v\leq\mu^{2}n^{-1}\mathbf{tr}\left(A^{2}\right)$, and so,

So it suffices to choose $\sigma_{a}^{2}=\mu^{4}\left\|A\right\|_{F}^{2}$, as desired. ∎

and by Lemma 19, $u^{T}Au\leq\mu^{2}n^{-1}\mathbf{tr}\left(A\right)$, and so

So it suffices to choose $\sigma_{r}^{2}=\mu^{4}\mathbf{tr}\left(A\right)^{2}$, as desired. ∎

D.1.2 Rectangular Entrywise Sampling

We recall that the rectangular entrywise samples are of the form

Now, since $(x+y)^{2}\leq 2(x^{2}+y^{2})$, if we let $P$ be the projection matrix onto the first $m$ basis vectors, then $\mathbf{E}\left[e_{i}e_{i}^{T}\right]=m^{-1}P$ and $\mathbf{E}\left[e_{m+j}e_{m+j}^{T}\right]=n^{-1}(I-P)$, and so,

Since this is true for any $y$ and $z$, it is true in particular for $z$ being an eigenvector of $A$. Therefore, it suffices to pick $\sigma_{a}^{2}=2\xi\left\|M\right\|_{F}^{2}$. Similarly, it is true in particular for $z=y$, and therefore it suffices to pick $\sigma_{r}^{2}=2\xi\left\|M\right\|_{F}^{2}$. This proves the lemma. ∎

D.1.3 Trace Sampling

In order to prove our second moment lemma for the trace sampling case, we must first derive some lemmas about the way this distribution behaves.

If we let $u$ denote $y^{T}x$, and $z$ denote the components of $x$ orthogonal to $y$, then $\left\|x\right\|^{2}=u^{2}+\left\|z\right\|^{2}$. Furthermore, by the properties of the normal distribution, $u$ and $z$ are independent. Therefore,

Now, $\mathbf{E}\left[u^{4}\right]$ is the fourth moment of the normal distribution, which is known to be $3$. Furthermore, $\mathbf{E}\left[\left\|z\right\|^{-4}\right]$ is the second moment of an inverse-chi-squared distribution with parameter $n-1$, which is also a known result. Substituting these in,

This quantity has the asymptotic properties we want. In particular, applying the constraint that $n>50$,

be the eigendecomposition of $W$. Then for any unit vector $z$,

By the Cauchy-Schwarz inequality applied to the expectation,

By Lemma 20, $\mathbf{E}\left[(z^{T}v)^{4}\right]\leq 4n^{-2}$, and so

Since this is true for any unit vector $z$, by the definition of the positive semidefinite relation,

Now, we prove the AVC lemma for this distribution.

Evaluating the expression we want to bound,

So it suffices to pick $\sigma_{a}^{2}=16\left\|A\right\|_{F}^{2}$, as desired. ∎

Evaluating the expression we want to bound,

So it suffices to pick $\sigma_{r}^{2}=16\left\|A\right\|_{F}^{2}$, as desired. ∎

D.1.4 Subspace Sampling

Recall that, in subspace sampling, our samples are of the form

where $Q$ and $R$ are independent projection matrices that select $m$ entries uniformly at random, and $v$ is uniformly and independently selected from the column space of $A$. Using this, we first prove some lemmas, then prove our bounds.

If $Q$ is a projection matrix that projects onto a subspace spanned by $m$ random standard basis vectors, and $v$ is a member of a subspace that is incoherent with parameter $\mu$, then for any vector $x$,

As a corollary, for any symmetric matrix $W\succeq 0$,

Let $\lambda_{i}$ be $1$ in the event that $e_{i}$ is in the column space of $Q$, and otherwise. Then an eigendecomposition of $Q$ is

Taking the expected value, and noting that $\lambda_{i}$ and $\lambda_{j}$ are independent, and have expected value $\mathbf{E}\left[\lambda_{i}\right]=mn^{-1}$,

Since $v$ is part of a subspace that is incoherent,

Evaluating the expression we want to bound,

So, we can choose $\sigma_{a}^{2}=r^{2}(1+\mu rm^{-1})^{2}$, as desired. ∎

Evaluating the expression we want to bound,

So, we can choose $\sigma_{r}^{2}=r^{2}(1+\mu rm^{-1})^{2}$, as desired. ∎

Appendix E Lower Bound on Alecton Rate

In this section, we prove a rough lower bound on the rate of convergence of an Alecton-like algorithm for bounded sampling distributions. Specifically, we analyze the case where, rather than choosing a constant $\eta$, we allow the step size to vary at each timestep. Our result shows that we can’t hope for a better step size rule that improves the convergence rate of Alecton to, for example, a linear rate.

To show this lower bound, we assume we run Alecton with $p=1$ for some sampling distribution such that for all $\eta$ and all $y$, for some constant $C$,

Further assume that for some eigenvector $u$ (with eigenvalue $\lambda\geq 0$) that is not global solution, the sample variance in the direction of $u$ satisfies

This quantity measures the error of the iterate at timestep $k$ in the direction of $u$. We will show that the expected value of $\rho_{k}$ can only decrease with at best a $\Omega\left(\frac{1}{K+1}\right)$ rate.

For any $a\geq 0$, $b\geq 0$, and $0\leq x\leq 1$,

Under the above conditions, regardless of how we choose the step size in the Alecton algorithm, even if we are able to choose a different step size each iteration, the expected error will still satisfy

Using the Alecton update rule with a time-varying step size $\eta_{k}$,

Since, by symmetry, $\mathbf{E}\left[\rho_{0}\right]=n^{-1}$, we have

Appendix F Handling Constraints

Alecton can easily be adapted to solve the problem of finding a low-rank approximation to a matrix under a spectahedral constraint. That is, we want to solve the problem

This is equivalent to the decomposed problem

This will have a minimum when $y=u_{1}$. We can therefore solve the problem using only the angular phase of Alecton, which recovers the vector $u_{1}$. The same convergence analysis described above still applies.

For an example of a constrained problem that Alecton cannot handle, because it is NP-hard, see the elliptope-constrained MAXCUT embedding in Appendix A. This shows that constrained problems can’t be solved efficiently by SGD algorithms in all cases.

Appendix G Towards a Linear Rate

is $X=A$. Performing a rank-$p$ quadratic substitution on this problem results in:

The specific case we will be looking at is where the operator $\Omega$ satisfies the $p$-RIP constraint.

This definition encodes the notion that $\Omega$ preserves the norm of low-rank matrices under its transformation. We can prove a simple lemma that extends this to the inner product.

If $\Omega$ is $(p+q)$-RIP with parameter $\delta$, then for any symmetric matrices $X$ and $Y$ of rank at most $p$ and $q$ respectively,

Since $\mathbf{rank}\left(X-aY\right)\leq\mathbf{rank}\left(X\right)+\mathbf{rank}\left(Y\right)\leq p+q$, we can apply our RIP inequalities, which produces

Substituting $a=\frac{\left\|X\right\|_{F}}{\left\|Y\right\|_{F}}$ results in

Finally, we prove our main theorem that shows that the quadratically transformed objective function is strongly convex in a ball about the solution.

and $\Omega$ is $3p$-RIP with parameter $\delta$, then for all $Y$, if we let $\lambda_{p}$ denote the smallest positive eigenvalue of $A$ then

The directional derivative of $f$ along some direction $V$ will be, by the product rule,

The second derivative along this same direction will be

To this, we can apply the definition of RIP, and the corollary lemma, which results in

Now, since at the optimum, $\lambda_{\text{min}}(Y^{T}Y)=\lambda_{p}$, it follows that for general $Y$,

Substituting this in to the previous expression,

Since this is true for an arbitrary direction vector $U$, it follows that

This theorem shows that there is a region of size $O(1)$ (i.e. not dependent on $n$) within which the above problem is strongly convex. So, if we start within this region, any standard convex descent method will converge at a linear rate. In particular, coordinate descent will do so. Therefore, we can imagine doing the following:

First, use Alecton to, with high probability, recover an estimate $Y$ that for which $\left\|YY^{T}-A\right\|_{F}$ is sufficiently small for the objective function to be strongly convex with some probability. This will only require $O(n\log n)$ steps of the angular phase of the algorithm per iteration of Alecton, as stated in the main body of the paper. We will need $p$ iterations of the algorithm to recover a rank-$p$ estimate, so a total $O(np\log n)$ iterations will be required.

Use a descent method, such as coordinate descent, to recover additional precision of the estimate. This method is necessarily more heavyweight than an SGD scheme (see Section E for the reason why an SGD scheme cannot achieve a linear rate), but it will converge monotonically at a linear rate to the exact solution matrix $A$.

This hybrid method is in some sense a best-of-both worlds approach. We use fast SGD steps when we can afford to, and then switch to slower coordinate descent steps when we need additional precision.