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 an upper triangular matrix. The LU decomposition algorithm can be implemented in software to achieve the best possible accuracy from a general elimination based algorithm [1]. It is used internally by MATLAB for computing inverses, solving square matrix systems of equations, and calculating 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, which is the permutation matrix that orders the rows to match the row exchanges. [2] Note that the parentheses are important to ensure that substitution is applied in the correct order.

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

(6.16)¶ $\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 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 first work part of a problem manually and then turn to MATLAB to help with the computation. We will use the same example that was 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}$ and will form $\bf{L}$ from a record of the row operations needed to convert $\bf{A}$ into $\bf{U}$ . We will focus on how the $\bf{U}$ and $\bf{L}$ matrices are found. We will not make row exchanges in this example, but will use row exchanges in

LU Decomposition.

$\spalignsys{2x - 3y \+ \. = 3; 4x - 5y + z = 9; 2x - y - 3z = -1}$

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.

Add $-1$ of row 1 to row 3, called $\mathbf{E}_{1}$ .
Add $-2$ of row 1 to row 2, called $\mathbf{E}_{2}$ .
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 these 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

Another way to find E and U is with elimination of an augmented matrix. 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].

$\spalignaugmathalf{2 -3 0 1 0 0; 4 -5 1 0 1 0; 2 -1 -3 0 0 1}$

$\spalignaugmathalf{2 -3 0 1 0 0; 0 1 1 -2 1 0; 0 0 -5 -1 -2 1}$

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 just changing the sign of the additional nonzero value to 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 $\bf{L}$ could be found 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 calculations 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}$ . When presented with upper or lower triangular matrices, the left-divide operator will skip directly to substitution. So either another vector, $\bm{y}$ is needed, or parenthesis may be used to dictate the ordering of the operations.

6.6.2. Example with Row Exchanges¶

We find the LU decomposition with partial pivoting (row exchanges) in this example. 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.

$\mathbf{A} = \mat{ 0 12 -3; 8 -4 -6; -4 -2 12}$

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.

$\mathbf{P} = \mat{0 1 0; 1 0 0; 0 0 1}$

$\begin{array}{ll} \mathbf{L} = \mat{1 0 0; 0 1 0; 0 0 1} &\mathbf{U} = \mat{8 -4 -6; 0 12 -3; -4 -2 12} \end{array}$

Add $1/2$ of row 1 to row 3.

$\begin{array}{ll} \mathbf{L} = \mat{1 0 0; 0 1 0; -1/2 0 1} &\mathbf{U} = \mat{8 -4 -6; 0 12 -3; 0 -4 9} \end{array}$

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{array}{ll} \mathbf{L} = \mat{1 0 0; 0 1 0; -1/2 -1/3 1} &\mathbf{U} = \mat{8 -4 -6; 0 12 -3; 0 0 8} \end{array}$

6.6.3. MATLAB Examples of LU¶

First, we ask MATLAB to just 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, but $\bf{L}$ and $\bf{U}$ are different than they were with the manual solution because of the permutation matrix. So the solution comes from equation (6.16), or in MATLAB, x = U\(L\(P*b)). Note here that the parenthesis are important. The left-divide operators will not use the correct matrices without the parenthesis.

% 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¶

The MATLAB documentation for lu states that “The LU factorization is computed using a variant of Gaussian elimination”. 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 call 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].

A 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 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}$ .

Turing’s algorithm has good performance and produces accurate results. Unfortunately, the algorithm is specific to LU decomposition and not usable with traditional elimination.

It is a little difficult from a first reading of the code 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}$ an $\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. This equates to recording the row operations in the $\bf{L}$ matrix. There is no need to set those values to zero 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. This operation is preparation 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.

Let’s work through a step-by-step elimination example with a $2 {\times} 2$ matrix to see how the algorithm produces the same result as achieved by Gaussian elimination. We will use letters for the matrix values so that we don’t let fraction reductions cause us to lose track of the numeric operations.

$\mathbf{A} = \mat{a, b; c d}$

The elimination step is to add $\frac{-c}{a}$ of row 1 to row 2. In the $\bf{L}$ matrix, $\frac{c}{a}$ should be put in the (2, 1) position. The $\bf{U}$ matrix now reflects what we would see in the Gaussian elimination matrix. We show what the matrix looks like below in the LU algorithm and the resulting $\bf{L}$ and $\bf{U}$ matrices.

$\mathbf{A}_{\text{final}} = \mat{a, b; \frac{c}{a}, {d - \frac{c\,d}{a}}} \quad \mathbf{L} = \mat{1, 0; \frac{c}{a}, 1} \quad \mathbf{U} = \mat{a, b; 0, {d - \frac{c\,d}{a}}}$

Experimentation with this function found that it reliably returns the same results as MATLAB’s lu function. Please experiment with it yourself. Set breakpoints for the debugging tool or add printing statements to learn how the algorithm works. But for applications other than education, use the lu function that comes with MATLAB, which runs faster, and has more flexible usage options. Note that this function could be implemented with vectorized array indices, but nested for loops are used for clarity of reading the code.

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 triangular, positive definite, Hermitian, or upper Hessenberg have more efficient solvers that are briefly mentioned in Left-Divide Operator. The LU decomposition is a general algorithm because it is applicable to all non-singular, square matrix systems.

[2]	Equation (6.16) 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))`).

6.6. LU Decomposition¶

6.6.1. LU Example¶

6.6.2. Example with Row Exchanges¶

6.6.3. MATLAB Examples of LU¶

6.6.4. Turing’s LU Elimination Algorithm¶

6.6.5. Determinant Shortcut¶

Table of Contents

Previous topic

Next topic

This Page