diff options
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 46 |
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 | |||
65 | import Data.List(foldl1') | 65 | import 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. |
68 | class (Normed (Matrix t), Linear Matrix t) => GenMat t where | 68 | class (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 | ||
106 | instance GenMat Double where | 106 | instance 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 | ||
119 | instance GenMat (Complex Double) where | 119 | instance 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@ |
135 | eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t) | 135 | eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) |
136 | eigSH m | m `equal` ctrans m = eigSH' m | 136 | eigSH 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@. |
142 | chol :: GenMat t => Matrix t -> Matrix t | 142 | chol :: Field t => Matrix t -> Matrix t |
143 | chol m | m `equal` ctrans m = cholSH m | 143 | chol 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 | ||
146 | square m = rows m == cols m | 146 | square m = rows m == cols m |
147 | 147 | ||
148 | det :: GenMat t => Matrix t -> t | 148 | det :: Field t => Matrix t -> t |
149 | det m | square m = s * (product $ toList $ takeDiag $ u) | 149 | det 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. |
154 | inv :: GenMat t => Matrix t -> Matrix t | 154 | inv :: Field t => Matrix t -> Matrix t |
155 | inv m | square m = m `linearSolve` ident (rows m) | 155 | inv 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. |
159 | pinv :: GenMat t => Matrix t -> Matrix t | 159 | pinv :: Field t => Matrix t -> Matrix t |
160 | pinv m = linearSolveSVD m (ident (rows m)) | 160 | pinv 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@. |
165 | full :: Field t | 165 | full :: 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) |
167 | full svd m = (u, d ,v) where | 167 | full 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@. |
176 | economy :: Field t | 176 | economy :: 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) |
178 | economy svd m = (u', subVector 0 d s, v') where | 178 | economy 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 |
201 | mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t | 201 | mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t |
202 | mXm = multiply | 202 | mXm = multiply |
203 | 203 | ||
204 | -- matrix - vector product | 204 | -- matrix - vector product |
205 | mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t | 205 | mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t |
206 | mXv m v = flatten $ m `mXm` (asColumn v) | 206 | mXv m v = flatten $ m `mXm` (asColumn v) |
207 | 207 | ||
208 | -- vector - matrix product | 208 | -- vector - matrix product |
209 | vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t | 209 | vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t |
210 | vXm v m = flatten $ (asRow v) `mXm` m | 210 | vXm 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. |
267 | nullspacePrec :: GenMat t | 267 | nullspacePrec :: 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@). |
279 | nullVector :: GenMat t => Matrix t -> Vector t | 279 | nullVector :: Field t => Matrix t -> Vector t |
280 | nullVector = last . nullspacePrec 1 | 280 | nullVector = 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 | ||
319 | haussholder :: (GenMat a) => a -> Vector a -> Matrix a | 319 | haussholder :: (Field a) => a -> Vector a -> Matrix a |
320 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) | 320 | haussholder 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 | |||
328 | zt k v = join [subVector 0 (dim v - k) v, constant 0 k] | 328 | zt k v = join [subVector 0 (dim v - k) v, constant 0 k] |
329 | 329 | ||
330 | 330 | ||
331 | unpackQR :: (GenMat t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) | 331 | unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) |
332 | unpackQR (pq, tau) = (q,r) | 332 | unpackQR (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 | ||
342 | unpackHess :: (GenMat t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) | 342 | unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) |
343 | unpackHess hf m | 343 | unpackHess 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. |
360 | rcond :: GenMat t => Matrix t -> Double | 360 | rcond :: Field t => Matrix t -> Double |
361 | rcond m = last s / head s | 361 | rcond 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. |
366 | rank :: GenMat t => Matrix t -> Int | 366 | rank :: Field t => Matrix t -> Int |
367 | rank m | pnorm PNorm1 m < eps = 0 | 367 | rank 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. |
384 | matFunc :: GenMat t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) | 384 | matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) |
385 | matFunc f m = case diagonalize (complex m) of | 385 | matFunc 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 | -} |
423 | expm :: GenMat t => Matrix t -> Matrix t | 423 | expm :: Field t => Matrix t -> Matrix t |
424 | expm = expGolub | 424 | expm = expGolub |