# 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.Create(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.Create(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 = lu.ApplyP(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.