\(\newcommand{\horz}{\rule[.5ex]{2.5ex}{0.5pt}}\)
3.4. Power iteration#
There is in general no exact method for computing SVDs. Instead we must rely on iterative methods, that is, methods that progressively approach the solution. We describe in this section the power iteration method. This method is behind an effective numerical approach for computing SVDs.
The focus here is on numerical methods and we will not spend much time computing SVDs by hand. But note that the connection between the SVD and the spectral decomposition of \(A^T A\) can be used for this purpose on small examples.
3.4.1. Key lemma#
We now derive the main idea behind an algorithm to compute singular vectors.
Eigenvectors in semidefinite case We start with eigenvectors. For our purposes, it will suffice to (1) restrict ourselves to positive semidefinite matrices and (2) focus on the largest eigenvalue.
Let \(A\) be a symmetric, positive semindefinite matrix in \(\mathbb{R}^{d \times d}\). By the Spectral Theorem, it has an eigenvector decomposition
where further \(0 \leq \lambda_d \leq \cdots \leq \lambda_1\) by the Characterization of Positive Semidefiniteness. To see how one might go about computing it, we start with the following observation.
Because of the orthogonality of \(Q\), the powers of \(A\) have a simple representation. The square gives
Repeating, we obtain
Further
When \(\lambda_1 > \lambda_2, \ldots, \lambda_d\), we get that \(\lambda_1^{k} \gg \lambda_2^{k}, \ldots, \lambda_d^{k}\) when \(k\) is large. Then, we get the approximation
This leads to the following idea to compute \(\mathbf{q}_1\):
LEMMA (Power Iteration in the Positive Semidefinite Case) Let \(A\) be a symmetric, positive semindefinite matrix in \(\mathbb{R}^{d \times d}\) with eigenvector decomposition \(A= Q \Lambda Q^T\) where the eigenvalues satisfy \(0 \leq \lambda_d \leq \cdots \leq \lambda_2 < \lambda_1\). Assume that \(\mathbf{x} \in \mathbb{R}^d\) is a vector such that \(\langle \mathbf{q}_1, \mathbf{x} \rangle > 0\). Then
as \(k \to +\infty\). If instead \(\langle \mathbf{q}_1, \mathbf{x} \rangle < 0\), then the limit is \(- \mathbf{q}_1\). \(\flat\)
The proof is similar to the case of singular vectors, which is detailed next.
SVD Let \(U \Sigma V^T\) be a (compact) SVD of \(A\). To see how one might go about computing it, we start with the following observation which relies on a different way to state the connection between the SVD and the spectral decomposition of \(A^T A\). Because of the orthogonality of \(U\) and \(V\), the powers of \(A^T A\) have a simple representation. Indeed
so that
Repeating
Further,
So
When \(\sigma_1 > \sigma_2, \ldots, \sigma_r\), which is often the case with real datasets, we get that \(\sigma_1^{2k} \gg \sigma_2^{2k}, \ldots, \sigma_r^{2k}\) when \(k\) is large. Then, we get the approximation
Finally, we arrive at:
LEMMA (Power Iteration) Let \(A \in \mathbb{R}^{n\times m}\) be a matrix and let \(U \Sigma V^T\) be a (compact) SVD of \(A\) such that \(\sigma_1 > \sigma_2 > 0\). Define \(B = A^T A\) and assume that \(\mathbf{x} \in \mathbb{R}^m\) is a vector satisfying \(\langle \mathbf{v}_1, \mathbf{x} \rangle > 0\). Then
as \(k \to +\infty\). If instead \(\langle \mathbf{v}_1, \mathbf{x} \rangle < 0\), then the limit is \(- \mathbf{v}_1\). \(\flat\)
Proof idea: We use the approximation above and divide by the norm to get a unit norm vector in the direction of \(\mathbf{v}_1\).
Proof: We have
So
This goes to \(\mathbf{v}_1\) as \(k\to +\infty\) if the expression in the first curly brackets goes to \(1\) and the one in the second curly brackets goes to \(0\). We prove this in the next claim.
LEMMA As \(k\to +\infty\),
\(\flat\)
Proof: Because the \(\mathbf{v}_j\)s are an orthonormal basis,
So, as \(k\to +\infty\), using the fact that \(\langle \mathbf{v}_1, \mathbf{x} \rangle \neq 0\) by assumption
since \(\sigma_j < \sigma_1\) for all \(j =2,\ldots,r\). That implies the first part of the claim by taking a square root and using \(\langle \mathbf{v}_1, \mathbf{x} \rangle > 0\). The second part of the claim follows essentially from the same argument. \(\square\) \(\square\)
EXAMPLE: We revisit the example
We previously compute its SVD and found that
This time we use the Power Iteration Lemma. Here
Taking powers of this matrix is easy
Let’s choose an arbitrary initial vector \(\mathbf{x}\), say \((-1, 2)^T\). Then
So
as \(k \to +\infty\). In fact, in this case, convergence occurs after one step. \(\lhd\)
3.4.2. Computing the top singular vector#
Power iteration gives us a way to compute \(\mathbf{v}_1\) – at least approximately if we use a large enough \(k\). But how do we find an appropriate vector \(\mathbf{x}\), as required by the Power Iteration Lemma? It turns out that a random vector will do. For instance, let \(\mathbf{X}\) be an \(m\)-dimensional spherical Gaussian with mean \(0\) and variance \(1\). Then, \(\mathbb{P}[\langle \mathbf{v}_1, \mathbf{X} \rangle = 0] = 0\). We will show this when we discuss multivariate Gaussians. Note that, if \(\langle \mathbf{v}_1, \mathbf{X} \rangle < 0\), we will instead converge to \(-\mathbf{v}_1\) which is also a right singular vector.
NUMERICAL CORNER: We implement the algorithm suggested by the Power Iteration Lemma. That is, we compute \(B^{k} \mathbf{x}\), then normalize it. To obtain the corresponding singular value and left singular vector, we use that \(\sigma_1 = \|A \mathbf{v}_1\|\) and \(\mathbf{u}_1 = A \mathbf{v}_1/\sigma_1\).
def topsing(A, maxiter=10):
x = rng.normal(0,1,np.shape(A)[1])
B = A.T @ A
for _ in range(maxiter):
x = B @ x
v = x / LA.norm(x)
s = LA.norm(A @ v)
u = A @ v / s
return u, s, v
We will apply it to our previous two-cluster example. The necessary functions are in mmids.py, which is available on the GitHub of the notes.
d, n, w = 10, 100, 3.
X = mmids.two_mixed_clusters(d, n, w)
Show code cell source
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(X[:,0], X[:,1], marker='x')
plt.show()
Let’s compute the top singular vector.
u, s, v = topsing(X)
print(v)
[-0.99403346 0.06169445 0.04646685 -0.00370871 0.02909716 -0.02980242
-0.02583767 0.04676195 -0.03515519 -0.00966285]
This is approximately \(-\mathbf{e}_1\). We get roughly the same answer (possibly up to sign) from Python’s numpy.linalg.svd function.
u, s, vh = LA.svd(X)
print(vh.T[:,0])
[ 0.99403346 -0.06169444 -0.04646684 0.00370871 -0.02909717 0.02980244
0.02583767 -0.04676197 0.03515518 0.00966286]
Recall that, when we applied \(k\)-means clustering to this example with \(d=1000\) dimension, we obtained a very poor clustering.
d, n, w = 1000, 100, 3.
X = mmids.two_mixed_clusters(d, n, w)
Show code cell source
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(X[:,0], X[:,1], marker='x')
plt.show()
assign = mmids.kmeans(X, 2)
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
99641.677051717
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(X[:,0], X[:,1], c=assign)
plt.show()
Let’s try again, but after projecting on the top singular vector. Recall that this corresponds to finding the best one-dimensional approximating subspace. The projection can be computed using the truncated SVD \(Z= U_{(1)} \Sigma_{(1)} V_{(1)}^T\). We can interpret the rows of \(U_{(1)} \Sigma_{(1)}\) as the coefficients of each data point in the basis \(\mathbf{v}_1\). We will work in that basis. We need one small hack: because our implementation of \(k\)-means clustering expects data points in at least \(2\) dimension, we add a column of \(0\)’s.
u, s, v = topsing(X)
Xproj = np.stack((u*s, np.zeros(np.shape(X)[0])), axis=-1)
Show code cell source
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(Xproj[:,0], Xproj[:,1], alpha=0.25)
plt.ylim([-3,3])
plt.show()
There is a small – yet noticeable – gap around 0.
A histogram of the first component of Xproj gives a better sense of the density of points.
Show code cell source
plt.hist(Xproj[:,0])
plt.show()
We run \(k\)-means clustering on the projected data.
assign = mmids.kmeans(Xproj, 2)
1309.1441138296427
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
386.8878320256463
Show code cell source
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(X[:,0], X[:,1], c=assign)
plt.show()
Much better. We give a more formal explanation of this outcome in an optional section. In essence, quoting [BHK, Section 7.5.1]:
[…] let’s understand the central advantage of doing the projection to [the top \(k\) right singular vectors]. It is simply that for any reasonable (unknown) clustering of data points, the projection brings data points closer to their cluster centers.
Finally, looking at the top right singular vector (or its first ten entries for lack of space), we see that it does align quite well (but not perfectly) with the first dimension.
print(v[:10])
[ 0.62787593 0.02260667 -0.0020277 0.01371159 0.01163581 0.01610368
0.00903951 -0.01153729 -0.00906383 -0.00990353]
\(\unlhd\)
3.4.3. Computing more singular vectors#
We have shown how to compute the first singular vector. How do we compute more singular vectors? One approach is to first compute \(\mathbf{v}_1\) (or \(-\mathbf{v}_1\)), then find a vector \(\mathbf{y}\) orthogonal to it, and proceed as above. And then we repeat until we have all \(m\) right singular vectors.
We are often interested only in the top, say \(\ell < m\), singular vectors. An alternative approach in that case is to start with \(\ell\) random vectors and, first, find an orthonormal basis for the space they span. Then to quote [BHK, Section 3.7.1]:
Then compute \(B\) times each of the basis vectors, and find an orthonormal basis for the space spanned by the resulting vectors. Intuitively, one has applied \(B\) to a subspace rather than a single vector. One repeatedly applies \(B\) to the subspace, calculating an orthonormal basis after each application to prevent the subspace collapsing to the one dimensional subspace spanned by the first singular vector.
We will not prove here that this approach, known as orthogonal iteration, works. The proof is similar to that of the Power Iteration Lemma.
NUMERICAL CORNER: We implement this last algorithm. We will need our previous implementation of Gram-Schimdt.
def svd(A, l, maxiter=100):
V = rng.normal(0,1,(np.size(A,1),l))
for _ in range(maxiter):
W = A @ V
Z = A.T @ W
V, R = mmids.gramschmidt(Z)
W = A @ V
S = [LA.norm(W[:, i]) for i in range(np.size(W,1))]
U = np.stack([W[:,i]/S[i] for i in range(np.size(W,1))],axis=-1)
return U, S, V
Note that above we avoided forming the matrix \(A^T A\). With a small number of iterations, that approach potentially requires fewer arithmetic operations overall and it allows to take advantage of the possible sparsity of \(A\) (i.e. the fact that it may have many zeros).
We apply it again to our two-cluster example.
d, n, w = 1000, 100, 3.
X = mmids.two_mixed_clusters(d, n, w)
Let’s try again, but after projecting on the top two singular vectors. Recall that this corresponds to finding the best two-dimensional approximating subspace. The projection can be computed using the truncated SVD \(Z= U_{(2)} \Sigma_{(2)} V_{(2)}^T\). We can interpret the rows of \(U_{(2)} \Sigma_{(2)}\) as the coefficients of each data point in the basis \(\mathbf{v}_1,\mathbf{v}_2\). We will work in that basis.
U, S, V = svd(X, 2)
Xproj = np.stack((U[:,0]*S[0], U[:,1]*S[1]), axis=-1)
assign = mmids.kmeans(Xproj, 2)
3329.139169940858
2525.8080945260817
2513.7186850891817
2499.8370677155963
2498.548904664579
2498.548904664579
2498.548904664579
2498.548904664579
2498.548904664579
2498.548904664579
Show code cell source
fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.scatter(X[:,0], X[:,1], c=assign)
plt.show()
Finally, looking at the first two right singular vectors, we see that the first one does align quite well with the first dimension.
print(np.stack((V[:,0], V[:,1]), axis=-1))
[[ 0.65108828 0.04380735]
[-0.00804941 0.032942 ]
[-0.03645334 -0.00097783]
...
[-0.00579195 -0.01518118]
[ 0.01448518 -0.02433161]
[-0.05940409 -0.01091726]]
\(\unlhd\)