Syntax

Description

The lu function expresses a matrix A as the product of two essentially triangular matrices, one of them a permutation of a lower triangular matrix and the other an upper triangular matrix. The factorization is often called the LU, or sometimes the LR, factorization. A can be rectangular.

Y = lu(A) returns matrix Y that contains the strictly lower triangular L, i.e., without its unit diagonal, and the upper triangular U as submatrices. That is, if [L,U,P] = lu(A), then Y = U+L-eye(size(A)). The permutation matrix P is not returned.

[L,U] = lu(A) returns an upper triangular matrix in U and a permuted lower triangular matrix in L such that A = L*U. Return value L is a product of lower triangular and permutation matrices.

[L,U,P] = lu(A) returns an upper triangular matrix in U, a lower triangular matrix L with a unit diagonal, and a permutation matrix P, such that L*U = P*A. The statement lu(A,'matrix') returns identical output values.

[L,U,P,Q] = lu(A) for sparse nonempty A, returns a unit lower triangular matrix L, an upper triangular matrix U, a row permutation matrix P, and a column reordering matrix Q, so that P*A*Q = L*U. If A is empty or not sparse, lu displays an error message. The statement lu(A,'matrix') returns identical output values.

[L,U,P,Q,R] = lu(A) returns unit lower triangular matrix L, upper triangular matrix U, permutation matrices P and Q, and a diagonal scaling matrix R so that P*(R\A)*Q = L*U for sparse non-empty A. Typically, but not always, the row-scaling leads to a sparser and more stable factorization. The statement lu(A,'matrix') returns identical output values.

[...] = lu(A,'vector') returns the permutation information in two row vectors p and q. You can specify from 1 to 5 outputs. Output p is defined as A(p,:)=L*U, output q is defined as A(p,q)=L*U, and output R is defined as R(:,p)\A(:,q)=L*U.

[...] = lu(A,thresh) controls pivoting. This syntax applies to sparse matrices only. The thresh input is a one- or two-element vector of type single or double that defaults to [0.1, 0.001]. If A is a square matrix with a mostly symmetric structure and mostly nonzero diagonal, MATLAB^® uses a symmetric pivoting strategy. For this strategy, the diagonal where

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

is selected. If the diagonal entry fails this test, a pivot entry below the diagonal is selected, using thresh(1). In this case, L has entries with absolute value 1/min(thresh) or less.

If A is not as described above, MATLAB uses an asymmetric strategy. In this case, the sparsest row i where

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

is selected. A value of 1.0 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 when MATLAB uses an asymmetric strategy.

Smaller values of thresh(1) and thresh(2) 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. The statement lu(A,thresh,'matrix') returns identical output values.

[...] = lu(A,thresh,'vector') controls the pivoting strategy and also returns the permutation information in row vectors, as described above. The thresh input must precede 'vector' in the input argument list.

Note   In rare instances, incorrect factorization results in P*A*Q ≠ L*U. Increase thresh, to a maximum of 1.0 (regular partial pivoting), and try again.

Arguments

Examples

Example 1

Start with

A = [ 1    2    3
      4    5    6
      7    8    0 ];

To see the LU factorization, call lu with two output arguments.

[L1,U] = lu(A)

L1 =
    0.1429    1.0000         0
    0.5714    0.5000    1.0000
    1.0000         0         0

U =
    7.0000    8.0000         0
         0    0.8571    3.0000
         0         0    4.5000

Notice that L1 is a permutation of a lower triangular matrix: if you switch rows 2 and 3, and then switch rows 1 and 2, the resulting matrix is lower triangular and has 1s on the diagonal. Notice also that U is upper triangular. To check that the factorization does its job, compute the product

L1*U

which returns the original A. The inverse of the example matrix, X = inv(A), is actually computed from the inverses of the triangular factors

X = inv(U)*inv(L1)

Using three arguments on the left side to get the permutation matrix as well,

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

returns a truly lower triangular L2, the same value of U, and the permutation matrix P.

L2 =

    1.0000         0         0
    0.1429    1.0000         0
    0.5714    0.5000    1.0000

U =
    7.0000    8.0000         0
         0    0.8571    3.0000
         0         0    4.5000

P =
    0    0    1
    1    0    0
    0    1    0

Note that L2 = P*L1.

P*L1

ans =

    1.0000         0         0
    0.1429    1.0000         0
    0.5714    0.5000    1.0000

To verify that L2*U is a permuted version of A, compute L2*U and subtract it from P*A:

P*A - L2*U

ans =
     0     0     0
     0     0     0
     0     0     0

In this case, inv(U)*inv(L) results in the permutation of inv(A) given by inv(P)*inv(A).

The determinant of the example matrix is

d = det(A)

d =    27

It is computed from the determinants of the triangular factors

d = det(L)*det(U)

The solution to Ax = b is obtained with matrix division

x = A\b

The solution is actually computed by solving two triangular systems

y = L\b
x = U\y

Example 2

The 1-norm of their difference is within roundoff error, indicating that L*U = P*B*Q.

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

B = bucky;

Use the sparse matrix syntax with four outputs to get the row and column permutation matrices.

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

Apply the permutation matrices to B, and subtract the product of the lower and upper triangular matrices.

Z = P*B*Q - L*U;
norm(Z,1)

ans =
  7.9936e-015

Example 3

This example illustrates the benefits of using the 'vector' option. Note how much memory is saved by using the lu(F,'vector') syntax.

F = gallery('uniformdata',[1000 1000],0);
g = sum(F,2);
[L,U,P] = lu(F);
[L,U,p] = lu(F,'vector');
whos P p

Name          Size           Bytes  Class      Attributes
P          1000x1000       8000000  double
p             1x1000          8000  double

The following two statements are equivalent. The first typically requires less time:

x = U\(L\(g(p,:)));
y = U\(L\(P*g));

See Also

cond | det | ilu | inv | qr | rref