# Solving Systems Of Linear Equations

**Numerics.NET**
supports a simple, uniform interface for solving systems of simultaneous linear equations
and related operations such as computing the inverse of a matrix. All matrix
types and Matrix Decompositions support this interface.

## Linear operators

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

If you want to... | Use this method: |
---|---|

Solve a system of simultaneous linear equations with one right-hand side | |

Solve a system of simultaneous linear equations with multiple right-hand sides. | |

Solve a system of simultaneous linear equations with one right-hand side, optionally overwriting the right-hand side with the solution. | |

Solve a system of simultaneous linear equations with multiple right-hand sides, optionally overwriting the right-hand side with the solution. | |

Solve an over-determined system of simultaneous linear equations with one right-hand side in the least squares sense. | |

Solve an over-determined system of simultaneous linear equations with multiple right-hand sides in the least squares sense. | |

Determine whether a matrix is singular. | |

Calculate the numerical rank of a matrix. | |

Calculate the numerical rank of a matrix using the specified tolerance. | |

Calculate the exact condition number. | |

Estimate the condition number of a matrix. | |

Calculate the inverse of a square matrix. | |

Calculate the Moore-Penrose pseudo-inverse of a matrix. | |

Calculate the determinant of a square matrix. |

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 Matrix<T> 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.CreateFromArray(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();
```

## Solving over-determined systems

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.CreateFromArray(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);
```

## Solving under-determined systems

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.CreateFromArray(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;
```

## A note on performance

The methods defined by the LinearOperator<T> 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 Matrix<T> class uses an The LU Decomposition for its implementation. For badly conditioned matrices, using a The QR Decomposition may be more appropriate.

## References

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