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