summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs46
1 files changed, 23 insertions, 23 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 794ef69..b7e208a 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -51,7 +51,7 @@ module Numeric.LinearAlgebra.Algorithms (
51-- * Util 51-- * Util
52 haussholder, 52 haussholder,
53 unpackQR, unpackHess, 53 unpackQR, unpackHess,
54 GenMat(linearSolveSVD,lu,eigSH',cholSH) 54 Field(linearSolveSVD,lu,eigSH',cholSH)
55) where 55) where
56 56
57 57
@@ -65,7 +65,7 @@ import Numeric.LinearAlgebra.Linear
65import Data.List(foldl1') 65import Data.List(foldl1')
66 66
67-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 67-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
68class (Normed (Matrix t), Linear Matrix t) => GenMat t where 68class (Normed (Matrix t), Linear Matrix t) => Field t where
69 -- | Singular value decomposition using lapack's dgesvd or zgesvd. 69 -- | Singular value decomposition using lapack's dgesvd or zgesvd.
70 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) 70 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t)
71 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) 71 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t)
@@ -103,7 +103,7 @@ class (Normed (Matrix t), Linear Matrix t) => GenMat t where
103 ctrans :: Matrix t -> Matrix t 103 ctrans :: Matrix t -> Matrix t
104 104
105 105
106instance GenMat Double where 106instance Field Double where
107 svd = svdR 107 svd = svdR
108 lu = GSL.luR 108 lu = GSL.luR
109 linearSolve = linearSolveR 109 linearSolve = linearSolveR
@@ -116,7 +116,7 @@ instance GenMat Double where
116 hess = unpackHess hessR 116 hess = unpackHess hessR
117 schur = schurR 117 schur = schurR
118 118
119instance GenMat (Complex Double) where 119instance Field (Complex Double) where
120 svd = svdC 120 svd = svdC
121 lu = GSL.luC 121 lu = GSL.luC
122 linearSolve = linearSolveC 122 linearSolve = linearSolveC
@@ -132,37 +132,37 @@ instance GenMat (Complex Double) where
132-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev. 132-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
133-- 133--
134-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ 134-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
135eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) 135eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
136eigSH m | m `equal` ctrans m = eigSH' m 136eigSH m | m `equal` ctrans m = eigSH' m
137 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" 137 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
138 138
139-- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. 139-- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
140-- 140--
141-- If @c = chol m@ then @m == c \<> ctrans c@. 141-- If @c = chol m@ then @m == c \<> ctrans c@.
142chol :: GenMat t => Matrix t -> Matrix t 142chol :: Field t => Matrix t -> Matrix t
143chol m | m `equal` ctrans m = cholSH m 143chol m | m `equal` ctrans m = cholSH m
144 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" 144 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
145 145
146square m = rows m == cols m 146square m = rows m == cols m
147 147
148det :: GenMat t => Matrix t -> t 148det :: Field t => Matrix t -> t
149det m | square m = s * (product $ toList $ takeDiag $ u) 149det m | square m = s * (product $ toList $ takeDiag $ u)
150 | otherwise = error "det of nonsquare matrix" 150 | otherwise = error "det of nonsquare matrix"
151 where (_,u,_,s) = lu m 151 where (_,u,_,s) = lu m
152 152
153-- | Inverse of a square matrix using lapacks' dgesv and zgesv. 153-- | Inverse of a square matrix using lapacks' dgesv and zgesv.
154inv :: GenMat t => Matrix t -> Matrix t 154inv :: Field t => Matrix t -> Matrix t
155inv m | square m = m `linearSolve` ident (rows m) 155inv m | square m = m `linearSolve` ident (rows m)
156 | otherwise = error "inv of nonsquare matrix" 156 | otherwise = error "inv of nonsquare matrix"
157 157
158-- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. 158-- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss.
159pinv :: GenMat t => Matrix t -> Matrix t 159pinv :: Field t => Matrix t -> Matrix t
160pinv m = linearSolveSVD m (ident (rows m)) 160pinv m = linearSolveSVD m (ident (rows m))
161 161
162-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. 162-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
163-- 163--
164-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. 164-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@.
165full :: Field t 165full :: Element t
166 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) 166 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
167full svd m = (u, d ,v) where 167full svd m = (u, d ,v) where
168 (u,s,v) = svd m 168 (u,s,v) = svd m
@@ -173,7 +173,7 @@ full svd m = (u, d ,v) where
173-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. 173-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
174-- 174--
175-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. 175-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@.
176economy :: Field t 176economy :: Element t
177 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 177 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
178economy svd m = (u', subVector 0 d s, v') where 178economy svd m = (u', subVector 0 d s, v') where
179 (u,s,v) = svd m 179 (u,s,v) = svd m
@@ -198,15 +198,15 @@ i = 0:+1
198 198
199 199
200-- matrix product 200-- matrix product
201mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t 201mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t
202mXm = multiply 202mXm = multiply
203 203
204-- matrix - vector product 204-- matrix - vector product
205mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t 205mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t
206mXv m v = flatten $ m `mXm` (asColumn v) 206mXv m v = flatten $ m `mXm` (asColumn v)
207 207
208-- vector - matrix product 208-- vector - matrix product
209vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t 209vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t
210vXm v m = flatten $ (asRow v) `mXm` m 210vXm v m = flatten $ (asRow v) `mXm` m
211 211
212 212
@@ -264,7 +264,7 @@ instance Normed (Matrix (Complex Double)) where
264----------------------------------------------------------------------- 264-----------------------------------------------------------------------
265 265
266-- | The nullspace of a matrix from its SVD decomposition. 266-- | The nullspace of a matrix from its SVD decomposition.
267nullspacePrec :: GenMat t 267nullspacePrec :: Field t
268 => Double -- ^ relative tolerance in 'eps' units 268 => Double -- ^ relative tolerance in 'eps' units
269 -> Matrix t -- ^ input matrix 269 -> Matrix t -- ^ input matrix
270 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 270 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
@@ -276,7 +276,7 @@ nullspacePrec t m = ns where
276 ns = drop rank $ toRows $ ctrans v 276 ns = drop rank $ toRows $ ctrans v
277 277
278-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). 278-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
279nullVector :: GenMat t => Matrix t -> Vector t 279nullVector :: Field t => Matrix t -> Vector t
280nullVector = last . nullspacePrec 1 280nullVector = last . nullspacePrec 1
281 281
282------------------------------------------------------------------------ 282------------------------------------------------------------------------
@@ -316,7 +316,7 @@ pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
316 316
317-- many thanks, quickcheck! 317-- many thanks, quickcheck!
318 318
319haussholder :: (GenMat a) => a -> Vector a -> Matrix a 319haussholder :: (Field a) => a -> Vector a -> Matrix a
320haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) 320haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
321 where w = asColumn v 321 where w = asColumn v
322 322
@@ -328,7 +328,7 @@ zt 0 v = v
328zt k v = join [subVector 0 (dim v - k) v, constant 0 k] 328zt k v = join [subVector 0 (dim v - k) v, constant 0 k]
329 329
330 330
331unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) 331unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
332unpackQR (pq, tau) = (q,r) 332unpackQR (pq, tau) = (q,r)
333 where cs = toColumns pq 333 where cs = toColumns pq
334 m = rows pq 334 m = rows pq
@@ -339,7 +339,7 @@ unpackQR (pq, tau) = (q,r)
339 hs = zipWith haussholder (toList tau) vs 339 hs = zipWith haussholder (toList tau) vs
340 q = foldl1' mXm hs 340 q = foldl1' mXm hs
341 341
342unpackHess :: (GenMat t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) 342unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
343unpackHess hf m 343unpackHess hf m
344 | rows m == 1 = ((1><1)[1],m) 344 | rows m == 1 = ((1><1)[1],m)
345 | otherwise = (uH . hf) m 345 | otherwise = (uH . hf) m
@@ -357,13 +357,13 @@ uH (pq, tau) = (p,h)
357-------------------------------------------------------------------------- 357--------------------------------------------------------------------------
358 358
359-- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. 359-- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD.
360rcond :: GenMat t => Matrix t -> Double 360rcond :: Field t => Matrix t -> Double
361rcond m = last s / head s 361rcond m = last s / head s
362 where (_,s',_) = svd m 362 where (_,s',_) = svd m
363 s = toList s' 363 s = toList s'
364 364
365-- | Number of linearly independent rows or columns. 365-- | Number of linearly independent rows or columns.
366rank :: GenMat t => Matrix t -> Int 366rank :: Field t => Matrix t -> Int
367rank m | pnorm PNorm1 m < eps = 0 367rank m | pnorm PNorm1 m < eps = 0
368 | otherwise = dim s where (_,s,_) = economy svd m 368 | otherwise = dim s where (_,s,_) = economy svd m
369 369
@@ -381,7 +381,7 @@ diagonalize m = if rank v == n
381 else eig m 381 else eig m
382 382
383-- | Generic matrix functions for diagonalizable matrices. 383-- | Generic matrix functions for diagonalizable matrices.
384matFunc :: GenMat t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) 384matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double)
385matFunc f m = case diagonalize (complex m) of 385matFunc f m = case diagonalize (complex m) of
386 Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v 386 Just (l,v) -> v `mXm` diag (liftVector f l) `mXm` inv v
387 Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix" 387 Nothing -> error "Sorry, matFunc requieres a diagonalizable matrix"
@@ -420,5 +420,5 @@ expGolub m = iterate msq f !! j
420{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, 420{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
421 based on a scaled Pade approximation. 421 based on a scaled Pade approximation.
422-} 422-}
423expm :: GenMat t => Matrix t -> Matrix t 423expm :: Field t => Matrix t -> Matrix t
424expm = expGolub 424expm = expGolub