# sparse-linear-algebra

Numerical computation in native Haskell https://github.com/ocramz/sparse-linear-algebra

LTS Haskell 9.18: | 0.2.9.8 |

Stackage Nightly 2017-12-15: | 0.2.9.8 |

Latest on Hackage: | 0.2.9.8 |

**Marco Zocca**

**zocca.marco gmail**

#### Module documentation for 0.2.9.8

- Control
- Control.Exception

- Data
- Numeric
- Numeric.Eps
- Numeric.LinearAlgebra

# sparse-linear-algebra

Numerical computation in native Haskell

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.

## News

Oct 7., 2017: The library is evolving in a number of ways, to reflect performance observations and user requests:

typeclasses and instances for primitive types will become

`sparse-linear-algebra-core`

, along with a typeclass-oriented reformulation of the numerical algorithms that used to depend on the nested IntMap representation. This will let other developers build on top of this library, in the spirit of`vector-space`

and`linear`

.The

`vector`

-based backend is being reworked.An

`accelerate`

-based backend is under development.

## 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`

)

- LU-based (
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`

)

- Rayleigh quotient iteration (
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 ( 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
( 3 elements ) , 3 NZ ( density 100.000 % )
1.50 , -2.00 , 1.00
```

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 ( density 100.000 % )
3.00 , 2.00 , 5.00
```

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 LU factors defined previously we can cross-verify the two solution methods:

```
λ> x' <- luSolve l u b
λ> prd x'
( 3 elements ) , 3 NZ ( density 100.000 % )
1.50 , -2.00 , 1.00
```

## License

GPL3, see LICENSE

## Credits

Inspired by

`linear`

: https://hackage.haskell.org/package/linear`vector-space`

: https://hackage.haskell.org/package/vector-space`sparse-lin-alg`

: https://github.com/laughedelic/sparse-lin-alg

## 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] T.A. Davis, Direct Methods for Sparse Linear Systems, 2006

[4] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997

[5] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992

## Changes

```
0.2.9.7
------
Improved pretty printer:
* Fixed display precision (e.g. 2 decimal digits), fixed column width output for vectors and matrices
* Small and large values (wrt fixed precision) switch to scientific notation
0.2.9.4
-------
Exceptions constructors are exported by Numeric.LinearAlgebra.Sparse
0.2.9.1
-------
* Uses classes from `vector-space` such as AdditiveGroup, VectorSpace and InnerSpace
* QuickCheck tests for algebraic properties, such as matrix-vector products and soon more abstract ones e.g. positive semi-definite matrices
* Getting rid of `error` in favor of MonadThrow exceptions for high-level operations such as matrix algorithms
```