summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-16 12:36:52 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-16 12:36:52 +0200
commita2d99e7d0e83fcedf3a856cdb927309e28a8eddd (patch)
treec6192d21e4b59e5527aac5d25d8b8baefc2b3052 /packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
parent5d15d765c0204854b587b03232b20336ddf91ced (diff)
container and algorithms moved to base
Diffstat (limited to 'packages/base/src/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Algorithms.hs744
1 files changed, 744 insertions, 0 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
new file mode 100644
index 0000000..6f40683
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -0,0 +1,744 @@
1{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
2{-# LANGUAGE CPP #-}
3{-# LANGUAGE MultiParamTypeClasses #-}
4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE TypeFamilies #-}
6
7-----------------------------------------------------------------------------
8{- |
9Module : Numeric.LinearAlgebra.Algorithms
10Copyright : (c) Alberto Ruiz 2006-9
11License : GPL-style
12
13Maintainer : Alberto Ruiz (aruiz at um dot es)
14Stability : provisional
15Portability : uses ffi
16
17High level generic interface to common matrix computations.
18
19Specific functions for particular base types can also be explicitly
20imported from "Numeric.LinearAlgebra.LAPACK".
21
22-}
23-----------------------------------------------------------------------------
24
25module Numeric.LinearAlgebra.Algorithms (
26-- * Supported types
27 Field(),
28-- * Linear Systems
29 linearSolve,
30 luSolve,
31 cholSolve,
32 linearSolveLS,
33 linearSolveSVD,
34 inv, pinv, pinvTol,
35 det, invlndet,
36 rank, rcond,
37-- * Matrix factorizations
38-- ** Singular value decomposition
39 svd,
40 fullSVD,
41 thinSVD,
42 compactSVD,
43 singularValues,
44 leftSV, rightSV,
45-- ** Eigensystems
46 eig, eigSH, eigSH',
47 eigenvalues, eigenvaluesSH, eigenvaluesSH',
48 geigSH',
49-- ** QR
50 qr, rq, qrRaw, qrgr,
51-- ** Cholesky
52 chol, cholSH, mbCholSH,
53-- ** Hessenberg
54 hess,
55-- ** Schur
56 schur,
57-- ** LU
58 lu, luPacked,
59-- * Matrix functions
60 expm,
61 sqrtm,
62 matFunc,
63-- * Nullspace
64 nullspacePrec,
65 nullVector,
66 nullspaceSVD,
67 orth,
68-- * Norms
69 Normed(..), NormType(..),
70 relativeError,
71-- * Misc
72 eps, peps, i,
73-- * Util
74 haussholder,
75 unpackQR, unpackHess,
76 ranksv
77) where
78
79
80import Data.Packed.Development hiding ((//))
81import Data.Packed
82import Numeric.LinearAlgebra.LAPACK as LAPACK
83import Data.List(foldl1')
84import Data.Array
85import Data.Packed.Numeric
86
87
88{- | Generic linear algebra functions for double precision real and complex matrices.
89
90(Single precision data can be converted using 'single' and 'double').
91
92-}
93class (Product t,
94 Convert t,
95 Container Vector t,
96 Container Matrix t,
97 Normed Matrix t,
98 Normed Vector t,
99 Floating t,
100 RealOf t ~ Double) => Field t where
101 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
102 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
103 sv' :: Matrix t -> Vector Double
104 luPacked' :: Matrix t -> (Matrix t, [Int])
105 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
106 linearSolve' :: Matrix t -> Matrix t -> Matrix t
107 cholSolve' :: Matrix t -> Matrix t -> Matrix t
108 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
109 linearSolveLS' :: Matrix t -> Matrix t -> Matrix t
110 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
111 eigSH'' :: Matrix t -> (Vector Double, Matrix t)
112 eigOnly :: Matrix t -> Vector (Complex Double)
113 eigOnlySH :: Matrix t -> Vector Double
114 cholSH' :: Matrix t -> Matrix t
115 mbCholSH' :: Matrix t -> Maybe (Matrix t)
116 qr' :: Matrix t -> (Matrix t, Vector t)
117 qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t
118 hess' :: Matrix t -> (Matrix t, Matrix t)
119 schur' :: Matrix t -> (Matrix t, Matrix t)
120
121
122instance Field Double where
123 svd' = svdRd
124 thinSVD' = thinSVDRd
125 sv' = svR
126 luPacked' = luR
127 luSolve' (l_u,perm) = lusR l_u perm
128 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
129 cholSolve' = cholSolveR
130 linearSolveLS' = linearSolveLSR
131 linearSolveSVD' = linearSolveSVDR Nothing
132 eig' = eigR
133 eigSH'' = eigS
134 eigOnly = eigOnlyR
135 eigOnlySH = eigOnlyS
136 cholSH' = cholS
137 mbCholSH' = mbCholS
138 qr' = qrR
139 qrgr' = qrgrR
140 hess' = unpackHess hessR
141 schur' = schurR
142
143instance Field (Complex Double) where
144#ifdef NOZGESDD
145 svd' = svdC
146 thinSVD' = thinSVDC
147#else
148 svd' = svdCd
149 thinSVD' = thinSVDCd
150#endif
151 sv' = svC
152 luPacked' = luC
153 luSolve' (l_u,perm) = lusC l_u perm
154 linearSolve' = linearSolveC
155 cholSolve' = cholSolveC
156 linearSolveLS' = linearSolveLSC
157 linearSolveSVD' = linearSolveSVDC Nothing
158 eig' = eigC
159 eigOnly = eigOnlyC
160 eigSH'' = eigH
161 eigOnlySH = eigOnlyH
162 cholSH' = cholH
163 mbCholSH' = mbCholH
164 qr' = qrC
165 qrgr' = qrgrC
166 hess' = unpackHess hessC
167 schur' = schurC
168
169--------------------------------------------------------------
170
171square m = rows m == cols m
172
173vertical m = rows m >= cols m
174
175exactHermitian m = m `equal` ctrans m
176
177--------------------------------------------------------------
178
179-- | Full singular value decomposition.
180svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
181svd = {-# SCC "svd" #-} svd'
182
183-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
184--
185-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
186thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
187thinSVD = {-# SCC "thinSVD" #-} thinSVD'
188
189-- | Singular values only.
190singularValues :: Field t => Matrix t -> Vector Double
191singularValues = {-# SCC "singularValues" #-} sv'
192
193-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
194--
195-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
196fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
197fullSVD m = (u,d,v) where
198 (u,s,v) = svd m
199 d = diagRect 0 s r c
200 r = rows m
201 c = cols m
202
203-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
204compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
205compactSVD m = (u', subVector 0 d s, v') where
206 (u,s,v) = thinSVD m
207 d = rankSVD (1*eps) m s `max` 1
208 u' = takeColumns d u
209 v' = takeColumns d v
210
211
212-- | Singular values and all right singular vectors.
213rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
214rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
215 | otherwise = let (_,s,v) = svd m in (s,v)
216
217-- | Singular values and all left singular vectors.
218leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
219leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
220 | otherwise = let (u,s,_) = thinSVD m in (u,s)
221
222
223--------------------------------------------------------------
224
225-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
226luPacked :: Field t => Matrix t -> (Matrix t, [Int])
227luPacked = {-# SCC "luPacked" #-} luPacked'
228
229-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
230luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
231luSolve = {-# SCC "luSolve" #-} luSolve'
232
233-- | 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'.
234-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
235linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
236linearSolve = {-# SCC "linearSolve" #-} linearSolve'
237
238-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
239cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
240cholSolve = {-# SCC "cholSolve" #-} cholSolve'
241
242-- | 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.
243linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
244linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'
245
246
247-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
248linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
249linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
250
251--------------------------------------------------------------
252
253-- | Eigenvalues and eigenvectors of a general square matrix.
254--
255-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
256eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
257eig = {-# SCC "eig" #-} eig'
258
259-- | Eigenvalues of a general square matrix.
260eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
261eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
262
263-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
264eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
265eigSH' = {-# SCC "eigSH'" #-} eigSH''
266
267-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
268eigenvaluesSH' :: Field t => Matrix t -> Vector Double
269eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH
270
271-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
272--
273-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
274eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
275eigSH m | exactHermitian m = eigSH' m
276 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
277
278-- | Eigenvalues of a complex hermitian or real symmetric matrix.
279eigenvaluesSH :: Field t => Matrix t -> Vector Double
280eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
281 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
282
283--------------------------------------------------------------
284
285-- | QR factorization.
286--
287-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
288qr :: Field t => Matrix t -> (Matrix t, Matrix t)
289qr = {-# SCC "qr" #-} unpackQR . qr'
290
291qrRaw m = qr' m
292
293{- | generate a matrix with k orthogonal columns from the output of qrRaw
294-}
295qrgr n (a,t)
296 | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)"
297 | otherwise = qrgr' n (a,t)
298
299-- | RQ factorization.
300--
301-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
302rq :: Field t => Matrix t -> (Matrix t, Matrix t)
303rq m = {-# SCC "rq" #-} (r,q) where
304 (q',r') = qr $ trans $ rev1 m
305 r = rev2 (trans r')
306 q = rev2 (trans q')
307 rev1 = flipud . fliprl
308 rev2 = fliprl . flipud
309
310-- | Hessenberg factorization.
311--
312-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
313-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
314hess :: Field t => Matrix t -> (Matrix t, Matrix t)
315hess = hess'
316
317-- | Schur factorization.
318--
319-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
320-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
321-- upper triangular in 2x2 blocks.
322--
323-- \"Anything that the Jordan decomposition can do, the Schur decomposition
324-- can do better!\" (Van Loan)
325schur :: Field t => Matrix t -> (Matrix t, Matrix t)
326schur = schur'
327
328
329-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
330mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
331mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'
332
333-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
334cholSH :: Field t => Matrix t -> Matrix t
335cholSH = {-# SCC "cholSH" #-} cholSH'
336
337-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
338--
339-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@.
340chol :: Field t => Matrix t -> Matrix t
341chol m | exactHermitian m = cholSH m
342 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
343
344
345-- | Joint computation of inverse and logarithm of determinant of a square matrix.
346invlndet :: Field t
347 => Matrix t
348 -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
349invlndet m | square m = (im,(ladm,sdm))
350 | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
351 where
352 lp@(lup,perm) = luPacked m
353 s = signlp (rows m) perm
354 dg = toList $ takeDiag $ lup
355 ladm = sum $ map (log.abs) dg
356 sdm = s* product (map signum dg)
357 im = luSolve lp (ident (rows m))
358
359
360-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
361det :: Field t => Matrix t -> t
362det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
363 | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
364 where (lup,perm) = luPacked m
365 s = signlp (rows m) perm
366
367-- | Explicit LU factorization of a general matrix.
368--
369-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
370-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
371lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
372lu = luFact . luPacked
373
374-- | Inverse of a square matrix. See also 'invlndet'.
375inv :: Field t => Matrix t -> Matrix t
376inv m | square m = m `linearSolve` ident (rows m)
377 | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"
378
379
380-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave).
381pinv :: Field t => Matrix t -> Matrix t
382pinv = pinvTol 1
383
384{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
385
386@
387m = (3><3) [ 1, 0, 0
388 , 0, 1, 0
389 , 0, 0, 1e-10] :: Matrix Double
390@
391
392>>> pinv m
3931. 0. 0.
3940. 1. 0.
3950. 0. 10000000000.
396
397>>> pinvTol 1E8 m
3981. 0. 0.
3990. 1. 0.
4000. 0. 1.
401
402-}
403
404pinvTol :: Field t => Double -> Matrix t -> Matrix t
405pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where
406 (u,s,v) = thinSVD m
407 sl@(g:_) = toList s
408 s' = real . fromList . map rec $ sl
409 rec x = if x <= g*tol then x else 1/x
410 tol = (fromIntegral (max r c) * g * t * eps)
411 r = rows m
412 c = cols m
413 d = dim s
414 u' = takeColumns d u
415 v' = takeColumns d v
416
417
418-- | Numeric rank of a matrix from the SVD decomposition.
419rankSVD :: Element t
420 => Double -- ^ numeric zero (e.g. 1*'eps')
421 -> Matrix t -- ^ input matrix m
422 -> Vector Double -- ^ 'sv' of m
423 -> Int -- ^ rank of m
424rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
425
426-- | Numeric rank of a matrix from its singular values.
427ranksv :: Double -- ^ numeric zero (e.g. 1*'eps')
428 -> Int -- ^ maximum dimension of the matrix
429 -> [Double] -- ^ singular values
430 -> Int -- ^ rank of m
431ranksv teps maxdim s = k where
432 g = maximum s
433 tol = fromIntegral maxdim * g * teps
434 s' = filter (>tol) s
435 k = if g > teps then length s' else 0
436
437-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
438eps :: Double
439eps = 2.22044604925031e-16
440
441
442-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1
443peps :: RealFloat x => x
444peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
445
446
447-- | The imaginary unit: @i = 0.0 :+ 1.0@
448i :: Complex Double
449i = 0:+1
450
451-----------------------------------------------------------------------
452
453-- | The nullspace of a matrix from its precomputed SVD decomposition.
454nullspaceSVD :: Field t
455 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
456 -- or Right \"theoretical\" matrix rank.
457 -> Matrix t -- ^ input matrix m
458 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
459 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
460nullspaceSVD hint a (s,v) = vs where
461 tol = case hint of
462 Left t -> t
463 _ -> eps
464 k = case hint of
465 Right t -> t
466 _ -> rankSVD tol a s
467 vs = drop k $ toRows $ ctrans v
468
469
470-- | The nullspace of a matrix. See also 'nullspaceSVD'.
471nullspacePrec :: Field t
472 => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
473 -> Matrix t -- ^ input matrix
474 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
475nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
476
477-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
478nullVector :: Field t => Matrix t -> Vector t
479nullVector = last . nullspacePrec 1
480
481orth :: Field t => Matrix t -> [Vector t]
482-- ^ Return an orthonormal basis of the range space of a matrix
483orth m = take r $ toColumns u
484 where
485 (u,s,_) = compactSVD m
486 r = ranksv eps (max (rows m) (cols m)) (toList s)
487
488------------------------------------------------------------------------
489
490-- many thanks, quickcheck!
491
492haussholder :: (Field a) => a -> Vector a -> Matrix a
493haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
494 where w = asColumn v
495
496
497zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
498 where xs = toList v
499
500zt 0 v = v
501zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k]
502
503
504unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
505unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r)
506 where cs = toColumns pq
507 m = rows pq
508 n = cols pq
509 mn = min m n
510 r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
511 vs = zipWith zh [1..mn] cs
512 hs = zipWith haussholder (toList tau) vs
513 q = foldl1' mXm hs
514
515unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
516unpackHess hf m
517 | rows m == 1 = ((1><1)[1],m)
518 | otherwise = (uH . hf) m
519
520uH (pq, tau) = (p,h)
521 where cs = toColumns pq
522 m = rows pq
523 n = cols pq
524 mn = min m n
525 h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs
526 vs = zipWith zh [2..mn] cs
527 hs = zipWith haussholder (toList tau) vs
528 p = foldl1' mXm hs
529
530--------------------------------------------------------------------------
531
532-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
533rcond :: Field t => Matrix t -> Double
534rcond m = last s / head s
535 where s = toList (singularValues m)
536
537-- | Number of linearly independent rows or columns.
538rank :: Field t => Matrix t -> Int
539rank m = rankSVD eps m (singularValues m)
540
541{-
542expm' m = case diagonalize (complex m) of
543 Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
544 Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
545 where exp = vectorMapC Exp
546-}
547
548diagonalize m = if rank v == n
549 then Just (l,v)
550 else Nothing
551 where n = rows m
552 (l,v) = if exactHermitian m
553 then let (l',v') = eigSH m in (real l', v')
554 else eig m
555
556-- | Generic matrix functions for diagonalizable matrices. For instance:
557--
558-- @logm = matFunc log@
559--
560matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
561matFunc f m = case diagonalize m of
562 Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v
563 Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"
564
565--------------------------------------------------------------
566
567golubeps :: Integer -> Integer -> Double
568golubeps p q = a * fromIntegral b / fromIntegral c where
569 a = 2^^(3-p-q)
570 b = fact p * fact q
571 c = fact (p+q) * fact (p+q+1)
572 fact n = product [1..n]
573
574epslist :: [(Int,Double)]
575epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
576
577geps delta = head [ k | (k,g) <- epslist, g<delta]
578
579
580{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
581 based on a scaled Pade approximation.
582-}
583expm :: Field t => Matrix t -> Matrix t
584expm = expGolub
585
586expGolub :: Field t => Matrix t -> Matrix t
587expGolub m = iterate msq f !! j
588 where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m
589 a = m */ fromIntegral ((2::Int)^j)
590 q = geps eps -- 7 steps
591 eye = ident (rows m)
592 work (k,c,x,n,d) = (k',c',x',n',d')
593 where k' = k+1
594 c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
595 x' = a <> x
596 n' = n |+| (c' .* x')
597 d' = d |+| (((-1)^k * c') .* x')
598 (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q
599 f = linearSolve df nf
600 msq x = x <> x
601
602 (<>) = multiply
603 v */ x = scale (recip x) v
604 (.*) = scale
605 (|+|) = add
606
607--------------------------------------------------------------
608
609{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
610It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@.
611
612@m = (2><2) [4,9
613 ,0,4] :: Matrix Double@
614
615>>> sqrtm m
616(2><2)
617 [ 2.0, 2.25
618 , 0.0, 2.0 ]
619
620-}
621sqrtm :: Field t => Matrix t -> Matrix t
622sqrtm = sqrtmInv
623
624sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
625 where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a
626 | otherwise = fixedPoint (b:rest)
627 fixedPoint _ = error "fixedpoint with impossible inputs"
628 f (y,z) = (0.5 .* (y |+| inv z),
629 0.5 .* (inv y |+| z))
630 (.*) = scale
631 (|+|) = add
632 (|-|) = sub
633
634------------------------------------------------------------------
635
636signlp r vals = foldl f 1 (zip [0..r-1] vals)
637 where f s (a,b) | a /= b = -s
638 | otherwise = s
639
640swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s)
641 | otherwise = (arr,s)
642
643fixPerm r vals = (fromColumns $ elems res, sign)
644 where v = [0..r-1]
645 s = toColumns (ident r)
646 (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
647
648triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
649 where el p q = if q-p>=h then v else 1 - v
650
651luFact (l_u,perm) | r <= c = (l ,u ,p, s)
652 | otherwise = (l',u',p, s)
653 where
654 r = rows l_u
655 c = cols l_u
656 tu = triang r c 0 1
657 tl = triang r c 0 0
658 l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r
659 u = l_u |*| tu
660 (p,s) = fixPerm r perm
661 l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c
662 u' = takeRows c (l_u |*| tu)
663 (|+|) = add
664 (|*|) = mul
665
666---------------------------------------------------------------------------
667
668data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
669
670class (RealFloat (RealOf t)) => Normed c t where
671 pnorm :: NormType -> c t -> RealOf t
672
673instance Normed Vector Double where
674 pnorm PNorm1 = norm1
675 pnorm PNorm2 = norm2
676 pnorm Infinity = normInf
677 pnorm Frobenius = norm2
678
679instance Normed Vector (Complex Double) where
680 pnorm PNorm1 = norm1
681 pnorm PNorm2 = norm2
682 pnorm Infinity = normInf
683 pnorm Frobenius = pnorm PNorm2
684
685instance Normed Vector Float where
686 pnorm PNorm1 = norm1
687 pnorm PNorm2 = norm2
688 pnorm Infinity = normInf
689 pnorm Frobenius = pnorm PNorm2
690
691instance Normed Vector (Complex Float) where
692 pnorm PNorm1 = norm1
693 pnorm PNorm2 = norm2
694 pnorm Infinity = normInf
695 pnorm Frobenius = pnorm PNorm2
696
697
698instance Normed Matrix Double where
699 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
700 pnorm PNorm2 = (@>0) . singularValues
701 pnorm Infinity = pnorm PNorm1 . trans
702 pnorm Frobenius = pnorm PNorm2 . flatten
703
704instance Normed Matrix (Complex Double) where
705 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
706 pnorm PNorm2 = (@>0) . singularValues
707 pnorm Infinity = pnorm PNorm1 . trans
708 pnorm Frobenius = pnorm PNorm2 . flatten
709
710instance Normed Matrix Float where
711 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
712 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
713 pnorm Infinity = pnorm PNorm1 . trans
714 pnorm Frobenius = pnorm PNorm2 . flatten
715
716instance Normed Matrix (Complex Float) where
717 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
718 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
719 pnorm Infinity = pnorm PNorm1 . trans
720 pnorm Frobenius = pnorm PNorm2 . flatten
721
722-- | Approximate number of common digits in the maximum element.
723relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
724relativeError x y = dig (norm (x `sub` y) / norm x)
725 where norm = pnorm Infinity
726 dig r = round $ -logBase 10 (realToFrac r :: Double)
727
728----------------------------------------------------------------------
729
730-- | Generalized symmetric positive definite eigensystem Av = lBv,
731-- for A and B symmetric, B positive definite (conditions not checked).
732geigSH' :: Field t
733 => Matrix t -- ^ A
734 -> Matrix t -- ^ B
735 -> (Vector Double, Matrix t)
736geigSH' a b = (l,v')
737 where
738 u = cholSH b
739 iu = inv u
740 c = ctrans iu <> a <> iu
741 (l,v) = eigSH' c
742 v' = iu <> v
743 (<>) = mXm
744