Subsection 5.2.1 Gaussian elimination
¶ VIDEO
Homework 5.2.1.1 .
Solve
\begin{equation*}
\left(\begin{array}{rrr}
2 \amp -1 \amp 1\\
-4 \amp 0 \amp 1\\
4 \amp 0 \amp -2
\end{array}\right)
\left( \begin{array}{c}
\chi_0 \\
\chi_1 \\
\chi_2
\end{array} \right)
=
\left( \begin{array}{r}
-6 \\
2 \\
0
\end{array} \right) .
\end{equation*}
Answer
\begin{equation*}
\left( \begin{array}{c}
\chi_0 \\
\chi_1 \\
\chi_2
\end{array} \right)
=
\left( \begin{array}{r}
-1 \\
2 \\
-2
\end{array} \right) .
\end{equation*}
Solution
We employ Gaussian elimination applied to an appended system:
\begin{equation*}
\left(\begin{array}{rrr| r}
2 \amp -1 \amp 1 \amp -6\\
-4 \amp 0 \amp 1 \amp 2\\
4 \amp 0 \amp -2 \amp 0
\end{array}\right)
\end{equation*}
Compute the multiplier \(\lambda_{10} = (-4)/(2) =
-2 \)
Subtract \(\lambda_{10} = -2 \) times the first row from the second row, yielding
\begin{equation*}
\left(\begin{array}{rrr| r}
2 \amp -1 \amp 1 \amp -6\\
0 \amp -2 \amp 3 \amp -10\\
4 \amp 0 \amp -2 \amp 0
\end{array}\right)
\end{equation*}
Compute the multiplier \(\lambda_{20} = (4)/(2) =
2 \)
Subtract \(\lambda_{20} = 2 \) times the first row from the third row, yielding
\begin{equation*}
\left(\begin{array}{rrr| r}
2 \amp -1 \amp 1 \amp -6\\
0 \amp -2 \amp 3 \amp -10\\
0 \amp 2 \amp -4 \amp 12
\end{array}\right)
\end{equation*}
Compute the multiplier \(\lambda_{21} = (2)/(-2) =
-1 \)
Subtract \(\lambda_{21} = -1 \) times the second row from the third row, yielding
\begin{equation*}
\left(\begin{array}{rrr| r}
2 \amp -1 \amp 1 \amp -6\\
0 \amp -2 \amp 3 \amp -10\\
0 \amp 0 \amp -1 \amp 2
\end{array}\right)
\end{equation*}
Solve the triangular system
\begin{equation*}
\left(\begin{array}{rrr}
2 \amp -1 \amp 1 \\
0 \amp -2 \amp 3 \\
0 \amp 0 \amp -1
\end{array}\right)
\left(\begin{array}{c}
\chi_0 \\
\chi_1 \\
\chi_2
\end{array}
\right)
=
\left(\begin{array}{r}
-6 \\
-10\\
2
\end{array}\right)
\end{equation*}
to yield
\begin{equation*}
\left( \begin{array}{c}
\chi_0 \\
\chi_1 \\
\chi_2
\end{array} \right)
=
\left( \begin{array}{r}
-1 \\
2 \\
-2
\end{array} \right) .
\end{equation*}
The exercise in Homework 5.2.1.1 motivates the following algorithm, which reduces the linear system \(A x = b \) stored in \(n \times n \) matrix \(A \) and right-hand side vector \(b \) of size \(n \) to an upper triangular system.
\begin{equation*}
\begin{array}{l}
{\bf for~} j:=0, \ldots , n-1 \\
~~~ {\bf for~} i:=j+1, \ldots , n-1 \\
~~~ ~~~ \lambda_{i,j} := \alpha_{i,j} / \alpha_{j,j} \\
~~~ ~~~ \alpha_{i,j} := 0 \\
~~~ ~~~
\left.
\begin{array}{l}
{\bf for~} k=j+1, \ldots , n-1 \\
~~~ \alpha_{i,k} := \alpha_{i,k} - \lambda_{i,j}
\alpha_{j,k} \\
{\bf endfor} \\
\beta_i := \beta_i - \lambda_{i,j} \beta_j
\end{array}
\right\} \mbox{ subtract } \lambda_{i,j} \mbox{ times row }
j \mbox{ from row } k
\\
~~~ {\bf endfor} \\
{\bf endfor}
\end{array}
\end{equation*}
This algorithm completes as long as no divide by zero is encountered.
Let us manipulate this a bit. First, we notice that we can first reduce the matrix to an upper triangular matrix, and then update the right-hand side using the multipliers that were computed along the way (if these are stored):
\begin{equation*}
\begin{array}{l}
\mbox{reduce } A \mbox{ to upper triangular form} \\
{\bf for~} j:=0, \ldots , n-1 \\
~~~ {\bf for~} i:=j+1, \ldots , n-1 \\
~~~ ~~~ \lambda_{i,j} := \alpha_{i,j} / \alpha_{j,j} \\
~~~ ~~~ \alpha_{i,j} := 0 \\
~~~ ~~~
\left. \begin{array}{l}
{\bf for~} k=j+1, \ldots , n-1 \\
~~~ \alpha_{i,k} := \alpha_{i,k} - \lambda_{i,j}
\alpha_{j,k} \\
{\bf endfor}
\end{array}
\right\} \mbox{ subtract } \lambda_{i,j} \mbox{ times row }
j \mbox{ from row } k
\\
~~~ {\bf endfor} \\
{\bf endfor} \\
\\
\mbox{update } b \mbox{ using multipliers (forward substitution)} \\
{\bf for~} j:=0, \ldots , n-1 \\
~~~ {\bf for~} i:=j+1, \ldots , n-1 \\
~~~ ~~~ \beta_i := \beta_i - \lambda_{i,j} \beta_j \\
~~~ {\bf endfor} \\
{\bf endfor}
\end{array}
\end{equation*}
Ignoring the updating of the right-hand side (a process known as forward substitution), for each iteration we can first compute the multipliers and then update the matrix:
\begin{equation*}
\begin{array}{l}
{\bf for~} j:=0, \ldots , n-1 \\
~~~
\left.
\begin{array}{l}
{\bf for~} i:=j+1, \ldots , n-1 \\
~~~ \lambda_{i,j} := \alpha_{i,j} / \alpha_{j,j} \\
~~~ \alpha_{i,j} := 0 \\
{\bf endfor}
\end{array}
\right\} \mbox{ compute multipliers} \\
~~~ {\bf for~} i:=j+1, \ldots , n-1 \\
~~~ ~~~
\left.
\begin{array}{l}
{\bf for~} k=j+1, \ldots , n-1 \\
~~~ \alpha_{i,k} := \alpha_{i,k} - \lambda_{i,j}
\alpha_{j,k} \\
{\bf endfor}
\end{array}
\right\} \mbox{ subtract } \lambda_{i,j} \mbox{ times row }
j \mbox{ from row } k
\\
~~~ {\bf endfor} \\
{\bf endfor}
\end{array}
\end{equation*}
Since we know that \(\alpha_{i,j} \) is set to zero, we can use its location to store the multiplier:
\begin{equation*}
\begin{array}{l}
{\bf for~} j:=0, \ldots , n-1 \\
~~~
\left.
\begin{array}{l}
{\bf for~} i:=j+1, \ldots , n-1 \\
~~~ \alpha_{i,j} := \lambda_{i,j} = \alpha_{i,j} / \alpha_{j,j} \\
{\bf endfor}
\end{array}
\right\} \mbox{ compute all multipliers} \\
~~~ {\bf for~} i:=j+1, \ldots , n-1 \\
~~~ ~~~
\left.
\begin{array}{l}
{\bf for~} k=j+1, \ldots , n-1 \\
~~~ \alpha_{i,k} := \alpha_{i,k} - \alpha_{i,j}
\alpha_{j,k} \\
{\bf endfor}
\end{array}
\right\} \mbox{ subtract } \lambda_{i,j} \mbox{ times row }
j \mbox{ from row } k
\\
~~~ {\bf endfor} \\
{\bf endfor}
\end{array}
\end{equation*}
Finally, we can cast the computation in terms of operations with vectors and submatrices:
\begin{equation*}
\begin{array}{l}
{\bf for~} j:=0, \ldots , n-1 \\
~~~
\left( \begin{array}{c}
\alpha_{j+1,j} \\
\vdots \\
\alpha_{n-1,j}
\end{array}
\right)
:=
\left( \begin{array}{c}
\lambda_{j+1,j} \\
\vdots \\
\lambda_{n-1,j}
\end{array}
\right) / \alpha_{j,j}
\\
~~~
\left( \begin{array}{c c c}
\alpha_{j+1,j+1} \amp \cdots \amp \alpha_{j+1,n-1} \\
\vdots \amp \amp \vdots \\
\alpha_{n-1,j+1} \amp \cdots \amp \alpha_{n-1,n-1} \\
\end{array}
\right)
:= \\
~~~ ~~~ ~~~
\left( \begin{array}{c c c}
\alpha_{j+1,j+1} \amp \cdots \amp \alpha_{j+1,n-1} \\
\vdots \amp \amp \vdots \\
\alpha_{n-1,j+1} \amp \cdots \amp \alpha_{n-1,n-1} \\
\end{array}
\right)
-
\left( \begin{array}{c}
\alpha_{j+1,j} \\
\vdots \\
\alpha_{n-1,j}
\end{array}
\right)
\left( \begin{array}{c c c}
\alpha_{j,j+1} \amp \cdots \amp \alpha_{j,n-1} \\
\end{array} \right) \\
{\bf endfor}
\end{array}
\end{equation*}
In Figure 5.2.1.1 this algorithm is presented with our FLAME notation.
\begin{equation*}
\newcommand{\routinename}{A = \mbox{GE}( A )}
\newcommand{\guard}{
n( A_{TL} ) \lt n( A )
}
\newcommand{\partitionings}{
A \rightarrow
\FlaTwoByTwo{A_{TL}}{A_{TR}}
{A_{BL}}{A_{BR}}
}
\newcommand{\partitionsizes}{
A_{TL} {\rm ~is~} 0 \times 0
}
\newcommand{\repartitionings}{
\FlaTwoByTwo{A_{TL}}{A_{TR}}
{A_{BL}}{A_{BR}}
\rightarrow
\FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}}
{a_{10}^T}{\alpha_{11}}{a_{12}^T}
{A_{20}}{a_{21}}{A_{22}}
}
\newcommand{\moveboundaries}{
\FlaTwoByTwo{A_{TL}}{A_{TR}}
{A_{BL}}{A_{BR}}
\leftarrow
\FlaThreeByThreeTL{A_{00}}{a_{01}}{A_{02}}
{a_{10}^T}{\alpha_{11}}{a_{12}^T}
{A_{20}}{a_{21}}{A_{22}}
}
\newcommand{\update}{
\begin{array}{l}
a_{21} := l_{21} = a_{21} / \alpha_{11} \\
A_{22} := A_{22} - a_{21} a_{12}^T
\end{array}
}
\FlaAlgorithm
\end{equation*}
Figure 5.2.1.1. Gaussian elimination algorithm that reduced a matrix \(A\) to upper triangular form, storing the multipliers below the diagonal.
Homework 5.2.1.2 .
Apply the algorithm Figure 5.2.1.1 to the matrix
\begin{equation*}
\left(\begin{array}{rrr}
2 \amp -1 \amp 1\\
-4 \amp 0 \amp 1\\
4 \amp 0 \amp -2
\end{array}\right)
\end{equation*}
and report the resulting matrix. Compare the contents of that matrix to the upper triangular matrix computed in the solution of Homework 5.2.1.1 .
Answer
\begin{equation*}
\left(\begin{array}{r r r}
2 \amp -1 \amp 1 \\
-2 \amp -2 \amp 3 \\
2 \amp -1 \amp -1
\end{array}\right)
\end{equation*}
Solution
Partition:
\begin{equation*}
\left(\begin{array}{|rrr} \hline
2 \amp -1 \amp 1\\
-4 \amp 0 \amp 1\\
4 \amp 0 \amp -2
\end{array}\right)
\end{equation*}
First iteration:
\(\alpha_{21} := \lambda_{21} = \alpha_{21} /
\alpha_{11} \text{:}\)
\begin{equation*}
\left(\begin{array}{rrr} \hline
2 \amp -1 \amp 1\\
-2 \amp 0 \amp 1\\
2 \amp 0 \amp -2
\end{array}\right)
\end{equation*}
\(A_{22} := A_{22} - a_{21} a_{12}^T \text{:}\)
\begin{equation*}
\left(\begin{array}{| r rr} \hline
2 \amp -1 \amp 1\\
-2 \amp -2 \amp 3\\
2 \amp 2 \amp -4
\end{array}\right)
\end{equation*}
State at bottom of iteration:
\begin{equation*}
\left(\begin{array}{r | r r}
2 \amp -1 \amp 1\\ \hline
-2 \amp -2 \amp 3\\
2 \amp 2 \amp -4
\end{array}\right)
\end{equation*}
Second iteration:
\(\alpha_{21} := \lambda_{21} = \alpha_{21} /
\alpha_{11} \text{:}\)
\begin{equation*}
\left(\begin{array}{r| rr}
2 \amp -1 \amp 1\\ \hline
-2 \amp -2 \amp 3\\
2 \amp -1 \amp -4
\end{array}\right)
\end{equation*}
\(A_{22} := A_{22} - a_{21} a_{12}^T \text{:}\)
\begin{equation*}
\left(\begin{array}{r | rr}
2 \amp -1 \amp 1\\ \hline
-2 \amp -2 \amp 3\\
2 \amp -1 \amp -1
\end{array}\right)
\end{equation*}
State at bottom of iteration:
\begin{equation*}
\left(\begin{array}{r r| r}
2 \amp -1 \amp 1\\
-2 \amp -2 \amp 3\\ \hline
2 \amp -1 \amp -1
\end{array}\right)
\end{equation*}
Third iteration:
\(\alpha_{21} := \lambda_{21} = \alpha_{21} /
\alpha_{11} \text{:}\)
\begin{equation*}
\left(\begin{array}{r r| r}
2 \amp -1 \amp 1\\
-2 \amp -2 \amp 3\\ \hline
2 \amp -1 \amp -1
\end{array}\right)
\end{equation*}
(computation with empty vector).
\(A_{22} := A_{22} - a_{21} a_{12}^T \text{:}\)
\begin{equation*}
\left(\begin{array}{rr|r}
2 \amp -1 \amp 1\\
-2 \amp -2 \amp 3\\ \hline
2 \amp -1 \amp -1
\end{array}\right)
\end{equation*}
(update of empty matrix)
State at bottom of iteration:
\begin{equation*}
\left(\begin{array}{r r r |}
2 \amp -1 \amp 1\\
-2 \amp -2 \amp 3\\
2 \amp -1 \amp -1 \\ \hline
\end{array}\right)
\end{equation*}
The upper triangular matrix computed in Homework 5.2.1.1 was
\begin{equation*}
\left(\begin{array}{rrr| r}
2 \amp -1 \amp 1\\
0 \amp -2 \amp 3\\
0 \amp 0 \amp -1 \\ \hline
\end{array}\right)
\end{equation*}
which can be found in the upper triangular part of the updated matrix \(A \text{.}\)
Homework 5.2.1.3 .
Applying Figure 5.2.1.1 to the matrix
\begin{equation*}
A =
\left(\begin{array}{rrr}
2 \amp -1 \amp 1\\
-4 \amp 0 \amp 1\\
4 \amp 0 \amp -2
\end{array}\right)
\end{equation*}
yielded
\begin{equation*}
\left(\begin{array}{r r r }
2 \amp -1 \amp 1 \\
-2 \amp -2 \amp 3 \\
2 \amp -1 \amp -1
\end{array}\right) .
\end{equation*}
This can be thought of as an array that stores the unit lower triangular matrix \(L \) below the diagonal (with implicit ones on its diagonal) and upper triangular matrix \(U \) on and above its diagonal:
\begin{equation*}
L =
\left(\begin{array}{r r r }
1 \amp 0 \amp 0 \\
-2 \amp 1 \amp 0 \\
2 \amp -1 \amp 1
\end{array}\right)
\quad \mbox{and} \quad
U =
\left(\begin{array}{r r r }
2 \amp -1 \amp 1 \\
0 \amp -2 \amp 3 \\
0 \amp 0 \amp -1
\end{array}\right)
\end{equation*}
Compute \(B = L U \) and compare it to \(A \text{.}\)
Answer Magic! \(B = A \text{!}\)
Solution
\begin{equation*}
B = L U =
\left(\begin{array}{r r r }
1 \amp 0 \amp 0 \\
-2 \amp 1 \amp 0 \\
2 \amp -1 \amp 1
\end{array}\right)
\left(\begin{array}{r r r }
2 \amp -1 \amp 1 \\
0 \amp -2 \amp 3 \\
0 \amp 0 \amp -1
\end{array}\right) =
\left(\begin{array}{rrr}
2 \amp -1 \amp 1\\
-4 \amp 0 \amp 1\\
4 \amp 0 \amp -2
\end{array}\right)
= A.
\end{equation*}
