6.6. LU Decomposition

A variant of Gaussian elimination called LU decomposition (for Lower–Upper) factors a matrix into products of two matrices. The first is a lower triangular matrix, and the second is an upper triangular matrix. The LU decomposition algorithm can be implemented in software to achieve the best possible accuracy from an elimination-based algorithm. [1] It is used internally by MATLAB to solve square matrix systems of equations and to calculate determinants. MATLAB’s lu function finds a matrix’s LU decomposition.

Recall that to solve a system of equations, we can use Gaussian elimination with an augmented matrix so that changes made to matrix \(\bf{A}\) are also applied to vector \(\bm{b}\).

\[\begin{array}{lcl} \mathbf{A}\,\bm{x} = \bm{b} & \mapsto & \mathbf{U}\,\bm{x} = \bm{c} \end{array}\]

In the modified equation, \(\bf{U}\) is an upper triangular matrix for which simple back substitution may be used to solve for the unknown vector \(\bm{x}\). The LU decomposition method differs in that it operates only on the \(\bf{A}\) matrix to find a matrix factoring, \(\mathbf{A} = \bf{L\,U}\), where \(\bf{L}\) is lower triangular and \(\bf{U}\) is upper triangular. Thus the substitution phase for solving \(\mathbf{A}\,\bm{x} = \bm{b}\) has two steps—solving \(\mathbf{L}\,\bm{y} = \bm{b}\) and then solving \(\mathbf{U}\,\bm{x} = \bm{y}\).

LU decomposition changes the equations used to solve systems of equations as follows. In addition to \(\bf{L}\) and \(\bf{U}\), we have an elimination matrix, \(\bf{E}\). Its inverse is \(\bf{L}\). We also introduce the \(\bf{P}\) matrix, the permutation matrix that orders the rows to match the row exchanges. [2]

\[\begin{split}\begin{array}{lcl} \mathbf{A} & \mapsto & \mathbf{E\,P\,A} = \mathbf{U} \\ \mathbf{P\,A} = \mathbf{E}^{-1}\mathbf{U} & \mapsto & \mathbf{P\,A} = \mathbf{L\,U} \\ \mathbf{A}\,\bm{x} = \bm{b} & \mapsto & \mathbf{P}^T\mathbf{L\,U}\,\bm{x} = \bm{b} \\ \bm{x} = \mathbf{A}^{-1}\bm{b} & \mapsto & \bm{x} =\mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{P}\,\bm{b} \end{array}\end{split}\]

(6.13)\[\boxed{\bm{x} =\mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{P}\,\bm{b}}\]

Which algorithm is faster?

There are some applications where there is a need to solve several systems of equations with the same \(\mathbf{A}\) matrix, but with different \(\bm{b}\) vectors. It is faster to find the LU factors once and then use forward and backward substitution with each \(\bm{b}\) vector than to do full elimination each time. The combined forward and backward substitution for LU decomposition has a run time of \(\mathcal{O}(n^2)\).

Elimination with an augmented matrix and the LU decomposition algorithm have the same order of magnitude run time performance, which, depending on the exact implementation, is approximately \(\mathcal{O}(n^3)\). But more precisely, the cost of computation for LU decomposition is \((2/3)n^3 + n^2\) floating point operations (flops) [HIGHAM11], which explains why experimentation shows that LU decomposition is faster than the rref function, especially for larger matrices. Here is a quick experiment performed on a personal computer.

>> A = rand(100); b = rand(100,1);

>> tic; [L,U,P] = lu(A);  x = U\(L\(P*b)); toc;
Elapsed time is 0.000605 seconds.

>> tic; x = rref([A b]); toc;
Elapsed time is 0.069932 seconds.

6.6.1. LU Example

Before demonstrating how to use the lu function in MATLAB, we will work on part of a problem and then turn to MATLAB to help with the computation. We will use the same example solved in

Elimination. For purposes of understanding how the \(\bf{U}\) and \(\bf{L}\) matrices are factors of the \(\bf{A}\) matrix, we will follow a traditional elimination strategy to find \(\bf{U}\). We will form \(\bf{L}\) from a record of the row operations needed to convert \(\bf{A}\) into \(\bf{U}\). We will not make row exchanges in this example, but will use row exchanges in

LU Decomposition.

\[\begin{split}\left\{\begin{aligned} 2x - 3y {\quad} {\quad} &= 3 \\ 4x - 5y + z {\:} &= 9 \\ 2x - y - 3z {\,} &= -1 \end{aligned}\right.\end{split}\]

The needed row operations are the same as noted in Elimination. Each row operation is represented by a matrix \(\mathbf{E}_{i,j}\) that can multiply with \(\bf{A}\) to affect the change. These matrices start with an identity matrix and then change one element to represent each operation. We call these matrices Elementary Matrices. Then the product of the elementary matrices is called the elimination matrix.

1. Add \(-1\) of row 1 to row 3, called \(\mathbf{E}_{1}\).

2. Add \(-2\) of row 1 to row 2, called \(\mathbf{E}_{2}\).

3. Add \(-2\) of row 2 to row 3, called \(\mathbf{E}_{3}\).

>> E1
E1 =
    1     0     0
    0     1     0
   -1     0     1
>> E2
E2 =
    1     0     0
   -2     1     0
    0     0     1
>> E3
E3 =
    1     0     0
    0     1     0
    0    -2     1

The combined row operations can be found by multiplying them together. The order of matrices must be such that the first operation applied is next to \(\bf{A}\) in the equation \(\mathbf{U} = \bf{E\,A}\).

>> E = E3*E2*E1
E =
    1     0     0
   -2     1     0
   -1    -2     1
>> U = E*A
U =
    2    -3     0
    0     1     1
    0     0    -5

We can also use an augmented matrix to find E and U. Start with \(\left[\,\mathbf{A}\,|\, \mathbf{I}\, \right]\) and perform row operations until the matrix in the position of \(\bf{A}\) is upper triangular. The resulting augmented matrix is then \(\left[\,\mathbf{U}\,|\, \mathbf{E}\, \right]\) [KUTZ13].

\[\begin{split}\left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 4 & -5 & 1 & 0 & 1 & 0 \\ 2 & -1 & -3 & 0 & 0 & 1 \end{array}\right] % \end{equation*}% % \begin{equation*} \; \mapsto \; \left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & -2 & 1 & 0 \\ 0 & 0 & -5 & -1 & -2 & 1 \end{array}\right]\end{split}\]

The \(\bf{L}\) matrix is the inverse of \(\bf{E}\). Recall that to take the inverse of a product of matrices, we reverse the order of the matrices—\((\mathbf{A\,B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}\).

Note that each elementary matrix only differs from the identity matrix at one element. The inverse of each elementary matrix is found by changing the sign of the additional nonzero value in the identity matrix.

>> inv(E1)
ans =
    1     0     0
    0     1     0
    1     0     1

>> inv(E2)
ans =
    1     0     0
    2     1     0
    0     0     1

>> inv(E3)
ans =
    1     0     0
    0     1     0
    0     2     1
>> L = inv(E1)*inv(E2)*inv(E3)
L =
    1     0     0
    2     1     0
    1     2     1

>> inv(E)
ans =
    1.0000         0         0
    2.0000    1.0000         0
    1.0000    2.0000    1.0000

Also note that the \(\bf{L}\) matrix does the reverse of \(\bf{E}\), so we could find \(\bf{L}\) directly from the list of row operations.

Now, with \(\bf{L}\) and \(\bf{U}\) found, we can find the solution to \(\mathbf{L\,U}\,\bm{x} = \bm{b}\). Since these are lower and upper triangular matrices, the left-divide operator will quickly solve them. If several solutions are needed, the one-time calculation of \(\bf{L}\) and \(\bf{U}\) will substantially expedite the calculations.

>> b
b =
    3
    9
   -1
>> x = U\(L\b)
x =
    3.0000
    1.0000
    2.0000

Note

The substitution phase has two steps: solving \(\mathbf{L}\,\bm{y} = \bm{b}\) and \(\mathbf{U}\,\bm{x} = \bm{y}\). Either another vector, \(\bm{y}\), is needed, or parentheses may be used to dictate the ordering of the operations.

6.6.2. Example with Row Exchanges

This example shows the LU decomposition with partial pivoting (row exchanges). We will show the row exchanges with a permutation matrix. A permutation matrix has the rows of an identity matrix, but the order of the rows shows what row exchanges were performed.

\[\begin{split}\mathbf{A} = \begin{bmatrix} 0 & 12 & -3 \\ 8 & -4 & -6 \\ -4 & -2 & 12 \end{bmatrix}\end{split}\]

Looking at the first column, we will use the largest value, 8, as our first pivot, so we exchange rows one and two and define a permutation matrix. The \(\mathbf{U}\) matrix will reflect the row exchanges. The \(\mathbf{L}\) matrix begins as an identity matrix, and the ones on the diagonal will remain. We will document our elimination row operations directly in the \(\mathbf{L}\) matrix as the negative values of what we would add to the elementary matrices.

\[\begin{split}\mathbf{P} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

\[\begin{split}\begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ -4 & -2 & 12 \end{bmatrix} \end{array}\end{split}\]

Add \(1/2\) of row 1 to row 3.

\[\begin{split}\begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1/2 & 0 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ 0 & -4 & 9 \end{bmatrix} \end{array}\end{split}\]

We will keep the pivot of \(12\) in row 2, so no change to the permutation matrix is needed. Add \(1/3\) of row 2 to row 3.

\[\begin{split}\begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1/2 & -1/3 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ 0 & 0 & 8 \end{bmatrix} \end{array}\end{split}\]

6.6.3. MATLAB Examples of LU

First, we ask MATLAB to compute \(\bf{L}\) and \(\bf{U}\). We can see that \(\bf{L}\) is not lower triangular as we think it should be. That is because we did not accept a permutation matrix in the output.

>> [L, U] = lu(A)
L =
    0.5000   -0.3333    1.0000
    1.0000         0         0
    0.5000    1.0000         0
U =
    4.0000   -5.0000    1.0000
         0    1.5000   -3.5000
         0         0   -1.6667

Now we will accept a permutation matrix as an output matrix. The relationship of the matrices is \(\mathbf{P\,A} = \mathbf{L\,U}\). Note that the permutation matrix merely swaps the rows of whatever vector or matrix it is multiplied by. Like the identity matrix, permutation matrices are orthogonal, thus their inverse is just their transpose, so \(\mathbf{A} = \mathbf{P}^{T}\bf{L\,U}\). Now, \(\bf{L}\) is lower triangular as expected. So the solution comes from equation (6.13), or in MATLAB, x = U\(L\(P*b)).

% Returning to
% our previous example ...
>> [L, U, P] = lu(A)
L =
    1.0000         0         0
    0.5000    1.0000         0
    0.5000   -0.3333    1.0000
U =
    4.0000   -5.0000    1.0000
         0    1.5000   -3.5000
         0         0   -1.6667
P =
    0     1     0
    0     0     1
    1     0     0
>> x = U\(L\(P*b))
x =
    3.0000
    1.0000
    2.0000

Another option is to have MATLAB return the row exchange information as an ordering of the rows rather than as a permutation matrix.

>> [L, U, P] = ...
    lu(A, 'vector');
>> P
P =
    2     3     1
>> x = U\(L\(b(P)))
x =
    3.0000
    1.0000
    2.0000

6.6.4. Turing’s LU Elimination Algorithm

There are several implementations of LU decomposition. In this section, we present a simple, yet efficient, and numerically stable LU algorithm that was developed and reported in a 1948 paper by British mathematician Alan Turing, who is considered to be a founding father of computer science [TURING48], [DOPICO13]. Modern numerical linear algebra literature calls Turing’s algorithm the “kij” version of the outer product LU algorithm. Golub and Van Loan’s text on Matrix Computations gives a proof showing that potential round-off errors are minimized by this algorithm [GOLUB13].

The simple function turingLU implements Turing’s kij algorithm with a few added lines of code for partial pivoting.

function [L, U, P] = turingLU(A)
% [L,U,P] = TURINGLU(A)
% A simplified LU factorization of a matrix using Alan Turing's
% kij variant of the Gaussian elimination algorithm.
%
% Input:  A - square matrix
% Output: L - lower triangular matrix
%         U - upper triangular matrix
%         P - permutation matrix
% P*A = L*U, A = P'*L*U
%
% This code is for educational purposes only.
% For any other purpose, use MATLAB's lu function.

[m, n] = size(A);
if m ~= n
    error('Matrix is not square')
end
P = eye(n);
% The kij nested elimination loop
for k = 1:(n - 1)
    [A(k:n,:), idx] = sortrows(A(k:n,:), k, 'descend', ...
        'ComparisonMethod', 'abs');
    I = P(k:n,:);
    P(k:n,:) = I(idx,:); % Permutation matrix
    for i = k + 1:n
        A(i, k) = A(i, k)/A(k, k);
        for j = (k + 1):n
            A(i, j) = A(i, j) - A(i, k) * A(k, j);
        end
    end
end
% L is now in the lower triangular of A.
% U is now in the upper triangular of A.
L = tril(A, -1) + eye(n); % extract lower
U = triu(A);              % and upper triangular

The \(\bf{L}\) and \(\bf{U}\) matrices are taken directly from the modified \(\bf{A}\) matrix after the algorithm is completed. After the elimination steps are finished, the \(\bf{L}\) matrix is in the lower triangular of \(\bf{A}\) and the \(\bf{U}\) matrix is in the upper triangular of \(\bf{A}\).

From a first reading of the code, it is a little difficult to see how it achieves the same result as found from normal Gaussian elimination. With elimination, row operations zero out the lower triangular part of the matrix, and then \(\bf{L}\) is found from the action of the row operations. Whereas, in Turing’s algorithm, the row operations directly yield both \(\bf{L}\) and \(\bf{U}\).

The variable k references the pivot variable. Row exchanges are performed each time the pivot is advanced. There are only two statements that alter the values of the matrix. First, the elements in the column below the pivot are divided by the pivot, which equates to recording the row operations in the \(\bf{L}\) matrix. Setting those values to zero is unnecessary since they will become part of the \(\bf{L}\) matrix. Secondly, row operations are performed between the pivot row and the rows below the pivot row, preparing for future steps that set the \(\bf{L}\) matrix values at positions left of the diagonal. It also sets the needed \(\bf{U}\) matrix values for elements on or to the right of the diagonal.

6.6.5. Determinant Shortcut

Another application of LU decomposition is a faster algorithm for computing the determinant of a matrix. The product of the determinants of a matrix factorization is equal to the determinant of the original matrix. Thus,

\[|\mathbf{A}| = |\mathbf{P}^T|\,|\mathbf{L}|\,|\mathbf{U}|.\]

Because of the zeros in all three of these matrices, their determinants are quick to compute. It is not necessary to take the transpose of \(\bf{P}\) before finding its determinant because the determinant of a square matrix is the same as the determinant of its transpose, det(P) = det(P’). The determinant of \(\bf{P}\) is either 1 or -1. Since each row has only a single one, the determinant is quick to calculate. The determinant of \(\bf{L}\) is always 1, which is the product of the ones along the diagonal. The determinant of \(\bf{U}\) is the product of the values on its diagonal.

Refer to the normal method for calculating the determinant in Determinant and take a moment to see why the determinant of an upper or lower triangular matrix is the product of the diagonal values.

The MATLAB det function computes the determinant of a matrix as follows.

>> [L, U, P] = lu(A)
>> determinant = det(P)*prod(diag(U))

Footnotes

[1]

Matrices in a special form, such as diagonal, triangular, or positive definite, have more efficient solvers.

[2]

Equation (6.13) is expressed with a pure mathematics notation using the inverse of matrices. The expression that is preferred with MATLAB is an expression that directly uses back substitution on triangular matrices (x = U\(L\(P*b))). Note that the parentheses are important to ensure that substitution is applied in the correct order.