Previously, we used the SVD to find the best rank-\(k\) approximation to a matrix \(\mathbf{X}\). For this application, notice that we didn’t really need to find all the singular vectors and values.

We can save time by computing an approximate solution to the SVD. In particular, we will only compute the top \(k\) singular vectors and values. We can do this with iterative algorithms that achieve time complexity \(O(ndk)\) instead of \(O(nd^2)\). There are many algorithms for this problem:

  • Krylov subspace methods like the Lanczos method are most commonly used in practice.

  • Power method is the simplest Krylov subspace method and still works very well.

The power method computes the top right singular vector, \(\mathbf{v}_1\), of a matrix \(\mathbf{X} \in \mathbb{R}^{n\times d}\) with singular value decomposition \(\mathbf{X} = \mathbf{U\Sigma V}^\top\). The method is as follows:

Power Method:

  • Choose \(\mathbf{z}^{(0)}\) randomly. E.g. \(\mathbf{z}^{(0)} \sim \mathcal{N}(0,1)\).
  • \(\mathbf{z}^{(0)} = \mathbf{z}^{(0)} /\|\mathbf{z}^{(0)}\|_2\)
  • For \(i = 1,\ldots, q\)
    • \(\mathbf{z}^{(i)} = \mathbf{X}^\top \mathbf{X} \mathbf{z}^{(i-1)} )\)
    • \(n_i = \|\mathbf{z}^{(i)}\|_2\)
    • \(\mathbf{z}^{(i)} = \mathbf{z}^{(i)}/n_i\)
  • Return \(\mathbf{z}^{(q)}\)

In other words, each of our iterates \(\mathbf{z}^{(i)}\) is simply a scaling of a column in the following matrix: \[ K = \begin{bmatrix} \mathbf{z}^{(0)} & \mathbf{A}\mathbf{z}^{(0)} & \mathbf{A}^2 \mathbf{z}^{(0)} & \mathbf{A}^3\mathbf{z}^{(0)} \ldots \mathbf{A}^q \mathbf{z}^{(0)}\end{bmatrix}, \] where \(\mathbf{A} = \mathbf{X}^\top\mathbf{X}\). Typically we run \(q \ll d\) iterations, so \(K\) is a tall, narrow matrix. The span of \(K\)’s columns is called a Krylov subspace because it is generated by starting with a single vector \(\mathbf{z}^{(0)}\) and repeatedly multiplying by a fixed matrix \(A\). We will return to \(K\) shortly.

Analysis Preliminaries

Let \(\mathbf{v}_1, \ldots, \mathbf{v}_d\) be the right singular vectors of \(\mathbf{X}\) (the columns of \(\mathbf{V}\)). We write each iterate in terms of this basis of vectors: \[ \begin{align*} \mathbf{z}^{(0)} &= c_1^{(0)}\mathbf{v}_1 + {c}_2^{(0)}\mathbf{v}_2 + \ldots + c_d^{(0)} \mathbf{v}_d \\ \mathbf{z}^{(1)} &= c_1^{(1)}\mathbf{v}_1 + {c}_2^{(1)}\mathbf{v}_2 + \ldots + c_d^{(1)} \mathbf{v}_d\\ &\vdots\\ \mathbf{z}^{(q)} &= c_1^{(q)}\mathbf{v}_1 + {c}_2^{(i)}\mathbf{v}_2 + \ldots + c_d^{(q)} \mathbf{v}_d \end{align*} \]

Note that \(c_j^{(i)} = \langle \mathbf{z}^{(i)}, v_i\rangle\). Intuitively, we expect the relative size of \(c_1^{(i)}\) to increase in comparison to the other coefficients as \(i\) increases, which is good because then \(\mathbf{z}^{(q)}\) will start looking more and more like \(\mathbf{v}_1\). Why do we expect this? Because \(\mathbf{z}^{(i)} = \frac{1}{n_i}\mathbf{A}^q \mathbf{z}^{(i-1)}\) and thus we can check that \(c_j^{(i)} = \frac{1}{n_i}\sigma_j^2 c^{(i-1)}_j\). The \(\frac{1}{n_i}\) term is fixed for all \(j\), but the \(\sigma_j^2\) term will be largest for \(j=1\), and thus the \(c_1^{(i)}\) will increase in size more than any other term. We will analyze this formally soon.

First however, let’s see what it suffices to prove.

Claim 1: If \(c_j^{(q)}/c_1^{(q)} \leq \sqrt{\epsilon/d}\) for all \(j\neq 1\) than either \(\|\mathbf{v}_1 -\mathbf{z}^{(q)}\|_2^2\leq 2\epsilon\) or \(\|-\mathbf{v}_1 - \mathbf{z}^{(q)}\|_2^2\leq 2\epsilon\) .

Proof. Since \(c_1^{(q)} \leq 1\), the above would imply that \(c_j^{(q)} \leq \sqrt{\epsilon/d}\). Since \(\sum_{k=1}^d \left(c_k^{(q)}\right)^2 = 1\) it follows that \(\left(c_1^{(q)}\right)^2 \geq (1-\epsilon)\) and thus \(\left|c_1^{(q)}\right| \geq 1-\epsilon\). Then for any unit vector \(\mathbf{x}\) we have \(\|\mathbf{x} -\mathbf{z}^{(q)}\|_2^2 = 2 - 2\langle\mathbf{x},\mathbf{z}^{(q)}\rangle\). Since \(\langle\mathbf{v}_1,\mathbf{z}^{(q)}\rangle = c_1^{(q)}\) we conclude that either \(\langle\mathbf{v}_1,\mathbf{z}^{(q)}\rangle \geq 1-\epsilon\) or \(\langle- \mathbf{v}_1,\mathbf{z}^{(q)}\rangle \geq 1-\epsilon\). Suppose without loss of generality it’s the first: \(\langle\mathbf{v}_1,\mathbf{z}^{(q)}\rangle \geq 1-\epsilon\). Then we would conclude that \(\|\mathbf{v}_1 -\mathbf{z}^{(q)}\|_2 = 2 - 2\langle\mathbf{v}_1,\mathbf{z}^{(q)}\rangle \leq 2 - 2\cdot (1-\epsilon) = 2\epsilon\). \(\square\)

So, to show power methods converges, our goal is to prove that \(\left|c_j^{(q)}/c_1^{(q)}\right| \leq \sqrt{\epsilon/d}\) for all \(j\neq 1\).

Heart of the Analysis

To prove this bound, we note that \(c_j^{(q)} = S\cdot \sigma_j^{2q} c_j^{(0)}\) for any \(j\) and \(c_1^{(q)} = S\cdot \sigma_1^{2q} c_1^{(0)}\) where \(S = 1/\prod_{i=1}^q n_i\) is some fixed scaling. So: \[ \left|\frac{c_j^{(q)} }{c_1^{(q)} }\right| = (\sigma_j/\sigma_1)^{2q} \left|\frac{c_j^{(0)} }{c_1^{(0)} }\right|. \] As claimed in class (and can be prove as an exercise – try using rotational invariance of the Gaussian ) if \(z^{(0)}\) is a randomly initialized start vector, \(\left|\frac{c_j^{(0)} }{c_1^{(0)} }\right| \leq d^{3}\) with high probability. That may seem like a lot, but \((\sigma_j/\sigma_1)^{2q}\) is going to be a tiny number, so will easily cancel that out. In particular, \[ (\sigma_j/\sigma_1)^{2q} = \left(1-\frac{\sigma_1-\sigma_j}{\sigma_1}\right)^{2q} \leq (1-\gamma)^{2q}, \] where \(\gamma = \frac{\sigma_1-\sigma_2}{\sigma_1}\) is our spectral gap parameter. As long as we set \[ q = \frac{\log(d^3\sqrt{d/\epsilon})}{\gamma} = O\left(\frac{\log(d/\epsilon)}{\gamma}\right), \] then \((1-\gamma)^{2q} \leq \frac{\sqrt{\epsilon/d}}{d^3}\), and thus \[ \left|\frac{c_j^{(q)} }{c_1^{(q)} }\right| = (\sigma_j/\sigma_1)^{2q} \left|\frac{c_j^{(0)} }{c_1^{(0)} }\right| \leq d^3\cdot \frac{\sqrt{\epsilon/d}}{d^3} \leq \sqrt{\epsilon/d}, \] as desired.

Alternative Guarantee

In machine learning applications, we care less about actually approximating \(\mathbf{v}_1\), and more that the \(\mathbf{z}\) is a “good top singular vector” in that it offers a good rank-1 approximation to \(\mathbf{X}\). Tha tis, we want to prove that \[\|\mathbf{X} - \mathbf{X}\mathbf{z}\mathbf{z}^\top\|_F^2 = \|\mathbf{X}\|_F^2 - \|\mathbf{X}\mathbf{z}\mathbf{z}^\top\|_F^2\] is small, or equivalently, that \[\|\mathbf{X}\mathbf{z}\mathbf{z}^\top\|_F^2\] is large. The largest it could be is \(\|\mathbf{X}\mathbf{v}_1\mathbf{v}_1^\top\|_F^2 = \sigma_1^2\). And from the same argument above, where we claimed that either \(\langle\mathbf{v}_1,\mathbf{z}^{(q)}\rangle \geq 1-\epsilon\) or \(\langle- \mathbf{v}_1,\mathbf{z}^{(q)}\rangle \geq 1-\epsilon\), it is not hard to check that \[\|\mathbf{X}\mathbf{z}\mathbf{z}^\top\|_F^2 \geq (1-\epsilon)^2 \sigma_1^2,\] so after \(O\left(\frac{\log(d/\epsilon)}{\gamma}\right)\) iterations we get a near optimal low-rank approximation.

The Lanczos Method

We will now see how to improve on power method using what is known as the Lanczos method. Like power method, Lanczos is considered a Krylov subspace method because it will return a solution in the span of the Krylov subspace that we introduced before. Power method clearly does this – it returns a scaling of the last column of \(K\). The whole idea behind Lanczos is to avoid “throwing away’” information from earlier columns like power method, but instead to take advantage of the whole space. It turns at that doing so can be very helpful – we will get a bound that depends on \(1/\sqrt{\gamma}\) instead of \(1/\gamma\).

Specifically, to define the Lanczos method, we will let \(\mathbf{Q}\in \mathbb{R}^{d\times k}\) be a matrix with orthonormal columns that spans \(K\). In practice, you need to be careful about how this is computed for numerical stability reasons, but we won’t worry about that for now. Imagine your computer has infinite precision and we just compute \(\mathbf{Q}\) by orthonormalizing \(K\).

Lanczos Method

  • Compute an orthonormal span \(\mathbf{Q}\) for the degree \(q\) Krylov subspace.
  • Let \(\mathbf{z}\) be the top eigenvector of \(\mathbf{Q}^\top \mathbf{A}\mathbf{Q} = \mathbf{Q}^\top \mathbf{X}^\top \mathbf{X} \mathbf{Q}\)
  • Return \(\mathbf{Q} \mathbf{z}\).

Importantly, the first step only requires \(q\) matrix-vector multiplications with \(\mathbf{X}^\top\mathbf{X}\), each of which can be implemented in \(O(nd)\) time, just like we did for power method. The second step might look a bit circular at first glance. We want an approximation algorithm for computing the top eigenvector of \({\mathbf{X}}^\top\mathbf{X}\) and the above method uses a top eigenvector algorithm as a subroutine. But note that \(\mathbf{Q}^\top {\mathbf{AQ}}\) only has size \(q\times q\), where \(q\ll d\) is our iteration count. So even if it’s too expensive to compute a direct eigendecomposition of \(\mathbf{X}^\top\mathbf{X}\), \({\mathbf{Q}}^\top{\mathbf{X}}^\top{\mathbf{X}}{\mathbf{Q}}\) can be computed in \(O(ndq) + O(dq^2) + O(q^3) + O(ndq)\) time (the first part is the time to construct \(\mathbf{Q}\) and \(\mathbf{X}\mathbf{Q}\), the second to compute \({\mathbf{Q}}^\top{\mathbf{X}}^\top{\mathbf{X}}{\mathbf{Q}}\), and the third to find its top eigenvector using a direct method).

Analysis Preliminaries

Our first claim is that Lanczos returns the best approximate singular vector in the span of the Krylov subspace. Then we will argue that there always exists some vector in the span of the subspace that is significantly better than what power method returns, so the Lanczos solution must be significantly better as well.

Claim 2: Amongst all vectors in the span of the Krylov subspace (i.e., any vector \(\mathbf{y}\) that can be written as \(\mathbf{y} = \mathbf{Q}\mathbf{x}\) for some \(x\in \mathbb{R}^k\)), \(\mathbf{y}^* = \mathbf{Q}\mathbf{z}\) minimizes the low-rank approximation error \(\|\mathbf{X} - \mathbf{X}\mathbf{y}\mathbf{y}^\top\|_F^2\).

Proof. First, we can prove that \(\mathbf{y}\) should always be chosen to have unit norm, so that \(\mathbf{y}\mathbf{y}^\top\) is a projection. Accordingly, it must also be that \(\mathbf{x} = \mathbf{Q}^\top\mathbf{y}\) has unit norm. Then, proving the claim above is equivalent to proving that \(\mathbf{y}^* = {\mathbf{Q}}\mathbf{z}\) maximizes \(\|\mathbf{X}\mathbf{y}\mathbf{y}^\top\|_F^2 = \|\mathbf{X}\mathbf{y}\|_F^2 = \|\mathbf{X}\mathbf{y}\|_2^2 = \|\mathbf{X}\mathbf{Q}\mathbf{x}\|_2^2\). It is not hard to check that setting \(\mathbf{x}\) to be the maximum right singular vector of \({\mathbf{X}}{\mathbf{Q}}\) maximizes this \(\|\mathbf{X}\mathbf{Q}\mathbf{x}\|_2^2\) amongst all unit vectors, and this right singular vector is the same at the top eigenvector of \(\mathbf{Q}^\top \mathbf{X}^\top \mathbf{X} \mathbf{Q}\).

Claim 3: If we run the Lanczos method for \(O\left(\frac{\log(d/\epsilon)}{\sqrt{\gamma}}\right)\) There is some vector \(\mathbf{w}\) of the form \(\mathbf{w} = \mathbf{Q}\mathbf{x}\) such that either \(\langle\mathbf{v}_1,\mathbf{w}\rangle \geq 1-\epsilon\) or \(\langle- \mathbf{v}_1,\mathbf{w}\rangle \geq 1-\epsilon\).

In other words, there is some \(\mathbf{w}\) in the Krylov subspace that has a large inner product with the true top singular vector \(\mathbf{v}_1\). As seen earlier for power method, this can be used to prove that e.g. \(\|\mathbf{X}\mathbf{w}\mathbf{w}^\top\|_F^2\) is large, and from Claim 2, we know that the \(\mathbf{v}\) returned by Lanczos can only do better. So, we focus on proving Claim 3.

Heart of the Analysis

The key idea is to observe that any vector in the span of the Krylov subspace of size \(q\) – i.e. any vector \(\mathbf{w}\) that can be written \(\mathbf{w} = \mathbf{Q}\mathbf{x}\) is equal to: \[ \mathbf{w} = p(\mathbf{A})\mathbf{z}^{(0)}, \] for some degree q polynomial \(p\). For example, we might have \(p(\mathbf{A}) = 2\mathbf{\mathbf{A}}^2 - 4\mathbf{A}^3 + \mathbf{A}^6\) or \(p(\mathbf{A}) = I - \mathbf{A} - 10\mathbf{A}^5\). And moreover, for any degree \(q\) polynomial \(p\), there is some \(\mathbf{x}\) such that \(\mathbf{Q}\mathbf{x} = p(\mathbf{A})\mathbf{z}^{(0)}\). To see that this is true, we note that all of the monomials \(\mathbf{z}^{(0)}, A\mathbf{z}^{(0)}, \ldots, \mathbf{A}^q\mathbf{z}^{(0)}\) lie in the span of \(\mathbf{Q}\), so any linear combination does as well.

This means that showing there is a good approximate top singular vector \(\mathbf{w}\) in the span of the Krylov subspace can be reduced to finding a good polynomial \(p(\mathbf{A})\). Specifically, notice that, if we write \(p(\mathbf{A})\mathbf{z}^{(0)}\) in the span of the singular vectors we have: \[ p(\mathbf{A})\mathbf{z}^{(0)} = g_1 \mathbf{v}_1 + g_2 \mathbf{v}_2 + \ldots + g_d \mathbf{v}_d \] where \[ g_j = c_j^{(0)}p(\sigma_j^2). \] Our goal, which is exactly the same as when we analyzed power method, is to show that \(g_1\) is much larger than \(g_j\) for all \(j \neq 1\). In other words we want to find a polynomial such that \(p(t)\) is small for all values of \(0\leq t< \sigma_1^2\), but then the polynomial should suddenly jump to be large at \(\sigma_1\). The simplest degree \(q\) polynomial that does this is \(p(t) = t^q\). However, it turns out there are more sharply jumping polynomials, which can be obtained by shifting/scaling a type of polynomials known as a Chebyshev polynomial. An example that jumps at \(\sigma_1^2 = 1\) is shown below. The difference from \(t^q\) is pretty remarkable – even though the polynomial \(p\) shown below is nearly as small for all \(0\leq t< 1\), it is much larger at \(t=1\).

Concretely we can claim the following, which is a bit tricky to prove, but well known (see Lemma 5 here for a full proof).

Claim 4: There is a degree \(O\left(\sqrt{\frac{1}{\gamma}}\log\frac{1}{\epsilon}\right)\) degree polynomial \(\hat{p}\) such that \(\hat{p}(1) = 1\) and \(|p(t)| \leq \epsilon\) for \(0 \leq t \leq 1-\gamma\).

In contrast, to ensure that \(t^q\) satisfies the above, we would need degree q = \(O\left({\frac{1}{\gamma}\log\frac{1}{\epsilon}}\right)\) – a quadratically worse bound. This is what will account for the quadratic difference in performance between the Lanczos and Power Methods.

The fact that there exist much steeper low-degree polynomails than \(t^q\) is the sort of thing studied in the mathematical field known as Approximation Theory. That might seem pretty obscure, but steep polynomials are surprisingly useful in computer science, appearing everywhere from classic, to Quantum complexity theory, to learning theory. If you’re interested in learning more about this, check out slides for this talk.

Finishing Up

Finally, we use Claim 4 to finish the analysis of Lanczos. Consider the vector \(\hat{p}\left(\frac{1}{\sigma_1^2}A\right)\mathbf{z}^{(0)}\), which as argued above lies in the Krylov subspace. As discussed before, our job is to prove that: \[ \frac{c_j^{(0)}\hat{p}\left(\frac{1}{\sigma_1^2}\sigma_j^2\right)}{c_1^{(0)}\hat{p}\left(\frac{1}{\sigma_1^2}\sigma_1^2\right)} \leq \sqrt{\epsilon/d}, \] for all \(j \neq i\), as long as \(\hat{p}\) is chosen to have degree \(q = O\left(\sqrt{\frac{1}{\gamma}}\log\frac{d}{\epsilon}\right)\)

Consider the numerator first. \(\frac{\sigma_j^2}{\sigma_1^2} = \left(\frac{\sigma_j}{\sigma_1}\right)^2 = \left(1 - \left(1-\frac{\sigma_j}{\sigma_1}\right)\right)^2 = \left(1 - \frac{\sigma_1 - \sigma_j}{\sigma_1}\right)^2 \leq (1-\gamma)^2 \leq (1-\gamma)\) .

Accordingly, if we set \(q = O\left(\sqrt{\frac{1}{\gamma}}\log\frac{1}{\epsilon'}\right)\) where \(\epsilon' = \sqrt{\epsilon/d}/d^3\), then \(\hat{p}\left(\frac{1}{\sigma_1^2}\sigma_j^2\right) \leq \epsilon'\). And the denomenator is simply equal to \(c_1^{(0)}\hat{p}(1) = c_1^{(0)}\cdot 1\). So, since \(c_j^{(0)}/c_1^{(0)}\leq d^3\) with high probability, as argued earlier, the equation holds. And thus Claim 3 is proven by setting \(\mathbf{w} = \hat{p}\left(\frac{1}{\sigma_1^2} \mathbf{A} \right)\mathbf{z}^{(0)}\).