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:
If A and b are known, this is a system of three linear equations for the three unknown components of x.
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 = 1.0; b = 2.0; b = 3.0;
// Solve A x = b as x =A^(-1) b
ColumnVector x = A.Inverse() * b;
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
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();
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.