Numeric.LinearAlgebra.HMatrix

Basic types and data processing

module Numeric.LinearAlgebra.Data

Arithmetic and numeric classes

The standard numeric classes are defined elementwise:

>>> vector [1,2,3] * vector [3,0,-2]
fromList [3.0,0.0,-6.0]

>>> matrix 3 [1..9] * ident 3
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 5.0, 0.0
 , 0.0, 0.0, 9.0 ]

In arithmetic operations single-element vectors and matrices (created from numeric literals or using scalar) automatically expand to match the dimensions of the other operand:

>>> 5 + 2*ident 3 :: Matrix Double
(3><3)
 [ 7.0, 5.0, 5.0
 , 5.0, 7.0, 5.0
 , 5.0, 5.0, 7.0 ]

>>> matrix 3 [1..9] + matrix 1 [10,20,30]
(3><3)
 [ 11.0, 12.0, 13.0
 , 24.0, 25.0, 26.0
 , 37.0, 38.0, 39.0 ]

Products

dot

dot :: Numeric t => Vector t -> Vector t -> t

<·> :: Numeric t => Vector t -> Vector t -> t infixr 8

infix synonym for dot

>>> vector [1,2,3,4] <·> vector [-2,0,1,1]
5.0

>>> let 𝑖 = 0:+1 :: ℂ
>>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
2.0 :+ 0.0

(the dot symbol "·" is obtained by Alt-Gr .)

matrix-vector

app :: Numeric t => Matrix t -> Vector t -> Vector t

dense matrix-vector product

(#>) :: Numeric t => Matrix t -> Vector t -> Vector t infixr 8

infix synonym for app

>>> let m = (2><3) [1..]
>>> m
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]

>>> let v = vector [10,20,30]

>>> m #> v
fromList [140.0,320.0]

(!#>) :: GMatrix -> Vector Double -> Vector Double infixr 8

general matrix - vector product

>>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]
>>> m !#> vector [1..2000]
fromList [1000.0,4000.0]

matrix-matrix

mul :: Numeric t => Matrix t -> Matrix t -> Matrix t

dense matrix product

(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t infixr 8

infix synonym of mul

>>> let a = (3><5) [1..]
>>> a
(3><5)
 [  1.0,  2.0,  3.0,  4.0,  5.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0 ]

>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
>>> b
(5><2)
 [  1.0, 3.0
 ,  0.0, 2.0
 , -1.0, 5.0
 ,  7.0, 7.0
 ,  6.0, 0.0 ]

>>> a <> b
(3><2)
 [  56.0,  50.0
 , 121.0, 135.0
 , 186.0, 220.0 ]

The matrix product is also implemented in the Data.Monoid instance, where single-element matrices (created from numeric literals or using scalar) are used for scaling.

>>> import Data.Monoid as M
>>> let m = matrix 3 [1..6]
>>> m M.<> 2 M.<> diagl[0.5,1,0]
(2><3)
 [ 1.0,  4.0, 0.0
 , 4.0, 10.0, 0.0 ]

mconcat uses optimiseMult to get the optimal association order.

other

outer :: Product t => Vector t -> Vector t -> Matrix t

Outer product of two vectors.

>>> fromList [1,2,3] `outer` fromList [5,2,3]
(3><3)
 [  5.0, 2.0, 3.0
 , 10.0, 4.0, 6.0
 , 15.0, 6.0, 9.0 ]

kronecker :: Product t => Matrix t -> Matrix t -> Matrix t

Kronecker product of two matrices.

m1=(2><3)
 [ 1.0,  2.0, 0.0
 , 0.0, -1.0, 3.0 ]
m2=(4><3)
 [  1.0,  2.0,  3.0
 ,  4.0,  5.0,  6.0
 ,  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0 ]

>>> kronecker m1 m2
(8><9)
 [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
 ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
 ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
 , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
 ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
 ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
 ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
 ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]

cross :: Vector Double -> Vector Double -> Vector Double

cross product (for three-element real vectors)

scale :: Container c e => e -> c e -> c e

multiplication by scalar

sumElements :: Container c e => c e -> e

the sum of elements

prodElements :: Container c e => c e -> e

the product of elements

Linear Systems

(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t infixl 7

Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSVD)

a = (3><2)
 [ 1.0,  2.0
 , 2.0,  4.0
 , 2.0, -1.0 ]

v = vector [13.0,27.0,1.0]

>>> let x = a <\> v
>>> x
fromList [3.0799999999999996,5.159999999999999]

>>> a #> x
fromList [13.399999999999999,26.799999999999997,1.0]

It also admits multiple right-hand sides stored as columns in a matrix.

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

Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use linearSolveLS or linearSolveSVD.

a = (2><2)
 [ 1.0, 2.0
 , 3.0, 5.0 ]

b = (2><3)
 [  6.0, 1.0, 10.0
 , 15.0, 3.0, 26.0 ]

>>> linearSolve a b
Just (2><3)
 [ -1.4802973661668753e-15,     0.9999999999999997, 1.999999999999997
 ,       3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ]

>>> let Just x = it
>>> disp 5 x
2x3
-0.00000  1.00000  2.00000
 3.00000  0.00000  4.00000

>>> a <> x
(2><3)
 [  6.0, 1.0, 10.0
 , 15.0, 3.0, 26.0 ]

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.

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.

cgSolve

Arguments

:: Bool

is symmetric

-> GMatrix

coefficient matrix

-> Vector Double

right-hand side

-> Vector Double

solution

cgSolve'

Arguments

:: Bool

symmetric

-> R

relative tolerance for the residual (e.g. 1E-4)

-> R

relative tolerance for δx (e.g. 1E-3)

-> Int

maximum number of iterations

-> GMatrix

coefficient matrix

-> V

initial solution

-> V

right-hand side

-> [CGState]

solution

Inverse and pseudoinverse

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.

m = (3><3) [ 1, 0,    0
           , 0, 1,    0
           , 0, 0, 1e-10] :: Matrix Double

>>> pinv m
1. 0.           0.
0. 1.           0.
0. 0. 10000000000.

>>> pinvTol 1E8 m
1. 0. 0.
0. 1. 0.
0. 0. 1.

Determinant and rank

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

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

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

Number of linearly independent rows or columns. See also ranksv

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.

Norms

class Normed a where

Methods

norm_0 :: a -> ℝ

norm_1 :: a -> ℝ

norm_2 :: a -> ℝ

norm_Inf :: a -> ℝ

Instances

Normed (Vector ℂ)

Normed (Vector ℝ)

Normed (Matrix ℂ)

Normed (Matrix ℝ)

norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ

norm_nuclear :: Field t => Matrix t -> ℝ

Nullspace and range

orth :: Field t => Matrix t -> Matrix t

return an orthonormal basis of the range space of a matrix. See also orthSVD.

nullspace :: Field t => Matrix t -> Matrix t

return an orthonormal basis of the null space of a matrix. See also nullspaceSVD.

null1 :: Matrix Double -> Vector Double

solution of overconstrained homogeneous linear system

null1sym :: Matrix Double -> Vector Double

solution of overconstrained homogeneous symmetric linear system

SVD

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

Full singular value decomposition.

a = (5><3)
 [  1.0,  2.0,  3.0
 ,  4.0,  5.0,  6.0
 ,  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0
 , 13.0, 14.0, 15.0 ] :: Matrix Double

>>> let (u,s,v) = svd a

>>> disp 3 u
5x5
-0.101   0.768   0.614   0.028  -0.149
-0.249   0.488  -0.503   0.172   0.646
-0.396   0.208  -0.405  -0.660  -0.449
-0.543  -0.072  -0.140   0.693  -0.447
-0.690  -0.352   0.433  -0.233   0.398

>>> s
fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]

>>> disp 3 v
3x3
-0.519  -0.751   0.408
-0.576  -0.046  -0.816
-0.632   0.659   0.408

>>> let d = diagRect 0 s 5 3
>>> disp 3 d
5x3
35.183  0.000  0.000
 0.000  1.477  0.000
 0.000  0.000  0.000
 0.000  0.000  0.000

>>> disp 3 $ u <> d <> tr v
5x3
 1.000   2.000   3.000
 4.000   5.000   6.000
 7.000   8.000   9.000
10.000  11.000  12.000
13.000  14.000  15.000

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 <> tr v.

a = (5><3)
 [  1.0,  2.0,  3.0
 ,  4.0,  5.0,  6.0
 ,  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0
 , 13.0, 14.0, 15.0 ] :: Matrix Double

>>> let (u,s,v) = thinSVD a

>>> disp 3 u
5x3
-0.101   0.768   0.614
-0.249   0.488  -0.503
-0.396   0.208  -0.405
-0.543  -0.072  -0.140
-0.690  -0.352   0.433

>>> s
fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]

>>> disp 3 v
3x3
-0.519  -0.751   0.408
-0.576  -0.046  -0.816
-0.632   0.659   0.408

>>> disp 3 $ u <> diag s <> tr v
5x3
 1.000   2.000   3.000
 4.000   5.000   6.000
 7.000   8.000   9.000
10.000  11.000  12.000
13.000  14.000  15.000

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.

a = (5><3)
 [  1.0,  2.0,  3.0
 ,  4.0,  5.0,  6.0
 ,  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0
 , 13.0, 14.0, 15.0 ] :: Matrix Double

>>> let (u,s,v) = compactSVD a

>>> disp 3 u
5x2
-0.101   0.768
-0.249   0.488
-0.396   0.208
-0.543  -0.072
-0.690  -0.352

>>> s
fromList [35.18264833189422,1.4769076999800903]

>>> disp 3 u
5x2
-0.101   0.768
-0.249   0.488
-0.396   0.208
-0.543  -0.072
-0.690  -0.352

>>> disp 3 $ u <> diag s <> tr v
5x3
 1.000   2.000   3.000
 4.000   5.000   6.000
 7.000   8.000   9.000
10.000  11.000  12.000
13.000  14.000  15.000

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 (as columns).

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

Singular values and all right singular vectors (as columns).

Eigensystems

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

Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.

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

a = (3><3)
 [ 3, 0, -2
 , 4, 5, -1
 , 3, 1,  0 ] :: Matrix Double

>>> let (l, v) = eig a

>>> putStr . dispcf 3 . asRow $ l
1x3
1.925+1.523i  1.925-1.523i  4.151

>>> putStr . dispcf 3 $ v
3x3
-0.455+0.365i  -0.455-0.365i   0.181
        0.603          0.603  -0.978
 0.033+0.543i   0.033-0.543i  -0.104

>>> putStr . dispcf 3 $ complex a <> v
3x3
-1.432+0.010i  -1.432-0.010i   0.753
 1.160+0.918i   1.160-0.918i  -4.059
-0.763+1.096i  -0.763-1.096i  -0.433

>>> putStr . dispcf 3 $ v <> diag l
3x3
-1.432+0.010i  -1.432-0.010i   0.753
 1.160+0.918i   1.160-0.918i  -4.059
-0.763+1.096i  -0.763-1.096i  -0.433

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

Eigenvalues and eigenvectors (as columns) of a complex hermitian or real symmetric matrix, in descending order.

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

a = (3><3)
 [ 1.0, 2.0, 3.0
 , 2.0, 4.0, 5.0
 , 3.0, 5.0, 6.0 ]

>>> let (l, v) = eigSH a

>>> l
fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]

>>> disp 3 $ v <> diag l <> tr v
3x3
1.000  2.000  3.000
2.000  4.000  5.000
3.000  5.000  6.000

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 (not ordered) of a general square matrix.

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

Eigenvalues (in descending order) 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.

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

qrgr :: Field t => Int -> (Matrix t, Vector t) -> Matrix t

generate a matrix with k orthogonal columns from the output of qrRaw

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.

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

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

For diagonalizable matrices you can try matFunc sqrt:

>>> matFunc sqrt ((2><2) [1,0,0,-1])
(2><2)
 [ 1.0 :+ 0.0, 0.0 :+ 0.0
 , 0.0 :+ 0.0, 0.0 :+ 1.0 ]

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

Generic matrix functions for diagonalizable matrices. For instance:

logm = matFunc log

Correlation and convolution

corr

Arguments

:: (Container Vector t, Product t)

=> Vector t

kernel

-> Vector t

source

-> Vector t

correlation

>>> corr (fromList[1,2,3]) (fromList [1..10])
fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]

conv :: (Container Vector t, Product t, Num t) => Vector t -> Vector t -> Vector t

convolution (corr with reversed kernel and padded input, equivalent to polynomial product)

>>> conv (fromList[1,1]) (fromList [-1,1])
fromList [-1.0,0.0,1.0]

corrMin :: (Container Vector t, RealElement t, Product t) => Vector t -> Vector t -> Vector t

similar to corr, using min instead of (*)

corr2 :: Product a => Matrix a -> Matrix a -> Matrix a

2D correlation (without padding)

>>> disp 5 $ corr2 (konst 1 (3,3)) (ident 10 :: Matrix Double)
8x8
3  2  1  0  0  0  0  0
2  3  2  1  0  0  0  0
1  2  3  2  1  0  0  0
0  1  2  3  2  1  0  0
0  0  1  2  3  2  1  0
0  0  0  1  2  3  2  1
0  0  0  0  1  2  3  2
0  0  0  0  0  1  2  3

conv2

Arguments

:: (Num (Matrix a), Product a, Container Vector a)

=> Matrix a

kernel

-> Matrix a

-> Matrix a

2D convolution

>>> disp 5 $ conv2 (konst 1 (3,3)) (ident 10 :: Matrix Double)
12x12
1  1  1  0  0  0  0  0  0  0  0  0
1  2  2  1  0  0  0  0  0  0  0  0
1  2  3  2  1  0  0  0  0  0  0  0
0  1  2  3  2  1  0  0  0  0  0  0
0  0  1  2  3  2  1  0  0  0  0  0
0  0  0  1  2  3  2  1  0  0  0  0
0  0  0  0  1  2  3  2  1  0  0  0
0  0  0  0  0  1  2  3  2  1  0  0
0  0  0  0  0  0  1  2  3  2  1  0
0  0  0  0  0  0  0  1  2  3  2  1
0  0  0  0  0  0  0  0  1  2  2  1
0  0  0  0  0  0  0  0  0  1  1  1

Random arrays

type Seed = Int

data RandDist

Constructors

Uniform

uniform distribution in [0,1)

Gaussian

normal distribution with mean zero and standard deviation one

Instances

Enum RandDist

randomVector

Arguments

:: Seed

-> RandDist

distribution

-> Int

vector size

-> Vector Double

Obtains a vector of pseudorandom elements (use randomIO to get a random seed).

rand :: Int -> Int -> IO (Matrix Double)

pseudorandom matrix with uniform elements between 0 and 1

randn :: Int -> Int -> IO (Matrix Double)

pseudorandom matrix with normal elements

>>> disp 3 =<< randn 3 5
3x5
0.386  -1.141   0.491  -0.510   1.512
0.069  -0.919   1.022  -0.181   0.745
0.313  -0.670  -0.097  -1.575  -0.583

gaussianSample

Arguments

:: Seed

-> Int

number of rows

-> Vector Double

mean vector

-> Matrix Double

covariance matrix

-> Matrix Double

result

Obtains a matrix whose rows are pseudorandom samples from a multivariate Gaussian distribution.

uniformSample

Arguments

:: Seed

-> Int

number of rows

-> [(Double, Double)]

ranges for each column

-> Matrix Double

result

Obtains a matrix whose rows are pseudorandom samples from a multivariate uniform distribution.

Misc

meanCov :: Matrix Double -> (Vector Double, Matrix Double)

Compute mean vector and covariance matrix of the rows of a matrix.

>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
(fromList [4.010341078059521,5.0197204699640405],
(2><2)
 [     1.9862461923890056, -1.0127225830525157e-2
 , -1.0127225830525157e-2,     3.0373954915729318 ])

rowOuters :: Matrix Double -> Matrix Double -> Matrix Double

outer products of rows

>>> a
(3><2)
 [   1.0,   2.0
 ,  10.0,  20.0
 , 100.0, 200.0 ]
>>> b
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]

>>> rowOuters a (b ||| 1)
(3><8)
 [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
 ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
 , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]

peps :: RealFloat x => x

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

relativeError :: (Normed c t, Num (c t)) => NormType -> c t -> c t -> Double

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

optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t

udot :: Product e => Vector e -> Vector e -> e

unconjugated dot product

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

-> Matrix t

nullspace

The nullspace of a matrix from its precomputed SVD decomposition.

orthSVD

Arguments

:: Field t

=> Either Double Int

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

-> Matrix t

input matrix m

-> (Matrix t, Vector Double)

leftSV of m

-> Matrix t

orth

The range space a matrix from its precomputed SVD decomposition.

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.

type ℝ = Double

type ℂ = Complex Double

iC :: ℂ

imaginary unit

Auxiliary classes

class Storable a => Element a

Supported matrix elements.

This class provides optimized internal operations for selected element types. It provides unoptimised defaults for any Storable type, so you can create instances simply as:

instance Element Foo

Instances

Element Double

Element Float

Element (Complex Double)

Element (Complex Float)

class (Complexable c, Fractional e, Element e) => Container c e

Basic element-by-element functions for numeric containers

Minimal complete definition

size', scalar', conj', scale', scaleRecip, addConstant, add, sub, mul, divide, equal, arctan2', cmap', konst', build', atIndex', minIndex', maxIndex', minElement', maxElement', sumElements', prodElements', step', cond', find', assoc', accum'

Instances

Container Vector Double

Container Vector Float

Container Vector a => Container Matrix a

Container Vector (Complex Double)

Container Vector (Complex Float)

class (Num e, Element e) => Product e

Matrix product and related functions

Minimal complete definition

multiply, absSum, norm1, norm2, normInf

Instances

Product Double

Product Float

Product (Complex Double)

Product (Complex Float)

class (Container Vector t, Container Matrix t, Konst t Int Vector, Konst t (Int, Int) Matrix, Product t) => Numeric t

Instances

Numeric Double

Numeric Float

Numeric (Complex Double)

Numeric (Complex Float)

class LSDiv c

Minimal complete definition

linSolve

Instances

LSDiv Vector

LSDiv Matrix

class Complexable c

Structures that may contain complex numbers

Minimal complete definition

toComplex', fromComplex', comp', single', double'

Instances

Complexable Vector

Complexable Matrix

class (Element t, Element (Complex t), RealFloat t) => RealElement t

Supported real types

Instances

RealElement Double

RealElement Float

type family RealOf x

Instances

type RealOf Double = Double

type RealOf Float = Float

type RealOf (Complex Double) = Double

type RealOf (Complex Float) = Float

type family ComplexOf x

Instances

type ComplexOf Double = Complex Double

type ComplexOf Float = Complex Float

type ComplexOf (Complex Double) = Complex Double

type ComplexOf (Complex Float) = Complex Float

type family SingleOf x

Instances

type SingleOf Double = Float

type SingleOf Float = Float

type SingleOf (Complex a) = Complex (SingleOf a)

type family DoubleOf x

Instances

type DoubleOf Double = Double

type DoubleOf Float = Double

type DoubleOf (Complex a) = Complex (DoubleOf a)

type family IndexOf c

Instances

type IndexOf Vector = Int

type IndexOf Matrix = (Int, Int)

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

Generic linear algebra functions for double precision real and complex matrices.

(Single precision data can be converted using single and double).

Minimal complete definition

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

Instances

Field Double

Field (Complex Double)

class Transposable m mt | m -> mt, mt -> m where

Methods

tr :: m -> mt

(conjugate) transpose

Instances

Transposable GMatrix GMatrix

Container Vector t => Transposable (Matrix t) (Matrix t)

(KnownNat n, KnownNat m) => Transposable (M m n) (M n m)

(KnownNat n, KnownNat m) => Transposable (L m n) (L n m)

data CGState

Constructors

CGState

Fields

cgp :: V

conjugate gradient

cgr :: V

residual

cgr2 :: R

squared norm of residual

cgx :: V

current solution

cgdx :: R

normalized size of correction

class Testable t where

Minimal complete definition

checkT

Methods

checkT :: t -> (Bool, IO ())

ioCheckT :: t -> IO (Bool, IO ())

Instances

Testable GMatrix

(KnownNat n', KnownNat m') => Testable (L n' m')