hmatrix-0.16.0.2: src/Numeric/LinearAlgebra/Algorithms.hs
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}
-----------------------------------------------------------------------------
{- |
Module : Numeric.LinearAlgebra.Algorithms
Copyright : (c) Alberto Ruiz 2006-14
License : BSD3
Maintainer : Alberto Ruiz
Stability : provisional
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,
mbLinearSolve,
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,
orthSVD,
orth,
-- * Norms
Normed(..), NormType(..),
relativeError', relativeError,
-- * Misc
eps, peps, i,
-- * Util
haussholder,
unpackQR, unpackHess,
ranksv
) where
import Data.Packed
import Numeric.LinearAlgebra.LAPACK as LAPACK
import Data.List(foldl1')
import Data.Array
import Data.Packed.Internal.Numeric
import Data.Packed.Internal(shSize)
{- | 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
mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (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) ??
mbLinearSolve' = mbLinearSolveR
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
mbLinearSolve' = mbLinearSolveC
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.
@
a = (5><3)
[ 1.0, 2.0, 3.0
, 4.0, 5.0, 6.0
, 7.0, 8.0, 9.0
, 10.0, 11.0, 12.0
, 13.0, 14.0, 15.0 ] :: Matrix Double
@
>>> let (u,s,v) = svd a
>>> disp 3 u
5x5
-0.101 0.768 0.614 0.028 -0.149
-0.249 0.488 -0.503 0.172 0.646
-0.396 0.208 -0.405 -0.660 -0.449
-0.543 -0.072 -0.140 0.693 -0.447
-0.690 -0.352 0.433 -0.233 0.398
>>> s
fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
>>> disp 3 v
3x3
-0.519 -0.751 0.408
-0.576 -0.046 -0.816
-0.632 0.659 0.408
>>> let d = diagRect 0 s 5 3
>>> disp 3 d
5x3
35.183 0.000 0.000
0.000 1.477 0.000
0.000 0.000 0.000
0.000 0.000 0.000
>>> disp 3 $ u <> d <> tr v
5x3
1.000 2.000 3.000
4.000 5.000 6.000
7.000 8.000 9.000
10.000 11.000 12.000
13.000 14.000 15.000
-}
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 \<> tr v@.
@
a = (5><3)
[ 1.0, 2.0, 3.0
, 4.0, 5.0, 6.0
, 7.0, 8.0, 9.0
, 10.0, 11.0, 12.0
, 13.0, 14.0, 15.0 ] :: Matrix Double
@
>>> let (u,s,v) = thinSVD a
>>> disp 3 u
5x3
-0.101 0.768 0.614
-0.249 0.488 -0.503
-0.396 0.208 -0.405
-0.543 -0.072 -0.140
-0.690 -0.352 0.433
>>> s
fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
>>> disp 3 v
3x3
-0.519 -0.751 0.408
-0.576 -0.046 -0.816
-0.632 0.659 0.408
>>> disp 3 $ u <> diag s <> tr v
5x3
1.000 2.000 3.000
4.000 5.000 6.000
7.000 8.000 9.000
10.000 11.000 12.000
13.000 14.000 15.000
-}
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 \<> tr 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.
@
a = (5><3)
[ 1.0, 2.0, 3.0
, 4.0, 5.0, 6.0
, 7.0, 8.0, 9.0
, 10.0, 11.0, 12.0
, 13.0, 14.0, 15.0 ] :: Matrix Double
@
>>> let (u,s,v) = compactSVD a
>>> disp 3 u
5x2
-0.101 0.768
-0.249 0.488
-0.396 0.208
-0.543 -0.072
-0.690 -0.352
>>> s
fromList [35.18264833189422,1.4769076999800903]
>>> disp 3 u
5x2
-0.101 0.768
-0.249 0.488
-0.396 0.208
-0.543 -0.072
-0.690 -0.352
>>> disp 3 $ u <> diag s <> tr v
5x3
1.000 2.000 3.000
4.000 5.000 6.000
7.000 8.000 9.000
10.000 11.000 12.000
13.000 14.000 15.000
-}
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 (as columns).
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 (as columns).
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 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 = {-# 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 (not ordered) and eigenvectors (as columns) of a general square matrix.
If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
@
a = (3><3)
[ 3, 0, -2
, 4, 5, -1
, 3, 1, 0 ] :: Matrix Double
@
>>> let (l, v) = eig a
>>> putStr . dispcf 3 . asRow $ l
1x3
1.925+1.523i 1.925-1.523i 4.151
>>> putStr . dispcf 3 $ v
3x3
-0.455+0.365i -0.455-0.365i 0.181
0.603 0.603 -0.978
0.033+0.543i 0.033-0.543i -0.104
>>> putStr . dispcf 3 $ complex a <> v
3x3
-1.432+0.010i -1.432-0.010i 0.753
1.160+0.918i 1.160-0.918i -4.059
-0.763+1.096i -0.763-1.096i -0.433
>>> putStr . dispcf 3 $ v <> diag l
3x3
-1.432+0.010i -1.432-0.010i 0.753
1.160+0.918i 1.160-0.918i -4.059
-0.763+1.096i -0.763-1.096i -0.433
-}
eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
eig = {-# SCC "eig" #-} eig'
-- | Eigenvalues (not ordered) 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 (as columns) of a complex hermitian or real symmetric matrix, in descending order.
If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@
@
a = (3><3)
[ 1.0, 2.0, 3.0
, 2.0, 4.0, 5.0
, 3.0, 5.0, 6.0 ]
@
>>> let (l, v) = eigSH a
>>> l
fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
>>> disp 3 $ v <> diag l <> tr v
3x3
1.000 2.000 3.000
2.000 4.000 5.000
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"
-- | 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"
--------------------------------------------------------------
-- | 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
-> Matrix t -- ^ 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 = conj (dropColumns k 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 = toColumns $ 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
-- | The range space a matrix from its precomputed SVD decomposition.
orthSVD :: Field t
=> Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
-- or Right \"theoretical\" matrix rank.
-> Matrix t -- ^ input matrix m
-> (Matrix t, Vector Double) -- ^ 'leftSV' of m
-> Matrix t -- ^ orth
orthSVD hint a (v,s) = vs where
tol = case hint of
Left t -> t
_ -> eps
k = case hint of
Right t -> t
_ -> rankSVD tol a s
vs = takeColumns k v
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. See also 'ranksv'
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<delta]
{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
based on a scaled Pade approximation.
-}
expm :: Field t => 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><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
where el p q = if q-p>=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)
relativeError :: (Normed c t, Num (c t)) => NormType -> c t -> c t -> Double
relativeError t a b = realToFrac r
where
norm = pnorm t
na = norm a
nb = norm b
nab = norm (a-b)
mx = max na nb
mn = min na nb
r = if mn < peps
then mx
else nab/mx
----------------------------------------------------------------------
-- | 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