summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Numeric')
-rw-r--r--packages/base/src/Numeric/Container.hs242
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Algorithms.hs744
2 files changed, 986 insertions, 0 deletions
diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs
new file mode 100644
index 0000000..b7d3b80
--- /dev/null
+++ b/packages/base/src/Numeric/Container.hs
@@ -0,0 +1,242 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Numeric.Container
11-- Copyright : (c) Alberto Ruiz 2010-14
12-- License : GPL
13--
14-- Maintainer : Alberto Ruiz <aruiz@um.es>
15-- Stability : provisional
16-- Portability : portable
17--
18-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
19--
20-- The 'Container' class is used to define optimized generic functions which work
21-- on 'Vector' and 'Matrix' with real or complex elements.
22--
23-- Some of these functions are also available in the instances of the standard
24-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
25--
26-----------------------------------------------------------------------------
27{-# OPTIONS_HADDOCK hide #-}
28
29module Numeric.Container (
30 -- * Basic functions
31 module Data.Packed,
32 konst, build,
33 linspace,
34 diag, ident,
35 ctrans,
36 -- * Generic operations
37 Container(..),
38 -- * Matrix product
39 Product(..), udot, dot, (◇),
40 Mul(..),
41 Contraction(..),
42 optimiseMult,
43 mXm,mXv,vXm,LSDiv(..),
44 outer, kronecker,
45 -- * Element conversion
46 Convert(..),
47 Complexable(),
48 RealElement(),
49
50 RealOf, ComplexOf, SingleOf, DoubleOf,
51
52 IndexOf,
53 module Data.Complex
54) where
55
56import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
57import Data.Packed.Numeric
58import Data.Complex
59import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
60import Data.Monoid(Monoid(mconcat))
61
62------------------------------------------------------------------
63
64{- | Creates a real vector containing a range of values:
65
66>>> linspace 5 (-3,7::Double)
67fromList [-3.0,-0.5,2.0,4.5,7.0]@
68
69>>> linspace 5 (8,2+i) :: Vector (Complex Double)
70fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
71
72Logarithmic spacing can be defined as follows:
73
74@logspace n (a,b) = 10 ** linspace n (a,b)@
75-}
76linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
77linspace 0 (a,b) = fromList[(a+b)/2]
78linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
79 where s = (b-a)/fromIntegral (n-1)
80
81--------------------------------------------------------
82
83class Contraction a b c | a b -> c
84 where
85 infixl 7 <.>
86 {- | Matrix product, matrix vector product, and dot product
87
88Examples:
89
90>>> let a = (3><4) [1..] :: Matrix Double
91>>> let v = fromList [1,0,2,-1] :: Vector Double
92>>> let u = fromList [1,2,3] :: Vector Double
93
94>>> a
95(3><4)
96 [ 1.0, 2.0, 3.0, 4.0
97 , 5.0, 6.0, 7.0, 8.0
98 , 9.0, 10.0, 11.0, 12.0 ]
99
100matrix × matrix:
101
102>>> disp 2 (a <.> trans a)
1033x3
104 30 70 110
105 70 174 278
106110 278 446
107
108matrix × vector:
109
110>>> a <.> v
111fromList [3.0,11.0,19.0]
112
113dot product:
114
115>>> u <.> fromList[3,2,1::Double]
11610
117
118For complex vectors the first argument is conjugated:
119
120>>> fromList [1,i] <.> fromList[2*i+1,3]
1211.0 :+ (-1.0)
122
123>>> fromList [1,i,1-i] <.> complex a
124fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
125
126-}
127 (<.>) :: a -> b -> c
128
129
130instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
131 u <.> v = conj u `udot` v
132
133instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
134 (<.>) = mXv
135
136instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
137 (<.>) v m = (conj v) `vXm` m
138
139instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
140 (<.>) = mXm
141
142
143--------------------------------------------------------------------------------
144
145class Mul a b c | a b -> c where
146 infixl 7 <>
147 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
148 (<>) :: Product t => a t -> b t -> c t
149
150instance Mul Matrix Matrix Matrix where
151 (<>) = mXm
152
153instance Mul Matrix Vector Vector where
154 (<>) m v = flatten $ m <> asColumn v
155
156instance Mul Vector Matrix Vector where
157 (<>) v m = flatten $ asRow v <> m
158
159--------------------------------------------------------------------------------
160
161class LSDiv c where
162 infixl 7 <\>
163 -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
164 (<\>) :: Field t => Matrix t -> c t -> c t
165
166instance LSDiv Vector where
167 m <\> v = flatten (linearSolveSVD m (reshape 1 v))
168
169instance LSDiv Matrix where
170 (<\>) = linearSolveSVD
171
172--------------------------------------------------------------------------------
173
174class Konst e d c | d -> c, c -> d
175 where
176 -- |
177 -- >>> konst 7 3 :: Vector Float
178 -- fromList [7.0,7.0,7.0]
179 --
180 -- >>> konst i (3::Int,4::Int)
181 -- (3><4)
182 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
183 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
184 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
185 --
186 konst :: e -> d -> c e
187
188instance Container Vector e => Konst e Int Vector
189 where
190 konst = konst'
191
192instance Container Vector e => Konst e (Int,Int) Matrix
193 where
194 konst = konst'
195
196--------------------------------------------------------------------------------
197
198class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
199 where
200 -- |
201 -- >>> build 5 (**2) :: Vector Double
202 -- fromList [0.0,1.0,4.0,9.0,16.0]
203 --
204 -- Hilbert matrix of order N:
205 --
206 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
207 -- >>> putStr . dispf 2 $ hilb 3
208 -- 3x3
209 -- 1.00 0.50 0.33
210 -- 0.50 0.33 0.25
211 -- 0.33 0.25 0.20
212 --
213 build :: d -> f -> c e
214
215instance Container Vector e => Build Int (e -> e) Vector e
216 where
217 build = build'
218
219instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
220 where
221 build = build'
222
223--------------------------------------------------------------------------------
224
225{- | alternative operator for '(\<.\>)'
226
227x25c7, white diamond
228
229-}
230(◇) :: Contraction a b c => a -> b -> c
231infixl 7 ◇
232(◇) = (<.>)
233
234-- | dot product: @cdot u v = 'udot' ('conj' u) v@
235dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
236dot u v = udot (conj u) v
237
238--------------------------------------------------------------------------------
239
240optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
241optimiseMult = mconcat
242
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-07 21:09:24 +0000