From a2d99e7d0e83fcedf3a856cdb927309e28a8eddd Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 16 May 2014 12:36:52 +0200 Subject: container and algorithms moved to base --- packages/base/hmatrix-base.cabal | 2 + packages/base/src/Numeric/Container.hs | 242 +++++++ .../base/src/Numeric/LinearAlgebra/Algorithms.hs | 744 ++++++++++++++++++++ packages/hmatrix/hmatrix.cabal | 4 +- packages/hmatrix/src/Numeric/Container.hs | 275 -------- packages/hmatrix/src/Numeric/HMatrix/Data.hs | 1 + .../src/Numeric/LinearAlgebra/Algorithms.hs | 746 --------------------- packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs | 24 +- 8 files changed, 1013 insertions(+), 1025 deletions(-) create mode 100644 packages/base/src/Numeric/Container.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra/Algorithms.hs delete mode 100644 packages/hmatrix/src/Numeric/Container.hs delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 8eeb97e..9a84b7e 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal @@ -37,10 +37,12 @@ library Data.Packed.Development, Numeric.Conversion Numeric.LinearAlgebra.LAPACK + Numeric.LinearAlgebra.Algorithms Data.Packed.Numeric Numeric.Vectorized Numeric.Vector Numeric.Matrix + Numeric.Container other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs new file mode 100644 index 0000000..b7d3b80 --- /dev/null +++ b/packages/base/src/Numeric/Container.hs @@ -0,0 +1,242 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Container +-- Copyright : (c) Alberto Ruiz 2010-14 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. +-- +-- The 'Container' class is used to define optimized generic functions which work +-- on 'Vector' and 'Matrix' with real or complex elements. +-- +-- Some of these functions are also available in the instances of the standard +-- numeric Haskell classes provided by "Numeric.LinearAlgebra". +-- +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + +module Numeric.Container ( + -- * Basic functions + module Data.Packed, + konst, build, + linspace, + diag, ident, + ctrans, + -- * Generic operations + Container(..), + -- * Matrix product + Product(..), udot, dot, (◇), + Mul(..), + Contraction(..), + optimiseMult, + mXm,mXv,vXm,LSDiv(..), + outer, kronecker, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex +) where + +import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) +import Data.Packed.Numeric +import Data.Complex +import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) +import Data.Monoid(Monoid(mconcat)) + +------------------------------------------------------------------ + +{- | Creates a real vector containing a range of values: + +>>> linspace 5 (-3,7::Double) +fromList [-3.0,-0.5,2.0,4.5,7.0]@ + +>>> linspace 5 (8,2+i) :: Vector (Complex Double) +fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] + +Logarithmic spacing can be defined as follows: + +@logspace n (a,b) = 10 ** linspace n (a,b)@ +-} +linspace :: (Container Vector e) => Int -> (e, e) -> Vector e +linspace 0 (a,b) = fromList[(a+b)/2] +linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] + where s = (b-a)/fromIntegral (n-1) + +-------------------------------------------------------- + +class Contraction a b c | a b -> c + where + infixl 7 <.> + {- | Matrix product, matrix vector product, and dot product + +Examples: + +>>> let a = (3><4) [1..] :: Matrix Double +>>> let v = fromList [1,0,2,-1] :: Vector Double +>>> let u = fromList [1,2,3] :: Vector Double + +>>> a +(3><4) + [ 1.0, 2.0, 3.0, 4.0 + , 5.0, 6.0, 7.0, 8.0 + , 9.0, 10.0, 11.0, 12.0 ] + +matrix × matrix: + +>>> disp 2 (a <.> trans a) +3x3 + 30 70 110 + 70 174 278 +110 278 446 + +matrix × vector: + +>>> a <.> v +fromList [3.0,11.0,19.0] + +dot product: + +>>> u <.> fromList[3,2,1::Double] +10 + +For complex vectors the first argument is conjugated: + +>>> fromList [1,i] <.> fromList[2*i+1,3] +1.0 :+ (-1.0) + +>>> fromList [1,i,1-i] <.> complex a +fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] + +-} + (<.>) :: a -> b -> c + + +instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where + u <.> v = conj u `udot` v + +instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where + (<.>) = mXv + +instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where + (<.>) v m = (conj v) `vXm` m + +instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where + (<.>) = mXm + + +-------------------------------------------------------------------------------- + +class Mul a b c | a b -> c where + infixl 7 <> + -- | Matrix-matrix, matrix-vector, and vector-matrix products. + (<>) :: Product t => a t -> b t -> c t + +instance Mul Matrix Matrix Matrix where + (<>) = mXm + +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> asColumn v + +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ asRow v <> m + +-------------------------------------------------------------------------------- + +class LSDiv c where + infixl 7 <\> + -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) + (<\>) :: Field t => Matrix t -> c t -> c t + +instance LSDiv Vector where + m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + +instance LSDiv Matrix where + (<\>) = linearSolveSVD + +-------------------------------------------------------------------------------- + +class Konst e d c | d -> c, c -> d + where + -- | + -- >>> konst 7 3 :: Vector Float + -- fromList [7.0,7.0,7.0] + -- + -- >>> konst i (3::Int,4::Int) + -- (3><4) + -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] + -- + konst :: e -> d -> c e + +instance Container Vector e => Konst e Int Vector + where + konst = konst' + +instance Container Vector e => Konst e (Int,Int) Matrix + where + konst = konst' + +-------------------------------------------------------------------------------- + +class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f + where + -- | + -- >>> build 5 (**2) :: Vector Double + -- fromList [0.0,1.0,4.0,9.0,16.0] + -- + -- Hilbert matrix of order N: + -- + -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double + -- >>> putStr . dispf 2 $ hilb 3 + -- 3x3 + -- 1.00 0.50 0.33 + -- 0.50 0.33 0.25 + -- 0.33 0.25 0.20 + -- + build :: d -> f -> c e + +instance Container Vector e => Build Int (e -> e) Vector e + where + build = build' + +instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e + where + build = build' + +-------------------------------------------------------------------------------- + +{- | alternative operator for '(\<.\>)' + +x25c7, white diamond + +-} +(◇) :: Contraction a b c => a -> b -> c +infixl 7 ◇ +(◇) = (<.>) + +-- | dot product: @cdot u v = 'udot' ('conj' u) v@ +dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t +dot u v = udot (conj u) v + +-------------------------------------------------------------------------------- + +optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t +optimiseMult = mconcat + diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs new file mode 100644 index 0000000..6f40683 --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs @@ -0,0 +1,744 @@ +{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} + +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Algorithms +Copyright : (c) Alberto Ruiz 2006-9 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +High level generic interface to common matrix computations. + +Specific functions for particular base types can also be explicitly +imported from "Numeric.LinearAlgebra.LAPACK". + +-} +----------------------------------------------------------------------------- + +module Numeric.LinearAlgebra.Algorithms ( +-- * Supported types + Field(), +-- * Linear Systems + linearSolve, + luSolve, + cholSolve, + linearSolveLS, + linearSolveSVD, + inv, pinv, pinvTol, + det, invlndet, + rank, rcond, +-- * Matrix factorizations +-- ** Singular value decomposition + svd, + fullSVD, + thinSVD, + compactSVD, + singularValues, + leftSV, rightSV, +-- ** Eigensystems + eig, eigSH, eigSH', + eigenvalues, eigenvaluesSH, eigenvaluesSH', + geigSH', +-- ** QR + qr, rq, qrRaw, qrgr, +-- ** Cholesky + chol, cholSH, mbCholSH, +-- ** Hessenberg + hess, +-- ** Schur + schur, +-- ** LU + lu, luPacked, +-- * Matrix functions + expm, + sqrtm, + matFunc, +-- * Nullspace + nullspacePrec, + nullVector, + nullspaceSVD, + orth, +-- * Norms + Normed(..), NormType(..), + relativeError, +-- * Misc + eps, peps, i, +-- * Util + haussholder, + unpackQR, unpackHess, + ranksv +) where + + +import Data.Packed.Development hiding ((//)) +import Data.Packed +import Numeric.LinearAlgebra.LAPACK as LAPACK +import Data.List(foldl1') +import Data.Array +import Data.Packed.Numeric + + +{- | Generic linear algebra functions for double precision real and complex matrices. + +(Single precision data can be converted using 'single' and 'double'). + +-} +class (Product t, + Convert t, + Container Vector t, + Container Matrix t, + Normed Matrix t, + Normed Vector t, + Floating 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) + sv' :: Matrix t -> Vector Double + luPacked' :: Matrix t -> (Matrix t, [Int]) + luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t + linearSolve' :: Matrix t -> Matrix t -> Matrix t + cholSolve' :: Matrix t -> Matrix t -> Matrix t + linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t + linearSolveLS' :: Matrix t -> Matrix t -> Matrix t + eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) + eigSH'' :: Matrix t -> (Vector Double, Matrix t) + eigOnly :: Matrix t -> Vector (Complex Double) + eigOnlySH :: Matrix t -> Vector Double + cholSH' :: Matrix t -> Matrix t + mbCholSH' :: Matrix t -> Maybe (Matrix t) + qr' :: Matrix t -> (Matrix t, Vector t) + qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t + hess' :: Matrix t -> (Matrix t, Matrix t) + schur' :: Matrix t -> (Matrix t, Matrix t) + + +instance Field Double where + svd' = svdRd + thinSVD' = thinSVDRd + sv' = svR + luPacked' = luR + luSolve' (l_u,perm) = lusR l_u perm + linearSolve' = linearSolveR -- (luSolve . luPacked) ?? + cholSolve' = cholSolveR + linearSolveLS' = linearSolveLSR + linearSolveSVD' = linearSolveSVDR Nothing + eig' = eigR + eigSH'' = eigS + eigOnly = eigOnlyR + eigOnlySH = eigOnlyS + cholSH' = cholS + mbCholSH' = mbCholS + qr' = qrR + qrgr' = qrgrR + hess' = unpackHess hessR + schur' = schurR + +instance Field (Complex Double) where +#ifdef NOZGESDD + svd' = svdC + thinSVD' = thinSVDC +#else + svd' = svdCd + thinSVD' = thinSVDCd +#endif + sv' = svC + luPacked' = luC + luSolve' (l_u,perm) = lusC l_u perm + linearSolve' = linearSolveC + cholSolve' = cholSolveC + linearSolveLS' = linearSolveLSC + linearSolveSVD' = linearSolveSVDC Nothing + eig' = eigC + eigOnly = eigOnlyC + eigSH'' = eigH + eigOnlySH = eigOnlyH + cholSH' = cholH + mbCholSH' = mbCholH + qr' = qrC + qrgr' = qrgrC + hess' = unpackHess hessC + schur' = schurC + +-------------------------------------------------------------- + +square m = rows m == cols m + +vertical m = rows m >= cols m + +exactHermitian m = m `equal` ctrans m + +-------------------------------------------------------------- + +-- | Full singular value decomposition. +svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +svd = {-# SCC "svd" #-} svd' + +-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@. +-- +-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@. +thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +thinSVD = {-# SCC "thinSVD" #-} thinSVD' + +-- | Singular values only. +singularValues :: Field t => Matrix t -> Vector Double +singularValues = {-# SCC "singularValues" #-} sv' + +-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. +-- +-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. +fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) +fullSVD m = (u,d,v) where + (u,s,v) = svd m + d = diagRect 0 s r c + r = rows m + c = cols m + +-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors. +compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) +compactSVD m = (u', subVector 0 d s, v') where + (u,s,v) = thinSVD m + d = rankSVD (1*eps) m s `max` 1 + u' = takeColumns d u + v' = takeColumns d v + + +-- | Singular values and all right singular vectors. +rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) +rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) + | otherwise = let (_,s,v) = svd m in (s,v) + +-- | Singular values and all left singular vectors. +leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) +leftSV m | vertical m = let (u,s,_) = svd m in (u,s) + | otherwise = let (u,s,_) = thinSVD m in (u,s) + + +-------------------------------------------------------------- + +-- | 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' + +-- | 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' + +-- | 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 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 = {-# 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. +linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t +linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD' + + +-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'. +linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t +linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS' + +-------------------------------------------------------------- + +-- | Eigenvalues and eigenvectors of a general square matrix. +-- +-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ +eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) +eig = {-# SCC "eig" #-} eig' + +-- | Eigenvalues of a general square matrix. +eigenvalues :: Field t => Matrix t -> Vector (Complex Double) +eigenvalues = {-# SCC "eigenvalues" #-} eigOnly + +-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. +eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) +eigSH' = {-# SCC "eigSH'" #-} eigSH'' + +-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. +eigenvaluesSH' :: Field t => Matrix t -> Vector Double +eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH + +-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix. +-- +-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ +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" + +-- | Eigenvalues 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" + +-------------------------------------------------------------- + +-- | QR factorization. +-- +-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. +qr :: Field t => Matrix t -> (Matrix t, Matrix t) +qr = {-# SCC "qr" #-} unpackQR . qr' + +qrRaw m = qr' m + +{- | generate a matrix with k orthogonal columns from the output of qrRaw +-} +qrgr n (a,t) + | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)" + | otherwise = qrgr' n (a,t) + +-- | RQ factorization. +-- +-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular. +rq :: Field t => Matrix t -> (Matrix t, Matrix t) +rq m = {-# SCC "rq" #-} (r,q) where + (q',r') = qr $ trans $ rev1 m + r = rev2 (trans r') + q = rev2 (trans q') + rev1 = flipud . fliprl + rev2 = fliprl . flipud + +-- | Hessenberg factorization. +-- +-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary +-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal). +hess :: Field t => Matrix t -> (Matrix t, Matrix t) +hess = hess' + +-- | Schur factorization. +-- +-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary +-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is +-- upper triangular in 2x2 blocks. +-- +-- \"Anything that the Jordan decomposition can do, the Schur decomposition +-- can do better!\" (Van Loan) +schur :: Field t => Matrix t -> (Matrix t, Matrix t) +schur = schur' + + +-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. +mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) +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' + +-- | Cholesky factorization of a positive definite hermitian or symmetric matrix. +-- +-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans 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" + + +-- | Joint computation of inverse and logarithm of determinant of a square matrix. +invlndet :: Field t + => Matrix t + -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) +invlndet m | square m = (im,(ladm,sdm)) + | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" + where + lp@(lup,perm) = luPacked m + s = signlp (rows m) perm + dg = toList $ takeDiag $ lup + ladm = sum $ map (log.abs) dg + sdm = s* product (map signum dg) + im = luSolve lp (ident (rows m)) + + +-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'. +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 + +-- | Explicit LU factorization of a general matrix. +-- +-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, +-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. +lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) +lu = luFact . luPacked + +-- | Inverse of a square matrix. See also 'invlndet'. +inv :: Field t => Matrix t -> Matrix t +inv m | square m = m `linearSolve` ident (rows m) + | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" + + +-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave). +pinv :: Field t => Matrix t -> Matrix t +pinv = pinvTol 1 + +{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value. + +@ +m = (3><3) [ 1, 0, 0 + , 0, 1, 0 + , 0, 0, 1e-10] :: Matrix Double +@ + +>>> pinv m +1. 0. 0. +0. 1. 0. +0. 0. 10000000000. + +>>> pinvTol 1E8 m +1. 0. 0. +0. 1. 0. +0. 0. 1. + +-} + +pinvTol :: Field t => Double -> Matrix t -> Matrix t +pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where + (u,s,v) = thinSVD m + sl@(g:_) = toList s + s' = real . fromList . map rec $ sl + rec x = if x <= g*tol then x else 1/x + tol = (fromIntegral (max r c) * g * t * eps) + r = rows m + c = cols m + d = dim s + u' = takeColumns d u + v' = takeColumns d v + + +-- | Numeric rank of a matrix from the SVD decomposition. +rankSVD :: Element t + => Double -- ^ numeric zero (e.g. 1*'eps') + -> Matrix t -- ^ input matrix m + -> Vector Double -- ^ 'sv' of m + -> Int -- ^ rank of m +rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s) + +-- | Numeric rank of a matrix from its singular values. +ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') + -> Int -- ^ maximum dimension of the matrix + -> [Double] -- ^ singular values + -> Int -- ^ rank of m +ranksv teps maxdim s = k where + g = maximum s + tol = fromIntegral maxdim * g * teps + s' = filter (>tol) s + k = if g > teps then length s' else 0 + +-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave). +eps :: Double +eps = 2.22044604925031e-16 + + +-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 +peps :: RealFloat x => x +peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) + + +-- | The imaginary unit: @i = 0.0 :+ 1.0@ +i :: Complex Double +i = 0:+1 + +----------------------------------------------------------------------- + +-- | The nullspace of a matrix from its precomputed SVD decomposition. +nullspaceSVD :: Field t + => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), + -- or Right \"theoretical\" matrix rank. + -> Matrix t -- ^ input matrix m + -> (Vector Double, Matrix t) -- ^ 'rightSV' of m + -> [Vector t] -- ^ list of unitary vectors spanning the nullspace +nullspaceSVD hint a (s,v) = vs where + tol = case hint of + Left t -> t + _ -> eps + k = case hint of + Right t -> t + _ -> rankSVD tol a s + vs = drop k $ toRows $ ctrans v + + +-- | The nullspace of a matrix. See also 'nullspaceSVD'. +nullspacePrec :: Field t + => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps') + -> Matrix t -- ^ input matrix + -> [Vector t] -- ^ list of unitary vectors spanning the nullspace +nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) + +-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. +nullVector :: Field t => Matrix t -> Vector t +nullVector = last . nullspacePrec 1 + +orth :: Field t => Matrix t -> [Vector t] +-- ^ Return an orthonormal basis of the range space of a matrix +orth m = take r $ toColumns u + where + (u,s,_) = compactSVD m + r = ranksv eps (max (rows m) (cols m)) (toList s) + +------------------------------------------------------------------------ + +-- many thanks, quickcheck! + +haussholder :: (Field a) => a -> Vector a -> Matrix a +haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) + where w = asColumn v + + +zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs) + where xs = toList v + +zt 0 v = v +zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k] + + +unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) +unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r) + where cs = toColumns pq + m = rows pq + n = cols pq + mn = min m n + r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs + vs = zipWith zh [1..mn] cs + hs = zipWith haussholder (toList tau) vs + q = foldl1' mXm hs + +unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) +unpackHess hf m + | rows m == 1 = ((1><1)[1],m) + | otherwise = (uH . hf) m + +uH (pq, tau) = (p,h) + where cs = toColumns pq + m = rows pq + n = cols pq + mn = min m n + h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs + vs = zipWith zh [2..mn] cs + hs = zipWith haussholder (toList tau) vs + p = foldl1' mXm hs + +-------------------------------------------------------------------------- + +-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values. +rcond :: Field t => Matrix t -> Double +rcond m = last s / head s + where s = toList (singularValues m) + +-- | Number of linearly independent rows or columns. +rank :: Field t => Matrix t -> Int +rank m = rankSVD eps m (singularValues m) + +{- +expm' m = case diagonalize (complex m) of + Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v + Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices" + where exp = vectorMapC Exp +-} + +diagonalize m = if rank v == n + then Just (l,v) + else Nothing + where n = rows m + (l,v) = if exactHermitian m + then let (l',v') = eigSH m in (real l', v') + else eig m + +-- | Generic matrix functions for diagonalizable matrices. For instance: +-- +-- @logm = matFunc log@ +-- +matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +matFunc f m = case diagonalize m of + Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v + Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" + +-------------------------------------------------------------- + +golubeps :: Integer -> Integer -> Double +golubeps p q = a * fromIntegral b / fromIntegral c where + a = 2^^(3-p-q) + b = fact p * fact q + c = fact (p+q) * fact (p+q+1) + fact n = product [1..n] + +epslist :: [(Int,Double)] +epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] + +geps delta = head [ k | (k,g) <- epslist, g Matrix t -> Matrix t +expm = expGolub + +expGolub :: Field t => Matrix t -> Matrix t +expGolub m = iterate msq f !! j + where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m + a = m */ fromIntegral ((2::Int)^j) + q = geps eps -- 7 steps + eye = ident (rows m) + work (k,c,x,n,d) = (k',c',x',n',d') + where k' = k+1 + c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k) + x' = a <> x + n' = n |+| (c' .* x') + d' = d |+| (((-1)^k * c') .* x') + (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q + f = linearSolve df nf + msq x = x <> x + + (<>) = multiply + v */ x = scale (recip x) v + (.*) = scale + (|+|) = add + +-------------------------------------------------------------- + +{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. +It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@. + +@m = (2><2) [4,9 + ,0,4] :: Matrix Double@ + +>>> sqrtm m +(2><2) + [ 2.0, 2.25 + , 0.0, 2.0 ] + +-} +sqrtm :: Field t => Matrix t -> Matrix t +sqrtm = sqrtmInv + +sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) + where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a + | otherwise = fixedPoint (b:rest) + fixedPoint _ = error "fixedpoint with impossible inputs" + f (y,z) = (0.5 .* (y |+| inv z), + 0.5 .* (inv y |+| z)) + (.*) = scale + (|+|) = add + (|-|) = sub + +------------------------------------------------------------------ + +signlp r vals = foldl f 1 (zip [0..r-1] vals) + where f s (a,b) | a /= b = -s + | otherwise = s + +swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) + | otherwise = (arr,s) + +fixPerm r vals = (fromColumns $ elems res, sign) + where v = [0..r-1] + s = toColumns (ident r) + (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals) + +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) + where + r = rows l_u + c = cols l_u + tu = triang r c 0 1 + tl = triang r c 0 0 + l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r + u = l_u |*| tu + (p,s) = fixPerm r perm + l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c + u' = takeRows c (l_u |*| tu) + (|+|) = add + (|*|) = mul + +--------------------------------------------------------------------------- + +data NormType = Infinity | PNorm1 | PNorm2 | Frobenius + +class (RealFloat (RealOf t)) => Normed c t where + pnorm :: NormType -> c t -> RealOf t + +instance Normed Vector Double where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = norm2 + +instance Normed Vector (Complex Double) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + +instance Normed Vector Float where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + +instance Normed Vector (Complex Float) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 + + +instance Normed Matrix Double where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Double) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix Float where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Float) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +-- | Approximate number of common digits in the maximum element. +relativeError :: (Normed c t, Container c t) => c t -> c t -> Int +relativeError x y = dig (norm (x `sub` y) / norm x) + where norm = pnorm Infinity + dig r = round $ -logBase 10 (realToFrac r :: Double) + +---------------------------------------------------------------------- + +-- | Generalized symmetric positive definite eigensystem Av = lBv, +-- for A and B symmetric, B positive definite (conditions not checked). +geigSH' :: Field t + => Matrix t -- ^ A + -> Matrix t -- ^ B + -> (Vector Double, Matrix t) +geigSH' a b = (l,v') + where + u = cholSH b + iu = inv u + c = ctrans iu <> a <> iu + (l,v) = eigSH' c + v' = iu <> v + (<>) = mXm + diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal index 369c412..9719c96 100644 --- a/packages/hmatrix/hmatrix.cabal +++ b/packages/hmatrix/hmatrix.cabal @@ -18,7 +18,7 @@ cabal-version: >=1.8 build-type: Simple -extra-source-files: Config.hs THANKS.md INSTALL.md CHANGELOG +extra-source-files: THANKS.md INSTALL.md CHANGELOG extra-source-files: examples/deriv.hs examples/integrate.hs @@ -90,9 +90,7 @@ library Numeric.GSL.Fitting, Numeric.GSL.ODE, Numeric.GSL, - Numeric.Container, Numeric.LinearAlgebra, - Numeric.LinearAlgebra.Algorithms, Numeric.LinearAlgebra.Util, Graphics.Plot, Numeric.HMatrix, diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs deleted file mode 100644 index e7f23d4..0000000 --- a/packages/hmatrix/src/Numeric/Container.hs +++ /dev/null @@ -1,275 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE UndecidableInstances #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Container --- Copyright : (c) Alberto Ruiz 2010-14 --- License : GPL --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. --- --- The 'Container' class is used to define optimized generic functions which work --- on 'Vector' and 'Matrix' with real or complex elements. --- --- Some of these functions are also available in the instances of the standard --- numeric Haskell classes provided by "Numeric.LinearAlgebra". --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Numeric.Container ( - -- * Basic functions - module Data.Packed, - konst, build, - linspace, - diag, ident, - ctrans, - -- * Generic operations - Container(..), - -- * Matrix product - Product(..), udot, dot, (◇), - Mul(..), - Contraction(..), - optimiseMult, - mXm,mXv,vXm,LSDiv(..), - outer, kronecker, - -- * Random numbers - RandDist(..), - randomVector, - gaussianSample, - uniformSample, - meanCov, - -- * Element conversion - Convert(..), - Complexable(), - RealElement(), - - RealOf, ComplexOf, SingleOf, DoubleOf, - - IndexOf, - module Data.Complex, - -- * IO - dispf, disps, dispcf, vecdisp, latexFormat, format, - loadMatrix, saveMatrix, fromFile, fileDimensions, - readMatrix, - fscanfVector, fprintfVector, freadVector, fwriteVector, -) where - -import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) -import Data.Packed.Numeric -import Numeric.IO -import Data.Complex -import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) -import Numeric.Random -import Data.Monoid(Monoid(mconcat)) - ------------------------------------------------------------------- - -{- | Creates a real vector containing a range of values: - ->>> linspace 5 (-3,7::Double) -fromList [-3.0,-0.5,2.0,4.5,7.0]@ - ->>> linspace 5 (8,2+i) :: Vector (Complex Double) -fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] - -Logarithmic spacing can be defined as follows: - -@logspace n (a,b) = 10 ** linspace n (a,b)@ --} -linspace :: (Container Vector e) => Int -> (e, e) -> Vector e -linspace 0 (a,b) = fromList[(a+b)/2] -linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] - where s = (b-a)/fromIntegral (n-1) - --------------------------------------------------------- - -class Contraction a b c | a b -> c - where - infixl 7 <.> - {- | Matrix product, matrix vector product, and dot product - -Examples: - ->>> let a = (3><4) [1..] :: Matrix Double ->>> let v = fromList [1,0,2,-1] :: Vector Double ->>> let u = fromList [1,2,3] :: Vector Double - ->>> a -(3><4) - [ 1.0, 2.0, 3.0, 4.0 - , 5.0, 6.0, 7.0, 8.0 - , 9.0, 10.0, 11.0, 12.0 ] - -matrix × matrix: - ->>> disp 2 (a <.> trans a) -3x3 - 30 70 110 - 70 174 278 -110 278 446 - -matrix × vector: - ->>> a <.> v -fromList [3.0,11.0,19.0] - -dot product: - ->>> u <.> fromList[3,2,1::Double] -10 - -For complex vectors the first argument is conjugated: - ->>> fromList [1,i] <.> fromList[2*i+1,3] -1.0 :+ (-1.0) - ->>> fromList [1,i,1-i] <.> complex a -fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] - --} - (<.>) :: a -> b -> c - - -instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where - u <.> v = conj u `udot` v - -instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where - (<.>) = mXv - -instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where - (<.>) v m = (conj v) `vXm` m - -instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where - (<.>) = mXm - - --------------------------------------------------------------------------------- - -class Mul a b c | a b -> c where - infixl 7 <> - -- | Matrix-matrix, matrix-vector, and vector-matrix products. - (<>) :: Product t => a t -> b t -> c t - -instance Mul Matrix Matrix Matrix where - (<>) = mXm - -instance Mul Matrix Vector Vector where - (<>) m v = flatten $ m <> asColumn v - -instance Mul Vector Matrix Vector where - (<>) v m = flatten $ asRow v <> m - --------------------------------------------------------------------------------- - -class LSDiv c where - infixl 7 <\> - -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) - (<\>) :: Field t => Matrix t -> c t -> c t - -instance LSDiv Vector where - m <\> v = flatten (linearSolveSVD m (reshape 1 v)) - -instance LSDiv Matrix where - (<\>) = linearSolveSVD - --------------------------------------------------------------------------------- - -class Konst e d c | d -> c, c -> d - where - -- | - -- >>> konst 7 3 :: Vector Float - -- fromList [7.0,7.0,7.0] - -- - -- >>> konst i (3::Int,4::Int) - -- (3><4) - -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 - -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 - -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] - -- - konst :: e -> d -> c e - -instance Container Vector e => Konst e Int Vector - where - konst = konst' - -instance Container Vector e => Konst e (Int,Int) Matrix - where - konst = konst' - --------------------------------------------------------------------------------- - -class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f - where - -- | - -- >>> build 5 (**2) :: Vector Double - -- fromList [0.0,1.0,4.0,9.0,16.0] - -- - -- Hilbert matrix of order N: - -- - -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double - -- >>> putStr . dispf 2 $ hilb 3 - -- 3x3 - -- 1.00 0.50 0.33 - -- 0.50 0.33 0.25 - -- 0.33 0.25 0.20 - -- - build :: d -> f -> c e - -instance Container Vector e => Build Int (e -> e) Vector e - where - build = build' - -instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e - where - build = build' - --------------------------------------------------------------------------------- - -{- | Compute mean vector and covariance matrix of the rows of a matrix. - ->>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) -(fromList [4.010341078059521,5.0197204699640405], -(2><2) - [ 1.9862461923890056, -1.0127225830525157e-2 - , -1.0127225830525157e-2, 3.0373954915729318 ]) - --} -meanCov :: Matrix Double -> (Vector Double, Matrix Double) -meanCov x = (med,cov) where - r = rows x - k = 1 / fromIntegral r - med = konst k r `vXm` x - meds = konst 1 r `outer` med - xc = x `sub` meds - cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) - --------------------------------------------------------------------------------- - -{- | alternative operator for '(\<.\>)' - -x25c7, white diamond - --} -(◇) :: Contraction a b c => a -> b -> c -infixl 7 ◇ -(◇) = (<.>) - --- | dot product: @cdot u v = 'udot' ('conj' u) v@ -dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t -dot u v = udot (conj u) v - --------------------------------------------------------------------------------- - -optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t -optimiseMult = mconcat - diff --git a/packages/hmatrix/src/Numeric/HMatrix/Data.hs b/packages/hmatrix/src/Numeric/HMatrix/Data.hs index 568dc05..5d7ce4f 100644 --- a/packages/hmatrix/src/Numeric/HMatrix/Data.hs +++ b/packages/hmatrix/src/Numeric/HMatrix/Data.hs @@ -64,6 +64,7 @@ module Numeric.HMatrix.Data( import Data.Packed.Vector import Data.Packed.Matrix import Numeric.Container +import Numeric.IO import Numeric.LinearAlgebra.Util import Data.Complex diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs deleted file mode 100644 index b4d0c5b..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs +++ /dev/null @@ -1,746 +0,0 @@ -{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE TypeFamilies #-} - ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Algorithms -Copyright : (c) Alberto Ruiz 2006-9 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : uses ffi - -High level generic interface to common matrix computations. - -Specific functions for particular base types can also be explicitly -imported from "Numeric.LinearAlgebra.LAPACK". - --} ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - - -module Numeric.LinearAlgebra.Algorithms ( --- * Supported types - Field(), --- * Linear Systems - linearSolve, - luSolve, - cholSolve, - linearSolveLS, - linearSolveSVD, - inv, pinv, pinvTol, - det, invlndet, - rank, rcond, --- * Matrix factorizations --- ** Singular value decomposition - svd, - fullSVD, - thinSVD, - compactSVD, - singularValues, - leftSV, rightSV, --- ** Eigensystems - eig, eigSH, eigSH', - eigenvalues, eigenvaluesSH, eigenvaluesSH', - geigSH', --- ** QR - qr, rq, qrRaw, qrgr, --- ** Cholesky - chol, cholSH, mbCholSH, --- ** Hessenberg - hess, --- ** Schur - schur, --- ** LU - lu, luPacked, --- * Matrix functions - expm, - sqrtm, - matFunc, --- * Nullspace - nullspacePrec, - nullVector, - nullspaceSVD, - orth, --- * Norms - Normed(..), NormType(..), - relativeError, --- * Misc - eps, peps, i, --- * Util - haussholder, - unpackQR, unpackHess, - ranksv -) where - - -import Data.Packed.Development hiding ((//)) -import Data.Packed -import Numeric.LinearAlgebra.LAPACK as LAPACK -import Data.List(foldl1') -import Data.Array -import Data.Packed.Numeric - - -{- | Generic linear algebra functions for double precision real and complex matrices. - -(Single precision data can be converted using 'single' and 'double'). - --} -class (Product t, - Convert t, - Container Vector t, - Container Matrix t, - Normed Matrix t, - Normed Vector t, - Floating 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) - sv' :: Matrix t -> Vector Double - luPacked' :: Matrix t -> (Matrix t, [Int]) - luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t - linearSolve' :: Matrix t -> Matrix t -> Matrix t - cholSolve' :: Matrix t -> Matrix t -> Matrix t - linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t - linearSolveLS' :: Matrix t -> Matrix t -> Matrix t - eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) - eigSH'' :: Matrix t -> (Vector Double, Matrix t) - eigOnly :: Matrix t -> Vector (Complex Double) - eigOnlySH :: Matrix t -> Vector Double - cholSH' :: Matrix t -> Matrix t - mbCholSH' :: Matrix t -> Maybe (Matrix t) - qr' :: Matrix t -> (Matrix t, Vector t) - qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t - hess' :: Matrix t -> (Matrix t, Matrix t) - schur' :: Matrix t -> (Matrix t, Matrix t) - - -instance Field Double where - svd' = svdRd - thinSVD' = thinSVDRd - sv' = svR - luPacked' = luR - luSolve' (l_u,perm) = lusR l_u perm - linearSolve' = linearSolveR -- (luSolve . luPacked) ?? - cholSolve' = cholSolveR - linearSolveLS' = linearSolveLSR - linearSolveSVD' = linearSolveSVDR Nothing - eig' = eigR - eigSH'' = eigS - eigOnly = eigOnlyR - eigOnlySH = eigOnlyS - cholSH' = cholS - mbCholSH' = mbCholS - qr' = qrR - qrgr' = qrgrR - hess' = unpackHess hessR - schur' = schurR - -instance Field (Complex Double) where -#ifdef NOZGESDD - svd' = svdC - thinSVD' = thinSVDC -#else - svd' = svdCd - thinSVD' = thinSVDCd -#endif - sv' = svC - luPacked' = luC - luSolve' (l_u,perm) = lusC l_u perm - linearSolve' = linearSolveC - cholSolve' = cholSolveC - linearSolveLS' = linearSolveLSC - linearSolveSVD' = linearSolveSVDC Nothing - eig' = eigC - eigOnly = eigOnlyC - eigSH'' = eigH - eigOnlySH = eigOnlyH - cholSH' = cholH - mbCholSH' = mbCholH - qr' = qrC - qrgr' = qrgrC - hess' = unpackHess hessC - schur' = schurC - --------------------------------------------------------------- - -square m = rows m == cols m - -vertical m = rows m >= cols m - -exactHermitian m = m `equal` ctrans m - --------------------------------------------------------------- - --- | Full singular value decomposition. -svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) -svd = {-# SCC "svd" #-} svd' - --- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@. --- --- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@. -thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) -thinSVD = {-# SCC "thinSVD" #-} thinSVD' - --- | Singular values only. -singularValues :: Field t => Matrix t -> Vector Double -singularValues = {-# SCC "singularValues" #-} sv' - --- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. --- --- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. -fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) -fullSVD m = (u,d,v) where - (u,s,v) = svd m - d = diagRect 0 s r c - r = rows m - c = cols m - --- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors. -compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) -compactSVD m = (u', subVector 0 d s, v') where - (u,s,v) = thinSVD m - d = rankSVD (1*eps) m s `max` 1 - u' = takeColumns d u - v' = takeColumns d v - - --- | Singular values and all right singular vectors. -rightSV :: Field t => Matrix t -> (Vector Double, Matrix t) -rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v) - | otherwise = let (_,s,v) = svd m in (s,v) - --- | Singular values and all left singular vectors. -leftSV :: Field t => Matrix t -> (Matrix t, Vector Double) -leftSV m | vertical m = let (u,s,_) = svd m in (u,s) - | otherwise = let (u,s,_) = thinSVD m in (u,s) - - --------------------------------------------------------------- - --- | 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' - --- | 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' - --- | 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 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 = {-# 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. -linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t -linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD' - - --- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'. -linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t -linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS' - --------------------------------------------------------------- - --- | Eigenvalues and eigenvectors of a general square matrix. --- --- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ -eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) -eig = {-# SCC "eig" #-} eig' - --- | Eigenvalues of a general square matrix. -eigenvalues :: Field t => Matrix t -> Vector (Complex Double) -eigenvalues = {-# SCC "eigenvalues" #-} eigOnly - --- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. -eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) -eigSH' = {-# SCC "eigSH'" #-} eigSH'' - --- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. -eigenvaluesSH' :: Field t => Matrix t -> Vector Double -eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH - --- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix. --- --- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ -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" - --- | Eigenvalues 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" - --------------------------------------------------------------- - --- | QR factorization. --- --- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. -qr :: Field t => Matrix t -> (Matrix t, Matrix t) -qr = {-# SCC "qr" #-} unpackQR . qr' - -qrRaw m = qr' m - -{- | generate a matrix with k orthogonal columns from the output of qrRaw --} -qrgr n (a,t) - | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)" - | otherwise = qrgr' n (a,t) - --- | RQ factorization. --- --- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular. -rq :: Field t => Matrix t -> (Matrix t, Matrix t) -rq m = {-# SCC "rq" #-} (r,q) where - (q',r') = qr $ trans $ rev1 m - r = rev2 (trans r') - q = rev2 (trans q') - rev1 = flipud . fliprl - rev2 = fliprl . flipud - --- | Hessenberg factorization. --- --- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary --- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal). -hess :: Field t => Matrix t -> (Matrix t, Matrix t) -hess = hess' - --- | Schur factorization. --- --- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary --- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is --- upper triangular in 2x2 blocks. --- --- \"Anything that the Jordan decomposition can do, the Schur decomposition --- can do better!\" (Van Loan) -schur :: Field t => Matrix t -> (Matrix t, Matrix t) -schur = schur' - - --- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. -mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) -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' - --- | Cholesky factorization of a positive definite hermitian or symmetric matrix. --- --- If @c = chol m@ then @c@ is upper triangular and @m == ctrans 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" - - --- | Joint computation of inverse and logarithm of determinant of a square matrix. -invlndet :: Field t - => Matrix t - -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) -invlndet m | square m = (im,(ladm,sdm)) - | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" - where - lp@(lup,perm) = luPacked m - s = signlp (rows m) perm - dg = toList $ takeDiag $ lup - ladm = sum $ map (log.abs) dg - sdm = s* product (map signum dg) - im = luSolve lp (ident (rows m)) - - --- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'. -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 - --- | Explicit LU factorization of a general matrix. --- --- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, --- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. -lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) -lu = luFact . luPacked - --- | Inverse of a square matrix. See also 'invlndet'. -inv :: Field t => Matrix t -> Matrix t -inv m | square m = m `linearSolve` ident (rows m) - | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" - - --- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave). -pinv :: Field t => Matrix t -> Matrix t -pinv = pinvTol 1 - -{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value. - -@ -m = (3><3) [ 1, 0, 0 - , 0, 1, 0 - , 0, 0, 1e-10] :: Matrix Double -@ - ->>> pinv m -1. 0. 0. -0. 1. 0. -0. 0. 10000000000. - ->>> pinvTol 1E8 m -1. 0. 0. -0. 1. 0. -0. 0. 1. - --} - -pinvTol :: Field t => Double -> Matrix t -> Matrix t -pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where - (u,s,v) = thinSVD m - sl@(g:_) = toList s - s' = real . fromList . map rec $ sl - rec x = if x <= g*tol then x else 1/x - tol = (fromIntegral (max r c) * g * t * eps) - r = rows m - c = cols m - d = dim s - u' = takeColumns d u - v' = takeColumns d v - - --- | Numeric rank of a matrix from the SVD decomposition. -rankSVD :: Element t - => Double -- ^ numeric zero (e.g. 1*'eps') - -> Matrix t -- ^ input matrix m - -> Vector Double -- ^ 'sv' of m - -> Int -- ^ rank of m -rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s) - --- | Numeric rank of a matrix from its singular values. -ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') - -> Int -- ^ maximum dimension of the matrix - -> [Double] -- ^ singular values - -> Int -- ^ rank of m -ranksv teps maxdim s = k where - g = maximum s - tol = fromIntegral maxdim * g * teps - s' = filter (>tol) s - k = if g > teps then length s' else 0 - --- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave). -eps :: Double -eps = 2.22044604925031e-16 - - --- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 -peps :: RealFloat x => x -peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x) - - --- | The imaginary unit: @i = 0.0 :+ 1.0@ -i :: Complex Double -i = 0:+1 - ------------------------------------------------------------------------ - --- | The nullspace of a matrix from its precomputed SVD decomposition. -nullspaceSVD :: Field t - => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), - -- or Right \"theoretical\" matrix rank. - -> Matrix t -- ^ input matrix m - -> (Vector Double, Matrix t) -- ^ 'rightSV' of m - -> [Vector t] -- ^ list of unitary vectors spanning the nullspace -nullspaceSVD hint a (s,v) = vs where - tol = case hint of - Left t -> t - _ -> eps - k = case hint of - Right t -> t - _ -> rankSVD tol a s - vs = drop k $ toRows $ ctrans v - - --- | The nullspace of a matrix. See also 'nullspaceSVD'. -nullspacePrec :: Field t - => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps') - -> Matrix t -- ^ input matrix - -> [Vector t] -- ^ list of unitary vectors spanning the nullspace -nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) - --- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision. -nullVector :: Field t => Matrix t -> Vector t -nullVector = last . nullspacePrec 1 - -orth :: Field t => Matrix t -> [Vector t] --- ^ Return an orthonormal basis of the range space of a matrix -orth m = take r $ toColumns u - where - (u,s,_) = compactSVD m - r = ranksv eps (max (rows m) (cols m)) (toList s) - ------------------------------------------------------------------------- - --- many thanks, quickcheck! - -haussholder :: (Field a) => a -> Vector a -> Matrix a -haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w)) - where w = asColumn v - - -zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs) - where xs = toList v - -zt 0 v = v -zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k] - - -unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) -unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r) - where cs = toColumns pq - m = rows pq - n = cols pq - mn = min m n - r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs - vs = zipWith zh [1..mn] cs - hs = zipWith haussholder (toList tau) vs - q = foldl1' mXm hs - -unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) -unpackHess hf m - | rows m == 1 = ((1><1)[1],m) - | otherwise = (uH . hf) m - -uH (pq, tau) = (p,h) - where cs = toColumns pq - m = rows pq - n = cols pq - mn = min m n - h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs - vs = zipWith zh [2..mn] cs - hs = zipWith haussholder (toList tau) vs - p = foldl1' mXm hs - --------------------------------------------------------------------------- - --- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values. -rcond :: Field t => Matrix t -> Double -rcond m = last s / head s - where s = toList (singularValues m) - --- | Number of linearly independent rows or columns. -rank :: Field t => Matrix t -> Int -rank m = rankSVD eps m (singularValues m) - -{- -expm' m = case diagonalize (complex m) of - Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v - Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices" - where exp = vectorMapC Exp --} - -diagonalize m = if rank v == n - then Just (l,v) - else Nothing - where n = rows m - (l,v) = if exactHermitian m - then let (l',v') = eigSH m in (real l', v') - else eig m - --- | Generic matrix functions for diagonalizable matrices. For instance: --- --- @logm = matFunc log@ --- -matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) -matFunc f m = case diagonalize m of - Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v - Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" - --------------------------------------------------------------- - -golubeps :: Integer -> Integer -> Double -golubeps p q = a * fromIntegral b / fromIntegral c where - a = 2^^(3-p-q) - b = fact p * fact q - c = fact (p+q) * fact (p+q+1) - fact n = product [1..n] - -epslist :: [(Int,Double)] -epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] - -geps delta = head [ k | (k,g) <- epslist, g Matrix t -> Matrix t -expm = expGolub - -expGolub :: Field t => Matrix t -> Matrix t -expGolub m = iterate msq f !! j - where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m - a = m */ fromIntegral ((2::Int)^j) - q = geps eps -- 7 steps - eye = ident (rows m) - work (k,c,x,n,d) = (k',c',x',n',d') - where k' = k+1 - c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k) - x' = a <> x - n' = n |+| (c' .* x') - d' = d |+| (((-1)^k * c') .* x') - (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q - f = linearSolve df nf - msq x = x <> x - - (<>) = multiply - v */ x = scale (recip x) v - (.*) = scale - (|+|) = add - --------------------------------------------------------------- - -{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. -It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@. - -@m = (2><2) [4,9 - ,0,4] :: Matrix Double@ - ->>> sqrtm m -(2><2) - [ 2.0, 2.25 - , 0.0, 2.0 ] - --} -sqrtm :: Field t => Matrix t -> Matrix t -sqrtm = sqrtmInv - -sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) - where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a - | otherwise = fixedPoint (b:rest) - fixedPoint _ = error "fixedpoint with impossible inputs" - f (y,z) = (0.5 .* (y |+| inv z), - 0.5 .* (inv y |+| z)) - (.*) = scale - (|+|) = add - (|-|) = sub - ------------------------------------------------------------------- - -signlp r vals = foldl f 1 (zip [0..r-1] vals) - where f s (a,b) | a /= b = -s - | otherwise = s - -swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s) - | otherwise = (arr,s) - -fixPerm r vals = (fromColumns $ elems res, sign) - where v = [0..r-1] - s = toColumns (ident r) - (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals) - -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) - where - r = rows l_u - c = cols l_u - tu = triang r c 0 1 - tl = triang r c 0 0 - l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r - u = l_u |*| tu - (p,s) = fixPerm r perm - l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c - u' = takeRows c (l_u |*| tu) - (|+|) = add - (|*|) = mul - ---------------------------------------------------------------------------- - -data NormType = Infinity | PNorm1 | PNorm2 | Frobenius - -class (RealFloat (RealOf t)) => Normed c t where - pnorm :: NormType -> c t -> RealOf t - -instance Normed Vector Double where - pnorm PNorm1 = norm1 - pnorm PNorm2 = norm2 - pnorm Infinity = normInf - pnorm Frobenius = norm2 - -instance Normed Vector (Complex Double) where - pnorm PNorm1 = norm1 - pnorm PNorm2 = norm2 - pnorm Infinity = normInf - pnorm Frobenius = pnorm PNorm2 - -instance Normed Vector Float where - pnorm PNorm1 = norm1 - pnorm PNorm2 = norm2 - pnorm Infinity = normInf - pnorm Frobenius = pnorm PNorm2 - -instance Normed Vector (Complex Float) where - pnorm PNorm1 = norm1 - pnorm PNorm2 = norm2 - pnorm Infinity = normInf - pnorm Frobenius = pnorm PNorm2 - - -instance Normed Matrix Double where - pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns - pnorm PNorm2 = (@>0) . singularValues - pnorm Infinity = pnorm PNorm1 . trans - pnorm Frobenius = pnorm PNorm2 . flatten - -instance Normed Matrix (Complex Double) where - pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns - pnorm PNorm2 = (@>0) . singularValues - pnorm Infinity = pnorm PNorm1 . trans - pnorm Frobenius = pnorm PNorm2 . flatten - -instance Normed Matrix Float where - pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns - pnorm PNorm2 = realToFrac . (@>0) . singularValues . double - pnorm Infinity = pnorm PNorm1 . trans - pnorm Frobenius = pnorm PNorm2 . flatten - -instance Normed Matrix (Complex Float) where - pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns - pnorm PNorm2 = realToFrac . (@>0) . singularValues . double - pnorm Infinity = pnorm PNorm1 . trans - pnorm Frobenius = pnorm PNorm2 . flatten - --- | Approximate number of common digits in the maximum element. -relativeError :: (Normed c t, Container c t) => c t -> c t -> Int -relativeError x y = dig (norm (x `sub` y) / norm x) - where norm = pnorm Infinity - dig r = round $ -logBase 10 (realToFrac r :: Double) - ----------------------------------------------------------------------- - --- | Generalized symmetric positive definite eigensystem Av = lBv, --- for A and B symmetric, B positive definite (conditions not checked). -geigSH' :: Field t - => Matrix t -- ^ A - -> Matrix t -- ^ B - -> (Vector Double, Matrix t) -geigSH' a b = (l,v') - where - u = cholSH b - iu = inv u - c = ctrans iu <> a <> iu - (l,v) = eigSH' c - v' = iu <> v - (<>) = mXm - diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs index fb2d54c..e4f21b0 100644 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs @@ -27,6 +27,7 @@ module Numeric.LinearAlgebra.Util( unitary, mt, pairwiseD2, + meanCov, rowOuters, null1, null1sym, @@ -55,6 +56,7 @@ module Numeric.LinearAlgebra.Util( ) where import Numeric.Container +import Numeric.IO import Numeric.LinearAlgebra.Algorithms hiding (i) import Numeric.Matrix() import Numeric.Vector() @@ -196,7 +198,27 @@ size m = (rows m, cols m) mt :: Matrix Double -> Matrix Double mt = trans . inv ----------------------------------------------------------------------- +-------------------------------------------------------------------------------- + +{- | Compute mean vector and covariance matrix of the rows of a matrix. + +>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) +(fromList [4.010341078059521,5.0197204699640405], +(2><2) + [ 1.9862461923890056, -1.0127225830525157e-2 + , -1.0127225830525157e-2, 3.0373954915729318 ]) + +-} +meanCov :: Matrix Double -> (Vector Double, Matrix Double) +meanCov x = (med,cov) where + r = rows x + k = 1 / fromIntegral r + med = konst k r `vXm` x + meds = konst 1 r `outer` med + xc = x `sub` meds + cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) + +-------------------------------------------------------------------------------- -- | Matrix of pairwise squared distances of row vectors -- (using the matrix product trick in blog.smola.org) -- cgit v1.2.3