LUDecomposition
Usage Message: LUDecomposition[m] generates a representation of the LU decomposition of a matrix m.
Attributes[LUDecomposition] = {Protected} Related Symbols: LUBackSubstitution
LinearSolve
QRDecomposition
SchurDecomposition
LUFactor
(defined in the LinearAlgebra`GaussianElimination` package)
LUSolve
(defined in the LinearAlgebra`GaussianElimination` package)
Notes: LUDecomposition[m] returns a result {lu, perm, cond}, where lu contains the LU factorization of m, perm is a list of integers that specify the row-permutation of the result, and cond gives an estimate of a condition number of m. The result from LUDecomposition is intended for use in the LUBackSubstitution function. Here is an example showing the use of LUDecomposition and LUBackSubstitution to solve a system of three equations.
In[1]:= m = {{1.2, 4.2, 1.7}, {6.3, 1.7, 5.5}, {3.8, 4.5, 2.3}};
In[2]:= b = {1.6, 4.1, 3.1};
In[3]:= data = LUDecomposition[m]
Out[3]= {{{6.3, 1.7, 5.5}, {0.190476, 3.87619, 0.652381},
> {0.603175, 0.896396, -1.60225}}, {2, 1, 3}, 15.158}
In[4]:= LUBackSubstitution[data, b]
Out[4]= {0.538401, 0.200041, 0.0669103}
Here, for comparison, is the same solution computed using LinearSolve:
In[5]:= LinearSolve[m, b]
Out[5]= {0.538401, 0.200041, 0.0669103}
LUDecomposition is conceptually based on factoring a matrix into a lower-triangular factor L and an upper-triangular factor U such that the matrix product of L and U gives the original matrix. As is typical for most LU decomposition algorithms, the result from LUDecomposition differs from this simple description in two respects. First, to save memory, both L and U are stored in the same matrix. This matrix is returned as the element lu in the result from LUDecomposition. The non-zero elements in L are given by the elements below the diagonal in lu (with diagonal elements equal to 1), and the remaining elements in lu (including the diagonal elements) give the non-zero elements in U. Second, for reasons of numerical stability, the matrix product of L and U gives a row-permuted form of the original matrix. The permutation is specified by a list of integers which is returned as the element perm in the result from LUDecomposition. If the original matrix is m, then the matrix product of L and U is m[[perm]]. Here are two functions, Lower and Upper, which are useful for constructing the L and U matrices from the result that is returned by LUDecomposition.
In[6]:= Lower[m_?MatrixQ] :=
Table[Which[i > j, m[[i,j]], i == j, 1, True, 0],
{i, Length[m]}, {j, Length[m]}]
In[7]:= Upper[m_?MatrixQ] :=
Table[If[i < = j, m[[i, j]], 0], {i, Length[m]}, {j, Length[m]}]
Here is an example to illustrate the relationship between the original matrix and the result from LUDecomposition.
In[8]:= m = {{6,2,9,0},{1,5,7,1},{2,3,2,2},{1,5,2,5}};
In[9]:= {lu, perm, cond} = LUDecomposition[m]; MatrixForm[lu]
Out[9]//MatrixForm= 1 5 7 1
2 -7 -12 0
1 0 -5 4
6 4 -3 6
In[10]:= L = Lower[lu]; MatrixForm[L]
Out[10]//MatrixForm= 1 0 0 0
2 1 0 0
1 0 1 0
6 4 -3 1
In[11]:= U = Upper[lu]; MatrixForm[U]
Out[11]//MatrixForm= 1 5 7 1
0 -7 -12 0
0 0 -5 4
0 0 0 6
In[12]:= MatrixForm[L . U]
Out[12]//MatrixForm= 1 5 7 1
2 3 2 2
1 5 2 5
6 2 9 0
This last result can be seen to be a row-permuted form of the original matrix.
In[13]:= MatrixForm[m]
Out[13]//MatrixForm= 6 2 9 0
1 5 7 1
2 3 2 2
1 5 2 5
The result from L.U is equivalent to m[[perm]]:
In[14]:= perm
Out[14]= {2, 3, 4, 1}
In[15]:= MatrixForm[m[[perm]]]
Out[15]//MatrixForm= 1 5 7 1
2 3 2 2
1 5 2 5
6 2 9 0
Alternatively, the permutation perm can be used to construct a permutation matrix P=IdentityMatrix[Length[m]][[perm]]:
In[16]:= P = IdentityMatrix[Length[m]][[perm]]; MatrixForm[P]
Out[16]//MatrixForm= 0 1 0 0
0 0 1 0
0 0 0 1
1 0 0 0
With this definition of P, L.U is equivalent to P.m, and m is given by Transpose[P].L.U:
In[17]:= P . m == L . U
Out[17]= True
In[18]:= Transpose[P] . L . U == m
Out[18]= True
Known Bugs: LUDecomposition in Version 3.0 can go into an infinite loop for exact matrices in which the first few columns are identical. For example, the following calculation should finish in a fraction of a second, but instead will never finish:
In[1]:= TimeConstrained[
LUDecomposition[{{1, 2, 1}, {1, 2, 2}, {1, 2, 3}}], 60]
Out[1]= $Aborted
You can work around this problem either by using inexact numbers, or by rearranging the columns in the matrix.
Additional Online Documentation:
Mathematica 3.0
http://documents.wolfram.com/v3/RefGuide/LUDecomposition.html
Mathematica 4.0
http://documents.wolfram.com/v4/RefGuide/LUDecomposition.html
Questions or comments? Send email to support@wolfram.com.
| |
