Extreme Optimization > Mathematics Library for .NET > User's Guide > Current Page > Solving Systems Of Linear Equations

Extreme Optimization Mathematics Library for .NET

Solving Systems Of Linear Equations

The Extreme Optimization Mathematics Library for .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 decompositions support this interface.

The LinearTransformation class

All the functionality for solving systems of simultaneous linear equations and related operations are defined by the LinearTransformation abstract base 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(Vector)
Solve a system of simultaneous linear equations with multiple right-hand sides.	Solve(Matrix)
Solve a system of simultaneous linear equations with one right-hand side, optionally overwriting the right-hand side with the solution.	Solve(Vector, Boolean)
Solve a system of simultaneous linear equations with multiple right-hand sides, optionally overwriting the right-hand side with the solution.	Solve(Matrix, Boolean)
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(Double)
Calculate the exact condition number.	GetConditionNumber()
Estimate the condition number of a matrix.	EstimateConditionNumber()
Calculate theinverse of a square matrix.	Inverse()
Calculate the Moore-Penrose pseudo-inverse of a matrix.	GetPseudoInverse()
Calculate the determinant of a square matrix.	Determinant()

The Solve method is overloaded. To solve for one right-hand side, pass this value as a Vector as the first parameter. To solve for multiple right-hand sides, store the vectors inthe columns of a GeneralMatrix. An optional second parameter specifies whether the right-hand sides should be overwritten with the solution. By default, a new vector or matrix instance is created of the same dimensions as the right-hand side.

A matrix is singular when not all of its rows or columns are linearly independent. The IsSingular method returns a Boolean value that indicates this condition.

The inverse matrix is returned by the GetInverse method. The inverse matrix is only defined for square matrices. If this method is applied to a rectangular matrix, an exception of type DimensionMismatchException is thrown.

The condition number of a matrix is defined as the ratio of its largest to its smallest singular value. Because the calculation of singular values is a very expensive operation, an estimate that is cheaper to calculate is usually preferred. The EstimateConditionNumber method returns such an estimate.

The condition number gives an indication ofthe worst case loss of precision when solving a system of simultaneous linear equations. The condition number of a singular matrix is infinite, which is returned as Double.PositiveInfinity. Note that, due to round-off error, most singular matrices will not appear to be exactly singular. For these matrices, the condition number is very large, of the order of 1/MachineConstants.Epsilon.

The determinant of a matrix is returned by theGetDeterminant method. The determinant is only defined for square matrices. If this method is applied to a rectangular matrix, an exception of type DimensionMismatchException is thrown.

An example

The following example uses the Solve method of GeneralMatrix to solve a system of three linear equations:

// Remember matrices are stored in column major order:
GeneralMatrix A = new GeneralMatrix(3, 3, new double[]
    {1, 3, 1,
     1, -2, 2,
     1, 3, -1}
Vector b = new Vector(0.0, 10.0, 2.0);
// Use the Solve method on a to get the solution:
Vector x = A.Solve(b);
// Print the results:
Console.WriteLine("x = {0}", x[0]);
Console.WriteLine("y = {0}", x[1]);
Console.WriteLine("z = {0}", x[2]);
// Now show the inverse is correctly calculated:
GeneralMatrix invA = A.Inverse();
Console.WriteLine("A * A^-1 = {0}", A * invA);
// Print an estimate for the condition number
Console.WriteLine("cond(A) = ", A.ConditionNumber());

' Remember matrices are stored in column major order:
Dim A As GeneralMatrix = New GeneralMatrix(3, 3, New Double() _
    {1, 3, 1, _
     1, -2, 2, _
     1, 3, -1}
Dim b As Vector = new Vector(0.0, 10.0, 2.0)
' Use the Solve method on a to get the solution:
Dim x As Vector = A.Solve(b)
' Print the results:
Console.WriteLine("x = {0}", x(0))
Console.WriteLine("y = {0}", x(1))
Console.WriteLine("z = {0}", x(2))
' Now show the inverse is correctly calculated:
Dim invA As GeneralMatrix = A.Inverse()
Console.WriteLine("A * A^-1 = {0}", A * invA)
' Print an estimate for the condition number
Console.WriteLine("cond(A) = ", A.ConditionNumber())

In the above example, 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 Solve method of the matrix A, with b as its only parameter. 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.

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:

// Remember matrices are stored in column major order:
GeneralMatrix A = new GeneralMatrix(3, 3, new double[]
    {1, 3, 1,
     1, -2, 2});
Vector b = new Vector(0.0, 10.0);
// Calculate the pseudo-inverse of A:
Matrix Aplus = A.GetPseudoInverse();
// And multiply b to get the solution:
Vector x = Aplus * b;
// Print the results:
Console.WriteLine("x = {0}", x[0]);
Console.WriteLine("y = {0}", x[1]);
Console.WriteLine("z = {0}", x[2]);

' Remember matrices are stored in column major order:
Dim A As GeneralMatrix = New GeneralMatrix(3, 3, New Double() _
    {1, 3, 1, _
     1, -2, 2})
Dim b As Vector = new Vector(0.0, 10.0)
' Use the Solve method on a to get the solution:
Dim Aplus As Matrix = A.GetPseudoInverse()
' And multiply by b to get the solution:
Dim x As Vector = Matrix.Multiply(Aplus, b)
' Print the results:
Console.WriteLine("x = {0}", x(0))
Console.WriteLine("y = {0}", x(1))
Console.WriteLine("z = {0}", x(2))

Solving an over-determined system corresponds to solving a least squares problem, which is covered in the next section.

A note on performance

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

References

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

Copyright 2004-2008, Extreme Optimization. All rights reserved.
Extreme Optimization, Complexity made simple, M#, and M Sharp  are trademarks of ExoAnalytics Inc.
Microsoft, Visual C#, Visual Basic, Visual Studio, Visual Studio.NET, and the Visual Studio Logo are registered trademarks of Microsoft Corporation