summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric/LinearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Numeric/LinearAlgebra')
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Algorithms.hs209
1 files changed, 193 insertions, 16 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index 21faf6e..25700bc 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -180,13 +180,97 @@ exactHermitian m = m `equal` ctrans m
180 180
181-------------------------------------------------------------- 181--------------------------------------------------------------
182 182
183-- | Full singular value decomposition. 183{- | Full singular value decomposition.
184
185@
186a = (5><3)
187 [ 1.0, 2.0, 3.0
188 , 4.0, 5.0, 6.0
189 , 7.0, 8.0, 9.0
190 , 10.0, 11.0, 12.0
191 , 13.0, 14.0, 15.0 ] :: Matrix Double
192@
193
194>>> let (u,s,v) = svd a
195
196>>> disp 3 u
1975x5
198-0.101 0.768 0.614 0.028 -0.149
199-0.249 0.488 -0.503 0.172 0.646
200-0.396 0.208 -0.405 -0.660 -0.449
201-0.543 -0.072 -0.140 0.693 -0.447
202-0.690 -0.352 0.433 -0.233 0.398
203
204>>> s
205fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
206
207>>> disp 3 v
2083x3
209-0.519 -0.751 0.408
210-0.576 -0.046 -0.816
211-0.632 0.659 0.408
212
213>>> let d = diagRect 0 s 5 3
214>>> disp 3 d
2155x3
21635.183 0.000 0.000
217 0.000 1.477 0.000
218 0.000 0.000 0.000
219 0.000 0.000 0.000
220
221>>> disp 3 $ u <> d <> tr v
2225x3
223 1.000 2.000 3.000
224 4.000 5.000 6.000
225 7.000 8.000 9.000
22610.000 11.000 12.000
22713.000 14.000 15.000
228
229-}
184svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 230svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
185svd = {-# SCC "svd" #-} svd' 231svd = {-# SCC "svd" #-} svd'
186 232
187-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@. 233{- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
188-- 234
189-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@. 235If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> tr v@.
236
237@
238a = (5><3)
239 [ 1.0, 2.0, 3.0
240 , 4.0, 5.0, 6.0
241 , 7.0, 8.0, 9.0
242 , 10.0, 11.0, 12.0
243 , 13.0, 14.0, 15.0 ] :: Matrix Double
244@
245
246>>> let (u,s,v) = thinSVD a
247
248>>> disp 3 u
2495x3
250-0.101 0.768 0.614
251-0.249 0.488 -0.503
252-0.396 0.208 -0.405
253-0.543 -0.072 -0.140
254-0.690 -0.352 0.433
255
256>>> s
257fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
258
259>>> disp 3 v
2603x3
261-0.519 -0.751 0.408
262-0.576 -0.046 -0.816
263-0.632 0.659 0.408
264
265>>> disp 3 $ u <> diag s <> tr v
2665x3
267 1.000 2.000 3.000
268 4.000 5.000 6.000
269 7.000 8.000 9.000
27010.000 11.000 12.000
27113.000 14.000 15.000
272
273-}
190thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 274thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
191thinSVD = {-# SCC "thinSVD" #-} thinSVD' 275thinSVD = {-# SCC "thinSVD" #-} thinSVD'
192 276
@@ -196,7 +280,7 @@ singularValues = {-# SCC "singularValues" #-} sv'
196 280
197-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. 281-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
198-- 282--
199-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. 283-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> tr v@.
200fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) 284fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
201fullSVD m = (u,d,v) where 285fullSVD m = (u,d,v) where
202 (u,s,v) = svd m 286 (u,s,v) = svd m
@@ -204,7 +288,47 @@ fullSVD m = (u,d,v) where
204 r = rows m 288 r = rows m
205 c = cols m 289 c = cols m
206 290
207-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors. 291{- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
292
293@
294a = (5><3)
295 [ 1.0, 2.0, 3.0
296 , 4.0, 5.0, 6.0
297 , 7.0, 8.0, 9.0
298 , 10.0, 11.0, 12.0
299 , 13.0, 14.0, 15.0 ] :: Matrix Double
300@
301
302>>> let (u,s,v) = compactSVD a
303
304>>> disp 3 u
3055x2
306-0.101 0.768
307-0.249 0.488
308-0.396 0.208
309-0.543 -0.072
310-0.690 -0.352
311
312>>> s
313fromList [35.18264833189422,1.4769076999800903]
314
315>>> disp 3 u
3165x2
317-0.101 0.768
318-0.249 0.488
319-0.396 0.208
320-0.543 -0.072
321-0.690 -0.352
322
323>>> disp 3 $ u <> diag s <> tr v
3245x3
325 1.000 2.000 3.000
326 4.000 5.000 6.000
327 7.000 8.000 9.000
32810.000 11.000 12.000
32913.000 14.000 15.000
330
331-}
208compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 332compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
209compactSVD m = (u', subVector 0 d s, v') where 333compactSVD m = (u', subVector 0 d s, v') where
210 (u,s,v) = thinSVD m 334 (u,s,v) = thinSVD m
@@ -213,12 +337,12 @@ compactSVD m = (u', subVector 0 d s, v') where
213 v' = takeColumns d v 337 v' = takeColumns d v
214 338
215 339
216-- | Singular values and all right singular vectors. 340-- | Singular values and all right singular vectors (as columns).
217rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) 341rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
218rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) 342rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
219 | otherwise = let (_,s,v) = svd m in (s,v) 343 | otherwise = let (_,s,v) = svd m in (s,v)
220 344
221-- | Singular values and all left singular vectors. 345-- | Singular values and all left singular vectors (as columns).
222leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) 346leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
223leftSV m | vertical m = let (u,s,_) = svd m in (u,s) 347leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
224 | otherwise = let (u,s,_) = thinSVD m in (u,s) 348 | otherwise = let (u,s,_) = thinSVD m in (u,s)
@@ -258,13 +382,46 @@ linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
258 382
259-------------------------------------------------------------- 383--------------------------------------------------------------
260 384
261-- | Eigenvalues and eigenvectors of a general square matrix. 385{- | Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.
262-- 386
263-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ 387If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
388
389@
390a = (3><3)
391 [ 3, 0, -2
392 , 4, 5, -1
393 , 3, 1, 0 ] :: Matrix Double
394@
395
396>>> let (l, v) = eig a
397
398>>> putStr . dispcf 3 . asRow $ l
3991x3
4001.925+1.523i 1.925-1.523i 4.151
401
402>>> putStr . dispcf 3 $ v
4033x3
404-0.455+0.365i -0.455-0.365i 0.181
405 0.603 0.603 -0.978
406 0.033+0.543i 0.033-0.543i -0.104
407
408>>> putStr . dispcf 3 $ complex a <> v
4093x3
410-1.432+0.010i -1.432-0.010i 0.753
411 1.160+0.918i 1.160-0.918i -4.059
412-0.763+1.096i -0.763-1.096i -0.433
413
414>>> putStr . dispcf 3 $ v <> diag l
4153x3
416-1.432+0.010i -1.432-0.010i 0.753
417 1.160+0.918i 1.160-0.918i -4.059
418-0.763+1.096i -0.763-1.096i -0.433
419
420-}
264eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 421eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
265eig = {-# SCC "eig" #-} eig' 422eig = {-# SCC "eig" #-} eig'
266 423
267-- | Eigenvalues of a general square matrix. 424-- | Eigenvalues (not ordered) of a general square matrix.
268eigenvalues :: Field t => Matrix t -> Vector (Complex Double) 425eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
269eigenvalues = {-# SCC "eigenvalues" #-} eigOnly 426eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
270 427
@@ -276,14 +433,34 @@ eigSH' = {-# SCC "eigSH'" #-} eigSH''
276eigenvaluesSH' :: Field t => Matrix t -> Vector Double 433eigenvaluesSH' :: Field t => Matrix t -> Vector Double
277eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH 434eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH
278 435
279-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix. 436{- | Eigenvalues and eigenvectors (as columns) of a complex hermitian or real symmetric matrix, in descending order.
280-- 437
281-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ 438If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@
439
440@
441a = (3><3)
442 [ 1.0, 2.0, 3.0
443 , 2.0, 4.0, 5.0
444 , 3.0, 5.0, 6.0 ]
445@
446
447>>> let (l, v) = eigSH a
448
449>>> l
450fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
451
452>>> disp 3 $ v <> diag l <> tr v
4533x3
4541.000 2.000 3.000
4552.000 4.000 5.000
4563.000 5.000 6.000
457
458-}
282eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) 459eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
283eigSH m | exactHermitian m = eigSH' m 460eigSH m | exactHermitian m = eigSH' m
284 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" 461 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
285 462
286-- | Eigenvalues of a complex hermitian or real symmetric matrix. 463-- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
287eigenvaluesSH :: Field t => Matrix t -> Vector Double 464eigenvaluesSH :: Field t => Matrix t -> Vector Double
288eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m 465eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
289 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix" 466 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"