lu

LU matrix factorization

Syntax

[L,U] = lu(A)

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

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

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

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

[___] = lu(S,thresh)

[___] = lu(___,outputForm)

Description

[L,U] = lu(A) factorizes the full or sparse matrix A into an upper triangular matrix U and a permuted lower triangular matrix L such that A = L*U.

[L,U,P] = lu(A) also returns a permutation matrix P such that A = P'*L*U. With this syntax, L is unit lower triangular and U is upper triangular.

[L,U,P] = lu(A,outputForm) returns P in the form specified by outputForm. Specify outputForm as 'vector' to return P as a permutation vector such that A(P,:) = L*U.

[L,U,P,Q] = lu(S) factorizes sparse matrix S into a unit lower triangular matrix L, an upper triangular matrix U, a row permutation matrix P, and a column permutation matrix Q, such that P*S*Q = L*U.

[L,U,P,Q,D] = lu(S) also returns a diagonal scaling matrix D such that P*(D\S)*Q = L*U. Typically, the row-scaling leads to a sparser and more stable factorization.

[___] = lu(S,thresh) specifies thresholds for the pivoting strategy employed by lu using any of the previous output argument combinations. Depending on the number of output arguments specified, the default value and requirements for the thresh input are different. See the thresh argument description for details.

[___] = lu(___,outputForm) returns P and Q in the form specified by outputForm. Specify outputForm as 'vector' to return P and Q as permutation vectors. You can use any of the input argument combinations in previous syntaxes.

Examples

LU Factorization of Matrix

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    -7     0
    -3     2     6
     5    -1     5

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    -7     0
    -3     2     6
     5    -1     5

Solve Linear System with LU Factorization

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

LU Factorization of Sparse Matrix

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 round-off error, indicating that L*U = P*S*Q.

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

ans = 
2.2204e-16

Save Memory with Permutation Vectors

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 round-off), the solution using the permutation vector typically requires a little less time.

Reduce Fill-in of Sparse Matrix Factorization

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')

Input Arguments

A — Input matrix

matrix

Input matrix. A can be full or sparse as well as square or rectangular in size.

Data Types: single | double
Complex Number Support: Yes

S — Sparse input matrix

sparse matrix

Sparse input matrix. S can be square or rectangular in size.

Data Types: double
Complex Number Support: Yes

thresh — Pivoting thresholds for sparse matrices

scalar | two-element vector

Pivoting thresholds for sparse matrices, specified as a scalar or two-element vector. Valid values are in the interval [0 1]. The way you specify thresh depends on how many outputs are specified in the call to lu:

• For three outputs or less, thresh must be a scalar, and the default value is 1.0.

• For four outputs or more, thresh can be a scalar or a two element vector. The default value is [0.1 0.001]. If you specify thresh as a scalar, then that only replaces the first value in the vector.

At a high level, this input enables you to make trade-offs between accuracy and total execution time. Smaller values of thresh tend to lead to sparser LU factors, but the solution can become inaccurate. Larger values can lead to a more accurate solution (but not always), and usually an increase in the total work and memory usage.

lu selects a pivoting strategy based first on the number of output arguments and second on the properties of the matrix being factorized. In all cases, setting the threshold value(s) to 1.0 results in partial pivoting, while setting them to 0 causes the pivots to be chosen only based on the sparsity of the resulting matrix. All values of L have an absolute value of 1/min(thresh) or less.

• Three or fewer output arguments — The algorithm selects the diagonal pivot if it satisfies the equation

  A(j,j) >= thresh * max(abs(A(j:m,j)))

  Otherwise, it selects the row that contains the element of largest absolute value.

• Symmetric Pivoting Strategy — If S is a square sparse matrix with a mostly symmetric structure and mostly nonzero diagonal, then lu uses a symmetric pivoting strategy. For this strategy, the algorithm selects the diagonal pivot j if it satisfies the inequality:

  A(i,j) >= thresh(2) * max(abs(A(j:m,j)))

  If the diagonal entry fails this test, then lu selects the sparsest row i satisfying the inequality:

  A(i,j) >= thresh(1) * max(abs(A(j:m,j)))

• Nonsymmetric Pivoting Strategy — If S does not satisfy the requirements for the symmetric pivoting strategy, then lu uses a nonsymmetric strategy. In this case, lu selects the sparsest row i satisfying the inequality:

  A(i,j) >= thresh(1) * max(abs(A(j:m,j)))

  A value of 1.0 for thresh(1) results in conventional partial pivoting. Entries in L have an absolute value of 1/thresh(1) or less. The second element of the thresh input vector is not used with the nonsymmetric strategy.

Note

In some rare cases, an incorrect factorization results in P*S*Q ≠ L*U. If this occurs, increase thresh to a maximum of 1.0 (regular partial pivoting), and try again.

outputForm — Shape of permutation outputs

'matrix' (default) | 'vector'

Shape of permutation outputs, specified as 'matrix' or 'vector'. This flag controls whether lu returns the row permutations P and column permutations Q as permutation matrices or permutation vectors.

As matrices, the outputs P and Q satisfy these identities:

• Three outputs — P satisfies P*A = L*U.

• Four outputs — P and Q satisfy P*S*Q = L*U.

• Five outputs — P, Q, and D satisfy P*(D\S)*Q = L*U.

As vectors, the outputs P and Q satisfy these identities:

• Three outputs — P satisfies A(P,:) = L*U

• Four outputs — P and Q satisfy S(P,Q) = L*U

• Five outputs — P, Q, and D satisfy D(:,P)\S(:,Q) = L*U.

Example: [L,U,P] = lu(A,'vector')

Output Arguments

L — Lower triangular factor

matrix

Lower triangular factor, returned as a matrix. The form of L depends on whether the row permutations P are returned in a separate output:

• If the third output P is specified, then L is returned as a unit lower triangular matrix (that is, a lower triangular matrix with 1s on the main diagonal).

• If the third output P is not specified, then L is returned as a row-permutation of a unit lower triangular matrix. Specifically, it is the product P'*L of the outputs P and L returned in the three output case.

U — Upper triangular factor

matrix

Upper triangular factor, returned as an upper triangular matrix.

P — Row permutation

vector | matrix

Row permutation, returned as a permutation matrix or, if the 'vector' option is specified, as a permutation vector. Use this output to improve the numerical stability of the calculation.

See outputForm for a description of the identities that this output satisfies.

Q — Column permutation

vector | matrix

Column permutation, returned as a permutation matrix or, if the 'vector' option is specified, as a permutation vector. Use this output to reduce the fill-in (number of nonzeros) in the factors of a sparse matrix.

See outputForm for a description of the identities that this output satisfies.

D — Row scaling

diagonal matrix

Row scaling, returned as a diagonal matrix. D is used to scale the values in S such that P*(D\S)*Q = L*U. Typically, but not always, the row scaling leads to a sparser and more stable factorization.

Algorithms

The LU factorization is computed using a variant of Gaussian elimination. Computing an accurate solution is dependent upon the value of the condition number of the original matrix cond(A). If the matrix has a large condition number (it is nearly singular), then the computed factorization might not be accurate.

The LU factorization is a key step in obtaining the inverse with inv and the determinant with det. It is also the basis for the linear equation solution or matrix division obtained with the operators \ and /. This necessarily means that the numerical limitations of lu are also present in these dependent functions.

References

[1] Gilbert, John R., and Tim Peierls. “Sparse Partial Pivoting in Time Proportional to Arithmetic Operations.” SIAM Journal on Scientific and Statistical Computing 9, no. 5 (September 1988): 862–874. https://doi.org/10.1137/0909058.

[2] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[3] Davis, Timothy A. "Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontal method." ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 196–199. https://doi.org/10.1145/992200.992206.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• Code generation does not support sparse matrix inputs for this function.

GPU Code Generation

Generate CUDA® code for NVIDIA® GPUs using GPU Coder™.

Usage notes and limitations:

• Code generation does not support sparse matrix inputs for this function.

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

This function fully supports thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The lu function supports GPU array input with these usage notes and limitations:

• A must be nonsparse.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Distributed Arrays

Partition large arrays across the combined memory of your cluster using Parallel Computing Toolbox™.

Usage notes and limitations:

• A must be nonsparse.

For more information, see Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox).

Version History

Introduced before R2006a