QR Decomposition

QR Decomposition is a useful tool for solving linear systems. Given the matrix A, it writes A = Q * R, where Q is orthogonal (Q^-1 = Q^T) and R is upper-right-triangular.

With a QR decomposition, it's easy to solve Ax = b: just write it as R x = Q^T b, find Q^T b by multiplication, and solve for x by back-substitution. QR decomposition is a little slower than LU decomposition, but is often a little more stable. It's my go-to all-round technique for solving linear systems.

Here's code to create a sample matrix A that we will use in the following examples:

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 QR Decompose a Matrix?

Easy:

SquareQRDecomposition qrd = A.QRDecomposition();

QR Decomposition is an O(N³) operation, so if your want to solve large systems, you should investigate how your performance scales with problem size.

Note that you can also QR decompose a non-square matrix.

How do I get a solution to Ax = b?

Also easy:

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

Here's some more code to prove the solution:

PrintMatrix("A x", A * x);

How to I get Q and R?

They are available as properties of the QRDecomposition. Here's some code to prove that Q is orthogonal and that the QR decomposition works:

SquareMatrix Q = qrd.QMatrix;
SquareMatrix R = qrd.RMatrix;
PrintMatrix("Q  Q^T", Q * Q.Transpose());
PrintMatrix("QR", Q * R);

Note that QMatrix and RMatrix are properties that return read-only matrices, which makes them fast operations. If you want to modify one of these matrices, just call Copy() to get a write-able copy.

Can I get A^-1?

Yes! Here's some code to get it and prove it does its job:

SquareMatrix AI = qrd.Inverse();
PrintMatrix("A^{-1}", AI);
PrintMatrix("A^{-1} A", AI * A);

Indulge us in another 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 if my matrix is singular?

If A is exactly singular, QRDecompose() will fail with a DivideByZeroException. If A is nearly singular, QRDecompose() can succeed but give solutions dominated by round-off errors. Singular and nearly singular systems can be analyzed by Singular Value Decomposition.