Solving Sparse Systems
The most common applications of large, sparse matrices involve the solution of systems of equations. There are two basic methods:
A direct method works the same as solving dense systems: the first step is to decompose the matrix into simpler matrices which can then be used to solve the system more quickly.
An iterative method works like algorithms for solving nonlinear equations: an approximate solution is improved iteratively until it is close enough to the actual solution by some measure.
Both methods have their place. Direct methods are generally used for smaller dense systems, and where the same system has to be solved with different right-hand sides. Iterative methods really stand out for large systems, where they can be orders of magnitude faster than direct methods.
Direct sparse solvers
A direct solver computes the decomposition of the matrix in order to calculate the solution directly in one step. Direct methods are most useful small to medium-size sparse matrices.
The LU decomposition of a matrix expresses the matrix as the product of a lower triangular matrix and an upper triangular matrix. The matrix can be of any type, and of any shape.The LU decomposition is usually written as
AQ = PLU,
where P and Q are permutation matrices, L is a lower-triangular matrix, and U is an upper-triangular matrix. If A has more rows than columns, then L is rectangular, and R is square. If A has more columns than rows, then U is rectangular, and R is square. The algorithm uses row pivoting to improve stability.
For sparse matrices, many zero elements will become nonzero. This is unavoidable. However, it can be mitigated by applying a suitable column permutation. The sparse LU decomposition uses an approximate minimum degree (AMD) ordering. The matrix Q is not available explicitly. Instead, the ColumnPermutation property returns the Permutation object, which can then be used to apply the permutation, or its inverse, directly to vectors or matrices.
Iterative sparse solvers
Iterative sparse solvers produce a sequence of approximations that converges to the solution of the linear system. Iterative solvers have several advantages over direct solvers. First, there is no need to explicitly form the matrix. All that is needed is a function to multiply a vector by the matrix. This makes it possible to solve very large systems with only relatively modest memory requirements. Second, preconditioners may be used to accelerate convergence. See below for more details.
The classes that implement iterative solvers live in the Numerics.NET.LinearAlgebra.IterativeSolvers namespace.
Types of iterative solvers
There are two basic types of iterative solvers. So-called stationary methods use an iteration formula that is independent of the iteration. Stationary methods do not use a preconditioner. non-stationary methods use information that changes at each iteration. In general, non-stationary methods are more effective. Another aspect that distinguishes the algorithms is whether they are suitable for non-symmetric systems. The different algorithms that are available are listed below.
Class | Description |
|---|---|
Implements simple Jacobi iteration. | |
Implements Gauss-Seidel iteration. | |
Implements Successive Over-Relaxation. Requires a relaxation factor between 0 and 2. |
Class | Description |
|---|---|
Implements the Conjugate Gradient (CG) method, suitable for symmetric systems. Convergence can be erratic, and may sometimes break down completely. | |
Implements the Conjugate Gradient Squared (CGS) method, suitable for symmetric and non-symmetric systems. Convergence can be erratic, and may sometimes break down completely. | |
Implements the Bi-Conjugate Gradient (BiCG) method, suitable for symmetric and non-symmetric systems. For symmetric systems, it takes twice as long as the standard Conjugate Gradient method. | |
Implements the Bi-Conjugate Gradient Stabalized (BiCGStab) method, suitable for symmetric and non-symmetric systems. Convergence is more stable than the Conjugate Gradient Squared method. | |
Implements the Quasi-Minimal Residual (QMR) method, suitable for symmetric and non-symmetric systems. This method is more robust than the Bi-Conjugate Gradient method. | |
Implements the Generalized Minimal Residual (GMRES) method, suitable for symmetric and non-symmetric systems. |
Preconditioners
A preconditioner is a transformation of a matrix that puts it in a form with better convergence properties. A preconditioner takes the form of a matrix that approximates the matrix of the system, but that is easier to solve than the original system. Preconditioners inherit from Preconditioner<T> and must implement the Solve method and in some cases SolveTranspose.
Several preconditioners have been predefined. The classes reside in the Numerics.NET.LinearAlgebra.IterativeSolvers.Preconditioners namespace. They are as follows:
Class | Description |
|---|---|
Implements the identity preconditioner. It has no effect and is equivalent to not using a preconditioner at all. | |
Implements the Jacobi preconditioner, which approximates the matrix by its diagonal. This is the default preconditioner for matrices whose elements can be individually accessed. | |
Implements an ILU(0) preconditioner. The matrix is approximated by an LU decomposition where during the decomposition process entries that are not present in the original matrix are ignored. This is often the most effective but also the most expensive of the pre-defined preconditioners. |
Solving equations
All classes that implement iterative sparse solvers inherit from the IterativeSparseSolver<T> class. They take the matrix of the linear system as the first parameter in the constructor(s). The matrix must be supplied as an object that inherits from LinearOperator<T>, and requires only that the Multiply method be implemented.
Once the solver object is constructed, a preconditioner can optionally be set using the Preconditioner property. The right-hand side may be set through the RightHandSide property. An initial guess can be supplied by setting the InitialGuess property.
The Solve method performs the actual calculation. It has three overloads. The first overload takes no arguments and solves the system using the current values of the RightHandSide and InitialGuess properties. The second overload takes one argument which specifies the right-hand side. The third overload takes two arguments. The first argument is once again the right-hand side. The second parameter is a corresponding initial guess. An exception is thrown if the length of the right-hand side or the initial guess does not match the size of the matrix.
Convergence is determined by a combination of two tests. The SolutionTest succeeds when the difference between two approximations is less than the tolerance. The ResidualTest succeeds when the residuals are below the tolerance. The properties of these convergence test objects can be modified to fine tune the convergence criteria.
If the algorithm fails to converge, no exception is thrown unless the ThrowExceptionOnFailure property is set to true. The Status property indicates how the algorithm terminated. Its relevant values and their meaning are listed below.
Value | Description |
|---|---|
NoResult | The algorithm has not been executed. |
Converged | The algorithm ended normally. The desired accuracy has been achieved. |
IterationLimitExceeded | The number of iterations needed to achieve the desired accuracy is greater than MaxIterations. |
RoundOffError | Round-off prevented the algorithm from achieving the desired accuracy. |
It is a good idea to always check the value of the Status property after calling Solve. The Residuals property returns the residuals.
The following example creates a BiCG solver for a non-symmetric matrix, selects an incomplete LU preconditioner, solves the system and prints the number of iterations and the estimated error:
var solver = new BiConjugateGradientSolver<double>(A);
solver.Preconditioner = new IncompleteLUPreconditioner<double>(A);
solver.Solve(b);
Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded);
Console.WriteLine("Estimated error: {0}", solver.SolutionReport.Error);