The Singular Value Decomposition (SVD)¶

In this notebook we take a first look at one of the most important matrix decompositions, the SVD. The SVD gives what is arguably the "best" basis for the row and column spaces of a matrix, and reveals the "true nature" of a matrix in a unique way. It is heavily used in nearly all applications of linear algebra, especially to non-square matrices, especially for data analysis.

Which orthonormal basis?¶

Orthonormal bases are "nice", and make it especially easy to do things like projections or changes of basis. However, for any matrix $A$ with rank > 1, there are infinitely many possible orthonormal bases for $C(A)$ and $C(A^T)$. Are they all equally good?

Let's look at a simple $2\times2$ matrix $$A = \begin{pmatrix} 1 & 1 \\ -1 & \frac{1}{4} \end{pmatrix} .$$ This matrix is rank 2, so $C(A) = C(A^T) = \mathbb{R}^2$. The obvious orthonormal basis would be the unit vectors (1,0) and (0,1). But is there a better choice? A choice that comes from $A$ itself?

Let's see what $A$ does to two orthonormal vectors $v_1$ and $v_2$ as we rotate the "input" basis.

Some key observations:¶

• $A$ takes points $v_1, v_2$ on the unit circle and maps them to an ellipse. In general ($n$ dimensions), linear tranformations map the unit sphere to an ellipsoid.

• In general, even if $v_1 \perp v_2$, $Av_1$ and $Av_2$ are not orthogonal.

• However, for a particular choice of orthonormal input basis, we get orthonormal outputs.

The “right” bases for A¶

In fact, it turns out that we can always do this for any $m \times n$ matrix $A$ of rank $r$:

• We can choose an orthonormal basis $v_1, v_2, \ldots , v_r \in \mathbb{R}^n$ of $C(A^T) = N(A)^\perp$ such that $Av_1, Av_2, \ldots, Av_r$ are an orthogonal basis for $C(A)$!!!

• To make an orthonormal basis for $C(A)$, we simply divide $Av_k$ etcetera by their lengths $\sigma_k = \Vert Av_k\Vert$ to get orthonormal vectors $u_k = Av_k /\sigma_k$.

• Geometrically, $A$ transforms spheres to ellipsoids with axes oriented along $u_k$ (the "principal axes"), and the "semi-axes" (half-diameters) of the ellipsoid are the $\sigma_k$.

(This is not obvious. We won't have the tools to prove this, or to find these vectors, until later in 18.06.)

Equivalently: $$ Av_k = \sigma_k u_k $$ The matrix $A$ transforms orthormal vectors $v_k$ into orthonormal vectors $u_k$ "stretched" by the factors $\sigma_k$.

The $u_k$ and $v_k$ are called the left and right singular vectors and the $\sigma_k$ are called the singular values. Together, they form the singular value decomposition (the SVD).

In particular, if we put the above relation in matrix form, we have: $$ A \underbrace{ \begin{pmatrix} v_1 & \cdots & v_r \end{pmatrix} }_V = \underbrace{ \begin{pmatrix} u_1 & \cdots & u_r \end{pmatrix} }_U \underbrace{ \begin{pmatrix} \sigma_1 & & & \\ & \sigma_2 & & \\ & & \ddots & \\ & & & \sigma_r \end{pmatrix} }_\Sigma $$

In fact, since $C(V)=C(A^T)=N(A)^\perp$, so that $N(V^T) = N(A)$, we can write: $$ \boxed{A = U\Sigma V^T} = AVV^T = A[VV^T + \underbrace{(I-VV^T)}_{\mbox{proj. to }N(A)}] $$ which is in the form of a matrix factorization, the SVD.

In Julia, you can get these vectors and values with the svd function:

Modulo an arbitrary choice of sign ($v_1$ is flipped in sign), this is exactly the case of $\theta \approx 29^\circ$ that we found above which gave $Av_1 \perp Av_2$. Julia has now computed a much more precise angle where this happens:

Any matrix = sum of rank-1¶

The SVD $A = U\Sigma V^T$ can be thought of as: to compute $Ax$ for any $x$ you (1) compute $x$'s components in the $V$ basis with $V^Tx$, then (2) multiply each coefficient by the $\sigma_k$, and (3) add up in the $U$ basis. In the SVD basis, any matrix $A$ acts like a bunch of scalars $\sigma_k$ once you have decomposed $x$ into the right input and output basis!

• Note: are we missing anything when we compute $V^Tx$? $V$ is only a basis for $C(A^T)$, so what happens if $x$ does not lie entirely in the row space of $A$? Any remaining components $x$ must be in $C(A^T)^\perp = N(A)$, so $A$ would send those other components to zero! So, to find out what $A$ does to any vector $x$, we only need to know the components of x in the row space.

This can be also written $$ \boxed{A = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T + \cdots + \sigma_r u_r v_r^T} $$ That is, any matrix $A$ of rank $r$ can be written as a sum of r rank-1 matrices with orthogonal vectors, weighted by the singular values $\sigma$.

Data “Compression”¶

It is conventional to sort the singular values in descending order $$ \boxed{\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r} $$ In this way, the earlier singular vectors describe the more "important" part of the matrix $A$: the part that has a bigger effect on the output of $Ax$.

Because the SVD decomposes A into orthogonal components, sorted by "importance" it is extremely useful in data analysis for pulling out the "most important" parts of a data. There is a whole statistical technique, principal component analysis (PCA) that is based on this.

If we take only the first k (< r) singular values of $A$, we can think of it as an "approximate" version of A that throws out the "less important" parts. This is called a low-rank approximation of a matrix: $$ A \approx \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T + \cdots + \sigma_k u_k v_k^T $$

A great way to visualize this is to look at images, which can be thought of as matrices of numbers representing the red, green, and blue (RGB) components of the image.

Moral of the story: if k is large enough,

U[:,1:k] * Diagonal(σ[1:k]) * V[:,1:k]'

is a good approximation to A often requiring less storage and less computation, revealing the "most important" parts of A.

Nearly rank-deficient matrices¶

Consider the matrix $$ B = \begin{pmatrix} 1 & 2 \\ 1 & 2.01 \end{pmatrix} $$ The rows (and columns) of $B$ are nearly linearly dependent. Maybe they were "supposed" to be exactly dependent but there was some error (e.g. measurement error) in constructing $B$.

The SVD of $B$ reveals that fact, because one of the singular values is much smaller than the other:

We say that the matrix $B$ is ill-conditioned. The condition number of a full-rank matrix is the ratio $\sigma_1 / \sigma_r = \sigma_\max / \sigma_\min$ of largest to smallest singular value. It is computed in Julia by the cond function:

If we plot what $B$ does to a circle, it almost projects it to a line. The output ellipse is so thin we can't really see it is an ellipse.

So, $B$ is almost the rank-1 matrix $\tilde{B} = \sigma_1 u_1 v_1^T$:

What happens if we blindly try to solve $Bx = b$ without knowing this fact?

First, because $B$ is nearly singular, $x = B^{-1} b$ tends to be pretty large for most $b$: we almost divided by zero:

Second, the result can be hugely sensitive to small changes in b:

This is a general problem with "ill-conditioned" (nearly singular) systems of equations: the results can be very sensitive to tiny errors in the inputs (or to tiny roundoff errors during the computation).