Singular value decomposition

Joshua R. Davis, Carleton College, Math 134, Fall 2024

Introduction

When we speak of “decomposing” a matrix A, we mean that we write A as a product of other matrices --- hopefully, matrices that are easier to understand. For example, diagonalizing A means finding an invertible matrix P and a diagonal matrix D such that A = P D P^-1. Another name for diagonalization is the eigenvalue decomposition, and it’s useful because D is easy to understand, has the same eigenvalues as A, and has eigenvectors that are closely related to A’s.

The singular value decomposition (SVD) is similar to, but different from, the eigenvalue decomposition. The SVD applies to any m x n matrix A., and it writes A as A = U S V^T, where U is m x m orthogonal, V is n x n orthogonal, and S is m x n and diagonal-ish. To compute the SVD, we orthogonally diagonalize A^T A as A^T A = V D V^T. Then we set S to the square root of D. Then we compute U from the columns of V and the diagonal entries of S.

This notebook shows some of the preliminary discussion leading up to the SVD. To start, execute the following code cell. (You aren’t expected to understand it. It lets us generate random points in the unit sphere in three dimensional space.)

Clear[rayUniform];
rayUniform[] :=
Module[{v},
v = Table[RandomVariate[NormalDistribution[]], {3}];
v / Norm[v]];
rayUniform[n_]:= Table[rayUniform[], {n}];

First example

Here is a 2 x 3 matrix A. It represents a linear transformation from 3D space to 2D space.

a = {
{4, -1, 2}, 
{0, 2, 1}
};
MatrixForm[a]

Here are 1000 points on the unit sphere in 3D space. And here they are transformed to 2D space by A. Intuitively, A flattens the sphere down to a plane while smearing it out along that plane.

data = rayUniform[1000];
ListPointPlot3D[data, BoxRatios->{1, 1, 1}]
ListPlot[Map[Function[v, a.v], data], AspectRatio->1, PlotRange->{{-5, 5}, {-5, 5}}]

Here is the eigensystem of A^T A.

{{lam1, lam2, lam3}, {v1, v2, v3}} = Eigensystem[Transpose[a].a]

The eigenvectors are guaranteed to be orthogonal, but let’s ensure that they’re orthonormal.

v1 = v1 / Norm[v1]
v2 = v2 / Norm[v2]
v3 = v3 / Norm[v3]

Here is the long axis of the 2D ellipse, along with its length.

a.v1
Sqrt[lam1]

Here is the short axis of the 2D ellipse and its length.

a.v2
Sqrt[lam2]

Here is the analogous third thing. It doesn’t really appear in the picture, because it’s the stuff that got crunched by A. (More precisely, v3 got crunched down to A v3 = 0.)

a.v3
Sqrt[lam3]

Second example

Here is another example. Figure out what each chunk of results is telling you, referring back to the previous section as needed.

a = {
{2, -1, 2}, 
{1, 2, -2}
};
MatrixForm[a]

data = rayUniform[1000];
ListPointPlot3D[data, BoxRatios->{1, 1, 1}]
ListPlot[Map[Function[v, a.v], data], AspectRatio->1, PlotRange->{{-5, 5}, {-5, 5}}]

{{lam1, lam2, lam3}, {v1, v2, v3}} = Eigensystem[Transpose[a].a]

v1 = v1 / Norm[v1]
v2 = v2 / Norm[v2]
v3 = v3 / Norm[v3]

a.v1
Sqrt[lam1]

a.v2
Sqrt[lam2]

a.v3
Sqrt[lam3]

SVD

Let A be either of the 2 x 3 matrices above. To get at the SVD of A, we first compute the diagonalization of A^T A.

MatrixForm[Transpose[a].a]
{{lam1, lam2, lam3}, {v1, v2, v3}} = Eigensystem[Transpose[a].a]

The eigenvectors are guaranteed to be orthogonal, but let’s make sure they’re orthonormal. And let’s put them into the columns of a matrix V.

v1 = v1 / Norm[v1]
v2 = v2 / Norm[v2]
v3 = v3 / Norm[v3]
v = Transpose[{v1, v2, v3}];
MatrixForm[v]

We compute the singular values and form the matrix S out of them.

sig1 = Sqrt[lam1]
sig2 = Sqrt[lam2]
sig3 = Sqrt[lam3]
s = {
{sig1, 0, 0}, 
{0, sig2, 0}
};
MatrixForm[s]

We compute each vector u from the corresponding singular value and v. We put the u vectors into the columns of a matrix U. (Along the way, the results get a little complicated, so we ask Mathematica to simplify.)

u1 = Simplify[a.v1 / sig1]
u2 = Simplify[a.v2 / sig2]
u = Transpose[{u1, u2}];
MatrixForm[u]

So in summary here are U, S, and V. And let’s check that U S V^T = A. Yippee.

MatrixForm[u]
MatrixForm[s]
MatrixForm[v]
MatrixForm[u.s.Transpose[v]]
MatrixForm[a]

The SVD is so important that of course Mathematica has it built in. Here’s how you’d compute it, if you weren’t trying to learn it at the same time. Compared to our manual answers above, the first column of U is negated, and the first column of V is negated, but that’s okay.

{u, s, v} = SingularValueDecomposition[a];
MatrixForm[u]
MatrixForm[s]
MatrixForm[v]

Numerical robustness

When computations with real numbers are done on a computer, the real numbers are almost always approximated as floating-point numbers. Put simply, this means that their decimal expansions are truncated. In other words, the number that you think that you’re using may not be exactly the number that the computer is really using. This issue is called numerical imprecision or round-off error.

A classic example is 1/3. Mathematica, unlike most computer programs, can represent this number exactly. But you can also ask Mathematica for a numerical approximation, which is close to, but not exactly equal to, the number that you wanted.

1/3
N[1/3]

(Mathematica keeps more than six decimal digits in its memory. It’s showing us a shortened form of what it is tracking. But the point still stands, that the numerical approximation has to be truncated somewhere.)

Another example is that numbers close to 0 might get rounded to 0. Or they might act like 0 in computations, even though they’re not 0. The following code asks Mathematica: How small does a number have to be, for you to start making bad results with it? Then it shows an example of a bad result.

$MachineEpsilon
1.0 + $MachineEpsilon

Approximating real numbers with floating-point numbers means that numerical calculations on a computer are almost always slightly wrong. Good algorithms are able to take slightly wrong inputs and produce slightly wrong outputs; we say that they are numerically robust. Bad algorithms, when given slightly wrong inputs, can produce wildly wrong outputs.

Algorithms that rely heavily on division tend not to be robust. Why? Because small changes to a number x can produce big changes to its reciprocal 1 / x. In fact, if a small number x gets rounded to 0, then suddenly 1 / x isn’t a number at all.

Unfortunately, many of the algorithms that we learn in an introductory linear algebra course are not robust. That’s because they’re designed for educational value and for being done by hand. Robust alternatives tend to be more complicated, so they’re not suitable for introductory students.

For example, the process of row-reducing a matrix relies heavily on division. It is not robust. Row reduction is almost never used on actual computers.

For another example, the process of diagonalizing a matrix is not robust, because slight changes to a matrix can produce wildly different eigenspace structures. If you set epsilon to 0 in the following code, then you get a classic non-diagonalizable matrix. If you set epsilon to a really small floating-point number, such as 0.000000001, then you get a nearly identical matrix that is diagonalizable. But notice that the eigenvectors are really close to each other. Intuitively, the two eigenvectors collapse down to one eigenvector, as epsilon goes to 0. This is a nightmare for robustness.

epsilon = 0.000001;
a = {
{1 + epsilon, 1}, 
{0, 1 - epsilon}
};
MatrixForm[a]
Eigensystem[a]

Why is the SVD so important?

The SVD is important because:
	A. Many linear algebra problems can be reframed in terms of the SVD.
	B. The SVD is numerically robust.
Why is the SVD robust? Because the only complicated step in it is diagonalizing a symmetric matrix A^T A. But isn’t diagonalization non-robust, because eigenspaces can collapse? In general yes, but not for symmetric matrices, because their eigenspaces are always orthogonal and of full dimension. So the collapsing cannot happen. But doesn’t the process of diagonalizing a symmetric matrix still require us to do a bunch of row reductions, to find the eigenvectors? Yes, but there are other approaches to finding eigenvectors, that avoid row reduction. Diagonalization of symmetric matrices can be made numerically robust, if we pick the details wisely.

We’re not going to study the details, because they take a while to learn. For example, professional-level computer code computes the SVD of A without comping A^T A at all. Why? Because numerical imprecision in A gets magnified when we compute A^T A. Experts in this subject strive to remove every bit of non-robustness that they can.

Example: Rank

We usually compute the rank of a matrix A by row-reducing and counting the pivot 1s. But row-reducing is not robust. A more robust alternative is to compute the SVD of A and count how many singular values are positive. In theory, that’s equal to the rank. Nice.

In practice, things get a little confusing when there are singular values close to 0. Maybe they’re non-zero just because of numerical imprecision? So we set a threshold, based on the precision of the computer, and we regard singular values smaller than that threshold to be 0. In this way, we obtain an approximate rank. It’s the best we can do, because fundamentally the entire concept of rank is not robust. I mean, two matrices, whose entries are close to each other, can have wildly different ranks. But at least the SVD doesn’t introduce additional imprecision into the problem.

Example: Condition number

Any square matrix A has a condition number, which measures “how non-invertible” A is. The condition number of A gives important insight into how robust computations with A are going to be.

Like rank, the condition number can be defined in a couple of different ways. A robust way to define it is as the ratio of the greatest singular value to the least singular value. The following code generates a random square matrix A and computes its condition number.

a = Table[RandomVariate[NormalDistribution[0, 8]], {3}, {3}];
MatrixForm[a]
{u, s, v} = SingularValueDecomposition[a];
MatrixForm[s]
s[[1, 1]] / s[[3, 3]]

Question: What is the condition number for the identity matrix? What is the condition number for a non-invertible matrix? (Hint: The SVD is easy to compute for diagonal square matrices.)

Example: Pseudo-inverse

The reduced SVD is used to construct the pseudoinverse of a matrix, which in practice is what actually powers the least-squares machinery and hence linear regression. I don’t have code for it right now, but we’ll discuss it in class.