sparse-linear-algebra-0.3.1: Numerical computing in native Haskell

Numeric.LinearAlgebra.Sparse

Description

This module exposes the user interface to the library.

Synopsis

# Linear solvers

## Iterative methods

Arguments

 :: (LinearSystem v, MonadIO m, MonadThrow m) => MatrixType v System matrix -> v Right-hand side -> m v Result

Solve a linear system; uses GMRES internally as default method

## Moore-Penrose pseudoinverse

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

Least-squares approximation of a rectangular system of equations.

## Direct methods

Arguments

 :: (Scalar (SpVector t) ~ t, MonadThrow m, Elt t, InnerSpace (SpVector t), Epsilon t, PrintDense (SpVector t), MonadIO m) => 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 :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t), PrintDense (SpVector t), Epsilon t, MonadThrow m, MonadIO m) => SpMatrix t -> SpVector t -> m (SpVector t) Source #

Forward substitution solver

### Backward substitution

triUpperSolve :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t), PrintDense (SpVector t), Epsilon t, MonadThrow m, MonadIO m) => SpMatrix t -> SpVector t -> m (SpVector t) Source #

Backward substitution solver

# Eigensolvers

Arguments

 :: (MonadThrow m, MonadIO m, Num' a, Normed (SpVector a), MatrixRing (SpMatrix a), Typeable (Magnitude (SpVector a)), PrintDense (SpVector a), PrintDense (SpMatrix 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.

Arguments

 :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t, Elt t, V (SpVector t), Epsilon t, PrintDense (SpMatrix t), MatrixRing (SpMatrix t), MonadThrow m, MonadIO 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

Arguments

 :: (Elt a, MatrixRing (SpMatrix a), PrintDense (SpMatrix a), Epsilon a, MonadThrow m, MonadIO 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

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

Arguments

 :: (Elt a, Epsilon a, MonadThrow m, MonadIO m, PrintDense (SpMatrix a)) => 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

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 ith iteration, it finds (i + 1) coefficients (the ith column of the Hessenberg matrix H) and the (i + 1)th Krylov vector.

# Utilities

## Givens' rotation

givens :: (Elt a, MonadThrow m) => SpMatrix a -> IxRow -> IxCol -> m (SpMatrix a) Source #

Givens matrix : a planar rotation embedded in R^n.

>>> aa = fromListSM
>>> g <- givens aa 1 0


Row version of the method: given a matrix element below the diagonal, indexed (i,j), choose a row index i' that is below the diagonal as well and distinct from i such that the corresponding element is nonzero.

To zero out entry A(i, j) we must find row i' such that A(i', j) is non-zero but A has zeros in row i' for all column indices < j.

NB: The Givens' matrix differs from Identity in 4 entries

NB2: The form of a Complex rotation matrix in R^2 is as follows (* indicates complex conjugation):

    ( c    s )
G =(        )
( -s*  c*)


## Condition number

conditionNumberSM :: (MonadThrow m, MonadIO m, MatrixRing (SpMatrix a), PrintDense (SpMatrix a), Num' a) => SpMatrix a -> m a Source #

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

## Householder reflection

hhRefl :: (Num a, AdditiveGroup a) => SpVector a -> SpMatrix a Source #

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

## Matrix partitioning

Arguments

 :: SpMatrix a -> (SpMatrix a, SpMatrix a, SpMatrix a) Subdiagonal, diagonal, superdiagonal partitions

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

# Creation and conversion of sparse data

## SpVector

### Sparse

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

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

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

Populate a list with SpVector contents

### Dense

#### " from a list of entries

vr :: [Double] -> SpVector Double Source #

Create a dense SpVector from a list of Double's

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

#### " from/to a Vector of entries

Populate a SpVector with the contents of a Vector.

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

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

#### " having constant elements

constv :: Int -> a -> SpVector a Source #

DENSE vector with constant elements

## SpMatrix

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

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

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

Populate list with SpMatrix contents

## Packing and unpacking, rows or columns of a sparse matrix

### ", using lists as container

fromRowsL :: [SpVector a] -> SpMatrix a Source #

Pack a list of SpVectors as rows of an SpMatrix

toRowsL :: SpMatrix a -> [SpVector a] Source #

Unpack the rows of an SpMatrix into a list of SpVectors

fromColsL :: [SpVector a] -> SpMatrix a Source #

Pack a list of SpVectors as columns an SpMatrix

toColsL :: SpMatrix a -> [SpVector a] Source #

Unpack the columns of an SpMatrix into a list of SpVectors

### ", using Vector as container

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

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

## Block operations

(-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a Source #

Vertical stacking of matrix blocks

(-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a Source #

Horizontal stacking of matrix blocks

fromBlocksDiag :: [SpMatrix a] -> SpMatrix a Source #

Assemble a block-diagonal square matrix from a list of square matrices, arranging these along the main diagonal

## Special matrices

eye :: Num a => Int -> SpMatrix a Source #

eye n : identity matrix of rank n

mkDiagonal :: Int -> [a] -> SpMatrix a Source #

mkDiagonal n ll : create a diagonal matrix of size n from a list ll of elements

mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a Source #

mkSubDiagonal n o xx creates a square SpMatrix of size n with xx on the oth subdiagonal

permutationSM :: Num a => Int -> [IxRow] -> SpMatrix a Source #

Permutation matrix from a (possibly incomplete) list of row swaps starting from row 0 e.g. permutationSM 5 [1,3] first swaps rows (0, 1) and then rows (1, 3) :

>>> prd (permutationSM 5 [1,3] :: SpMatrix Double)

( 5 rows, 5 columns ) , 5 NZ ( density 20.000 % )

_      , 1.00   , _      , _      , _
_      , _      , _      , 1.00   , _
_      , _      , 1.00   , _      , _
1.00   , _      , _      , _      , _
_      , _      , _      , _      , 1.00


permutPairsSM :: Num a => Int -> [(IxRow, IxRow)] -> SpMatrix a Source #

Permutation matrix from a (possibly incomplete) list of row pair swaps e.g. permutPairs 5 [(2,4)] swaps rows 2 and 4 :

>>> prd (permutPairsSM 5 [(2,4)] :: SpMatrix Double)

( 5 rows, 5 columns ) , 5 NZ ( density 20.000 % )

1.00   , _      , _      , _      , _
_      , 1.00   , _      , _      , _
_      , _      , _      , _      , 1.00
_      , _      , _      , 1.00   , _
_      , _      , 1.00   , _      , _


# Predicates

isOrthogonalSM :: (Eq a, Epsilon a, MatrixRing (SpMatrix a)) => SpMatrix a -> Bool Source #

Is the matrix orthogonal? i.e. Q^t ## Q == I

Is the matrix diagonal?

# Manipulation of sparse data

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

Filter

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

Indexed filter

# Sparsity-related predicates

nearZero :: Epsilon a => a -> Bool Source #

Determine if a quantity is near zero.

nearOne :: Epsilon a => a -> Bool Source #

Is this quantity close to 1 ?

isNz :: Epsilon a => a -> Bool Source #

Is this quantity distinguishable from 0 ?

# Operators

## Scaling

(.*) :: VectorSpace v => Scalar v -> v -> v infixr 7 Source #

Scale a vector

(./) :: (VectorSpace v, s ~ Scalar v, Fractional s) => v -> s -> v infixr 7 Source #

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

## Inner product

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

Inner/dot product

## Matrix-vector product

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

Matrix-vector action

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

Dual matrix-vector action

## Matrix-matrix product

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

Matrix-matrix product

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

Matrix transpose times matrix (A^T B)

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

Matrix times matrix transpose (A B^T)

### Sparsifying matrix-matrix product

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

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

Sparsifying A B^T

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

Sparsifying A^T B

## Vector outer product

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

Outer product

# Common operations

dim :: FiniteDim f => f -> FDSize f Source #

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

nnz :: HasData f => f -> Int Source #

Number of nonzeros

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

Sparsity (fraction of nonzero elements)

## Vector spaces

cvx :: VectorSpace v => Scalar v -> v -> v -> v Source #

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

### Norms and normalization

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

Lp norm (p > 0)

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

Euclidean (L2) norm

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

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

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

Normalize w.r.t. Lp norm

normalize2 :: Normed v => v -> v Source #

Normalize w.r.t. L2 norm

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

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

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

L1 norm

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

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

## Matrix-related

transpose :: MatrixRing m => m -> m Source #

Matrix transpose (Hermitian conjugate in the Complex case)

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

Matrix trace

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

Frobenius norm

# Pretty-printing

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

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

# Iteration combinators

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

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

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

data IterationConfig a b Source #

Constructors

 IterConf print function for type b FieldsnumIterationsMax :: IntMax.# of iterationsprintDebugInfo :: BoolPrint iteration info to stdout iterationView :: a -> bProject state to a type b`printDebugIO :: b -> IO ()

Instances

 Show (IterationConfig a b) Source # MethodsshowsPrec :: 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 Source #

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

# Internal

linSolve0 :: ((~#) * * (Scalar b) b, (~#) * * (Magnitude b) (RealScalar b), (~#) * * (RealScalar b) (Scalar b), Floating (Scalar b), Epsilon b, Elt b, PrintDense (SpMatrix b), PrintDense (SpVector b), MatrixRing (SpMatrix b), Normed b, MonadThrow m, MonadIO m, Typeable * (Magnitude b), Typeable * b, Show b) => LinSolveMethod -> SpMatrix b -> SpVector b -> SpVector b -> m (SpVector b) Source #

Interface method to individual linear solvers

Iterative methods for linear systems

Constructors

Instances

 Source # Methods Source # MethodsshowList :: [LinSolveMethod] -> ShowS #

# Exceptions

Errors associated with partial functions

Instances

 Source # Methods Source # Methods Source # Methods

Input errors

Instances

 Source # Methods Source # MethodsshowList :: [InputError] -> ShowS # Source # Methods

Out of bounds index errors

Instances

 Eq i => Eq (OutOfBoundsIndexError i) Source # Methods Source # MethodsshowList :: [OutOfBoundsIndexError i] -> ShowS # (Show i, Typeable * i) => Exception (OutOfBoundsIndexError i) Source # Methods

Operand size mismatch errors

Instances

 Source # Methods Source # Methods Source # Methods

Matrix exceptions

Instances

 Eq i => Eq (MatrixException i) Source # Methods(==) :: MatrixException i -> MatrixException i -> Bool #(/=) :: MatrixException i -> MatrixException i -> Bool # Show i => Show (MatrixException i) Source # MethodsshowList :: [MatrixException i] -> ShowS # (Show i, Typeable * i) => Exception (MatrixException i) Source # Methods

Numerical iteration errors

Instances

 Eq a => Eq (IterationException a) Source # Methods Show a => Show (IterationException a) Source # MethodsshowList :: [IterationException a] -> ShowS # (Show a, Typeable * a) => Exception (IterationException a) Source # Methods

# Orphan instances

 Source # <> uses the GMRES method as default Methods(<\>) :: (MonadIO m, MonadThrow m) => MatrixType (SpVector Double) -> SpVector Double -> m (SpVector Double) Source #