CodePlexProject Hosting for Open Source Software

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:

A_{11} x_{1} + A_{12} x_{2} + A_{13} x_{3} = b_{1}

A_{21} x_{1} + A_{22} x_{2} + A_{23} x_{3} = b_{2}

A_{31} x_{1} + A_{32} x_{2} + A_{33} x_{3} = b_{3}

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:

### 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.

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.

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.

### 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.

A

A

A

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

// 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;

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();

Last edited Apr 2, 2011 at 12:45 AM by ichbin, version 2