Singular value decomposition (SVD)

Theorem Every (real) m×n matrix A can be written as a product A = USV^T. Here, U and V are orthogonal matrices, with U being m×m and V being n×n. S is m×n, and has the form S =

s₁	0	0	...	0	...	0
0	s₂	0	...	0	...	0
...	...	...	...	...	...	...
0	0	0	...	s_r	...	0
0	0	0	...	0	...	0
...	...	...	...	...	...	...
0	0	0	...	0	...	0

The diagonal entries are positive, and are ordered from greatest to least; r is the rank of A. The s_k's are called the singular values of A.

Applications Gilbert Strang has called this theorem the "Fundamental Theorem of Linear Algebra." The SVD contains very explicit information concerning everything one would want to know about a matrix.

Condition number The condition number of a square, invertible matrix A is defined by
cond(A) = s₁/s_n.
It measures how many significant digits are preserved when one tries to solve Ax = b. For example, if b is known to 6 digits and cond(A) = 10³, then x is known to 6 -3 = 3 digits.
Numerical rank The rank of A is r, the number of (positive) singular values. The numerical rank of A is
rank(A,t) = # of singular values greater than a tolerance t.
Again, this is useful in working with problems having finite precision.
Least squares The solution to finding a minimum of $\| Ax - b \|$ is easily done with the SVD. First we rewrite the problem using the SVD for A.
$\| Ax - b\| = \| USV^Tx -UU^T b\|=\|Sz - c\|$, where $z = V^Tx$ and $c = U^T b$
Then we note that
$\|Sz - c \|^2 = \sum_{k=1}^r (s_kz_k-c_k)^2 +\sum_{k=r+1}^n c_k^2$
Choosing $z_k = c_k/s_k$ for $ k = 1 ... r$ and $z_k = 0$ for $k = r+1 ... n$ then yields $\|Sz - c \|^2= \sum_{k=r+1}^n c_k^2$, which not only solves the problem, but also gives the solution $x = Vz$ with smallest length $\| x \| = \|z\|$.

Finding the SVD of a matrix The matrix S above is just the matrix of A relative to new bases for both the input and output spaces.

The input space basis The matrix A^TA is real and symmetric, so there is an orthonormal basis of eigenvectors relative to which it is diagonal. Let the eigenvalues of A^TA be z₁, ..., z_n, and the corresponding basis of eigenvectors be {v₁, ..., v_n}. The eigenvalues are nonnegative, as this calculation shows:
A^TA v_k = z_k v_k
v_k^T A^T A v_k = z_k v_k^T v_k
|| A v_k ||² = z_k || v_k ||²
|| A v_k ||² = z_k (|| v_k || =1 )
This calculation also shows that a vector v is in the null space of A if and only if it is an eigenvecor corresponding to the eigenvalue 0. If r = rank(A), then "rank + nullity = # of columns" tells us that the the nullity(A) = n - r. This means that there are r eigenvectors for the remaining eigenvalues. List these as z₁ >= z₂ >= ... >= z_r > 0. Our input basis is now chosen as {v₁, ..., v_r, v_r+1, ..., v_n }. The numbering is the same as that for the eigenvalues. We now define the matrix V via
V = [ v₁ ... v_r v_r+1 ... v_n ].
The output space basis For k = 1, ..., r, let
u_k = A v_k / || A v_k ||
u_k = A v_k / z₁^1/2
We can also write this as the following equation:
A v_k = z_k^1/2 u_k .
The u_k's are orthonormal, for k = 1, ..., r. Again, we see this from these equations.
u_j^T u_k = v_j^T A^T A v_k / z_j^1/2 z_k^1/2
u_j^T u_k = z_k v_j^T v_k / z_j^1/2 z_k^1/2
The orthonormality of the v's implies that the right side above is 0 unless j = k, in which case it is 1. Thus, the u's are orthonormal. Fill this set out with m - r vectors to form an orthonormal basis for the output space, R^m. This gives us the basis {u₁, .., u_r, u_r+1, .., u_m}. As before, we define the m×m orthogonal matrix
U = [u₁, .., u_r, u_r+1, .., u_m].
The matrix of A relative to the new bases We let S be M_A. We compute it via the formula for the matrix of a linear transformation.
S = [ [Av₁]_U [Av₂]_U ... [Av_r]_U [Av_r+1]_U ... [Av_n]_U ]
The v's with k = r+1,...,n, are all in the null space of A. Thus the last n - r vectors are 0, and so are their corresponding columns relative to the basis of u's. The other vectors we get from the definition of the u's for k = 1, ..., r. The end result is this.
S = [ [z₁^1/2 u₁]_U [z₂^1/2 u₂]_U ... [z_k^1/2 u_r]_U [ 0 ]_U ... [ 0 ]_U ]
S = [ z₁^1/2 [u₁]_U ... z_r^1/2 [u_r]_U 0 ... 0 ]
S = [ z₁^1/2 e₁ ... z_r^1/2 e_r 0 ... 0 ].
If we let s_k = z_k^1/2 for k = 1, ..., r, we get the same S as the one given in the statement of the theorem. These s_k's are the singular values of A. The matrix S is related to A via multiplication by change-of-basis matrices. The matrix U changes from new output to old output bases, and V changes from new input to old input bases. Since V^T = V^-1, we have that V^T changes from old input to new input bases. In the end, this gives us A =USV^T

Example Consider the matrix A =

 2  -2
 1   1
-2   2

Here, A^TA =

 9  -7
-7   9

The eigenvalues of this matrix are z₁ = 16 and z₂ = 2. The singular values are s₁ = 4 and s₂ = 2^1/2. We can immediately write out what S is. We have S =

 4 0
 0 2^1/2
 0 0

The eigenvector corresponding to 16 is v₁ = 2^-1/2(1,-1)^T, and the one corresponding to 4 is v₂ = 2^-1/2(1,1)^T. Hence, we see that V =

 2^-1/2   2^-1/2
-2^-1/2   2^-1/2

Next, we find the u's.

u₁ = A v₁ / z₁^1/2
u₁ = 2^-1/2 (4, 0, -4)^T/4
u₁ =( 2^-1/2 , 0, - 2^-1/2 )^T.

A simlar calculation gives us
u₂ =( 0, 1, 0 )^T.

We now have to add to these to a "fill" vector
u₃ =( 2^-1/2 , 0, 2^-1/2 )^T
to complete the new output basis. This finally yields U =

 2^-1/2 0  2^-1/2
 0   1  0
-2^-1/2 0  2^-1/2