Solving Systems Of Linear Equations

All the functionality for solving systems of simultaneous linear equations and related operations are defined by the LinearOperatorT abstract base class. Both matrices and matrix decompositions inherit from this class. The operations defined by this class are listed in the following table:

In the example below, the matrix A of the system is first created. Note that matrix elements are stored in column major order by default. The first three components are the coefficients of x, followed by the coefficients for y and z. The vector b contains the right-hand sides of the equations.

The system is solved simply by calling the MatrixT method of the matrix A, with b as its only argument. The results are then printed.

As a little add-on, the inverse of A is calculated. The result of multiplying a matrix with its inverse must be the identity matrix. This value is printed as well as an estimate for the condition number of A.

var A = Matrix.Create(3, 3, new double[]
    { 1, 1, 1,
      3,-2, 3,
      1, 2,-1 }, MatrixElementOrder.RowMajor);
var b = Vector.Create(0.0, 10.0, 2.0);
// Use the Solve method on a to get the solution:
var x = A.Solve(b);
var invA = A.GetInverse();

var cond = A.GetConditionNumber();
var condEstimate = A.EstimateConditionNumber();

Dim A = Matrix.Create(3, 3, New Double() {
    1, 1, 1,
    3, -2, 3,
    1, 2, -1}, MatrixElementOrder.RowMajor)
Dim b = Vector.Create(0.0, 10.0, 2.0)
' Use the Solve method on a to get the solution:
Dim x = A.Solve(b)
Dim invA = A.GetInverse()

Dim cond = A.GetConditionNumber()
Dim condEstimate = A.EstimateConditionNumber()

let A = Matrix.Create(3, 3, 
          [| 1.0; 1.0; 1.0;
             3.0;-2.0; 3.0;
             1.0; 2.0;-1.0 |], MatrixElementOrder.RowMajor)
let b = Vector.Create(0.0, 10.0, 2.0)
// Use the Solve method on a to get the solution:
let x = A.Solve(b)
let invA = A.GetInverse()

let cond = A.GetConditionNumber()
let condEstimate = A.EstimateConditionNumber()

An over-determined system is a system with more equations than unknowns. Solving an over-determined system corresponds to solving a least squares problem. The simplest way to find the least-squares solution is to call the LeastSquaresSolve method or one of its overloads. This method can be called directly on a matrix for a one-off solution. Several decomposition classes also support solving least squares problems, most commonly QR or the SVD.

var A = Matrix.Create(6, 4, new double[] {
    1, 1, 1, 1, 1, 1,
    1, 2, 3, 4, 5, 6,
    1, 4, 9, 16, 25, 36,
    1, 2, 1, 2, 1, 2 }, MatrixElementOrder.ColumnMajor);
var b = Vector.Create(1.0, 3.0, 6.0, 11.0, 15.0, 21.0);
var x = A.LeastSquaresSolve(b);

Dim A = Matrix.Create(6, 4, New Double() {
    1, 1, 1, 1, 1, 1,
    1, 2, 3, 4, 5, 6,
    1, 4, 9, 16, 25, 36,
    1, 2, 1, 2, 1, 2}, MatrixElementOrder.ColumnMajor)
Dim b = Vector.Create(1.0, 3.0, 6.0, 11.0, 15.0, 21.0)
Dim x = A.LeastSquaresSolve(b)

let A = Matrix.Create(6, 4, 
          [|
            1.0; 1.0; 1.0; 1.0; 1.0; 1.0;
            1.0; 2.0; 3.0; 4.0; 5.0; 6.0;
            1.0; 4.0; 9.0; 16.0; 25.0; 36.0;
            1.0; 2.0; 1.0; 2.0; 1.0; 2.0 
          |], MatrixElementOrder.ColumnMajor)
let b = Vector.Create(1.0, 3.0, 6.0, 11.0, 15.0, 21.0)
let x = A.LeastSquaresSolve(b)

An under-determined system is a system with less equations than unknowns. There are no methods that directly solve an under-determined system. However, solving such a system is very easy using the pseudo-inverse. Multiplying the right-hand side of a system of equations by the pseudo-inverse of the matrix of the system produces a solution to the system of equations. Moreover, this is the solution with smallest norm.

The example below solves a system of 2 equations in three unknowns:

var A = Matrix.Create(2, 3, new double[]
    {1, 3, 1,
     1, -2, 2}, MatrixElementOrder.RowMajor);
var b = Vector.Create(0.0, 10.0);
var Aplus = A.GetPseudoInverse();
var x = Aplus * b;

Dim A = Matrix.Create(2, 3, New Double() {
     1, 3, 1,
     1, -2, 2}, MatrixElementOrder.RowMajor)
Dim b = Vector.Create(0.0, 10.0)
Dim Aplus = A.GetPseudoInverse()
Dim x = Aplus * b

let A = Matrix.Create(2, 3, 
            [|
                 1.0; 3.0; 1.0;
                 1.0;-2.0; 2.0
            |], MatrixElementOrder.RowMajor)
let b = Vector.Create(0.0, 10.0)
let Aplus = A.GetPseudoInverse()
let x = Aplus * b

The methods defined by the LinearOperatorT class are connected because they all rely on some shared preliminary calculations. Nearly always, some kind of matrix decomposition is performed first. The classes that represent Matrix Decompositions have been optimized to take advantage of this situation. Matrix classes, on the other hand, always perform all calculations from scratch.

As a result, it is much more desirable from a performance perspective to first create a matrix decomposition and call these methods on the decomposition rather than on the matrix directly. This also gives you more control on how the calculations are performed. For instance, the MatrixT class uses an The LU Decomposition for its implementation. For badly conditioned matrices, using a The QR Decomposition may be more appropriate.

G. H. Golub, C. F. Van Loan, Matrix Computations (3rd Ed), Johns Hopkins University Press, 1996.