svdNotes.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{svdNotes}
%\renewcommand{\dirname}{notes/ece1254/}
%%\newcommand{\dateintitle}{}
%\newcommand{\keywords}{Singular Value Decomposition, SVD, image space, orthonormal basis}
%
%\input{../peeter_prologue_print2.tex}
%
%\beginArtNoToc
%
%\generatetitle{Singular Value Decomposition}
%\label{chap:svdNotes}
%\section{Motivation}
%
In \cref{dfn:multiphysicsL6:1} \textAndIndex{singular value decomposition} (SVD) was presented without any mention of how to compute it.
%\ref{dfn:svdNotes:1}
%
% FIXME: omit on integration
%Recall the definition \cref{dfn:svdNotes:1}
%
%%%\makedefinition{Singular value decomposition (SVD)}{dfn:svdNotes:1}{
%%%...
%
Here I'll review some of the key ideas from the MIT OCW SVD lecture by Prof.\ Gilbert Strang \citep{ocw:svd}.  This is largely to avoid having to rewatch this again in a few years as I just did.
%
The idea behind the SVD is to find an orthogonal basis that relates the image space of the transformation, as well as the basis for the vectors that the transformation applies to.  That is a relation of the form
%
\begin{equation}\label{eqn:svdNotes:281}
\Bu_i = M \Bv_j
\end{equation}
%
Where \( \Bv_j \) are orthonormal basis vectors, and \( \Bu_i \) are orthonormal basis vectors for the image space.  %Strang's lectures call these the row and column spaces (or reverse?)
%
Suppose that such a set of vectors can be computed and that \( M \) has a representation of the desired form
%
\begin{equation}\label{eqn:svdNotes:301}
M = U \Sigma V^\conj
\end{equation}
%
where
\begin{equation}\label{eqn:svdNotes:321}
U =
\begin{bmatrix}
   \Bu_1 & \Bu_2 & \cdots & \Bu_m
\end{bmatrix},
\end{equation}
%
and
%
\begin{equation}\label{eqn:svdNotes:341}
V =
\begin{bmatrix}
   \Bv_1 & \Bv_2 & \cdots & \Bv_n
\end{bmatrix}.
\end{equation}
%
By left or right multiplication of \( M \) with its (conjugate) transpose,
the decomposed form of these products can be observed to have a very simple form
%
\begin{subequations}
\begin{equation}\label{eqn:svdNotes:361}
M^\conj M
=
V \Sigma^\conj U^\conj U \Sigma V^\conj
=
V \Sigma^\conj \Sigma V^\conj,
\end{equation}
\begin{equation}\label{eqn:svdNotes:381}
M M^\conj
=
U \Sigma V^\conj V \Sigma^\conj U^\conj
=
U \Sigma \Sigma^\conj U^\conj.
\end{equation}
\end{subequations}
%
Because \( \Sigma \) is diagonal, the products \( \Sigma^\conj \Sigma \) and \( \Sigma \Sigma^\conj \) are also both diagonal, and populated with the absolute squares of the singular values that have been presumed to exist.  Because \( \Sigma \Sigma^\conj \) is an \( m \times m \) matrix, whereas \( \Sigma^\conj \Sigma \) is an \( n \times n \) matrix, so the numbers of zeros in each of these will differ, but each will have the structure
%
\begin{equation}\label{eqn:svdNotes:401}
\Sigma^\conj \Sigma \sim \Sigma \Sigma^\conj \sim
\begin{bmatrix}
   \Abs{\sigma_1}^2 &     &   &     &    &    &\\
      & \Abs{ \sigma_2 }^2 &   &     &    &    &\\
      &     & \ddots &     &    &    &\\
      &     &   & \Abs{ \sigma_r }^2 &    &    &\\
      &     &     &   &       0 &    &    \\
      &     &     &   &    & \ddots  &    \\
      &     &     &   &    &    & 0       \\
\end{bmatrix}
\end{equation}
%
This shows one method of computing the singular value decomposition (for full rank systems).  Solution of the eigensystem of either \( M M^\conj \) or \( M^\conj M \) is required to find both the singular values, and one of \( U \) or \( V \).  \( U \) can be found from the \( r \) orthonormal eigenvectors of \( M M^\conj \) supplemented with a mutually orthonormal set of vectors from the Null space of \( M \).  \( \Sigma \) can be found by taking the square roots of the eigenvalues of \( M M^\conj \).  With both \( \Sigma \) and \( U \) computed, \( V \) can be found by inversion.  Alternatively, it is possible to solve for \( \Sigma \) and \( V \) by computing the eigensystem of \( M^\conj M \), also supplementing those eigenvectors with vectors from the null space, and then compute \( U \) by inversion.
%
Working the examples from the lecture is instructive to gain some insight about how this works.
%
\makeexample{Full rank \( 2 \times 2 \) matrix.}{example:svdNotes:1}{
In the lecture the SVD decomposition is computed for
%
\begin{equation}\label{eqn:svdNotes:421}
M =
\begin{bmatrix}
   4 & 4 \\
   3 & -3
\end{bmatrix}
\end{equation}
%
The matrix product with its conjugate is
%
\begin{subequations}
\begin{equation}\label{eqn:svdNotes:441}
M M^\conj
=
\begin{bmatrix}
32 &    0 \\
0 &  18
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:svdNotes:461}
M^\conj M
=
\begin{bmatrix}
25 & 7 \\
7 & 25
\end{bmatrix}
\end{equation}
\end{subequations}
%
The first is already diagonalized so \( U = I \), and the singular values are found by inspection \( \{ \sqrt{32}, \sqrt{18} \} \), or
%
\begin{equation}\label{eqn:svdNotes:541}
\Sigma =
\begin{bmatrix}
\sqrt{32} & 0 \\
0 & \sqrt{18}
\end{bmatrix}
\end{equation}
%
Because the system is full rank, \( V \) follows by inversion
%
\begin{equation}\label{eqn:svdNotes:481}
\Sigma^{-1} U^\conj M  = V^\conj,
\end{equation}
%
or
\begin{equation}\label{eqn:svdNotes:501}
V = M^\conj U \lr{ \Sigma^{-1} }^\conj.
\end{equation}
%
In this case that is
%
\begin{equation}\label{eqn:svdNotes:521}
V = \inv{\sqrt{2}}
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
\end{equation}
%
\( V \) could also be computed directly by diagonalizing \( M^\conj M \).  The eigenvectors are \( \lr{ 1, \pm 1 }/\sqrt{2} \), with respective eigenvalues \( \{ 32, 18 \} \).
%
This gives \cref{eqn:svdNotes:541} and \cref{eqn:svdNotes:521}.  Again, because the system is full rank, \( U \) can be computed by inversion
%
\begin{equation}\label{eqn:svdNotes:561}
U = M V \Sigma^{-1}.
\end{equation}
%
Carrying out this calculation recovers \( U = I \) as expected.  Looks like I used a different matrix than Prof.\ Strang used in his lecture (alternate signs on the 3's).  He had some trouble that arrived from independent calculation of the respective eigenspaces.  Calculating one from the other avoids that trouble since there are different signed eigenvalues that can be chosen.  Because specific mappings between the eigenspaces that satisfy the \( \Bu_i = M \Bv_j \) constraints encoded by the relationship \( M = U \Sigma V^\conj \) are desired, these \( U \) and \( V \) matrices cannot be selected independently.
}
%
A non-full rank example, as in the lecture, is also useful
%
\makeexample{A \( 2 \times 2 \) matrix without full rank.}{example:svdNotes:2}{
%
How about
%
\begin{equation}\label{eqn:svdNotes:581}
M =
\begin{bmatrix}
1 & 1 \\
2 & 2
\end{bmatrix}.
\end{equation}
%
The matrix and conjugate product is
%
\begin{subequations}
\begin{equation}\label{eqn:svdNotes:601}
M^\conj M =
\begin{bmatrix}
5 & 5 \\
5 & 5
\end{bmatrix}
\end{equation}
\begin{equation}\label{eqn:svdNotes:621}
M M^\conj =
\begin{bmatrix}
2 & 4 \\
4 & 8
\end{bmatrix}
\end{equation}
\end{subequations}
%
For which the non-zero eigenvalue is \( 10 \) and the corresponding eigenvalue is
\begin{equation}\label{eqn:svdNotes:641}
\Bv =
\inv{\sqrt{2}}
\begin{bmatrix}
1 \\
1
\end{bmatrix}.
\end{equation}
%
This gives
%
\begin{equation}\label{eqn:svdNotes:661}
\Sigma =
\begin{bmatrix}
\sqrt{10} & 0 \\
0 & 0
\end{bmatrix}
\end{equation}
%
Since \( V \) must be orthonormal, there is only one choice (up to a sign) for the vector from the null space.
%
Try
%
\begin{equation}\label{eqn:svdNotes:681}
V =
\inv{\sqrt{2}}
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
\end{equation}
%
\( M M^\conj \) has eigenvalues \( \{ 10, 0 \} \) as expected.  The eigenvector for the non-zero eigenvalue is found to be
%
\begin{equation}\label{eqn:svdNotes:701}
\Bu = \inv{\sqrt{5}}
\begin{bmatrix}
1 \\
2
\end{bmatrix}.
\end{equation}
%
It is easy to expand this to an orthonormal basis.  Is it required to pick a specific sign relative to the choice made for \( V \)?
%
Try
%
\begin{equation}\label{eqn:svdNotes:721}
U = \inv{\sqrt{5}}
\begin{bmatrix}
1 & 2 \\
2 & -1
\end{bmatrix}.
\end{equation}
%
Multiplying out \( U \Sigma V^\conj \) gives
%
\begin{equation}\label{eqn:svdNotes:741}
\begin{aligned}
\inv{\sqrt{5}}
\begin{bmatrix}
1 & 2 \\
2 & -1
\end{bmatrix}
\begin{bmatrix}
\sqrt{10} & 0 \\
0 & 0
\end{bmatrix}
\inv{\sqrt{2}}
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
&=
\begin{bmatrix}
1 & 2 \\
2 & -1
\end{bmatrix}
\begin{bmatrix}
1 & 0 \\
0 & 0
\end{bmatrix}
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
\\ &=
\begin{bmatrix}
1 & 0 \\
2 & 0
\end{bmatrix}
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
\\ &=
\begin{bmatrix}
1 & 1 \\
2 & 2
\end{bmatrix}.
\end{aligned}
\end{equation}
%
It appears that this works.  This has not demonstrated why that should be, and could be due to luck with signs.  To fully understand how to do this type of computation, more theoretical work would be required.  That is probably a good reason to leave such a task to computational software.
}
%
%\EndArticle