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