The sparse-linear-algebra package

[Tags:gpl, library, test]

Overview

The sparse-linear-algebra library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. The user interface is provided by the top-level module Numeric.LinearAlgebra.Sparse:

import           Numeric.LinearAlgebra.Sparse

Please refer to the README file for usage examples.


[Skip to Readme]

Properties

Versions 0.1.0.0, 0.1.0.1, 0.1.0.2, 0.1.0.3, 0.2.0.0, 0.2.0.1, 0.2.0.2, 0.2.0.3, 0.2.0.4, 0.2.0.5, 0.2.0.7, 0.2.0.8, 0.2.0.9, 0.2.1.0, 0.2.1.1, 0.2.2.0, 0.2.9, 0.2.9.1, 0.2.9.2, 0.2.9.3, 0.2.9.4, 0.2.9.5, 0.2.9.6, 0.2.9.7 (info)
Change log CHANGELOG.markdown
Dependencies base (>=4.7 && <5), containers, exceptions, mtl (>=2.2.1), transformers, vector, vector-algorithms, vector-space [details]
License GPL-3
Copyright 2016-2017 Marco Zocca
Author Marco Zocca
Maintainer zocca.marco gmail
Category Numeric
Home page https://github.com/ocramz/sparse-linear-algebra
Source repository head: git clone https://github.com/ocramz/sparse-linear-algebra
Uploaded Sat Feb 25 22:31:52 UTC 2017 by ocramz
Distributions LTSHaskell:0.2.9.7, NixOS:0.2.9.7, Stackage:0.2.9.7, Tumbleweed:0.2.9.7
Downloads 534 total (21 in the last 30 days)
Votes
1 []
Status Docs available [build log]
Last success reported on 2017-02-25 [all 1 reports]

Modules

[Index]

Downloads

Maintainer's Corner

For package maintainers and hackage trustees

Readme for sparse-linear-algebra

Readme for sparse-linear-algebra-0.2.9.7

sparse-linear-algebra

Numerical computation in native Haskell

Hackage Build Status

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) (non-Hermitian systems)

    • BiConjugate Gradient (BCG)

    • Conjugate Gradient Squared (CGS)

    • BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems)

    • 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)

  • Eigenvalue algorithms

    • QR (eigsQR)

    • QR-Arnoldi (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)
  • Matrix factorization algorithms

    • Golub-Kahan-Lanczos bidiagonalization (gklBidiag)

    • Singular value decomposition (SVD)

  • Iterative linear solvers

    • Transpose-Free Quasi-Minimal 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 pretty-printed 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 non-zero 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 matrix-matrix 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 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 Floats and 10^-12 for Doubles).

λ> 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 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 ( 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 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
( 3 elements ) ,  3 NZ ( sparsity 1.0 )

1.4999999999999998 -1.9999999999999998 0.9999999999999998

The result can be verified by computing the matrix-vector action amat #> x, which should (ideally) be very close to the right-hand side b :

λ> prd $ amat #> x
( 3 elements ) ,  3 NZ ( sparsity 1.0 )

2.9999999999999996 1.9999999999999996 4.999999999999999

The library also provides a forward-backward 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 cross-verify 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

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