diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-06-20 11:31:24 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-06-20 11:31:24 +0200 |
commit | 21a70f76d47d4851f27161eab350a31ab0651047 (patch) | |
tree | faf81be978dc0789b47c43115c633b2b2e4b14c5 | |
parent | 2869503ccb552238e330562a62a076e48d567f79 (diff) |
documentation
-rw-r--r-- | packages/base/hmatrix.cabal | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Algorithms.hs | 209 |
2 files changed, 194 insertions, 17 deletions
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index c6c421b..0eeb6f0 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal | |||
@@ -7,7 +7,7 @@ Maintainer: Alberto Ruiz | |||
7 | Stability: provisional | 7 | Stability: provisional |
8 | Homepage: https://github.com/albertoruiz/hmatrix | 8 | Homepage: https://github.com/albertoruiz/hmatrix |
9 | Synopsis: Numeric Linear Algebra | 9 | Synopsis: Numeric Linear Algebra |
10 | Description: Dense linear algebra based on BLAS and LAPACK. | 10 | Description: Linear algebra based on BLAS and LAPACK. |
11 | Category: Math | 11 | Category: Math |
12 | tested-with: GHC==7.8 | 12 | tested-with: GHC==7.8 |
13 | 13 | ||
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 | @ | ||
186 | a = (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 | ||
197 | 5x5 | ||
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 | ||
205 | fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15] | ||
206 | |||
207 | >>> disp 3 v | ||
208 | 3x3 | ||
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 | ||
215 | 5x3 | ||
216 | 35.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 | ||
222 | 5x3 | ||
223 | 1.000 2.000 3.000 | ||
224 | 4.000 5.000 6.000 | ||
225 | 7.000 8.000 9.000 | ||
226 | 10.000 11.000 12.000 | ||
227 | 13.000 14.000 15.000 | ||
228 | |||
229 | -} | ||
184 | svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | 230 | svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) |
185 | svd = {-# SCC "svd" #-} svd' | 231 | svd = {-# 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@. | 235 | If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> tr v@. |
236 | |||
237 | @ | ||
238 | a = (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 | ||
249 | 5x3 | ||
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 | ||
257 | fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15] | ||
258 | |||
259 | >>> disp 3 v | ||
260 | 3x3 | ||
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 | ||
266 | 5x3 | ||
267 | 1.000 2.000 3.000 | ||
268 | 4.000 5.000 6.000 | ||
269 | 7.000 8.000 9.000 | ||
270 | 10.000 11.000 12.000 | ||
271 | 13.000 14.000 15.000 | ||
272 | |||
273 | -} | ||
190 | thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | 274 | thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) |
191 | thinSVD = {-# SCC "thinSVD" #-} thinSVD' | 275 | thinSVD = {-# 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@. |
200 | fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) | 284 | fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) |
201 | fullSVD m = (u,d,v) where | 285 | fullSVD 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 | @ | ||
294 | a = (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 | ||
305 | 5x2 | ||
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 | ||
313 | fromList [35.18264833189422,1.4769076999800903] | ||
314 | |||
315 | >>> disp 3 u | ||
316 | 5x2 | ||
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 | ||
324 | 5x3 | ||
325 | 1.000 2.000 3.000 | ||
326 | 4.000 5.000 6.000 | ||
327 | 7.000 8.000 9.000 | ||
328 | 10.000 11.000 12.000 | ||
329 | 13.000 14.000 15.000 | ||
330 | |||
331 | -} | ||
208 | compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) | 332 | compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) |
209 | compactSVD m = (u', subVector 0 d s, v') where | 333 | compactSVD 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). |
217 | rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) | 341 | rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) |
218 | rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) | 342 | rightSV 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). |
222 | leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) | 346 | leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) |
223 | leftSV m | vertical m = let (u,s,_) = svd m in (u,s) | 347 | leftSV 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@ | 387 | If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ |
388 | |||
389 | @ | ||
390 | a = (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 | ||
399 | 1x3 | ||
400 | 1.925+1.523i 1.925-1.523i 4.151 | ||
401 | |||
402 | >>> putStr . dispcf 3 $ v | ||
403 | 3x3 | ||
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 | ||
409 | 3x3 | ||
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 | ||
415 | 3x3 | ||
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 | -} | ||
264 | eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) | 421 | eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) |
265 | eig = {-# SCC "eig" #-} eig' | 422 | eig = {-# SCC "eig" #-} eig' |
266 | 423 | ||
267 | -- | Eigenvalues of a general square matrix. | 424 | -- | Eigenvalues (not ordered) of a general square matrix. |
268 | eigenvalues :: Field t => Matrix t -> Vector (Complex Double) | 425 | eigenvalues :: Field t => Matrix t -> Vector (Complex Double) |
269 | eigenvalues = {-# SCC "eigenvalues" #-} eigOnly | 426 | eigenvalues = {-# SCC "eigenvalues" #-} eigOnly |
270 | 427 | ||
@@ -276,14 +433,34 @@ eigSH' = {-# SCC "eigSH'" #-} eigSH'' | |||
276 | eigenvaluesSH' :: Field t => Matrix t -> Vector Double | 433 | eigenvaluesSH' :: Field t => Matrix t -> Vector Double |
277 | eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH | 434 | eigenvaluesSH' = {-# 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@ | 438 | If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@ |
439 | |||
440 | @ | ||
441 | a = (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 | ||
450 | fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575] | ||
451 | |||
452 | >>> disp 3 $ v <> diag l <> tr v | ||
453 | 3x3 | ||
454 | 1.000 2.000 3.000 | ||
455 | 2.000 4.000 5.000 | ||
456 | 3.000 5.000 6.000 | ||
457 | |||
458 | -} | ||
282 | eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) | 459 | eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) |
283 | eigSH m | exactHermitian m = eigSH' m | 460 | eigSH 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. |
287 | eigenvaluesSH :: Field t => Matrix t -> Vector Double | 464 | eigenvaluesSH :: Field t => Matrix t -> Vector Double |
288 | eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m | 465 | eigenvaluesSH 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" |