LU Decomposition

LU Decomposition is a tool for solving linear systems. Given a square matrix A, it finds matrices L and U such that A = L U, where L is lower-left triangular and U is upper-right triangular.

With an LU decomposition in hand, it's easy to solve Ax = b; just write it as L(Ux) = b. Use back-substitution to solve Ly = b and then use back-substitution again to solve Ux=y. Solving a linear system in this way is both faster and much more stable than finding A^-1 and then multiplying out A^-1b. LU decomposition is the fastest known technique for solving a linear system.

For the examples that follow, we will use the matrix initialized by the following code:

using Meta.Numerics.Matrices;

SquareMatrix A = new SquareMatrix(new double[,] {
    { 1, -2, 3 },
    { 2, -5, 12 },
    { 0, 2, -10 }
});

We will also use the PrintMatrix method defined earlier.

How do I LU Decompose a matrix?

Easy:

LUDecomposition lud = A.LUDecomposition();

Here's some code to prove that our decomposition actually does its job:

PrintMatrix("LU", lud.LMatrix() * lud.UMatrix());
PrintMatrix("PA", lud.PMatrix() * A);

Notice that our LU isn't actually a decomposition of A, but rather the decomposition of a row-wise permutation of A. This slight complication is a necessary consequence of the partial pivoting that is used to make LU decomposition stable. Since permutations are trivially invert-able, this complication makes our LU decomposition slightly less pretty, but no less useful or performant.

How do I solve Ax = b?

Also easy:

ColumnVector b = new ColumnVector(2, 8, -4);
ColumnVector x = lud.Solve(b);
PrintMatrix("x", x);

Here's some code to prove that the solution actually works:

PrintMatrix("Ax", A * x);

Can I get det A and A^-1?

Yes. Here is some code to print the determinant and show that A-inverse does what's it's supposed to:

Console.WriteLine($"det(a) = {lud.Determinant()}");
PrintMatrix("AI A", lud.Inverse() * A);

Indulge us in a reminder that it's both faster and more accurate to solve Ax=b by using the Solve() method than by calling Inverse() to get AI and then multiplying AI * b.

What other decompositions are possible?

Meta.Numerics can also do QR Decomposition and Singular Value Decomposition. It can also do eigenvalue decompositions, but those are not useful for solving linear systems.