Numeric.LinearAlgebra.Algorithms

High level generic interface to common matrix computations.

Specific functions for particular base types can also be explicitly imported from Numeric.LinearAlgebra.LAPACK.

Supported types

class (Product t, Convert t, Container Vector t, Container Matrix t, Normed Matrix t, Normed Vector t, Floating t, RealOf t ~ Double) => Field t

Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can transform single precision objects using single and double).

Minimal complete definition

svd', thinSVD', sv', luPacked', luSolve', linearSolve', cholSolve', linearSolveSVD', linearSolveLS', eig', eigSH'', eigOnly, eigOnlySH, cholSH', mbCholSH', qr', hess', schur'

Instances

Field Double

Field (Complex Double)

Linear Systems

linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t

Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use linearSolveLS or linearSolveSVD. It is similar to luSolve . luPacked, but linearSolve raises an error if called on a singular system.

luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t

Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by luPacked.

cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t

Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by chol.

linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t

Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use linearSolveSVD.

linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t

Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than linearSolveLS. The effective rank of A is determined by treating as zero those singular valures which are less than eps times the largest singular value.

inv :: Field t => Matrix t -> Matrix t

Inverse of a square matrix. See also invlndet.

pinv :: Field t => Matrix t -> Matrix t

Pseudoinverse of a general matrix with default tolerance (pinvTol 1, similar to GNU-Octave).

pinvTol :: Field t => Double -> Matrix t -> Matrix t

pinvTol r computes the pseudoinverse of a matrix with tolerance tol=r*g*eps*(max rows cols), where g is the greatest singular value.

> let m = fromLists [[1,0,    0]
                    ,[0,1,    0]
                    ,[0,0,1e-10]]
  --
> pinv m
1. 0.           0.
0. 1.           0.
0. 0. 10000000000.
  --
> pinvTol 1E8 m
1. 0. 0.
0. 1. 0.
0. 0. 1.

det :: Field t => Matrix t -> t

Determinant of a square matrix. To avoid possible overflow or underflow use invlndet.

invlndet

Arguments

:: Field t

=> Matrix t

-> (Matrix t, (t, t))

(inverse, (log abs det, sign or phase of det))

Joint computation of inverse and logarithm of determinant of a square matrix.

rank :: Field t => Matrix t -> Int

Number of linearly independent rows or columns.

rcond :: Field t => Matrix t -> Double

Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.

Matrix factorizations

Singular value decomposition

svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)

Full singular value decomposition.

fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)

A version of svd which returns an appropriate diagonal matrix with the singular values.

If (u,d,v) = fullSVD m then m == u <> d <> trans v.

thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)

A version of svd which returns only the min (rows m) (cols m) singular vectors of m.

If (u,s,v) = thinSVD m then m == u <> diag s <> trans v.

compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)

Similar to thinSVD, returning only the nonzero singular values and the corresponding singular vectors.

singularValues :: Field t => Matrix t -> Vector Double

Singular values only.

leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)

Singular values and all left singular vectors.

rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)

Singular values and all right singular vectors.

Eigensystems

eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))

Eigenvalues and eigenvectors of a general square matrix.

If (s,v) = eig m then m <> v == v <> diag s

eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)

Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.

If (s,v) = eigSH m then m == v <> diag s <> ctrans v

eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)

Similar to eigSH without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.

eigenvalues :: Field t => Matrix t -> Vector (Complex Double)

Eigenvalues of a general square matrix.

eigenvaluesSH :: Field t => Matrix t -> Vector Double

Eigenvalues of a complex hermitian or real symmetric matrix.

eigenvaluesSH' :: Field t => Matrix t -> Vector Double

Similar to eigenvaluesSH without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.

geigSH'

Arguments

:: Field t

=> Matrix t

A

-> Matrix t

B

-> (Vector Double, Matrix t)

Generalized symmetric positive definite eigensystem Av = lBv, for A and B symmetric, B positive definite (conditions not checked).

QR

qr :: Field t => Matrix t -> (Matrix t, Matrix t)

QR factorization.

If (q,r) = qr m then m == q <> r, where q is unitary and r is upper triangular.

rq :: Field t => Matrix t -> (Matrix t, Matrix t)

RQ factorization.

If (r,q) = rq m then m == r <> q, where q is unitary and r is upper triangular.

Cholesky

chol :: Field t => Matrix t -> Matrix t

Cholesky factorization of a positive definite hermitian or symmetric matrix.

If c = chol m then c is upper triangular and m == ctrans c <> c.

cholSH :: Field t => Matrix t -> Matrix t

Similar to chol, without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.

mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)

Similar to cholSH, but instead of an error (e.g., caused by a matrix not positive definite) it returns Nothing.

Hessenberg

hess :: Field t => Matrix t -> (Matrix t, Matrix t)

Hessenberg factorization.

If (p,h) = hess m then m == p <> h <> ctrans p, where p is unitary and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).

Schur

schur :: Field t => Matrix t -> (Matrix t, Matrix t)

Schur factorization.

If (u,s) = schur m then m == u <> s <> ctrans u, where u is unitary and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is upper triangular in 2x2 blocks.

"Anything that the Jordan decomposition can do, the Schur decomposition can do better!" (Van Loan)

LU

lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)

Explicit LU factorization of a general matrix.

If (l,u,p,s) = lu m then m == p <> l <> u, where l is lower triangular, u is upper triangular, p is a permutation matrix and s is the signature of the permutation.

luPacked :: Field t => Matrix t -> (Matrix t, [Int])

Obtains the LU decomposition of a matrix in a compact data structure suitable for luSolve.

Matrix functions

expm :: Field t => Matrix t -> Matrix t

Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, based on a scaled Pade approximation.

sqrtm :: Field t => Matrix t -> Matrix t

Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try matFunc sqrt.

m = (2><2) [4,9
           ,0,4] :: Matrix Double

>sqrtm m
(2><2)
 [ 2.0, 2.25
 , 0.0,  2.0 ]

matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)

Generic matrix functions for diagonalizable matrices. For instance:

logm = matFunc log

Nullspace

nullspacePrec

Arguments

:: Field t

=> Double

relative tolerance in eps units (e.g., use 3 to get 3*eps)

-> Matrix t

input matrix

-> [Vector t]

list of unitary vectors spanning the nullspace

The nullspace of a matrix. See also nullspaceSVD.

nullVector :: Field t => Matrix t -> Vector t

The nullspace of a matrix, assumed to be one-dimensional, with machine precision.

nullspaceSVD

Arguments

:: Field t

=> Either Double Int

Left "numeric" zero (eg. 1*eps), or Right "theoretical" matrix rank.

-> Matrix t

input matrix m

-> (Vector Double, Matrix t)

rightSV of m

-> [Vector t]

list of unitary vectors spanning the nullspace

The nullspace of a matrix from its SVD decomposition.

orth :: Field t => Matrix t -> [Vector t]

Return an orthonormal basis of the range space of a matrix

Norms

class RealFloat (RealOf t) => Normed c t where

Methods

pnorm :: NormType -> c t -> RealOf t

Instances

Normed Vector Double

Normed Vector Float

Normed Matrix Double

Normed Matrix Float

Normed Vector (Complex Double)

Normed Vector (Complex Float)

Normed Matrix (Complex Double)

Normed Matrix (Complex Float)

data NormType

Constructors

Infinity

PNorm1

PNorm2

Frobenius

relativeError :: (Normed c t, Container c t) => c t -> c t -> Int

Approximate number of common digits in the maximum element.

Misc

eps :: Double

The machine precision of a Double: eps = 2.22044604925031e-16 (the value used by GNU-Octave).

peps :: RealFloat x => x

1 + 0.5*peps == 1, 1 + 0.6*peps /= 1

i :: Complex Double

The imaginary unit: i = 0.0 :+ 1.0

Util

haussholder :: Field a => a -> Vector a -> Matrix a

unpackQR :: Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)

unpackHess :: Field t => (Matrix t -> (Matrix t, Vector t)) -> Matrix t -> (Matrix t, Matrix t)

ranksv

Arguments

:: Double

numeric zero (e.g. 1*eps)

-> Int

maximum dimension of the matrix

-> [Double]

singular values

-> Int

rank of m

Numeric rank of a matrix from its singular values.