## Documentation Center |

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] = 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.

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 *A**x* = *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));

Was this topic helpful?