diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 90 |
1 files changed, 44 insertions, 46 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 6500e0b..4823cec 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -31,7 +31,7 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
31 | cholSolve, | 31 | cholSolve, |
32 | linearSolveLS, | 32 | linearSolveLS, |
33 | linearSolveSVD, | 33 | linearSolveSVD, |
34 | inv, pinv, | 34 | inv, pinv, pinvTol, |
35 | det, invlndet, | 35 | det, invlndet, |
36 | rank, rcond, | 36 | rank, rcond, |
37 | -- * Matrix factorizations | 37 | -- * Matrix factorizations |
@@ -73,7 +73,6 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
73 | -- * Util | 73 | -- * Util |
74 | haussholder, | 74 | haussholder, |
75 | unpackQR, unpackHess, | 75 | unpackQR, unpackHess, |
76 | pinvTol, | ||
77 | ranksv | 76 | ranksv |
78 | ) where | 77 | ) where |
79 | 78 | ||
@@ -95,7 +94,9 @@ class (Product t, | |||
95 | Container Vector t, | 94 | Container Vector t, |
96 | Container Matrix t, | 95 | Container Matrix t, |
97 | Normed Matrix t, | 96 | Normed Matrix t, |
98 | Normed Vector t) => Field t where | 97 | Normed Vector t, |
98 | Floating t, | ||
99 | RealOf t ~ Double) => Field t where | ||
99 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 100 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
100 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 101 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
101 | sv' :: Matrix t -> Vector Double | 102 | sv' :: Matrix t -> Vector Double |
@@ -330,9 +331,9 @@ chol m | exactHermitian m = cholSH m | |||
330 | 331 | ||
331 | 332 | ||
332 | -- | Joint computation of inverse and logarithm of determinant of a square matrix. | 333 | -- | Joint computation of inverse and logarithm of determinant of a square matrix. |
333 | invlndet :: (Floating t, Field t) | 334 | invlndet :: Field t |
334 | => Matrix t | 335 | => Matrix t |
335 | -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) | 336 | -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) |
336 | invlndet m | square m = (im,(ladm,sdm)) | 337 | invlndet m | square m = (im,(ladm,sdm)) |
337 | | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" | 338 | | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" |
338 | where | 339 | where |
@@ -363,9 +364,42 @@ inv :: Field t => Matrix t -> Matrix t | |||
363 | inv m | square m = m `linearSolve` ident (rows m) | 364 | inv m | square m = m `linearSolve` ident (rows m) |
364 | | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" | 365 | | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" |
365 | 366 | ||
366 | -- | Pseudoinverse of a general matrix. | 367 | |
368 | -- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave). | ||
367 | pinv :: Field t => Matrix t -> Matrix t | 369 | pinv :: Field t => Matrix t -> Matrix t |
368 | pinv m = linearSolveSVD m (ident (rows m)) | 370 | pinv = pinvTol 1 |
371 | |||
372 | {- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value. | ||
373 | |||
374 | @\> let m = 'fromLists' [[1,0, 0] | ||
375 | ,[0,1, 0] | ||
376 | ,[0,0,1e-10]] | ||
377 | \ -- | ||
378 | \> 'pinv' m | ||
379 | 1. 0. 0. | ||
380 | 0. 1. 0. | ||
381 | 0. 0. 10000000000. | ||
382 | \ -- | ||
383 | \> pinvTol 1E8 m | ||
384 | 1. 0. 0. | ||
385 | 0. 1. 0. | ||
386 | 0. 0. 1.@ | ||
387 | |||
388 | -} | ||
389 | |||
390 | pinvTol :: Field t => Double -> Matrix t -> Matrix t | ||
391 | pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where | ||
392 | (u,s,v) = thinSVD m | ||
393 | sl@(g:_) = toList s | ||
394 | s' = real . fromList . map rec $ sl | ||
395 | rec x = if x <= g*tol then x else 1/x | ||
396 | tol = (fromIntegral (max r c) * g * t * eps) | ||
397 | r = rows m | ||
398 | c = cols m | ||
399 | d = dim s | ||
400 | u' = takeColumns d u | ||
401 | v' = takeColumns d v | ||
402 | |||
369 | 403 | ||
370 | -- | Numeric rank of a matrix from the SVD decomposition. | 404 | -- | Numeric rank of a matrix from the SVD decomposition. |
371 | rankSVD :: Element t | 405 | rankSVD :: Element t |
@@ -439,39 +473,6 @@ orth m = take r $ toColumns u | |||
439 | 473 | ||
440 | ------------------------------------------------------------------------ | 474 | ------------------------------------------------------------------------ |
441 | 475 | ||
442 | {- Pseudoinverse of a real matrix with the desired tolerance, expressed as a | ||
443 | multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). | ||
444 | |||
445 | @\> let m = 'fromLists' [[1,0, 0] | ||
446 | ,[0,1, 0] | ||
447 | ,[0,0,1e-10]] | ||
448 | \ -- | ||
449 | \> 'pinv' m | ||
450 | 1. 0. 0. | ||
451 | 0. 1. 0. | ||
452 | 0. 0. 10000000000. | ||
453 | \ -- | ||
454 | \> pinvTol 1E8 m | ||
455 | 1. 0. 0. | ||
456 | 0. 1. 0. | ||
457 | 0. 0. 1.@ | ||
458 | |||
459 | -} | ||
460 | --pinvTol :: Double -> Matrix Double -> Matrix Double | ||
461 | pinvTol t m = v' `mXm` diag s' `mXm` trans u' where | ||
462 | (u,s,v) = thinSVDRd m | ||
463 | sl@(g:_) = toList s | ||
464 | s' = fromList . map rec $ sl | ||
465 | rec x = if x < g*tol then 1 else 1/x | ||
466 | tol = (fromIntegral (max r c) * g * t * eps) | ||
467 | r = rows m | ||
468 | c = cols m | ||
469 | d = dim s | ||
470 | u' = takeColumns d u | ||
471 | v' = takeColumns d v | ||
472 | |||
473 | --------------------------------------------------------------------- | ||
474 | |||
475 | -- many thanks, quickcheck! | 476 | -- many thanks, quickcheck! |
476 | 477 | ||
477 | haussholder :: (Field a) => a -> Vector a -> Matrix a | 478 | haussholder :: (Field a) => a -> Vector a -> Matrix a |
@@ -545,7 +546,7 @@ diagonalize m = if rank v == n | |||
545 | matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | 546 | matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) |
546 | matFunc f m = case diagonalize m of | 547 | matFunc f m = case diagonalize m of |
547 | Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v | 548 | Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v |
548 | Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" | 549 | Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" |
549 | 550 | ||
550 | -------------------------------------------------------------- | 551 | -------------------------------------------------------------- |
551 | 552 | ||
@@ -556,6 +557,7 @@ golubeps p q = a * fromIntegral b / fromIntegral c where | |||
556 | c = fact (p+q) * fact (p+q+1) | 557 | c = fact (p+q) * fact (p+q+1) |
557 | fact n = product [1..n] | 558 | fact n = product [1..n] |
558 | 559 | ||
560 | epslist :: [(Int,Double)] | ||
559 | epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] | 561 | epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] |
560 | 562 | ||
561 | geps delta = head [ k | (k,g) <- epslist, g<delta] | 563 | geps delta = head [ k | (k,g) <- epslist, g<delta] |
@@ -567,11 +569,7 @@ geps delta = head [ k | (k,g) <- epslist, g<delta] | |||
567 | expm :: Field t => Matrix t -> Matrix t | 569 | expm :: Field t => Matrix t -> Matrix t |
568 | expm = expGolub | 570 | expm = expGolub |
569 | 571 | ||
570 | expGolub :: ( Fractional t, Element t, Field t | 572 | expGolub :: Field t => Matrix t -> Matrix t |
571 | , Normed Matrix t | ||
572 | , RealFrac (RealOf t) | ||
573 | , Floating (RealOf t) | ||
574 | ) => Matrix t -> Matrix t | ||
575 | expGolub m = iterate msq f !! j | 573 | expGolub m = iterate msq f !! j |
576 | where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m | 574 | where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m |
577 | a = m */ fromIntegral ((2::Int)^j) | 575 | a = m */ fromIntegral ((2::Int)^j) |