The Extreme Optimization Numerical Libraries for .NET
support 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.
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:
If you want to... | Use this method: |
---|
Solve a system of simultaneous linear equations with one right-hand side | Solve(VectorT) |
Solve a system of simultaneous linear equations with multiple right-hand sides. | Solve(MatrixT) |
Solve a system of simultaneous linear equations with one right-hand side,
optionally overwriting the right-hand side with the solution.
| Solve(DenseVectorT, Boolean) |
Solve a system of simultaneous linear equations with multiple right-hand sides, optionally overwriting the
right-hand side with the solution.
| Solve(DenseMatrixT, Boolean) |
Solve an over-determined system of simultaneous linear equations
with one right-hand side in the least squares sense. | LeastSquaresSolve(VectorT) |
Solve an over-determined system of simultaneous linear equations
with multiple right-hand sides in the least squares sense.
| LeastSquaresSolve(MatrixT) |
Determine whether a matrix is singular. | IsSingular |
Calculate the numerical rank of a matrix. | Rank |
Calculate the numerical rank of a matrix using the specified tolerance. | Rank(T) |
Calculate the exact condition number. | GetConditionNumber |
Estimate the condition number of a matrix. | EstimateConditionNumber |
Calculate the inverse of a square matrix. | GetInverse |
Calculate the Moore-Penrose pseudo-inverse of a matrix. | GetPseudoInverse |
Calculate the determinant of a square matrix. | GetDeterminant |
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);
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)
Dim x = A.Solve(b)
Dim invA = A.GetInverse()
Dim cond = A.GetConditionNumber()
Dim condEstimate = A.EstimateConditionNumber()
No code example is currently available or this language may not be supported.
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)
let x = A.Solve(b)
let invA = A.GetInverse()
let cond = A.GetConditionNumber()
let 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.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)
No code example is currently available or this language may not be supported.
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)
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.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
No code example is currently available or this language may not be supported.
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.