sparse-linear-algebra: Numerical computation in native Haskell

[ gpl, library, numeric ] [ Propose Tags ]
This version is deprecated.

Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details


[Skip to Readme]

Downloads

Maintainer's Corner

Package maintainers

For package maintainers and hackage trustees

Candidates

  • No Candidates
Versions [RSS] 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, 0.2.9.8, 0.2.9.9, 0.3, 0.3.1 (info)
Dependencies base (>=4.7 && <5), containers, hspec, mtl (>=2.2.1), mwc-random, primitive (>=0.6.1.0), QuickCheck [details]
License GPL-3.0-only
Copyright 2016 Marco Zocca
Author Marco Zocca
Maintainer zocca.marco gmail
Category Numeric
Home page https://github.com/ocramz/sparse-linear-algebra
Source repo head: git clone https://github.com/ocramz/sparse-linear-algebra
Uploaded by ocramz at 2016-10-23T23:24:43Z
Distributions
Reverse Dependencies 8 direct, 0 indirect [details]
Executables sparse-linear-algebra
Downloads 12120 total (73 in the last 30 days)
Rating 2.0 (votes: 1) [estimated by Bayesian average]
Your Rating
  • λ
  • λ
  • λ
Status Docs available [build log]
Last success reported on 2016-11-19 [all 1 reports]

Readme for sparse-linear-algebra-0.2.0.8

[back to package description]

sparse-linear-algebra

Numerical computation in native Haskell

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

Algorithms :

  • Iterative linear solvers

    • BiConjugate Gradient (BCG)

    • Conjugate Gradient Squared (CGS)

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

    • Transpose-Free Quasi-Minimal Residual (TFQMR)

  • Matrix decompositions

    • QR factorization

    • LU factorization

  • Eigenvalue algorithms

    • QR algorithm

    • Rayleigh quotient iteration

  • 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 and pretty-printing

To create a sparse matrix from an array of its entries we use fromListSM :

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

e.g.

> amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)]

And similarly for sparse vectors : fromListSV :: Int -> [(Int, a)] -> SpVector a.

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

> prd amat
( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )

[2,0,0]
[4,3,2]
[0,0,5]

The zeros are just added at pretty printing time; sparse vectors and matrices should only contain non-zero entries.

Matrix operations

Matrix factorizations are available as lu and qr respectively, and are straightforward to verify by using the matrix product ## :

> (l, u) = lu amat
> prd $ l ## u
( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )

[2.0,0.0,0.0]
[4.0,3.0,2.0]
[0.0,0.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, is fixed at 10^-8.

> prd $ l #~# u
( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )

[2.0,0.0,0.0]
[4.0,3.0,2.0]
[0.0,0.0,5.0]

Linear systems

Linear systems can be solved with either linSolve (which also requires choosing a method) or with <\> (which uses BiCGSTAB as default) :

> b = fromListSV 3 [(0,3),(1,2),(2,5)]
> 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]

This is also an experiment in principled scientific programming :

  • set the stage by declaring typeclasses and some useful generic operations (normed linear vector spaces, i.e. finite-dimensional spaces equipped with an inner product that induces a distance function),

  • define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via Data.IntMap, are made instances of the VectorSpace and Additive classes respectively). This allows to decouple the algorithms from the actual implementation of the backend,

  • implement the algorithms, following 1:1 the textbook [1]

License

GPL3, see LICENSE

Credits

Inspired by

References

[1] : Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000