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
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 |
A | Rectangular matrix to be factored. |
thresh | Pivot threshold for sparse matrices. Valid values are
in the interval |
L | Factor of |
U | Upper triangular matrix that is a factor of |
P | Row permutation matrix satisfying the equation |
Q | Column permutation matrix satisfying the equation |
R | Row-scaling matrix |
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 1
s
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
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));