Solving Sparse Systems

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 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 Extreme.Mathematics.LinearAlgebra.IterativeSolvers namespace.

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.

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 PreconditionerT and must implement the Solve method and in some cases SolveTranspose.

Several preconditioners have been predefined. The classes reside in the Extreme.Mathematics.LinearAlgebra.IterativeSolvers.Preconditioners namespace. They are as follows:

All classes that implement iterative sparse solvers inherit from the IterativeSparseSolverT 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 LinearOperatorT, 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.

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);

Dim solver = New BiConjugateGradientSolver(Of Double)(A)
solver.Preconditioner = New IncompleteLUPreconditioner(Of Double)(A)
solver.Solve(b)
Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded)
Console.WriteLine("Estimated error: {0}", solver.SolutionReport.Error)

let solver = new BiConjugateGradientSolver<float>(A)
solver.Preconditioner <- new IncompleteLUPreconditioner<float>(A)
let x = solver.Solve(b)
printfn "Solved in %d iterations." solver.IterationsNeeded
printfn "Estimated error: %f" solver.SolutionReport.Error