\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

\textheight=22.5cm
\textwidth=15cm
\topmargin=-1cm
\oddsidemargin=0.5cm
\evensidemargin=0.5cm
\parindent=0pt
\parskip=6pt

\setcounter{section}{1}
\setcounter{subsection}{4}
\setcounter{equation}{1}

\begin{document}

\begin{center}
  {\large \textbf{Mathematical Tripos Part IB: Easter 2001}}

  {\Large \textbf{Numerical Analysis -- Lecture
      2\footnote{Corrections and suggestions to these notes should
        be emailed to \texttt{a.shadrin@damtp.cam.ac.uk}. All handouts
        are available on the WWW at the URL 
        \texttt{http://www.damtp.cam.ac.uk/user/na/PartIB/Handouts.html}.}}}
\end{center}

\section{Factorization of structured matrices}

\subsection{Symmetric matrices}

Let $A$ be an $n\times n$ symmetric matrix (i.e.,
$A_{k,\ell}=A_{\ell,k}$). An analogue of LU factorization takes advantage
of symmetry: we express $A$ in the form of the product $LDL^{\CC{T}}$,
where $L$ is $n\times n$ lower triangular, with ones on its diagonal,
whereas $D$ is a diagonal matrix. Subject to its existence,
we can write this factorization as
$$A=\left[\begin{array}{cccc}\MM{l}_1 & \MM{l}_2 & \cdots \MM{l}_n
\end{array}\right]\left[\begin{array}{cccc}
D_{1,1} & 0 & \cdots & 0\\ 
0 & D_{2,2} & \ddots & \vdots\\ 
\vdots & \ddots & \ddots & 0\\ 
0 & \cdots & 0 & D_{n,n} \end{array}\right]
\left[\begin{array}{c}\MM{l}_1^{\CC{T}}\\\MM{l}_2^{\CC{T}}
    \\\vdots\\\MM{l}_n^{\CC{T}} 
\end{array}\right]=\sum_{k=1}^n D_{k,k}\MM{l}_k\MM{l}_k^{\CC{T}}$$
where, as before, $\MM{l}_k$ is the $k$th column of $L$.\\
The analogy with the algorithm of Section 1.2 becomes obvious by letting
$U=DL^{\CC{T}}$, but the present form lends itself better to exploitation of
symmetry. Specifically, to compute this factorization, we let $A_0=A$
and for $k=1,2,\ldots,n$ let $\MM{l}_k$ be the multiple of the $k$th
column of $A_{k-1}$ such that $L_{k,k}=1$. Set
$D_{k,k}=(A_{k-1})_{k,k}$ and form $A_k=A_{k-1}-D_{k,k}\MM{l}_k
\MM{l}_k^{\CC{T}}$. 

{\bf Example} Let
$A=A_0=\left[\begin{array}{cc}2&4\\4&11\end{array}\right]$. Hence
$\MM{l}_1=\left[\begin{array}{c}1\\2\end{array}\right]$, $D_{1,1}=2$
and
$$A_1=A_0-D_{1,1}\MM{l}_1\MM{l}_1^{\CC{T}}=\left[\begin{array}{cc}2&4\\4&11
\end{array}\right]-2\left[\begin{array}{cc}1&2\\2&4 \end{array}\right]
=\left[\begin{array}{cc}0&0\\0&3\end{array}\right].$$
We deduce that
$\MM{l}_2=\left[\begin{array}{c}0\\1\end{array}\right]$, $D_{2,2}=3$
and 
$A=\left[\begin{array}{cc} 1&0\\2&1\end{array}\right] \left[\begin{array}{cc}
2&0\\0&3\end{array}\right]\left[\begin{array}{cc}1&2\\0&1\end{array}\right]$.

\subsection{Symmetric positive definite matrices}

Recall that $A$ is positive definite if $\MM{x}^{\CC{T}}A\MM{x}>0$ for
all $\MM{x}\neq{\bf 0}$. 

{\bf Theorem} Let $A$ be a real $n\times n$ symmetric matrix. It is
positive definite if and only if it has an $LDL^{\CC{T}}$
factorization in which the diagonal elements of $D$ are all positive.\\
\hspace*{8pt}{\bf Proof.} Suppose that $A=LDL^{\CC{T}}$ and let
$\MM{x}\in\BB{R}^n\setminus\{\MM{0}\}$. Since $L$ is nonsingular,
$\MM{y}:=L^{\CC{T}}\MM{x}\neq\MM{0}$. Then
$\MM{x}^{\CC{T}}A\MM{x}=\MM{y}^{\CC{T}}D\MM{y}=\sum_{k=1}^n D_{k,k}y_k^2>0$, hence
$A$ is positive definite.\\ 
Conversely, suppose that $A$ is positive definite. We wish to
demonstrate that an $LDL^{\CC{T}}$ factorization exists. We denote by
$\MM{e}_k\in\BB{R}^n$ the $k$th unit (a.k.a.\ coordinate)
vector. Hence $\MM{e}_1^{\CC{T}}A\MM{e}_1=A_{1,1}>0$ and $\MM{l}_1$ \&
$D_{1,1}$ are well 
defined. We now show that $(A_{k-1})_{k,k}>0$ for $k=1,2,\ldots$. The
result is true for $k=1$ and we continue by induction (hence may
assume that $A_{k-1}=A-\sum_{j=1}^{k-1}D_{j,j}\MM{l}_j\MM{l}_j^{\CC{T}}$
has been computed successfully).\\
We define $\MM{x}\in\BB{R}^n$ as follows. The bottom $n-k$ components
are zero, $x_k=1$ and $x_1,x_2,\ldots,x_{k-1}$ are calculated in a
reverse order, each $x_j$ being chosen so that $\MM{l}_j^{\CC{T}}\MM{x}=0$
for $j=k-1,k-2,\ldots,1$. In other words, since
$0=\MM{l}_j^{\CC{T}}\MM{x}=\sum_{i=1}^n L_{i,j}x_i=\sum_{i=j}^k L_{i,j}x_i$,
we let $x_j=-\sum_{i=j+1}^k L_{i,j}x_i$, $j=k-1,k-2,\ldots,1$.\\
Since the first $k-1$ rows \& columns of $A_{k-1}$ vanish, our choice
implies that $(A_{k-1})_{k,k}=\MM{x}^{\CC{T}}A_{k-1}\MM{x}$. Thus, from the
definition of $A_{k-1}$ and the choice of $\MM{x}$,
$$(A_{k-1})_{k,k}=\MM{x}^{\CC{T}}A_{k-1}\MM{x}=\MM{x}^{\CC{T}}
\left(A-\sum_{j=1}^{k-1}D_{j,j} \MM{l}_j\MM{l}_j^{\CC{T}}\right)
  \MM{x}=\MM{x}^{\CC{T}}A\MM{x}-\sum_{j=1}^{k-1} D_{j,j} (\MM{l}_j^{\CC{T}}\MM{x})^2
  =\MM{x}^{\CC{T}}A\MM{x}>0,$$
as required. Hence $(A_{k-1})_{k,k}>0$, $k=1,2,\ldots,n$, and the
factorization exists. \QED

{\bf Conclusion} It is possible to check if a symmetric matrix is
positive definite by trying to form its $LDL^{\CC{T}}$ factorization.

{\bf Cholesky factorization} Define $D^{1/2}$ as the diagonal
matrix whose $(k,k)$ element is $D_{k,k}^{1/2}$, hence
$D^{1/2}D^{1/2}=D$. Then, $A$ being positive definite, we can
write
$$A=(LD^{1/2})(D^{1/2}L^{\CC{T}})=(LD^{1/2})(LD^{1/2})^{\CC{T}}.$$
In other words, letting $\tilde{L}:=LD^{1/2}$ by $L$, we obtain the {\em
Cholesky factorization\/} $A=\tilde{L}\tilde{L}^{\CC{T}}$.

\subsection{Sparse matrices}

Frequently it is required to solve {\em very\/} large systems
$A\MM{x}=\MM{b}$ ($n=10^5$ is considered small in this context!) where
nearly all the elements of $A$ are zero. Such a matrix is called {\em
sparse\/} and efficient solution of $A\MM{x}=\MM{b}$ should exploit
sparsity. In particular, we wish the matrices $L$ and $U$ to inherit
as much as possible of the sparsity of $A$.  The only tool at our
disposal at the moment is the freedom to exchange rows and columns to
minimise {\em fill-in.\/} To this end the following theorem is useful.

{\bf Theorem} Let $A=LU$ be an LU factorization (without pivoting) of a
sparse matrix. Then all leading zeros in the rows of $A$ to the left of the
diagonal are inherited by $L$ and all the leading zeros in the columns of $A$
above the diagonal are inherited by $U$.

{\bf Proof} Follows from {\sf Question 2} on {\sf Examples'
Sheet 1}. \QED

This theorem suggests that if one requires a factorization of a sparse
matrix then one might try to reorder its rows and columns by a preliminary
calculation so that many of the zero elements are leading zero elements in
rows and columns. This will reduce the fill-in.

\textbf{Example 1} The LU factorisation of
$$\left[\begin{array}{ccccc}
-3 & 1 & 1 & 2 & 0\\
1 & -3 & 0 & 0 & 1\\
1 & 0 & 2 & 0 & 0\\
2 & 0 & 0 & 3 & 0\\
0 & 1 & 0 & 0 & 3
        \end{array}\right]=\left[\begin{array}{ccccc}
1 & 0 & 0 & 0 & 0\\
-\frac13 & 1 & 0 & 0 & 0\\
-\frac13 & -\frac18 & 1 & 0 & 0\\
-\frac23 & -\frac14 & \frac{6}{19} & 1 & 0\\
0 & -\frac38 & \frac{1}{19} & \frac{4}{81} & 1  \end{array}\right]
\left[\begin{array}{ccccc}
-3 & 1 & 1 & 2 & 0\\
0 & -\frac83 & \frac13 & \frac23 & 1\\
0 & 0 & \frac{19}{8} & \frac34 & \frac18\\
0 & 0 & 0 & \frac{81}{19} & \frac{4}{19}\\
0 & 0 & 0 & 0 & \frac{272}{81}
      \end{array}\right],$$
has significant fill-in. However, reordering (symmetrically) rows and
columns yields
$$\left[\begin{array}{ccccc}
2 & 0 & 1 & 0 & 0\\
0 & 3 & 2 & 0 & 0\\
1 & 2 & -3 & 0 & 1\\
0 & 0 & 0 & 3 & 1\\
0 & 0 & 1 & 1 & -3
        \end{array}\right]=\left[\begin{array}{ccccc}
1 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
\frac12 & \frac23 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & -\frac{6}{29} & \frac13 & 1 \end{array}\right]
\left[\begin{array}{ccccc}
2 & 0 & 1 & 0 & 0\\
0 & 3 & 2 & 0 & 0\\
0 & 0 & -\frac{29}{6} & 0 & 1\\
0 & 0 & 0 & 3 & 1\\
0 & 0 & 0 & 0 & -\frac{272}{87}
      \end{array}\right].$$

\textbf{Example 2} If the nonzeros of $A$ occur only on the diagonal, in one
row and in one column, then the full row and column should
be placed at the bottom and on the right of~$A$, respectively.


\end{document}