Numeric.LinearAlgebra.Sparse

This module exposes the high-level functionality of the library.

Linear solvers

Iterative methods

linSolve0 :: ((~#) * * (Scalar (SpVector t)) t, (~#) * * (MatrixType (SpVector t)) (SpMatrix t), Show t, Epsilon t, MonadThrow m, MonadIO m, LinearVectorSpace (SpVector t), Typeable * (Magnitude (SpVector t)), Typeable * t, Normed (SpVector t), Elt t) => LinSolveMethod -> SpMatrix t -> SpVector t -> SpVector t -> m (SpVector t)

Interface method to individual linear solvers

data LinSolveMethod

Constructors

GMRES_

CGNE_

BCG_

CGS_

BICGSTAB_

Instances

Eq LinSolveMethod

Methods

(==) :: LinSolveMethod -> LinSolveMethod -> Bool

(/=) :: LinSolveMethod -> LinSolveMethod -> Bool

Show LinSolveMethod

Methods

showsPrec :: Int -> LinSolveMethod -> ShowS

show :: LinSolveMethod -> String

showList :: [LinSolveMethod] -> ShowS

(<\>) :: (LinearSystem v, MonadIO m, MonadThrow m) => MatrixType v -> v -> m v

Moore-Penrose pseudoinverse

pinv :: (MatrixType v ~ SpMatrix a, LinearSystem v, Epsilon a, MonadThrow m, MonadIO m) => SpMatrix a -> v -> m v

Least-squares approximation of a rectangular system of equaitons. Uses \ for the linear solve

Preconditioners

jacobiPre :: Fractional a => SpMatrix a -> SpMatrix a

The Jacobi preconditioner is just the reciprocal of the diagonal

ilu0

Arguments

:: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t), Epsilon t, MonadThrow m)

=> SpMatrix t

-> m (SpMatrix t, SpMatrix t)

L, U (with holes)

Used for Incomplete LU : remove entries in m corresponding to zero entries in m2 (this is called ILU(0) in the preconditioner literature)

mSsor

Arguments

:: (MatrixRing (SpMatrix b), Fractional b)

=> SpMatrix b

-> b

relaxation factor

-> (SpMatrix b, SpMatrix b)

Left, right factors

Symmetric Successive Over-Relaxation. `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner. When ω = 1, the SOR reduces to Gauss-Seidel; when ω > 1 and ω < 1, it corresponds to over-relaxation and under-relaxation, respectively.

Direct methods

luSolve

Arguments

:: (Scalar (SpVector t) ~ t, MonadThrow m, Elt t, InnerSpace (SpVector t), Epsilon t)

=> SpMatrix t

Lower triangular

-> SpMatrix t

Upper triangular

-> SpVector t

r.h.s.

-> m (SpVector t)

Direct solver based on a triangular factorization of the system matrix.

Forward substitution

triLowerSolve

Arguments

:: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t), Epsilon t, MonadThrow m)

=> SpMatrix t

Lower triangular

-> SpVector t

-> m (SpVector t)

Forward substitution solver

Backward substitution

triUpperSolve

Arguments

:: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t), Epsilon t, MonadThrow m)

=> SpMatrix t

Upper triangular

-> SpVector t

-> m (SpVector t)

Backward substitution solver

Eigensolvers

eigsQR

Arguments

:: (MonadThrow m, MonadIO m, Elt a, Normed (SpVector a), MatrixRing (SpMatrix a), Epsilon a, Typeable (Magnitude (SpVector a)), Typeable a, Show a)

=> Int

-> Bool

Print debug information

-> SpMatrix a

Operand matrix

-> m (SpVector a)

Eigenvalues

`eigsQR n mm` performs n iterations of the QR algorithm on matrix mm, and returns a SpVector containing all eigenvalues

eigRayleigh :: ((~#) * * (MatrixType b) (SpMatrix (Scalar b)), Show b, Show (Scalar b), Floating (Scalar b), LinearSystem b, MonadIO m, MonadThrow m, Normed b, Typeable * (Magnitude b), Typeable * b, Typeable * (Scalar b)) => Int -> Bool -> (b -> IO ()) -> SpMatrix (Scalar b) -> (b, Scalar b) -> m (b, Scalar b)

`eigsRayleigh n mm` performs n iterations of the Rayleigh algorithm on matrix mm and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair.

Matrix factorization algorithms

QR

qr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) => SpMatrix a -> m (SpMatrix a, SpMatrix a)

Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.

LU

lu

Arguments

:: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t), Epsilon t, MonadThrow m)

=> SpMatrix t

-> m (SpMatrix t, SpMatrix t)

L, U

Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A

Cholesky

chol

Arguments

:: (Elt a, Epsilon a, MonadThrow m)

=> SpMatrix a

-> m (SpMatrix a)

L

Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A . This is an implementation of the Cholesky–Banachiewicz algorithm, i.e. proceeding row by row from the upper-left corner.

Arnoldi iteration

arnoldi

Arguments

:: (MatrixType (SpVector a) ~ SpMatrix a, V (SpVector a), Scalar (SpVector a) ~ a, Epsilon a, MonadThrow m)

=> SpMatrix a

System matrix

-> SpVector a

r.h.s.

-> Int

Max. # of iterations

-> m (SpMatrix a, SpMatrix a)

Q, H

Given a matrix A, a vector b and a positive integer n, this procedure finds the basis of an order n Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q. At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.

Matrix partitioning

diagPartitions

Arguments

:: SpMatrix a

-> (SpMatrix a, SpMatrix a, SpMatrix a)

Subdiagonal, diagonal, superdiagonal partitions

Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts

Utilities

Givens' rotation

givens :: (Elt a, MonadThrow m) => SpMatrix a -> Int -> Int -> m (SpMatrix a)

Givens method, row version: choose other row index i' s.t. i' is : * below the diagonal * corresponding element is nonzero

QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is non-zero but A has zeros in row k for all columns less than j.

NB: the current version is quite inefficient in that: 1. the Givens' matrix G_i is different from Identity only in 4 entries 2. at each iteration i we multiply G_i by the previous partial result M. Since this corresponds to a rotation, and the givensCoef function already computes the value of the resulting non-zero component (output r), `G_i ## M` can be simplified by just changing two entries of M (i.e. zeroing one out and changing the other into r).

Condition number

conditionNumberSM :: (MonadThrow m, MatrixRing (SpMatrix a), Num' a, Typeable a) => SpMatrix a -> m a

uses the R matrix from the QR factorization

Householder reflection

hhMat :: Num a => a -> SpVector a -> SpMatrix a

hhRefl :: Num a => SpVector a -> SpMatrix a

Householder reflection: a vector x uniquely defines an orthogonal plane; the Householder operator reflects any point v with respect to this plane: v' = (I - 2 x >< x) v

From/to SpVector

fromListSV :: Foldable t => Int -> t (Int, a) -> SpVector a

toListSV :: SpVector a -> [(Key, a)]

From/to SpMatrix

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

Create new SpMatrix using data from list (row, col, value)

toListSM :: SpMatrix t -> [(IxRow, IxCol, t)]

Populate list with SpMatrix contents

Iteration combinators

untilConvergedG0

Arguments

:: (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s)

=> String

-> IterationConfig s v

-> v

Known value

-> (s -> s)

-> s

-> m s

untilConvergedG0 is a special case of untilConvergedG that assesses convergence based on the L2 distance to a known solution xKnown

untilConvergedG :: (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s) => String -> IterationConfig s v -> (v -> Bool) -> (s -> s) -> s -> m s

This function makes some default choices on the modifyInspectGuarded machinery: convergence is assessed using the squared L2 distance between consecutive states, and divergence is detected when this function is increasing between pairs of measurements.

untilConvergedGM :: (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s) => String -> IterationConfig s v -> (v -> Bool) -> (s -> m s) -> s -> m s

", monadic version

modifyInspectGuarded

Arguments

:: (MonadThrow m, MonadIO m, Typeable s, Typeable a, Show s, Show a)

=> String

Calling function name

-> IterationConfig s v

Configuration

-> ([v] -> a)

State summary array projection

-> (a -> Bool)

Convergence criterion

-> (a -> a -> Bool)

Divergence criterion

-> (v -> Bool)

Final state acceptance criterion

-> (s -> s)

State evolution

-> s

Initial state

-> m s

Final state

modifyInspectGuarded is a high-order abstraction of a numerical iterative process. It accumulates a rolling window of 3 states and compares a summary q of the latest 2 with that of the previous two in order to assess divergence (e.g. if `q latest2 > q prev2` then it). The process ends when either we hit an iteration budget or relative convergence is verified. The function then assesses the final state with a predicate qfinal (e.g. against a known solution; if this is not known, the user can just supply `const True`)

modifyInspectGuardedM :: (MonadThrow m, MonadIO m, Typeable s, Show s, Typeable a, Show a) => String -> IterationConfig s v -> ([v] -> a) -> (a -> Bool) -> (a -> a -> Bool) -> (v -> Bool) -> (s -> m s) -> s -> m s

", monadic version

data IterationConfig a b

Constructors

IterConf

Fields

• numIterationsMax :: Int
• printDebugInfo :: Bool
• iterationView :: a -> b
• printDebugIO :: b -> IO ()

Instances

Show (IterationConfig a b)

Methods

showsPrec :: Int -> IterationConfig a b -> ShowS

show :: IterationConfig a b -> String

showList :: [IterationConfig a b] -> ShowS

modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s

modifyUntilM :: MonadState s m => (s -> Bool) -> (s -> m s) -> m s

Orphan instances

LinearSystem (SpVector Double)

Methods

(<\>) :: (MonadIO m, MonadThrow m) => MatrixType (SpVector Double) -> SpVector Double -> m (SpVector Double)

LinearSystem (SpVector (Complex Double))

Methods

(<\>) :: (MonadIO m, MonadThrow m) => MatrixType (SpVector (Complex Double)) -> SpVector (Complex Double) -> m (SpVector (Complex Double))