luAlgorithm.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{luAlgorithm}
%\renewcommand{\dirname}{notes/FIXMEwheretodirname/}
%%\newcommand{\dateintitle}{}
%%\newcommand{\keywords}{}
%
%\input{../peeter_prologue_print2.tex}
%
%\beginArtNoToc
%
%\generatetitle{Illustrating the LU algorithm with pivots by example}
%%\chapter{Illustrating the LU algorthm with pivots by example}
%%\label{chap:luAlgorithm}
%
Two previous examples of LU factorizations were given.  I found one more to be the key to understanding how to implement this as a Matlab algorithm, required for \cref{multiphysics:problemSet1:2}.
%
A matrix that contains both pivots and elementary matrix operations is
%
\begin{equation}\label{eqn:luAlgorithm:20}
M=
\begin{bmatrix}
0 & 0 & 2 & 1 \\
0 & 0 & 1 & 1 \\
2 & 0 & 2 & 0 \\
1 & 1 & 1 & 1
\end{bmatrix}
\end{equation}
%
The objective is to apply a sequence of row permutations or elementary row operations to \( M \) that put \( M \) into upper triangular form, while also tracking all the inverse operations.  When no permutations were required to produce \( U \), then a factorization \( M = L' U \) is produced where \( L' \) is lower triangular.
%
The row operations to be applied to \( M \) are
%
\begin{equation}\label{eqn:luAlgorithm:40}
U =
L_k^{-1}
L_{k-1}^{-1} \cdots
L_2^{-1}
L_1^{-1}
M,
\end{equation}
%
with
%
\begin{equation}\label{eqn:luAlgorithm:60}
L' = L_0 L_1 L_2 \cdots L_{k-1} L_k
\end{equation}
%
Here \( L_0 = I \), the identity matrix, and \( L_i^{-1} \) is either a permutation matrix interchanging two rows of the identity matrix, or it is an elementary row operation encoding the operation \( r_j \rightarrow r_j - M r_i \), where \( r_i \) is the pivot row, and \( r_j, j > i \) are the rows that the Gaussian elimination operations are applied to.
%
For this example matrix, the \( M_{1 1} \) value cannot be used as the pivot element since it is zero.  In general,
the row with the biggest absolute value in the column should be used.
In this case that is row 3.  The first row operation is therefore a \( 1,3 \) permutation.
For numeric stability, use
%
\begin{equation}\label{eqn:luAlgorithm:80}
L_1^{-1} =
\begin{bmatrix}
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix},
\end{equation}
%
which gives
%
\begin{equation}\label{eqn:luAlgorithm:100}
M \rightarrow
L_1^{-1}
M
=
\begin{bmatrix}
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
0 & 0 & 2 & 1 \\
0 & 0 & 1 & 1 \\
2 & 0 & 2 & 0 \\
1 & 1 & 1 & 1 \\
\end{bmatrix}
=
\begin{bmatrix}
2 & 0 & 2 & 0 \\
0 & 0 & 1 & 1 \\
0 & 0 & 2 & 1 \\
1 & 1 & 1 & 1 \\
\end{bmatrix}.
\end{equation}
%
Computationally, avoiding the matrix multiplication that achieve this permutation is desired.  Instead just swap the two rows in question.
%
The inverse of this operation is the same permutation, so for \( L' \) the first stage computation is
%
\begin{equation}\label{eqn:luAlgorithm:120}
L \sim L_0 L_1 = L_1.
\end{equation}
%
As before, a matrix operation would be very expensive.  When the application of the permutation matrix is from the right, it results in an exchange of columns \(1,3\) of the \( L_0 \) matrix (which happens to be identity at this point).  So the matrix operation can be done as a column exchange directly using submatrix notation.
%
Now proceed down the column, doing all the non-zero row elimination operations required.  In this case, there is only one operation todo
%
\begin{equation}\label{eqn:luAlgorithm:140}
r_4 \rightarrow r_4 - \frac{1}{2} r_1.
\end{equation}
%
This has the matrix form
%
\begin{equation}\label{eqn:luAlgorithm:160}
L_2^{-1} =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
-1/2 & 0 & 0 & 1 \\
\end{bmatrix}.
\end{equation}
%
The next stage of the \( U \) computation is
%
\begin{equation}\label{eqn:luAlgorithm:180}
\begin{aligned}
M
&\rightarrow L_2^{-1} L_1^{-1} M \\
&=
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
-1/2 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
2 & 0 & 2 & 0 \\
0 & 0 & 1 & 1 \\
0 & 0 & 2 & 1 \\
1 & 1 & 1 & 1 \\
\end{bmatrix} \\
&=
\begin{bmatrix}
2 & 0 & 2 & 0 \\
0 & 0 & 1 & 1 \\
0 & 0 & 2 & 1 \\
0 & 1 & 0 & 1 \\
\end{bmatrix}.
\end{aligned}
\end{equation}
%
Again, this is not an operation that should be done as a matrix operation.  Instead act directly on the rows in question with \cref{eqn:luAlgorithm:140}.
%
Note that the inverse of this matrix operation is very simple.
An amount \( r_1/2 \) has been subtracted from \( r_4 \), so to invert this all that is required is adding back \( r_1/2 \).  That is
%
\begin{equation}\label{eqn:luAlgorithm:200}
L_2
=
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
1/2 & 0 & 0 & 1 \\
\end{bmatrix}.
\end{equation}
%
Observe that when this is applied from the right to \( L_0 L_1 \rightarrow L_0 L_1 L_2\), the interpretation is a column operation
%
\begin{equation}\label{eqn:luAlgorithm:220}
c_1 \rightarrow c_1 + m c_4,
\end{equation}
%
In general, if application of the row operation
%
\begin{equation}\label{eqn:luAlgorithm:240}
r_j \rightarrow r_j - m r_i,
\end{equation}
%
to the current state of the matrix \( U \), requires application of the operation
%
\begin{equation}\label{eqn:luAlgorithm:260}
r_i \rightarrow r_i + m r_j,
\end{equation}
%
to the current state of the matrix \( L' \).
%
The next step is to move on to reduction of column 2, and for that only a permutation operation is required
%
\begin{equation}\label{eqn:luAlgorithm:280}
L_3 =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
\end{bmatrix},
\end{equation}
%
Application of a \( 2,4 \) row interchange to U, and a \( 2,4 \) column interchange to \( L' \) gives
%
\begin{equation}\label{eqn:luAlgorithm:300}
M \rightarrow
\begin{bmatrix}
2 & 0 & 2 & 0 \\
0 & 1 & 0 & 1 \\
0 & 0 & 2 & 1 \\
0 & 0 & 1 & 1 \\
\end{bmatrix}.
\end{equation}
%
The final operation is a regular row operation
%
\begin{equation}\label{eqn:luAlgorithm:320}
r_4 \rightarrow r_4 - \inv{2} r_3,
\end{equation}
%
with matrix
\begin{equation}\label{eqn:luAlgorithm:340}
L_4^{-1} =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & -1/2 & 1 \\
\end{bmatrix}
\end{equation}
%
The composite permutation performed so far is
%
\begin{equation}\label{eqn:luAlgorithm:360}
P = L_3 L_1 I.
\end{equation}
%
This should also be computed by performing row interchanges, not matrix multiplication.
%
To solve the system
%
\begin{equation}\label{eqn:luAlgorithm:380}
M \Bx = L' U \Bx = \Bb,
\end{equation}
%
solve the equivalent problem
%
\begin{equation}\label{eqn:luAlgorithm:400}
P L' U \Bx = P \Bb,
\end{equation}
%
To do this let \( \By = U \Bx \), for
%
\begin{equation}\label{eqn:luAlgorithm:420}
P L' \By = P \Bb.
\end{equation}
%
The matrix \( L = P L' \) is lower triangular, as \( P \) contained all the permutations that applied along the way (FIXME: this is a statement, not a proof, and not obvious).
%
The system
%
\begin{equation}\label{eqn:luAlgorithm:440}
L \By = P \Bb,
\end{equation}
%
can now be solved using forward substitution, after which the upper triangular system
%
\begin{equation}\label{eqn:luAlgorithm:460}
\By = U \Bx,
\end{equation}
%
can be solved using only back substitution.
%
%\EndArticle
%\EndNoBibArticle