\magnification=\magstep1
\centerline{Mathematics 372 -- Numerical Linear Algebra}
\centerline{Lab2/Problem Set 3 -- {\it due:} Friday, February 9}
\bigskip
\noindent
{\it General Information}
\bigskip
In this lab/problem set we will study the {\it Cholesky} decomposition
of positive-definite symmetric matrices, using MATLAB.  A particular
item of interest in one of the later questions will be the relation of 
the ``sparsity pattern'' of $A$ (that is, the locations of zero and 
nonzero entries in the matrix) and the sparsity pattern of the Cholesky 
factor $R$.
\bigskip
\noindent
{\it MATLAB Practice}
\bigskip
The command for computing Cholesky factorizations in MATLAB
is {\tt chol} (logically enough!).  To use it, you can enter
the matrix $A$ first as we discussed before, then enter a command
of the form 
\bigskip
\centerline{\tt chol(A)}
\bigskip
\noindent
The output will be the Cholesky factor $R$.  For example enter the 
matrix
$$A = \pmatrix{16 & 4 & 8 & 4\cr 4 & 10 & 8 & 4\cr 8 & 8 & 12 & 10\cr 4 & 4 & 10 & 12\cr},\leqno(1)$$
then compute the Cholesky factor and call the result $R$.  Of course we should have
$A = R^t R$ (the Cholesky factor is the ``matrix square root'').  To check that the 
result is correct, you can compute $R^t R$ in MATLAB with 
\bigskip
\centerline{\tt R' * R}
\bigskip
\noindent
(the ``prime'' is transpose in MATLAB, $*$ is matrix multiplication as we saw before).
Try it and compare with the original matrix $A$.  The matrix (1) is simpler than many 
because the Cholesky factor is also an integer matrix.  That is not always true since
as we know the computation of the entries in $R$ involves taking square roots.  
To see a more typical example, let's add a $4\time 4$ identity matrix $I$ to this $A$:
\bigskip
\centerline{\tt AA = A + eye(4)}
\bigskip
\noindent
(someone had a sense of humor here!).  Then compute the Cholesky factor and check
that it is correct.  Call your computed Cholesky factor $RR$.

The Cholesky factorization is only defined for positive definite symmetric matrices.
If you enter a matrix that is symmetric but not positive definite, 
then MATLAB will produce an error message.  For instance try {\tt chol(B)} for 
$$B = \pmatrix{-1 & 1 \cr 1 & -1}.$$
$B$ is not positive definite (why not?)  

Something more subtle happens if you enter a matrix that is not symmetric.  
In that case, MATLAB doesn't actually check to see whether the matrix is symmetric, 
it just {\it assumes that it is} and works with the entries above the diagonal 
as in the formulas for Cholesky factorization that we wrote down in class.  For 
instance, try 
$$T = \pmatrix{12 & 1 \cr 2 & 11}$$
and compute {\tt chol(T)}.  If the output is $R$, what matrix is $R'*R$ here?

Recall that we got into the whole ``matrix factorization game'' and the idea of the 
Cholesky factorization by way of the observation that solving systems of linear equations
for upper- or lower-triangular coefficient matrices is {\it especially simple}
(no row-reductions are necessary if we use backward or forward substitution).  
For instance suppose we want to solve $A x = b$ for some $b$.  Since
$A = R^t R$, as we know, this can be done in two steps.  By associativity
of matrix multiplication, $R^t(Rx) = b$ First we solve the system
$R^t y = b$, then we solve $Rx = y$ for $x$.  Carry out this process in MATLAB
for the system $AA x = b$ with $AA$ as above, and a vector $b$ of your choice.  Check your
final answer.  

Finally, for some of the problems below, we will need to introduce the idea of 
{\it sparse} matrices.  A somewhat imprecise, but suggestive, way to understand
this idea is that a sparse matrix is one where some ``large'' proportion of the total
entries are {\it zeroes}.  One could also quantify this, of course, 
by computing the actual proportion of zero entries in the matrix.  
Matrices with ``almost all'' entries nonzero are called ``dense'' matrices.  
MATLAB has an interesting command called {\tt spy} that produces
a {\it plot} indicating the sparsity structure of a matrix.  For instance, enter
\bigskip
\centerline{\tt A = [1 0 -1; 2 3 0; -1 -1 2]}
\centerline{\tt spy(A)}
\bigskip
\noindent
What does the plot mean?  Of course this is more interesting for ``big''
matrices (where we might not even want to try to list all the entries!)
For example, there are a number of ``builtin'' big example matrices in 
MATLAB, including one $479 \times 479$ matrix called {\tt west0479}.  
To access it and generate
the spy plot, enter
\bigskip
\centerline{\tt load west0479}
\centerline{\tt spy(west0479)}
\bigskip
\noindent
Even though we don't want to try to see the whole $479 \times 479$ matrix, we
can visualize where it has nonzero entries and see some things
about it.  For instance, is {\tt west0479} a symmetric matrix?
The number $nnz = 1887$ at the bottom of the plot is the actual
number of {\it nonzero} entries.  What fraction of the total entries
in this matrix are zero?

If a matrix is very sparse in this intuitive sense, then 
storing the whole matrix (with all the zeroes explicitly filled in) is an
inefficient use of the computer memory (especially if the matrix is {\it big},
say $10000\times 10000$).  So an alternate way to think about storing the information
in a matrix is to give the total numbers of rows and columns, then 
{\it list only the nonzero entries}, together with their row 
and column numbers.  If a pair of a row number and a column number 
{\it doesn't} appear, then that entry is automatically {\it zero}.  
A matrix stored in this way is properly called a matrix with a {\it sparse 
representation}, but we will also often refer to it as a {\it sparse matrix} 
(even if the proportion of zero entries is not that large).  
A moment's thought will convince you that it is possible to compute 
matrix sums, products, etc.  using this sort of representation of matrices 
(though it would be awkward by hand!).   And MATLAB has two different ways of 
dealing with matrices:
\bigskip
\item{$\bullet$} as {\it dense} matrices (all the matrices we have entered
manually to this point have been dense), or
\item{$\bullet$} as {\it sparse} (i.e. sparsely represented) matrices
({\tt west0479} above is a sparse matrix).  For instance try entering the
command
\bigskip
\centerline{\tt west0479}
\bigskip
\noindent
Two commands {\tt sparse} and {\tt dense} convert from one representation
to the other (and are inverse functions of one another!).  
With $A$ the $3\times 3$ matrix above 
\bigskip
\centerline{\tt AS = sparse(A)}
\bigskip
\noindent
The output is the sparse representation.  Note that there are exactly seven items
in the list (each item being a pair of a row and column number, together with 
the corresponding entry in $A$).  The command {\tt dense(AS)} returns
to the dense representation.

The MATLAB matrix commands we have seen to this point are set up so that 
either dense or sparse matrices can be supplied as the input ({\it as long as
the same form is used for all inputs!}), and the 
output will be a matrix of the same type.   Some other commands in the future
will {\it only} take sparse matrix inputs.

\bigskip
\noindent
{\it Problems}
\bigskip
\noindent
{\it Note: Several problems below deal with rather large matrices.
You probably won't want to print them out, either in dense or in
sparse representations!  Recall that a semicolon at the end of 
a MATLAB statement does the calculation but suppresses the output.}
\bigskip
\noindent
I.  Recall the systems of linear equations we obtained by discretizing
the 2nd order ODE boundary value problem in Lab 1.  If in the ODE 
$$ay'' + by' + cy = f,$$
the constants are 
$a = -1$, $b = 0$, $c = 1$ (as before), 
then the coefficient matrix $A_N$ is positive definite for all $N$.  
\medskip
\item{A)} Use MATLAB to compute the Cholesky factors $R_N$ for the coefficient
matrices of these systems for $N = 6, 8, 20$.   Show how you could recover
the approximate solutions of the boundary value problem
you had in the first lab using the Cholesky factors $R_N$.
\item{B)} What do you notice about the locations of nonzero entries 
in $R_N$ in these examples?  
\item{C)} (non-lab problem) An $n\times n$ matrix $A = (a_{ij})$ is said to be
{\it tridiagonal} if $a_{ij} = 0$ for all entries with $j > i + 1$ and also $i > j + 1$
(fancier way to say the same thing:  $a_{ij} = 0$ whenever $|i - j| > 1$).  The 
entries in the ``band'' consisting of the main diagonal and the upper and lower
diagonals on either side of it may or may not be zero.
Give a full mathematical proof using the formulas we developed for computing the 
Cholesky factor that if $A$ is tridiagonal, then $R = (r_{ij})$ is an upper triangular 
matrix with $r_{ij} = 0$ whenver $j > i+1$.  In other words $R$ is upper triangular
{\it and also} tridiagonal.   Suggestion:  To see the ideas you need
for this proof, it will probably help to work out one or two $4\times 4$
or $5\times 5$ examples by hand just to see which entries are actually 
used where in the computation.  The matrix
$$A = \pmatrix{4 & 1 & 0 & 0\cr
               1 & 4 & 1 & 0\cr
               0 & 1 & 4 & 1\cr
               0 & 0 & 1 & 4\cr}$$
should be instructive.  Why do the zero entries in $R$ end up being zero?

\bigskip
\noindent
II.  The ``bucky ball'' is a polyhedron with 60 vertices, 32 faces (12 pentagons 
and 20 hexagons), and 30 edges.  If you imagine taking one and ``inflating'' it
so the surface becomes spherical, you will see the edges forming the pattern of 
seams on a {\it soccer ball}.  This shape also appears in chemistry.  There 
is a carbon-hydrogen chemical compound called ``fullerene'' where carbon atoms 
are arranged at the vertices of one of these polyhedra, and the carbon-carbon 
bonds are the edges.  To see a picture, look up the Wikipedia entry for ``Fullerene''.
The names ``bucky ball'' and ``Fullerene'' 
refer to the architect and futurist Buckminster Fuller who proposed 
constructing buildings called geodesic domes with these shapes.

We can get a {\it matrix} from this shape in the following way.  First we label
the vertices $1 \ldots 60$ in some fixed way.  Then we make a $60\times 60$
matrix $B = (b_{ij})$ where 
$$
b_{ij} = \cases{ 1 & if there is an edge vertex $i$ to vertex $j$\cr
                 0 & if there is no edge\cr}$$
(There are no edges from vertex $i$ to vertex $i$ for any $i$.)
This gives a pretty sparse $60\times 60$ matrix which is supplied
in MATLAB under the name {\tt bucky} (you don't {\tt load} this one, 
it's a predeclared name in MATLAB).  Note:  This example is also 
discussed in the text, but try to figure this out on your own before
you look there to check your work.  
\bigskip
\item{A)} (non-lab) From the picture of fullerene in Wikipedia above, or from 
your experience of soccer balls, (or if absolutely necessary 
by examining the entries in {\tt bucky} directly), how many nonzero 
entries are there on each column or row of $B$?  Why?
\item{B)} (non-lab) Prove that $B$ is {\it not positive definite} (no matter how
we label the vertices).  
\item{C)} To make a positive definite matrix, we can add a positive
multiple of the identity matrix $I_{60}$.  Since {\tt bucky} is a MATLAB
{\it sparse} matrix, we want a {\it sparse} identity matrix.  Compute 
the Cholesky factors for the matrices
$$B + 4I_{60}, \quad B + 10 I_{60}, \quad B + 20 I_{60}$$
and examine their {\tt spy} plots. 
Does the spy plot of the Cholesky factor of $B + c I_{60}$ appear
to depend on the value of $c$?  
\bigskip
\noindent
III.  We will consider the relation between
the spy plots of Cholesky factors $R$ and the spy plots of the original
positive definite matrices $A = R^tR$ in this question.
\bigskip
The {\it envelope} of a symmetric matrix $A$ or an upper triangular
matrix $U$
is a set of ordered pairs $(k,\ell)$ with $\ell \ge k$ representing the  
locations of the {\it nonzero} entries in the upper triangle of the matrix.
It is defined like this:  For each column $j = 1, \ldots, n$ we locate
the {\it first nonzero entry in that column} (counting down from the top 
of the matrix).  Say this happens in location $(m,j)$.
Then the rest of the pairs representing locations in that column 
down to the diagonal go into the envelope: 
$(m-1,j), (m-2,j), \ldots, (j-1,j)$.  
\bigskip
\item{A)}  On a printout spy plot of the bucky matrix $B$ from problem II, indicate
(by hand) where the envelope of $B$ is.  Do the same on the spy plot of the
Cholesky factor of $B + 4 I_{60}$.  
\item{B)}  (non-lab) The tridiagonal examples in question I and the example from
part A should make the following theorem plausible.

\proclaim Theorem. If $A$ is positive definite symmetric and $R$
is its Cholesky factor, then $A$ and $R$ have the same envelope.\par

\item{}  There is a proof of this in the text in Theorem 1.5.7 (page 57-58).
It depends on the ``bordered form'' of Cholesky factorization on page 43.
Your job in this part is to work through the details of that proof, and 
write out the complete argument in your own words.  In the course of this,
you will want to supply a proof for Exercise 1.4.35 in the text.  
\bye
