The LU Decomposition
The LU decomposition or LU factorization of a matrix expresses the matrix as the product of a lower triangular matrix and an upper triangular matrix. The matrix can be of any type, and of any shape. The LU decomposition is usually written as
A = PLU,
where P is a permutation matrix, L is a lower-triangular matrix, and U is an upper-triangular matrix. If A has more rows than columns, then L is rectangular, and R is square. If A has more columns than rows, then U is rectangular, and R is square. The algorithm uses row pivoting.
The LU decomposition is most commonly used in the solution of systems of simultaneous linear equations. The method is at least twice as fast as other decomposition algorithms. Unfortunately, the method has less desirable numerical problems. For ill-conditioned problems, or when high accuracy is needed, QR would be preferred.
Working with LU Decompositions
The LU decomposition is implemented by the LUDecomposition<T> class. It has no constructors. Instead, it is created by calling the GetLUDecomposition method on the matrix. This method has two overloads. The first overload has no arguments. The second overload takes a Boolean value that specifies whether the contents of the matrix may be overwritten by the LU decomposition. The default is false.
var A = Matrix.CreateFromArray(4, 4, new double[] {
1.80, 5.25, 1.58, -1.11,
2.88,-2.95, -2.69, -0.66,
2.05,-0.95, -2.90, -0.59,
-0.89,-3.80,-1.04, 0.80
}, MatrixElementOrder.ColumnMajor);
var lu = A.GetLUDecomposition();The Decompose method performs the actual decomposition. This method copies the matrix if necessary. It then calls the appropriate LAPACK routine to perform the actual decomposition. This method is called by other methods as needed. You will rarely need to call it explicitly.
Once the decomposition is computed, a number of operations can be performed in much less time. You can repeatedly solve a system of linear equations with different right-hand sides. You can also calculate the determinant and the inverse of the base matrix:
var b = Matrix.CreateFromArray(4, 2, new double[] {
9.52,24.35,0.77,-6.22,
18.47,2.25,-13.28,-6.21
}, MatrixElementOrder.ColumnMajor);
var x = lu.Solve(b);
var invA = lu.GetInverse();
var detA = lu.GetDeterminant();The LowerTriangularFactor and UpperTriangularFactor properties return a TriangularMatrix<T> containing the lower and upper triangular matrices, L and U, of the decomposition:
var L = lu.LowerTriangularFactor;
var U = lu.UpperTriangularFactor;
var P = lu.RowPermutation;The permutation matrix P is not available explicitly. Instead, the RowPermutation property returns the permutation object itself. You can then apply this to a vector or the rows of a matrix:
var Px = P.Apply(x);
x.Rows.PermuteInPlace(P);Band LU Decomposition
The LU decomposition without pivoting of a band matrix is made up of a lower band matrix with lower bandwidth the same as the original matrix and an upper band matrix with upper bandwidth the same as the original matrix. However, pivoting destroys this band structure to a large degree. The bandwidth of the upper triangular matrix is the total bandwidth of the original matrix, and the lower triangular matrix generally doesn't have any band structure, although the number of nonzero elements in each column is limited.
As a result, the factors of the decomposition can't be isolated easily, but the other advantages of the LU decomposition remain. It is much more efficient to repeatedly solve a system of linear equations using the LU decomposition than repeating all the calculations each time.
If a band matrix is to be overwritten with its LU decomposition, the matrix must have enough storage allocated to allow for the upper band to fill up to its maximum width. This can be achieved by setting the allocateForLUDecomposition in the CreateBanded<T>(Int32, Int32, Int32, Int32, Boolean) method to true. The other limitation is that the lower triangular factor of the decomposition is not available.