\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

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

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

  {\Large \textbf{Numerical Analysis -- Lecture
      1\footnote{Please email all corrections and suggestions to these
      notes 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{LU factorization of matrices}

\subsection{Definition and applications}

Let $A$ be a real $n\times n$ matrix. We say that the $n\times n$ matrices
$L$ and $U$ are an \textit{LU factorization\/} of $A$ if
\textbf{\textit{(1)\/}} $L$ is lower-triangular (i.e., $L_{i,j}=0$,
$i<j$); \textbf{\textit{(2)\/}} $U$ is upper-triangular, $U_{i,j}=0$,
$i>j$; and \textbf{\textit{(3)\/}} $A=LU$. Therefore the factorization
takes the form 
\vspace*{-8pt}
$$\left[ \begin{picture}(30,20)
\put (0,-13) {\line(1,0){30}}
\put (30,-13) {\line(0,1){30}}
\put (30,17) {\line(-1,0){30}}
\put (0,17) {\line(0,-1){30}}
         \end{picture}\right]= 
\left[ \begin{picture}(30,20)
\put (0,-13) {\line(1,0){30}}
\put (30,-13) {\line(-1,1){30}}
\put (0,17) {\line(0,-1){30}}
       \end{picture}\right] \times 
\left[\begin{picture}(30,20)
\put (0,17) {\line(1,0){30}}
\put (0,17) {\line(1,-1){30}}
\put (30,17) {\line(0,-1){30}}
      \end{picture}\right].$$

\vspace*{-4pt}
{\bf Application 1} Calculation of a determinant: $\det A=\det L
\det U=\left(\prod_{k=1}^n L_{k,k}\right)\cdot\left( \prod_{k=1}^n
  U_{k,k}\right)$. 

{\bf Application 2} Testing nonsingularity: $A=LU$ is nonsingular iff
all the diagonal elements of $L$ and $U$ are nonzero.

{\bf Application 3} Solution of linear systems: Let $A=LU$ and suppose
we wish to solve $A\MM{x}=\MM{b}$. This is the same as
$L(U\MM{x})=\MM{b}$, which we decompose into $L\MM{y}=\MM{b}$,
$U\MM{x}=\MM{y}$. Both latter systems are triangular and can be
calculated easily. Thus, $L_{1,1}y_1=b_1$ gives $y_1$, next
$L_{2,1}y_1+L_{2,2}y_2=b_2$ yields $y_2$ etc. Having found $\MM{y}$, we
solve for $\MM{x}$ by reversing the order: $U_{n,n}x_n=y_n$ gives
$x_n$, $U_{n-1,n-1}x_{n-1}+U_{n-1,n}x_n=y_{n-1}$ produces $x_{n-1}$
and so on. This requires $\mathcal{O}(n^2)$ computational
operations (usually we only bother to count multiplications/divisions). 

{\bf Application 4} The inverse of $A$: It is straightforward to
devise a direct way of calculating the inverse of triangular matrices,
subsequently forming $A^{-1}=U^{-1}L^{-1}$.

\vspace*{-4pt}
\subsection{The calculation of LU factorization}

We denote the {\em columns\/} of $L$ by
$\MM{l}_1,\MM{l}_2,\ldots,\MM{l}_n$ and the {\em rows\/} of $U$ by
$\MM{u}_1^{\CC{T}},\MM{u}_2^{\CC{T}},\ldots,\MM{u}_n^{\CC{T}}$. Hence
\begin{equation}\label{1.1}
  A=LU=[
  \begin{array}{cccc}
    \MM{l}_1 & \MM{l}_2 & \cdots & \MM{l}_n
  \end{array}
  ] \left[
    \begin{array}{c}
      \MM{u}_1^{\CC{T}}\\ \MM{u}_2^{\CC{T}}\\\vdots\\
      \MM{u}_n^{\CC{T}}  
    \end{array}
  \right]=\sum_{k=1}^n \MM{l}_k \MM{u}_k^{\CC{T}}.
\end{equation}
Since the first $k-1$ components of $\MM{l}_k$ and $\MM{u}_k$ are all
zero, each rank-one matrix $\MM{l}_k \MM{u}_k^{\CC{T}}$ has zeros
in its first $k$ rows and columns.

Assume that the factorization exists (hence the diagonal elements of
$L$ are nonzero) and that $A$ is nonsingular.
Since $\MM{l}_k\MM{u}_k^{\CC{T}}$ stays the same if we replace
$\MM{l}_k\rightarrow\alpha\MM{l}_k$, $\MM{u}_k\rightarrow
\alpha^{-1}\MM{u}_k$, where $\alpha\neq0$, we may assume w.l.o.g.\
that all diagonal 
elements of $L$ equal one. In  other words, the $k$th row of $\MM{l}_k
\MM{u}_k^{\CC{T}}$ is $\MM{u}_k^{\CC{T}}$ and its $k$th column is $U_{k,k}$ times 
$\MM{l}_k$.

We begin our calculation by extracting $\MM{l}_1$ and
$\MM{u}_1^{\CC{T}}$ from $A$, and then proceed similarly to extract
$\MM{l}_2$ and $\MM{u}_2^{\CC{T}}$, etc.

First we note that since the leading $k-1$ elements of $\MM{l}_k$ and
$\MM{u}_k$ are zero for $k\geq2$, it follows from \R{1.1} that
$\MM{u}_1^{\CC{T}}$ is the first row of $A$ and $\MM{l}_1$ is the
first column of $A$, divided by $A_{1,1}$ (so that $L_{1,1}=1$).

Next, having found $\MM{l}_1$ and $\MM{u}_1$, we form the matrix $A-\MM{l}_1
\MM{u}_1^{\CC{T}}=\sum_{k=2}^n \MM{l}_k \MM{u}_k^{\CC{T}}$. The first row \& 
column of $A$ are zero and it follows that $\MM{u}_2^{\CC{T}}$ is the second
row of $A-\MM{l}_1\MM{u}_1^{\CC{T}}$, while $\MM{l}_2$ is its second
column, scaled so that $L_{2,2}=1$.

{\bf The LU algorithm:} Set $A_0:=A$. For all $k=1,2,\ldots,n$ set
$\MM{u}_k^{\CC{T}}$ to the $k$th row of $A_{k-1}$ and $\MM{l}_k$ to the $k$th
column of $A_{k-1}$, scaled so that $L_{k,k}=1$. Further, calculate
$A_k:=A_{k-1} -\MM{l}_k\MM{u}_k^{\CC{T}}$ before incrementing $k$.

Note that all elements in the first $k$ rows \& columns of $A_k$ are
zero. Hence, we can use the storage of the original $A$ to accumulate
$L$ and $U$. The full LU factorization requires $\mathcal{O}(n^3)$
computational operations.

\vspace*{-4pt}
\subsection{Relation to Gaussian elimination}

The equation $A_k=A_{k-1}-\MM{l}_k\MM{u}_k^{\CC{T}}$ has the property that
the $j$th row of $A_k$ is the $j$th row of $A_{k-1}$ minus $L_{j,k}$ times
$\MM{u}_k^{\CC{T}}$ (the $k$th row of $A_{k-1}$). Moreover, the multipliers
$L_{k,k},L_{k+1,k},\ldots,L_{n,k}$ are chosen so that the outcome of
this {\em   elementary row operation\/} is that the $k$th column of
$A_k$ is zero. This construction is analogous to Gaussian 
elimination for solving $A\MM{x}=\MM{b}$. An important difference is
that in LU we do not consider the right hand side $\MM{b}$ until the
factorization is complete. This is useful e.g.\ when there are many
right hand sides, in particular if not all the $\MM{b}$'s are known
at the outset: in Gaussian
elimination the solution for each new $\MM{b}$ would require
$\mathcal{O}(n^3)$ computational operations, whereas with LU factorization
$\mathcal{O}(n^3)$ operations are required for the initial factorization, but
then the solution for each new $\MM{b}$ only requires $\mathcal{O}(n^2)$
operations.

\vspace*{-4pt}
\subsection{Pivoting}

Naive LU factorization fails when, for example, $A_{1,1}=0$. The
remedy is to exchange rows of $A$, a technique called {\em column
pivoting.\/} This is equivalent to picking a suitable equation for
eliminating the first unknown in Gaussian elimination. 
Specifically, column pivoting means that, having obtained
$A_{k-1}$, we exchange two rows of $A_{k-1}$ so that the element of
largest magnitude in the $k$th column is in the `pivotal position'
$(k,k)$. In other words, 
$$|(A_{k-1})_{k,k}|=\max\{ |(A_{k-1})_{j,k}|\, :\, j=1,2,\ldots,n\}.$$
Of course, the same exchange is required in the portion of $L$ that has
been formed already (i.e., the first $k-1$ columns). Also, we need to
record the permutation of rows to solve for the right hand side and/or
to compute the determinant. (The exchange of rows can be regarded as the
pre-multiplication of the relevant matrix by a permutation matrix.)

Column pivoting copes with zeros at the pivot position, except when the
whole $k$th column of $A_{k-1}$ is zero -- in that case it is usual to let
$\MM{l}_k$ be the $k$th unit vector while, as before, choose
$\MM{u}_k^{\CC{T}}$ as the $k$th row of $A_k$). Such a choice
preserves the condition that the matrix $\MM{l}_k \MM{u}_k^{\CC{T}}$
has the same $k$th row and column as $A_{k-1}$.  Thus $A_k := A_{k-1}
- \MM{l}_k \MM{u}_k^{\CC{T}}$ still has zeros in its $k$th row and
column as required. 

An important advantage of column pivoting is that every element of $L$
has magnitude at most one. This avoids not just division by zero but
also tends to reduce the chance of very large numbers occuring during
the factorization, a phenomenon that might lead to{\em ill
  conditioning\/} and to accumulation of {\em roundoff error.\/}

In {\em row pivoting\/} one exchanges columns of $A_{k-1}$, rather
than rows (sic!), whereas {\em total pivoting\/} corresponds to
exchange of both rows and columns, so that the modulus of the
pivotal element $(A_{k-1})_{k,k}$ is maximised.

\end{document}