LU matrix factorization
Y = lu(A)
[L,U] = lu(A)
[L,U,P] = lu(A)
[L,U,P,Q] = lu(A)
[L,U,P,Q,R] = lu(A)
[...] = lu(A,'vector')
[...] = lu(A,thresh)
[...] = lu(A,thresh,'vector')
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, for sparse A, 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)). For nonsparse A, Y is the output from the LAPACK dgetrf or zgetrf routine. The permutation matrix P is not returned.
[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.
Rectangular matrix to be factored.
Pivot threshold for sparse matrices. Valid values are in the interval [0,1]. If you specify the fourth output Q, the default is 0.1. Otherwise, the default is 1.0.
Factor of A. Depending on the form of the function, L is either a unit lower triangular matrix, or else the product of a unit lower triangular matrix with P'.
Upper triangular matrix that is a factor of A.
Row permutation matrix satisfying the equation L*U = P*A, or L*U = P*A*Q. Used for numerical stability.
Column permutation matrix satisfying the equation P*A*Q = L*U. Used to reduce fill-in in the sparse case.
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
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
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
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));