lu - LU matrix factorization - MATLAB (original) (raw)

Syntax

Description

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

example

[[L](#d126e1027346),[U](#d126e1027389),[P](#d126e1027406)] = lu([A](#d126e1027033)) 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.

example

[[L](#d126e1027346),[U](#d126e1027389),[P](#d126e1027406)] = lu([A](#d126e1027033),[outputForm](#mw%5Fc7ac10f5-4b1d-4b24-aac5-cb08fc302896)) returns P in the form specified byoutputForm. Specify outputForm as'vector' to return P as a permutation vector such that A(P,:) = L*U.

example

[[L](#d126e1027346),[U](#d126e1027389),[P](#d126e1027406),[Q](#d126e1027434)] = lu([S](#d126e1027066)) 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 matrixQ, such that P*S*Q = L*U.

example

[[L](#d126e1027346),[U](#d126e1027389),[P](#d126e1027406),[Q](#d126e1027434),[D](#d126e1027462)] = lu([S](#d126e1027066)) also returns a diagonal scaling matrix D such thatP*(D\S)*Q = L*U. Typically, the row-scaling leads to a sparser and more stable factorization.

[___] = lu([S](#d126e1027066),[thresh](#d126e1027099)) specifies thresholds for the pivoting strategy employed bylu 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](#mw%5Fc7ac10f5-4b1d-4b24-aac5-cb08fc302896)) returns P and Q in the form specified byoutputForm. Specify outputForm as'vector' to return P andQ as permutation vectors. You can use any of the input argument combinations in previous syntaxes.

example

Examples

collapse all

Compute the LU factorization of a matrix and examine the resulting factors. LU factorization is a way of decomposing a matrix A into an upper triangular matrix U, a lower triangular matrix L, and a permutation matrix P such that PA=LU. These matrices describe the steps needed to perform Gaussian elimination on the matrix until it is in reduced row echelon form. The L matrix contains all of the multipliers, and the permutation matrix 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 = 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.

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

ans = 3×3

10    -7     0
-3     2     6
 5    -1     5

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 Ax=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

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:

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

x = 5×1

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

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.

Compute the LU factorization of S using the sparse matrix syntax with four outputs to return the row and column permutation matrices.

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 = PSQ - L*U; norm(e,1)

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.

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

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

Figure contains 2 axes objects. Axes object 1 with title L factor, xlabel nz = 10351 contains a line object which displays its values using only markers. Axes object 2 with title U factor, xlabel nz = 6046 contains a line object which displays its values using only markers.

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

Figure contains 2 axes objects. Axes object 1 with title L factor, xlabel nz = 1854 contains a line object which displays its values using only markers. Axes object 2 with title U factor, xlabel nz = 2391 contains a line object which displays its values using only markers.

Input Arguments

collapse all

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

Data Types: single | double
Complex Number Support: Yes

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

Data Types: single | double
Complex Number Support: Yes

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 ofthresh 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) to1.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 — IfS is a square sparse matrix with a mostly symmetric structure and mostly nonzero diagonal, thenlu uses a symmetric pivoting strategy. For this strategy, the algorithm selects the diagonal pivotj 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 inL have an absolute value of1/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 inP*S*Q ≠ L*U. If this occurs, increase thresh to a maximum of1.0 (regular partial pivoting), and try again.

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

As matrices, the outputsP 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, andD satisfy P*(D\S)*Q = L*U.

As vectors, the outputsP 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, andD satisfy D(:,P)\S(:,Q) = L*U.

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

Output Arguments

collapse all

Lower triangular factor, returned as a matrix. The form ofL depends on whether the row permutationsP are returned in a separate output:

If the third output P is specified, thenL 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 productP'*L of the outputs P and L returned in the three output case.

Upper triangular factor, returned as an upper triangular 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.

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.

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 withinv 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

expand all

Usage notes and limitations:

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

Usage notes and limitations:

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

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

Version History

Introduced before R2006a

expand all

You can specify the sparse input matrix S as single precision. If S is single, the output arguments are also single except those that are related to indexing, such as ordering and permutation vectors.