From a273fdb74b04db6d57d5c9b15e676d83357e71fd Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 16 Jul 2015 20:16:59 +0200 Subject: Her, LU, LDL, Linear, Additive --- packages/base/src/Internal/Algorithms.hs | 138 ++++++++++++++++----- packages/base/src/Internal/CG.hs | 20 +-- packages/base/src/Internal/Modular.hs | 10 +- packages/base/src/Internal/Numeric.hs | 50 +++++--- packages/base/src/Internal/Util.hs | 14 ++- packages/base/src/Numeric/LinearAlgebra.hs | 47 +++---- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 3 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 12 +- packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 48 +++---- .../src/Numeric/LinearAlgebra/Tests/Instances.hs | 18 +-- .../src/Numeric/LinearAlgebra/Tests/Properties.hs | 12 +- 11 files changed, 239 insertions(+), 133 deletions(-) (limited to 'packages') diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index 3d25491..d2f17f4 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs @@ -4,6 +4,12 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} + ----------------------------------------------------------------------------- {- | Module : Internal.Algorithms @@ -32,6 +38,7 @@ import Data.List(foldl1') import qualified Data.Array as A import Internal.ST import Internal.Vectorized(range) +import Control.DeepSeq {- | Generic linear algebra functions for double precision real and complex matrices. @@ -43,6 +50,10 @@ class (Numeric t, Normed Matrix t, Normed Vector t, Floating t, + Linear t Vector, + Linear t Matrix, + Additive (Vector t), + Additive (Matrix t), RealOf t ~ Double) => Field t where svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) @@ -306,25 +317,38 @@ leftSV m | vertical m = let (u,s,_) = svd m in (u,s) -------------------------------------------------------------- +-- | LU decomposition of a matrix in a compact format. +data LU t = LU (Matrix t) [Int] deriving Show + +instance (NFData t, Numeric t) => NFData (LU t) + where + rnf (LU m _) = rnf m + -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. -luPacked :: Field t => Matrix t -> (Matrix t, [Int]) -luPacked = {-# SCC "luPacked" #-} luPacked' +luPacked :: Field t => Matrix t -> LU t +luPacked x = {-# SCC "luPacked" #-} LU m p + where + (m,p) = luPacked' x -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'. -luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t -luSolve = {-# SCC "luSolve" #-} luSolve' +luSolve :: Field t => LU t -> Matrix t -> Matrix t +luSolve (LU m p) = {-# SCC "luSolve" #-} luSolve' (m,p) -- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t linearSolve = {-# SCC "linearSolve" #-} linearSolve' --- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. +-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. mbLinearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t) mbLinearSolve = {-# SCC "linearSolve" #-} mbLinearSolve' -- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'. -cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t +cholSolve + :: Field t + => Matrix t -- ^ Cholesky decomposition of the coefficient matrix + -> Matrix t -- ^ right hand sides + -> Matrix t -- ^ solution cholSolve = {-# SCC "cholSolve" #-} cholSolve' -- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value. @@ -338,20 +362,28 @@ linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS' -------------------------------------------------------------------------------- +-- | LDL decomposition of a complex Hermitian or real symmetric matrix in a compact format. +data LDL t = LDL (Matrix t) [Int] deriving Show + +instance (NFData t, Numeric t) => NFData (LDL t) + where + rnf (LDL m _) = rnf m + -- | Similar to 'ldlPacked', without checking that the input matrix is hermitian or symmetric. It works with the lower triangular part. -ldlPackedSH :: Field t => Matrix t -> (Matrix t, [Int]) -ldlPackedSH = {-# SCC "ldlPacked" #-} ldlPacked' +ldlPackedSH :: Field t => Matrix t -> LDL t +ldlPackedSH x = {-# SCC "ldlPacked" #-} LDL m p + where + (m,p) = ldlPacked' x -- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'. -ldlPacked :: Field t => Matrix t -> (Matrix t, [Int]) -ldlPacked m - | exactHermitian m = {-# SCC "ldlPacked" #-} ldlPackedSH m - | otherwise = error "ldlPacked requires complex Hermitian or real symmetrix matrix" - +ldlPacked :: Field t => Her t -> LDL t +ldlPacked (Her m) = ldlPackedSH m --- | Solution of a linear system (for several right hand sides) from the precomputed LDL factorization obtained by 'ldlPacked'. -ldlSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t -ldlSolve = {-# SCC "ldlSolve" #-} ldlSolve' +-- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'. +-- +-- Note: this can be slower than the general solver based on the LU decomposition. +ldlSolve :: Field t => LDL t -> Matrix t -> Matrix t +ldlSolve (LDL m p) = {-# SCC "ldlSolve" #-} ldlSolve' (m,p) -------------------------------------------------------------- @@ -429,14 +461,12 @@ fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575] 3.000 5.000 6.000 -} -eigSH :: Field t => Matrix t -> (Vector Double, Matrix t) -eigSH m | exactHermitian m = eigSH' m - | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" +eigSH :: Field t => Her t -> (Vector Double, Matrix t) +eigSH (Her m) = eigSH' m -- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix. -eigenvaluesSH :: Field t => Matrix t -> Vector Double -eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m - | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix" +eigenvaluesSH :: Field t => Her t -> Vector Double +eigenvaluesSH (Her m) = eigenvaluesSH' m -------------------------------------------------------------- @@ -490,14 +520,18 @@ mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH' -- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. cholSH :: Field t => Matrix t -> Matrix t -cholSH = {-# SCC "cholSH" #-} cholSH' +cholSH = cholSH' -- | Cholesky factorization of a positive definite hermitian or symmetric matrix. -- -- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@. -chol :: Field t => Matrix t -> Matrix t -chol m | exactHermitian m = cholSH m - | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" +chol :: Field t => Her t -> Matrix t +chol (Her m) = {-# SCC "chol" #-} cholSH' m + +-- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. +mbChol :: Field t => Her t -> Maybe (Matrix t) +mbChol (Her m) = {-# SCC "mbChol" #-} mbCholSH' m + -- | Joint computation of inverse and logarithm of determinant of a square matrix. @@ -507,7 +541,7 @@ invlndet :: Field t invlndet m | square m = (im,(ladm,sdm)) | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" where - lp@(lup,perm) = luPacked m + lp@(LU lup perm) = luPacked m s = signlp (rows m) perm dg = toList $ takeDiag $ lup ladm = sum $ map (log.abs) dg @@ -519,8 +553,9 @@ invlndet m | square m = (im,(ladm,sdm)) det :: Field t => Matrix t -> t det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix" - where (lup,perm) = luPacked m - s = signlp (rows m) perm + where + LU lup perm = luPacked m + s = signlp (rows m) perm -- | Explicit LU factorization of a general matrix. -- @@ -720,7 +755,7 @@ diagonalize m = if rank v == n else Nothing where n = rows m (l,v) = if exactHermitian m - then let (l',v') = eigSH m in (real l', v') + then let (l',v') = eigSH (trustSym m) in (real l', v') else eig m -- | Generic matrix functions for diagonalizable matrices. For instance: @@ -835,8 +870,9 @@ fixPerm' s = res $ mutable f s0 triang r c h v = (r>=h then v else 1 - v -luFact (l_u,perm) | r <= c = (l ,u ,p, s) - | otherwise = (l',u',p, s) +luFact (LU l_u perm) + | r <= c = (l ,u ,p, s) + | otherwise = (l',u',p, s) where r = rows l_u c = cols l_u @@ -929,7 +965,13 @@ relativeError norm a b = r ---------------------------------------------------------------------- -- | Generalized symmetric positive definite eigensystem Av = lBv, --- for A and B symmetric, B positive definite (conditions not checked). +-- for A and B symmetric, B positive definite. +geigSH :: Field t + => Her t -- ^ A + -> Her t -- ^ B + -> (Vector Double, Matrix t) +geigSH (Her a) (Her b) = geigSH' a b + geigSH' :: Field t => Matrix t -- ^ A -> Matrix t -- ^ B @@ -943,3 +985,33 @@ geigSH' a b = (l,v') v' = iu <> v (<>) = mXm +-------------------------------------------------------------------------------- + +-- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric. +-- +-- It can be created using 'sym', 'xTx', or 'trustSym', and the matrix can be extracted using 'her'. +data Her t = Her (Matrix t) deriving Show + +-- | Extract the general matrix from a 'Her' structure, forgetting its symmetric or Hermitian property. +her :: Her t -> Matrix t +her (Her x) = x + +-- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@). +sym :: Field t => Matrix t -> Her t +sym x = Her (scale 0.5 (tr x `add` x)) + +-- | Compute the contraction @tr x <> x@ of a general matrix. +xTx :: Numeric t => Matrix t -> Her t +xTx x = Her (tr x `mXm` x) + +instance Field t => Linear t Her where + scale x (Her m) = Her (scale x m) + +instance Field t => Additive (Her t) where + add (Her a) (Her b) = Her (a `add` b) + +-- | At your own risk, declare that a matrix is complex Hermitian or real symmetric +-- for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used. +trustSym :: Matrix t -> Her t +trustSym x = (Her x) + diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs index f0142cd..cc10ad8 100644 --- a/packages/base/src/Internal/CG.hs +++ b/packages/base/src/Internal/CG.hs @@ -32,11 +32,11 @@ v /// b = debugMat b 2 asRow v type V = Vector R data CGState = CGState - { cgp :: V -- ^ conjugate gradient - , cgr :: V -- ^ residual - , cgr2 :: R -- ^ squared norm of residual - , cgx :: V -- ^ current solution - , cgdx :: R -- ^ normalized size of correction + { cgp :: Vector R -- ^ conjugate gradient + , cgr :: Vector R -- ^ residual + , cgr2 :: R -- ^ squared norm of residual + , cgx :: Vector R -- ^ current solution + , cgdx :: R -- ^ normalized size of correction } cg :: Bool -> (V -> V) -> (V -> V) -> CGState -> CGState @@ -89,23 +89,25 @@ takeUntil q xs = a++ take 1 b where (a,b) = break q xs +-- | Solve a sparse linear system using the conjugate gradient method with default parameters. cgSolve :: Bool -- ^ is symmetric -> GMatrix -- ^ coefficient matrix - -> Vector Double -- ^ right-hand side - -> Vector Double -- ^ solution + -> Vector R -- ^ right-hand side + -> Vector R -- ^ solution cgSolve sym a b = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0 where n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double)) +-- | Solve a sparse linear system using the conjugate gradient method with default parameters. cgSolve' :: Bool -- ^ symmetric -> R -- ^ relative tolerance for the residual (e.g. 1E-4) -> R -- ^ relative tolerance for δx (e.g. 1E-3) -> Int -- ^ maximum number of iterations -> GMatrix -- ^ coefficient matrix - -> V -- ^ initial solution - -> V -- ^ right-hand side + -> Vector R -- ^ initial solution + -> Vector R -- ^ right-hand side -> [CGState] -- ^ solution cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 64ed2bb..a3421a8 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -33,7 +33,7 @@ import Internal.Element import Internal.Container import Internal.Vectorized (prodI,sumI,prodL,sumL) import Internal.LAPACK (multiplyI, multiplyL) -import Internal.Algorithms(luFact) +import Internal.Algorithms(luFact,LU(..)) import Internal.Util(Normed(..),Indexable(..), gaussElim, gaussElim_1, gaussElim_2, luST, luSolve', luPacked', magnit, invershur) @@ -169,7 +169,7 @@ instance forall m . KnownNat m => Container Vector (Mod m I) size' = dim scale' s x = vmod (scale (unMod s) (f2i x)) addConstant c x = vmod (addConstant (unMod c) (f2i x)) - add a b = vmod (add (f2i a) (f2i b)) + add' a b = vmod (add' (f2i a) (f2i b)) sub a b = vmod (sub (f2i a) (f2i b)) mul a b = vmod (mul (f2i a) (f2i b)) equal u v = equal (f2i u) (f2i v) @@ -209,7 +209,7 @@ instance forall m . KnownNat m => Container Vector (Mod m Z) size' = dim scale' s x = vmod (scale (unMod s) (f2i x)) addConstant c x = vmod (addConstant (unMod c) (f2i x)) - add a b = vmod (add (f2i a) (f2i b)) + add' a b = vmod (add' (f2i a) (f2i b)) sub a b = vmod (sub (f2i a) (f2i b)) mul a b = vmod (mul (f2i a) (f2i b)) equal u v = equal (f2i u) (f2i v) @@ -371,7 +371,9 @@ test = (ok, info) checkLU okf t = norm_Inf $ flatten (l <> u <> p - t) where - (l,u,p,_ :: Int) = luFact $ mutable (luST okf) t + (l,u,p,_ :: Int) = luFact (LU x' p') + where + (x',p') = mutable (luST okf) t checkSolve aa = norm_Inf $ flatten (aa <> x - bb) where diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs index a8ae2bb..e8c7440 100644 --- a/packages/base/src/Internal/Numeric.hs +++ b/packages/base/src/Internal/Numeric.hs @@ -49,7 +49,7 @@ class Element e => Container c e scalar' :: e -> c e scale' :: e -> c e -> c e addConstant :: e -> c e -> c e - add :: c e -> c e -> c e + add' :: c e -> c e -> c e sub :: c e -> c e -> c e -- | element by element multiplication mul :: c e -> c e -> c e @@ -100,7 +100,7 @@ instance Container Vector I size' = dim scale' = vectorMapValI Scale addConstant = vectorMapValI AddConstant - add = vectorZipI Add + add' = vectorZipI Add sub = vectorZipI Sub mul = vectorZipI Mul equal u v = dim u == dim v && maxElement' (vectorMapI Abs (sub u v)) == 0 @@ -139,7 +139,7 @@ instance Container Vector Z size' = dim scale' = vectorMapValL Scale addConstant = vectorMapValL AddConstant - add = vectorZipL Add + add' = vectorZipL Add sub = vectorZipL Sub mul = vectorZipL Mul equal u v = dim u == dim v && maxElement' (vectorMapL Abs (sub u v)) == 0 @@ -179,7 +179,7 @@ instance Container Vector Float size' = dim scale' = vectorMapValF Scale addConstant = vectorMapValF AddConstant - add = vectorZipF Add + add' = vectorZipF Add sub = vectorZipF Sub mul = vectorZipF Mul equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 @@ -216,7 +216,7 @@ instance Container Vector Double size' = dim scale' = vectorMapValR Scale addConstant = vectorMapValR AddConstant - add = vectorZipR Add + add' = vectorZipR Add sub = vectorZipR Sub mul = vectorZipR Mul equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 @@ -253,7 +253,7 @@ instance Container Vector (Complex Double) size' = dim scale' = vectorMapValC Scale addConstant = vectorMapValC AddConstant - add = vectorZipC Add + add' = vectorZipC Add sub = vectorZipC Sub mul = vectorZipC Mul equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 @@ -289,7 +289,7 @@ instance Container Vector (Complex Float) size' = dim scale' = vectorMapValQ Scale addConstant = vectorMapValQ AddConstant - add = vectorZipQ Add + add' = vectorZipQ Add sub = vectorZipQ Sub mul = vectorZipQ Mul equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 @@ -327,7 +327,7 @@ instance (Num a, Element a, Container Vector a) => Container Matrix a size' = size scale' x = liftMatrix (scale' x) addConstant x = liftMatrix (addConstant x) - add = liftMatrix2 add + add' = liftMatrix2 add' sub = liftMatrix2 sub mul = liftMatrix2 mul equal a b = cols a == cols b && flatten a `equal` flatten b @@ -387,9 +387,6 @@ scalar = scalar' conj :: Container c e => c e -> c e conj = conj' --- | multiplication by scalar -scale :: Container c e => e -> c e -> c e -scale = scale' arctan2 :: (Fractional e, Container c e) => c e -> c e -> c e arctan2 = arctan2' @@ -581,6 +578,10 @@ class ( Container Vector t , Konst t (Int,Int) Matrix , CTrans t , Product t + , Additive (Vector t) + , Additive (Matrix t) + , Linear t Vector + , Linear t Matrix ) => Numeric t instance Numeric Double @@ -912,11 +913,30 @@ instance (CTrans t, Container Vector t) => Transposable (Matrix t) (Matrix t) tr = ctrans tr' = trans -class Linear t v +class Additive c where - scalarL :: t -> v - addL :: v -> v -> v - scaleL :: t -> v -> v + add :: c -> c -> c + +class Linear t c + where + scale :: t -> c t -> c t + + +instance Container Vector t => Linear t Vector + where + scale = scale' + +instance Container Matrix t => Linear t Matrix + where + scale = scale' + +instance Container Vector t => Additive (Vector t) + where + add = add' + +instance Container Matrix t => Additive (Matrix t) + where + add = add' class Testable t diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 4123e6c..36b7855 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -458,12 +458,12 @@ rowOuters a b = a' * b' -------------------------------------------------------------------------------- -- | solution of overconstrained homogeneous linear system -null1 :: Matrix Double -> Vector Double +null1 :: Matrix R -> Vector R null1 = last . toColumns . snd . rightSV -- | solution of overconstrained homogeneous symmetric linear system -null1sym :: Matrix Double -> Vector Double -null1sym = last . toColumns . snd . eigSH' +null1sym :: Her R -> Vector R +null1sym = last . toColumns . snd . eigSH -------------------------------------------------------------------------------- @@ -712,7 +712,9 @@ luST ok (r,_) x = do , 0, 0, 0, 0, 1 ] -} -luPacked' x = mutable (luST (magnit 0)) x +luPacked' x = LU m p + where + (m,p) = mutable (luST (magnit 0)) x -------------------------------------------------------------------------------- @@ -782,7 +784,7 @@ forwSust' lup rhs = foldl' f (rhs?[]) ls (b - l<>x) -luSolve'' (lup,p) b = backSust' lup (forwSust' lup pb) +luSolve'' (LU lup p) b = backSust' lup (forwSust' lup pb) where pb = b ?? (Pos (fixPerm' p), All) @@ -827,7 +829,7 @@ backSust lup rhs = fst $ mutable f rhs , 7, 10, 6 ] -} -luSolve' (lup,p) b = backSust lup (forwSust lup pb) +luSolve' (LU lup p) b = backSust lup (forwSust lup pb) where pb = b ?? (Pos (fixPerm' p), All) diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index 9a924e0..7be2600 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -53,11 +53,11 @@ module Numeric.LinearAlgebra ( -- -- * Products - -- ** dot + -- ** Dot dot, (<.>), - -- ** matrix-vector + -- ** Matrix-vector (#>), (<#), (!#>), - -- ** matrix-matrix + -- ** Matrix-matrix (<>), -- | The matrix product is also implemented in the "Data.Monoid" instance, where -- single-element matrices (created from numeric literals or using 'scalar') @@ -73,20 +73,25 @@ module Numeric.LinearAlgebra ( -- 'mconcat' uses 'optimiseMult' to get the optimal association order. - -- ** other + -- ** Other outer, kronecker, cross, - scale, + scale, add, sumElements, prodElements, -- * Linear systems + -- ** General (<\>), - linearSolve, linearSolveLS, linearSolveSVD, - luSolve, - luSolve', + -- ** Determined + linearSolve, + luSolve, luPacked, + luSolve', luPacked', + -- ** Symmetric indefinite + ldlSolve, ldlPacked, + -- ** Positive definite cholSolve, - ldlSolve, + -- ** Sparse cgSolve, cgSolve', @@ -113,21 +118,18 @@ module Numeric.LinearAlgebra ( leftSV, rightSV, -- * Eigendecomposition - eig, eigSH, eigSH', - eigenvalues, eigenvaluesSH, eigenvaluesSH', - geigSH', + eig, eigSH, + eigenvalues, eigenvaluesSH, + geigSH, -- * QR qr, rq, qrRaw, qrgr, -- * Cholesky - chol, cholSH, mbCholSH, + chol, mbChol, -- * LU - lu, luPacked, luPacked', luFact, - - -- * LDL - ldlPacked, ldlPackedSH, + lu, luFact, -- * Hessenberg hess, @@ -150,14 +152,16 @@ module Numeric.LinearAlgebra ( -- * Misc meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, - iC, + iC, sym, xTx, trustSym, her, -- * Auxiliary classes - Element, Container, Product, Numeric, LSDiv, + Element, Container, Product, Numeric, LSDiv, Her, Complexable, RealElement, RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf, - Field, + Field, Linear(), Additive(), Transposable, + LU(..), + LDL(..), CGState(..), Testable(..) ) where @@ -169,7 +173,7 @@ import Numeric.Vector() import Internal.Matrix import Internal.Container hiding ((<>)) import Internal.Numeric hiding (mul) -import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve') +import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve',ldlPacked') import qualified Internal.Algorithms as A import Internal.Util import Internal.Random @@ -246,4 +250,3 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m) -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'. orth m = orthSVD (Left (1*eps)) m (leftSV m) - diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index bac1c0c..5ce529c 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -13,11 +13,12 @@ compatibility with previous version, to be removed module Numeric.LinearAlgebra.HMatrix ( module Numeric.LinearAlgebra, - (¦),(——),ℝ,ℂ,(<·>),app,mul + (¦),(——),ℝ,ℂ,(<·>),app,mul, cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH' ) where import Numeric.LinearAlgebra import Internal.Util +import Internal.Algorithms(cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH') infixr 8 <·> (<·>) :: Numeric t => Vector t -> Vector t -> t diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index 0dab0e6..ded69fa 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -63,9 +63,9 @@ import GHC.TypeLits import Numeric.LinearAlgebra hiding ( (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||), row,col,vector,matrix,linspace,toRows,toColumns, - (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', - eigenvalues,eigenvaluesSH,eigenvaluesSH',build, - qr,size,dot,chol,range,R,C) + (<\>),fromList,takeDiag,svd,eig,eigSH, + eigenvalues,eigenvaluesSH,build, + qr,size,dot,chol,range,R,C,Her,her,sym) import qualified Numeric.LinearAlgebra as LA import Data.Proxy(Proxy) import Internal.Static @@ -292,10 +292,10 @@ her m = Her $ (m + LA.tr m)/2 instance KnownNat n => Eigen (Sym n) (R n) (L n n) where - eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH' $ m + eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH . LA.trustSym $ m eigensystem (Sym (extract -> m)) = (mkR l, mkL v) where - (l,v) = LA.eigSH' m + (l,v) = LA.eigSH . LA.trustSym $ m instance KnownNat n => Eigen (Sq n) (C n) (M n n) where @@ -305,7 +305,7 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n) (l,v) = LA.eig m chol :: KnownNat n => Sym n -> Sq n -chol (extract . unSym -> m) = mkL $ LA.cholSH m +chol (extract . unSym -> m) = mkL $ LA.chol $ LA.trustSym m -------------------------------------------------------------------------------- diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 2ff1580..30480d7 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs @@ -127,8 +127,8 @@ expmTest2 = expm nd2 :~15~: (2><2) mbCholTest = utest "mbCholTest" (ok1 && ok2) where m1 = (2><2) [2,5,5,8 :: Double] m2 = (2><2) [3,5,5,9 :: Complex Double] - ok1 = mbCholSH m1 == Nothing - ok2 = mbCholSH m2 == Just (chol m2) + ok1 = mbChol (trustSym m1) == Nothing + ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2) --------------------------------------------------------------------- @@ -403,8 +403,8 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3 -------------------------------------------------------------------------------- sliceTest = utest "slice test" $ and - [ testSlice chol (gen 5 :: Matrix R) - , testSlice chol (gen 5 :: Matrix C) + [ testSlice (chol . trustSym) (gen 5 :: Matrix R) + , testSlice (chol . trustSym) (gen 5 :: Matrix C) , testSlice qr (rec :: Matrix R) , testSlice qr (rec :: Matrix C) , testSlice hess (agen 5 :: Matrix R) @@ -420,12 +420,12 @@ sliceTest = utest "slice test" $ and , testSlice eig (agen 5 :: Matrix R) , testSlice eig (agen 5 :: Matrix C) - , testSlice eigSH (gen 5 :: Matrix R) - , testSlice eigSH (gen 5 :: Matrix C) + , testSlice (eigSH . trustSym) (gen 5 :: Matrix R) + , testSlice (eigSH . trustSym) (gen 5 :: Matrix C) , testSlice eigenvalues (agen 5 :: Matrix R) , testSlice eigenvalues (agen 5 :: Matrix C) - , testSlice eigenvaluesSH (gen 5 :: Matrix R) - , testSlice eigenvaluesSH (gen 5 :: Matrix C) + , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R) + , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C) , testSlice svd (rec :: Matrix R) , testSlice thinSVD (rec :: Matrix R) @@ -489,10 +489,10 @@ sliceTest = utest "slice test" $ and , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5) , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5) - , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ gen 5) - , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ gen 5) - , testSlice (cholSolve (chol $ gen 5:: Matrix R)) (agen 5) - , testSlice (cholSolve (chol $ gen 5:: Matrix C)) (agen 5) + , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5) + , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5) + , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5) + , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5) , ok_qrgr (rec :: Matrix R) , ok_qrgr (rec :: Matrix C) @@ -515,8 +515,8 @@ sliceTest = utest "slice test" $ and test_lus m = testSlice f lup where - f x = luSolve (x,p) m - (lup,p) = luPacked m + f x = luSolve (LU x p) m + (LU lup p) = luPacked m gen :: Numeric t => Int -> Matrix t gen n = diagRect 1 (konst 5 n) n n @@ -588,11 +588,11 @@ runTests n = do test (linearSolveProp (luSolve.luPacked) . rSqWC) test (linearSolveProp (luSolve.luPacked) . cSqWC) putStrLn "------ ldlSolve" - test (linearSolveProp (ldlSolve.ldlPacked) . rSymWC) - test (linearSolveProp (ldlSolve.ldlPacked) . cSymWC) + test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC) + test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC) putStrLn "------ cholSolve" - test (linearSolveProp (cholSolve.chol) . rPosDef) - test (linearSolveProp (cholSolve.chol) . cPosDef) + test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef) + test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef) putStrLn "------ luSolveLS" test (linearSolveProp linearSolveLS . rSqWC) test (linearSolveProp linearSolveLS . cSqWC) @@ -865,8 +865,8 @@ eigBench = do let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) s = m + tr m m `seq` s `seq` putStrLn "" - time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) - time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) + time "eigenvalues symmetric 1000x1000" (eigenvaluesSH (trustSym m)) + time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m)) time "eigenvalues general 1000x1000" (eigenvalues m) time "eigenvectors general 1000x1000" (snd $ eig m) @@ -893,12 +893,14 @@ solveBenchN n = do time ("svd solve " ++ show n) (linearSolveSVD a b) time (" ls solve " ++ show n) (linearSolveLS a b) time (" solve " ++ show n) (linearSolve a b) - time ("cholSolve " ++ show n) (cholSolve (chol a) b) +-- time (" LU solve " ++ show n) (luSolve (luPacked a) b) + time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b) + time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b) solveBench = do solveBenchN 500 solveBenchN 1000 - -- solveBenchN 1500 + solveBenchN 1500 -------------------------------- @@ -906,7 +908,7 @@ cholBenchN n = do let x = uniformSample 777 (2*n) (replicate n (-1,1)) a = tr x <> x a `seq` putStr "" - time ("chol " ++ show n) (chol a) + time ("chol " ++ show n) (chol $ trustSym a) cholBench = do putStrLn "" diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 7c54535..4704989 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs @@ -14,7 +14,7 @@ Arbitrary instances for vectors, matrices. module Numeric.LinearAlgebra.Tests.Instances( Sq(..), rSq,cSq, Rot(..), rRot,cRot, - Her(..), rHer,cHer, + rHer,cHer, WC(..), rWC,cWC, SqWC(..), rSqWC, cSqWC, rSymWC, cSymWC, PosDef(..), rPosDef, cPosDef, @@ -81,12 +81,12 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where -- a complex hermitian or real symmetric matrix -newtype (Her a) = Her (Matrix a) deriving Show +--newtype (Her a) = Her (Matrix a) deriving Show instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where arbitrary = do Sq m <- arbitrary let m' = m/2 - return $ Her (m' + tr m') + return $ sym m' class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a @@ -125,9 +125,9 @@ newtype (PosDef a) = PosDef (Matrix a) deriving Show instance (Numeric a, ArbitraryField a, Num (Vector a)) => Arbitrary (PosDef a) where arbitrary = do - Her m <- arbitrary + m <- arbitrary let (_,v) = eigSH m - n = rows m + n = rows (her m) l <- replicateM n (choose (0,100)) let s = diag (fromList l) p = v <> real s <> tr v @@ -161,8 +161,8 @@ fM m = m :: FM zM m = m :: ZM -rHer (Her m) = m :: RM -cHer (Her m) = m :: CM +rHer m = her m :: RM +cHer m = her m :: CM rRot (Rot m) = m :: RM cRot (Rot m) = m :: CM @@ -176,8 +176,8 @@ cWC (WC m) = m :: CM rSqWC (SqWC m) = m :: RM cSqWC (SqWC m) = m :: CM -rSymWC (SqWC m) = m + tr m :: RM -cSymWC (SqWC m) = m + tr m :: CM +rSymWC (SqWC m) = sym m :: Her R +cSymWC (SqWC m) = sym m :: Her C rPosDef (PosDef m) = m :: RM cPosDef (PosDef m) = m :: CM diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs index 207a303..2ac3588 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs @@ -39,7 +39,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( expmDiagProp, multProp1, multProp2, subProp, - linearSolveProp, linearSolveProp2 + linearSolveProp, linearSolvePropH, linearSolveProp2 ) where import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary) @@ -209,11 +209,11 @@ eigProp m = complex m <> v |~| v <> diag s eigSHProp m = m <> v |~| v <> real (diag s) && unitary v && m |~| v <> real (diag s) <> tr v - where (s, v) = eigSH m + where (s, v) = eigSH' m eigProp2 m = fst (eig m) |~| eigenvalues m -eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m +eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m ------------------------------------------------------------------ @@ -246,9 +246,9 @@ schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme where (u,s) = schur m cholProp m = m |~| tr c <> c && upperTriang c - where c = chol m + where c = chol (trustSym m) -exactProp m = chol m == chol (m+0) +exactProp m = chol (trustSym m) == chol (trustSym (m+0)) expmDiagProp m = expm (logm m) :~ 7 ~: complex m where logm = matFunc log @@ -263,6 +263,8 @@ multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a) linearSolveProp f m = f m m |~| ident (rows m) +linearSolvePropH f m = f m (her m) |~| ident (rows (her m)) + linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) where q = min (rows a) (cols a) b = a <> x -- cgit v1.2.3