Extreme Optimization™: Complexity made simple.

Math and Statistics
Libraries for .NET

  • Home
  • Features
    • Math Library
    • Vector and Matrix Library
    • Statistics Library
    • Performance
    • Usability
  • Documentation
    • Introduction
    • Math Library User's Guide
    • Vector and Matrix Library User's Guide
    • Data Analysis Library User's Guide
    • Statistics Library User's Guide
    • Reference
  • Resources
    • Downloads
    • QuickStart Samples
    • Sample Applications
    • Frequently Asked Questions
    • Technical Support
  • Blog
  • Order
  • Company
    • About us
    • Testimonials
    • Customers
    • Press Releases
    • Careers
    • Partners
    • Contact us
Introduction
Deployment Guide
Nuget packages
Configuration
Using Parallelism
Expand Mathematics Library User's GuideMathematics Library User's Guide
Expand Vector and Matrix Library User's GuideVector and Matrix Library User's Guide
Expand Data Analysis Library User's GuideData Analysis Library User's Guide
Expand Statistics Library User's GuideStatistics Library User's Guide
Expand Data Access Library User's GuideData Access Library User's Guide
Expand ReferenceReference
  • Extreme Optimization
    • Features
    • Solutions
    • Documentation
    • QuickStart Samples
    • Sample Applications
    • Downloads
    • Technical Support
    • Download trial
    • How to buy
    • Blog
    • Company
    • Resources
  • Documentation
    • Introduction
    • Deployment Guide
    • Nuget packages
    • Configuration
    • Using Parallelism
    • Mathematics Library User's Guide
    • Vector and Matrix Library User's Guide
    • Data Analysis Library User's Guide
    • Statistics Library User's Guide
    • Data Access Library User's Guide
    • Reference
  • Vector and Matrix Library User's Guide
    • Basic Concepts
    • Vectors
    • Matrices
    • Structured Matrix Types
    • Matrix Decompositions
    • Sparse Vectors and Matrices
    • Complex Linear Algebra
    • Single-Precision Linear Algebra
    • Distributed and GPU Computing
  • Sparse Vectors and Matrices
    • Sparse Vectors
    • Sparse Matrices
    • Solving Sparse Systems
  • Solving Sparse Systems

Solving Sparse Systems

Extreme Optimization Numerical Libraries for .NET Professional

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 Extreme.Mathematics.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.

Stationary iterative solver classes

Class

Description

JacobiSolverT

Implements simple Jacobi iteration.

GaussSeidelSolverT

Implements Gauss-Seidel iteration.

SuccessiveOverRelaxationSolverT

Implements Successive Over-Relaxation. Requires a relaxation factor between 0 and 2.

Non-Stationary iterative solver classes

Class

Description

ConjugateGradientSolverT

Implements the Conjugate Gradient (CG) method, suitable for symmetric systems. Convergence can be erratic, and may sometimes break down completely.

ConjugateGradientSquaredSolverT

Implements the Conjugate Gradient Squared (CGS) method, suitable for symmetric and non-symmetric systems. Convergence can be erratic, and may sometimes break down completely.

BiConjugateGradientSolverT

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.

BiConjugateGradientStabilizedSolverT

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.

QuasiMinimalResidualSolverT

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.

QuasiMinimalResidualSolverT

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 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:

Class

Description

IdentityPreconditionerT

Implements the identity preconditioner. It has no effect and is equivalent to not using a preconditioner at all.

JacobiPreconditionerT

Implements the Jacobi preconditioner, which approximates the matrix by its diagonal. This is the default preconditioner for matrices whose elements can be individually accessed.

IncompleteLUPreconditionerT

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 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.

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:

C#
VB
C++
F#
Copy
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)

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

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

Copyright (c) 2004-2021 ExoAnalytics Inc.

Send comments on this topic to support@extremeoptimization.com

Copyright © 2004-2021, 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 Optimized for Visual Studio logo
are registered trademarks of Microsoft Corporation.