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.

Choose between direct and iterative methods according to matrix structure, factor fill, available memory, accuracy and the number of right-hand sides. Preconditioner setup can be amortized over repeated solves; iterative convergence depends on both the operator and the preconditioner.

A direct solver computes the decomposition of the matrix in order to calculate the solution directly in one step. Direct methods are useful for small to medium-size sparse matrices and repeated right-hand sides.

The LU decomposition of a matrix expresses the matrix as the product of a lower triangular matrix and an upper triangular matrix. The matrix may be rectangular. 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 U is square. If A has more columns than rows, then U is rectangular, and L 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 solvers use matrix-vector products and can accept a matrix-free LinearOperator<T>. Some methods additionally require transpose or conjugate-transpose products. Choose the method from the operator's mathematical properties and the preconditioner's action, then check the residual of the original system.

For complex systems, Hermitian means equal to the conjugate transpose. Symmetry alone does not make a matrix suitable for conjugate gradient. Positive definiteness is an additional requirement.

IterationsBetweenRestarts and IterationsBetweenRestarts control restart length. Both default to 30. Restarting limits storage but can slow or prevent convergence. MINRES and FGMRES do not apply the inherited LeftPreconditioner property.

Stationary alternatives are JacobiSolver<T>, GaussSeidelSolver<T> and SuccessiveOverRelaxationSolver<T>. Their convergence depends on the matrix and, for SOR, the relaxation parameter. They do not use the Krylov preconditioner property.

A preconditioner approximates the operator with a system that is cheaper to solve. Stronger factors can reduce iterations while increasing setup, storage and application costs. No choice is uniformly best; measure total cost for the matrices and right-hand sides used by your application.

Matrix-free operators use identity preconditioning by default. A user-supplied preconditioner derives from Preconditioner<T> and implements the required solve actions. BiCG and QMR require the corresponding transpose/adjoint actions as well. FGMRES permits varying actions; CG, MINRES and ordinary GMRES require the action to remain fixed during each solve.

Incomplete LU and Cholesky offer explicit relative pivot protection, diagonal shifts and equilibration. A shift changes only the preconditioner target; the solver still solves the original system. Incomplete Cholesky, including its modified and threshold variants, may break down even on an SPD matrix. Successful factorization does not certify that the original operator is positive definite. Choose a shift, ordering or different preconditioner when factorization fails.

ILUT's fill cap applies separately to strict lower and upper rows. Threshold IC's cap applies to strict lower columns. Lower drop tolerances and larger caps can improve an approximation at a greater cost. Factor nonzero counts, factorization counts and minimum relative pivots help assess the resulting factors; threshold IC also reports discarded candidates.

GetReverseCuthillMcKeeOrdering() returns an explicit reverse Cuthill-McKee permutation. It can reduce bandwidth but does not promise reduced factor fill or faster convergence. A permutation maps factor coordinates to original coordinates and is copied by the factor constructor.

Incomplete-factor constructors accept ordering and equilibration settings. The factor handles these transformations internally. Keep the original matrix and right-hand side on the solver; solutions and residual checks remain in the original coordinates.

Incomplete Cholesky uses matching row and column scaling to preserve Hermitian structure. ILU and ILUT offer row equilibration for general systems. These choices can help a factorization, but do not guarantee improved conditioning or convergence. CG and MINRES still require a fixed, Hermitian positive definite preconditioner.

Create a solver with a square matrix or linear operator, set Preconditioner if needed, and call Solve(Vector<T>). The overload with an initial guess starts from that vector; otherwise the default guess is zero. Set the inherited maximum iteration count and tolerance for the application's accuracy requirements.

ResidualTest is enabled by default and SolutionTest is disabled. The residual test measures the original system. A small update alone does not establish an accurate solution. Configure absolute and relative tolerances and the norm explicitly when comparing applications with different scales.

Check the inherited status after each solve. Converged indicates accepted convergence; IterationLimitExceeded indicates exhausted iterations, and RoundOffError can indicate numerical breakdown or an unattainable tolerance. The inherited throw-on-failure option controls whether algorithm failure throws. Invalid arguments and factorization failures can still throw.

Residuals contains the true original-system residual at termination. The error reported by SolutionReport is a residual measure, not a bound on solution error. Ill-conditioned systems can have a small residual and an inaccurate solution.

The CG example reads an SPD Matrix Market file using MatrixMarketFile, constructs threshold IC with explicit ordering, scaling and a small factor-only shift, and checks status. The caller supplies a matching RHS. The tuning values illustrate the API and require adjustment for other matrices.

SparseMatrix<double> matrix = (SparseMatrix<double>)MatrixMarketFile.ReadMatrix<double>(matrixPath);
Permutation ordering = matrix.GetReverseCuthillMcKeeOrdering();
ThresholdIncompleteCholeskyPreconditioner<double> factor =
    new ThresholdIncompleteCholeskyPreconditioner<double>(
        matrix, 1e-4, 16, 1e-14, .01, equilibrate: true, ordering: ordering);
ConjugateGradientSolver<double> solver = new ConjugateGradientSolver<double>(matrix);
solver.Preconditioner = factor;
solver.MaxIterations = 2000;
solver.ConvergenceTests.AbsoluteTolerance = 0;
solver.ConvergenceTests.RelativeTolerance = 1e-12;
Vector<double> solution = solver.Solve(rhs);
if (solver.Status != AlgorithmStatus.Converged)
    throw new InvalidOperationException("CG did not converge: " + solver.Status);
Console.WriteLine("Original residual norm: " + solver.Residuals.Norm());

Dim matrix = DirectCast(MatrixMarketFile.ReadMatrix(Of Double)(matrixPath), SparseMatrix(Of Double))
Dim ordering = matrix.GetReverseCuthillMcKeeOrdering()
Dim factor = New ThresholdIncompleteCholeskyPreconditioner(Of Double)(
    matrix, 1.0E-4, 16, 1.0E-14, 0.01, equilibrate:=True, ordering:=ordering)
Dim solver = New ConjugateGradientSolver(Of Double)(matrix)
solver.Preconditioner = factor
solver.MaxIterations = 2000
solver.ConvergenceTests.AbsoluteTolerance = 0
solver.ConvergenceTests.RelativeTolerance = 1.0E-12
Dim solution = solver.Solve(rhs)
If solver.Status <> AlgorithmStatus.Converged Then
    Throw New InvalidOperationException("CG did not converge: " & solver.Status.ToString())
End If
Console.WriteLine("Original residual norm: " & solver.Residuals.Norm().ToString())

No code example is currently available or this language may not be supported.

let matrix = MatrixMarketFile.ReadMatrix<float>(matrixPath) :?> SparseMatrix<float>
let ordering = matrix.GetReverseCuthillMcKeeOrdering()
let factor = new ThresholdIncompleteCholeskyPreconditioner<float>(
    matrix, 1e-4, 16, 1e-14, 0.01, equilibrate = true, ordering = ordering)
let solver = new ConjugateGradientSolver<float>(matrix)
solver.Preconditioner <- factor
solver.MaxIterations <- 2000
solver.ConvergenceTests.AbsoluteTolerance <- 0.0
solver.ConvergenceTests.RelativeTolerance <- 1e-12
let solution = solver.Solve(rhs)
if solver.Status <> AlgorithmStatus.Converged then
    invalidOp ("CG did not converge: " + string solver.Status)
printfn "Original residual norm: %g" (solver.Residuals.Norm())

This MINRES example solves a small symmetric indefinite system using its default identity preconditioner.

SparseCompressedColumnMatrix<double> matrix = Matrix.CreateSparse<double>(2, 2);
matrix[0, 0] = 2; matrix[0, 1] = 1;
matrix[1, 0] = 1; matrix[1, 1] = -1;
Vector<double> rhs = Vector.Create(4.0, -1.0);
MinimumResidualSolver<double> solver = new MinimumResidualSolver<double>(matrix);
// Identity is the default; an alternative must be fixed, linear and HPD.
solver.ConvergenceTests.AbsoluteTolerance = 0;
solver.ConvergenceTests.RelativeTolerance = 1e-12;
Vector<double> solution = solver.Solve(rhs);
if (solver.Status != AlgorithmStatus.Converged)
    throw new InvalidOperationException("MINRES did not converge: " + solver.Status);

Dim matrix = Global.Numerics.NET.Matrix.CreateSparse(Of Double)(2, 2)
matrix(0, 0) = 2 : matrix(0, 1) = 1
matrix(1, 0) = 1 : matrix(1, 1) = -1
Dim rhs = Vector.Create(4.0, -1.0)
Dim solver = New MinimumResidualSolver(Of Double)(matrix)
' Identity is the default; an alternative must be fixed, linear and HPD.
solver.ConvergenceTests.AbsoluteTolerance = 0
solver.ConvergenceTests.RelativeTolerance = 1.0E-12
Dim solution = solver.Solve(rhs)
If solver.Status <> AlgorithmStatus.Converged Then
    Throw New InvalidOperationException("MINRES did not converge: " & solver.Status.ToString())
End If

No code example is currently available or this language may not be supported.

let matrix = Matrix.CreateSparse<float>(2, 2)
matrix.[0, 0] <- 2.0
matrix.[0, 1] <- 1.0
matrix.[1, 0] <- 1.0
matrix.[1, 1] <- -1.0
let rhs = Vector.Create(4.0, -1.0)
let solver = new MinimumResidualSolver<float>(matrix)
// Identity is the default; an alternative must be fixed, linear and HPD.
solver.ConvergenceTests.AbsoluteTolerance <- 0.0
solver.ConvergenceTests.RelativeTolerance <- 1e-12
let solution = solver.Solve(rhs)
if solver.Status <> AlgorithmStatus.Converged then
    invalidOp ("MINRES did not converge: " + string solver.Status)

The FGMRES example reads a general sparse Matrix Market file and maps an ILUT factor from equilibrated coordinates. This example uses a fixed factor, which FGMRES also supports. Choose FGMRES when a custom preconditioner varies or uses an inner nonlinear solve; ordinary GMRES is an alternative for a fixed linear action.

SparseMatrix<double> matrix = (SparseMatrix<double>)MatrixMarketFile.ReadMatrix<double>(matrixPath);
Permutation ordering = matrix.GetReverseCuthillMcKeeOrdering();
ThresholdIncompleteLUPreconditioner<double> factor =
    new ThresholdIncompleteLUPreconditioner<double>(matrix, 1e-4, 24, 1e-14, 0,
        equilibrateRows: true, ordering: ordering);
FlexibleGeneralizedMinimalResidualSolver<double> solver =
    new FlexibleGeneralizedMinimalResidualSolver<double>(matrix);
solver.Preconditioner = factor;
solver.IterationsBetweenRestarts = 60;
solver.MaxIterations = 2000;
solver.ConvergenceTests.AbsoluteTolerance = 0;
solver.ConvergenceTests.RelativeTolerance = 1e-12;
Vector<double> solution = solver.Solve(rhs);
if (solver.Status != AlgorithmStatus.Converged)
    throw new InvalidOperationException("FGMRES did not converge: " + solver.Status);
// The factor handles ordering and row scaling internally.
Console.WriteLine("Original residual norm: " + solver.Residuals.Norm());

Dim matrix = DirectCast(MatrixMarketFile.ReadMatrix(Of Double)(matrixPath), SparseMatrix(Of Double))
Dim ordering = matrix.GetReverseCuthillMcKeeOrdering()
Dim factor = New ThresholdIncompleteLUPreconditioner(Of Double)(
    matrix, 1.0E-4, 24, 1.0E-14, 0, equilibrateRows:=True, ordering:=ordering)
Dim solver = New FlexibleGeneralizedMinimalResidualSolver(Of Double)(matrix)
solver.Preconditioner = factor
solver.IterationsBetweenRestarts = 60
solver.MaxIterations = 2000
solver.ConvergenceTests.AbsoluteTolerance = 0
solver.ConvergenceTests.RelativeTolerance = 1.0E-12
Dim solution = solver.Solve(rhs)
If solver.Status <> AlgorithmStatus.Converged Then
    Throw New InvalidOperationException("FGMRES did not converge: " & solver.Status.ToString())
End If
' The factor handles ordering and row scaling internally.
Console.WriteLine("Original residual norm: " & solver.Residuals.Norm().ToString())

No code example is currently available or this language may not be supported.

let matrix = MatrixMarketFile.ReadMatrix<float>(matrixPath) :?> SparseMatrix<float>
let ordering = matrix.GetReverseCuthillMcKeeOrdering()
let factor = new ThresholdIncompleteLUPreconditioner<float>(
    matrix, 1e-4, 24, 1e-14, 0.0, equilibrateRows = true, ordering = ordering)
let solver = new FlexibleGeneralizedMinimalResidualSolver<float>(matrix)
solver.Preconditioner <- factor
solver.IterationsBetweenRestarts <- 60
solver.MaxIterations <- 2000
solver.ConvergenceTests.AbsoluteTolerance <- 0.0
solver.ConvergenceTests.RelativeTolerance <- 1e-12
let solution = solver.Solve(rhs)
if solver.Status <> AlgorithmStatus.Converged then
    invalidOp ("FGMRES did not converge: " + string solver.Status)
// The factor handles ordering and row scaling internally.
printfn "Original residual norm: %g" (solver.Residuals.Norm())

Reuse a solver for sequential solves to amortize workspace allocation. Solver instances contain mutable iteration state and must not run concurrent solves. The returned solution and residual vectors belong to the solver; copy them before a later solve if they must be retained.

The sparse ILU, Cholesky and SSOR factors snapshot their coefficients and retain their source for explicit refresh. Changing the source matrix does not automatically update an existing factor. Their refresh methods rebuild factors at unchanged dimensions, accept changed sparsity where documented, and preserve the previous usable factor if refresh fails.

These sparse factors support concurrent applications with distinct destinations. Refresh requires exclusive access. Transformations snapshot ordering, scaling and matrix coefficients; after changing the original matrix, rebuild the transformation and mapped factor together. Refreshing a factor on an old transformed matrix cannot incorporate changes to the original matrix.