% FILE antonNM.tex

% AR6 Chapter 9. AR7 Section 9.9
\section*{\hNM\ Introduction to Numerical Methods of Linear Algebra}

%Further notes on Maple and numerical linear algebra
%(and other notes on ``Maple and floating point numerical computation")
%are available from G. Keady.
%One version of the numerical linear algebra notes is organised
%as an accompaniment to the relevant chapter in
%E. Kreyszig {\it Advanced Engineering Mathematics} (Wiley).

In this Chapter we basically confine ourselves to a few of the
items treated in \cite{AR}, 
and relating to Gaussian elimination methods
of solving linear systems.
See also the chapter on numerical linear algebra in
\cite{K} for another account suitable for engineering maths
students.

The `numerical methods' refer to methods for machine-precision
(SmallFloat)
floating point computation. 
Except for small problems, in practice one does not use Maple
for this sort of task.

\ifodd\ARed\setcounter{subsection}{7}\else\setcounter{subsection}{0}\fi
\subsection{Comparison of Procedures for Solving Linear Systems}

In the problem of SmallFloat approximate solution of 
linear systems $Ax=b$, one does not find $A^{-1}$, and nor
does one use Cramer's rule. The best methods are essentially
Gaussian elimination.


\subsection{$LU$-Decompositions}

\headq{Definitions.} A square matrix is said to be
\begin{enumerate}
\item[] {\it upper triangular} if all entries below the main diagonal are zero
  (i.e. it is in row echelon form).
\item[] {\it lower triangular} if all entries above the main diagonal are zero. 
\item[] {\it unit (lower or upper) triangular} if it is (lower or upper)
  triangular matrix and its diagonal entries are all one.
\end{enumerate}
For example, the elementary matrices $E_{ij}$ (of Chapter 1)
corresponding to the Maple command {\tt addrow} are unit lower triangular.
\medskip

A result, which we treated in Chapter 1 merely as an example of
matrix arithmetic, is as follows:

%\setcounter{theorem}{55}
%\begin{theorem}
%\headq{THEOREM }
\theoremNN{\it
  {\em (a)} The product of two lower triangular matrices is lower triangular.

  \qquad {\em (b)} The product of two unit lower triangular matrices is unit
  lower triangular.
} %\end{theorem}
\theoremNNEnd

% EXERCISE BELOW IS IN Chapter 1
%\headq{EXERCISE.} Verify the preceding theorem for 
%general $3\times 3$ lower triangular matrices and\newline
%\quad for general $3\times 3$ unit lower triangular matrices.
\medskip

%\setcounter{theorem}{57}
%\begin{theorem}[LU-Decomposition]
%\headq{THEOREM }
\theoremNN{\bf (LU-Decomposition).}{\it
  If the Gaussian elimination of the square matrix $A$ requires no row
  interchanges, yielding nonzero pivots, then there are unique matrices
  $L$, $U$ such that
    $$A=LU$$ 
  where $U$ is an upper triangular matrix and $L$ is a unit lower triangular
  matrix.
} %\end{theorem}
\theoremNNEnd

\label{lufac}
\begin{maple}
# Example 2 - AR6 Section 9.2 p442; AR7 Section 9.9 p512
# CAUTION. Anton makes U unit upper triangular;
# below we make L unit lower triangular.
> with(linalg):
# Find the LU decomposition of the matrix
> A:= matrix([ [ 2,  6,  2],
               [-3, -8,  0],
               [ 4,  9,  2] ]):
# Set the diagonal entries of  L  equal to 1 and  U  equal to A
> L:= diag(1,1,1):                               U:= A:
#
# Now reduce  U  to row-echelon form, storing the INVERSES of the row
# operations used in the matrix  L
\end{maple}
\goodbreak
\begin{maple}
> L[2,1]:= U[2,1]/U[1,1]:               U:= addrow( U,1,2,-L[2,1]):
  L[3,1]:= U[3,1]/U[1,1]:               U:= addrow( U,1,3,-L[3,1]):
  L[3,2]:= U[3,2]/U[2,2]:               U:= addrow( U,2,3,-L[3,2]):
  print( L,U );
                        [   1    0  0 ]  [ 2  6  2 ]
                        [             ]  [         ]
                        [ -3/2   1  0 ], [ 0  1  3 ]
                        [             ]  [         ]
                        [   2   -3  1 ]  [ 0  0  7 ]

# Verify that    A = L*U
> equal( A,evalm(L&*U) );
                                      true
\end{maple}
\goodbreak

\medskip
Not every matrix $A$ has a factorization $A=LU$.\quad
Consider, for example, the matrix
  $$ A = \pmatrix{0& a12\cr
                a21& a22\cr }\qquad\qquad {\rm where}~a21,a12\ne 0. $$
However systems of equations involving this last $2\times 2$ matrix are
still easy to solve using row interchanges (swaps) and permutation matrices.
\smallskip

While discussing $LU$-factorisation, we choose to look ahead to
material which in \cite{AR} is discussed in 
Exercise 17 of this section 
and is related to [AR6] Section 9.4.
\smallskip

\headq {Definition.} An {\it elementary permutation matrix} is a matrix
$P_{ij}$ corresponding to the Maple command {\tt swaprow}, and a
{\it permutation matrix} is a product of such matrices.
\smallskip

\headq{EXERCISE.} Show that the product of two permutation matrices 
is a permutation matrix. Also show that any permutation matrix is invertible
and its inverse is also a permutation matrix.
\medskip
 
%\setcounter{theorem}{62}
%\begin{theorem} 
%\headq{THEOREM }
\theoremNN{\it 
  Let $A$ be an invertible matrix.
  Then there is a permutation matrix $P$ such that
    $$ P A = L U ,$$
  where $L$ is unit lower triangular, and $U$ is upper triangular
  (both depending on $P$).
} %\end{theorem}
\theoremNNEnd

\begin{maple}
# Exercise 18 - AR6 Section 9.2 p449; AR7 Section 9.9 p519
# An example of the  PA = LU  decomposition of a matrix
> with(linalg):
> A:= matrix([ [3,-1, 0],
               [3,-1, 1],
               [0, 2, 1] ]):
# The following row operations will reduce  A  to Echelon form:
L:= diag(1,1,1):                               U:= swaprow(A,2,3):
L[3,1]:= U[3,1]/U[1,1]:               U:= addrow( U,1,3,-L[3,1]);
                                   [ 3  -1  0 ]
                                   [          ]
                              U := [ 0   2  1 ]
                                   [          ]
                                   [ 0   0  1 ]


# Alternatively, perform the above row interchanges on the identity matrix
> L:= diag(1,1,1): P23:= swaprow(L,2,3);
# Now perform an LU decomposition on PA, as in LU Example above
> PA:= evalm( P23&*A );
# Looking at the result on the screen, and working the steps
# one gets to L
   L:=matrix([[1,0,0],[0,1,0],[1,0,1]]):
   equal(A,P23 &* L &* U);
                                     true

> equal( PA,evalm(L&*U) );
                                     true
\end{maple}
%\headq{EXERCISES.}
%\begin{enumerate}
%\item Complete the above example 7.
%\item Write Maple procedures to
%  \begin{enumerate}
%  \item do forward substitution;
%  \item find the LU-decomposition of a matrix $A$;
%  \item replace a matrix $A$ by $PA$ and do an LU-decomposition of $PA$.
%  \end{enumerate}
%\end{enumerate}

The point of $LU$-decomposition is that it speeds up the solution of many linear
systems $A\b x =\b b_1$, $A\b x =\b b_2$,\ldots $A\b x =\b b_m$, where the
coefficient matrix is the same each time.


\subsubsection*{Transposes, Symmetry, and Band Matrices}
There are various applications of symmetric matrices.
It is also worth noting that the symmetry of a matrix
may be used when storing it, since only the values on and below the diagonal
need be entered.  Maple takes advantage of this by allowing you to type,
for example
\begin{maple}
> A:= array( symmetric,1..3,1..3 ):
> with(linalg):
> equal( A,transpose(A) );
                                      true
\end{maple}

A matrix $A$ is said to be {\it banded with bandwidth $2k+1$} if $a_{ij}=0$ 
whenever $\|i-j\|>k$.
Banded matrices occur in practice.
To avoid conflict with Maple's {\tt band} function, we use the word
{\em banded} where many linear algebra texts use {\em band}.  

The simplest examples (after diagonal matrices) of banded matrices
are tridiagonal matrices.
A matrix $A$ is said to be {\it tridiagonal} if $a_{ij}=0$ 
whenever $\|i-j\|>1$.

%In the Maple and numerical linear algebra document,
%the derivation of the 
There is a special structure to the Choleski factorization 
$L D L^T$ of a tridiagonal matrix $A$:
the $L$ matrix is tridiagonal (as well as unit lower triangular)


\medskip
Again we stress that, in solving large sets of
linear equations in floats, ONE DOES NOT USE {\tt det} OR {\tt inverse},
but uses matrix factorization techniques.

This is especially obvious in the case of banded matrices.  The
number of operations for factoring such an $n\times n$ banded matrix is $O(n)$,
and for solving the equations $L\b y = \b b$ and $L^T\b x =\b y$ is also $O(n)$.
The inverse matrix is generally not a banded matrix. Forming $A^{-1}\b b$ is
at least $O(n^2)$ -- i.e. many, many more operations.

\smallskip
Interestingly, for exact solution in rationals the situation is different
in that rationals themselves have different complexities, and this means
that multiplications can have different costs.  At this stage,
let's agree to leave it to the package designers and other specialists.
It is rare in practice to want to solve very large systems exactly. 

\medskip\goodbreak
\subsubsection*{Numerical and Programming Considerations}
Maple is an example of a Computer Algebra System, which gives exact solutions
to problems.
However, most problems cannot be solved exactly, and even if they could,
in practice one only wants approximations to answers or graphics.
For this reason, the advances in computing have led to
`numerical analysis' taking an increasingly important role in science. 
In practice numerical libraries (such as NAG) are very important.

In situations where greater accuracy is needed than double-precision, 
the arbitrary-precision capabilities of Maple can be useful.  For example,
the ability to reset the value of the Maple variable {\tt Digits}
allows the effect of round-off error to be explored.

There are two sorts of error in floating-point numeric computation: truncation
error and roundoff error. The first is the error associated with truncating
some infinite process to a finite one.
Truncation errors do arise in connection with the
iterative methods of solving linear systems such as the
Gauss-Siedel and Jacobi methods.

