diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-16 12:36:52 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-16 12:36:52 +0200 |
commit | a2d99e7d0e83fcedf3a856cdb927309e28a8eddd (patch) | |
tree | c6192d21e4b59e5527aac5d25d8b8baefc2b3052 /packages/base | |
parent | 5d15d765c0204854b587b03232b20336ddf91ced (diff) |
container and algorithms moved to base
Diffstat (limited to 'packages/base')
-rw-r--r-- | packages/base/hmatrix-base.cabal | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/Container.hs | 242 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Algorithms.hs | 744 |
3 files changed, 988 insertions, 0 deletions
diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 8eeb97e..9a84b7e 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal | |||
@@ -37,10 +37,12 @@ library | |||
37 | Data.Packed.Development, | 37 | Data.Packed.Development, |
38 | Numeric.Conversion | 38 | Numeric.Conversion |
39 | Numeric.LinearAlgebra.LAPACK | 39 | Numeric.LinearAlgebra.LAPACK |
40 | Numeric.LinearAlgebra.Algorithms | ||
40 | Data.Packed.Numeric | 41 | Data.Packed.Numeric |
41 | Numeric.Vectorized | 42 | Numeric.Vectorized |
42 | Numeric.Vector | 43 | Numeric.Vector |
43 | Numeric.Matrix | 44 | Numeric.Matrix |
45 | Numeric.Container | ||
44 | 46 | ||
45 | other-modules: Data.Packed.Internal, | 47 | other-modules: Data.Packed.Internal, |
46 | Data.Packed.Internal.Common, | 48 | Data.Packed.Internal.Common, |
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 | |||
29 | module 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 | |||
56 | import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) | ||
57 | import Data.Packed.Numeric | ||
58 | import Data.Complex | ||
59 | import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) | ||
60 | import Data.Monoid(Monoid(mconcat)) | ||
61 | |||
62 | ------------------------------------------------------------------ | ||
63 | |||
64 | {- | Creates a real vector containing a range of values: | ||
65 | |||
66 | >>> linspace 5 (-3,7::Double) | ||
67 | fromList [-3.0,-0.5,2.0,4.5,7.0]@ | ||
68 | |||
69 | >>> linspace 5 (8,2+i) :: Vector (Complex Double) | ||
70 | fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] | ||
71 | |||
72 | Logarithmic spacing can be defined as follows: | ||
73 | |||
74 | @logspace n (a,b) = 10 ** linspace n (a,b)@ | ||
75 | -} | ||
76 | linspace :: (Container Vector e) => Int -> (e, e) -> Vector e | ||
77 | linspace 0 (a,b) = fromList[(a+b)/2] | ||
78 | linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] | ||
79 | where s = (b-a)/fromIntegral (n-1) | ||
80 | |||
81 | -------------------------------------------------------- | ||
82 | |||
83 | class Contraction a b c | a b -> c | ||
84 | where | ||
85 | infixl 7 <.> | ||
86 | {- | Matrix product, matrix vector product, and dot product | ||
87 | |||
88 | Examples: | ||
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 | |||
100 | matrix × matrix: | ||
101 | |||
102 | >>> disp 2 (a <.> trans a) | ||
103 | 3x3 | ||
104 | 30 70 110 | ||
105 | 70 174 278 | ||
106 | 110 278 446 | ||
107 | |||
108 | matrix × vector: | ||
109 | |||
110 | >>> a <.> v | ||
111 | fromList [3.0,11.0,19.0] | ||
112 | |||
113 | dot product: | ||
114 | |||
115 | >>> u <.> fromList[3,2,1::Double] | ||
116 | 10 | ||
117 | |||
118 | For complex vectors the first argument is conjugated: | ||
119 | |||
120 | >>> fromList [1,i] <.> fromList[2*i+1,3] | ||
121 | 1.0 :+ (-1.0) | ||
122 | |||
123 | >>> fromList [1,i,1-i] <.> complex a | ||
124 | fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] | ||
125 | |||
126 | -} | ||
127 | (<.>) :: a -> b -> c | ||
128 | |||
129 | |||
130 | instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where | ||
131 | u <.> v = conj u `udot` v | ||
132 | |||
133 | instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where | ||
134 | (<.>) = mXv | ||
135 | |||
136 | instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where | ||
137 | (<.>) v m = (conj v) `vXm` m | ||
138 | |||
139 | instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where | ||
140 | (<.>) = mXm | ||
141 | |||
142 | |||
143 | -------------------------------------------------------------------------------- | ||
144 | |||
145 | class 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 | |||
150 | instance Mul Matrix Matrix Matrix where | ||
151 | (<>) = mXm | ||
152 | |||
153 | instance Mul Matrix Vector Vector where | ||
154 | (<>) m v = flatten $ m <> asColumn v | ||
155 | |||
156 | instance Mul Vector Matrix Vector where | ||
157 | (<>) v m = flatten $ asRow v <> m | ||
158 | |||
159 | -------------------------------------------------------------------------------- | ||
160 | |||
161 | class 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 | |||
166 | instance LSDiv Vector where | ||
167 | m <\> v = flatten (linearSolveSVD m (reshape 1 v)) | ||
168 | |||
169 | instance LSDiv Matrix where | ||
170 | (<\>) = linearSolveSVD | ||
171 | |||
172 | -------------------------------------------------------------------------------- | ||
173 | |||
174 | class 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 | |||
188 | instance Container Vector e => Konst e Int Vector | ||
189 | where | ||
190 | konst = konst' | ||
191 | |||
192 | instance Container Vector e => Konst e (Int,Int) Matrix | ||
193 | where | ||
194 | konst = konst' | ||
195 | |||
196 | -------------------------------------------------------------------------------- | ||
197 | |||
198 | class 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 | |||
215 | instance Container Vector e => Build Int (e -> e) Vector e | ||
216 | where | ||
217 | build = build' | ||
218 | |||
219 | instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e | ||
220 | where | ||
221 | build = build' | ||
222 | |||
223 | -------------------------------------------------------------------------------- | ||
224 | |||
225 | {- | alternative operator for '(\<.\>)' | ||
226 | |||
227 | x25c7, white diamond | ||
228 | |||
229 | -} | ||
230 | (◇) :: Contraction a b c => a -> b -> c | ||
231 | infixl 7 ◇ | ||
232 | (◇) = (<.>) | ||
233 | |||
234 | -- | dot product: @cdot u v = 'udot' ('conj' u) v@ | ||
235 | dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t | ||
236 | dot u v = udot (conj u) v | ||
237 | |||
238 | -------------------------------------------------------------------------------- | ||
239 | |||
240 | optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t | ||
241 | optimiseMult = 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 | {- | | ||
9 | Module : Numeric.LinearAlgebra.Algorithms | ||
10 | Copyright : (c) Alberto Ruiz 2006-9 | ||
11 | License : GPL-style | ||
12 | |||
13 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
14 | Stability : provisional | ||
15 | Portability : uses ffi | ||
16 | |||
17 | High level generic interface to common matrix computations. | ||
18 | |||
19 | Specific functions for particular base types can also be explicitly | ||
20 | imported from "Numeric.LinearAlgebra.LAPACK". | ||
21 | |||
22 | -} | ||
23 | ----------------------------------------------------------------------------- | ||
24 | |||
25 | module 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 | |||
80 | import Data.Packed.Development hiding ((//)) | ||
81 | import Data.Packed | ||
82 | import Numeric.LinearAlgebra.LAPACK as LAPACK | ||
83 | import Data.List(foldl1') | ||
84 | import Data.Array | ||
85 | import 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 | -} | ||
93 | class (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 | |||
122 | instance 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 | |||
143 | instance 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 | |||
171 | square m = rows m == cols m | ||
172 | |||
173 | vertical m = rows m >= cols m | ||
174 | |||
175 | exactHermitian m = m `equal` ctrans m | ||
176 | |||
177 | -------------------------------------------------------------- | ||
178 | |||
179 | -- | Full singular value decomposition. | ||
180 | svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | ||
181 | svd = {-# 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@. | ||
186 | thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | ||
187 | thinSVD = {-# SCC "thinSVD" #-} thinSVD' | ||
188 | |||
189 | -- | Singular values only. | ||
190 | singularValues :: Field t => Matrix t -> Vector Double | ||
191 | singularValues = {-# 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@. | ||
196 | fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) | ||
197 | fullSVD 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. | ||
204 | compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | ||
205 | compactSVD 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. | ||
213 | rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) | ||
214 | rightSV 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. | ||
218 | leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) | ||
219 | leftSV 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'. | ||
226 | luPacked :: Field t => Matrix t -> (Matrix t, [Int]) | ||
227 | luPacked = {-# SCC "luPacked" #-} luPacked' | ||
228 | |||
229 | -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'. | ||
230 | luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t | ||
231 | luSolve = {-# 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. | ||
235 | linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t | ||
236 | linearSolve = {-# SCC "linearSolve" #-} linearSolve' | ||
237 | |||
238 | -- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'. | ||
239 | cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t | ||
240 | cholSolve = {-# 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. | ||
243 | linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t | ||
244 | linearSolveSVD = {-# 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'. | ||
248 | linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t | ||
249 | linearSolveLS = {-# 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@ | ||
256 | eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) | ||
257 | eig = {-# SCC "eig" #-} eig' | ||
258 | |||
259 | -- | Eigenvalues of a general square matrix. | ||
260 | eigenvalues :: Field t => Matrix t -> Vector (Complex Double) | ||
261 | eigenvalues = {-# 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. | ||
264 | eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) | ||
265 | eigSH' = {-# 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. | ||
268 | eigenvaluesSH' :: Field t => Matrix t -> Vector Double | ||
269 | eigenvaluesSH' = {-# 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@ | ||
274 | eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) | ||
275 | eigSH 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. | ||
279 | eigenvaluesSH :: Field t => Matrix t -> Vector Double | ||
280 | eigenvaluesSH 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. | ||
288 | qr :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
289 | qr = {-# SCC "qr" #-} unpackQR . qr' | ||
290 | |||
291 | qrRaw m = qr' m | ||
292 | |||
293 | {- | generate a matrix with k orthogonal columns from the output of qrRaw | ||
294 | -} | ||
295 | qrgr 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. | ||
302 | rq :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
303 | rq 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). | ||
314 | hess :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
315 | hess = 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) | ||
325 | schur :: Field t => Matrix t -> (Matrix t, Matrix t) | ||
326 | schur = schur' | ||
327 | |||
328 | |||
329 | -- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. | ||
330 | mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) | ||
331 | mbCholSH = {-# 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. | ||
334 | cholSH :: Field t => Matrix t -> Matrix t | ||
335 | cholSH = {-# 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@. | ||
340 | chol :: Field t => Matrix t -> Matrix t | ||
341 | chol 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. | ||
346 | invlndet :: Field t | ||
347 | => Matrix t | ||
348 | -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) | ||
349 | invlndet 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'. | ||
361 | det :: Field t => Matrix t -> t | ||
362 | det 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. | ||
371 | lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) | ||
372 | lu = luFact . luPacked | ||
373 | |||
374 | -- | Inverse of a square matrix. See also 'invlndet'. | ||
375 | inv :: Field t => Matrix t -> Matrix t | ||
376 | inv 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). | ||
381 | pinv :: Field t => Matrix t -> Matrix t | ||
382 | pinv = 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 | @ | ||
387 | m = (3><3) [ 1, 0, 0 | ||
388 | , 0, 1, 0 | ||
389 | , 0, 0, 1e-10] :: Matrix Double | ||
390 | @ | ||
391 | |||
392 | >>> pinv m | ||
393 | 1. 0. 0. | ||
394 | 0. 1. 0. | ||
395 | 0. 0. 10000000000. | ||
396 | |||
397 | >>> pinvTol 1E8 m | ||
398 | 1. 0. 0. | ||
399 | 0. 1. 0. | ||
400 | 0. 0. 1. | ||
401 | |||
402 | -} | ||
403 | |||
404 | pinvTol :: Field t => Double -> Matrix t -> Matrix t | ||
405 | pinvTol 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. | ||
419 | rankSVD :: 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 | ||
424 | rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s) | ||
425 | |||
426 | -- | Numeric rank of a matrix from its singular values. | ||
427 | ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') | ||
428 | -> Int -- ^ maximum dimension of the matrix | ||
429 | -> [Double] -- ^ singular values | ||
430 | -> Int -- ^ rank of m | ||
431 | ranksv 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). | ||
438 | eps :: Double | ||
439 | eps = 2.22044604925031e-16 | ||
440 | |||
441 | |||
442 | -- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 | ||
443 | peps :: RealFloat x => x | ||
444 | peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) | ||
445 | |||
446 | |||
447 | -- | The imaginary unit: @i = 0.0 :+ 1.0@ | ||
448 | i :: Complex Double | ||
449 | i = 0:+1 | ||
450 | |||
451 | ----------------------------------------------------------------------- | ||
452 | |||
453 | -- | The nullspace of a matrix from its precomputed SVD decomposition. | ||
454 | nullspaceSVD :: 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 | ||
460 | nullspaceSVD 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'. | ||
471 | nullspacePrec :: 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 | ||
475 | nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) | ||
476 | |||
477 | -- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. | ||
478 | nullVector :: Field t => Matrix t -> Vector t | ||
479 | nullVector = last . nullspacePrec 1 | ||
480 | |||
481 | orth :: Field t => Matrix t -> [Vector t] | ||
482 | -- ^ Return an orthonormal basis of the range space of a matrix | ||
483 | orth 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 | |||
492 | haussholder :: (Field a) => a -> Vector a -> Matrix a | ||
493 | haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) | ||
494 | where w = asColumn v | ||
495 | |||
496 | |||
497 | zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs) | ||
498 | where xs = toList v | ||
499 | |||
500 | zt 0 v = v | ||
501 | zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k] | ||
502 | |||
503 | |||
504 | unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) | ||
505 | unpackQR (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 | |||
515 | unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) | ||
516 | unpackHess hf m | ||
517 | | rows m == 1 = ((1><1)[1],m) | ||
518 | | otherwise = (uH . hf) m | ||
519 | |||
520 | uH (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. | ||
533 | rcond :: Field t => Matrix t -> Double | ||
534 | rcond m = last s / head s | ||
535 | where s = toList (singularValues m) | ||
536 | |||
537 | -- | Number of linearly independent rows or columns. | ||
538 | rank :: Field t => Matrix t -> Int | ||
539 | rank m = rankSVD eps m (singularValues m) | ||
540 | |||
541 | {- | ||
542 | expm' 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 | |||
548 | diagonalize 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 | -- | ||
560 | matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
561 | matFunc 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 | |||
567 | golubeps :: Integer -> Integer -> Double | ||
568 | golubeps 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 | |||
574 | epslist :: [(Int,Double)] | ||
575 | epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] | ||
576 | |||
577 | geps 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 | -} | ||
583 | expm :: Field t => Matrix t -> Matrix t | ||
584 | expm = expGolub | ||
585 | |||
586 | expGolub :: Field t => Matrix t -> Matrix t | ||
587 | expGolub 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. | ||
610 | It 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 | -} | ||
621 | sqrtm :: Field t => Matrix t -> Matrix t | ||
622 | sqrtm = sqrtmInv | ||
623 | |||
624 | sqrtmInv 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 | |||
636 | signlp r vals = foldl f 1 (zip [0..r-1] vals) | ||
637 | where f s (a,b) | a /= b = -s | ||
638 | | otherwise = s | ||
639 | |||
640 | swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) | ||
641 | | otherwise = (arr,s) | ||
642 | |||
643 | fixPerm 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 | |||
648 | triang 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 | |||
651 | luFact (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 | |||
668 | data NormType = Infinity | PNorm1 | PNorm2 | Frobenius | ||
669 | |||
670 | class (RealFloat (RealOf t)) => Normed c t where | ||
671 | pnorm :: NormType -> c t -> RealOf t | ||
672 | |||
673 | instance Normed Vector Double where | ||
674 | pnorm PNorm1 = norm1 | ||
675 | pnorm PNorm2 = norm2 | ||
676 | pnorm Infinity = normInf | ||
677 | pnorm Frobenius = norm2 | ||
678 | |||
679 | instance Normed Vector (Complex Double) where | ||
680 | pnorm PNorm1 = norm1 | ||
681 | pnorm PNorm2 = norm2 | ||
682 | pnorm Infinity = normInf | ||
683 | pnorm Frobenius = pnorm PNorm2 | ||
684 | |||
685 | instance Normed Vector Float where | ||
686 | pnorm PNorm1 = norm1 | ||
687 | pnorm PNorm2 = norm2 | ||
688 | pnorm Infinity = normInf | ||
689 | pnorm Frobenius = pnorm PNorm2 | ||
690 | |||
691 | instance Normed Vector (Complex Float) where | ||
692 | pnorm PNorm1 = norm1 | ||
693 | pnorm PNorm2 = norm2 | ||
694 | pnorm Infinity = normInf | ||
695 | pnorm Frobenius = pnorm PNorm2 | ||
696 | |||
697 | |||
698 | instance 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 | |||
704 | instance 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 | |||
710 | instance 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 | |||
716 | instance 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. | ||
723 | relativeError :: (Normed c t, Container c t) => c t -> c t -> Int | ||
724 | relativeError 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). | ||
732 | geigSH' :: Field t | ||
733 | => Matrix t -- ^ A | ||
734 | -> Matrix t -- ^ B | ||
735 | -> (Vector Double, Matrix t) | ||
736 | geigSH' 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 | |||