Matrices are commonly used to represent systems of linear equations. For example, A x = b, written out in terms of components is the system of linear equations:

A11 x1 + A12 x2 + A13 x3 = b1
A21 x1 + A22 x2 + A23 x3 = b2
A31 x1 + A32 x2 + A33 x3 = b3

If A and b are known, this is a system of three linear equations for the three unknown components of x.

Matrix Inverse

Naively, the matrix equation can be solved by using the inverse of A to write x = A-1 b. Meta.Numerics can certainly do this for you:

// Define a square matrix
SquareMatrix A = new SquareMatrix(3);
A[0,0] = 1.0; A[0,1] = 2.0; A[0,2] = 3.0;
A[1,0] = 4.0; A[1,1] = 5.0; A[1,2] = 6.0;
A[2,0] = 7.0; A[2,1] = 8.0; A[2,2] = 9.0;
// Define a column vector
ColumnVector b = new ColumnVector(3);
b[0] = 1.0; b[1] = 2.0; b[2] = 3.0;
// Solve A x = b as x =A^(-1) b
ColumnVector x = A.Inverse() * b;

Matrix Decomposition

Finding a matrix inverse is almost never the best way to solve a system of equations, though. Instead, you should form one of any number of matrix decompositions, such as an LU decomposition or QR decomposition. Such decompositions represent the original matrix as a product of matrices with special forms that make solving sytems of equations very efficient. Using a decomposition to solve a system of equations is usually faster and more accurate that finding the matrix inverse and multiplying it by the right-hand-side vector.

The following code solves the same system of equation using LU decomposition, which the fastest known decomposition of an arbitrary squre matrix for this purpose.

// Decompose A = L U
SquareLUDecomposition LU = A.LUDecomposition();
// Solve A x = b via x = U^(-1) L^(-1) b
ColumnVector x= LU.Solve(b);

As with the matrix inverse, once you have obtained a decomposition, you can use it to solve for any many different right-hand-side vectors.

// Solve A y = c using existing LU decomposition
ColumnVector y = LU.Solve(c);

The decomposition classes not only allow you to solve the system for different right-hand-side vectors, but also to quickly determine other properties of the original matrix, such as its determinant and, if you want, its matrix inverse.

SquareMatrix AI = LU.Inverse();

Cholesky Decomposition

If A has some special symmetry, such a being symmetric or tri-diagonal, using the matrix class for that symmetry will provide you with methods of decomposition optimized for that kind of matrix. For example, the Cholesky decomposition is a decomposition that exists only for symmetric, positive definite matrices. It is provided by the CholeskyDecomposition() method of the SymmetricMatrix class, and is about twice as fast as LU decomposition.

Last edited Apr 1, 2011 at 11:45 PM by ichbin, version 2


No comments yet.