% Hill1.tex

\pagenumbering{arabic}

\section{Introduction to Linear Equations and Matrices}

\subsection{Introduction to Linear Systems and Matrices}
The central problem of linear algebra is the solution of linear equations.
You should already have some experience of these, and in particular you should
know that any system of linear equations can have 0, 1 or an infinite number of
solutions.

Maple has a command {\tt linsolve}\index{linsolve} for solving linear systems.
To illustrate this, as well as seeing an application of matrix arithmetic,
we look at Example 4 of Hill Section 1.1, page 4.
\medskip

\begin{verbatim}
# Hill Section 1.1 Example 4, p4
# An example of a Markov process, allowing determination of long-term trends
with(linalg):
# Suppose that the population of U.S.A. remains stable at 270 million. Also
# each year     20 per cent of people living in Texas move out of Texas
#               10 per cent of those outside Texas move into Texas.
#
# (a) If y1, y2 are the populations inside and outside Texas respectively
#     at the beginning of a year,  and  u1, u2 are the populations inside and
#     outside Texas respectively at the end of that year, then we can write
#                       u = T y,      where  u = [u1,u2],  y = [y1,y2]  and
T:= matrix([ [80/100, 10/100],
             [20/100, 90/100] ]):
#
# (b) At the beginning of 1991, the populations are given by y where
y:= vector( [50000000,220000000] ):
#
#     At the end of 1991, the populations are u1, u2 where u = [u1,u2] and
> u:= evalm( T&*y );
                          u := [ 62000000, 208000000 ]

# (c) To find the populations at the beginning of 1990, solve  T x = y:
> x:= linsolve(T,y);
                       x := [ 230000000/7, 1660000000/7 ]

# (d) If the populations  z  remain the same for two consecutive years, then
#                        z = T z
#               (Id - T) z = 0,         where
> Id:= diag( 1,1):
> z:= linsolve( Id-T, vector([0,0]) );
                            z := [ _t[1], 2 _t[1] ]

#    (In Maple V,  _t[1]  was just  t1 ).
# Since the total population is 270 million  =  _t[1] + 2*_t[1]  =  3*_t[1],
#    the populations are _t[1] = 90 million and 2*_t[1] = 180 million.
#
# (e) To find out what happens over the six year period from 1990 to 1995,
#     repeat (b) six times  or use the following loop:
> for i from 1 to 6 do    x:= evalm( T&*x );    od;
\end{verbatim}
\goodbreak
\begin{verbatim}
                          x := [ 50000000, 220000000 ]

                          x := [ 62000000, 208000000 ]

                          x := [ 70400000, 199600000 ]

                          x := [ 76280000, 193720000 ]

                          x := [ 80396000, 189604000 ]

                          x := [ 83277200, 186722800 ]

# These answers are tending (slowly!) toward the answer  z  in (d).
\end{verbatim}

Maple also has a command {\tt augment}\index{augment} for forming the
augmented matrix of a system of equations, and {\tt backsub}\index{backsub}
for solving a system in row-echelon form:

\begin{verbatim}
# Hill Section 1.1 Example 5, p5
> with(linalg):
# The coefficient and augmented matrices of the following linear system
#        2 x + 3 y - 2 z = -7,
#              3 y + 2 z =  3,
#                    5 z = 15
# (which is in row-echelon form)  are given by
Coeff:= matrix([ [2, 3,-2],
                 [0, 3, 2],
                 [0, 0, 5] ]):
linsys:= augment( Coeff, vector([-7,3,15]) );
                                     [ 2  3  -2  -7 ]
                                     [              ]
                           linsys := [ 0  3   2   3 ]
                                     [              ]
                                     [ 0  0   5  15 ]

# and the solution is given by
> backsub(");
                                  [ 1, -1, 3 ]
\end{verbatim}

\goodbreak
\subsection{Gaussian Elimination}
In the last section, we mentioned Maple's command, {\tt linsolve}.
However, it remains worthwhile to know that row operations  
underpin many algorithms for solving linear systems. 
Accordingly we follow Hill in looking at row operations.
In fact, for floating point numeric computation on ordinary serial
computers, Gaussian elimination is the best method for large systems.

\begin{verbatim}
# Hill Section 1.2 Example 4, p14
# This concerns the "Gaussian elimination" method of solving linear systems. 
> with(linalg):
\end{verbatim}
\goodbreak
\index{addrow}
\begin{verbatim}
# Consider the linear system of equations
#  2*x1 -   x2 -   x3 =  3 ,
# -6*x1 + 6*x2 + 5*x3 = -3 ,
#  4*x1 + 4*x2 + 7*x3 =  3 .
# This system can be written using an augmented matrix A_b = [A | b] where
> A:= matrix([ [ 2,-1,-1],
>              [-6, 6, 5],
>              [ 4, 4, 7] ]):
> b:= vector( [3,-3,3] ):                                 A_b:= augment(A,b);
                                  [  2  -1  -1   3 ]
                                  [                ]
                           A_b := [ -6   6   5  -3 ]
                                  [                ]
                                  [  4   4   7   3 ]

# The example begins by reducing A_b to row-echelon form with row operations.
# It should be emphasised that these operations are likely to be more useful
# in learning some matrix manipulative processes than they are in actually
# solving equations, as the latter is now built into Maple's  linsolve.
# 
> pivot:= A[1,1]:          # The pivot is the first non-zero entry in the row
# 
# The command to add  -A[2,1]/pivot  times the first row to the second row is
> A1:= addrow( A_b,1,2,-A[2,1]/pivot ):
# and similarly to add  -A[3,1]/pivot  times the first row to the third row:
> A2:= addrow( A1, 1,3,-A[3,1]/pivot ); 
                                  [ 2  -1  -1   3 ]
                                  [               ]
                            A2 := [ 0   3   2   6 ]
                                  [               ]
                                  [ 0   6   9  -3 ]

> pivot:= A2[2,2]:                     Ech:= addrow( A2, 2,3,-A2[3,2]/pivot ); 
                                  [ 2  -1  -1   3  ]
                                  [                ]
                           Ech := [ 0   3   2   6  ]
                                  [                ]
                                  [ 0   0   5  -15 ]

# If you have had some experience with Maple procedures,
# you could put the above row operations into a procedure.
# (If you are just beginning your Maple you should ignore this next part.)
#
> noswaps:=proc(A)
>   local pivot,A1,A2,Ech;
>   pivot:= A[1,1]:
>   A1:=addrow(  A,1,2,-A[2,1]/pivot ):   A2:=addrow( A1,1,3,-A1[3,1]/pivot ):
>   pivot:= A2[2,2]:          Ech:= addrow( A2, 2,3,-A2[3,2]/pivot ); 
> end:
\end{verbatim}
\goodbreak
\begin{verbatim}
#
> equal( noswaps(A_b),Ech );
                                      true

# showing that the result of  noswaps  is  Ech.  The answer to the system is
> backsub(Ech);
                                  [ 2, 4, -3 ]
\end{verbatim}

\noindent The case where row swaps are required is given in Example 5 of Hill
p.16 below.\index{swaprow}

\begin{verbatim}
# Hill Section 1.2 Example 5, p16
> with(linalg):
# Consider the linear system of equations
#               6*x3 + 2*x4 - 4*x5 - 8*x6 = 8 ,
#               3*x3 +   x4 - 2*x5 - 4*x6 = 4 ,
# 2*x1 - 3*x2 +   x3 + 4*x4 - 7*x5 +   x6 = 2 ,
# 6*x1 - 9*x2        +11*x4 -19*x5 + 3*x6 = 1 .
# The associated augmented matrix is
> A_b:= matrix([ [ 0, 0, 6,  2, -4,-8, 8],
>                [ 0, 0, 3,  1, -2,-4, 4],
>                [ 2,-3, 1,  4, -7, 1, 2],
>                [ 6,-9, 0, 11,-19, 3, 1] ]):
# In this case, the pivot is in the third row, and so start by
# interchanging the first and third rows:
> A1:= swaprow(A_b,1,3);
                             [ 2  -3  1   4   -7   1  2 ]
                             [                          ]
                             [ 0   0  3   1   -2  -4  4 ]
                       A1 := [                          ]
                             [ 0   0  6   2   -4  -8  8 ]
                             [                          ]
                             [ 6  -9  0  11  -19   3  1 ]

# and then proceed as before:
> A2:= addrow(A1,1,4,-A1[4,1]/A1[1,1]):   A3:= addrow(A2,2,3,-A2[3,3]/A2[2,3]):
> A4:= addrow(A3,2,4,-A3[4,3]/A3[2,3]);
                             [ 2  -3  1  4  -7   1   2 ]
                             [                         ]
                             [ 0   0  3  1  -2  -4   4 ]
                       A4 := [                         ]
                             [ 0   0  0  0   0   0   0 ]
                             [                         ]
                             [ 0   0  0  0   0  -4  -1 ]

# Since all the zero rows of an Echelon form must be at the bottom,
# we finally interchange the last two rows:
\end{verbatim}
\goodbreak
\begin{verbatim}
> Ech:= swaprow(A4,3,4);
                              [ 2  -3  1  4  -7   1   2 ]
                              [                         ]
                              [ 0   0  3  1  -2  -4   4 ]
                       Ech := [                         ]
                              [ 0   0  0  0   0  -4  -1 ]
                              [                         ]
                              [ 0   0  0  0   0   0   0 ]

# The  gausselim  function in the linalg package also gives this matrix
> equal( gausselim(A_b),Ech );
                                      true
\end{verbatim}

\headq{An optional remark.} There are some theoretical
(but not practical numerical) uses in continuing the elimination process
until the matrix is in a still simpler {\bf row-reduced echelon form}.
The difference is that all pivots are normalized to 1 by dividing each row
by its pivot, and zeroes are produced not only below but also above each pivot.
See Appendix B of Hill for more details.

\noindent The Maple command for this is\quad {\tt rref}\index{rref}\quad or
\quad {\tt gaussjord}\index{gaussjord}.

\goodbreak
\subsection{The Algebra of Matrices}
These have been considered in the opening section on the  {\tt linalg} package.
\medskip

\headq{EXERCISES}
\nopagebreak
\begin{enumerate}
\item Let $R$, $S$ be two arbitrary $2\times 2$ matrices, i.e.

  \verb$ R:= matrix(2,2);                        S:= matrix(2,2);$
  \begin{enumerate}
    \item Use Maple to calculate the products $RS$ and $SR$. \\
      \quad (Check that you could find them without Maple too!)
    \item Use the function {\tt equal} to discover if these products are the
      same.
  \end{enumerate}
  (The systematic study of matrices, their noncommutative multiplication, 
  etc., started to interest people relatively recently, only last century.)

\item Use the {\tt linalg} multiplication to CHECK the answers in Hill
  \S1.3 Example 7 (p27).
\end{enumerate}

\noindent The following example illustrates the different ways of viewing matrix
multiplication:

\begin{verbatim}
# Hill Section 1.3 Examples 8 - 10 (pp. 29-32).
> with(linalg):
# These examples deal with the system of equations
#        2x - 3y + 4z =  5
#        4x      -  z = -2
# whose coefficient matrix is
> A:= matrix([ [2,-3, 4],
>              [4, 0,-1] ]):
# Let
> X:= vector( [x,y,z] ):       AX:= evalm( A&*X ):       B:= vector( [5,-2] ):
> evalm(AX) = evalm(B);
                    [ 2 x - 3 y + 4 z, 4 x - z ] = [ 5, -2 ]

# Thus  A*X = B  is the matrix equation of the system of equations.
\end{verbatim} 
\goodbreak
\begin{verbatim}
#
# Now make the substitution
> lprint();  subs( x = 3*s-t, y = 2*s+5*t, z = -s+2*t, " );

                    [ - 4 s - 9 t, 13 s - 6 t ] = [ 5, -2 ]

# The coefficient matrix of this substitution is
> C:= matrix([ [ 3,-1],
>              [ 2, 5],
>              [-1, 2] ]):
# and the coefficient matrix of the new system (in s and t) is
> AC:= evalm( A&*C );
                                      [ -4  -9 ]
                                AC := [        ]
                                      [ 13  -6 ]

#
# The relationship between matrix multiplication and the columns of a matrix
#
# Example 9 (p 31).
> A1:= col(A,1):              A2:= col(A,2):              A3:= col(A,3):
> evalm( x*A1 + y*A2 + z*A3 - AX );
                                    [ 0, 0 ]

# This shows that  AX  is the linear combination  x*A1 + y*A2 + z*A3.
#
# Example 10 (p 32).
> equal( evalm( augment( A&*col(C,1), A&*col(C,2) ) ), AC );
                                      true

# This shows that the j'th column of A*C equals A times (j'th column of C)
#
> row1:= submatrix(A,1..1,1..3):                 row2:= submatrix(A,2..2,1..3):
> equal( evalm( stack(row1&*C,row2&*C) ), AC );
                                      true

# This shows that the i'th row of A*C equals (i'th row of A) times C
\end{verbatim} 
\medskip

\headq{EXERCISE} Let $R$ be the matrix in the previous exercise, and
\begin{enumerate}
\item[] \verb$ C:= diag(d1,d2);$
  \begin{enumerate}
  \item Show that the first column of $RC$ is $d1$ times the first column of
    $R$ and the second column of $RC$ is $d2$ times the second column of $R$.
  \item Show that the first row of $CR$ is $d1$ times the first row of
    $R$ and the second row of $CR$ is $d2$ times the second row of $R$.
  \item Use the Maple command\quad \verb$evalm$\quad to evaluate $C^n$
    for $n=2,3,4$.\quad Do you notice any pattern?
  \end{enumerate}
\end{enumerate}
This exercise illustrates that diagonal matrices are particularly easy to
manipulate, as will become increasingly clear later.

\goodbreak
\subsection{Inverses and Elementary Matrices}
\headq{EXERCISE}\quad Use the\quad {\tt inverse}\quad command in the\quad
{\tt linalg}\quad package to
\begin{enumerate}
\item[]
  \begin{enumerate}
  \item check the matrix inverses given in Hill \S1.4 Examples 2 and 3
    (pp. 38,39);
  \item verify that\quad $(RS)^{-1} = S^{-1}R^{-1}$\quad
    but\quad $(RS)^{-1} \ne R^{-1}S^{-1}$,\quad
    where $R$ and $S$ are the matrices of the previous exercises,
  \item find the inverse of the diagonal matrix $C$ of the previous exercise.
  \end{enumerate}
\end{enumerate}

\setcounter{theorem}{46}
\begin{theorem} Let $C$ be the $n\times n$ diagonal matrix
  $$ C = \pmatrix{d_1 \cr & d_2\cr & & \ddots\cr & & & d_n\cr }. $$
\begin{enumerate} \begin{enumerate}
  \item If $A$ is an $n\times p$ matrix and $B$ is an $m\times n$ matrix, then
    the $i$'th row of the matrix $CA$ is $d_i$ times the $i$'th row of $A$ and
    the $i$'th column of $BC$ is $d_i$ times the $i$'th column of $B$.
  \item $C$ is invertible if and only if every $d_i$ is non-zero. In this case
 $$C^{-1}=\pmatrix{d_1^{-1}\cr & d_2^{-1}\cr & & \ddots\cr & & & d_n^{-1}\cr}.$$
  \end{enumerate} \end{enumerate}
\end{theorem}

\medskip
The elementary matrices can be obtained in Maple by using the row operations
{\tt addrow, mulrow, swaprow} on the identity matrix:\index{mulrow}

\begin{verbatim}
# Hill Section 1.4 Example 5, p41
> with(linalg):
> I3:= diag( 1,1,1):                                   M1:= mulrow(I3,1,c):
> P23:= swaprow(I3,2,3):                             E13:= addrow(I3,1,3,m):
> print( M1,P23,E13 );
                     [ c  0  0 ]  [ 1  0  0 ]  [ 1  0  0 ]
                     [         ]  [         ]  [         ]
                     [ 0  1  0 ], [ 0  0  1 ], [ 0  1  0 ]
                     [         ]  [         ]  [         ]
                     [ 0  0  1 ]  [ 0  1  0 ]  [ m  0  1 ]
\end{verbatim}
\headq{Note} The matrices {\tt M1, P23, E13}
will be printed as $M_1$, $P_{23}$, $E_{13}$ respectively below.
\smallskip

\headq{EXERCISE} Perform the same row operations on the matrix
\begin{verbatim}
A:= matrix([ [ 2,-2,3,1],
             [ 1, 0,2,4],
             [-6, 2,1,5] ]);
\end{verbatim}
and show that the results are the same as $M_1 A$, $P_{23} A$ and $E_{13} A$
respectively. 

\setcounter{theorem}{42}
\begin{theorem}
  Each elementary matrix is invertible, and its inverse is an elementary matrix
  of the same type.
\end{theorem}

\goodbreak
\headq{EXERCISES}
\nopagebreak
\begin{enumerate}
\item Verify the preceding theorem for the matrices $M_1$, $P_{23}$ and $E_{13}$.
\item Verify, for the matrix $P_{23}$ above, that\quad $P_{23}^2 = I$\quad
  and so\quad $P_{23}^{-1} = P_{23}$. 
\item Let $A$ be the matrix in the previous exercise, and $E_{12}$,
  $E_{13}$, $E_{23}$ be the (elementary) matrices
  \begin{verbatim}
  E12:= addrow( I3,1,2,-1/2):
  E13:= addrow( I3,1,3,3):
  E23:= addrow( I3,2,3,4):
  \end{verbatim}
  Show that $E_{23} E_{13} E_{12} A$ is in row echelon form.
\end{enumerate}

\setcounter{theorem}{49}
\begin{theorem}
  Let $A$ be an $n\times n$ matrix. Then the following are equivalent:
  \begin{enumerate}
  \item[]
    \begin{enumerate}
    \item $A$ is invertible;
    \item $AX=B$ has a unique solution;
    \item $AX=\b 0$ has the unique solution $X=\b 0$;
    \item $A$ is row-equivalent to the identity matrix;
    \item $A$ is a product of elementary matrices;
    \end{enumerate}
  \end{enumerate}
\end{theorem}
\bigskip

The fact that any matrix can be reduced to Echelon form by pre-multiplying it
by a sequence of (invertible) elementary matrices suggests a method for
inverting a matrix by hand illustrated by Hill \S1.4 Examples 9 and 10 on
pages 46-48 (although Maple has a built-in function to do this).
\smallskip

\begin{verbatim}
# Hill Section 1.4 Example 9, p46/47
> with(linalg):
# Use the function  augment  to construct the matrix  [A : I]
> A:= matrix([ [ 2, 1, 0],
>              [-4,-1,-3],
>              [ 3, 1, 2] ]):
> Id:= diag( 1,1,1):                       A_I:= augment(A,Id):
#
# Reduce this matrix to Echelon form
> AI1:= gausselim(A_I);
                              [ 2  1   0     1    0   0 ]
                              [                         ]
                       AI1 := [ 0  1   -3    2    1   0 ]
                              [                         ]
                              [ 0  0  1/2  -1/2  1/2  1 ]

# Do a formal back-substitution on this using Maple's row operations 
> AI2:= addrow(AI1,3,2,6):                        AI3:= addrow(AI2,2,1,-1);
                             [ 2  0   0     2    -4  -6 ]
                             [                          ]
                      AI3 := [ 0  1   0    -1    4    6 ]
                             [                          ]
                             [ 0  0  1/2  -1/2  1/2   1 ]
\end{verbatim}
\begin{verbatim}
# Now divide each row by the diagonal element.
> AI4:= mulrow(AI3,1,1/2):                         AI5:= mulrow(AI4,3,2);
                                [ 1  0  0   1  -2  -3 ]
                                [                     ]
                         AI5 := [ 0  1  0  -1   4   6 ]
                                [                     ]
                                [ 0  0  1  -1   1   2 ]

# The inverse of the matrix A is now the right hand half of this matrix:
> Inverse:= submatrix( AI5,1..3,4..6 ):        equal( evalm(Inverse&*A), Id );
                                      true

# Hill Section 1.4 Example 10, p48
> A_I[3,3]:= 3/2:     # Otherwise, the entries are the same as in Example 9
#
# Reduce the matrix to Echelon form as above
> Ech:= gausselim(A_I);
                              [ 2  1   0    1    0   0 ]
                              [                        ]
                       Ech := [ 0  1  -3    2    1   0 ]
                              [                        ]
                              [ 0  0   0  -1/2  1/2  1 ]

# As the left hand half of the last row of this matrix consists of
#   a row of zeroes, the matrix A has no inverse.
\end{verbatim}
\headq{EXERCISE} At each stage of Example 9 above, let\par
  
  \verb$ A.i:= submatrix(AI.i,1..3,1..3);$\par
  \verb$ E.i:= submatrix(AI.i,1..3,4..6);$

Verify that\quad $E_i A = A_i$,\quad and hence show that\quad
$E_5 E_4 E_3 E_2 E_1 = A^{-1}$.

\goodbreak
\subsection{Gaussian Elimination as a Matrix Factorization}
\headq{Definitions.} A square matrix is said to be
\begin{enumerate}
\item[] {\bf upper triangular} if all entries below the main diagonal are zero
  (i.e. row echelon form).
\item[] {\bf lower triangular} if all entries above the main diagonal are zero. 
\item[] {\bf unit (lower or upper) triangular} if it is (lower or upper)
  triangular matrix and its diagonal entries are all one.
\end{enumerate}
For example, all elementary matrices are lower triangular, and
the elementary matrices $E_{ij}$ of the previous section
(corresponding to the Maple command {\tt addrow}) are unit lower triangular.

\setcounter{theorem}{55}
\begin{theorem}
  {\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}

\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]
  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}

\label{lufac}
\begin{verbatim}
# Hill Section 1.5 Example 3, p59/60
> with(linalg):
# Find the LU decomposition of the matrix
> A:= matrix([ [ 2,  4, -4],
>              [ 1, -4,  3],
>              [-6, -9, 10] ]):
# Start with 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
> 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   4   -4 ]
                        [              ]  [            ]
                        [ 1/2    1   0 ], [ 0  -6   5  ]
                        [              ]  [            ]
                        [  -3  -1/2  1 ]  [ 0   0  1/2 ]

# Verify that    A = L*U
> equal( A,evalm(L&*U) );
                                      true

# Hill Section 1.5 Example 4, p61
# Solve the equation  A*X = B  (i.e. L*U*X = B)  where
> B:= vector( [12,-21,-24] ):
#
# Let  U*X = Y  (and so  L*Y = B)  where Y is found by forward substitution:
> y1:= B[1]:                                 y2:= B[2] - L[2,1]*y1:
> y3:= B[3] - L[3,1]*y1 - L[3,2]*y2:         Y:= vector( [y1,y2,y3] );
                             Y := [ 12, -27, -3/2 ]

# The solution is now given by
> X:= backsub( augment(U,Y) );
                               X := [ -4, 2, -3 ]
\end{verbatim}
\goodbreak

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

\headq {Definition.} An {\bf elementary permutation matrix} is a matrix
$P_{ij}$ corresponding to the maple command {\tt swaprow}, and a
{\bf 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} 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 only on $P$ and $A$).
\end{theorem}

\begin{verbatim}
# Hill Section 1.5 Example 6, pp 63-65
# An example of the  PA = LU  decomposition of a matrix
> with(linalg):
> A:= matrix([ [0,0, 3],
>              [2,4,-1],
>              [6,5, 5] ]):
# The following row operations will reduce  A  to Echelon form:
> U:= swaprow(A,1,2):
> U:= addrow(U,1,3,-U[3,1]/U[1,1]):
> U:= swaprow(U,2,3):         # (The value of this matrix U is given later)
#
# Alternatively, perform the above row interchanges on the identity matrix
> P:= diag(1,1,1):           P:= swaprow(P,1,2):        P:= swaprow(P,2,3);
                                     [ 0  1  0 ]
                                     [         ]
                                P := [ 0  0  1 ]
                                     [         ]
                                     [ 1  0  0 ]

> PA:= evalm( P&*A );
                                     [ 2  4  -1 ]
                                     [          ]
                               PA := [ 6  5   5 ]
                                     [          ]
                                     [ 0  0   3 ]

# Now perform an LU decomposition on PA, as in Example 3 above
> L:= diag(1,1,1):           L[2,1]:= PA[2,1]/PA[1,1]:
> print( L,U );
                           [ 1  0  0 ]  [ 2   4  -1 ]
                           [         ]  [           ]
                           [ 3  1  0 ], [ 0  -7   8 ]
                           [         ]  [           ]
                           [ 0  0  1 ]  [ 0   0   3 ]

> equal( PA,evalm(L&*U) );
                                      true

# Hill Section 1.5 Example 7, p 65
# The above can be used (as in example 4) to solve the equation  A*X = B where
> B:= vector( [6,-10,-7] ):
# by solving  P*A*X = P*B  where
> PB:= evalm( P&*B );
                              PB := [ -10, -7, 6 ]
\end{verbatim}
\goodbreak
\headq{EXERCISES}
\nopagebreak
\begin{enumerate}
\item Complete the above example 7.
\item If you have some experience with Maple procedures, try writing
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. For details, see Hill's Appendices
A and B.
\bigskip
\goodbreak

\subsection{Transposes, Symmetry, and Band Matrices}
The importance of symmetric matrices will be seen when you reach Hill \S 5.4
and \S 5.7.  For the present, it is 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{verbatim}
# Hill Section 1.6 Theorem 1.69, p70.
> A:= array( symmetric,1..3,1..3 ):
> with(linalg):
> equal( A,transpose(A) );
                                      true
\end{verbatim}

To avoid conflict with Maple's {\tt band} function, we will use the word
{\bf banded} where Hill uses {\bf band}.  Banded matrices have a nice
factorization as can be seen in the following example.

\begin{verbatim}
# Hill Section 1.6 Example 5, p73.
> with(linalg):
# The LU decomposition of  the matrix
> A:= matrix([ [ 1, 1,  0, 0],
>              [ 2, 4,  1, 0],
>              [ 0, 6,  7, 1],
>              [ 0, 0, 12, 2] ]):
# is given by
> L:= diag(1,1,1,1):                               U:= A:
> L[2,1]:= U[2,1]/U[1,1]:               U:= addrow( U,1,2,-L[2,1]):
> L[3,2]:= U[3,2]/U[2,2]:               U:= addrow( U,2,3,-L[3,2]):
> L[4,3]:= U[4,3]/U[3,3]:               U:= addrow( U,3,4,-L[4,3]):
#
> print( L,U );          # these are both banded matrices.
                        [ 1  0  0  0 ]  [ 1  1  0   0 ]
                        [            ]  [             ]
                        [ 2  1  0  0 ]  [ 0  2  1   0 ]
                        [            ], [             ]
                        [ 0  3  1  0 ]  [ 0  0  4   1 ]
                        [            ]  [             ]
                        [ 0  0  3  1 ]  [ 0  0  0  -1 ]

> inverse(A);             # this is NOT a banded matrix.
                           [ 1/2   1/4  -1/4   1/8 ]
                           [                       ]
                           [ 1/2  -1/4   1/4  -1/8 ]
                           [                       ]
                           [  -3   3/2  -1/2   1/4 ]
                           [                       ]
                           [  18   -9     3    -1  ]
\end{verbatim}

\medskip
Your texts and lectures will 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 (See the comments on page 95 of Hill).

This is especially obvious in the case of banded matrices such as above
where only 3 steps were required in finding the LU decomposition of the
matrix $A$, whereas a general $4\times 4$ matrix may involve
finding\quad $L[2,1], L[3,1], L[4,1], L[3,2], L[4,2], L[4,3]$\quad and so
require 6 steps. In general, the LU decomposition of a banded $n\times n$
matrix with bands of length $\ell$ requires $(\ell-2)(n-1)$ steps, whereas
that of an arbitrary $n\times n$ matrix may require\quad
$(n-1)+(n-2)+\cdots +1 = n(n-1)/2$\quad steps.
Also the number of steps for solving the equations $L\b y=\b b$ and
$L^T\b x=\b y$ is each at most $n$. However, the inverse matrix is generally
not a banded matrix, and forming $A^{-1}\b b$ requires at least $n^2$ steps.
This is clearly a great saving for large matrices as Hill's truss example shows.

\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
\subsection{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 are very important.
For more details of Numerical Computing packages, see Appendix D of Hill.

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}\index{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, and is discussed in Section 5.8 of Hill.

\subsubsection*{Partial Pivotting}
Roundoff error is the error associated with the storage or arithmetic of
floating-point numbers on computers with a finite precision.
An example is given in the following.

\begin{verbatim}
# Hill Section 1.7 Example 1, p75
# We set the value of Digits to 5, and try to solve the equations
> eqn1:= 10^(-8)*x+y=1;                                   eqn2:= x+y=2;
                         eqn1 := 1/100000000 x + y = 1

                               eqn2 := x + y = 2

> vars:= {x,y}:            Digits:= 5:            fsolve( {eqn1,eqn2},vars );
                            {y = 1.0000, x = 1.0000}

# If we pivot on the first coefficient of  eqn2,  we get
> eqn:= evalf(eqn1-eqn2/10^8);                      fsolve( {eqn,eqn2},vars );
                           eqn := 1.00000 y = 1.00000

                            {y = 1.0000, x = 1.0000}

# which is the correct answer to five significant figures.
#
# However, if we try to pivot on the first coefficient of  eqn1,  we get
> eqn:= evalf(eqn2-10^8*eqn1);                     fsolve( {eqn1,eqn},vars );
                                         9                9
                     eqn :=  - .100000*10  y = -.100000*10

                              {x = 0, y = 1.0000}

# which is not the correct answer.
# Note that Maple gives the correct answer (in both cases) with the value
> Digits:= 10:
\end{verbatim}
\headq{Strategy 1.73 (Partial Pivotting)} {\em When doing Gaussian
elimination, a computer using floating-point arithmetic should choose,
as the pivot, the largest (in absolute value) matrix entry in the column
being considered. This is done by interchanging the corresponding rows.}

\subsubsection*{Overwriting Matrices}
Situations often arise in which we do not need to store all the entries of a
matrix. One example is in the LU decomposition of a matrix as in Example~3 of
Hill Section 1.5 (see page~\pageref{lufac} of these notes). When the matrix
$(L,U)$ is stored, the zeroes in $L$ may be over-written by the entries of $U$.

In order to save space under such circumstances, Maple allows an argument
{\tt sparse}\index{sparse} when defining a matrix using the function {\tt array}
(but {\em not} when using the function {\tt matrix}). This causes only those
entries actually assigned to be stored.
  
Another example of this economy of space is given in the symmetric matrices
mentioned in Hill section~1.6. In this case, only just over half of
the entries need to be stored, and so the Maple command (with $n$ assigned to a
positive integer)
\begin{verbatim}
A:= array(symmetric,1..n,1..n);
\end{verbatim} 
stores the $n(n+1)/2$ entries on or above the diagonal, together with the rule
\quad {\tt A[j,i]:= A[i,j]}.

\subsubsection*{Ill-conditioned Systems}
Example 1 above, on p75 of Hill, shows how a problem of round-off can be
dealt with using partial pivotting. However there are some systems of
linear equations where this cannot be avoided. In the next example, note how a
small change in the second equation causes a large change in the answer.

\begin{verbatim}
# Example 6 of Hill Section 1.7, p79:      Solve
> eqn1:= x + 2*y = 3;                             eqn2:= 1.001*x + 2*y = 3+f;
                              eqn1 := x + 2 y = 3

                         eqn2 := 1.001 x + 2 y = 3 + f

> vars:= {x,y}:                  Digits:= 5:
> for f from 1/1000 by 1/1000 to 3/1000 do
>     print( evalf(eqn2),solve( {eqn1,eqn2},vars) );
> od:
               1.001 x + 2. y = 3.0010, {x = 1.0000, y = 1.0000}

                1.001 x + 2. y = 3.0020, {x = 2.0000, y = .5000}

                  1.001 x + 2. y = 3.0030, {y = 0, x = 3.0000}
\end{verbatim} 
  
\subsubsection*{Condition Number}
This is a measure of how ill-conditioned a system is.
Maple can also calculate matrix condition numbers,\quad {\tt cond}\index{cond}.
\quad For example, using the matrix of the above equations,
\begin{verbatim}
> with(linalg):
> A:= matrix([ [ 1,2], [ 1.001,2] ]):
> cond(A);
                                  6002.000000

# Note that Maple's function  cond  may have a second parameter
# which yields different answers.
\end{verbatim} 
For more details on the different definitions of the condition number, see
Section 6.4 of Hill.

\vfill
\eject

