6.2. LU Decomposition with Partial Pivoting
The formula for the elements of \(L\) in LU decomposition
\[ \begin{align*}
\ell_{ij} &= \dfrac{1}{u_{jj}} \left(a_{ij} - \displaystyle \sum_{k=1}^{j-1} \ell_{ik}u_{kj}\right),
\end{align*} \]
requires division by the pivot \(u_{jj}\). If a pivot is zero, the LU decomposition cannot proceed. Even when the pivot is non-zero, a very small pivot can lead to large multipliers \(\ell_{ij}\), making the computation sensitive to rounding errors. This can produce inaccurate solutions even when the exact solution is perfectly well behaved.
For example, consider the linear system
\[\begin{split} \begin{align*}
\begin{pmatrix}
0.0001 & 1 \\
1 & 1
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &=
\begin{pmatrix} 1 \\ 2 \end{pmatrix},
\end{align*}\end{split}\]
which has the solution
\[ \begin{align*}
x_1 &= \frac{10000}{9999} \approx 1.0001, &
x_2 &= \frac{9998}{9999} \approx 0.9999.
\end{align*} \]
Computing the LU decomposition of the coefficient matrix gives
\[\begin{split} \begin{align*}
L &= \begin{pmatrix}
1 & 0 \\
10000 & 1
\end{pmatrix}, &
U &= \begin{pmatrix}
0.0001 & 1 \\
0 & -9999
\end{pmatrix},
\end{align*} \end{split}\]
To illustrate the effect of finite-precision arithmetic, suppose all calculations are rounded to three decimal places at each stage. Solving \(L\mathbf{y} = \mathbf{b}\) by forward substitution
\[\begin{split} \begin{align*}
\begin{pmatrix}
1 & 0 \\
10000 & 1
\end{pmatrix}
\begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &=
\begin{pmatrix} 1 \\ 2 \end{pmatrix}
\quad \implies \quad
\begin{aligned} y_1 &= 1 \\ y_2 &= 2 - 10000(1) = -9998, \end{aligned}
\end{align*} \end{split}\]
and solving \(U \mathbf{x} = \mathbf{y}\) by back substitution
\[\begin{split} \begin{align*}
\begin{pmatrix}
0.0001 & 1 \\
0 & -9999
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &=
\begin{pmatrix} 1 \\ -9998 \end{pmatrix}
\quad \implies \quad
\begin{aligned}
x_2 &= \frac{-9998}{-9999} \approx 1.000 \\
x_1 &\approx \frac{1}{0.0001}(1 - 1) = 0.
\end{aligned}
\end{align*} \end{split}\]
So we have a computed solution \(x_1 = 0\) whereas the exact solution is \(x_1 \approx 1\). In partial pivoting, the current pivot is chosen to be the entry with the largest absolute value in the current column among the rows that have not yet been eliminated. This is achieved by swapping rows before the elimination step.
\[\begin{split} \begin{align*}
\begin{pmatrix}
1 & 1 \\
0.0001 & 1
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &=
\begin{pmatrix} 2 \\ 1 \end{pmatrix}.
\end{align*}\end{split}\]
Computing the LU decomposition of the coefficient matrix gives
\[\begin{split} \begin{align*}
L &= \begin{pmatrix}
1 & 0 \\
0.0001 & 1
\end{pmatrix}, &
U &= \begin{pmatrix}
1 & 1 \\
0 & 0.9999
\end{pmatrix}.
\end{align*} \end{split}\]
Solving \(L \mathbf{y} = \mathbf{b}\) using the same precision as before
\[\begin{split} \begin{align*}
\begin{pmatrix}
1 & 0 \\
0.0001 & 1
\end{pmatrix}
\begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &=
\begin{pmatrix} 2 \\ 1 \end{pmatrix}
\quad \implies \quad
\begin{aligned}
y_1 &= 2, \\
y_2 &= 1 - 0.0001(2) \approx 1,
\end{aligned}
\end{align*} \end{split}\]
Solving \(U \mathbf{x} = \mathbf{y}\)
\[\begin{split} \begin{align*}
\begin{pmatrix}
1 & 1 \\
0 & 0.9999
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &=
\begin{pmatrix} 2 \\ 0.9998 \end{pmatrix}
\quad \implies \quad
\begin{aligned}
x_2 &= \frac{0.9998}{0.9999} \approx 1, \\
x_1 &\approx 2 - 1 = 1.
\end{aligned}
\end{align*} \end{split}\]
Now we have the solution \(x_1 \approx 1\) and \(x_2 \approx 1\) which is closer to the exact solution of \(x_1 \approx 1.0001\) and \(x_2 \approx 0.9999\). Note that swapping rows of the coefficient matrix does not change the solution to the system as long as the same swap is applied to the constant vector.
In practice, partial pivoting is performed automatically during Gaussian elimination and LU decomposition because it usually improves numerical stability and avoids division by zero.
Definition 6.1 (Partial pivoting)
At the \(k\)-th step of Gaussian elimination, partial pivoting selects as the pivot the entry with the largest absolute values in column \(k\) among rows \(k, k+1, \ldots, n\). If this entry is not already in row \(k\), the corresponding rows are swapped before elimination proceeds.
Formally, for column \(k\), we choose the row \(p\) such that
\[ |a_{pk}| = \max_{i=k, \ldots, n} |a_{ik}|. \]
Rows \(k\) and \(p\) are then swapped. These row swaps can be represented using a permutation matrix \(P\). A permutation matrix is obtained by performing the same row swaps on the identity matrix as are applied to \(A\). The matrix \(PA\) is the matrix obtained after performing all row swaps determined by partial pivoting.
Example 6.3
Apply partial pivoting to the following matrix and determine the permutation matrix \(P\)
\[\begin{split} A =\begin{pmatrix}
0 & 1 & -2\\
1 & 0 & 2\\
3 & -2 & 2
\end{pmatrix}. \end{split}\]
Solution
Applying partial pivoting to \(A\)
\[\begin{split} \begin{align*}
&\begin{pmatrix}
0 & 1 & -2\\
1 & 0 & 2\\
3 & -2 & 2
\end{pmatrix}
\begin{array}{l} R_1 \leftrightarrow R_3 \\ \phantom{x} \\ \phantom{x} \end{array}
\longrightarrow
&\begin{pmatrix}
3 & -2 & 2\\
1 & 0 & 2\\
0 & 1 & -2
\end{pmatrix}
\begin{array}{l} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array}\\ \\
\longrightarrow
&\begin{pmatrix}
3 & -2 & 2\\
0 & 1 & -2\\
1 & 0 & 2
\end{pmatrix}
\begin{array}{l} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array}.
\end{align*} \end{split}\]
Performing the same row operations on the identity matrix
\[\begin{split} \begin{align*}
&\begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{array}{l} R_1 \leftrightarrow R_3 \\ \phantom{x} \\ \phantom{x} \end{array}
\longrightarrow
&\begin{pmatrix}
0 & 0 & 1\\
0 & 1 & 0\\
1 & 0 & 0
\end{pmatrix}
\begin{array}{l} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array}\\ \\
\longrightarrow
&\begin{pmatrix}
0 & 0 & 1\\
1 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}.
\end{align*} \end{split}\]
Therefore, the permutation matrix is
\[\begin{split} P = \begin{pmatrix}
0 & 0 & 1\\
1 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}. \end{split}\]
6.2.1. LUP decomposition
The standard form of an LU decomposition with partial pivoting is
\[ PA = LU, \]
where \(P\) is a permutation matrix. Since partial pivoting consists of row swaps, the same row swaps must be applied to both sides of the system \(A \mathbf{x} = \mathbf{b}.\) Pre-multiplying by \(P\) gives
\[ P A \mathbf{x} = P \mathbf{b}. \]
Using the factorisation \(PA = LU\), we obtain
\[ LU \mathbf{x} = P \mathbf{b}. \]
Introducing an intermediate vector \(\mathbf{y}\) such that
\[ U \mathbf{x} = \mathbf{y}, \]
the system becomes
\[ L \mathbf{y} = P \mathbf{b}. \]
Since \(L\) is lower triangular, \(\mathbf{y}\) is computed using forward substitution. Once \(\mathbf{y}\) has been determined, the solution vector \(\mathbf{x}\) is obtained by solving
\[ U\mathbf{x}=\mathbf{y}, \]
using back substitution.
Algorithm 6.3 (Solving system of linear equations using LU decomposition with partial pivoting)
Inputs: A system of linear equations expressed using the coefficient matrix \(A\) and variable vector \(\mathbf{b}\) such that \(A \mathbf{x} = \mathbf{b}\).
Outputs: The solution vector \(\mathbf{x}\).
Perform partial pivoting on \(A\) and determine the permutation matrix \(P\)
Calculate the LU decomposition of \(PA\)
\(\mathbf{b} \gets P \mathbf{b}\)
\(y_i \gets b_i - \displaystyle\sum_{j=1}^{i - 1} \ell_{ij} y_j, \qquad i = 1, \ldots, n\)
\(x_n \gets \dfrac{y_n}{u_{nn}}\)
\(x_i \gets \dfrac{1}{u_{ii}} \left( y_i - \displaystyle \sum_{j=i+1}^{n} u_{ij} x_j \right), \qquad i = n-1, \ldots, 1\)
Return \(\mathbf{x}\).
Example 6.4
Solve the following system of linear equations using LU decomposition with partial pivoting
\[\begin{split} \begin{align*}
\begin{pmatrix}
0 & 1 & -2\\
1 & 0 & 2\\
3 & -2 & 2
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} =
\begin{pmatrix} 10 \\ -4 \\ -8 \end{pmatrix}.
\end{align*} \end{split}\]
Solution
We have seen from the Example 6.3 that applying partial pivoting to the coefficient matrix results in
\[\begin{split} \begin{align*}
P &= \begin{pmatrix}
0 & 0 & 1\\
1 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}, &
PA &= \begin{pmatrix}
3 & -2 & 2\\
0 & 1 & -2\\
1 & 0 & 2
\end{pmatrix}, &
\end{align*} \end{split}\]
Calculating the LU decomposition of \(PA\)
\[\begin{split} \begin{align*}
j &= 1: & u_{11} &= a_{11} = 3, \\
&& \ell_{21} &= \frac{1}{u_{11}} a_{21} = \frac{1}{3}(0) = 0, \\
&& \ell_{31} &= \frac{1}{u_{11}} a_{31} = \frac{1}{3}(1) = \frac{1}{3}, \\
j &= 2: & u_{12} &= a_{12} = -2, \\
&& u_{22} &= a_{22} - \ell_{21}u_{12} = 1 - 0(-2) = 1, \\
&& \ell_{32} &= \frac{1}{u_{22}}(a_{32} - \ell_{31}u_{12}) = \frac{1}{1}\left(0 - \frac{1}{3}(-2)\right) = \frac{2}{3}, \\
j &= 3: & u_{13} &= a_{13} = 2, \\
&& u_{23} &= a_{23} - \ell_{21}u_{13} = -2 - 0(2) = -2, \\
&& u_{33} &= a_{33} - \ell_{31}u_{13} - \ell_{32}u_{23} = 2 - \frac{1}{3}(2) - \frac{2}{3}(-2) = \frac{8}{3},
\end{align*} \end{split}\]
therefore
\[\begin{split} \begin{align*}
L &= \begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
\frac{1}{3} & \frac{2}{3} & 1
\end{pmatrix}, &
U &= \begin{pmatrix}
3 & -2 & 2\\
0 & 1 & -2\\
0 & 0 & \frac{8}{3}
\end{pmatrix}.
\end{align*} \end{split}\]
Solving \(L \mathbf{y} = P \mathbf{b}\) using forward substitution
\[\begin{split} \begin{align*}
\begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
\frac{1}{3} & \frac{2}{3} & 1
\end{pmatrix}
\begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} =
\begin{pmatrix}
0 & 0 & 1\\
1 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}
\begin{pmatrix} 10 \\ -4 \\ -8 \end{pmatrix} =
\begin{pmatrix} -8 \\ 10 \\ -4 \end{pmatrix},
\end{align*} \end{split}\]
gives
\[\begin{split} \begin{align*}
y_1 &= -8, \\
y_2 &= 10 - 0(8)=10, \\
y_3 &= -4 - \frac{1}{3}(-8) + \frac{2}{3}(10) = -8.
\end{align*} \end{split}\]
Solving \(U \mathbf{x} = \mathbf{y}\) using back substitution
\[\begin{split} \begin{align*}
\begin{pmatrix}
3 & -2 & 2\\
0 & 1 & -2\\
0 & 0 & 8/3
\end{pmatrix}
\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} =
\begin{pmatrix} -8 \\ 10 \\ -8 \end{pmatrix},
\end{align*} \end{split}\]
gives
\[\begin{split} \begin{align*}
x_3 &= \frac{3}{8}(-8) = -3, \\
x_2 &= \frac{1}{1}(10+2(-3)) = 4, \\
x_1 &= \frac{1}{3}(-8+2(4)-2(-3)) = 2.
\end{align*} \end{split}\]
So the solution vector is \(\mathbf{x}=(2,4,-3)\).
