\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}

\DeclareMathOperator*{\argmin}{arg\,min}

\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Python,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={(*@}{@*)}          % if you want to add a comment within your code
}

\begin{document}


\begin{center}
  \Huge{Principal Components}\\
  \Large{MATH1900: Machine Learning}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/ml\_2019/principal\_components/principal\_components.pdf
\hrule

\begin{center}
  \includegraphics[scale=0.90]{components.png}
  \vskip 0.1in
  {\it{What are the top two trends in my data?}}
\end{center}
\vskip 0.1in

\begin{tcolorbox}[colback=blue!5,title=Principal Component Analysis]
{\it{Find a small set of orthogonal basis vectors that approximate a large data set.}}
\end{tcolorbox}

\vskip 0.1in
If we ask 1,000 people to fill out a survey of 50 questions, it's likely
that every survey will be different.  However, it may also be the case
that distinct patterns can be observed, corresponding to differences in
men and women, young and old, conservative and liberal.  Identifying
these patterns will help us to understand the data, to replace the
raw data by a simpler model that explains much of the results.
\vskip 0.1in
We can use principal component analysis (PCA) to search for these
patterns.  We try to boil down our data to reveal the most information 
in the fewest number of components. 
\vskip 0.1in
If there was only one question on the survey, it would make sense to
compute the average answer, and then the variance to report how much
answers can deviate from it.  Now, however, our task is more complicated. 
\vskip 0.1in
We will think of our data as a big, rectangular, numerical
matrix, and we will see that the singular value decomposition (SVD) can
give us useful answers.

\section{A = U * S * V}

The SVD factorization of an $m \times n$ matrix A (as computed in Python) 
has the form $A = U S V$ where the matrix $S$ has the shape of $A$, 
but is nonzero only on the diagonal, while the $m \times m$ matrix $U$ 
and the $n \times n$ matrix $V$ are orthogonal.  
\vskip 0.1in
In Python, if $A$ is an {\it{np.array()}}, then we can request the SVD
factorization by:
\begin{lstlisting}
U, s, V = np.linalg.svd ( A );
\end{lstlisting}
Instead of returning the matrix $S$, we get back a vector $s$, of length
$mn = \min(m,n)$, containing the diagonal entries of $S$.  We can build $S$ by
\begin{lstlisting}[caption = Create matrix S from vector s]
S = np.zeros ( A.shape )
mn = min ( m, n )
for i in range ( 0, mn ):
  S[i,i]= s[i]
\end{lstlisting}
\vskip 0.1in
To verify that the factorization is correct, we can do the following test:
\begin{lstlisting}[caption = Compare A and U*S*V]
  SV = np.matmul ( S, V )
  USV = np.matmul ( U, SV )
  diff = np.linalg.norm ( A - USV )
\end{lstlisting}
\vskip 0.1in
{\bf{Exercise:}} Compute {\tt{A = np.random.randn ( 5, 6 )}}, compute
the SVD, and verify that the factors, when multiplied, are almost
exactly equal to {\tt{A}}.
\vskip 0.1in
Most texts and algorithms return the transpose of
the matrix {\tt{V}}.  Python is the only system that I know of where
this is not so!
\vskip 0.1in
You can look at the code {\it{svd\_product.py}}, which checks that
$A = U * S * V$ for several small matrices.

\section{Theory of approximating A with a partial SVD}

You know that if column vectors $u$ and $v$ both have length $n$ that
you can compute the vector dot product, or vector inner product,
\begin{displaymath}
  u {\text{ dot }} v = u' \, v =\sum_{i=1}^n u_i v_i
\end{displaymath}
\vskip 0.1in
You may be less familiar with the {\it{outer product}} of 
a column $m$-vector $u$ and a row $n$-vector $v$:
\begin{displaymath}
  u \otimes v = u \, v
\end{displaymath}
which results in an $m \times n$ matrix $A$ whose $(i,j)$ entry is:
\begin{displaymath}
  A_{i,j} = u_i \, v_j
\end{displaymath}
Thus, if we have
\begin{displaymath}
  u = 
  \left( 
    \begin{array}{r} 1\\2\\3\\ 
    \end{array} 
  \right),
  v = 
  \left( 
    \begin{array}{rr} 10 & -2\\ 
    \end{array} 
  \right),
\end{displaymath}
Then
\begin{displaymath}
  u \otimes v = u \, v = \left( 
  \begin{array}{rrrr}
   10&  -2 \\
   20&  -4 \\
   30&  -6 \\
  \end{array}
  \right)
\end{displaymath}
So, although normally a $1000 \times 50$ matrix requires 50,000 values
to describe it, we can easily create a matrix that large just by
an outer product that requires a total of 1000 + 50 values.  The 
question is, what good is such an ``inflated'' matrix?
\vskip 0.1in
The answer goes the other way.  If our 50,000 values can actually
be described by 1050 values, then we have compressed the data down to
the actual information.  Even if we have to create 3 or 4 such
outer products and add them together and we only get an approximation
to our original big matrix, such an approximation could be good enough
for our purposes, and might also reveal some of the underlying structure of
the raw data.
\vskip 0.1in
So now let's see how the SVD can be used in this way.
\vskip 0.1in
Let's suppose that $A$ is an $m \times n$ matrix.  Our survey data,
for instance, might be stored as a $1000 \times 50$ table of numbers.
Our computational puzzle is: can we approximate the large amount of data 
in $A$ by a small set of numbers, which capture most of the information?
\vskip 0.1in
We know that computing the average answer for each question gives us
a lot of information to start with, so we can start our analysis
by creating a matrix $A_0$ so that each row contains the same 50 
average values.  We can compute this as an outer product,
$A_0 = u_0 \otimes v_0$ where the $m$ vector $u_0$ and the
$n$-vector $v_0$ are:
\begin{displaymath}
  u_0 = 
  \left( 
    \begin{array}{r} 1\\ 1\\ ...\\ 1\\ 
    \end{array} 
  \right), \quad
  v_0 = 
  \left( 
    \begin{array}{rrrr} \mu_1& \mu_2& ...& \mu_n\\ 
    \end{array} 
  \right)
\end{displaymath}
\vskip 0.1in
Now let's compute the SVD of $A-A_0$, so that we have exactly:
\begin{displaymath}
  A = A_0 + U \, S \, V = u_0 \otimes {v_0} + U \, S \, V
\end{displaymath}
The SVD suggests that, if we wanted to use just one more
matrix outer product, we should use the first column of $U$,
the first row of V, and multiply it by the first entry in S.
This approximation would have the form:
\begin{displaymath}
  A \approx A_0 + A_1 = u_0 \otimes {v_0} + s_1 \, u_1 \otimes {v_1}
\end{displaymath}
Using one more column gives us a better approximation:
\begin{displaymath}
  A \approx A_0 + A_1 + A_2 = u_0 \otimes {v_0} + s_1 \, u_1 \otimes {v_1}
  + s_2 \, u_2 \otimes {v_2}
\end{displaymath}
\vskip 0.1in
We can tell how good our approximation is.  If $A_k$ is
 formed from the first $k$ columns of the SVD, then:
\begin{displaymath}
  ||A-A_0-A_1-...-A_k||_2^2 = 
  \frac
  {\sum_{i=k+1}^{mn} s_i^2 }
  {\sum_{i=1}^{mn} s_i^2 }
\end{displaymath}

\section{Exercise: approximating A with a partial SVD}

Let's test out the theory on a simple $5 \times 5$ matrix $A$.  Once we
subtract the mean, we expect to have just 4 nonzero singular values.
Let's add up the first few outer products and see how we do:
\begin{lstlisting}[caption=Define a matrix $A$]
  import numpy as np
  m = 5
  n = 5
  A = np.array ( [
    [17, 24,  1,  8, 15],
    [23,  5,  7, 14, 16],
    [ 4,  6, 13, 20, 22],
    [10, 12, 19, 21,  3],
    [11, 18, 25,  2,  9] ] )
\end{lstlisting}

\begin{lstlisting}[caption=Compute mean and create $A0$]
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
\end{lstlisting}

\begin{lstlisting}[caption=Request SVD of $A-A_0$]
  U, s, V = np.linalg.svd ( A - A0 )
  diff = np.linalg.norm ( A - A0 )
\end{lstlisting}

\begin{lstlisting}[caption=Create outer product $A_1$]
  A1 = s[0] * np.outer ( U[:,0], V[0,:] )
  diff = np.linalg.norm ( A - A0 - A1 )
\end{lstlisting}

\begin{lstlisting}[caption=Create outer product $A_2$]
  A2 = s[1] * np.outer ( U[:,1], V[1,:] )
  diff = np.linalg.norm ( A - A0 - A1 - A2 )
\end{lstlisting}

At this halfway point, we can compare A and our approximation A0+A1+A2:
\begin{verbatim}
   -----A-----            -----A0+A1+A2-----
[[17 24  1  8 15]     [[20.7 20.4  2.9  5.7 15.2]
 [23  5  7 14 16]      [16.0 10.2  4.6 15.8 18.1]
 [ 4  6 13 20 22]      [ 9.6  3.9 13.0 22.0 16.3]
 [10 12 19 21  3]      [ 7.8 10.1 21.3 15.7  9.9]
 [11 18 25  2  9]]     [10.7 20.2 23.0  5.5  5.2]]
\end{verbatim}

The SVD promises that our sequence of approximations is the best
that can be computed, and always gives us a report of how far
off we are by looking at the singular values.
\vskip 0.1in
Once we have computed $A_0$, we can jump to any approximation
in one step.  This requires setting up the full matrix $S$,
and doing a nested multiplication of $S*V$, then of $U*(S*V)$:
For instance, to compute the third approximation 
A123 = A1 + A2 + A3:
\begin{lstlisting}[caption=Create an outer product $A_{123}$ in one step]
  A123 = np.matmul ( U[:,0:3], np.matmul ( S[0:3,0:3], V[0:3,:] ) )
\end{lstlisting}
\vskip 0.1in
You can look at some of these computations in {\it{svd\_approx.py}}.

\section{PCA by Singular Value Decomposition of A}

As a small, but more realistic example, we are ready to consider
a dataset containing 214 records, describing 9 measurements of
chemical properties of samples of glass.  
\vskip 0.1in
We will read a data file {\it{glass\_data.txt}} into a $m \times n$ 
array $A$.  We compute the average values of the 9 columns and 
write $ A_0 = u_0 \otimes v_0$, where each row of $A_0$ stores
the average values.  We compute the SVD of $A-A_0$.  The singular 
values $s$ tell us the importance of each column of $U$ and row of $V$.
If we approximate $A$ with a part of the SVD, then the corresponding
singular values tell us how close our approximation is.
\begin{lstlisting}[caption = Read the data into $A$]
  import numpy as np
  A = np.loadtxt ( 'glass_data.txt' )
  m, n = A.shape
\end{lstlisting}

\begin{lstlisting}[caption=Create column average matrix $A_0$]
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
  print ( '||A-A0|| = ', np.linalg.norm ( A - A0 ) )
\end{lstlisting}

\begin{lstlisting}[caption = Get SVD of $A-A_0$]
  U, s, V = np.linalg.svd ( A - A0 )
\end{lstlisting}

\begin{lstlisting}[caption = Create $A_1$  first outer product]
  A1 = s[0] * np.outer ( U[:,0], V[0,:] ) )
  print ( '||A-A0-A1|| = ', np.linalg.norm ( A - A0 - A1 ) )
\end{lstlisting}

\begin{lstlisting}[caption = Create $A_2$ second outer product]
  A2 = s[1] * np.outer ( U[:,1], V[1,:] ) )
  print ( '||A-A0-A1-A2|| = ', np.linalg.norm ( A - A0 - A1 -A2 ) )
\end{lstlisting}

\begin{lstlisting}[caption = Plot the singular values]
  import matplotlib.pyplot as plt
  mn = min ( m, n )
  x = np.arange ( mn )
  plt.plot ( x, s )
\end{lstlisting}

\begin{lstlisting}[caption = Plot variance captured by first K vectors]
  x = np.arange ( mn + 1 )
  y = [ 0.0 ]
  y = np.append ( y, np.cumsum ( s**2 ) )
  y = y / np.sum ( s**2 )
  plt.plot ( x, y )
\end{lstlisting}

\begin{center}
  \includegraphics[scale=0.50]{svd_glass_singular_values.png}
  \includegraphics[scale=0.50]{svd_glass_percentage.png}
  {\it{Singular values and variance capture}}
\end{center}
The point of this exercise is to suggest how the SVD can be used to
handle datasets much larger than this example.  We have only looked at
the idea of approximating the data with a small set of outer products,
but with the SVD there are many other tools that we can build.
\vskip 0.1in
You can look at some of these computations in {\it{svd\_glass.py}}.

\section{PCA for Image Data}

Consider an black and white image as an $m \times n$ array $A$ of numbers.
Now imagine each column of this array as being a data item, so that we have
$n$ items, each of dimension $m$.  (You need the PIL library for this!) 
\vskip 0.1in
If we compute the SVD of A,  we can approximate the
image using combinations of a small set of singular vectors, and
compute how much information we have captured.
\begin{lstlisting}[caption = Start up and get the image data]
  import numpy as np
  import matplotlib.pyplot as plt
  from PIL import image
  I = Image.open ( 'casablanca.png' )
  I = I.convert ( 'LA' )
\end{lstlisting}
 
\begin{lstlisting}[caption = Convert image to numeric format and display]
  A = np.array ( list ( I.getdata(band=0)),float )
  A.shape = ( I.size[1], I.size[0] )
  plt.imshow ( A, cmap = 'gray' )
  plt.title ( 'Numeric version of image' )
  plt.show ( )
\end{lstlisting}

\begin{lstlisting}[caption = Compute column averages]
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
\end{lstlisting}

\begin{lstlisting}[caption = Compute SVD and build S]
  U, s, V = np.linalg.svd ( A - A0 )
  S = np.zeros ( A.shape )
  for i in range ( 0, min ( m, n ) ):
    S[i,i] = s[i]
\end{lstlisting}

\begin{lstlisting}[caption = Construct average using A0 and 5 outer products]
  r = 5
  A5 = A0 + np.matmul ( U[:,0:r], np.matmul ( S[0:r,0:r], V[0:r,:] ) )
  plt.imshow ( A5, cmap = 'gray' )
  plt.title ( 'Reconstruction using %d components' % ( r ) )
  plt.show ( )
  plt.clf ( )
\end{lstlisting}
The last block of code can be modified to use 10, 20 or 40 
outer products.
\vskip 0.1in
You can look at some of these computations in {\it{svd\_image.py}}.

\section { Computing Assignment \#3}

Write a Python program which
\begin{packed_itemize}
  \item{Sets {\tt{A}} to the data in {\it{concrete\_data.txt}} (on web page);}
  \item{Sets {\tt{A0}} to the matrix in which each row contains
  the column averages;}
  \item{Computes {\tt{U, s, V}} of the SVD of {\tt{A-A0}};}
  \item{Computes {\tt{A1}} through {\tt{A3}}, 
   the first three SVD outer products;}
  \item{Prints $||A||, ||A-A0||, ||A-A0-A1||, ||A-A0-A1-A2||, 
  ||A-A0-A1-A2-A3||$.}
\end{packed_itemize}
\vskip 0.1in
Email a copy of your program to Dr Schneier {\bf{mhs64@pitt.edu}} 
by 11:59pm Friday, 27 September.

\end{document}
