sparse-linear-algebra

KO | EN

gitlite — search

#react #rust #linux #docker #macos #nodejs #windows #cli #vue #awesome-list #pytorch #ios #openai #open-source

sparse-linear-algebra

(★ 95)

Numerical computation in native Haskell

sparse-linear-algebra Latest Version Download

Download Latest Version (.zip)

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

# sparse-linear-algebra

Numerical computation in native Haskell

[![CI-library](https://github.com/ocramz/sparse-linear-algebra/actions/workflows/ci.yml/badge.svg?branch=master)](https://github.com/ocramz/sparse-linear-algebra/actions/workflows/ci.yml)
[![CI-matrix-factorizations](https://github.com/ocramz/sparse-linear-algebra/actions/workflows/matrix_factorizations.yml/badge.svg)](https://github.com/ocramz/sparse-linear-algebra/actions/workflows/matrix_factorizations.yml)
[![Hackage](https://img.shields.io/hackage/v/sparse-linear-algebra.svg)](https://hackage.haskell.org/package/sparse-linear-algebra)  
[![sparse-linear-algebra](http://stackage.org/package/sparse-linear-algebra/badge/lts)](http://stackage.org/lts/package/sparse-linear-algebra)
[![sparse-linear-algebra](http://stackage.org/package/sparse-linear-algebra/badge/nightly)](http://stackage.org/nightly/package/sparse-linear-algebra)

This library provides common numerical analysis functionality, without requiring any external bindings. It aims to serve as an experimental platform for scientific computation in a purely functional setting.

## Project status

January 2026: Re-enabled BiCGSTAB (Biconjugate Gradient Stabilized) solver with comprehensive tests. BiCGSTAB is more numerically stable than CGS and is the recommended choice for non-Hermitian systems. Also fixed CGS (Conjugate Gradient Squared) numerical convergence issues. Both solvers now properly handle early termination and use appropriate tolerances for iterative methods.

**Note on property tests**: Property-based tests for iterative solvers (BiCGSTAB, CGS) guard against degenerate cases (tiny systems, nearly zero RHS/solution, very sparse matrices) which can cause flaky behavior with randomly generated SPD matrices. This is expected behavior and reflects the numerical limitations of iterative methods on ill-conditioned systems.

December 2025: The project sat unmaintained for a few years but I'm not comfortable with leaving projects incomplete. There are some hard design problems under the hood (inefficient data representation with lots of memory copies, incorrect algorithms, numerical instability) that make progress a slog.

Mar 14, 2018: The core linear algebra operations work, but there are still a few (documented) bugs such as in the matrix factorizations department. Complex number support is still incomplete, so the users are advised to not rely on that for the time being. The issues related to Complex number handling are tracked in #51, #12, #30.

Refer to the Changelog for detailed status updates.

## Contents

* "BLAS" levels 1 to 3 : vector, matrix operations including matrix products etc.

* Iterative linear solvers (`<\>`)

* Generalized Minimal Residual (GMRES) (non-Hermitian systems) 🚧

* BiConjugate Gradient (BCG) 🚧

* Conjugate Gradient Squared (CGS) ✅
      - Suitable for non-Hermitian systems
      - Fast convergence for well-conditioned problems
      - May exhibit irregular convergence on ill-conditioned systems
      - Use `cgsInit` and `cgsStep` for manual iteration control
      - Supports debug tracing with `cgsStepDebug`

* BiConjugate Gradient Stabilized (BiCGSTAB) ✅
      - Suitable for non-Hermitian systems
      - More numerically stable than CGS
      - Recommended for ill-conditioned problems
      - Smoother convergence behavior than CGS
      - Use `bicgsInit` and `bicgstabStep` for manual iteration control
      - Access via `linSolve0` interface with `BICGSTAB_` method

* Transpose-Free Quasi-Minimal Residual (TFQMR) 🚧

* Moore-Penrose pseudoinverse (`pinv`) (rectangular systems)

* Direct linear solvers
    * LU-based (`luSolve`); forward and backward substitution (`triLowerSolve`, `triUpperSolve`)
    
 
* Matrix factorization algorithms
    * QR (`qr`) 🚧
    * LU (`lu`) 🚧
    * Cholesky (`chol`) 🚧
    * Arnoldi iteration (`arnoldi`) 🚧
    * Golub-Kahan-Lanczos bidiagonalization (`gklBidiag`) 🚧
    * Singular value decomposition (SVD) 🚧

* Eigenvalue algorithms 
    * QR (`eigsQR`) 🚧
    * QR-Arnoldi (`eigsArnoldi`) 🚧
    * Rayleigh quotient iteration (`eigRayleigh`) 🚧

* Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection

* Predicates : Matrix orthogonality test (A^T A ~= I)

---------

## Examples

The module `Numeric.LinearAlgebra.Sparse` contains the user interface.

### Creation of sparse data

The `fromListSM` function creates a sparse matrix from a collection of its entries in (row, column, value) format. This is its type signature:

fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a

and, in case you have a running GHCi session (the terminal is denoted from now on by `λ>`), you can try something like this:

λ> amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)] :: SpMatrix Double

Similarly, `fromListSV` is used to create sparse vectors:

fromListSV :: Int -> [(Int, a)] -> SpVector a

Alternatively, the user can copy the contents of a list to a (dense) SpVector using

fromListDenseSV :: Int -> [a] -> SpVector a

### Displaying sparse data

Both sparse vectors and matrices can be pretty-printed using `prd`:

λ> prd amat

( 3 rows, 3 columns ) , 5 NZ ( density 55.556 % )

2.00   , _      , _      
    4.00   , 3.00   , 2.00   
    _      , _      , 5.00

*Note (sparse storage)*: sparse data should only contain non-zero entries not to waste memory and computation.

*Note (approximate output)*: `prd` rounds the results to two significant digits, and switches to scientific notation for large or small values. Moreover, values which are indistinguishable from 0 (see the `Numeric.Eps` module) are printed as `_`.

### Matrix factorizations, matrix product

There are a few common matrix factorizations available; in the following example we compute the LU factorization of matrix `amat` and verify it with the matrix-matrix product `##` of its factors :

λ> (l, u) <- lu amat
    λ> prd $ l ## u
    
    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % )

2.00   , _      , _      
    4.00   , 3.00   , 2.00   
    _      , _      , 5.00

Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!).
To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) depends on the numerical type used (e.g. it is 10^-6 for `Float`s and 10^-12 for `Double`s).

λ> prd $ l #~# u
    
    ( 3 rows, 3 columns ) , 5 NZ ( density 55.556 % )

2.00   , _      , _      
    4.00   , 3.00   , 2.00   
    _      , _      , 5.00

A matrix is transposed using the `transpose` function.

Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (i.e. matrix transpose * matrix) and `##^` (matrix * matrix transpose):

λ> amat' = amat #^# amat
    λ> prd amat'
    
    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % )

20.00  , 12.00  , 8.00   
    12.00  , 9.00   , 6.00   
    8.00   , 6.00   , 29.00

λ> lc <- chol amat'
    λ> prd $ lc ##^ lc
    
    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % )

20.00  , 12.00  , 8.00   
    12.00  , 9.00   , 10.80  
    8.00   , 10.80  , 29.00

In the last example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only defined for symmetric positive-definite matrices.

### Linear systems

Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `<\>` (inspired by Matlab's "backslash" function. Currently this method uses GMRES as default) :

λ> b = fromListDenseSV 3 [3,2,5] :: SpVector Double
    λ> x <- amat <\> b
    λ> prd x