## Subsection5.2.1Gaussian elimination

###### Homework5.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*}
\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} \alpha_{j+1,j} \\ \vdots \\ \alpha_{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.

###### Homework5.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.

\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{.}$

###### Homework5.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{.}$

Magic! $B = A \text{!}$