Dusty Wilson 
Singular Value Decomposition
Western Washington University
Masters Project
Winter 2001

Contents

(1.)  The Dream That Motivates Singular Value Decomposition
(2.)  Exploring the World of Real Valued Symmetric Matrices
(3.)  In a World Where Dreams Come True: the Singular Value Decomposition Theorem
(4.)  Matrix Reduction Through Singular Value Decomposition
(5.)  Numerical Rank Approximation
(6.)  Image Compression
Appendix (1.).  Basic definitions from linear algebra.
Appendix (2.)  Mathematica notebook which generates examples of the spectral theorem 
Appendix (3.)  Mathematica notebook which calculates the singular value decomposition of a matrix as well as generate an animation of the geometric interpretation for two by two case
Appendix (4.)  Matlab m-file which performs singular value decomposition for image compression
Appendix (5.)  Project Proposal

(1.)  The Dream That Motivates Singular Value Decomposition

Objectives:
(i.)  Describe the situation
(ii)  Consider simple linear transformations
(iii) Propose a decomposition of our situation into the product of simple linear transformations

Let (Subscript[M, m x n]) be a real valued matrix and let x=((((Subscript[x, 1]),...,(Subscript[x, n])))^T) be a column vector in (ℝ^n).  We think of M: (ℝ^n)↦(ℝ^m) as a linear transformation where x↦M x.

Example 1(a)
Let (Subscript[M, 3 x 2])=(({{1, 2}, {3, -1}, {-2, -3}})) and (Subscript[x, 2 x 1])=(({{1}, {2}})).  Then x↦M x implies that (({{1}, {2}}))↦(({{1, 2}, {3, -1}, {-2, -3}}))·(({{1}, {2}}))=(({{5}, {1}, {-8}}))

Our Dream.
We dream of finding the most conceptually simple representation of M and then exploring some of the applications of this conceptually simple representation.  The goal of finding the simplest representation of M leads us to ask what a conceptually transparent transformation would look like.  The most transparent representation of M would be one where M is factored (broken up) into a series of simple transformations.  So, what are the simple transformations into which we hope to decompose M?

Diagonal Matrices.
The simplest matrices are diagonal matrices (Subscript[Σ, r]):(ℝ^r)↦(ℝ^r) whose diagonal σ=((((Subscript[σ, 1]),...,(Subscript[σ, r])))^T) lives in (ℝ^r).  Alternately, we write (Subscript[Σ, r])=diag⟨σ⟩ and σ=diag((Subscript[Σ, r])).  When considered as a transformation, x↦(Subscript[Σ, r])x , (Subscript[Σ, r])x=((((Subscript[σ, 1])(Subscript[x, 1]),...,(Subscript[σ, r])(Subscript[x, r])))^T) where x∈(ℝ^r).  That is, (Subscript[Σ, r]) "stretches" the (Subscript[x, i]) component of x by a factor of (Subscript[σ, i]).  The concept of a diagonal matrix can be extended to include matrices which are not square.  

We will call (Subscript[Λ, m x n]) "diagonal" if Λ=(Subscript[(({{(Subscript[Σ, r]), (Subscript[0, r x (n-r)])}, {(Subscript[0, (m-r) x r]), (Subscript[0, (m-r) x (n-r)])}})), m x n])

Here x↦Λ x=((((Subscript[σ, 1])(Subscript[x, 1]),...,(Subscript[σ, r])(Subscript[x, r]),0·(Subscript[x, r+1]),...,0·(Subscript[x, n])))^T)=((((Subscript[σ, 1])(Subscript[x, 1]),...,(Subscript[σ, r])(Subscript[x, r]),0,...,0))^T).  In terms of transformations, Λ "stretches" the (Subscript[x, i]) component of x by a factor of (Subscript[σ, i]) for i=1,...,r and sends the remaining n-r components of x to zero.

Example 1(b)
Given x=(({{1}, {2}})) as above, let (Subscript[Σ, 2])=diag⟨x⟩=(({{1, 0}, {0, 2}})).  Suppose that we wanted (Subscript[Λ, 3 x 2]).  We could easily create this Λ=(({{Σ}, {(Subscript[0, 1 x 2])}}))=(({{1, 0}, {0, 2}, {0, 0}})).

Orthonormal Matrices.
Clearly all matrices do not have a "diagonal" representation, so we will consider a second type of conceptually simple transformation (matrix).  A transformation V of this type preserves length and preserves orientation.  So, if x↦V x, then length(x)=length(V x) and orientation(x,y)=orientation(V x,V y) for all x,y∈(ℝ^n).  (We describe this mathematically by saying that ||x||=||V x|| and ⟨x,y⟩=⟨V x,V y⟩).  It is intuitively obvious from the restrictions on length and orientation preservation that if x∈(ℝ^n), then V:(ℝ^n)↦(ℝ^n).  That is, V∈(ℝ^(n x n)).  A matrix that satisfies these requirements is said to be orthonormal.  In terms of transformations, V  rotates the vector x around zero.

Combining the two classes of linear transformations, we dream that M can be represented as the product of orthonormal and "diagonal" matrices.  Since orthogonal matrices are square, it is reasonable to look for U, V, and Λ such that M=U·Λ·(V^T) where (Subscript[U, m x m]) and (Subscript[V, n x n]) are orthogonal matrices and (Subscript[Λ, m x n]) is a "diagonal" matrix.  Thus, we want to represent the transformation M as a rotation, then a component-wise stretch, followed by a second rotation.  In order to do this, we are in search of a function - SVD - into which we input matrices and which outputs U, V, and Λ.  That is, SVD((Subscript[M, m x n]))={(Subscript[U, m x m]),(Subscript[Λ, m x n]),(Subscript[V, n x n])}.

In our search for the function SVD, we will look first for the right most member of our product - (Subscript[V, n x n]) - an n x n matrix.  Since our function SVD accepts the input (Subscript[M, m x n]), the most obvious n x n matrix SVD has to work with is A=(Subscript[((M^T)M), n x n]).  In general, A itself will not satisfy the requirements of orthonormalness.  However, we dream that we can construct an orthonormal matrix V from A.  

In order to do this, we must take stock of the tools with which we may work: (a.)  A is a square real valued matrix, and (b.)  A is symmetric.  Given the tools of realness and symmetry, we will try to find V.

(2.)  Exploring the World of Real Valued Symmetric Matrices

Objectives:
(i)   Prove non-zero symmetric matrices have an eigenvalue and non-zero eigenvector which live in ℂ and (ℂ^n) respectively
(ii)  Prove that the eigenvalues and eigenvectors of a symmetric matrix live in ℝ and (ℝ^n) respectively
(iii) Prove that symmetric linear transformations are invariant
(iv) Prove that every real valued symmetric matrix has an orthogonal basis of eigenvectors
(v)  Discuss that change of bases of matrices and the implications of the real valued spectral theorem

Let A be the symmetric real valued matrix with which we ended part one.

Example 1(c)
Recall that A=(M^T)M=(((({{1, 2}, {3, -1}, {-2, -3}})))^T)·(({{1, 2}, {3, -1}, {-2, -3}}))=(({{14, 5}, {5, 14}})).  Clearly A is real valued and symmetric.

Let 𝒱 be a vector space and 𝕂 be a scalar field.

Let x,y∈𝒱 such that x=((((Subscript[x, 1]),...,(Subscript[x, n])))^T) and y=((((Subscript[y, 1]),...,(Subscript[y, n])))^T) where each (Subscript[x, i]) ,(Subscript[y, i])∈𝕂.

From this point forward, we specify that the inner product we are using is ⟨ x,y⟩=(∑_(i=1)^m)(Subscript[x, i])(Subscript[y, i]).  Note: all of the bases with which we are working are orthonormal bases, therefore, this is a positive definite inner product.

Define the norm henceforth as: ||x||=(Sqrt[⟨x,y ⟩]).

We begin by showing that A (assuming A is non-zero) has at least one eigenvalue λ and one non-zero eigenvector.

Theorem:
All non-zero n x n matrices A have a non-zero eigenvector with complex components and an eigenvalue λ∈ℂ.

□ Proof.
Let A be a n x n matrix with entries in ℂ.  Then the eigenvalues of A are found by determining what values of λ are solutions to the  polynomial equation:

Det(A-λ I)=0

Since Det(A-λ I) is a non-constant polynomial with complex coefficients, it has at least one solution λ∈ℂ.  (For a proof of this fundamental theorem from algebra, see Lang 279).  This implies that there exists an eigenvector z∈(ℂ^n) corresponding to the eigenvalue λ that is a non-trivial solution to the system of homogeneous equations A-λ I=(Subscript[0, (ℂ^n)]) and thus z≠(Subscript[0, (ℂ^n)]). ■

The existence of eigenvectors and eigenvalues in ℂ is not sufficient for us, we must show the existence of a non-zero eigenvector and an eigenvalue over ℝ.  (The following proof relies on a basic understanding of complex variables).

Theorem: Eigenpair existence theorem
Let A be a real symmetric n x n matrix with the eigenvalue λ in ℂ.  Then λ is real.  If z≠(Subscript[0, (ℂ^n)]) is a complex eigenvector corresponding to the eigenvalue λ, and z=x+i y where x,y∈(ℝ^n), then both x and y are real eigenvectors of A with eigenvalue λ, and x and y are not both equal to zero.

□ Proof.
Let z=((((Subscript[z, 1]),...,(Subscript[z, n])))^T) where z has the complex coordinates (Subscript[z, i]).  Then

z·(OverBar[z])=(OverBar[z])·z=(((OverBar[z]) )^T)z=(Subscript[(OverBar[z]), 1])(Subscript[z, 1])+…+(Subscript[(OverBar[z]), n])(Subscript[z, n])=|(Subscript[z, 1])((|)^2)+…+|(Subscript[z, n])((|)^n)>0  since z is nonzero.  (here " · " refers to the ordinary scalar product)

By assumption we have A z = λ z.  This implies

((OverBar[z])^T) A z=((OverBar[z])^T) λ z=λ ((OverBar[z])^T) z

Since this quantity is a 1 x 1 matrix, i.e., a scalar, it is equal to its transpose .  Reversing our equality we have:

((OverBar[z])^T)A z=(z^T)(A^T)(OverBar[z])=(z^T)A (OverBar[z]) (since A is symmetric)

But since A  is real valued we have thatA(OverBar[z])=(OverBar[A]) (OverBar[z])=(OverBar[A z])=(OverBar[λ z])=(OverBar[λ]) (OverBar[z])  and

λ ((OverBar[z])^T)z=(z^T)A (OverBar[z])=(z^T)(OverBar[λ]) (OverBar[z])=(OverBar[λ]) (z^T)(OverBar[z])=(OverBar[λ]) ((OverBar[z])^T)z 

Since we already showed that ((OverBar[z])^T)z >0  , it follows that λ=(OverBar[λ]), so λ is real.

Now, from A z=λ z and z=x+i y we get that
 
A z=A(x+i y)=A x+i A y = λ(x+i y)=λ x+i λ y

Since A, x, and y are all real, it follows that A x=λ x and A y = λ y.  Thus x and y are real eigenvectors of A with the real eigenvalue λ. ■ (Lang 214).

Alternately, we can show that the eigenvector z is real by noticing that Det(A-λ I)=0 has a solution for a real value of λ and therefore the system of real equations A-λ I=0 has a non-trivial solution z.  Clearly, z must be real.

Note: I titled the previous theorem the "eigenpair existence theorem."  It is very humorous to note that Mathematica's spell checker suggested replacing the word eigenpair with the word teenybopper.

Example 1(d)
The preceding two theorems guarantee that A has at least one non-zero eigenvalue and one non-zero eigenvector.  In fact, the matrix A=(({{14, 5}, {5, 14}})) has the eigenvalues: (Subscript[λ, 1])=19 and (Subscript[λ, 2])=9 while the corresponding eigenvectors are (Subscript[v, 1]^( * ))=(({{1}, {1}})) and (Subscript[v, 2]^( * ))=(({{-1}, {1}})).  (We use (Subscript[v, 1]^( * )) to notate a temporary value for (Subscript[v, 1]).)

Since we want to consider what happens to subspaces of 𝒱 under the symmetric real valued A, we will consider the following:

Observation: Let v be an eigenvector of A.  The subspace 𝒲=span⟨v⟩ is invariant under A.  (By definition of an eigenvector).

Theorem:
Let A:𝒱↦𝒱 be a symmetric linear transformation where 𝒱 is a vector space over the scalar field 𝕂 with a positive definite inner product.  Then if 𝒲 is a subspace of 𝒱 which is invariant under A, then (𝒲^⊥) (the perpendicular subspace of 𝒲) is also invariant under A.

□ Proof.
Suppose 𝒲 is invariant under A.  Let u∈(𝒲^⊥).  Then for all w∈𝒲 we have

⟨A u,w⟩=⟨u,A w⟩=0

By assumption we had that A w∈𝒲.  Hence A u∈(𝒲^⊥), thus (𝒲^⊥) is also invariant under A. ■  (Lang 219)

Theorem: Real Valued Spectral Theorem
Let 𝒱 be an n dimensional vector space over ℝ with a positive definite inner product.  Let A:𝒱 ↦ 𝒱 be a linear transformation.  Then A is symmetric with respect to the inner product if and only if 𝒱 has an orthonormal basis B consisting of eigenvectors of the linear transformation A.

□ Proof.
◇ (⇒): By the eigenpair existence theorem, there exists a non-zero eigenvector (Subscript[v, 1]) for A which we have shown to be real.  Let (Subscript[𝒲, 1])=span⟨(Subscript[v, 1])⟩ be a one dimensional subspace of 𝒱.  Then 𝒲 is invariant under A.  Furthermore, when restricted to the subspace (𝒲^⊥) of dimension  n-1, the transformation A is invariant and symmetric.  We may now view A as a symmetric transformation A:(𝒲^⊥) ↦ (𝒲^⊥).  Repeating the process n-times, we find the set of eigenvectors (Subscript[B, V])={(Subscript[v, 1]),...,(Subscript[v, n])} that form a basis for 𝒱 and have the characteristic that ⟨(Subscript[v, i]),(Subscript[v, j])⟩>0 [i=j]. (Note: our process halts when (Subscript[𝒲, n]^⊥)={(Subscript[0, (ℝ^n)])}).  Scaling each (Subscript[v, i]) by its length, we have an orthonormal basis for 𝒱 as desired. ◆
◇ (⇐):  Assume that 𝒱 has an orthonormal basis (Subscript[B, V]) consisting of the eigenvectors of the linear transformation A.  Then, with respect to the basis (Subscript[B, V]), the transformation A is a diagonal transformation and so A=(A^T).  ◆
As was to be shown.  ■

Example 1(e)
Given the matrix A=(({{14, 5}, {5, 14}})), by the real valued spectral theorem, we can find an orthonormal basis for A made of eigenvectors for A.  

(({{1}, {1}}))↦(Subscript[v, 1])=(({{(1/(Sqrt[2]))}, {(1/(Sqrt[2]))}})) and (({{-1}, {1}}))↦(Subscript[v, 2])=(({{-(1/(Sqrt[2]))}, {(1/(Sqrt[2]))}})).  

Define V=(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})).  We can show that these vectors are orthonormal by verifying that (V^T)·V=(Subscript[I, 2]).  (Here (Subscript[I, 2]) denotes the 2 x 2 identity matrix.)  

(V^T)·V=(((({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})))^T)·(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}}))=(({{1, 0}, {0, 1}}))

Notation:
The transformation x↦M x where M: (ℝ^n)↦(ℝ^m) can be thought of in the following way.  Let W={(Subscript[w, 1]),...,(Subscript[w, n])} and (W^( '))={(Subscript[w, 1]^( ')),...,(Subscript[w, n]^( '))} be bases for (ℝ^n).  Then for x∈(ℝ^n), there exist (Subscript[α, 1]),...,(Subscript[α, n]) and (Subscript[β, 1]),...,(Subscript[β, n]) such that x=(Subscript[α, 1])(Subscript[w, 1])+...+(Subscript[α, n])(Subscript[w, n]) and x=(Subscript[β, 1])(Subscript[w, 1]^( '))+...+(Subscript[β, n])(Subscript[w, n]^( ')) Now we have the two equalities:

M x=M((Subscript[α, 1])(Subscript[w, 1])+...+(Subscript[α, n])(Subscript[w, n]))=(Subscript[α, 1])M (Subscript[w, 1])+...+(Subscript[α, n])M (Subscript[w, n]) 

and similarly

M x=M((Subscript[β, 1])(Subscript[w, 1]^( '))+...+(Subscript[β, n])(Subscript[w, n]^( ')))=(Subscript[β, 1])M (Subscript[w, 1]^( '))+...+(Subscript[β, n])M (Subscript[w, n]^( ')). 

Now if (Subscript[w, i])≠(Subscript[w, i]^( ')), then M (Subscript[w, i])≠M (Subscript[w, i]^( ')).  However, the sum is the same regardless of the basis that we use.  So, we can see that our transformation M will take on different forms depending upon the basis that we work with.  These representations are referred to as the coordinatization of x with respect to W and the coordinatization of x with respect to (W^( ')).  

While we have varied the basis of the domain, it is clear that we can similarly vary the basis of the codomain.

Define (Subscript[M, (Subscript[B, 2])]^(Subscript[B, 1]))(M) to be the transformation M where x∈(ℝ^n) is coordinatized with respect to the basis (Subscript[B, 1]) and M x∈(ℝ^m) is coordinatized with respect to the basis (Subscript[B, 2]).

Note:
We most often think about vectors in (ℝ^n) as being coordinitized with respect to the standard basis (Subscript[E, n])={(Subscript[e, 1]),...(Subscript[e, n])}.  (We use the subscript n to denote the dimension of our basis.  However, this additional subscript will be left off in cases where it does not raise ambiguity).  In this case, we would write M=(Subscript[M, (Subscript[E, m])]^(Subscript[E, n]))(M).

Interpretation: the "Real Valued Spectral Theorem" in terms of singular value decomposition.
Given the assumptions of the previous theorem, we begin with the transformation matrix A whose columns are linear combinations of the standard basis (Subscript[E, n]).  Now, by the real valued spectral theorem, there exists an orthonormal basis (Subscript[B, V]) such that the A takes on the form (Subscript[M, (Subscript[B, V])]^(Subscript[B, V]))(A)=diag⟨(Subscript[λ, 1]),...,(Subscript[λ, n])⟩.  Taking advantage of what we know of the standard basis and the known matrix diag⟨(Subscript[λ, 1]),...,(Subscript[λ, n])⟩, we have

A=(Subscript[M, (Subscript[E, n])]^(Subscript[E, n]))(A)=(Subscript[M, (Subscript[E, n])]^(Subscript[B, V]))((Subscript[I, n]))·diag⟨(Subscript[λ, 1]),...,(Subscript[λ, n])⟩·(Subscript[M, (Subscript[B, V])]^(Subscript[E, n]))((Subscript[I, n]))  (Here (Subscript[I, n]) refers to the identity transformation in (ℝ^n)).

However, the orthonormal nature of the bases (Subscript[E, n]) and (Subscript[B, V]) implies that ((((Subscript[M, (Subscript[E, n])]^(Subscript[B, V]))((Subscript[I, n]))))^T)=(Subscript[M, (Subscript[B, V])]^(Subscript[E, n]))((Subscript[I, n])) since the nature of the matrices clearly shows that they are inverses.  Rewriting (Subscript[M, (Subscript[E, n])]^(Subscript[B, V]))((Subscript[I, n])) as V, we are left with

A=(Subscript[M, (Subscript[E, n])]^(Subscript[E, n]))(A)=V·diag⟨(Subscript[λ, 1]),...,(Subscript[λ, m])⟩· (V^T)

This is known as the spectral decomposition of the symmetric real valued transformation A.

Example 1(f)
Let A=(({{14, 5}, {5, 14}})).  A is a real valued symmetric matrix.  Construct a diagonal matrix Λ whose (i^th) entry is the (i^th) largest eigenvalue of M and a matrix V whose column (i^th) is the normalized eigenvector of A corresponding to (Subscript[Λ, i, i])
V=(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})) and Λ=diag⟨19,9⟩=(({{19, 0}, {0, 9}})).  
Now we can verify that (({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})) ·(({{19, 0}, {0, 9}}))·(((({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})))^T)=(({{14, 5}, {5, 14}})) as expected.

So, we conclude by saying that we can find an orthonormal basis of eigenvectors for any real valued symmetric matrix.

(3.)  In a World Where Dreams Come True: the Singular Value Decomposition Theorem

Objectives:
(i)    Determine the singular values of our singular value decomposition
(ii)   Find the second orthonormal matrix U
(iii)  Construct a "diagonal" matrix from the singular values
(iv)  Prove the Singular Value Decomposition Theorem

Digging into the recesses of our subconsciousness, we remember that just before our digression into symmetric matrices, we had the statement, "Since our function SVD takes in the input (Subscript[M, m x n]), the most obvious n x n matrix SVD has to work with is A=(Subscript[((M^T)M), n x n])."

It has been shown that we can find an orthonormal basis of eigenvectors for any real valued symmetric transformation A.  It remains to be shown that any real valued transformation has a singular value decomposition.  That is, we need to show that a non-symmetric transformation has a singular value decomposition.

Let (Subscript[B, V])={(Subscript[v, 1]),…,(Subscript[v, n])} be the orthonormal basis of eigenvectors for A with the corresponding eigenvalues {(Subscript[λ, 1]),(Subscript[…λ, n])}.

Define V=[(Subscript[v, 1]),…,(Subscript[v, n])].  Then A=V·diag⟨(Subscript[λ, 1]),…,(Subscript[λ, n])⟩·(V^T).

Example 1(g)
Recall that M, A, V, and Λ are:
M=(({{1, 2}, {3, -1}, {-2, -3}})) , A=(({{14, 5}, {5, 14}})) , V=(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})) , and Λ=(({{19, 0}, {0, 9}}))

Theorem:
Let (Subscript[B, m x n]) and (Subscript[C, n x p]) be matrices where rank(B)=(Subscript[r, 1]) and rank(C)=(Subscript[r, 2]).  Then rank(B C)≤min{(Subscript[r, 1]),(Subscript[r, 2])}.  (Def. rank(A)).

□ Proof.
Let B and C be as above.  Now the domain of C is of dimension (Subscript[r, 2]).  Since the product B C can be seen as the composition B∘C, the rank(B C)≤(Subscript[r, 2]), the dimension of the domain of C.  The domain of (B^T) is of dimension (Subscript[r, 1]).  The product (((B C))^T)=(C^T)(B^T) can be seen as the composition (C^T)∘(B^T)which has rank((((C B))^T))≤(Subscript[r, 1]).  Since rank(B C)=rank((((B C))^T)), we have that the rank of a product is less than or equal to the rank of each factor in the product.  ■

Since rank(M)=rank((M^T))=r≤min{m,n}, we have that there are r nonzero eigenvalues of (M^T)M.  Thus, M (Subscript[v, 1]),…,M (Subscript[v, r])≠(Subscript[0, (ℝ^m)]) and M (Subscript[v, r+1]),…,M (Subscript[v, n])=(Subscript[0, (ℝ^m)]).  Alternately, we can write that Col(M)=Span{M (Subscript[v, 1]),…,M (Subscript[v, r])} .  (Def. Col(M)).

Note: as a direct consequence, we have that Null(M)=Span{(Subscript[v, r+1]),(Subscript[…v, n])}.  (Def. Null(M)).

To summarize, what we have so far is (a.) an orthonormal basis (Subscript[B, V]) for (ℝ^n) and (b.) a basis M (Subscript[v, 1]),...,M (Subscript[v, r]) for the column space of M.

This is almost the situation we want - all we need is (Subscript[B, U])={(Subscript[u, 1]),...,(Subscript[u, m])}, an orthonormal basis for (ℝ^m), such that M (Subscript[v, i])=(Sqrt[(Subscript[λ, i])])(Subscript[u, i]) , for i=1,2,3,...,n.  Since the column space of M lives in (ℝ^m), we will attempt to build our basis (Subscript[B, U]) from the r vectors M (Subscript[v, 1]),...,M (Subscript[v, r])∈(ℝ^m) which we already have.  We will attempt to construct the basis (Subscript[B, U]) from the (Subscript[u, i])=(1/(Sqrt[(Subscript[λ, i])]))M (Subscript[v, i]) for i=1,...,r as the only vectors we know of in (ℝ^m).  (Note: the scalar (1/(Sqrt[(Subscript[λ, i])])) is added in order to normalize the vector (Subscript[u, i]).  However, in order to anticipate this result, we would have to repeat the proof below with less generality.)

Construction: the Basis (Subscript[B, U])
◇ Let M, (Subscript[B, V]), and (Subscript[{(Subscript[λ, i])}, i=1]^n) be as above.  Now, construct the basis (Subscript[B, U])={(Subscript[u, 1]),...,(Subscript[u, m])} as follows.

Set (Subscript[u, i])=(1/(Sqrt[(Subscript[λ, i])]))M (Subscript[v, i])  for i=1,2,3,...,r where r is the number of non-zero eigenvalues of (M^T)M

⟨(Subscript[u, i]),(Subscript[u, j])⟩=⟨(1/(Sqrt[(Subscript[λ, i])]))M (Subscript[v, i]),(1/(Sqrt[(Subscript[λ, i])]))M (Subscript[v, j])⟩
=⟨(1/(Sqrt[(Subscript[λ, i])(Subscript[λ, j])]))(Subscript[v, i]),(M^T)M (Subscript[v, j])⟩
=(1/(Sqrt[(Subscript[λ, i])(Subscript[λ, j])]))⟨(Subscript[v, i]),(Subscript[λ, j])(Subscript[v, j])⟩
=((Subscript[λ, j])/(Sqrt[(Subscript[λ, i])(Subscript[λ, j])]))⟨(Subscript[v, i]),(Subscript[v, j])⟩
=1 [i=j] [i,j =1,2,3,...,r]  (because the(Subscript[v, i])'sare orthonormal)

Thus (Subscript[B, U])={(Subscript[u, 1]),...,(Subscript[u, r])} are orthonormal.  To complete our basis, we use the Gram-Schmidt process to find vectors {(Subscript[u, r+1]),...,(Subscript[u, m])} such that the vectors (Subscript[B, U]) form an orthonormal basis for (ℝ^m). ◆

Define the matrix U=[(Subscript[u, 1]),...,(Subscript[u, m])].

Corollary
It follows directly that ||M (Subscript[v, i])||=(Sqrt[(Subscript[λ, i])])≥0 and thus r is the number of positive eigenvalues.

For future purposes, we will define (Subscript[σ, i])=(Sqrt[(Subscript[λ, i])]) [i=1,...,n].  By renumbering, we may assume  (Subscript[σ, 1])≥(Subscript[σ, 2])≥...≥(Subscript[σ, n])≥0

Example 1(h)
Given the matrix M=(({{1, 2}, {3, -1}, {-2, -3}})) , we found the orthonormal matrix V=(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})) and the singular values {(Subscript[σ, 1]),(Subscript[σ, 2])}={(Sqrt[19]),(Sqrt[9])}.  Now, we can find (Subscript[B, U])={(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])} in the following manner.  For i=1,2, we will set (Subscript[u, i])=(1/(Subscript[σ, i]))M (Subscript[v, i]) where (Subscript[v, i]) is the (i^th) column of V.
(Subscript[u, 1])=(1/(Sqrt[19]))(({{1, 2}, {3, -1}, {-2, -3}}))·(({{(1/(Sqrt[2]))}, {(1/(Sqrt[2]))}}))=(({{(3/(Sqrt[38]))}, {(Sqrt[(2/19)])}, {-(5/(Sqrt[38]))}})) and (Subscript[u, 2])=(1/(Sqrt[9]))(({{1, 2}, {3, -1}, {-2, -3}}))·(({{-(1/(Sqrt[2]))}, {(1/(Sqrt[2]))}}))=(({{(1/(3 (Sqrt[2])))}, {-((2 (Sqrt[2]))/3)}, {-(1/(3 (Sqrt[2])))}})).  
Then we find (Subscript[u, 3]) by using the Gram-Schmidt orthogonalization.  So, (Subscript[u, 3])=(({{(11/(3 (Sqrt[19])))}, {(1/(3 (Sqrt[19])))}, {(7/(3 (Sqrt[19])))}})).

The following corollary is an interesting side effect of the preceding theorems

Corollary:
If (Subscript[v, i]) is an eigenvector of (M^T)M with the eigenvalue (Subscript[λ, i])≠0, then M (Subscript[v, i]) is an eigenvector of (Subscript[(M (M^T)), m x m]) with the same eigenvalue (Subscript[λ, i]).

□ Proof.
Let (Subscript[v, i]) be an eigenvector of (M^T)M with the eigenvalue (Subscript[λ, i])≠0.  Now consider M (M^T)(M (Subscript[v, i]))

M (M^T)(M (Subscript[v, i]))=M((M^T)M)(Subscript[v, i])=M((Subscript[λ, i])(Subscript[v, i]))=(Subscript[λ, i])M (Subscript[v, i])

Therefore M (Subscript[v, i]) is an eigenvector of M (M^T) with the corresponding eigenvalue (Subscript[λ, i]). ■

Having two orthonormal matrices U and V, we are desperately in need of a 'diagonal" matrix, so we will construct one as follows

Construction: the matrix Λ
◇ Let the "diagonal" matrix Λ be constructed as follows.
Begin by constructing a sub matrix  (Subscript[Σ, r ])=(Subscript[(({{(Subscript[σ, 1]), 0, ., ., 0}, {0, ., ., ., .}, {., ., ., ., .}, {., ., ., ., 0}, {0, ., ., 0, (Subscript[σ, r])}})), r x r])

Now define Λ to be Λ=(Subscript[(({{(Subscript[Σ, r]), (Subscript[0, r x (n-r)])}, {(Subscript[0, (m-r) x r]), (Subscript[0, (m-r) x (n-r)])}})), m x n]). ◆

Example 1(h)
Given the singular values (Subscript[σ, 1])=(Sqrt[19]) and (Subscript[σ, 2])=(Sqrt[9]) we have that (Subscript[Σ, 2])=(({{(Sqrt[19]), 0}, {0, (Sqrt[9])}})) and Λ=(({{(Subscript[Σ, 2])}, {1 x 2}}))=(({{(Sqrt[19]), 0}, {0, (Sqrt[9])}, {0, 0}}))

We are now ready to prove the singular value decomposition theorem.

Theorem: The Singular Value Decomposition (SVD) Theorem
Any m x n real valued matrix M can be decomposed into the product

M=U·Λ·(V^T)

where (Subscript[U, m x m]) and (Subscript[V, n x n]) are orthonormal matrices and (Subscript[Λ, m x n]) is a "diagonal" matrix whose entries are the square roots of the eigenvalues of the matrix (M^T)M.

□ Proof.
Consider the matrix (M^T)M.  We have already shown that

(M^T)M=V·diag⟨(Subscript[λ, 1]),...,(Subscript[λ, n])⟩·(V^T).  
 
Since the product of diagonal matrices is also a diagonal matrix, this is equal to 

V·diag⟨(Subscript[σ, 1]),...,(Subscript[σ, n])⟩·diag⟨(Subscript[σ, 1]),...,(Subscript[σ, n])⟩·(V^T)=V·(Λ^T)·Λ·(V^T).  (where we take r=n)

Define (Subscript[U, m x m]) to be an orthonormal matrix.  Then (U^T)U=I and we continue the equality with

V·(Λ^T)·I·Λ·(V^T)=V·(Λ^T)·(U^T)·U·Λ·(V^T)=(((U·Λ·(V^T)))^T)·U·Λ·(V^T)

By comparison, we have that M=U·Λ·(V^T) and our claim is proved.  ■

The preceding proof of the singular value theorem is beautiful in its simplicity, but does not go very far towards telling us how to find the matrix U.  So we will continue on towards a second proof of the singular value decomposition theorem - this one more computationally valuable.

We begin by considering the real valued  transformation matrix (Subscript[M, m x n])=(Subscript[M, (Subscript[B, D])]^(Subscript[B, C]))(M) with real entries when (ℝ^n) is coordinitized with respect to the orthonormal basis (Subscript[B, C]) and (ℝ^m) with respect to the orthonormal basis (Subscript[B, D]).

Our goal is to find bases (Subscript[B, V]) and (Subscript[B, U]) such that (Subscript[M, (Subscript[B, D])]^(Subscript[B, C]))(M) can be replaced by (Subscript[M, (Subscript[B, D])]^(Subscript[B, U]))(I)(Subscript[M, (Subscript[B, U])]^(Subscript[B, V]))(M)(Subscript[M, (Subscript[B, V])]^(Subscript[B, C]))(I) where (Subscript[M, (Subscript[B, U])]^(Subscript[B, V]))(M) is a "diagonal" matrix.  This amounts to changing our bases from (Subscript[B, C])↦(Subscript[B, V])↦(Subscript[B, U])↦(Subscript[B, D]) where the first and last change of basis is an identity transformation and the center change of basis is a "diagonal" transformation.

Note:  In general, we think about the bases (Subscript[B, C]) and (Subscript[B, D]) as the standard bases with the proper dimension.  That is, we think of (Subscript[B, C])=(Subscript[E, n]) and (Subscript[B, D])=(Subscript[E, m]).

We are now ready to give our second proof of the singular value decomposition theorem.

Theorem: Singular Value Decomposition (SVD) Theorem
Any m x n real valued matrix M can be decomposed into the product

M=U·Λ·(V^T)

where (Subscript[U, m x m]) and (Subscript[V, n x n]) are orthonormal matrices and (Subscript[Λ, m x n]) is a "diagonal" matrix whose entries are the square roots of the eigenvalues of the matrix (M^T)M.

□ Proof.
Let M, (Subscript[{(Subscript[σ, i])}, i=1]^n), (Subscript[B, V]), and (Subscript[B, U]) be as above.  Define U=[(Subscript[u, 1]),…,(Subscript[u, m])]  and V=[(Subscript[v, 1]),…,(Subscript[v, n])].  By construction, we know that U and V are orthonormal.  Also by construction, we have that 

M V=[M (Subscript[v, 1])…M (Subscript[v, r])   0…0]=[(Subscript[σ, 1])(Subscript[u, 1])(Subscript[…σ, r])(Subscript[u, r])  0…0]

By comparison, we have

U·Λ=[(Subscript[u, 1]),…,(Subscript[u, r]) ,(Subscript[u, r+1]),(Subscript[…u, m])]·

=[(Subscript[σ, 1])(Subscript[u, 1]),…,(Subscript[σ, r])(Subscript[u, r]) , 0,…,0]

=M·V

So U·Λ=M·V which implies that U·Λ·(V^T)=M·V·(V^T).  However, V is an orthonormal matrix, so M=U·Λ·(V^T). ■

Example 1(j)
Given the matrix M=(({{1, 2}, {3, -1}, {-2, -3}})) , we found the orthonormal matrix V=(({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})) and the "diagonal" matrix Λ=(({{(Sqrt[19]), 0}, {0, (Sqrt[9])}, {0, 0}})) and (Subscript[B, U])={(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])}.  
Define U=[(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])]=(({{(3/(Sqrt[38])), (1/(3 (Sqrt[2]))), (11/(3 (Sqrt[19])))}, {(Sqrt[(2/19)]), -((2 (Sqrt[2]))/3), (1/(3 (Sqrt[19])))}, {-(5/(Sqrt[38])), -(1/(3 (Sqrt[2]))), (7/(3 (Sqrt[19])))}}))

Now we may verify that U·Λ·(V^T)=M

U·Λ·(V^T)=(({{(3/(Sqrt[38])), (1/(3 (Sqrt[2]))), (11/(3 (Sqrt[19])))}, {(Sqrt[(2/19)]), -((2 (Sqrt[2]))/3), (1/(3 (Sqrt[19])))}, {-(5/(Sqrt[38])), -(1/(3 (Sqrt[2]))), (7/(3 (Sqrt[19])))}}))·(({{(Sqrt[19]), 0}, {0, (Sqrt[9])}, {0, 0}}))·(((({{(1/(Sqrt[2])), -(1/(Sqrt[2]))}, {(1/(Sqrt[2])), (1/(Sqrt[2]))}})))^T)=(({{1, 2}, {3, -1}, {-2, -3}}))=M

Before concluding example one, let us consider a second example which will help us visualize the conclusion of example one.

Example 2
Let the matrix M=(({{1, 3}, {-5, 7}})).  Since A is a real valued matrix, there is a singular value decomposition of A.  

The diagonal matrix "Λ = "(({{8.81913190966076`, 0.`}, {0.`, 2.4945765893240006`}})) contains the singular values associated with A.

The matrices V=" "(({{-0.5257311121191336`, 0.85065080835204`}, {0.8506508083520399`, 0.5257311121191336`}})) and U=(({{0.22975292054736118`, 0.9732489894677303`}, {0.9732489894677301`, -0.2297529205473613`}})) are orthonormal matrices.

Putting these matrices together, we find that U·Λ·(V^T)=""(({{1.`, 3.`}, {-5.`, 6.999999999999999`}})).

Visualizing the Singular Value Decomposition M
This decomposition can be thought of in the following way.  Begin with the unit circle and the rotate the unit circle according to the transformation (V^T).  Now stretch out the eigenvectors of (V^T) according to the singular values found in Λ.  Finally, rotate the ellipse resulting from Λ·(V^T) according to the transformation found in U.  Thus, M can be broken up into a rotation: (V^T), component-wise stretch: Λ, and a second rotation: U.  (see figure)

Figure: Visualizing the singular value decomposition in Example 2

Example 1(k)
We have found the singular value decomposition of M.  We find it convenient to begin and end our transformation with respect to the standard basis.  That is, M=(Subscript[M, (Subscript[E, 3])]^(Subscript[E, 2]))(M). We have found the orthonormal bases (Subscript[B, V]) and (Subscript[B, U]) such that (Subscript[M, (Subscript[B, U])]^(Subscript[B, V]))(M)=Λ.  Furthermore, we have that (Subscript[M, (Subscript[B, V])]^(Subscript[E, 2]))((Subscript[I, 2]))=(V^T) and (Subscript[M, (Subscript[E, 3])]^(Subscript[B, U]))((Subscript[I, 3]))=U.  Therefore,

M=(Subscript[M, (Subscript[E, 3])]^(Subscript[E, 2]))(M)=(Subscript[M, (Subscript[E, 3])]^(Subscript[B, U]))((Subscript[I, 3]))·(Subscript[M, (Subscript[B, U])]^(Subscript[B, V]))(M)·(Subscript[M, (Subscript[B, V])]^(Subscript[E, 2]))((Subscript[I, 2]))=U·Λ·(V^T) as desired.

This transformation may be visualized in the following manner:
Begin with the unit circle in (ℝ^2) and rotate it as described in (V^T).  Then project our rotated unit circle onto the x-y plane in (ℝ^3), but leave the circle centered at the origin.  From here, "stretch" the circle into an ellipse component-wise  Here we are stretching with respect to the basis G as described in Λ.  Finally, rotate and project our ellipse onto the plane formed by the span⟨(Subscript[u, 1]),(Subscript[u, 2])⟩ while leaving our ellipse centered at the origin as described by U.

And so we have found a conceptually simple representation of M.  Long ago, we sought a conceptually simply representation of polynomials and we found it by factoring - here we have factored matrices.  After finding the factorization of a polynomial, we used it to find the roots of the polynomial - here we use the factorization to do . . .

(4.)  Matrix Reduction Through Singular Value Decomposition

We begin our search for applications of singular value decomposition by looking at a way we can reduce a matrix.

Theorem: Matrix Reduction
Let (Subscript[M, m x n]) be as above.  Then all the information contained in M can be represented by storing two matrices (hopefully smaller that M) and one vector.

□ Proof.
M=U·Λ·(V^T)
=(({{(Subscript[U, r]), (Subscript[U, m-r])}}))·(({{(Subscript[Σ, r]), (Subscript[0, r x (n-r)])}, {(Subscript[0, (m-r) x r]), (Subscript[0, (m-r) x (n-r)])}}))·(({{(Subscript[V, r]^T)}, {(Subscript[V, n-r]^T)}}))
=(Subscript[U, r])·(Subscript[Σ, r])·(Subscript[V, r]^T)
Therefore M=(Subscript[U, r])·(Subscript[Σ, r])·(Subscript[V, r]^T).  ■

So, we have two options for storing M.  (a.)  Store the m x n entries of M.  Or (b.)  Store the (m x r)+r+(n x r) entries in (Subscript[U, r]) ,(Subscript[V, r]^T), and diag⟨(Subscript[Σ, r])⟩.  Clearly, we would want to choose which ever option required less storage.

An additional advantage of this new representation of M is that it has a known inverse (the Moore-Penrose Inverse) (SuperPlus[M])=(Subscript[V, r])((Subscript[Σ, r])^(-1))(Subscript[U, r]^T).  To see that this is the inverse, consider 

((Subscript[U, r])·(Subscript[Σ, r])·(Subscript[V, r]^T))·((Subscript[V, r])·(Subscript[Σ, r]^(-1))·(Subscript[U, r]^T))=(Subscript[U, r])·(Subscript[Σ, r])·(Subscript[V, r]^T)·(Subscript[V, r])·((Subscript[Σ, r])^(-1))·(Subscript[U, r]^T)
                                                           =(Subscript[U, r])·(Subscript[Σ, r])·(Subscript[I, r])·((Subscript[Σ, r])^(-1))·(Subscript[U, r]^T)
                                                      =(Subscript[U, r])·(Subscript[I, r])·(Subscript[U, r]^T)
                                                  =(Subscript[I, m]^(rank=r))

That is, we get out (Subscript[I, m]^(rank=r))=(Subscript[(({{(Subscript[I, r]), 0}, {0, 0}})), m x m]).

Thus, given an SVD of M, we can reduce the problem so as to require less storage, and have an pseudo-inverse thrown in for free.

Extending Matrix Reduction
Suppose that M=U·Λ·(V^T) is the singular value decomposition of M.  If the rank of M is r, then we have already shown that M≅(Subscript[U, r])·(Subscript[Σ, r])·(Subscript[V, r]^T).  Now the singular values of M are ordered so that (Subscript[σ, 1])≥(Subscript[σ, 2])≥...≥(Subscript[σ, r])>0.  So, it is no great leap to further approximate M by disregarding the r-k smallest singular values.  Thus, we have that M≅(Subscript[M, k])=(Subscript[U, k])·(Subscript[Σ, k])·(Subscript[V, k]^T) where (Subscript[Σ, k]) is defined as above.

(5.)  Numerical Rank Approximation

Theoretically, the singular value decomposition of a matrix M reveal the rank of a matrix.  This is because the rank of M is the number of non-zero singular values in the singular value decomposition of M.  However, in practice, the rank is not so readily apparent.  In practice, finding the singular value decomposition of M requires inputting M into a computer and thus introducing computer round-off to the problem.  Computer round-off in and of itself would not be a problem except that the space of matrices with full rank is an open and dense space.  That is, (OverTilde[M]) - the numerical representation of our matrix M - is most likely a matrix with full rank.  And so rank(M)≤rank((OverTilde[M])).

Definition:
The matrix two-norm is defined to be the norm induced by the vector two-norm:  Let A∈(ℝ^(m x n)), then

 ||| A |||2 = (max_(x≠0)) ((||A x (Subscript[||, 2]))/(|| x (Subscript[||, 2]))) where x∈(ℝ^n)

Theorem
Let A∈(ℝ^(m x n)) have singular values (Subscript[σ, 1])≥...≥0.  Then ||| A |||2 = (Subscript[σ, 1])

□ Proof.
We need to show that (max_(x≠0)) ((||A x (Subscript[||, 2]))/(|| x (Subscript[||, 2])))=(Subscript[σ, 1]).  Well, we begin by observing that if (Subscript[v, 1]) is the eigenvector associated with (Subscript[σ, 1]), then 

(max_(x≠0)) ((||A (Subscript[v, 1]) (Subscript[||, 2]))/(|| (Subscript[v, 1]) (Subscript[||, 2])))=(Subscript[σ, 1])((|| (Subscript[u, 1])(Subscript[||, 2]))/(|| (Subscript[v, 1])(Subscript[||, 2])))=(Subscript[σ, 1])

so (max_(x≠0)) ((||A x (Subscript[||, 2]))/(|| x (Subscript[||, 2])))≥(Subscript[σ, 1]).  Now we must show that no other vector is "stretched" by more than (Subscript[σ, 1]).  Let x∈(ℝ^n).  Then x is a linear combination of the (Subscript[v, i])'s such that x=(Subscript[c, 1])(Subscript[v, 1])+...+(Subscript[c, n])(Subscript[v, n]).  Since the (Subscript[v, i])'s are orthonormal, we have that ||x(Subscript[||, 2]^2)=|(Subscript[c, 1])((|)^2)+...+|(Subscript[c, n])((|)^2).  Now the product A x=A ((Subscript[c, 1])(Subscript[v, 1])+...+(Subscript[c, n])(Subscript[v, n])).  This in turn is equal to (Subscript[σ, 1])(Subscript[c, 1])(Subscript[u, 1])+...+(Subscript[σ, r])(Subscript[c, r])(Subscript[u, r]) where r is the number or non-zero singular values of A.  (Alternately, rank(A)=r).  Since the (Subscript[u, i])'s are orthonormal, we have that

||A x(Subscript[||, 2]^2)=|(Subscript[σ, 1])(Subscript[c, 1])((|)^2)+...+|(Subscript[σ, r])(Subscript[c, r])((|)^2)≤(Subscript[σ, 1]^2) (|(Subscript[c, 1])((|)^2)+...+|(Subscript[c, r])((|)^2))≤(Subscript[σ, 1]^2) (|(Subscript[c, 1])((|)^2)+...+|(Subscript[c, n])((|)^2))=(Subscript[σ, 1]^2) || x (Subscript[||, 2]^2)

Thus ((||A x (Subscript[||, 2]))/(|| x (Subscript[||, 2])))≥(Subscript[σ, 1]) as desired.  ■  (Watkins 409).

The two-norm gives us a way to measure the distance between two matrices.  So, with that in mind, consider the following lemma.

Lemma
Let A∈(ℝ^(m x n)) where rank(A)=r<min{m,n}.  Then for all ϵ>0, there is a matrix (Subscript[A, ϵ]) such that rank((Subscript[A, ϵ]))=min{m,n} and (Subscript[A, ϵ]) is within ϵ>0 of A as measured by the two-norm

□ Proof.
Let A be as above and let ϵ>0 be given such that (Subscript[σ, r])>ϵ>0.  Now, by the singular value decomposition theorem, we have that A=U·Λ·(V^T).  Define (Subscript[A, ϵ]) to be (Subscript[A, ϵ])=U·(Subscript[Λ, ϵ])·(V^T) where (Subscript[Λ, ϵ]) is a diagonal matrix whose diagonal entries correspond to the diagonal entries of Λ minus ϵ.  That is, (Subscript[Λ, i,i])↦(Subscript[Λ, i,i])-ϵ ,for i=1,2,3,...,min{m,n}.  Certainly (Subscript[A, ϵ]) has full rank.  Now, A-(Subscript[A, ϵ])=U·Λ·(V^T)-U·(Subscript[Λ, ϵ])·(V^T)=U·(Λ-(Subscript[Λ, ϵ]))·(V^T)has a spectral radius of ϵ>0 and so A and (Subscript[A, ϵ]) are within ϵ>0 as measured with the two-norm.  ■

Theorem
Let A∈(ℝ^(m x n)) where rank(A)=r>0.  Let A=U·Λ·(V^T) be the singular value decomposition of A where V=[(Subscript[v, 1]),...,(Subscript[v, n])] and U=[(Subscript[u, 1]),...,(Subscript[u, m])].  Furthermore, for k=1,2,3,...,r-1, define (Subscript[A, k]) be the extended matrix reduction of A.  Then rank((Subscript[A, k]))=k and (Subscript[σ, k+1])= ||| A-(Subscript[A, k])||(Subscript[|, 2]) = min{ ||| A-B ||(Subscript[|, 2]) : rank(B) = k}.  That is, (Subscript[A, k]) is the matrix with rank k that is closest to A when measured with the two norm.

□ Proof.
Clearly the rank of (Subscript[A, k]) is k  and the observation that A-(Subscript[A, k])=U·(Λ-(Subscript[Λ, k]))·(V^T) makes it equally clear that |||A-(Subscript[A, k])||(Subscript[|, 2]) = (Subscript[σ, k+1]).  It only remains to show that (Subscript[A, k]) is the matrix with rank k that is closest to A.  Let B∈(ℝ^(m x n)) be a matrix with rank k.  First note that dim(Null(B))=n-k.  Also, note that dim(Span⟨(Subscript[v, 1]),...,(Subscript[v, k+1])⟩)=k+1.  Now Null(B)⊂(ℝ^n) and each (Subscript[v, 1]),...,(Subscript[v, k+1])∈(ℝ^n).  Since dim(Null(B))+dim(Span⟨(Subscript[v, 1]),...,(Subscript[v, k+1])⟩)=(n-k)+(k+1)=n+1>n, there exists x∈Null(B)⋂ Span⟨(Subscript[v, 1]),...,(Subscript[v, k+1])⟩.  We assume that ||x(Subscript[||, 2])=1.  Since x∈Span⟨(Subscript[v, 1]),...,(Subscript[v, k+1])⟩, there exist scalars (Subscript[c, 1]),...,(Subscript[c, k+1]) such that x=(Subscript[c, 1])(Subscript[v, 1])+...+(Subscript[c, k+1])(Subscript[v, k+1]).  Since (Subscript[v, 1]),...,(Subscript[v, k+1]) are orthonormal, we have that |(Subscript[c, 1])((|)^2)+...+|(Subscript[c, k+1])((|)^2)=||x(Subscript[||, 2]^2)=1.  Since x∈Null(B), we have that B x=(Subscript[0, (ℝ^m)]).  Thus,

(A-B)x=A x
           =A (Subscript[c, 1])(Subscript[v, 1])+...+A (Subscript[c, k+1])(Subscript[v, k+1])
        =(Subscript[c, 1])A (Subscript[v, 1])+...+(Subscript[c, k+1])A (Subscript[v, k+1])

recalling that A (Subscript[v, i])=(Subscript[σ, i])(Subscript[u, i]) we continue our equality with

                   =(Subscript[c, 1])(Subscript[σ, 1])(Subscript[u, 1])+...(Subscript[c, k+1])(Subscript[σ, k+1])(Subscript[u, k+1])
             ≥(Subscript[σ, k+1])((Subscript[c, 1])(Subscript[u, 1])+...+(Subscript[c, k+1])(Subscript[u, k+1]))

Since the (Subscript[u, i]) are orthonormal, by taking the two-norm of both sides, we have that: 

|||(A-B)x||(Subscript[|, 2])≥(Subscript[σ, k+1])
 
Finally, since x is of unit length, we have that

(Subscript[σ, k+1])≤ ||| A - B |||2 || x ||2 = ||| A - B |||2

And thus the theorem is proved.  ■

As a direct consequence of the preceding theorem, we have the following corollary

Corollary
Suppose A∈(ℝ^(m x n)) is such that rank(A)=r=min{m,n}.  Let (Subscript[σ, 1])≥...≥(Subscript[σ, r])>0 be the singular values of A.  Then if B∈(ℝ^(m x n)) satisfies |||A-B|||2 < (Subscript[σ, r]), then B also has rank r.

Summarizing the conclusions of the preceding lemma, theorem, and corollary, we have that every rank deficient matrix has a matrix of full rank arbitrarily close to it.  Furthermore, every matrix with full rank is surrounded by matrices with full rank.  So, the set of full rank matrices is an open and dense subset of (ℝ^(m x n)).  Thus, we have shown that "most" matrices have full rank.

If a matrix is rank deficient, we have shown that an small perturbation will most likely transform our matrix into one with full rank.  It follows that any type of calculation or storage that introduces even the slightest amount of uncertainty will make the resulting matrix's theoretical rank virtually impossible to detect.  However, it may still be possible to find the approximate numerical rank of the matrix.

Suppose that for some very small ϵ>0, the matrix A has singular values such that (Subscript[σ, 1])≥...≥(Subscript[σ, k])≫ϵ>(Subscript[σ, k+1])≥...  Mathematically, this can only take place if (a.) there exists a matrix B such that rank(B)=k and |||A-B|||2 < ϵ and (b.) all matrices C that have rank k-1 are much farther from A than ϵ>0.  That is |||A-C|||2 ≫ϵ.  Given these conditions, then we say that the numerical rank of A is k.

Of course, in order to determine the numerical rank, we would need to have some idea of what ϵ>0 would be reasonable.  One simple way of determining ϵ is to look at the precision to which calculations are made.  (For example, see the IEEE standards).  If, for example, computer round-off takes place at eps=(10^(-16)), then it would be reasonable to choose ϵ=(10^(-15)).  That is, we can safely assume that any singular value smaller than (10^(-15)) is zero.

While cases do occur where there is a gap in the singular values, it is often the case that no obvious gap can be found, but rather that we will make our rank approximation based on how much precision we require out of the results.

Example 3(a)
Let (Subscript[P, 128 x 128]) be a real valued matrix whose significance will be discussed later.  The 128 singular values of P are as follows: 

({{3601.2`, 1157.6`, 830.9`, 569.5`, 366.20000000000005`, 300.70000000000005`, 255.2`, 243.7`}, {197.1`, 163.39999999999998`, 151.9`, 131.89999999999998`, 119.`, 103.5`, 92.39999999999999`, 79.2`}, {71.8`, 66.4`, 61.800000000000004`, 59.3`, 48.300000000000004`, 47.5`, 43.3`, 41.9`}, {38.800000000000004`, 38.199999999999996`, 33.4`, 32.099999999999994`, 29.4`, 28.7`, 28.`, 25.`}, {23.400000000000002`, 23.400000000000002`, 22.5`, 20.6`, 18.9`, 18.2`, 18.1`, 17.6`}, {16.8`, 15.7`, 15.`, 14.5`, 14.200000000000001`, 12.9`, 12.8`, 12.5`}, {12.1`, 11.700000000000001`, 11.5`, 11.1`, 10.6`, 10.200000000000001`, 9.9`, 9.6`}, {9.5`, 9.2`, 9.1`, 8.7`, 8.4`, 8.200000000000001`, 7.9`, 7.5`}, {7.2`, 7.2`, 6.8999999999999995`, 6.7`, 6.5`, 6.3`, 6.2`, 6.1000000000000005`}, {5.7`, 5.6`, 5.3`, 5.3`, 5.1000000000000005`, 5.`, 4.7`, 4.7`}, {4.4`, 4.4`, 4.3`, 4.`, 3.9`, 3.8`, 3.7`, 3.6`}, {3.5`, 3.5`, 3.2`, 3.2`, 3.1`, 3.`, 2.9`, 2.8`}, {2.7`, 2.6`, 2.5`, 2.5`, 2.4`, 2.3`, 2.1`, 2.`}, {1.8`, 1.8`, 1.8`, 1.6`, 1.5`, 1.4`, 1.3`, 1.2`}, {1.2`, 1.1`, 1.`, 0.9`, 0.9`, 0.8`, 0.7`, 0.7`}, {0.5`, 0.4`, 0.4`, 0.3`, 0.3`, 0.2`, 0.1`, 0.`}})

Note: The singular values have been chopped at one decimal place.

While there are no major gaps in the singular values, it is obvious that the first few rows of values are of much more importance than the later rows.  So, it would be reasonable to expect a low rank approximation of rank eight, sixteen, twenty-four, or thirty-two to give satisfactory results.

(6.)  Image Compression

The Scenario
We are constructing a new satellite to orbit around Mars.  The cameras and sensors on our satellite are to take one image per second.  Then the images are transmitted back to Earth for processing.  The problem we face is this: the images the satellite will collect are very large - say 10MB each, but we can transmit at a much slower rate - say 1MB per second.  That is, in the ten seconds it takes us to transmit the first image, we will have collected nine more images.  So, we must find a way of compressing the images to a tenth of their original size without losing the major characteristics of the original image.  At our disposal, we have the virtually limitless power and storage of the satellites onboard computers.  Yet we need an algorithm which will compress the images without destroying the original - a simple shrinking of the image will not suffice.  Enter singular value decomposition.

Our goal is to use singular value decomposition to make a close approximation of the original image.  To do this, we being by viewing the matrix representation of the image.  That is, we break the image into an m x n lattice with the corresponding matrix A.  That is, IMAGE ↦ (Subscript[A, m x n]).  To (Subscript[a, i, j]), we assign a value corresponding to the color of the (((i,j))^th) region of the image.  Now, we can find the singular value decomposition of A and then consider the approximations of A: (Subscript[A, 1]), (Subscript[A, 2]), (Subscript[A, 3]), ...  We hope that for some small k, (Subscript[A, k]) sufficiently approximates A.

Before going too far, let us consider a simple example.

Example 4
Lets begin with a very simple example.  Consider the gray-scale scheme between zero and one where 0↦ BLACK and 1↦WHITE and values between 0 and 1 are assigned a greys.  Now, it seems logical to start with a 𝒮uper example, so let the matrix A represent be the matrix representation of the letter S.



In order to recall A, we must store 9 x 6 entries, or 54 entries.  Now, we find that the singular value decomposition of A is: (rounding to two places)

A=U·Λ·(V^T)=(Subscript[(({{0, 0, 0, 0, 1, 0, 0, 0, 0}, {0.60, 0.41, 0, -0.40, 0, -0.29, -0.29, -0.39, 0}, {0.12, -0.41, 0.45, -0.51, 0, -0.18, -0.18, 0.541, 0}, {0.12, -0.41, 0.45, -0.14, 0, 0.34, 0.34, -0.61, 0}, {0.48, 0, -0.45, -0.25, 0, 0.45, 0.45, 0.33, 0}, {0.12, 0.441, 0.45, 0.20, 0, 0.64, -0.36, 0.20, 0}, {0.12, 0.41, 0.45, 0.20, 0, -0.36, 0.64, 0.20, 0}, {0.60, 0.41, 0, 0.65, 0, -0.16, -0.16, 0.06, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1}})), 9 x 9])·(Subscript[(({{2.65, 0, 0, 0, 0, 0}, {0, 1.73, 0, 0, 0, 0}, {0, 0, 1.41, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}})), 9 x 6])·(((({{0, 0, 0, 0, 1, 0}, {0.32, -0.71, 0.63, 0, 0, 0}, {0.63, 0, -0.32, 0.71, 0, 0}, {0.63, 0, -0.32, -0.71, 0, 0}, {0.32, 0.71, 0.63, 0, 0, 0}, {0, 0, 0, 0, 0, 1}})))^T)

Now, let us consider the images associated with (Subscript[A, 1]), (Subscript[A, 2]), (Subscript[A, 3]), (Subscript[A, 4]),(Subscript[A, 5]), and (Subscript[A, 6])
(Subscript[A, 1])=(Subscript[U, 1])·(Subscript[D, 1])·(Subscript[V, 1]^T)=(Subscript[(({{0}, {0.60}, {0.12}, {0.12}, {0.48}, {0.12}, {0.12}, {0.60}, {0}})), 9 x 1])·(Subscript[(2.65), 1 x 1])·(Subscript[(({{0}, {0.32}, {0.63}, {0.63}, {0.32}, {0}})), 6 x 1]^T)has dimension 9 x 6 and the associated image: 
In order to view this image, I must store the entries from (Subscript[U, 1]), (Subscript[D, 1]), and (Subscript[V, 1]) which means storing (9 x 1) + 1+(6 x 1)=16 entries.  This image certainly will not be mistaken for an 𝒮, so let us try (Subscript[A, 2]).

(Subscript[A, 2])=(Subscript[U, 2])·(Subscript[D, 2])·(Subscript[V, 2]^T)=(Subscript[(({{0, 0}, {0.60, 0.41}, {0.12, -0.41}, {0.12, -0.41}, {0.48, 0}, {0.12, 0.441}, {0.12, 0.41}, {0.60, 0.41}, {0, 0}})), 9 x 2])·(Subscript[(({{2.65, 0}, {0, 1.73}})), 2 x 2])·(Subscript[(({{0, 0}, {0.32, -0.71}, {0.63, 0}, {0.63, 0}, {0.32, 0.71}, {0, 0}})), 6 x 2]^T)again has dimension 9 x 6 and associated image: 
We must now store the entries from (Subscript[U, 2]), (Subscript[D, 2]), and (Subscript[V, 2]) which means storing (9 x 2) + 2+(6 x 2)=32 entries.  (Note: we only save two entries from (Subscript[D, 2]) since we know it is a diagonal matrix and thus all its off-diagonal elements are zero.)  The image representation of (Subscript[A, 2]) contains a recognizable 𝒮, so we could choose to accept this approximation of A or find (Subscript[A, 3]).  Suppose that we were not satisfied with (Subscript[A, 2])=(Subscript[U, 2])·(Subscript[D, 2])·(Subscript[V, 2]^T) as an approximation for A, we can try (Subscript[A, 3]).
(Subscript[(({{0, 0, 0}, {0.60, 0.41, 0}, {0.12, -0.41, 0.45}, {0.12, -0.41, 0.45}, {0.48, 0, -0.45}, {0.12, 0.441, 0.45}, {0.12, 0.41, 0.45}, {0.60, 0.41, 0}, {0, 0, 0}})), 9 x3])·(Subscript[(({{2.65, 0, 0}, {0, 1.73, 0}, {0, 0, 1.41}})), 3 x 3])·(Subscript[(({{0, 0, 0}, {0.32, -0.71, 0.63}, {0.63, 0, -0.32}, {0.63, 0, -0.32}, {0.32, 0.71, 0.63}, {0, 0, 0}})), 6 x 3]^T)has  associated image: 

This last image is the same image that we began with.  To get this result, we store the entries from (Subscript[U, 3]), (Subscript[D, 3]), and (Subscript[V, 3]) which means storing (9 x 3) +3+(6 x 3)=48 entries.  Thus we need only store 48 entries to maintain an image with 54 entries - a compression of 11%.  Furthermore, we could have approximated A with the recognizable image (Subscript[A, 2]) which would have been a compression of 40%.   

Now, it may not be clear that A=(Subscript[A, 3])=(Subscript[A, 4])=(Subscript[A, 5])=(Subscript[A, 5])=(Subscript[A, 6]).  In order to see this, consider (Subscript[A, 4])
(Subscript[A, 4])=(Subscript[U, 4])·(Subscript[D, 4])·(Subscript[V, 4]^T)=(({{(Subscript[U, 3]), (Subscript[u, 4])}}))·(({{(Subscript[D, 3]), (Subscript[0, 3 x 1])}, {(Subscript[0, 1 x 3]), (Subscript[0, 1 x 1])}}))·(({{(Subscript[V, 3]^T)}, {(Subscript[v, 4]^T)}})) where (Subscript[u, 4]) is the (4^th) column in U and (Subscript[v, 4]) is similarly defined.  By expanding the product out, we see that (Subscript[A, 4])=(Subscript[U, 3])·(Subscript[D, 3])·(Subscript[V, 3]^T)=(Subscript[A, 3]).  This makes the remainder of the equality A=(Subscript[A, 3])=(Subscript[A, 4])=(Subscript[A, 5])=(Subscript[A, 5])=(Subscript[A, 6]) obvious.  Mathematically, it is makes sense that it would only require a reduction to (Subscript[A, 3]) to return to a matrix of rank 3.

To summarize this example, given a matrix whose image representation was the letter 𝒮, we were able to approximate our image with a 40% reduction in storage and we were able to recapture our original image with a savings of 11% in storage.  These savings seems relatively significant, so we ask ourselves what the cost was?  In order to compress the image of my 𝒮, I had to begin with a singular value decomposition which is computationally intensive (the decomposition required 2496 flops in Matlab).  Then, I had to look at each image: (Subscript[A, 1]), (Subscript[A, 2]), and (Subscript[A, 3]) to decide whether I was satisfied with my result.  That is, I actually had to compute multiply out each (Subscript[U, i])·(Subscript[D, i])·(Subscript[V, i]^T) for i=1,2,3 which requires (i^2)(2 m i+n) flops.  Finally, I must actually view the image to determine whether it is satisfactory.  This requires more computational time.  All together, the results of image compression through singular value decomposition are good in this example, however, the results come at the cost of some serious computations.

Note: In calculating the number of entries which must be stored for each (Subscript[A, i]), I glossed over the fact that we must save the dimensions of each matrix as well.  This would have made a noticeable difference in the previous example, however, it would be negligible if the image we were working with was much larger.

The previous example captures the theory of image compression through singular value decompression, however its simplicity does not give use the confidence in the method or viability of the process needed in order to ship it to Mars.  So, extend the methodology in the following way.  First, we are no longer satisfied with grey-scale images, so let us break the color spectrum up into - say 64 - pieces.  Thus, we assign to each integer between zero and sixty-three a specific color and work with our images accordingly.  Suppose we now have an image, instead of finding its singular value decomposition, we will simply find its singular values.  Given these singular values, we can decide how much of U, Λ, and V we want to find.

Consider another example.

Example 3(b)
Let (Subscript[P, 128 x 128]) be the aforementioned matrix.  P is the matrix representation of a U.S. penny.  The singular values of P were as listed in example 3(a).

The image corresponding to P is: 

Looking at the important features of the original image, we can clearly read the U.S. motto, "IN GOD WE TRUST," as well as "IBERTY" and "1984,"  So, we will be satisfied with any approximation which leaves these features legible.

Choosing a few miscellaneous steps along the way, we have (Subscript[P, 1]),(Subscript[P, 8]), (Subscript[P, 16]), and (Subscript[P, 32])

({{, }, {, }})
 
The key features of the penny are discernible in (Subscript[P, 16]) and most coin-shops would pay money for the image in (Subscript[P, 32]).  What is the compression here - for each additional singular value, we must store 257 more entries.  This means that (Subscript[P, 16]) requires the storage of 4112 entries and (Subscript[P, 32]) requires the storage of 8224 entries.  Contrast this with P itself which contains 16384 entries.  Thus, a 75% compression would suffice for maintaining the basic outlines of the key features while a 50% compression is good enough to mint.

Now that we understand image compression, it is finally feasible to give an interpretation of the singular values of a matrix.  Consider the matrix M associated with some image and the corresponding singular values {(Subscript[σ, 1]),…,(Subscript[σ, r])}.  At each iteration of our image compression process, we are adding on another layer of information where the layer is formed from the matrix (Subscript[σ, i])(Subscript[u, i])(Subscript[v, i]^T).  Here (Subscript[σ, i]) gives the intensity of the layer.  So, the (i^th) approximation is 

(Subscript[M, i])=(Subscript[σ, 1])(Subscript[u, 1])(Subscript[v, 1]^T)+…+(Subscript[σ, i])(Subscript[u, i])(Subscript[v, i]^T) 

and this iteration contains (((∑_(j=1)^i)(Subscript[σ, j]))/((∑_(j=1)^r)(Subscript[σ, j]))) percent of the total information contained in the original image.

Appendix (1.)  Definitions

Definition: Scalar Field
A set 𝕂⊂ℂis called a scalar field if
(1.)  0,1∈𝕂
(2.)  α+β∈𝕂 and -α∈𝕂 whenever α,β∈𝕂
(3.)  α·β∈𝕂 whenever α,β∈𝕂 and (1/α)∈𝕂 whenever α∈𝕂∼{0}
(4.)  (OverBar[α])∈𝕂 whenever α∈𝕂

Definition: Vector Space
A set 𝒱 is a vector space over a scalar field 𝕂 if it satisfies the following nine conditions. 
(1.)  For all α∈𝕂 and all u,v∈𝒱, we have α u∈𝒱 and u+v∈𝒱
(2.)  For all u,v,w∈𝒱, we have (u+v)+w=u+(v+w)  (associative property of addition for vectors)
(3.)  For all u,v∈𝒱, we have u+v=v+u (commutative property of addition)
(4.)  There is an element of 𝒱, denoted by (Subscript[0, 𝒱]), such that (Subscript[0, 𝒱])+u=u+(Subscript[0, 𝒱])=u for all u∈𝒱
(5.)  For each element u∈𝒱, there exists an element in 𝒱, denoted by -u, such that u+(-u)=(Subscript[0, 𝒱])
(6.)  For all α∈𝕂 and all u,v∈𝒱, we have α(u+v)=α u+α v (distributive property for scalars)
(7.)  For all α,β∈𝕂 and all u∈𝒱, we have (α+β) u=α u+β u
(8.)  For all α,β∈𝕂 and all u∈𝒱, we have (α β) u=α (β u) (associativity of scalars)
(9.)  For all u∈𝒱, there exists a scalar in 𝕂, denoted as 1, such that 1·u=u

Definition: Linear Transformation
Let 𝒰 and 𝒱 be vector spaces over the same scalar field 𝕂.  A transformation S:𝒰 ↦ 𝒱 is linear if S((Subscript[α, 1])(Subscript[u, 1])+(Subscript[α, 2])(Subscript[u, 2]))=(Subscript[α, 1])S((Subscript[u, 1]))+(Subscript[α, 2])S((Subscript[u, 2])) for all (Subscript[α, 1]),(Subscript[α, 2])∈𝕂, (Subscript[u, 1]),(Subscript[u, 2])∈𝒰

Definition: Inner Product
Let 𝒰 be a vector space over the scalar field 𝕂.  Let the inner product ⟨·,·⟩:𝒰×𝒰 ↦ 𝕂 be a function defined by the following two conditions:
1.)  ⟨α u+β v,w⟩=α⟨u,w⟩+β⟨v,w⟩ for all α,β∈𝕂 and for all u,v,w∈𝒰
2.)  ⟨u,v⟩=(OverBar[⟨v,u⟩]) for all u,v∈𝒰, where (OverBar[⟨v,u⟩]) denotes the complex conjugate of ⟨v,u⟩
Additionally, for the purposes of this paper, we place the following two restrictions on our inner product
3.)  ⟨u,v⟩=0 for all v∈𝒰 implies u=(Subscript[0, 𝒰])
4.)  ⟨v,v⟩≥0 for all v ∈𝒰
An inner product satisfying these four conditions is said to be a "positive definite inner product" and an inner product satisfying conditions (1), (2), and (4) is said be a "non-negative inner product."

Definition: Vector Norm
Let u,v∈(𝕂^n).  A function ||•||:(𝕂^n)↦𝕂 satisfying the following four conditions is referred to as a vector norm:
1.)  ||u||≥0
2.)  ||u||=0 if and only if u=(Subscript[0, (𝕂^n)])
3.)  ||α u||=|α|||u|| for all α∈𝕂
4.)  ||u+v||≤||u||+||v||.  (triangle inequality)

Definition: Orthonormal
Let 𝒰 be a vector space with a positive definite inner product.  The vectors (Subscript[u, 1]),...,(Subscript[u, m])∈𝒰 are said to be orthonormal if they satisfy the following two conditions:
1.)  ⟨(Subscript[u, i]),(Subscript[u, j])⟩=0 if and only if i≠j  (orthogonal)
2.)  ||(Subscript[u, i])||=1, for i=1,2,3,...,m  (normal)
Note: most often we refer to orthonormal bases, that is and orthonormal set of vectors {(Subscript[u, 1]),...,(Subscript[u, n])} such that the span{(Subscript[u, 1]),...,(Subscript[u, n])}=𝒰.
An alternate, but sometimes more convenient definition of orthonormal is as follows.  Let Q be an n x n matrix with linearly independent columns.  If (Q^T)Q=(Q^(-1))Q=I, that is (Q^T)=(Q^(-1)), then the columns of Q are orthonormal

Definition: Transpose (alternate definition)
The transpose of a matrix A is alternately defined as a matrix (A^T) such that ⟨A v,w⟩=⟨v,(A^T)w⟩ for any v,w∈𝒰, where 𝒰 is a vector space with a non-degenerate inner product.

Definition: Symmetric
A matrix A is symmetric if A=(A^T).  Alternately, we refer to a transformation as symmetric if ⟨A v,w⟩=⟨v,A w⟩ for any v,w∈𝒰, where 𝒰 is a vector space with a non-degenerate inner product.

Definition: Invariant (alternately referred to as Stable)
Let 𝒲 be a subspace of 𝒰, and let S: 𝒰->𝒰 be a symmetric linear map.  𝒲 is invariant under S if S(𝒲) ⊂ 𝒲, that is for all u ∈𝒲 we have S u∈𝒲

Definition: Perpendicular Subspaces.
Let 𝒰 be a vector space with a inner product.  Let 𝒲 be a subspace of 𝒰.  Then the subspace of 𝒰 perpendicular to 𝒲, denoted (𝒲^⊥) (read "𝒲 perp"), is the set of all u∈𝒰 such that ⟨u,w⟩=0 for all w∈𝒲.

Definition: Column Space
The column space of an m x n matrix A, written as Col(A), is the set of all linear combinations of the columns of A.  If A=[(Subscript[a, 1]),...,(Subscript[a, n])], then Col(A)=Span{(Subscript[a, 1]),...,(Subscript[a, n])}.  Alternately, the column space may be thought of as Col(A)={b|b=A x for some x in (ℝ^n)}.  This alternate definition of the column space clearly shows that it is equal to the range of A and a subspace of (ℝ^m).

Definition:  Null Space
The null space of a matrix A is the set Null(A) of all solutions to the homogeneous equation A x=(Subscript[0, 𝒱]).

Definition: Rank
The rank of a matrix A, denoted by rank(A), is the dimension of the column space of A.  Alternately, the rank may be thought of as the number of linearly independent columns in A.

Definition: Matrix Norm
A matrix norm is a function |||•|||:(ℝ^(n x n))↦ℝsuch that when A,B∈(ℝ^(n x n)), then the following four properties are satisfied.
a.)  |||A|||>0  if A≠0
b.)  |||α A|||=|α| |||A|||
c.)  |||A+B|||≤|||A|||+|||B|||
d.)  |||A B|||≤|||A||| |||B|||(consistency)
One such matrix norm is |||•||(Subscript[|, 2]):(ℝ^(n x n))↦ℝ This is referred to as the two-norm.  ||| A |||2 = max((Sqrt[(Subscript[λ, i])])) where (Sqrt[(Subscript[λ, i])])is the (i^th) eigenvalue of (A^T)A.  Sometimes,  the two-norm is defined to be the square root of the spectral radius of (A^T)A.

Appendix (2.)  Mathematica notebook which generates examples of the spectral theorem

M=(({{1, 2}, {3, -1}, {-2, -3}}));

<<LinearAlgebra`Orthogonalization`
<<LinearAlgebra`MatrixManipulation`
Off[General::spell1];

(*M=FullSimplify[Transpose[M].M]*)
A=Transpose[Eigenvectors[M]];
diag=Eigenvalues[M];
sortdiag=Sort[diag,Greater];
Λ=DiagonalMatrix[sortdiag];

n=Dimensions[A,2][[2]];
A=Transpose[A];
temp={};
For[i=1,i≤n,i++,
  AppendTo[temp,Normalize[A[[Position[diag,sortdiag[[i]]][[1]][[1]],All]]]]]
A=Transpose[temp];

A//MatrixForm
Λ//MatrixForm

Simplify[A.Λ.Transpose[A]]//MatrixForm

Appendix (3.)  Mathematica notebook which calculates the singular value decomposition of a matrix as well as generate an animation of the geometric interpretation for two by two case

There is a small error in the code below which causes one of the rotations to go the wrong direction.  I will correct it when I find it.  Other than that, this code is pretty decent.  It is the result of four days of solid work, and so I include it even with its imperfections.

<<LinearAlgebra`Orthogonalization`
<<LinearAlgebra`MatrixManipulation`
Off[General::spell1];
Off[ParametricPlot::pptr];

A=(({{1, 3}, {-5, 7}}));

m=Dimensions[A][[1]];
n=Dimensions[A][[2]];

V=Eigenvectors[Transpose[A].A];
diag=Eigenvalues[Transpose[A].A]//N;
sortdiag=Sort[diag,Greater]//N;
n=Dimensions[V,2][[2]];
temp={};
For[i=1,i≤n,i++,
AppendTo[temp,Normalize[V[[Position[diag,sortdiag[[i]]][[1]][[1]],All]]]]]
V=Transpose[temp]//N;

U={};
rank=Length[diag]-Count[diag,0.];

For[i=1,i≤rank,i++,
AppendTo[U,(1/(Sqrt[sortdiag[[i]]]))A.V[[All,i]]]];
For[i=rank+1,i≤m,i++,
AppendTo[U,Flatten[ZeroMatrix[1,m]]]];
U=Transpose[U];

If[m==rank && n==rank,
     (*Print[1];*)
    Λ=DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],
    If[m==rank,
      (* Print[2];*)
       Λ=AppendRows[DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[rank,n-rank]],
       If[n==rank,
          (*Print[3];*)
          Λ=AppendColumns[DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[m-rank,rank]],
         (*Print[4];*)
          Λ=BlockMatrix[{{DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[rank,n-rank]}, 
 {ZeroMatrix[m-rank,rank],ZeroMatrix[m-rank,n-rank]}}]]]];

Print["A = ",A//MatrixForm]
Print["U = ",N[U]//MatrixForm]
Print["Λ = ",N[Λ]//MatrixForm]
Print["V = ",N[V]//MatrixForm]
Print["U·Λ·\!\(V\^T\) = ",N[U.Λ.Transpose[V]]//MatrixForm]

x01[t_]=t;
x02[t_]=t;
x03[t_]=Abs[t];
x04[t_]=0;
y01[t_]=(Sqrt[1-(t^2)]);
y02[t_]=-(Sqrt[1-(t^2)]);
y03[t_]=0;
y04[t_]=Abs[t];

ν=ArcTan[((Transpose[V][[2,1]])/(Transpose[V][[1,1]]))];
x11[t_,θ_]=x01[t] Cos[θ]-y01[t] Sin[θ];
x12[t_,θ_]=x02[t] Cos[θ]-y02[t] Sin[θ];
x13[t_,θ_]=x03[t] Cos[θ]-y03[t] Sin[θ];
x14[t_,θ_]=x04[t] Cos[θ]-y04[t] Sin[θ];
y11[t_,θ_]=x01[t] Sin[θ]+y01[t] Cos[θ];
y12[t_,θ_]=x02[t] Sin[θ]+y02[t] Cos[θ];
y13[t_,θ_]=x03[t] Sin[θ]+y03[t] Cos[θ];
y14[t_,θ_]=x04[t] Sin[θ]+y04[t] Cos[θ];
a=Λ[[1,1]];b=Λ[[2,2]];
x21[t_,α_,β_]=t;
x22[t_,α_,β_]=t;
x23[t_,α_,β_]=Abs[t];
x24[t_,α_,β_]=0;
y21[t_,α_,β_]=(β/α)(Sqrt[(α^2)-(t^2)]);
y22[t_,α_,β_]=-(β/α)(Sqrt[(α^2)-(t^2)]);
y23[t_,α_,β_]=0;
y24[t_,α_,β_]=(β/α) Abs[t];
x31[t_,θ_,α_,β_]=x21[t,α,β] Cos[θ]-y21[t,α,β] Sin[θ];
x32[t_,θ_,α_,β_]=x22[t,α,β] Cos[θ]-y22[t,α,β] Sin[θ];
x33[t_,θ_,α_,β_]=x23[t,α,β] Cos[θ]-y23[t,α,β] Sin[θ];
x34[t_,θ_,α_,β_]=x24[t,α,β] Cos[θ]-y24[t,α,β] Sin[θ];
y31[t_,θ_,α_,β_]=x21[t,α,β] Sin[θ]+y21[t,α,β] Cos[θ];
y32[t_,θ_,α_,β_]=x22[t,α,β] Sin[θ]+y22[t,α,β] Cos[θ];
y33[t_,θ_,α_,β_]=x23[t,α,β] Sin[θ]+y23[t,α,β] Cos[θ];
y34[t_,θ_,α_,β_]=x24[t,α,β] Sin[θ]+y24[t,α,β] Cos[θ];
υ=ArcTan[((U[[2,1]])/(U[[1,1]]))];
x41[t_,θ_,ω_,α_,β_]=x31[t,θ,α,β] Cos[ω]-y31[t,θ,α,β] Sin[ω];
x42[t_,θ_,ω_,α_,β_]=x32[t,θ,α,β] Cos[ω]-y32[t,θ,α,β] Sin[ω];
x43[t_,θ_,ω_,α_,β_]=x33[t,θ,α,β] Cos[ω]-y33[t,θ,α,β] Sin[ω];
x44[t_,θ_,ω_,α_,β_]=x34[t,θ,α,β] Cos[ω]-y34[t,θ,α,β] Sin[ω];
y41[t_,θ_,ω_,α_,β_]=x31[t,θ,α,β] Sin[ω]+y31[t,θ,α,β] Cos[ω];
y42[t_,θ_,ω_,α_,β_]=x32[t,θ,α,β] Sin[ω]+y32[t,θ,α,β] Cos[ω];
y43[t_,θ_,ω_,α_,β_]=x33[t,θ,α,β] Sin[ω]+y33[t,θ,α,β] Cos[ω];
y44[t_,θ_,ω_,α_,β_]=x34[t,θ,α,β] Sin[ω]+y34[t,θ,α,β] Cos[ω];
n=25;

RotateV={};
StretchΛ={};
RotateU={};

unitcircle=ParametricPlot[{{x11[t,0],y11[t,0]},{x12[t,0],y12[t,0]},{x13[t,0],y13[t,0]},{x14[t,0],y14[t,0]}},{t,-1,1},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"The Unit Circle"];

For[i=0,i≤n,i++,
ϕ=(ν/n)i;
ψ=(υ/n)i;
fx=1+((a-1)/n)i;
fy=1+((b-1)/n)i;

AppendTo[RotateV,ParametricPlot[{{x11[t,ϕ],y11[t,ϕ]},{x12[t,ϕ],y12[t,ϕ]},{x13[t,ϕ],y13[t,ϕ]},{x14[t,ϕ],y14[t,ϕ]}},{t,-1,1},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Rotation : \!\(\* StyleBox[\"V\",\nFontSlant->\"Italic\"]\^\(\(\(\\\\\)\(\\ \)\) T\)\)"]];

AppendTo[StretchΛ,ParametricPlot[{{x31[t,ν,fx,fy],y31[t,ν,fx,fy]},{x32[t,ν,fx,fy],y32[t,ν,fx,fy]},{x33[t,ν,fx,fy],y33[t,ν,fx,fy]},{x34[t,ν,fx,fy],y34[t,ν,fx,fy]}},{t,-fx,fx},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Componentwise  Stretch : Λ"]];

AppendTo[RotateU,ParametricPlot[{{x41[t,ν,ψ,a,b],y41[t,ν,ψ,a,b]},{x42[t,ν,ψ,a,b],y42[t,ν,ψ,a,b]},{x43[t,ν,ψ,a,b],y43[t,ν,ψ,a,b]},{x44[t,ν,ψ,a,b],y44[t,ν,ψ,a,b]}},{t,-a,a},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Rotation : \!\(\* StyleBox[\"U\",\nFontSlant->\"Italic\"]\)"]]]

SVDMovie=Flatten[{RotateV,StretchΛ,RotateU}];

For[i=1,i≤Length[SVDMovie],i++,
Show[SVDMovie[[i]],DisplayFunction->$DisplayFunction]];

SVDecomp={{unitcircle,RotateV[[n+1]]},{StretchΛ[[n+1]],RotateU[[n+1]]}};

Show[GraphicsArray[SVDecomp]];

Appendix (4.)  Matlab m-file which performs singular value decomposition for image compression

function [IMAGE] = picttomat(string)
%
% [IMAGE] = picttomat(string)
%
% Reads in jpegs and saves them as matrices. Remember
% that you must be in the same directory as the jpeg
% that you are trying to input.  You also must enter
% the titles as strings.  For example, example.jpg 
% becomes 'example.jpg'
%
% Also see..............................
% contentsvd.m
% -> pictomat.m
% rgbdecomp.m
% im.m
% dashow.m

IMAGE = imread(string,'jpg');

function [U,S,V] = rgbdecomp(A)

% [M] = im(U,S,V)[U,S,V] = rgbdecomp(A)
%
% Finds U*S*V', the SVD of and image A and returns the smallest 
% matrices U,S, and V required to recreate A.  
%
% Note: View the image approximations using the m-file im.m
%
% Also see..............................
% contentsvd.m
% pictomat.m
% -> rgbdecomp.m
% im.m
% dashow.m
tic;
A = double(A);
[r,c,h]=size(A);

U = zeros(r,r,h);
V = zeros(c,c,h);
TEMP = zeros(r,c,h);

[U(:,:,1),TEMP(:,:,1),V(:,:,1)]=svd(A(:,:,1));
[U(:,:,2),TEMP(:,:,2),V(:,:,2)]=svd(A(:,:,2));
[U(:,:,3),TEMP(:,:,3),V(:,:,3)]=svd(A(:,:,3));

S=zeros(min(r,c),3);
S(:,1)=diag(TEMP(:,:,1));
S(:,2)=diag(TEMP(:,:,2));
S(:,3)=diag(TEMP(:,:,3));

R = min(r,c);
U = U(:,1:R,:);
V = V(:,1:R,:);

toc;

function [M1] = im(U,S,V)

% [M1] = im(U,S,V)
% Creates a 50 frame movie from U*S*V'.  Here U,S, and V are 
% the SVD of some  image A.  M1 is the resulting movie.
% 
% Note: Create U,S, and V using the m-file rgbdecomp.m
%
% Also see..............................
% contentsvd.m
% pictomat.m
% rgbdecomp.m
% -> im.m
% dashow.m

m=size(U,1);
n=size(V,1);
h=size(S,2);

clf reset;
clc;
map = colormap;

C = zeros(m,n,3);

for num=1:50,
    
    C(:,:,1) = C(:,:,1) + U(:,num,1)*S(num,1)*V(:,num,1)'; %Red
    C(:,:,2) = C(:,:,2) + U(:,num,2)*S(num,2)*V(:,num,2)'; %Green
    C(:,:,3) = C(:,:,3) + U(:,num,3)*S(num,3)*V(:,num,3)'; %Blue
    p = 100*((m*num)+num+(n*num))/(m*n);
    title(p);
    axis off;
    M1(num) = im2frame(uint8(C));
   
end

close(1);

function [ ] = dashow(M);
%
% [ ] = dashow(M)
%
% The purpose of this program is to display a premade "movie"
% of slides.  The key to this program is that it allows certain
% frames of the movie to be skipped - that is, the user has 
% control over how much is shown WHILE THE MOVIE IS IN PROGRESS.
%
% The movie M is made up of n sequences (columns) where each sequence
% has m frames (rows).  dashow.m plays one sequence after another.
% It moves to the next sequence when one of two things happen.  
% 1.) All of the frames in the sequence are displayed.
% 2.) The figure (window) in which the sequence is being displayed
% is closed down.  
% In either case, dashow.m progresses to the next sequence of frames.
%
% Movies M may be created using the program im.m.  Im.m creates
% sequences of frames M1,M2,M3,... and displays the movie created by
% M=[m1' m2' m3' ...]
%
% Also see..............................
% contentsvd.m
% pictomat.m
% rgbdecomp.m
% im.m
% -> dashow.m

figure(1);
title('Background: Maximize for best effect');
axis off;
pause(5);
[m,n] = size(M);
tic;


for j=1:n,
    [r,c,h] = size(frame2im(M(1,j)));
    figure(j+1);
    for i=1:m,
        if gcf==j+1
            image(frame2im(M(i,j)));
            title(100*i*(h*(r+c+1))/(r*c*h));
            axis off;
            pause(1.25);
        else
            
            break;
        end
    end
end
close(1);

toc;

Appendix (5.)  Project Proposal

Dusty Wilson
Singular Value Decomposition
Masters Project Proposal
Advisor: Dr. John Woll

In Dr. Branko Curgus's "Abstract Linear Algebra" course of 1999, we concluded our study with the Real Valued Spectral Theorem.  Since this theorem was proved on the last day of class, I propose to revisit and prove it once again.  Then I propose to prove the Singular Value Decomposition Theorem which generalizes the Real Valued Spectral Theorem from symmetric transformations to any real valued linear transformation.  In doing so, I will give the geometric explanation of singular value decomposition.  Furthermore, I will demonstrate my  understanding of applications of singular value decomposition including matrix approximation, rank approximation, and image compression through quantitative examples.  The theorems and applications are as follows.

Theorem: Real Valued Spectral Theorem
Let 𝒰 be an n dimensional vector space over ℝ with a positive definite inner product.  Let S:𝒰 ↦ 𝒰 be a linear transformation.  Then S is  symmetric with respect to the inner product if and only if 𝒰 has an orthonormal basis B consisting of the eigenvectors of the linear transformation S.

Theorem: Singular Value Decomposition
Any m x n real valued matrix A can be decomposed into the product

A=U·Λ·(V^T)

where (Subscript[U, m x m]) and (Subscript[V, n x n]) are orthonormal matrices and (Subscript[Λ, m x n]) is a "diagonal" matrix whose entries are the square roots of the eigenvalues of the matrix (A^T)A.

Application: Matrix Approximation
Let (Subscript[A, m x n]) be a real valued matrix with the singular value decomposition A=U·Λ·(V^T)and the singular values (Sqrt[(Subscript[λ, 1])])≥... ≥(Sqrt[(Subscript[λ,  r])])>0 where r is the rank of A.  The matrix A may be approximated by disregarding the smallest r-k singular values.
 Define (Subscript[D, k x k])=
Now, we have that  M=U·Λ·(V^T)≅(({{(Subscript[U, k]), (Subscript[U, m-k])}}))·(({{(Subscript[D, k x k]), (Subscript[0, k x (n-k)])}, {(Subscript[0, (m-k) x k]), (Subscript[0, (m-k) x (n-k)])}}))·(({{(Subscript[V, k]^T)}, {(Subscript[V, n-k]^T)}}))=(Subscript[U, k])·(Subscript[D, k x k])·(Subscript[V, k]^T)=(Subscript[M, k])

Application: Rank Approximation
We find the rank of A by counting the number of non-zero singular values of A.  However, the space of matrices with full rank is open and dense and so the computer representation of A most likely is of full rank.  In order to determine a low rank approximation of A, we choose ϵ>0 such that for a relatively small k, (Sqrt[(Subscript[λ, 1])])≥...≥(Sqrt[(Subscript[λ, k])])≫ϵ>(Sqrt[(Subscript[λ, k+1])])≥...  Thus, we approximate the computer representation of A with the reduced matrix (Subscript[A, k]).

Application: Image Compression
An image may be represented as (Subscript[A, m x n]), a real valued matrix whose entries correspond to the color at each point on the image.  Using the matrix approximation and rank approximation, we find  the smallest k such that (Subscript[A, k]) acceptably approximates the image A.

References Include:
Curgus, Branko.  "Notes from Math 504."  1999
Lay, David C.  "Linear Algebra and its Applications."  (New York, Addison Wesley) 1997
Watkins, David.  "Fundamentals of Matrix Computation."  (New York, Wiley) 1991

Bibliography

Axler, Sheldon. "Linear Algebra Done Right (2^nd) ed."  (New York, Springer) 1997
Curgus, Branko.  "Notes from Math 504"  1999
Lang, Serge.  "Linear Algebra, (3^rd) ed."  (New York, Springer-Verlag) 1987
Lay, David C.  "Linear Algebra and its Applications."  (New York, Addison Wesley) 1997
Watkins, David.  "Fundamentals of Matrix Computation."  (New York, Wiley) 1991
Woll, John. "Singular Value Decomposition."