Numeric.LinearAlgebra.Sparse

This module exposes the user interface to the library.

Linear solvers

Iterative methods

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

Solve a linear system; uses GMRES internally as default method

Preconditioners

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

The Jacobi preconditioner is just the reciprocal of the diagonal

ilu0Pre

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 the output matrix corresponding to zero entries in the input matrix (this is called ILU(0) in the preconditioner literature)

mSsorPre

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.

Moore-Penrose pseudoinverse

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

Least-squares approximation of a rectangular system of equations.

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

r.h.s

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

r.h.s

-> m (SpVector t)

Backward substitution solver

Eigensolvers

eigsQR

Arguments

:: (MonadThrow m, MonadIO m, Num' a, Normed (SpVector a), MatrixRing (SpMatrix a), Typeable (Magnitude (SpVector a)))

=> Int

Max. # of iterations

-> Bool

Print debug information

-> SpMatrix a

Operand matrix

-> m (SpVector a)

Eigenvalues {λ_i}

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

eigsArnoldi

Arguments

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

=> Int

-> SpMatrix t

-> SpVector t

-> m (SpMatrix t, SpMatrix t, SpVector t)

Q, O, {λ_i}

`eigsArnoldi n aa b` computes at most n iterations of the Arnoldi algorithm to find a Krylov subspace of (A, b), denoted Q, along with a Hessenberg matrix of coefficients H. After that, it computes the QR decomposition of H, denoted (O, R) and the eigenvalues {λ_i} of A are listed on the diagonal of the R factor.

Matrix factorization algorithms

QR

qr

Arguments

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

=> SpMatrix a

-> m (SpMatrix a, SpMatrix a)

Q, R

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.

NB: at each iteration i we multiply the Givens matrix 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 updating two entries of M (i.e. zeroing one out and changing the other into r).

However, we must also accumulate the G_i in order to build Q, and the present formulation follows this definition closely.

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 . Implements the Doolittle algorithm, which sets the diagonal of the L matrix to ones and expects all diagonal entries of A to be nonzero. Apply pivoting (row or column permutation) to enforce a nonzero diagonal of the A matrix (the algorithm throws an appropriate exception otherwise).

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. | NB: The algorithm throws an exception if some diagonal element of A is zero.

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 -> IxRow -> IxCol -> m (SpMatrix a)

Condition number

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

Matrix condition number: computes the QR factorization and extracts the extremal eigenvalues from the R factor

Householder reflection

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

Householder reflection: a vector x uniquely defines an orthogonal (hyper)plane, i.e. an orthogonal subspace; the Householder operator reflects any point v through this subspace: v' = (I - 2 x >< x) v

Creation and conversion of sparse data

SpVector

Sparse

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

Create new SpVector using data from a Foldable (e.g. a list) in (index, value) form

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

Populate a list with SpVector contents

Dense

" from a list of entries

vr :: [Double] -> SpVector Double

Create a dense SpVector from a list of Double's

vc :: [Complex Double] -> SpVector (Complex Double)

Create a dense SpVector from a list of Complex Double's

" from/to a Vector of entries

fromVector :: Vector a -> SpVector a

Populate a SpVector with the contents of a Vector.

toVectorDense :: Num a => SpVector a -> Vector a

Populate a Vector with the entries of a SpVector, replacing the missing entries with 0

" having constant elements

constv :: Int -> a -> SpVector a

DENSE vector with constant elements

SpMatrix

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

Create new SpMatrix using data from a Foldable (e.g. a list) in (row, col, value) form

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

Populate list with SpMatrix contents

Packing and unpacking, rows or columns of a sparse matrix

", using lists as container

fromRowsL :: [SpVector a] -> SpMatrix a

Pack a list of SpVectors as rows of an SpMatrix

toRowsL :: SpMatrix a -> [SpVector a]

Unpack the rows of an SpMatrix into a list of SpVectors

fromColsL :: [SpVector a] -> SpMatrix a

Pack a list of SpVectors as columns an SpMatrix

toColsL :: SpMatrix a -> [SpVector a]

Unpack the columns of an SpMatrix into a list of SpVectors

", using Vector as container

fromRowsV :: Vector (SpVector a) -> SpMatrix a

Pack a V.Vector of SpVectors as rows of an SpMatrix

fromColsV :: Vector (SpVector a) -> SpMatrix a

Pack a V.Vector of SpVectors as columns of an SpMatrix

Manipulation of sparse data

filterSV :: (a -> Bool) -> SpVector a -> SpVector a

Filter

ifilterSV :: (Int -> a -> Bool) -> SpVector a -> SpVector a

Indexed filter

Sparsity-related

nearZero :: Epsilon a => a -> Bool

Determine if a quantity is near zero.

nearOne :: Epsilon a => a -> Bool

Is this quantity close to 1 ?

isNz :: Epsilon a => a -> Bool

Is this quantity distinguishable from 0 ?

Operators

Scaling

(.*) :: VectorSpace v => Scalar v -> v -> v

Scale a vector

(./) :: (VectorSpace v, Fractional (Scalar v)) => v -> Scalar v -> v

Scale a vector by the reciprocal of a number (e.g. for normalization)

Inner product

(<.>) :: InnerSpace v => v -> v -> Scalar v

Inner/dot product

Matrix-vector product

(#>) :: LinearVectorSpace v => MatrixType v -> v -> v

Matrix-vector action

(<#) :: LinearVectorSpace v => v -> MatrixType v -> v

Dual matrix-vector action

Matrix-matrix product

(##) :: MatrixRing m => m -> m -> m

Matrix-matrix product

(#^#) :: MatrixRing m => m -> m -> m

Matrix transpose times matrix (A^T B)

(##^) :: MatrixRing m => m -> m -> m

Matrix times matrix transpose (A B^T)

Sparsifying matrix-matrix product

(#~#) :: (MatrixRing (SpMatrix a), Epsilon a) => SpMatrix a -> SpMatrix a -> SpMatrix a

After multiplying the two matrices, all elements x for which `| x | <= eps` are removed.

(#~^#) :: (MatrixRing (SpMatrix a), Epsilon a) => SpMatrix a -> SpMatrix a -> SpMatrix a

Sparsifying A B^T

(#~#^) :: (MatrixRing (SpMatrix a), Epsilon a) => SpMatrix a -> SpMatrix a -> SpMatrix a

Sparsifying A^T B

Vector outer product

(><) :: Num a => SpVector a -> SpVector a -> SpMatrix a

Outer product

Common operations

dim :: FiniteDim f => f a -> FDSize f

Dimension (i.e. Int for SpVector, (Int, Int) for SpMatrix)

nnz :: HasData f a => f a -> Int

Number of nonzeros

spy :: (Sparse f a, Fractional b) => f a -> b

Sparsity (fraction of nonzero elements)

Vector spaces

cvx :: (VectorSpace e, Num (Scalar e)) => Scalar e -> e -> e -> e

Convex combination of two vectors (NB: 0 <= a <= 1).

Norms and normalization

norm :: (Normed v, Floating (Magnitude v)) => RealScalar v -> v -> Magnitude v

Lp norm (p > 0)

norm2 :: (Normed v, Floating (Magnitude v)) => v -> Magnitude v

Euclidean (L2) norm

norm2' :: (Normed v, Floating (Scalar v)) => v -> Scalar v

Euclidean (L2) norm; returns a Complex (norm :+ 0) for Complex-valued vectors

normalize :: Normed v => RealScalar v -> v -> v

Normalize w.r.t. Lp norm

normalize2 :: Normed v => v -> v

Normalize w.r.t. L2 norm

normalize2' :: (Normed v, Floating (Scalar v)) => v -> v

Normalize w.r.t. norm2' instead of norm2

norm1 :: Normed v => v -> Magnitude v

L1 norm

hilbertDistSq :: InnerSpace v => v -> v -> Scalar v

`hilbertDistSq x y = || x - y ||^2` computes the squared L2 distance between two vectors

Matrix-related

transpose :: MatrixRing m => m -> m

Matrix transpose (Hermitian conjugate in the Complex case)

trace :: Num b => SpMatrix b -> b

Matrix trace

normFrobenius :: MatrixRing m => m -> MatrixNorm m

Frobenius norm

Pretty-printing

prd :: PrintDense a => a -> IO ()

Pretty-print with a descriptive header

prd0 :: PrintDense a => a -> IO ()

Pretty-print with no header

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 the function throws an exception and terminates). The process ends by either hitting an iteration budget or by relative convergence, whichever happens first. After the iterations stop, the function then assesses the final state with a predicate qfinal (e.g. for comparing the final state with a known one; if this is not available, 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

print function for type b

Fields

• numIterationsMax :: Int

  Max.# of iterations

• printDebugInfo :: Bool

  Print iteration info to stdout

• iterationView :: a -> b

  Project state to a type 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

Internal

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

Iterative methods for linear systems

Constructors

GMRES_

Generalized Minimal RESidual

CGNE_

Conjugate Gradient on the Normal Equations

BCG_

BiConjugate Gradient

CGS_

Conjugate Gradient Squared

BICGSTAB_

BiConjugate Gradient Stabilized

Instances

Eq LinSolveMethod

Methods

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

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

Show LinSolveMethod

Methods

showsPrec :: Int -> LinSolveMethod -> ShowS

show :: LinSolveMethod -> String

showList :: [LinSolveMethod] -> ShowS

Exceptions

data PartialFunctionError

Errors associated with partial functions

Instances

Eq PartialFunctionError

Methods

(==) :: PartialFunctionError -> PartialFunctionError -> Bool

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

Show PartialFunctionError

Methods

showsPrec :: Int -> PartialFunctionError -> ShowS

show :: PartialFunctionError -> String

showList :: [PartialFunctionError] -> ShowS

Exception PartialFunctionError

Methods

toException :: PartialFunctionError -> SomeException

fromException :: SomeException -> Maybe PartialFunctionError

displayException :: PartialFunctionError -> String

data InputError

Input errors

Instances

Eq InputError

Methods

(==) :: InputError -> InputError -> Bool

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

Show InputError

Methods

showsPrec :: Int -> InputError -> ShowS

show :: InputError -> String

showList :: [InputError] -> ShowS

Exception InputError

Methods

toException :: InputError -> SomeException

fromException :: SomeException -> Maybe InputError

displayException :: InputError -> String

data OutOfBoundsIndexError i

Out of bounds index errors

Instances

Eq i => Eq (OutOfBoundsIndexError i)

Methods

(==) :: OutOfBoundsIndexError i -> OutOfBoundsIndexError i -> Bool

(/=) :: OutOfBoundsIndexError i -> OutOfBoundsIndexError i -> Bool

Show i => Show (OutOfBoundsIndexError i)

Methods

showsPrec :: Int -> OutOfBoundsIndexError i -> ShowS

show :: OutOfBoundsIndexError i -> String

showList :: [OutOfBoundsIndexError i] -> ShowS

(Show i, Typeable * i) => Exception (OutOfBoundsIndexError i)

Methods

toException :: OutOfBoundsIndexError i -> SomeException

fromException :: SomeException -> Maybe (OutOfBoundsIndexError i)

displayException :: OutOfBoundsIndexError i -> String

data OperandSizeMismatch

Operand size mismatch errors

Instances

Eq OperandSizeMismatch

Methods

(==) :: OperandSizeMismatch -> OperandSizeMismatch -> Bool

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

Show OperandSizeMismatch

Methods

showsPrec :: Int -> OperandSizeMismatch -> ShowS

show :: OperandSizeMismatch -> String

showList :: [OperandSizeMismatch] -> ShowS

Exception OperandSizeMismatch

Methods

toException :: OperandSizeMismatch -> SomeException

fromException :: SomeException -> Maybe OperandSizeMismatch

displayException :: OperandSizeMismatch -> String

data MatrixException i

Matrix exceptions

Instances

Eq i => Eq (MatrixException i)

Methods

(==) :: MatrixException i -> MatrixException i -> Bool

(/=) :: MatrixException i -> MatrixException i -> Bool

Show i => Show (MatrixException i)

Methods

showsPrec :: Int -> MatrixException i -> ShowS

show :: MatrixException i -> String

showList :: [MatrixException i] -> ShowS

(Show i, Typeable * i) => Exception (MatrixException i)

Methods

toException :: MatrixException i -> SomeException

fromException :: SomeException -> Maybe (MatrixException i)

displayException :: MatrixException i -> String

data IterationException a

Numerical iteration errors

Instances

Eq a => Eq (IterationException a)

Methods

(==) :: IterationException a -> IterationException a -> Bool

(/=) :: IterationException a -> IterationException a -> Bool

Show a => Show (IterationException a)

Methods

showsPrec :: Int -> IterationException a -> ShowS

show :: IterationException a -> String

showList :: [IterationException a] -> ShowS

(Show a, Typeable * a) => Exception (IterationException a)

Methods

toException :: IterationException a -> SomeException

fromException :: SomeException -> Maybe (IterationException a)

displayException :: IterationException a -> String

Orphan instances

LinearSystem (SpVector Double)

<> uses the GMRES method as default

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