summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs90
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.
333invlndet :: (Floating t, Field t) 334invlndet :: 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))
336invlndet m | square m = (im,(ladm,sdm)) 337invlndet 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
363inv m | square m = m `linearSolve` ident (rows m) 364inv 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).
367pinv :: Field t => Matrix t -> Matrix t 369pinv :: Field t => Matrix t -> Matrix t
368pinv m = linearSolveSVD m (ident (rows m)) 370pinv = 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
3791. 0. 0.
3800. 1. 0.
3810. 0. 10000000000.
382\ --
383\> pinvTol 1E8 m
3841. 0. 0.
3850. 1. 0.
3860. 0. 1.@
387
388-}
389
390pinvTol :: Field t => Double -> Matrix t -> Matrix t
391pinvTol 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.
371rankSVD :: Element t 405rankSVD :: 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
443multiplicative 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
4501. 0. 0.
4510. 1. 0.
4520. 0. 10000000000.
453\ --
454\> pinvTol 1E8 m
4551. 0. 0.
4560. 1. 0.
4570. 0. 1.@
458
459-}
460--pinvTol :: Double -> Matrix Double -> Matrix Double
461pinvTol 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
477haussholder :: (Field a) => a -> Vector a -> Matrix a 478haussholder :: (Field a) => a -> Vector a -> Matrix a
@@ -545,7 +546,7 @@ diagonalize m = if rank v == n
545matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 546matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
546matFunc f m = case diagonalize m of 547matFunc 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
560epslist :: [(Int,Double)]
559epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] 561epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
560 562
561geps delta = head [ k | (k,g) <- epslist, g<delta] 563geps delta = head [ k | (k,g) <- epslist, g<delta]
@@ -567,11 +569,7 @@ geps delta = head [ k | (k,g) <- epslist, g<delta]
567expm :: Field t => Matrix t -> Matrix t 569expm :: Field t => Matrix t -> Matrix t
568expm = expGolub 570expm = expGolub
569 571
570expGolub :: ( Fractional t, Element t, Field t 572expGolub :: Field t => Matrix t -> Matrix t
571 , Normed Matrix t
572 , RealFrac (RealOf t)
573 , Floating (RealOf t)
574 ) => Matrix t -> Matrix t
575expGolub m = iterate msq f !! j 573expGolub 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)