Sparse linear algebra

ALGLIB numerical analysis library provides a rich set of sparse matrix functions available from C++, C#, Java and several other programming languages. The library supports several sparse matrix storage formats, sparse BLAS (sparse GEMV and its variants), factorizations (sparse Cholesky, LDLT and LU), direct and iterative sparse linear solvers.

Being a dual-licensed library, ALGLIB includes a free and a commercial edition. A fully functional free edition is ideal for academic research or small non-commercial projects. The commercial edition is intended for those seeking high performance, priority support, or the ability to use the library in commercial applications.

Getting started

Sparse matrix storage formats and initialization

ALGLIB supports several sparse matrix storage formats optimized for various sparsity patterns and operations:

• CRS (Compressed Row Storage) - a format optimized for sparse linear algebra operations. This format efficiently stores matrices with any sparsity pattern. However, it lacks flexible initialization (matrices are always initialized from top to bottom and from left to right) and an ability to modify sparsity pattern on the fly.
• HTS (Hash Table Storage) - a format optimized for easy initialization. Matrix elements are stored in a hash table, and one can easily add or remove elements. However, since most operations (including linear algebra) are not supported for this format, so one has to convert HTS to CRS (preferably) in order to use a sparse matrix.
• SKS (Skyline Storage) - a format optimized for low-profile matrices, including ones with low bandwidth.

We recommend using HTS format first because of its easy and convenient initialization, followed by conversion to CRS as soon as possible, as it offers the best performance.

Using sparse matrices

The basis sparse matrix functionality is provided in a sparse subpackage of the LinAlg package. This functionality includes:

• A rich set of matrix initialization functions with the ability to create matrices in CRS, HTS, or SKS formats, either from scratch or by reusing previously allocated memory
• A set of conversion functions allowing one to perform either in-place or out-of-place conversions between various sparse storage formats
• Additional functions, such as transposition, permutation, serialization, enumeration, and so on

The ALGLIB Reference Manual includes several examples of sparse matrix initialization and usage, including sparse_d_1 and sparse_d_crs.

Sparse linear algebra functions

Sparse BLAS

The sparse subpackage of the LinAlg package also provides a rich set of fast and efficient sparse BLAS functions, including:

• Sparse-dense matrix-matrix products of various types, including (S·A and S'·A)
• Matrix-vector products, including S·x, S'·x, and x'·S·x
• Triangular products and solvers, i.e., products in the form of S·x and S^-1·x, with S being a sparse triangular matrix

Our sparse linear algebra implementation efficiently handles both small and extremely large sparse matrices. As long as one has enough RAM, it is possible to handle matrices with more than 1.000.000.000 elements.

Sparse Cholesky and other triangular factorizations

Sparse triangular factorization is an important part of any sparse matrix package. The following sparse factorizations are provided in the trfac subpackage:

• Sparse Cholesky decomposition (supernodal) with AMD (Approximate Minimum Degree) ordering, capable of handling matrices with more than 1.000.000 rows. ALGLIB provides a flexible and feature-rich API that allows one to perform either one-off factorization or efficient multiple factorizations with the same sparsity patterns ('analyze'/'factorize' approach).
• Sparse LU decomposition with complete, row, or no pivoting

The commercial version of the sparse Cholesky solver supports multithreading. However, specific speed-up heavily depends on a sparsity pattern of a matrix being factorized.

Sparse direct linear solvers

A direct linear solver is one that relies on some kind of (sparse) triangular factorization to solve a linear system. ALGLIB provides several direct sparse solvers in its directsparsesolvers subpackage:

• A Cholesky-based sparse solver for symmetric positive definite systems
• An LU-based sparse solver for general (nonsymmetric) linear systems

Direct solvers usually provide the most accurate results - at least 6 digits of precision (often, many more than that). The biggest drawback is that they have unpredictable memory/time requirements - it is hard to tell in advance how much memory and time will be needed to perform a factorization.

Sparse iterative linear solvers

An iterative linear solver usually solves a linear system through a sequence of matrix-vector products. Sparse iterative solvers are provided by the iterativesparse, linlsqr, and lincg subpackages, which include the following algorithms:

• GMRES - a well-known algorithm for sparse general (nonsymmetric) systems, with in-core and out-of-core versions
• LSQR - an algorithm for sparse linear least squares. Originally intended for least squares problems with rectangular MxN matrices, it can also be used as a solver for problems with square system matrices.
• CG - a well-known algorithm for symmetric positive definite systems

When compared with direct solvers, iterative solvers usually have much more predictable memory requirements but much less predictable accuracy. Usually, they allocate only a limited amount of additional memory bounded by O(max(M,N), which can easily be predicted in advance. What is difficult to predict is how many iterations will be necessary to converge to 6 digits of precision.

Sparse eigensolvers

Eigenproblems with sparse matrices are an important subset of linear algebra. evd subpackage includes a subspace iteration eigensolver which can be used to compute several top eigenvalues/vectors of a large sparse matrix.

Other sparsity-aware functions

Many other algorithms in ALGLIB can utilize the sparsity of their inputs, such as:

• Linear programming and quadratic programming solvers, which are heavily optimized for sparse problems
• Large-scale nonlinear programming solvers, which work much faster when constraints are sparse