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