15_LA_gaussian
Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli
Note: This material largely follows the text “Numerical Linear Algebra” by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there.
In [ ]:
%matplotlib notebook
%precision 16
import numpy
import matplotlib.pyplot as plt
Gaussian Elimination¶
Gaussian elimination is the process where one transforms a matrix or linear system through a series of operations into one that is at the very least upper triangular (it is often also presented as including the operation that transforms $A$ completely to diagonal). These series of operations can be written as a sequence of successive matrix-matrix multiplications by lower triangular matrices. Letting these matrices be $L_j$ and the resulting upper triangular matrix be $U$ we can write this as
$$
\overbrace{L_{m-1} \cdots L_2 L_1}^{L^{-1}} A = U.
$$
Labeling the successive operations as $L^{-1}$ allows us to move $L$ to the other side of the equation and we see that in fact what we have done is computed another factorization of the matrix $A$ called the $LU$ factorization.
Example¶
As an example of this process lets consider the matrix
$$
A = \begin{bmatrix}
2 & 1 & 1 & 0 \\
4 & 3 & 3 & 1 \\
8 & 7 & 9 & 5 \\
6 & 7 & 9 & 8
\end{bmatrix}
$$
The first step is to remove the values in the first column below the diagonal, to do this we can multiply by the matrix
$$
L_1 = \begin{bmatrix}
1 & & & \\
-2 & 1 & & \\
-4 & & 1 & \\
-3 & & & 1
\end{bmatrix} \text{ so that } L_1 A = \begin{bmatrix}
2 & 1 & 1 & 0 \\
& 1 & 1 & 1 \\
& 3 & 5 & 5 \\
& 4 & 6 & 8
\end{bmatrix}.
$$
The next step is to remove the values below the diagonal of the second column. This can be done with
$$
L_2 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& -3 & 1 & \\
& -4 & & 1
\end{bmatrix} \text{ so that } L_2 L_1 A = \begin{bmatrix}
2 & 1 & 1 & 0 \\
& 1 & 1 & 1 \\
& & 2 & 2 \\
& & 2 & 4
\end{bmatrix}.
$$
Finally we multiply $A$ by $L_3$ defined as
$$
L_3 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& & 1 & \\
& & -1 & 1
\end{bmatrix} \text{ so that } L_3 L_2 L_1 A = \begin{bmatrix}
2 & 1 & 1 & 0 \\
& 1 & 1 & 1 \\
& & 2 & 2 \\
& & & 2
\end{bmatrix}
$$
completing the factorization with
$$
L^{-1} = L_3 L_2 L_1 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& & 1 & \\
& & -1 & 1
\end{bmatrix}
\begin{bmatrix}
1 & & & \\
& 1 & & \\
& -3 & 1 & \\
& -4 & & 1
\end{bmatrix}
\begin{bmatrix}
1 & & & \\
-2 & 1 & & \\
-4 & & 1 & \\
-3 & & & 1
\end{bmatrix}
$$
and
$$
U = \begin{bmatrix}
2 & 1 & 1 & 0 \\
& 1 & 1 & 1 \\
& & 2 & 2 \\
& & & 2
\end{bmatrix}
$$
We can actually easily invert $L$ and finally can write $A$ as
$$
A =
\begin{bmatrix}
1 & & & \\
2 & 1 & & \\
4 & 3 & 1 & \\
3 & 4 & 1 & 1
\end{bmatrix}\begin{bmatrix}
2 & 1 & 1 & 0 \\
& 1 & 1 & 1 \\
& & 2 & 2 \\
& & & 2
\end{bmatrix}
$$
Pivoting¶
As you may recall from your linear algebra course pivoting of the rows and columns of $A$ is often an important addition to Gaussian elimination to ensure that we can in fact factorize the matrix. As a simple example take the matrix
$$
A = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}.
$$
Without switch the rows Gaussian elimination would fail at the first step! This is made more hazardous if we consider the matrix
$$
A = \begin{bmatrix} 10^{-17} & 1 \\ 1 & 1 \end{bmatrix}.
$$
on a finite precision machine. In principle any row and column can be pivoted so that at each step we have the maximum value being used (on the diagonal) to perform the operations that compose the matrices $L_j$. In practice however we restrict ourselves to only pivoting rows of the matrix (called partial pivoting).
Consider again the example from above and switch the 1st and 3rd rows using the criteria that we always want to use the largest value to do perform the reduction. Defining the first pivot matrix as
$$
P_1 = \begin{bmatrix}
& & 1 & \\
& 1 & & \\
1 & & & \\
& & & 1
\end{bmatrix} \quad \text{so that} \quad
P_1 A = \begin{bmatrix}
8 & 7 & 9 & 5 \\
4 & 3 & 3 & 1 \\
2 & 1 & 1 & 0 \\
6 & 7 & 9 & 8
\end{bmatrix}
$$
Now defining the first $L_1$ as
$$
L_1 = \begin{bmatrix}
1 & & & \\
-\frac{1}{2} & 1 & & \\
-\frac{1}{4} & & 1 & \\
-\frac{3}{4} & & & 1
\end{bmatrix} \quad \text{so that} \quad
L_1 P_1 A = \begin{bmatrix}
8 & 7 & 9 & 5 \\
& -\frac{1}{2} & -\frac{3}{2} & -\frac{3}{2} \\
& -\frac{3}{4} & -\frac{5}{4} & -\frac{5}{4} \\
& \frac{7}{4} & \frac{9}{4} & \frac{17}{4}
\end{bmatrix}.
$$
Again examining the remaining values in column 2 the maximum value lies in row 4 so we want to interchange this with the second row (note that we do not want to move the first row as that will bring non-zero values into the first column below the diagonal).
$$
P_2 = \begin{bmatrix}
1 & & & \\
& & & 1 \\
& & 1 & \\
& 1 & &
\end{bmatrix} \quad \text{and} \quad
L_2 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& \frac{3}{7} & 1 & \\
& \frac{2}{7} & & 1
\end{bmatrix} \quad \text{so that} \quad
L_2 P_2 L_1 P_1 A = \begin{bmatrix}
8 & 7 & 9 & 5 \\
& \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\
& & -\frac{2}{7} & \frac{4}{7} \\
& & -\frac{6}{7} & -\frac{2}{7}
\end{bmatrix}.
$$
Finally
$$
P_3 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& & & 1 \\
& & 1 &
\end{bmatrix} \quad \text{and} \quad
L_3 = \begin{bmatrix}
1 & & & \\
& 1 & & \\
& & 1 & \\
& & -\frac{1}{3} & 1
\end{bmatrix} \quad \text{so that} \quad
L_3 P_3 L_2 P_2 L_1 P_1 A = \begin{bmatrix}
8 & 7 & 9 & 5 \\
& \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\
& & -\frac{6}{7} & -\frac{2}{7} \\
& & & \frac{2}{3} \\
\end{bmatrix}.
$$
LU Factorization with Partial Pivoting¶
Due to the nature of the pivot matrices we can disentangle them from the matrices $L_j$. Right now we have
$$
L_3 P_3 L_2 P_2 L_1 P_1 A = U
$$
where what we want is
$$
(L’_3 L’_2 L’_1) P_3 P_2 P_1 A = U.
$$
It turns out we can easily do this by
$$
L’_3 = L_3, \quad L’_2 = P_3 L_2 P_3^{-1}, \quad \text{and} \quad L’_1 = P_3 P_2 L_1 P_2^{-1} P_3^{-1}.
$$
These new matrices $L’_j$ can easily be computed and it turns out that the $L’_j$ matrices have the same structure with the rows permuted accordingly.
In general then the $LU$ factorization with partial pivoting of the above example can be written as
$$
\underbrace{\begin{bmatrix}
& & 1 & \\
& & & 1 \\
& 1 & & \\
1 & & &
\end{bmatrix}}_{P = P_3 P_2 P_1} \underbrace{\begin{bmatrix}
2 & 1 & 1 & 0 \\
4 & 3 & 3 & 1 \\
8 & 7 & 9 & 5 \\
6 & 7 & 9 & 8
\end{bmatrix}}_{A} =
\underbrace{\begin{bmatrix}
1 & & & \\
\frac{3}{4} & 1 & & \\
\frac{1}{2} & -\frac{2}{7} & 1 & \\
\frac{1}{4} & -\frac{3}{7} & \frac{1}{3} & 1 \\
\end{bmatrix}}_{L}
\underbrace{\begin{bmatrix}
8 & 7 & 9 & 5 \\
& \frac{7}{4} & \frac{9}{4} & \frac{17}{4} \\
& & -\frac{6}{7} & -\frac{2}{7} \\
& & & \frac{2}{3} \\
\end{bmatrix}}_{U}
$$
The algorithm for the general factorization with partial pivoting then can be written as
U = A
L = I
P = I
for k = 1 to m – 1
Select i >= k to maximize |u_{i,k}|
U[k, k:m] <==> U[i, k:m]
L[k, 1:k-1] <==> L[i, 1:k-1]
P[k, :] <==> P[i, :]
for j = k + 1 to m
L[j, k] = U[j, k] / U[k, k]
U[j, k:m] = U[j, k:m] – L[j, k] * U[k, k:m]
where <==> represents swapping of the two rows indicated.
Solving $Ax = b$¶
To complete our discussion lets consider the solution of the linear system $Ax = b$. We now have a factorization of the matrix $PA = LU$ so that the new system is $LU x = b$ (note that pivoting is allowed as long as we also interchange the elements of $b$ via $P \cdot b$). We can do this in two steps:
Compute the pivoted $LU$ factorization $PA = LU$.
Compute the solution to $L y = P b$ using forward substitution.
Compute the solution to $U x = y$ using backward substitution.
Forward Substitution¶
For forward substitution we proceed from the first row and progress downwards through the matrix. We can then consider the general $i$th row with
$$
L_{i,1} y_1 + L_{i,2} y_2 + \cdots + L_{i,i-1} y_{i-1} + y_i = b_i
$$
noting that we are using the fact that the matrix $L$ has 1 on its diagonal. We can now solve for $y_i$ as
$$
y_i = b_i – \left( L_{i,1} y_1 + L_{i,2} y_2 + \cdots + L_{i,i-1} y_{i-1} \right ).
$$
Backward Substitution¶
Backwards substitution requires us to move from the last row of $U$ and move upwards. We can consider again the general $i$th row with
$$
U_{i,i} x_i + U_{i,i+1} x_{i+1} + \ldots + U_{i,m-1} x_{m-1} + U_{i,m} x_m = y_i.
$$
We can now solve for $x_i$ as
$$
x_i = \frac{1}{U_{i,i}} \left( y_i – ( U_{i,i+1} x_{i+1} + \ldots + U_{i,m-1} x_{m-1} + U_{i,m} x_m) \right )
$$