matrix 0.2.3.0 → 0.2.4.0
raw patch · 2 files changed
+115/−10 lines, 2 files
Files
- Data/Matrix.hs +114/−9
- matrix.cabal +1/−1
Data/Matrix.hs view
@@ -47,6 +47,8 @@ , switchCols -- * Decompositions , luDecomp+ , luDecomp'+ , cholDecomp -- * Properties , trace , diagProd -- ** Determinants@@ -63,6 +65,7 @@ -- Data import Control.Monad.Primitive (PrimMonad, PrimState) import Data.List (maximumBy)+import Data.Ord (comparing) import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV @@ -719,8 +722,8 @@ luDecomp :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,a) luDecomp a = recLUDecomp a i i 1 1 n where+ i = (identity $ nrows a) n = min (nrows a) (ncols a)- i = identity $ nrows a recLUDecomp :: (Ord a, Fractional a) => Matrix a -- ^ U@@ -731,17 +734,17 @@ -> Int -- ^ Total rows -> (Matrix a,Matrix a,Matrix a,a) recLUDecomp u l p d k n =- if k == n then (u,l,p,d)- else recLUDecomp u'' l'' p' d' (k+1) n+ if k > n then (u,l,p,d)+ else recLUDecomp u'' l'' p' d' (k+1) n where -- Pivot strategy: maximum value in absolute value below the current row. i = maximumBy (\x y -> compare (abs $ u ! (x,k)) (abs $ u ! (y,k))) [ k .. n ] -- Switching to place pivot in current row. u' = switchRows k i u- l' = M n n $+ l' = M (nrows l) (ncols l) $ V.modify (\mv -> mapM_ (\j -> do- msetElem (l ! (k,j)) n (i,j) mv- msetElem (l ! (i,j)) n (k,j) mv+ msetElem (l ! (k,j)) (ncols l) (i,j) mv+ msetElem (l ! (i,j)) (ncols l) (k,j) mv ) [1 .. k-1] ) $ mvect l p' = switchRows k i p -- Permutation determinant@@ -750,10 +753,112 @@ (u'',l'') = go u' l' (k+1) ukk = u' ! (k,k) go u_ l_ j =- if j > n then (u_,l_)- else let x = (u_ ! (j,k)) / ukk- in go (combineRows j (-x) k u_) (setElem x (j,k) l_) (j+1)+ if j > nrows u_+ then (u_,l_)+ else let x = (u_ ! (j,k)) / ukk+ in go (combineRows j (-x) k u_) (setElem x (j,k) l_) (j+1) +-- | Matrix LU decomposition with /complete pivoting/.+-- The result for a matrix /M/ is given in the format /(U,L,P,Q,d,e)/ where:+--+-- * /U/ is an upper triangular matrix.+--+-- * /L/ is an /unit/ lower triangular matrix.+--+-- * /P,Q/ is a permutation matrix.+--+-- * /d,e/ is the determinant of /P,Q/.+--+-- * /PMQ = LU/.+--+-- These properties are only guaranteed when the input matrix is invertible.+-- An additional property matches thanks to the strategy followed for pivoting:+--+-- * /L_(i,j)/ <= 1, for all /i,j/.+--+-- This follows from the maximal property of the selected pivots, which also+-- leads to a better numerical stability of the algorithm.+--+-- Example:+--+-- > ( 1 0 ) ( 2 1 ) ( 1 0 0 ) ( 0 0 1 )+-- > ( 0 2 ) ( 0 2 ) ( 0 1 0 ) ( 0 1 0 ) ( 1 0 )+-- > luDecomp' ( 2 1 ) = ( ( 0 0 ) , ( 1/2 -1/4 1 ) , ( 1 0 0 ) , ( 0 1 ) , -1 , 1 )+luDecomp' :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)+luDecomp' a = recLUDecomp' a i i (identity $ ncols a) 1 1 1 n+ where+ i = identity $ nrows a+ n = min (nrows a) (ncols a)++recLUDecomp' :: (Ord a, Fractional a)+ => Matrix a -- ^ U+ -> Matrix a -- ^ L+ -> Matrix a -- ^ P+ -> Matrix a -- ^ Q+ -> a -- ^ d+ -> a -- ^ e+ -> Int -- ^ Current row+ -> Int -- ^ Total rows+ -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)+recLUDecomp' u l p q d e k n =+ if k > n || u'' ! (k, k) == 0+ then (u,l,p,q,d,e)+ else recLUDecomp' u'' l'' p' q' d' e' (k+1) n+ where+ -- Pivot strategy: maximum value in absolute value below the current row & col.+ (i, j) = maximumBy (comparing (\(i0, j0) -> abs $ u ! (i0,j0)))+ [ (i0, j0) | i0 <- [k .. nrows u], j0 <- [k .. ncols u] ]+ -- Switching to place pivot in current row.+ u' = switchCols k j $ switchRows k i u+ l'0 = M (nrows l) (ncols l) $+ V.modify (\mv -> forM_ [1..k-1] $ \ h -> do+ msetElem (l ! (k,h)) (ncols l) (i,h) mv+ msetElem (l ! (i,h)) (ncols l) (k,h) mv+ )+ $ mvect l+ l' = M (nrows l) (ncols l) $+ V.modify (\mv -> forM_ [1..k-1] $ \h -> do+ msetElem (l'0 ! (h,k)) (ncols l) (h,i) mv+ msetElem (l'0 ! (h,i)) (ncols l) (h,k) mv+ )+ $ mvect l'0+ p' = switchRows k i p+ q' = switchCols k j q+ -- Permutation determinant+ d' = if i == k then d else negate d+ e' = if j == k then e else negate e+ -- Cancel elements below the pivot.+ (u'',l'') = go u' l' (k+1)+ ukk = u' ! (k,k)+ go u_ l_ h =+ if h > nrows u_+ then (u_,l_)+ else let x = (u_ ! (h,k)) / ukk+ in go (combineRows h (-x) k u_) (setElem x (h,k) l_) (h+1)++-- CHOLESKY DECOMPOSITION++-- | Simple Cholesky decomposition of a symmetric, positive definite matrix.+-- The result for a matrix /M/ is a lower triangular matrix /L/ such that:+--+-- * /M = LL^T/.+--+-- Example:+--+-- > ( 2 -1 0 ) ( 1.41 0 0 )+-- > ( -1 2 -1 ) ( -0.70 1.22 0 )+-- > cholDecomp ( 0 -1 2 ) = ( 0.00 -0.81 1.15 )+cholDecomp :: (Floating a) => Matrix a -> Matrix a+cholDecomp a+ | (nrows a == 1) && (ncols a == 1) = fmap sqrt a+ | otherwise = joinBlocks (l11, l12, l21, l22) where+ (a11, a12, a21, a22) = splitBlocks 1 1 a+ l11' = sqrt (a11 ! (1,1))+ l11 = fromList 1 1 [l11']+ l12 = zero (nrows a12) (ncols a12)+ l21 = scaleMatrix (1/l11') a21+ a22' = a22 - multStd l21 (transpose l21)+ l22 = cholDecomp a22' ------------------------------------------------------- ------------------------------------------------------- ---- PROPERTIES
matrix.cabal view
@@ -1,5 +1,5 @@ Name: matrix -Version: 0.2.3.0 +Version: 0.2.4.0 Author: Daniel Díaz Category: Math Build-type: Simple