From 3815bc25f62124063e02af83fe3c907336dc86f5 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 28 Sep 2007 12:25:56 +0000 Subject: algorithms interface reorganized --- lib/LinearAlgebra/Algorithms.hs | 125 ++++++++++++++++++---------------------- lib/LinearAlgebra/Interface.hs | 7 +++ lib/LinearAlgebra/Linear.hs | 43 +++----------- 3 files changed, 69 insertions(+), 106 deletions(-) (limited to 'lib/LinearAlgebra') diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs index 162426f..a67f822 100644 --- a/lib/LinearAlgebra/Algorithms.hs +++ b/lib/LinearAlgebra/Algorithms.hs @@ -9,106 +9,114 @@ Maintainer : Alberto Ruiz (aruiz at um dot es) Stability : provisional Portability : uses ffi -A simple interface to the available matrix computations. We have defined some generic functions -on both real and complex matrices in such a way that higher level algorithms and -testing properties are valid for both base types. +A generic interface for a number of essential functions. Using it some higher level algorithms +and testing properties can be written for both real and complex matrices. -In any case the specific functions for particular base types can also be explicitly +In any case, the specific functions for particular base types can also be explicitly imported from the LAPACK and GSL.Matrix modules. -} ----------------------------------------------------------------------------- module LinearAlgebra.Algorithms ( --- * Basic Linear Algebra - scale, - add, - multiply, dot, outer, - linearSolve, linearSolveSVD, +-- * Linear Systems + linearSolve, + inv, pinv, + pinvTol, det, -- * Matrix factorizations - svd, lu, eig, eigSH, - qr, chol, --- * Utilities - Normed(..), NormType(..), - det,inv,pinv,full,economy, - pinvTol, --- pinvTolg, +-- ** Singular value decomposition + svd, + full, economy, +-- ** Eigensystems + eig, LinearAlgebra.Algorithms.eigS, LinearAlgebra.Algorithms.eigH, +-- ** Other + LinearAlgebra.Algorithms.qr, chol, +-- * Nullspace nullspacePrec, nullVector, --- * Conversions - toComplex, fromComplex, - comp, real, complex, - conj, ctrans, -- * Misc eps, i, - scaleRecip, - addConstant, - sub, - mul, - divide, + ctrans, + Normed(..), NormType(..), + GenMat(linearSolveSVD,lu,eigSH) ) where import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj) -import Data.Packed.Matrix +import Data.Packed import GSL.Matrix import GSL.Vector import LAPACK import Complex import LinearAlgebra.Linear --- | Base types for which some optimized matrix computations are available -class (Linear Matrix t) => Optimized t where +-- | matrix computations available for both real and complex matrices +class (Linear Matrix t) => GenMat t where svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) linearSolve :: Matrix t -> Matrix t -> Matrix t linearSolveSVD :: Matrix t -> Matrix t -> Matrix t - ctrans :: Matrix t -> Matrix t eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) eigSH :: Matrix t -> (Vector Double, Matrix t) chol :: Matrix t -> Matrix t + -- | conjugate transpose + ctrans :: Matrix t -> Matrix t -instance Optimized Double where +instance GenMat Double where svd = svdR lu = luR linearSolve = linearSolveR linearSolveSVD = linearSolveSVDR Nothing ctrans = trans eig = eigR - eigSH = eigS + eigSH = LAPACK.eigS chol = cholR -instance Optimized (Complex Double) where +instance GenMat (Complex Double) where svd = svdC lu = luC linearSolve = linearSolveC linearSolveSVD = linearSolveSVDC Nothing ctrans = conjTrans eig = eigC - eigSH = eigH + eigSH = LAPACK.eigH chol = error "cholC not yet implemented" -- waiting for GSL-1.10 +-- | eigensystem of a symmetric matrix +eigS :: Matrix Double -> (Vector Double, Matrix Double) +eigS = LAPACK.eigS + +-- | eigensystem of a hermitian matrix +eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +eigH = LAPACK.eigH + +qr :: Matrix Double -> (Matrix Double, Matrix Double) +qr = GSL.Matrix.qr + square m = rows m == cols m -det :: Optimized t => Matrix t -> t +det :: GenMat t => Matrix t -> t det m | square m = s * (product $ toList $ takeDiag $ u) | otherwise = error "det of nonsquare matrix" where (_,u,_,s) = lu m -inv :: Optimized t => Matrix t -> Matrix t +inv :: GenMat t => Matrix t -> Matrix t inv m | square m = m `linearSolve` ident (rows m) | otherwise = error "inv of nonsquare matrix" -pinv :: Optimized t => Matrix t -> Matrix t +pinv :: GenMat t => Matrix t -> Matrix t pinv m = linearSolveSVD m (ident (rows m)) - +full :: Field t + => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) full svd m = (u, d ,v) where (u,s,v) = svd m d = diagRect s r c r = rows m c = cols m +economy :: Field t + => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) economy svd m = (u', subVector 0 d s, v') where (u,s,v) = svd m sl@(g:_) = toList (complex s) @@ -123,39 +131,25 @@ economy svd m = (u', subVector 0 d s, v') where v' = takeColumns d v -{- | Machine precision of a Double. - ->> eps -> 2.22044604925031e-16 - -(The value used by GNU-Octave) - --} +-- | The machine precision of a Double: @eps == 2.22044604925031e-16@ (the value used by GNU-Octave). eps :: Double eps = 2.22044604925031e-16 -{- | The imaginary unit - -@> 'ident' 3 \<\> i -1.i 0. 0. - 0. 1.i 0. - 0. 0. 1.i@ - --} +-- | The imaginary unit: @i == 0.0 :+ 1.0@ i :: Complex Double i = 0:+1 -- | matrix product -mXm :: (Num t, Optimized t) => Matrix t -> Matrix t -> Matrix t +mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t mXm = multiply -- | matrix - vector product -mXv :: (Num t, Optimized t) => Matrix t -> Vector t -> Vector t +mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t mXv m v = flatten $ m `mXm` (asColumn v) -- | vector - matrix product -vXm :: (Num t, Optimized t) => Vector t -> Matrix t -> Vector t +vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t vXm v m = flatten $ (asRow v) `mXm` m @@ -167,15 +161,6 @@ norm2 = toScalarR Norm2 norm1 :: Vector Double -> Double norm1 = toScalarR AbsSum -vectorMax :: Vector Double -> Double -vectorMax = toScalarR Max -vectorMin :: Vector Double -> Double -vectorMin = toScalarR Min -vectorMaxIndex :: Vector Double -> Int -vectorMaxIndex = round . toScalarR MaxIdx -vectorMinIndex :: Vector Double -> Int -vectorMinIndex = round . toScalarR MinIdx - data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int pnormRV PNorm2 = norm2 @@ -199,7 +184,7 @@ pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` const --pnormCM _ _ = error "p norm not yet defined" -- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'. ---pnorm :: (Container t, Optimized a) => Int -> t a -> Double +--pnorm :: (Container t, GenMat a) => Int -> t a -> Double --pnorm = pnormG class Normed t where @@ -222,7 +207,7 @@ instance Normed (Matrix (Complex Double)) where ----------------------------------------------------------------------- -- | The nullspace of a matrix from its SVD decomposition. -nullspacePrec :: Optimized t +nullspacePrec :: GenMat t => Double -- ^ relative tolerance in 'eps' units -> Matrix t -- ^ input matrix -> [Vector t] -- ^ list of unitary vectors spanning the nullspace @@ -235,12 +220,12 @@ nullspacePrec t m = ns where ns = drop rank $ toRows $ ctrans v -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). -nullVector :: Optimized t => Matrix t -> Vector t +nullVector :: GenMat t => Matrix t -> Vector t nullVector = last . nullspacePrec 1 ------------------------------------------------------------------------ -{- | Pseudoinverse of a real matrix with the desired tolerance, expressed as a +{- Pseudoinverse of a real matrix with the desired tolerance, expressed as a multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). @\> let m = 'fromLists' [[1,0, 0] @@ -258,7 +243,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). 0. 0. 1.@ -} -pinvTol :: Double -> Matrix Double -> Matrix Double +--pinvTol :: Double -> Matrix Double -> Matrix Double pinvTol t m = v' `mXm` diag s' `mXm` trans u' where (u,s,v) = svdR m sl@(g:_) = toList s diff --git a/lib/LinearAlgebra/Interface.hs b/lib/LinearAlgebra/Interface.hs index 3392d54..0c65a8b 100644 --- a/lib/LinearAlgebra/Interface.hs +++ b/lib/LinearAlgebra/Interface.hs @@ -16,6 +16,7 @@ Operators for frequent operations. module LinearAlgebra.Interface( (<>),(<.>), + (<\>), (.*),(*/), (<|>),(<->), ) where @@ -23,6 +24,7 @@ module LinearAlgebra.Interface( import LinearAlgebra.Linear import Data.Packed.Vector import Data.Packed.Matrix +import LinearAlgebra.Algorithms class Mul a b c | a b -> c where infixl 7 <> @@ -59,6 +61,11 @@ a .* x = scale a x infixl 7 */ v */ x = scale (recip x) v +-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). +(<\>) :: (GenMat a) => Matrix a -> Vector a -> Vector a +infixl 7 <\> +m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + ------------------------------------------------ class Joinable a b where diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs index 85daa4a..3e6c55d 100644 --- a/lib/LinearAlgebra/Linear.hs +++ b/lib/LinearAlgebra/Linear.hs @@ -9,10 +9,10 @@ Maintainer : Alberto Ruiz (aruiz at um dot es) Stability : provisional Portability : uses ffi +Basic optimized operations on vectors and matrices. -} ----------------------------------------------------------------------------- --- #hide module LinearAlgebra.Linear ( Linear(..), @@ -21,25 +21,22 @@ module LinearAlgebra.Linear ( import Data.Packed.Internal -import Data.Packed.Matrix +import Data.Packed import GSL.Vector import Complex -- | basic optimized operations -class (Field e) => Linear c e where +class (Container c e) => Linear c e where scale :: e -> c e -> c e - scaleRecip :: e -> c e -> c e addConstant :: 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 + -- | element by element division divide :: c e -> c e -> c e - toComplex :: RealFloat e => (c e, c e) -> c (Complex e) - fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) - comp :: RealFloat e => c e -> c (Complex e) - conj :: RealFloat e => c (Complex e) -> c (Complex e) - real :: c Double -> c e - complex :: c e -> c (Complex Double) + -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ + scaleRecip :: e -> c e -> c e instance Linear Vector Double where scale = vectorMapValR Scale @@ -49,12 +46,6 @@ instance Linear Vector Double where sub = vectorZipR Sub mul = vectorZipR Mul divide = vectorZipR Div - toComplex = Data.Packed.Internal.toComplex - fromComplex = Data.Packed.Internal.fromComplex - comp = Data.Packed.Internal.comp - conj = Data.Packed.Internal.conj - real = id - complex = LinearAlgebra.Linear.comp instance Linear Vector (Complex Double) where scale = vectorMapValC Scale @@ -64,12 +55,6 @@ instance Linear Vector (Complex Double) where sub = vectorZipC Sub mul = vectorZipC Mul divide = vectorZipC Div - toComplex = undefined -- can't match - fromComplex = undefined - comp = undefined - conj = undefined - real = LinearAlgebra.Linear.comp - complex = id instance Linear Matrix Double where scale x = liftMatrix (scale x) @@ -79,14 +64,6 @@ instance Linear Matrix Double where sub = liftMatrix2 sub mul = liftMatrix2 mul divide = liftMatrix2 divide - toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex - fromComplex z = (reshape c r, reshape c i) - where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z) - c = cols z - comp = liftMatrix Data.Packed.Internal.comp - conj = liftMatrix Data.Packed.Internal.conj - real = id - complex = LinearAlgebra.Linear.comp instance Linear Matrix (Complex Double) where scale x = liftMatrix (scale x) @@ -96,12 +73,6 @@ instance Linear Matrix (Complex Double) where sub = liftMatrix2 sub mul = liftMatrix2 mul divide = liftMatrix2 divide - toComplex = undefined - fromComplex = undefined - comp = undefined - conj = undefined - real = LinearAlgebra.Linear.comp - complex = id -------------------------------------------------- -- cgit v1.2.3