Compute the LU factorization of a matrix and examine the resulting factors. LU factorization is a way of decomposing a matrix $\mathit{A}$ into an upper triangular matrix $\mathit{U}$, a lower triangular matrix $\mathit{L}$, and a permutation matrix $\mathit{P}$ such that $\mathrm{PA}=\mathrm{LU}$. These matrices describe the steps needed to perform Gaussian elimination on the matrix until it is in reduced row echelon form. The $\mathit{L}$ matrix contains all of the multipliers, and the permutation matrix $\mathit{P}$ accounts for row interchanges.

Create a 3-by-3 matrix and calculate the LU factors.

A = [10 -7 0
     -3  2 6
      5 -1 5];

[L,U] = lu(A)

L = 3×3

    1.0000         0         0
   -0.3000   -0.0400    1.0000
    0.5000    1.0000         0

U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000

Multiply the factors to recreate A. With the two-input syntax, lu incorporates the permutation matrix P directly into the L factor, such that the L being returned is really P'*L and thus A = L*U.

L*U

ans = 3×3

   10.0000   -7.0000         0
   -3.0000    2.0000    6.0000
    5.0000   -1.0000    5.0000

You can specify three outputs to separate the permutation matrix from the multipliers in L.

[L,U,P] = lu(A)

L = 3×3

    1.0000         0         0
    0.5000    1.0000         0
   -0.3000   -0.0400    1.0000

U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000

P = 3×3

     1     0     0
     0     0     1
     0     1     0

P'*L*U

ans = 3×3

   10.0000   -7.0000         0
   -3.0000    2.0000    6.0000
    5.0000   -1.0000    5.0000

Solve a linear system by performing an LU factorization and using the factors to simplify the problem. Compare the results with other approaches using the backslash operator and decomposition object.

Create a 5-by-5 magic square matrix and solve the linear system $\mathrm{Ax}=\mathit{b}$ with all of the elements of b equal to 65, the magic sum. Since 65 is the magic sum for this matrix (all of the rows and columns add to 65), the expected solution for x is a vector of 1s.

A = magic(5);
b = 65*ones(5,1);
x = A\b

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

For generic square matrices, the backslash operator computes the solution of the linear system using LU decomposition. LU decomposition expresses A as the product of triangular matrices, and linear systems involving triangular matrices are easily solved using substitution formulas.

To recreate the answer computed by backslash, compute the LU decomposition of A. Then, use the factors to solve two triangular linear systems:

y = L\(P*b);
x = U\y;

This approach of precomputing the matrix factors prior to solving the linear system can improve performance when many linear systems will be solved, since the factorization occurs only once and does not need to be repeated.

[L,U,P] = lu(A)

L = 5×5

    1.0000         0         0         0         0
    0.7391    1.0000         0         0         0
    0.4783    0.7687    1.0000         0         0
    0.1739    0.2527    0.5164    1.0000         0
    0.4348    0.4839    0.7231    0.9231    1.0000

U = 5×5

   23.0000    5.0000    7.0000   14.0000   16.0000
         0   20.3043   -4.1739   -2.3478    3.1739
         0         0   24.8608   -2.8908   -1.0921
         0         0         0   19.6512   18.9793
         0         0         0         0  -22.2222

P = 5×5

     0     1     0     0     0
     1     0     0     0     0
     0     0     0     0     1
     0     0     1     0     0
     0     0     0     1     0

y = L\(P*b);
x = U\y

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

The decomposition object also is useful to solve linear systems using specialized factorizations, since you get many of the performance benefits of precomputing the matrix factors but you do not need to know how to use the factors. Use the decomposition object with the 'lu' type to recreate the same results.

dA = decomposition(A,'lu');
x = dA\b

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

Compute the LU factorization of a sparse matrix and verify the identity L*U = P*S*Q.

Create a 60-by-60 sparse adjacency matrix of the connectivity graph of the Buckminster-Fuller geodesic dome.

S = bucky;

Compute the LU factorization of S using the sparse matrix syntax with four outputs to return the row and column permutation matrices.

[L,U,P,Q] = lu(S);

Permute the rows and columns of S with P*S*Q and compare the result with multiplying the triangular factors L*U. The 1-norm of their difference is within roundoff error, indicating that L*U = P*S*Q.

e = P*S*Q - L*U;
norm(e,1)

ans = 2.4425e-15

Compute the LU factorization of a matrix. Save memory by returning the row permutations as a vector instead of a matrix.

Create a 1000-by-1000 random matrix.

A = rand(1000);

Compute the LU factorization with the permutation information stored as a matrix P. Compare the result with the permutation information stored as a vector p. The larger the matrix, the more memory efficient it is to use a permutation vector.

[L1,U1,P] = lu(A);
[L2,U2,p] = lu(A,'vector');
whos P p

  Name         Size                Bytes  Class     Attributes

  P         1000x1000            8000000  double              
  p            1x1000               8000  double              

Using a permutation vector also saves on execution time in subsequent operations. For instance, you can use the previous LU factorizations to solve a linear system $\mathrm{Ax}=\mathit{b}$. Although the solutions obtained from the permutation vector and permutation matrix are equivalent (up to roundoff), the solution using the permutation vector typically requires a little less time.

Compare the results of computing the LU factorization of a sparse matrix with and without column permutations.

Load the west0479 matrix, which is a real-valued 479-by-479 sparse matrix.

load west0479
A = west0479;

Calculate the LU factorization of A by calling lu with three outputs. Generate spy plots of the L and U factors.

[L,U,P] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')

Now, calculate the LU factorization of A using lu with four outputs, which permutes the columns of A to reduce the number of nonzeros in the factors. The resulting factors are much sparser than if column permutations are not used.

[L,U,P,Q] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')