sparselinearalgebra
Numerical computation in native Haskell
https://github.com/ocramz/sparselinearalgebra
Version on this page:  0.2.9.7 
LTS Haskell 17.11:  0.3.1 
Stackage Nightly 20210510:  0.3.1 
Latest on Hackage:  0.3.1 
sparselinearalgebra0.2.9.7@sha256:a108761785de5caa59d4eaab17d1ec3fc785d7ef78a5e3533f5fb479c6384d0f,4094
Module documentation for 0.2.9.7
 Control
 Control.Exception
 Data
 Numeric
 Numeric.Eps
 Numeric.LinearAlgebra
sparselinearalgebra
Numerical computation in native Haskell
This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting.
Contents :

Iterative linear solvers (
<\>
)
Generalized Minimal Residual (GMRES) (nonHermitian systems)

BiConjugate Gradient (BCG)

Conjugate Gradient Squared (CGS)

BiConjugate Gradient Stabilized (BiCGSTAB) (nonHermitian systems)

MoorePenrose pseudoinverse (
pinv
) (rectangular systems)


Direct linear solvers
 LUbased (
luSolve
); forward and backward substitution (triLowerSolve
,triUpperSolve
)
 LUbased (

Matrix factorization algorithms

QR (
qr
) 
LU (
lu
) 
Cholesky (
chol
) 
Arnoldi iteration (
arnoldi
)


Eigenvalue algorithms

QR (
eigsQR
) 
QRArnoldi (
eigsArnoldi
)


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

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

Eigenvalue algorithms
 Rayleigh quotient iteration (
eigRayleigh
)
 Rayleigh quotient iteration (

Matrix factorization algorithms

GolubKahanLanczos bidiagonalization (
gklBidiag
) 
Singular value decomposition (SVD)


Iterative linear solvers
 TransposeFree QuasiMinimal Residual (TFQMR)
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 prettyprinted using prd
:
λ> prd amat
( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )
2.0 _ _
4.0 3.0 2.0
_ _ 5.0
Note: sparse data should only contain nonzero entries not to waste memory and computation.
Matrix operations
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 matrixmatrix product ##
of its factors :
λ> (l, u) < lu amat
λ> prd $ l ## u
( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )
2.0 _ _
4.0 3.0 2.0
_ _ 5.0
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 matrixmatrix 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 ( sparsity 0.5555555555555556 )
2.0 _ _
4.0 3.0 2.0
_ _ 5.0
A matrix is transposed using the transpose
function.
Sometimes we need to compute matrixmatrix 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 ( sparsity 1.0 )
20.0 12.0 8.0
12.0 9.0 6.0
8.0 6.0 29.0
λ> l < chol amat'
λ> prd $ l ##^ l
( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )
20.000000000000004 12.0 8.0
12.0 9.0 10.8
8.0 10.8 29.0
In the last example we have also shown the Cholesky decomposition (M = L L^T where L is a lowertriangular matrix), which is only defined for symmetric positivedefinite matrices.
Linear systems
Large sparse linear systems are best solved with iterative methods. sparselinearalgebra
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
( 3 elements ) , 3 NZ ( sparsity 1.0 )
1.4999999999999998 1.9999999999999998 0.9999999999999998
The result can be verified by computing the matrixvector action amat #> x
, which should (ideally) be very close to the righthand side b
:
λ> prd $ amat #> x
( 3 elements ) , 3 NZ ( sparsity 1.0 )
2.9999999999999996 1.9999999999999996 4.999999999999999
The library also provides a forwardbackward substitution solver (luSolve
) based on a triangular factorization of the system matrix (usually LU). This should be the preferred for solving smaller, dense systems. Using the data defined above we can crossverify the two solution methods:
λ> x' < luSolve l u b
λ> prd x'
( 3 elements ) , 3 NZ ( sparsity 1.0 )
1.5 2.0 1.0
License
GPL3, see LICENSE
Credits
Inspired by
linear
: https://hackage.haskell.org/package/linearvectorspace
: https://hackage.haskell.org/package/vectorspacesparselinalg
: https://github.com/laughedelic/sparselinalg
References
[1] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000
[2] G.H. Golub and C.F. Van Loan, Matrix Computations, 3rd ed., 1996
[3] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997
[4] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992
Changes
0.2.9.4

Exceptions constructors are exported by Numeric.LinearAlgebra.Sparse
0.2.9.1

* Uses classes from `vectorspace` such as AdditiveGroup, VectorSpace and InnerSpace
* QuickCheck tests for algebraic properties, such as matrixvector products and soon more abstract ones e.g. positive semidefinite matrices
* Getting rid of `error` in favor of MonadThrow exceptions for highlevel operations
such as matrix algorithms