summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2015-06-05 16:43:45 +0200
committerAlberto Ruiz <aruiz@um.es>2015-06-05 16:43:45 +0200
commit11d7c37dc8b314338bc6382d80e74aaec2bb5620 (patch)
tree33d9d74c064ad2d9bbc4da0643404b85e41032c6 /packages/base/src/Numeric
parent13856329ed09cd9f70f03363895545a6ca83374c (diff)
move algorithms
Diffstat (limited to 'packages/base/src/Numeric')
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Algorithms.hs972
1 files changed, 0 insertions, 972 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
deleted file mode 100644
index 0e085c1..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ /dev/null
@@ -1,972 +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-14
11License : BSD3
12Maintainer : Alberto Ruiz
13Stability : provisional
14
15High level generic interface to common matrix computations.
16
17Specific functions for particular base types can also be explicitly
18imported from "Numeric.LinearAlgebra.LAPACK".
19
20-}
21{-# OPTIONS_HADDOCK hide #-}
22-----------------------------------------------------------------------------
23
24module Numeric.LinearAlgebra.Algorithms (
25-- * Supported types
26 Field(),
27-- * Linear Systems
28 linearSolve,
29 mbLinearSolve,
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 orthSVD,
68 orth,
69-- * Norms
70 Normed(..), NormType(..),
71 relativeError', relativeError,
72-- * Misc
73 eps, peps, i,
74-- * Util
75 haussholder,
76 unpackQR, unpackHess,
77 ranksv
78) where
79
80
81import Data.Packed
82import Numeric.LinearAlgebra.LAPACK as LAPACK
83import Data.List(foldl1')
84import Data.Array
85import Data.Packed.Internal.Numeric
86import Data.Packed.Internal(shSize)
87
88
89{- | Generic linear algebra functions for double precision real and complex matrices.
90
91(Single precision data can be converted using 'single' and 'double').
92
93-}
94class (Product t,
95 Convert t,
96 Container Vector t,
97 Container Matrix t,
98 Normed Matrix t,
99 Normed Vector t,
100 Floating t,
101 RealOf t ~ Double) => Field t where
102 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
103 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
104 sv' :: Matrix t -> Vector Double
105 luPacked' :: Matrix t -> (Matrix t, [Int])
106 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
107 mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (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 mbLinearSolve' = mbLinearSolveR
132 cholSolve' = cholSolveR
133 linearSolveLS' = linearSolveLSR
134 linearSolveSVD' = linearSolveSVDR Nothing
135 eig' = eigR
136 eigSH'' = eigS
137 eigOnly = eigOnlyR
138 eigOnlySH = eigOnlyS
139 cholSH' = cholS
140 mbCholSH' = mbCholS
141 qr' = qrR
142 qrgr' = qrgrR
143 hess' = unpackHess hessR
144 schur' = schurR
145
146instance Field (Complex Double) where
147#ifdef NOZGESDD
148 svd' = svdC
149 thinSVD' = thinSVDC
150#else
151 svd' = svdCd
152 thinSVD' = thinSVDCd
153#endif
154 sv' = svC
155 luPacked' = luC
156 luSolve' (l_u,perm) = lusC l_u perm
157 linearSolve' = linearSolveC
158 mbLinearSolve' = mbLinearSolveC
159 cholSolve' = cholSolveC
160 linearSolveLS' = linearSolveLSC
161 linearSolveSVD' = linearSolveSVDC Nothing
162 eig' = eigC
163 eigOnly = eigOnlyC
164 eigSH'' = eigH
165 eigOnlySH = eigOnlyH
166 cholSH' = cholH
167 mbCholSH' = mbCholH
168 qr' = qrC
169 qrgr' = qrgrC
170 hess' = unpackHess hessC
171 schur' = schurC
172
173--------------------------------------------------------------
174
175square m = rows m == cols m
176
177vertical m = rows m >= cols m
178
179exactHermitian m = m `equal` ctrans m
180
181--------------------------------------------------------------
182
183{- | Full singular value decomposition.
184
185@
186a = (5><3)
187 [ 1.0, 2.0, 3.0
188 , 4.0, 5.0, 6.0
189 , 7.0, 8.0, 9.0
190 , 10.0, 11.0, 12.0
191 , 13.0, 14.0, 15.0 ] :: Matrix Double
192@
193
194>>> let (u,s,v) = svd a
195
196>>> disp 3 u
1975x5
198-0.101 0.768 0.614 0.028 -0.149
199-0.249 0.488 -0.503 0.172 0.646
200-0.396 0.208 -0.405 -0.660 -0.449
201-0.543 -0.072 -0.140 0.693 -0.447
202-0.690 -0.352 0.433 -0.233 0.398
203
204>>> s
205fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
206
207>>> disp 3 v
2083x3
209-0.519 -0.751 0.408
210-0.576 -0.046 -0.816
211-0.632 0.659 0.408
212
213>>> let d = diagRect 0 s 5 3
214>>> disp 3 d
2155x3
21635.183 0.000 0.000
217 0.000 1.477 0.000
218 0.000 0.000 0.000
219 0.000 0.000 0.000
220
221>>> disp 3 $ u <> d <> tr v
2225x3
223 1.000 2.000 3.000
224 4.000 5.000 6.000
225 7.000 8.000 9.000
22610.000 11.000 12.000
22713.000 14.000 15.000
228
229-}
230svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
231svd = {-# SCC "svd" #-} g . svd'
232 where
233 g (u,s,v) = (u,s,tr v)
234
235{- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
236
237If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> tr v@.
238
239@
240a = (5><3)
241 [ 1.0, 2.0, 3.0
242 , 4.0, 5.0, 6.0
243 , 7.0, 8.0, 9.0
244 , 10.0, 11.0, 12.0
245 , 13.0, 14.0, 15.0 ] :: Matrix Double
246@
247
248>>> let (u,s,v) = thinSVD a
249
250>>> disp 3 u
2515x3
252-0.101 0.768 0.614
253-0.249 0.488 -0.503
254-0.396 0.208 -0.405
255-0.543 -0.072 -0.140
256-0.690 -0.352 0.433
257
258>>> s
259fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
260
261>>> disp 3 v
2623x3
263-0.519 -0.751 0.408
264-0.576 -0.046 -0.816
265-0.632 0.659 0.408
266
267>>> disp 3 $ u <> diag s <> tr v
2685x3
269 1.000 2.000 3.000
270 4.000 5.000 6.000
271 7.000 8.000 9.000
27210.000 11.000 12.000
27313.000 14.000 15.000
274
275-}
276thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
277thinSVD = {-# SCC "thinSVD" #-} g . thinSVD'
278 where
279 g (u,s,v) = (u,s,tr v)
280
281
282-- | Singular values only.
283singularValues :: Field t => Matrix t -> Vector Double
284singularValues = {-# SCC "singularValues" #-} sv'
285
286-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
287--
288-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> tr v@.
289fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
290fullSVD m = (u,d,v) where
291 (u,s,v) = svd m
292 d = diagRect 0 s r c
293 r = rows m
294 c = cols m
295
296{- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
297
298@
299a = (5><3)
300 [ 1.0, 2.0, 3.0
301 , 4.0, 5.0, 6.0
302 , 7.0, 8.0, 9.0
303 , 10.0, 11.0, 12.0
304 , 13.0, 14.0, 15.0 ] :: Matrix Double
305@
306
307>>> let (u,s,v) = compactSVD a
308
309>>> disp 3 u
3105x2
311-0.101 0.768
312-0.249 0.488
313-0.396 0.208
314-0.543 -0.072
315-0.690 -0.352
316
317>>> s
318fromList [35.18264833189422,1.4769076999800903]
319
320>>> disp 3 u
3215x2
322-0.101 0.768
323-0.249 0.488
324-0.396 0.208
325-0.543 -0.072
326-0.690 -0.352
327
328>>> disp 3 $ u <> diag s <> tr v
3295x3
330 1.000 2.000 3.000
331 4.000 5.000 6.000
332 7.000 8.000 9.000
33310.000 11.000 12.000
33413.000 14.000 15.000
335
336-}
337compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
338compactSVD m = (u', subVector 0 d s, v') where
339 (u,s,v) = thinSVD m
340 d = rankSVD (1*eps) m s `max` 1
341 u' = takeColumns d u
342 v' = takeColumns d v
343
344
345-- | Singular values and all right singular vectors (as columns).
346rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
347rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
348 | otherwise = let (_,s,v) = svd m in (s,v)
349
350-- | Singular values and all left singular vectors (as columns).
351leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
352leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
353 | otherwise = let (u,s,_) = thinSVD m in (u,s)
354
355
356--------------------------------------------------------------
357
358-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
359luPacked :: Field t => Matrix t -> (Matrix t, [Int])
360luPacked = {-# SCC "luPacked" #-} luPacked'
361
362-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
363luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
364luSolve = {-# SCC "luSolve" #-} luSolve'
365
366-- | 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'.
367-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
368linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
369linearSolve = {-# SCC "linearSolve" #-} linearSolve'
370
371-- | 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'.
372mbLinearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
373mbLinearSolve = {-# SCC "linearSolve" #-} mbLinearSolve'
374
375-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
376cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
377cholSolve = {-# SCC "cholSolve" #-} cholSolve'
378
379-- | 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.
380linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
381linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'
382
383
384-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
385linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
386linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
387
388--------------------------------------------------------------
389
390{- | Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.
391
392If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
393
394@
395a = (3><3)
396 [ 3, 0, -2
397 , 4, 5, -1
398 , 3, 1, 0 ] :: Matrix Double
399@
400
401>>> let (l, v) = eig a
402
403>>> putStr . dispcf 3 . asRow $ l
4041x3
4051.925+1.523i 1.925-1.523i 4.151
406
407>>> putStr . dispcf 3 $ v
4083x3
409-0.455+0.365i -0.455-0.365i 0.181
410 0.603 0.603 -0.978
411 0.033+0.543i 0.033-0.543i -0.104
412
413>>> putStr . dispcf 3 $ complex a <> v
4143x3
415-1.432+0.010i -1.432-0.010i 0.753
416 1.160+0.918i 1.160-0.918i -4.059
417-0.763+1.096i -0.763-1.096i -0.433
418
419>>> putStr . dispcf 3 $ v <> diag l
4203x3
421-1.432+0.010i -1.432-0.010i 0.753
422 1.160+0.918i 1.160-0.918i -4.059
423-0.763+1.096i -0.763-1.096i -0.433
424
425-}
426eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
427eig = {-# SCC "eig" #-} eig'
428
429-- | Eigenvalues (not ordered) of a general square matrix.
430eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
431eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
432
433-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
434eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
435eigSH' = {-# SCC "eigSH'" #-} eigSH''
436
437-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
438eigenvaluesSH' :: Field t => Matrix t -> Vector Double
439eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH
440
441{- | Eigenvalues and eigenvectors (as columns) of a complex hermitian or real symmetric matrix, in descending order.
442
443If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@
444
445@
446a = (3><3)
447 [ 1.0, 2.0, 3.0
448 , 2.0, 4.0, 5.0
449 , 3.0, 5.0, 6.0 ]
450@
451
452>>> let (l, v) = eigSH a
453
454>>> l
455fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
456
457>>> disp 3 $ v <> diag l <> tr v
4583x3
4591.000 2.000 3.000
4602.000 4.000 5.000
4613.000 5.000 6.000
462
463-}
464eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
465eigSH m | exactHermitian m = eigSH' m
466 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
467
468-- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
469eigenvaluesSH :: Field t => Matrix t -> Vector Double
470eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
471 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
472
473--------------------------------------------------------------
474
475-- | QR factorization.
476--
477-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
478qr :: Field t => Matrix t -> (Matrix t, Matrix t)
479qr = {-# SCC "qr" #-} unpackQR . qr'
480
481qrRaw m = qr' m
482
483{- | generate a matrix with k orthogonal columns from the output of qrRaw
484-}
485qrgr n (a,t)
486 | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)"
487 | otherwise = qrgr' n (a,t)
488
489-- | RQ factorization.
490--
491-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
492rq :: Field t => Matrix t -> (Matrix t, Matrix t)
493rq m = {-# SCC "rq" #-} (r,q) where
494 (q',r') = qr $ trans $ rev1 m
495 r = rev2 (trans r')
496 q = rev2 (trans q')
497 rev1 = flipud . fliprl
498 rev2 = fliprl . flipud
499
500-- | Hessenberg factorization.
501--
502-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
503-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
504hess :: Field t => Matrix t -> (Matrix t, Matrix t)
505hess = hess'
506
507-- | Schur factorization.
508--
509-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
510-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
511-- upper triangular in 2x2 blocks.
512--
513-- \"Anything that the Jordan decomposition can do, the Schur decomposition
514-- can do better!\" (Van Loan)
515schur :: Field t => Matrix t -> (Matrix t, Matrix t)
516schur = schur'
517
518
519-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
520mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
521mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'
522
523-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
524cholSH :: Field t => Matrix t -> Matrix t
525cholSH = {-# SCC "cholSH" #-} cholSH'
526
527-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
528--
529-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@.
530chol :: Field t => Matrix t -> Matrix t
531chol m | exactHermitian m = cholSH m
532 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
533
534
535-- | Joint computation of inverse and logarithm of determinant of a square matrix.
536invlndet :: Field t
537 => Matrix t
538 -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
539invlndet m | square m = (im,(ladm,sdm))
540 | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
541 where
542 lp@(lup,perm) = luPacked m
543 s = signlp (rows m) perm
544 dg = toList $ takeDiag $ lup
545 ladm = sum $ map (log.abs) dg
546 sdm = s* product (map signum dg)
547 im = luSolve lp (ident (rows m))
548
549
550-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
551det :: Field t => Matrix t -> t
552det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
553 | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
554 where (lup,perm) = luPacked m
555 s = signlp (rows m) perm
556
557-- | Explicit LU factorization of a general matrix.
558--
559-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
560-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
561lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
562lu = luFact . luPacked
563
564-- | Inverse of a square matrix. See also 'invlndet'.
565inv :: Field t => Matrix t -> Matrix t
566inv m | square m = m `linearSolve` ident (rows m)
567 | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"
568
569
570-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave).
571pinv :: Field t => Matrix t -> Matrix t
572pinv = pinvTol 1
573
574{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
575
576@
577m = (3><3) [ 1, 0, 0
578 , 0, 1, 0
579 , 0, 0, 1e-10] :: Matrix Double
580@
581
582>>> pinv m
5831. 0. 0.
5840. 1. 0.
5850. 0. 10000000000.
586
587>>> pinvTol 1E8 m
5881. 0. 0.
5890. 1. 0.
5900. 0. 1.
591
592-}
593
594pinvTol :: Field t => Double -> Matrix t -> Matrix t
595pinvTol t m = v' `mXm` diag s' `mXm` ctrans u' where
596 (u,s,v) = thinSVD m
597 sl@(g:_) = toList s
598 s' = real . fromList . map rec $ sl
599 rec x = if x <= g*tol then x else 1/x
600 tol = (fromIntegral (max r c) * g * t * eps)
601 r = rows m
602 c = cols m
603 d = dim s
604 u' = takeColumns d u
605 v' = takeColumns d v
606
607
608-- | Numeric rank of a matrix from the SVD decomposition.
609rankSVD :: Element t
610 => Double -- ^ numeric zero (e.g. 1*'eps')
611 -> Matrix t -- ^ input matrix m
612 -> Vector Double -- ^ 'sv' of m
613 -> Int -- ^ rank of m
614rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
615
616-- | Numeric rank of a matrix from its singular values.
617ranksv :: Double -- ^ numeric zero (e.g. 1*'eps')
618 -> Int -- ^ maximum dimension of the matrix
619 -> [Double] -- ^ singular values
620 -> Int -- ^ rank of m
621ranksv teps maxdim s = k where
622 g = maximum s
623 tol = fromIntegral maxdim * g * teps
624 s' = filter (>tol) s
625 k = if g > teps then length s' else 0
626
627-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
628eps :: Double
629eps = 2.22044604925031e-16
630
631
632-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1
633peps :: RealFloat x => x
634peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
635
636
637-- | The imaginary unit: @i = 0.0 :+ 1.0@
638i :: Complex Double
639i = 0:+1
640
641-----------------------------------------------------------------------
642
643-- | The nullspace of a matrix from its precomputed SVD decomposition.
644nullspaceSVD :: Field t
645 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
646 -- or Right \"theoretical\" matrix rank.
647 -> Matrix t -- ^ input matrix m
648 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
649 -> Matrix t -- ^ nullspace
650nullspaceSVD hint a (s,v) = vs where
651 tol = case hint of
652 Left t -> t
653 _ -> eps
654 k = case hint of
655 Right t -> t
656 _ -> rankSVD tol a s
657 vs = dropColumns k v
658
659
660-- | The nullspace of a matrix. See also 'nullspaceSVD'.
661nullspacePrec :: Field t
662 => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
663 -> Matrix t -- ^ input matrix
664 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
665nullspacePrec t m = toColumns $ nullspaceSVD (Left (t*eps)) m (rightSV m)
666
667-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
668nullVector :: Field t => Matrix t -> Vector t
669nullVector = last . nullspacePrec 1
670
671-- | The range space a matrix from its precomputed SVD decomposition.
672orthSVD :: Field t
673 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
674 -- or Right \"theoretical\" matrix rank.
675 -> Matrix t -- ^ input matrix m
676 -> (Matrix t, Vector Double) -- ^ 'leftSV' of m
677 -> Matrix t -- ^ orth
678orthSVD hint a (v,s) = vs where
679 tol = case hint of
680 Left t -> t
681 _ -> eps
682 k = case hint of
683 Right t -> t
684 _ -> rankSVD tol a s
685 vs = takeColumns k v
686
687
688orth :: Field t => Matrix t -> [Vector t]
689-- ^ Return an orthonormal basis of the range space of a matrix
690orth m = take r $ toColumns u
691 where
692 (u,s,_) = compactSVD m
693 r = ranksv eps (max (rows m) (cols m)) (toList s)
694
695------------------------------------------------------------------------
696
697-- many thanks, quickcheck!
698
699haussholder :: (Field a) => a -> Vector a -> Matrix a
700haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
701 where w = asColumn v
702
703
704zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
705 where xs = toList v
706
707zt 0 v = v
708zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k]
709
710
711unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
712unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r)
713 where cs = toColumns pq
714 m = rows pq
715 n = cols pq
716 mn = min m n
717 r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
718 vs = zipWith zh [1..mn] cs
719 hs = zipWith haussholder (toList tau) vs
720 q = foldl1' mXm hs
721
722unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
723unpackHess hf m
724 | rows m == 1 = ((1><1)[1],m)
725 | otherwise = (uH . hf) m
726
727uH (pq, tau) = (p,h)
728 where cs = toColumns pq
729 m = rows pq
730 n = cols pq
731 mn = min m n
732 h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs
733 vs = zipWith zh [2..mn] cs
734 hs = zipWith haussholder (toList tau) vs
735 p = foldl1' mXm hs
736
737--------------------------------------------------------------------------
738
739-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
740rcond :: Field t => Matrix t -> Double
741rcond m = last s / head s
742 where s = toList (singularValues m)
743
744-- | Number of linearly independent rows or columns. See also 'ranksv'
745rank :: Field t => Matrix t -> Int
746rank m = rankSVD eps m (singularValues m)
747
748{-
749expm' m = case diagonalize (complex m) of
750 Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
751 Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
752 where exp = vectorMapC Exp
753-}
754
755diagonalize m = if rank v == n
756 then Just (l,v)
757 else Nothing
758 where n = rows m
759 (l,v) = if exactHermitian m
760 then let (l',v') = eigSH m in (real l', v')
761 else eig m
762
763-- | Generic matrix functions for diagonalizable matrices. For instance:
764--
765-- @logm = matFunc log@
766--
767matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
768matFunc f m = case diagonalize m of
769 Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v
770 Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"
771
772--------------------------------------------------------------
773
774golubeps :: Integer -> Integer -> Double
775golubeps p q = a * fromIntegral b / fromIntegral c where
776 a = 2^^(3-p-q)
777 b = fact p * fact q
778 c = fact (p+q) * fact (p+q+1)
779 fact n = product [1..n]
780
781epslist :: [(Int,Double)]
782epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
783
784geps delta = head [ k | (k,g) <- epslist, g<delta]
785
786
787{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
788 based on a scaled Pade approximation.
789-}
790expm :: Field t => Matrix t -> Matrix t
791expm = expGolub
792
793expGolub :: Field t => Matrix t -> Matrix t
794expGolub m = iterate msq f !! j
795 where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m
796 a = m */ fromIntegral ((2::Int)^j)
797 q = geps eps -- 7 steps
798 eye = ident (rows m)
799 work (k,c,x,n,d) = (k',c',x',n',d')
800 where k' = k+1
801 c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
802 x' = a <> x
803 n' = n |+| (c' .* x')
804 d' = d |+| (((-1)^k * c') .* x')
805 (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q
806 f = linearSolve df nf
807 msq x = x <> x
808
809 (<>) = multiply
810 v */ x = scale (recip x) v
811 (.*) = scale
812 (|+|) = add
813
814--------------------------------------------------------------
815
816{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
817It only works with invertible matrices that have a real solution.
818
819@m = (2><2) [4,9
820 ,0,4] :: Matrix Double@
821
822>>> sqrtm m
823(2><2)
824 [ 2.0, 2.25
825 , 0.0, 2.0 ]
826
827For diagonalizable matrices you can try 'matFunc' @sqrt@:
828
829>>> matFunc sqrt ((2><2) [1,0,0,-1])
830(2><2)
831 [ 1.0 :+ 0.0, 0.0 :+ 0.0
832 , 0.0 :+ 0.0, 0.0 :+ 1.0 ]
833
834-}
835sqrtm :: Field t => Matrix t -> Matrix t
836sqrtm = sqrtmInv
837
838sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
839 where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a
840 | otherwise = fixedPoint (b:rest)
841 fixedPoint _ = error "fixedpoint with impossible inputs"
842 f (y,z) = (0.5 .* (y |+| inv z),
843 0.5 .* (inv y |+| z))
844 (.*) = scale
845 (|+|) = add
846 (|-|) = sub
847
848------------------------------------------------------------------
849
850signlp r vals = foldl f 1 (zip [0..r-1] vals)
851 where f s (a,b) | a /= b = -s
852 | otherwise = s
853
854swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s)
855 | otherwise = (arr,s)
856
857fixPerm r vals = (fromColumns $ elems res, sign)
858 where v = [0..r-1]
859 s = toColumns (ident r)
860 (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
861
862triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
863 where el p q = if q-p>=h then v else 1 - v
864
865luFact (l_u,perm) | r <= c = (l ,u ,p, s)
866 | otherwise = (l',u',p, s)
867 where
868 r = rows l_u
869 c = cols l_u
870 tu = triang r c 0 1
871 tl = triang r c 0 0
872 l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r
873 u = l_u |*| tu
874 (p,s) = fixPerm r perm
875 l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c
876 u' = takeRows c (l_u |*| tu)
877 (|+|) = add
878 (|*|) = mul
879
880---------------------------------------------------------------------------
881
882data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
883
884class (RealFloat (RealOf t)) => Normed c t where
885 pnorm :: NormType -> c t -> RealOf t
886
887instance Normed Vector Double where
888 pnorm PNorm1 = norm1
889 pnorm PNorm2 = norm2
890 pnorm Infinity = normInf
891 pnorm Frobenius = norm2
892
893instance Normed Vector (Complex Double) where
894 pnorm PNorm1 = norm1
895 pnorm PNorm2 = norm2
896 pnorm Infinity = normInf
897 pnorm Frobenius = pnorm PNorm2
898
899instance Normed Vector Float where
900 pnorm PNorm1 = norm1
901 pnorm PNorm2 = norm2
902 pnorm Infinity = normInf
903 pnorm Frobenius = pnorm PNorm2
904
905instance Normed Vector (Complex Float) where
906 pnorm PNorm1 = norm1
907 pnorm PNorm2 = norm2
908 pnorm Infinity = normInf
909 pnorm Frobenius = pnorm PNorm2
910
911
912instance Normed Matrix Double where
913 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
914 pnorm PNorm2 = (@>0) . singularValues
915 pnorm Infinity = pnorm PNorm1 . trans
916 pnorm Frobenius = pnorm PNorm2 . flatten
917
918instance Normed Matrix (Complex Double) where
919 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
920 pnorm PNorm2 = (@>0) . singularValues
921 pnorm Infinity = pnorm PNorm1 . trans
922 pnorm Frobenius = pnorm PNorm2 . flatten
923
924instance Normed Matrix Float where
925 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
926 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
927 pnorm Infinity = pnorm PNorm1 . trans
928 pnorm Frobenius = pnorm PNorm2 . flatten
929
930instance Normed Matrix (Complex Float) where
931 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
932 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
933 pnorm Infinity = pnorm PNorm1 . trans
934 pnorm Frobenius = pnorm PNorm2 . flatten
935
936-- | Approximate number of common digits in the maximum element.
937relativeError' :: (Normed c t, Container c t) => c t -> c t -> Int
938relativeError' x y = dig (norm (x `sub` y) / norm x)
939 where norm = pnorm Infinity
940 dig r = round $ -logBase 10 (realToFrac r :: Double)
941
942
943relativeError :: Num a => (a -> Double) -> a -> a -> Double
944relativeError norm a b = r
945 where
946 na = norm a
947 nb = norm b
948 nab = norm (a-b)
949 mx = max na nb
950 mn = min na nb
951 r = if mn < peps
952 then mx
953 else nab/mx
954
955
956----------------------------------------------------------------------
957
958-- | Generalized symmetric positive definite eigensystem Av = lBv,
959-- for A and B symmetric, B positive definite (conditions not checked).
960geigSH' :: Field t
961 => Matrix t -- ^ A
962 -> Matrix t -- ^ B
963 -> (Vector Double, Matrix t)
964geigSH' a b = (l,v')
965 where
966 u = cholSH b
967 iu = inv u
968 c = ctrans iu <> a <> iu
969 (l,v) = eigSH' c
970 v' = iu <> v
971 (<>) = mXm
972
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-28 20:03:04 +0000