From 11d7c37dc8b314338bc6382d80e74aaec2bb5620 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 5 Jun 2015 16:43:45 +0200 Subject: move algorithms --- .../base/src/Numeric/LinearAlgebra/Algorithms.hs | 972 --------------------- 1 file changed, 972 deletions(-) delete mode 100644 packages/base/src/Numeric/LinearAlgebra/Algorithms.hs (limited to 'packages/base/src/Numeric') diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs deleted file mode 100644 index 0e085c1..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs +++ /dev/null @@ -1,972 +0,0 @@ -{-# 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" #-} g . svd' - where - g (u,s,v) = (u,s,tr v) - -{- | 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" #-} g . thinSVD' - where - g (u,s,v) = (u,s,tr v) - - --- | 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 = 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 = 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 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. - -@m = (2><2) [4,9 - ,0,4] :: Matrix Double@ - ->>> sqrtm m -(2><2) - [ 2.0, 2.25 - , 0.0, 2.0 ] - -For diagonalizable matrices you can try 'matFunc' @sqrt@: - ->>> matFunc sqrt ((2><2) [1,0,0,-1]) -(2><2) - [ 1.0 :+ 0.0, 0.0 :+ 0.0 - , 0.0 :+ 0.0, 0.0 :+ 1.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) - - -relativeError :: Num a => (a -> Double) -> a -> a -> Double -relativeError norm a b = r - where - 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 - -- cgit v1.2.3