# Spectral and singular value decompositions¶

## 3 Singular value decomposition¶

Course: Math 535 - Mathematical Methods in Data Science (MMiDS)
Author: Sebastien Roch, Department of Mathematics, University of Wisconsin-Madison
Updated: Sep 21, 2020

### 3.1 Definition¶

First recall a matrix $D \in \mathbb{R}^{n \times m}$ is diagonal if its non-diagonal entries are zero. That is, $i \neq j$ implies that $D_{ij} =0$. We now come to our main definition.

Definition (Singular Value Decomposition): Let $A \in \mathbb{R}^{n\times m}$ be a matrix with $m \leq n$. A singular value decomposition (SVD) of $A$ is a matrix factorization

$$A = U \Sigma V^T = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$$

where the columns of $U \in \mathbb{R}^{n \times r}$ and those of $V \in \mathbb{R}^{m \times r}$ are orthonormal, and $\Sigma \in \mathbb{R}^{r \times r}$ is a diagonal matrix. Here the $\mathbf{u}_j$'s are the columns of $U$ and are referred to as left singular vectors. Similarly the $\mathbf{v}_j$'s are the columns of $V$ and are referred to as right singular vectors. The $\sigma_j$'s, which are non-negative and in decreasing order

$$\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0$$

are the diagonal elements of $\Sigma$ and are referred to as singular values. $\lhd$

Exercise: Let $A = U \Sigma V^T = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$ be an SVD of $A$.

(a) Show that $\{\mathbf{u}_j:j\in[r]\}$ is a basis of $\mathrm{col}(A)$ and that $\{\mathbf{v}_j:j\in[r]\}$ is a basis of $\mathrm{row}(A)$.

(b) Show that $r$ is the rank of $A$. $\lhd$

This type of matrix factorization is illustrated below (in its full form).

(Source)

Remarkably, any matrix has an SVD.

Theorem (Existence of SVD): Any matrix $A \in \mathbb{R}^{n\times m}$ has a singular value decomposition.

The construction works as follows. Compute the greedy sequence $\mathbf{v}_1,\ldots,\mathbf{v}_r$ from the previous section until the first $r$ such that

$$\max \{\|A \mathbf{v}\|:\|\mathbf{v}\| = 1, \ \langle \mathbf{v}, \mathbf{v}_j \rangle = 0, \forall j \leq r\} = 0$$

or, otherwise, until $r=m$. The $\mathbf{v}_j$'s are orthonormal by construction. For $j=1,\ldots,m$, let

$$\sigma_j = \|A \mathbf{v}_j\| \quad\text{and}\quad \mathbf{u}_j = \frac{1}{\sigma_j} A \mathbf{v}_j.$$

Observe that, by our choice of $r$, the $\sigma_j$'s are $> 0$. In particular, the $\mathbf{u}_j$'s have unit norm by definition. We show next that they are also orthogonal.

Lemma (Left Singular Vectors are Orthogonal): For all $1 \leq i \neq j \leq r$, $\langle \mathbf{u}_i, \mathbf{u}_j \rangle = 0$.

Proof idea: Quoting [BHK, Section 3.6]:

Intuitively if $\mathbf{u}_i$ and $\mathbf{u}_j$, $i < j$, were not orthogonal, one would suspect that the right singular vector $\mathbf{v}_j$ had a component of $\mathbf{v}_i$ which would contradict that $\mathbf{v}_i$ and $\mathbf{v}_j$ were orthogonal. Let $i$ be the smallest integer such that $\mathbf{u}_i$ is not orthogonal to all other $\mathbf{u}_j$. Then to prove that $\mathbf{u}_i$ and $\mathbf{u}_j$ are orthogonal, we add a small component of $\mathbf{v}_j$ to $\mathbf{v}_i$, normalize the result to be a unit vector $\mathbf{v}'_i \propto \mathbf{v}_i + \epsilon \mathbf{v}_j$ and show that $\|A \mathbf{v}'_i\| > \|A \mathbf{v}_i\|$, a contradiction.

Proof idea (Existence of SVD): We show that $\sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$ acts on vectors in the same way that $A$ does, by an orthogonal decomposition over the span of the $\mathbf{v}_j$'s and its complement.

Proof (Existence of SVD): Let

$$C = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T.$$

We claim that $A = C$ as matrices. It suffices to show that the left-hand side and right-hand side are the same linear map, that is, $A \mathbf{v} = C \mathbf{v}$ for all $\mathbf{v} \in \mathbb{R}^m$. Indeed, one then has in particular that the $i$-th column in both cases is $A \mathbf{e}_i = C \mathbf{e}_i$ for all $i$.

Let $\mathbf{v} \in \mathbb{R}^m$ be any vector and let $\mathcal{Z} = \mathrm{span}(\mathbf{v}_1,\ldots,\mathbf{v}_r)$. Decompose $\mathbf{v}$ into orthogonal components

$$\mathbf{v} = \mathrm{proj}_{\mathcal{Z}}(\mathbf{v}) + (\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v})) = \sum_{j=1}^r \langle \mathbf{v}, \mathbf{v}_j\rangle \,\mathbf{v}_j + (\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v})).$$

Applying $A$ we get

\begin{align} A \mathbf{v} &= \sum_{j=1}^r \langle \mathbf{v}, \mathbf{v}_j\rangle \, A\mathbf{v}_j + A (\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v}))\\ &= \sum_{j=1}^r \langle \mathbf{v}, \mathbf{v}_j\rangle \, \sigma_j \mathbf{u}_j + A (\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v}))\\ &= \sum_{j=1}^r \langle \mathbf{v}, \mathbf{v}_j\rangle \, \sigma_j \mathbf{u}_j \end{align}

where we used that, by definition of $r$, $A \mathbf{w} = 0$ for any $\mathbf{w}$ orthogonal to $\mathcal{Z}$, which includes $\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v})$ by the Orthogonal Projection Theorem.

Similarly,

\begin{align} \left(\sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T\right)\mathbf{v} &= \left(\sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T\right) \left(\sum_{k=1}^r \langle \mathbf{v}, \mathbf{v}_k\rangle \, \mathbf{v}_k + (\mathbf{v} - \mathrm{proj}_{\mathcal{Z}}(\mathbf{v}))\right)\\ &= \sum_{j=1}^r \sum_{k=1}^r \,\langle \mathbf{v}, \mathbf{v}_k\rangle \, \sigma_j \mathbf{u}_j \,(\mathbf{v}_j^T \mathbf{v}_k)\\ &= \sum_{j=1}^r \langle \mathbf{v}, \mathbf{v}_j\rangle \,\sigma_j \mathbf{u}_j, \end{align}

where, again, we used the Orthogonal Decomposition Lemma on the second line.

That establishes the claim. $\square$

The SVD also has a natural geometric interpretation. To quote [Sol, p. 133]:

The SVD provides a complete geometric characterization of the action of $A$. Since $U$ and $V$ are orthogonal, they have no effect on lengths and angles; as a diagonal matrix, $\Sigma$ scales individual coordinate axes. Since the SVD always exists, all matrices $A \in \mathbb{R}^{n \times m}$ are a composition of an isometry, a scale in each coordinate, and a second isometry.

This sequence of operations is illustrated below.

(Source)

### 3.2 Power iteration¶

There is in general no exact method for computing SVDs. Instead we must rely on iterative methods, that is, methods that approach the solution. Let $U \Sigma V^T$ be an SVD of $A$. To see how one might go about computing such an SVD, we start with the following observation. Because of the orthogonality of $U$ and $V$, the powers of $A^T A$ have a simple representation. Indeed

$$B = A^T A = (U \Sigma V^T)^T (U \Sigma V^T) = V \Sigma^T U^T U \Sigma V^T = V \Sigma^T \Sigma V^T,$$

so that

$$B^2 = (V \Sigma^T \Sigma V^T) (V \Sigma^T \Sigma V^T) = V (\Sigma^T \Sigma)^2 V^T$$

and repeating

$$B^{k} = V (\Sigma^T \Sigma)^{k} V^T.$$

Further,

$$\widetilde{\Sigma} = \Sigma^T \Sigma = \begin{pmatrix} \sigma_1^2 & 0 & \cdots & 0\\ 0 & \sigma_2^2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \sigma_r^2 \end{pmatrix}$$

and

$$\widetilde{\Sigma}^k = \begin{pmatrix} \sigma_1^{2k} & 0 & \cdots & 0\\ 0 & \sigma_2^{2k} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \sigma_r^{2k} \end{pmatrix}.$$

When $\sigma_1 > \sigma_2$, 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. In that case, we get the approximation

$$B^{k} = \sum_{j=1}^r \sigma_j^{2k} \mathbf{v}_j \mathbf{v}_j^T. \approx \sigma_1^{2k} \mathbf{v}_1 \mathbf{v}_1^T.$$

This leads to the following idea to compute $\mathbf{v}_1$:

Lemma (Power Iteration): Let $A \in \mathbb{R}^{n\times m}$ be a matrix with $m \leq n$. Let $U \Sigma V^T$ be an SVD of $A$ such that $\sigma_1 > \sigma_2$. 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

$$\frac{B^{k} \mathbf{x}}{\|B^{k} \mathbf{x}\|} \to \mathbf{v}_1$$

as $k \to +\infty$.

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

$$B^{k}\mathbf{x} = \sum_{j=1}^r \sigma_j^{2k} \mathbf{v}_j \mathbf{v}_j^T \mathbf{x}$$

and, because the $\mathbf{v}_j$'s are an orthonormal basis,

$$\frac{1}{\sigma_1^{4k} (\mathbf{v}_1^T \mathbf{x})^2} \|B^{k}\mathbf{x}\|^2 = \frac{1}{\sigma_1^{4k} (\mathbf{v}_1^T \mathbf{x})^2} \sum_{j=1}^r \left|\sigma_j^{2k} (\mathbf{v}_j^T \mathbf{x})\right|^2 = 1 + \sum_{j=2}^r \frac{\sigma_j^{4k} (\mathbf{v}_j^T \mathbf{x})^2}{\sigma_1^{4k} (\mathbf{v}_1^T \mathbf{x})^2} \to 1$$

as $k\to +\infty$. Hence

$$\frac{B^{k} \mathbf{x}}{\|B^{k} \mathbf{x}\|} = \sum_{j=1}^r \mathbf{v}_j \frac{\sigma_j^{2k} (\mathbf{v}_j^T \mathbf{x})} {\|B^{k} \mathbf{x}\|} = \mathbf{v}_1 \frac{\sigma_1^{2k} (\mathbf{v}_1^T \mathbf{x})} {\|B^{k} \mathbf{x}\|} + \sum_{j=2}^r \mathbf{v}_j \frac{\sigma_j^{2k} (\mathbf{v}_j^T \mathbf{x})} {\|B^{k} \mathbf{x}\|} \to \mathbf{v}_1$$

as $k\to +\infty$, by the above and the assumption $\sigma_1 > \sigma_2$. $\square$

That gives us a way to compute $\mathbf{v}_1$ (approximately). How do we find an appropriate vector $\mathbf{x}$? 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$. 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.

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.

Exercise: Let $A$, $U \Sigma V^T$, $B$ be as in the Power Iteration Lemma. Assume further that $\sigma_2 > \sigma_3$ and that $\mathbf{y} \in \mathbb{R}^m$ satisfies both $\langle \mathbf{v}_1, \mathbf{y} \rangle = 0$ and $\langle \mathbf{v}_2, \mathbf{y} \rangle > 0$. Show that $\frac{B^{k} \mathbf{y}}{\|B^{k} \mathbf{y}\|} \to \mathbf{v}_2$ as $k \to +\infty$. How would you find such a $\mathbf{y}$? $\lhd$

But 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.

In [2]:
function mmids_svd(A, l; maxiter=100)
V = randn(size(A,2),l) # random initialization
for _ = 1:maxiter
W = A * V
Z = A' * W
V, R = mmids_gramschmidt(Z)
end
W = A * V
S = [norm(W[:, i]) for i=1:size(W,2)] # singular values
U = reduce(hcat,[W[:,i]/S[i] for i=1:size(W,2)]) # left singular vectors
return U, S, V
end

Out[2]:
mmids_svd (generic function with 1 method)
In [3]:
d, n, offset = 10, 100, 3
X1, X2 = two_clusters(d, n, offset)
X = vcat(X1, X2)
scatter(X[:,1], X[:,2], legend=false)

Out[3]:
In [4]:
U, S, V = mmids_svd(X, 1)
V

Out[4]:
10×1 Array{Float64,2}:
0.9962548680669375
0.03175449922317592
-0.046035453904927505
0.025825801595636282
-0.04449037015490855
-0.029726316647995778
-0.011654005208952954
-0.015931727206384222
0.015283439352855857
-0.013978138721069698
In [6]:
d, n, offset = 1000, 100, 3
X1, X2 = two_clusters(d, n, offset)
X = vcat(X1, X2)
assign = mmids_kmeans(X, 2);

G = sum(dist .^ 2) = 200285.18553307335
G = sum(dist .^ 2) = 200226.174659831
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306
G = sum(dist .^ 2) = 200206.05494999306

In [7]:
scatter(X[:,1], X[:,2], marker_z=assign, legend=false)

Out[7]:
In [8]:
U, S, V = mmids_svd(X, 2)
Xproj = hcat(U[:,1].*S[1], U[:,2].*S[2])
assign = mmids_kmeans(Xproj, 2);

G = sum(dist .^ 2) = 4242.179608092223
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775
G = sum(dist .^ 2) = 2467.0066479604775

In [9]:
scatter(X[:,1], X[:,2], marker_z=assign, legend=false)

Out[9]:

### 3.3 Low-rank approximation in the induced norm¶

Let $A \in \mathbb{R}^{n \times m}$ be a matrix with SVD

$$A = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T.$$

For $k < r$, truncate the sum at the $k$-th term

$$A_k = \sum_{j=1}^k \sigma_j \mathbf{u}_j \mathbf{v}_j^T.$$

The rank of $A_k$ is exactly $k$. Indeed, by construction,

1. the vectors $\{\mathbf{u}_j\,:\,j = 1,\ldots,k\}$ are orthonormal, and

2. since $\sigma_j > 0$ for $j=1,\ldots,k$ and the vectors $\{\mathbf{v}_j\,:\,j = 1,\ldots,k\}$ are orthonormal, $\{\mathbf{u}_j\,:\,j = 1,\ldots,k\}$ spans the column space of $A_k$.

We have shown before that $A_k$ is the best approximation to $A$ among matrices of rank at most $k$ in Frobenius norm. Specifically, the Greedy Finds Best Fit Theorem implies that, for any matrix $B \in \mathbb{R}^{n \times m}$ of rank at most $k$,

$$\|A - A_k\|_F \leq \|A - B\|_F.$$

We show in this section that the same holds in the induced norm. First, some observations.

Lemma (Matrix Norms and Singular Values): Let $A \in \mathbb{R}^{n \times m}$ be a matrix with SVD

$$A = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$$

where recall that $\sigma_1 \geq \sigma_2 \geq \cdots \sigma_r > 0$ and let $A_k$ be the truncation defined above. Then

$$\|A - A_k\|^2_F = \sum_{j=k+1}^r \sigma_j^2$$

and

$$\|A - A_k\|^2_2 = \sigma_{k+1}^2.$$

Proof: For the first claim, by definition, summing over the columns of $A - A_k$

$$\|A - A_k\|^2_F = \left\|\sum_{j=k+1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T\right\|_F^2 = \sum_{i=1}^m \left\|\sum_{j=k+1}^r \sigma_j v_{j,i} \mathbf{u}_j \right\|^2.$$

Because the $\mathbf{u}_j$'s are orthonormal, this is

$$\sum_{i=1}^m \sum_{j=k+1}^r \sigma_j^2 v_{j,i}^2 = \sum_{j=k+1}^r \sigma_j^2 \left(\sum_{i=1}^m v_{j,i}^2\right) = \sum_{j=k+1}^r \sigma_j^2$$

where we used that the $\mathbf{v}_j$'s are also orthonormal.

For the second claim, recall that the induced norm is defined as

$$\|B\|_2 = \max_{\mathbf{x} \in \mathbb{S}^{m-1}} \|B \mathbf{x}\|.$$

For any $\mathbf{x} \in \mathbb{S}^{m-1}$

$$\left\|(A - A_k)\mathbf{x}\right\|^2 = \left\| \sum_{j=k+1}^r \sigma_j \mathbf{u}_j (\mathbf{v}_j^T \mathbf{x}) \right\|^2 = \sum_{j=k+1}^r \sigma_j^2 \langle \mathbf{v}_j, \mathbf{x}\rangle^2.$$

Because the $\sigma_j$'s are in decreasing order, this is maximized when $\langle \mathbf{v}_j, \mathbf{x}\rangle = 1$ if $j=k+1$ and $0$ otherwise. That is, we take $\mathbf{x} = \mathbf{v}_{k+1}$ and the norm is then $\sigma_{k+1}^2$, as claimed. $\square$

Theorem (Low-Rank Approximation in the Induced Norm): Let $A \in \mathbb{R}^{n \times m}$ be a matrix with SVD

$$A = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$$

and let $A_k$ be the truncation defined above with $k < r$. For any matrix $B \in \mathbb{R}^{n \times m}$ of rank at most $k$,

$$\|A - A_k\|_2 \leq \|A - B\|_2.$$

Proof idea: We know that $\|A - A_k\|_2^2 = \sigma_{k+1}^2$. So we want to lower bound $\|A - B\|_2^2$ by $\sigma_{k+1}^2$. For that, we have to find an appropriate $\mathbf{z}$ for any given $B$ of rank at most $k$. The idea is to take a vector $\mathbf{z}$ in the intersection of the null space of $B$ and the span of the singular vectors $\mathbf{v}_1,\ldots,\mathbf{v}_{k+1}$. By the former, the squared norm of $(A - B)\mathbf{z}$ is equal to the squared norm of $A\mathbf{z}$ which lower bounds $\|A\|_2^2$. By the latter, $\|A \mathbf{z}\|^2$ is at least $\sigma_{k+1}^2$.

### 3.4 Why project?¶

We return to $k$-means clustering and why projecting to a lower-dimensional subspace can produce better results. We prove a simple inequality that provides some insight. 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.

To elaborate, suppose we have $n$ data points in $d$ dimension in the form of the rows $\mathbf{a}_i^T$, $i=1\ldots, n$, of matrix $A \in \mathbb{A}^{n \times d}$, where we assume that $n > d$ and that $A$ has full column rank. Imagine these data points come from an unknown ground-truth $k$-clustering assignment $g(i) \in [k]$, $i = 1,\ldots, n$, with corresponding unknown centers $\mathbf{c}_j$, $j = 1,\ldots, k$, for $g(i) = j$. Let $C \in \mathbb{R}^{n \times d}$ be the corresponding matrix, that is, row $i$ of $C$ is $\mathbf{c}_j^T$ if $g(i) = j$. The $k$-means objective of the true clustering is then

\begin{align} \sum_{j\in [k]} \sum_{i:g(i)=j} \|\mathbf{a}_i - \mathbf{c}_{j}\|^2 &= \sum_{j\in [k]} \sum_{i:g(i)=j} \sum_{\ell=1}^d (a_{i,\ell} - c_{j,\ell})^2\\ &= \sum_{i=1}^n \sum_{\ell=1}^d (a_{i,\ell} - c_{g(i),\ell})^2\\ &= \|A - C\|_F^2. \end{align}

The matrix $A$ has an SVD

$$A = \sum_{j=1}^r \sigma_j \mathbf{u}_j \mathbf{v}_j^T$$

and for $k < r$ we have the truncation

$$A_k = \sum_{j=1}^k \sigma_j \mathbf{u}_j \mathbf{v}_j^T.$$

It corresponds to projecting each row of $A$ onto the linear subspace spanned by the first $k$ right singular vectors $\mathbf{v}_1, \ldots, \mathbf{v}_k$. To see this, note that the $i$-th row of $A$ is $\boldsymbol{\alpha}_i^T = \sum_{j=1}^r \sigma_j u_{j,i} \mathbf{v}_j^T$ and that, because the $\mathbf{v}_j$'s are linearly independent and in particular $\mathbf{v}_1,\ldots,\mathbf{v}_k$ is an orthonormal basis of its span, the projection of $\boldsymbol{\alpha}_i$ onto $\mathrm{span}(\mathbf{v}_1,\ldots,\mathbf{v}_k)$ is

$$\sum_{\ell=1}^k \mathbf{v}_\ell \left\langle\sum_{j=1}^r \sigma_j u_{j,i} \mathbf{v}_j,\mathbf{v}_\ell\right\rangle = \sum_{\ell=1}^k \sigma_\ell u_{\ell,i} \mathbf{v}_\ell$$

which is the $i$-th row of $A_k$. The $k$-means objective of $A_k$ with respect to the ground-truth centers $\mathbf{c}_j$, $j=1,\ldots,k$, is $\|A_k - C\|_F^2$.

One more observation: the rank of $C$ is at most $k$. Indeed, there are $k$ different rows in $C$ so its row rank is $k$ if these different rows are linearly independent and less than $k$ otherwise.

Exercise: Show that for any matrices $A, B \in \mathbb{R}^{n \times m}$, the rank of the sum is less or equal than the sum of the ranks, that is,

$$\mathrm{rk}(A+B) \leq \mathrm{rk}(A) + \mathrm{rk}(B).$$

[Hint: Show that $\mathrm{col}(A+B) \subseteq \mathrm{col}(A) \cup \mathrm{col}(B)$.] $\lhd$

Theorem (Why Project): Let $A \in \mathbb{A}^{n \times d}$ be a matrix and let $A_k$ be the truncation above. For any matrix $C \in \mathbb{R}^{n \times d}$ of rank $\leq k$,

$$\|A_k - C\|_F^2 \leq 8 k \|A - C\|_2^2.$$

The content of this inequality is the following. The quantity $\|A_k - C\|_F^2$ is the $k$-means objective of the projection $A_k$ with respect to the true centers, that is, the sum of the squared distances to the centers. By the Matrix Norms and Singular Values Lemma, the inequality above gives that

$$\|A_k - C\|_F^2 \leq 8 k \sigma_1(A - C)^2,$$

where $\sigma_j(A - C)$ is the $j$-th singular value of $A - C$. On the other hand, by the same lemma, the $k$-means objective of the un-projected data is

$$\|A - C\|_F^2 = \sum_{j=1}^{\mathrm{rk}(A-C)} \sigma_j(A - C)^2.$$

If the rank of $A-C$ is much larger than $k$ and the singular values of $A-C$ decay slowly, then the latter quantity may be much larger. In other words, projecting may bring the data points closer to their true centers, potentially making it easier to cluster them.

Proof (Why Project): By the exercise preceding the statement, the rank of the difference $A_k - C$ is at most the sum of the ranks

$$\mathrm{rk}(A_k - C) \leq \mathrm{rk}(A_k) + \mathrm{rk}(-C) \leq 2 k$$

where we used that the rank of $A_k$ is $k$ and the rank of $C$ is $\leq k$ since it has $k$ distinct rows. So by the Matrix Norms and Singular Values Lemma,

$$\|A_k - C\|_F^2 \leq 2k \|A_k - C\|_2^2.$$

By the triangle inequality for matrix norms,

$$\|A_k - C\|_2 \leq \|A_k - A\|_2 + \|A - C\|_2.$$

By the Low-Rank Approximation in the Induced Norm Theorem,

$$\|A - A_k\|_2 \leq \|A - C\|_2$$

since $C$ has rank at most $k$. Putting these three inequalities together,

$$\|A_k - C\|_F^2 \leq 2k (2 \|A - C\|_2)^2 = 8k \|A - C\|_2^2.$$

$\square$

In [11]:
d, n, offset = 1000, 100, 3
X1, X2 = two_clusters(d, n, offset)
X = vcat(X1, X2);


In reality, we cannot compute the matrix norms of $X-C$ and $X_k-C$ as the true centers are not known. But, because this is simulated data, we happen to know the truth and we can check the validity of our results in this case. The centers are:

In [12]:
C1 = reduce(hcat, [vcat(-offset, zeros(d-1)) for i=1:n])'
C2 = reduce(hcat, [vcat(offset, zeros(d-1)) for i=1:n])'
C = vcat(C1, C2);

In [13]:
Fc = svd(X-C)
plot(Fc.S, legend=false)

Out[13]:
In [14]:
frob = sum((Fc.S).^2)

Out[14]:
198953.15831164474
In [15]:
top_sval_sq = (Fc.S[1])^2

Out[15]:
2032.8205617726574
In [16]:
F = svd(X)
frob_proj = sum((F.S[1] * F.U[:,1] * F.Vt[1,:]' - C).^2)

Out[16]:
1477.471629781237