6.3. LU decomposition with partial pivoting#

In the previous section we saw that the elements of \(L\) in the LU decomposition of \(A\) is

\[ \begin{align*} \ell_{ij} &= \dfrac{1}{u_{jj}} \left(a_{ij} - \displaystyle \sum_{k=1}^{j-1} \ell_{ik}u_{kj}\right), & i &= j+1, \ldots, n. \end{align*} \]

A problem can occur if the value of \(u_{jj}\) in the expression is zero or some small number it will mean that it is either undefined or prone to computational rounding errors due to the value of \(\ell_{ij}\) being very large. For example, consider the solution of the following system using LU decomposition

\[\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}\]

Calculating the LU decomposition of the coefficient matrix we have

\[\begin{split} \begin{align*} L &= \begin{pmatrix} 1 & 0 \\ 10000 & 1 \end{pmatrix}, & U &= \begin{pmatrix} 1 & 1 \\ 0 & -9999 \end{pmatrix}, \end{align*} \end{split}\]

and solving using forward and back substitution and working to four decimal place accuracy

\[\begin{split} \begin{align*} L \vec{y} &= \vec{b} \\ \begin{pmatrix} 1 & 0 \\ 10000 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &= \begin{pmatrix} 1 \\ 2 \end{pmatrix}\\ \therefore y_1 &= 1, \\ y_2 &= 2 - 10000(1) = -9998 , \\ \\ U \vec{x} &= \vec{y} \\ \begin{pmatrix} 1 & 1 \\ 0 & -9999 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 1 \\ -9998 \end{pmatrix} \\ \therefore x_2 &= \frac{9998}{9999} = 0.9998\dot{9} \approx 1,\\ x_1 &= 1 - 1 = 0. \end{align*} \end{split}\]

So we have a solution \(x_1 = 0\) and \(x_2 = 1\), checking this is correct by substituting into the original system gives

\[\begin{split} \begin{align*} \vec{b} &= \begin{pmatrix} 0.0001 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \end{align*} \end{split}\]

which is not equal to the constant vector. This problem can be overcome by using partial pivoting where rows of the coefficient matrix are permuted (swapped) to ensure that the value on the main diagonal of \(A\) is larger than the other values in the column beneath it. For example, if we swap the rows of the system from before we have

\[\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}\]

which has the LU decomposition

\[\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 using forward and back substitution with the same decimal place accuracy used before

\[\begin{split} \begin{align*} L \vec{y} &= \vec{b} \\ \begin{pmatrix} 1 & 0 \\ 0.0001 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &= \begin{pmatrix} 2 \\ 1 \end{pmatrix}\\ \therefore y_1 &= 2, \\ y_2 &= 1 - 0.0001(2) = 0.9998 , \\ \\ U \vec{x} &= \vec{y} \\ \begin{pmatrix} 1 & 1 \\ 0 & 0.9999 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 2 \\ 0.9998 \end{pmatrix} \\ \therefore x_2 &= \frac{0.9998}{0.9999} = 0.9998\dot{9} \approx 1, \\ x_1 &= 2 - 1 = 1. \end{align*} \end{split}\]

Now we have the solution \(x_1 = 1\) and \(x_2 = 1\) which substituted into the original system gives

\[\begin{split} \begin{align*} \vec{b} &= \begin{pmatrix} 0.0001 & 1 \\ 1 & 1 \\ \end{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1.0001 \\ 2 \end{pmatrix} \approx \begin{pmatrix} 1 \\ 2 \end{pmatrix}, \end{align*} \end{split}\]

which is correct. 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. The permutations applied to the coefficient matrix are recorded in a matrix \(P\) which is determined by applying the same permutations to the identity matrix.

Theorem 6.2 (Permutation matrix)

A sequence of row swaps applied to a matrix \(A\) can be recorded in a permutation matrix \(P\) by performing the same sequence to the identity matrix. Furthermore, the result of applying the row swaps to \(A\) can be determined by \(PA\).

Algorithm 6.3 (Partial pivoting)

Inputs: An \(m \times n\) matrix \(A\)

Outputs: The permutation matrix \(P\)

  • \(P \gets I\)

  • Start with the pivot element as the element in the first row and column of \(A\).

  • For each column of \(A\)

    • Perform a row swap on \(A\) with the row beneath the pivot row that has the largest absolute value element in the pivot column which is greater than the absolute value fo the pivot element.

    • Perform the same row swap on \(P\).

    • If the pivot element is zero move to the next column.

  • Return \(P\)

Example 6.3

Apply partial pivoting 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 (click to show)

Using elementary row operations

\[\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} \phantom{x} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array} \\ \\ \longrightarrow &\begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}. \end{align*} \end{split}\]

Apply the same row operations to 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} \phantom{x} \\ 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}\]

So \(P=\begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix}\). Note that pre-multiplying \(A\) by \(P\) gives the matrix \(A\) after partial pivoting has been applied

\[\begin{split} \begin{align*} \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} = \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}. \end{align*} \end{split}\]

6.3.1. Code#

The code below defines a function called partial_pivot() which performs partial pivoting on the matrix A and returns the pivoted matrix A and the permutation matrix P.

def partial_pivot(A):
    n = A.shape[0]
    P = np.eye(n)
    for j in range(n):
        maxpivot, maxpivotrow = abs(A[j,j]), j
        for i in range(j + 1, n):
            if abs(A[i,j]) > maxpivot:
                maxpivot, maxpivotrow = abs(A[i,j]), i

        A[[j,maxpivotrow]] = A[[maxpivotrow,j]] 
        P[[j,maxpivotrow]] = P[[maxpivotrow,j]]

    return P

function P = partial_pivot(A)

n = size(A, 1);
P = eye(n);
for j = 1 : n
    maxpivot = abs(A(j,j));
    maxpivotrrow = j;
    for i = j + 1 : n
        if abs(A(i,j)) > maxpivot
            maxpivot = abs(A(i,j));
            maxpivotrow = i;
        end
    end
    A([j,maxpivotrow],:) = A([maxpivotrow,j],:);
    P([j,maxpivotrow],:) = P([maxpivotrow,j],:);
end

end

6.3.2. LUP decomposition#

To implement LU decomposition with partial pivoting (LUP decomposition) we apply partial pivoting to the coefficient matrix of a system to determine a permutation matrix \(P\) before calculating the LU decomposition of \(PA\), i.e.,

\[ LU = PA. \]

Since we have applied row swaps to the coefficient matrix we must also apply the same row swaps to the constant vector. So to solve the system \(A \vec{x} = \vec{b}\) using LUP decomposition we have

\[\begin{split} \begin{align*} PA \vec{x} &= P \vec{b} \\ \therefore LU \vec{x} &= P \vec{b}, \end{align*} \end{split}\]

and we solve \(L \vec{y} = P \vec{b}\) for \(\vec{y}\) using forward substitution and \(U \vec{x} = \vec{y}\) for \(\vec{x}\) using back substitution.

Algorithm 6.4 (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 \(\vec{b}\) such that \(A \vec{x} = \vec{b}\).

Inputs: The variable vector \(\vec{x}\).

  • Apply partial pivoting to calculated the permutation matrix \(P\).

  • Calculate the LU decomposition of \(PA\) to determine \(L\) and \(U\) such that \(PA = LU\).

  • \(\vec{b} \gets P \vec{b}\)

  • For \(i = 1, \ldots, n\) do

    • \(y_i \gets b_i - \displaystyle\sum_{j=1}^i \ell_{ij} y_j\)

  • \(x_n \gets \dfrac{y_i}{u_{n,n}}\)

  • For \(i = n - 1, \ldots, 1\) do

    • \(x_i \gets \dfrac{1}{u_{ii}} \left( y_i - \displaystyle \sum_{j=i+1}^{n} u_{ij} x_j \right)\)

  • Return \(\vec{x} = (x_1, x_2, \ldots, x_n)\).

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 (click to show)

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}, \\ \therefore PA &= \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} = \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}, & \end{align*} \end{split}\]

Calculating the LU decomposition of \(PA\) using equation (6.1)

\[\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\\ 1/3 & 2/3 & 1 \end{pmatrix}, & U &= \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 0 & 0 & 8/3 \end{pmatrix}. \end{align*} \end{split}\]

Solving \(L \vec{y} = P \vec{b}\) using forward substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 1/3 & 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 \vec{x} = \vec{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 is \(\vec{x}=(2,4,-3)\).

6.3.3. Python code#

The code below calculates the solution to the system of linear equations from Example 6.4 using LUP decomposition.

# Define linear system
A = np.array([[0, 1, -2],
              [1, 0, 2],
              [3, -2, 2]])
b = np.array([10, -4, -8])

# Calculate LUP decomposition
P = partial_pivot(A)
L, U = lu(A)
print(f"L = \n{L},\n\nU = \n{U},\n\nP = \n{P}")

# Solve linear system
y = forward_substitution(L, np.dot(P, b))
x = back_substitution(U, y)
print(f"\nx = {x}")

The NumPy command np.dot(A, B) calculates the matrix multiplication \(AB\).

% Define linear system
A = [0, 1, -2 ;
    1, 0, 2 ;
    3, -2, 2];
b = [10 ; -4 ; -8];

% Calculate LUP decomposition
P = partial_pivot(A)
[L, U] = lu(P * A);

% Solve linear system
y = forward_substitution(L, P * b);
x = back_substitution(U, y)