summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CHANGES8
-rw-r--r--hmatrix.cabal4
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs284
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs7
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs22
5 files changed, 191 insertions, 134 deletions
diff --git a/CHANGES b/CHANGES
index deca1d5..f786fdf 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,11 +1,15 @@
10.7.3.0 10.8.0.0
2======= 2=======
3 3
4- added sv, fullSVD, thinSVD, compactSVD, leftSV, rightSV 4- singularValues, fullSVD, thinSVD, compactSVD, leftSV, rightSV
5 and complete interface to [d|z]gesdd. 5 and complete interface to [d|z]gesdd.
6 Algorithms based on the SVD of large matrices can now be 6 Algorithms based on the SVD of large matrices can now be
7 significantly faster. 7 significantly faster.
8 8
9- eigenvalues, eigenvaluesSH
10
11- linearSolveLS, rq
12
90.7.2.0 130.7.2.0
10======= 14=======
11 15
diff --git a/hmatrix.cabal b/hmatrix.cabal
index a7dca00..1b077e3 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
1Name: hmatrix 1Name: hmatrix
2Version: 0.7.3.0 2Version: 0.8.0.0
3License: GPL 3License: GPL
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
@@ -11,7 +11,7 @@ Description: Purely functional interface to basic linear algebra
11 and other numerical computations, internally implemented using 11 and other numerical computations, internally implemented using
12 GSL, BLAS and LAPACK. 12 GSL, BLAS and LAPACK.
13Category: Math 13Category: Math
14tested-with: GHC ==6.12.1 14tested-with: GHC ==6.10.4, GHC ==6.12.1
15 15
16cabal-version: >=1.2 16cabal-version: >=1.2
17build-type: Custom 17build-type: Custom
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 3d4a209..de88dd9 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -27,6 +27,7 @@ module Numeric.LinearAlgebra.Algorithms (
27-- * Linear Systems 27-- * Linear Systems
28 linearSolve, 28 linearSolve,
29 luSolve, 29 luSolve,
30 linearSolveLS,
30 linearSolveSVD, 31 linearSolveSVD,
31 inv, pinv, 32 inv, pinv,
32 det, rank, rcond, 33 det, rank, rcond,
@@ -36,13 +37,13 @@ module Numeric.LinearAlgebra.Algorithms (
36 fullSVD, 37 fullSVD,
37 thinSVD, 38 thinSVD,
38 compactSVD, 39 compactSVD,
39 sv, 40 singularValues,
40 leftSV, rightSV, 41 leftSV, rightSV,
41-- ** Eigensystems 42-- ** Eigensystems
42 eig, eigSH, eigSH', 43 eig, eigSH, eigSH',
43 eigenvalues, eigenvaluesSH, eigenvaluesSH', 44 eigenvalues, eigenvaluesSH, eigenvaluesSH',
44-- ** QR 45-- ** QR
45 qr, 46 qr, rq,
46-- ** Cholesky 47-- ** Cholesky
47 chol, cholSH, 48 chol, cholSH,
48-- ** Hessenberg 49-- ** Hessenberg
@@ -91,6 +92,7 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
91 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t 92 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
92 linearSolve' :: Matrix t -> Matrix t -> Matrix t 93 linearSolve' :: Matrix t -> Matrix t -> Matrix t
93 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t 94 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
95 linearSolveLS' :: Matrix t -> Matrix t -> Matrix t
94 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 96 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
95 eigSH'' :: Matrix t -> (Vector Double, Matrix t) 97 eigSH'' :: Matrix t -> (Vector Double, Matrix t)
96 eigOnly :: Matrix t -> Vector (Complex Double) 98 eigOnly :: Matrix t -> Vector (Complex Double)
@@ -103,35 +105,134 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
103 multiply' :: Matrix t -> Matrix t -> Matrix t 105 multiply' :: Matrix t -> Matrix t -> Matrix t
104 106
105 107
106-- | Full singular value decomposition using lapack's dgesdd or zgesdd. 108instance Field Double where
109 svd' = svdRd
110 thinSVD' = thinSVDRd
111 sv' = svR
112 luPacked' = luR
113 luSolve' (l_u,perm) = lusR l_u perm
114 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
115 linearSolveLS' = linearSolveLSR
116 linearSolveSVD' = linearSolveSVDR Nothing
117 ctrans' = trans
118 eig' = eigR
119 eigSH'' = eigS
120 eigOnly = eigOnlyR
121 eigOnlySH = eigOnlyS
122 cholSH' = cholS
123 qr' = unpackQR . qrR
124 hess' = unpackHess hessR
125 schur' = schurR
126 multiply' = multiplyR
127
128instance Field (Complex Double) where
129 svd' = svdCd
130 thinSVD' = thinSVDCd
131 sv' = svC
132 luPacked' = luC
133 luSolve' (l_u,perm) = lusC l_u perm
134 linearSolve' = linearSolveC
135 linearSolveLS' = linearSolveLSC
136 linearSolveSVD' = linearSolveSVDC Nothing
137 ctrans' = conj . trans
138 eig' = eigC
139 eigOnly = eigOnlyC
140 eigSH'' = eigH
141 eigOnlySH = eigOnlyH
142 cholSH' = cholH
143 qr' = unpackQR . qrC
144 hess' = unpackHess hessC
145 schur' = schurC
146 multiply' = multiplyC
147
148--------------------------------------------------------------
149
150-- | Full singular value decomposition.
107svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 151svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
108svd = svd' 152svd = svd'
109 153
110-- | Thin singular value decomposition using lapack's dgesvd or zgesvd. 154-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
155--
156-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
111thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 157thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
112thinSVD = thinSVD' 158thinSVD = thinSVD'
113 159
114-- | Singular values using lapack's dgesvd or zgesvd. 160-- | Singular values only.
115sv :: Field t => Matrix t -> Vector Double 161singularValues :: Field t => Matrix t -> Vector Double
116sv = sv' 162singularValues = sv'
163
164-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
165--
166-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
167fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
168fullSVD m = (u,d,v) where
169 (u,s,v) = svd m
170 d = diagRect s r c
171 r = rows m
172 c = cols m
173
174-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
175compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
176compactSVD m = (u', subVector 0 d s, v') where
177 (u,s,v) = thinSVD m
178 d = rankSVD (1*eps) m s `max` 1
179 u' = takeColumns d u
180 v' = takeColumns d v
181
182vertical m = rows m >= cols m
183
184-- | Singular values and all right singular vectors.
185rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
186rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
187 | otherwise = let (_,s,v) = svd m in (s,v)
188
189-- | Singular values and all right singular vectors.
190leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
191leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
192 | otherwise = let (u,s,_) = thinSVD m in (u,s)
193
194
195{-# DEPRECATED full "use fullSVD instead" #-}
196full svdFun m = (u, d ,v) where
197 (u,s,v) = svdFun m
198 d = diagRect s r c
199 r = rows m
200 c = cols m
201
202{-# DEPRECATED economy "use compactSVD instead" #-}
203economy svdFun m = (u', subVector 0 d s, v') where
204 (u,s,v) = svdFun m
205 d = rankSVD (1*eps) m s `max` 1
206 u' = takeColumns d u
207 v' = takeColumns d v
208
209
210--------------------------------------------------------------
117 211
118-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. 212-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
119luPacked :: Field t => Matrix t -> (Matrix t, [Int]) 213luPacked :: Field t => Matrix t -> (Matrix t, [Int])
120luPacked = luPacked' 214luPacked = luPacked'
121 215
122-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization 216-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
123-- obtained by 'luPacked'.
124luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t 217luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
125luSolve = luSolve' 218luSolve = luSolve'
126 219
127-- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv. 220-- | 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'.
128-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system. 221-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
129-- See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
130linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t 222linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
131linearSolve = linearSolve' 223linearSolve = linearSolve'
224
225-- | 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.
132linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t 226linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
133linearSolveSVD = linearSolveSVD' 227linearSolveSVD = linearSolveSVD'
134 228
229
230-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
231linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
232linearSolveLS = linearSolveLS'
233
234--------------------------------------------------------------
235
135-- | Eigenvalues and eigenvectors of a general square matrix. 236-- | Eigenvalues and eigenvectors of a general square matrix.
136-- 237--
137-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@ 238-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
@@ -150,24 +251,45 @@ eigSH' = eigSH''
150eigenvaluesSH' :: Field t => Matrix t -> Vector Double 251eigenvaluesSH' :: Field t => Matrix t -> Vector Double
151eigenvaluesSH' = eigOnlySH 252eigenvaluesSH' = eigOnlySH
152 253
153-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric. 254-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
154cholSH :: Field t => Matrix t -> Matrix t 255--
155cholSH = cholSH' 256-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
257eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
258eigSH m | m `equal` ctrans m = eigSH' m
259 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
260
261-- | Eigenvalues of a complex hermitian or real symmetric matrix.
262eigenvaluesSH :: Field t => Matrix t -> Vector Double
263eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
264 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
265
266--------------------------------------------------------------
156 267
157-- | QR factorization using lapack's dgeqr2 or zgeqr2. 268-- | QR factorization.
158-- 269--
159-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. 270-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
160qr :: Field t => Matrix t -> (Matrix t, Matrix t) 271qr :: Field t => Matrix t -> (Matrix t, Matrix t)
161qr = qr' 272qr = qr'
162 273
163-- | Hessenberg factorization using lapack's dgehrd or zgehrd. 274-- | RQ factorization.
275--
276-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
277rq :: Field t => Matrix t -> (Matrix t, Matrix t)
278rq m = (r,q) where
279 (q',r') = qr $ trans $ rev1 m
280 r = rev2 (trans r')
281 q = rev2 (trans q')
282 rev1 = flipud . fliprl
283 rev2 = fliprl . flipud
284
285-- | Hessenberg factorization.
164-- 286--
165-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary 287-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
166-- and h is in upper Hessenberg form. 288-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
167hess :: Field t => Matrix t -> (Matrix t, Matrix t) 289hess :: Field t => Matrix t -> (Matrix t, Matrix t)
168hess = hess' 290hess = hess'
169 291
170-- | Schur factorization using lapack's dgees or zgees. 292-- | Schur factorization.
171-- 293--
172-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary 294-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
173-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is 295-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
@@ -187,137 +309,45 @@ multiply :: Field t => Matrix t -> Matrix t -> Matrix t
187multiply = multiply' 309multiply = multiply'
188 310
189 311
190instance Field Double where 312-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
191 svd' = svdRd 313cholSH :: Field t => Matrix t -> Matrix t
192 thinSVD' = thinSVDRd 314cholSH = cholSH'
193 sv' = svR
194 luPacked' = luR
195 luSolve' (l_u,perm) = lusR l_u perm
196 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
197 linearSolveSVD' = linearSolveSVDR Nothing
198 ctrans' = trans
199 eig' = eigR
200 eigSH'' = eigS
201 eigOnly = eigOnlyR
202 eigOnlySH = eigOnlyS
203 cholSH' = cholS
204 qr' = unpackQR . qrR
205 hess' = unpackHess hessR
206 schur' = schurR
207 multiply' = multiplyR
208
209instance Field (Complex Double) where
210 svd' = svdCd
211 thinSVD' = thinSVDCd
212 sv' = svC
213 luPacked' = luC
214 luSolve' (l_u,perm) = lusC l_u perm
215 linearSolve' = linearSolveC
216 linearSolveSVD' = linearSolveSVDC Nothing
217 ctrans' = conj . trans
218 eig' = eigC
219 eigOnly = eigOnlyC
220 eigSH'' = eigH
221 eigOnlySH = eigOnlyH
222 cholSH' = cholH
223 qr' = unpackQR . qrC
224 hess' = unpackHess hessC
225 schur' = schurC
226 multiply' = multiplyC
227
228
229-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
230--
231-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
232eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
233eigSH m | m `equal` ctrans m = eigSH' m
234 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
235
236-- | Eigenvalues of a complex hermitian or real symmetric matrix.
237eigenvaluesSH :: Field t => Matrix t -> Vector Double
238eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
239 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
240 315
241-- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. 316-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
242-- 317--
243-- If @c = chol m@ then @m == ctrans c \<> c@. 318-- If @c = chol m@ then @m == ctrans c \<> c@.
244chol :: Field t => Matrix t -> Matrix t 319chol :: Field t => Matrix t -> Matrix t
245chol m | m `equal` ctrans m = cholSH m 320chol m | m `equal` ctrans m = cholSH m
246 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" 321 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
247 322
323
324
325
248square m = rows m == cols m 326square m = rows m == cols m
249 327
250-- | determinant of a square matrix, computed from the LU decomposition. 328-- | Determinant of a square matrix.
251det :: Field t => Matrix t -> t 329det :: Field t => Matrix t -> t
252det m | square m = s * (product $ toList $ takeDiag $ lup) 330det m | square m = s * (product $ toList $ takeDiag $ lup)
253 | otherwise = error "det of nonsquare matrix" 331 | otherwise = error "det of nonsquare matrix"
254 where (lup,perm) = luPacked m 332 where (lup,perm) = luPacked m
255 s = signlp (rows m) perm 333 s = signlp (rows m) perm
256 334
257-- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf. 335-- | Explicit LU factorization of a general matrix.
258-- 336--
259-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular, 337-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
260-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation. 338-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
261lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) 339lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
262lu = luFact . luPacked 340lu = luFact . luPacked
263 341
264-- | Inverse of a square matrix using lapacks' dgesv and zgesv. 342-- | Inverse of a square matrix.
265inv :: Field t => Matrix t -> Matrix t 343inv :: Field t => Matrix t -> Matrix t
266inv m | square m = m `linearSolve` ident (rows m) 344inv m | square m = m `linearSolve` ident (rows m)
267 | otherwise = error "inv of nonsquare matrix" 345 | otherwise = error "inv of nonsquare matrix"
268 346
269-- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss. 347-- | Pseudoinverse of a general matrix.
270pinv :: Field t => Matrix t -> Matrix t 348pinv :: Field t => Matrix t -> Matrix t
271pinv m = linearSolveSVD m (ident (rows m)) 349pinv m = linearSolveSVD m (ident (rows m))
272 350
273
274{-# DEPRECATED full "use fullSVD instead" #-}
275full svdFun m = (u, d ,v) where
276 (u,s,v) = svdFun m
277 d = diagRect s r c
278 r = rows m
279 c = cols m
280
281-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
282--
283-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
284fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
285fullSVD m = (u,d,v) where
286 (u,s,v) = svd m
287 d = diagRect s r c
288 r = rows m
289 c = cols m
290
291
292{-# DEPRECATED economy "use compactSVD instead" #-}
293economy svdFun m = (u', subVector 0 d s, v') where
294 (u,s,v) = svdFun m
295 d = rankSVD (1*eps) m s `max` 1
296 u' = takeColumns d u
297 v' = takeColumns d v
298
299-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
300--
301-- If @(u,s,v) = compactSVD m@ then @m == u \<> diag s \<> trans v@.
302compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
303compactSVD m = (u', subVector 0 d s, v') where
304 (u,s,v) = thinSVD m
305 d = rankSVD (1*eps) m s `max` 1
306 u' = takeColumns d u
307 v' = takeColumns d v
308
309vertical m = rows m >= cols m
310
311-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
312rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
313rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
314 | otherwise = let (_,s,v) = svd m in (s,v)
315
316-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
317leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
318leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
319 | otherwise = let (u,s,_) = thinSVD m in (u,s)
320
321-- | Numeric rank of a matrix from the SVD decomposition. 351-- | Numeric rank of a matrix from the SVD decomposition.
322rankSVD :: Element t 352rankSVD :: Element t
323 => Double -- ^ numeric zero (e.g. 1*'eps') 353 => Double -- ^ numeric zero (e.g. 1*'eps')
@@ -379,12 +409,12 @@ pnormCV PNorm1 = norm1 . mapVector magnitude
379pnormCV Infinity = vectorMax . mapVector magnitude 409pnormCV Infinity = vectorMax . mapVector magnitude
380--pnormCV _ = error "pnormCV not yet defined" 410--pnormCV _ = error "pnormCV not yet defined"
381 411
382pnormRM PNorm2 m = sv m @> 0 412pnormRM PNorm2 m = singularValues m @> 0
383pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m 413pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
384pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) 414pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
385--pnormRM _ _ = error "p norm not yet defined" 415--pnormRM _ _ = error "p norm not yet defined"
386 416
387pnormCM PNorm2 m = sv m @> 0 417pnormCM PNorm2 m = singularValues m @> 0
388pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m 418pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m
389pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) 419pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)
390--pnormCM _ _ = error "p norm not yet defined" 420--pnormCM _ _ = error "p norm not yet defined"
@@ -433,12 +463,12 @@ nullspaceSVD hint a (s,v) = vs where
433 463
434-- | The nullspace of a matrix. See also 'nullspaceSVD'. 464-- | The nullspace of a matrix. See also 'nullspaceSVD'.
435nullspacePrec :: Field t 465nullspacePrec :: Field t
436 => Double -- ^ relative tolerance in 'eps' units 466 => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
437 -> Matrix t -- ^ input matrix 467 -> Matrix t -- ^ input matrix
438 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 468 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
439nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m) 469nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
440 470
441-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'. 471-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
442nullVector :: Field t => Matrix t -> Vector t 472nullVector :: Field t => Matrix t -> Vector t
443nullVector = last . nullspacePrec 1 473nullVector = last . nullspacePrec 1
444 474
@@ -522,11 +552,11 @@ uH (pq, tau) = (p,h)
522-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values. 552-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
523rcond :: Field t => Matrix t -> Double 553rcond :: Field t => Matrix t -> Double
524rcond m = last s / head s 554rcond m = last s / head s
525 where s = toList (sv m) 555 where s = toList (singularValues m)
526 556
527-- | Number of linearly independent rows or columns. 557-- | Number of linearly independent rows or columns.
528rank :: Field t => Matrix t -> Int 558rank :: Field t => Matrix t -> Int
529rank m = rankSVD eps m (sv m) 559rank m = rankSVD eps m (singularValues m)
530 560
531{- 561{-
532expm' m = case diagonalize (complex m) of 562expm' m = case diagonalize (complex m) of
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 7284cd9..8c55486 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -178,6 +178,11 @@ runTests n = do
178 putStrLn "------ luSolve" 178 putStrLn "------ luSolve"
179 test (linearSolveProp (luSolve.luPacked) . rSqWC) 179 test (linearSolveProp (luSolve.luPacked) . rSqWC)
180 test (linearSolveProp (luSolve.luPacked) . cSqWC) 180 test (linearSolveProp (luSolve.luPacked) . cSqWC)
181 putStrLn "------ luSolveLS"
182 test (linearSolveProp linearSolveLS . rSqWC)
183 test (linearSolveProp linearSolveLS . cSqWC)
184 test (linearSolveProp2 linearSolveLS . rConsist)
185 test (linearSolveProp2 linearSolveLS . cConsist)
181 putStrLn "------ pinv (linearSolveSVD)" 186 putStrLn "------ pinv (linearSolveSVD)"
182 test (pinvProp . rM) 187 test (pinvProp . rM)
183 test (pinvProp . cM) 188 test (pinvProp . cM)
@@ -220,6 +225,8 @@ runTests n = do
220 putStrLn "------ qr" 225 putStrLn "------ qr"
221 test (qrProp . rM) 226 test (qrProp . rM)
222 test (qrProp . cM) 227 test (qrProp . cM)
228 test (rqProp . rM)
229 test (rqProp . cM)
223 putStrLn "------ hess" 230 putStrLn "------ hess"
224 test (hessProp . rSq) 231 test (hessProp . rSq)
225 test (hessProp . cSq) 232 test (hessProp . cSq)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 7bb914f..405ce64 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -32,13 +32,13 @@ module Numeric.LinearAlgebra.Tests.Properties (
32 svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4, 32 svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,
33 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, 33 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
34 eigProp, eigSHProp, eigProp2, eigSHProp2, 34 eigProp, eigSHProp, eigProp2, eigSHProp2,
35 qrProp, 35 qrProp, rqProp,
36 hessProp, 36 hessProp,
37 schurProp1, schurProp2, 37 schurProp1, schurProp2,
38 cholProp, 38 cholProp,
39 expmDiagProp, 39 expmDiagProp,
40 multProp1, multProp2, 40 multProp1, multProp2,
41 linearSolveProp 41 linearSolveProp, linearSolveProp2
42) where 42) where
43 43
44import Numeric.LinearAlgebra 44import Numeric.LinearAlgebra
@@ -177,10 +177,11 @@ svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
177 (s',v') = rightSVC m 177 (s',v') = rightSVC m
178 (u',s'') = leftSVC m 178 (u',s'') = leftSVC m
179 179
180svdProp7 m = s |~| s' && u |~| u' && v |~| v' 180svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
181 where (u,s,v) = svd m 181 where (u,s,v) = svd m
182 (s',v') = rightSV m 182 (s',v') = rightSV m
183 (u',s'') = leftSV m 183 (u',s'') = leftSV m
184 s''' = singularValues m
184 185
185------------------------------------------------------------------ 186------------------------------------------------------------------
186 187
@@ -201,6 +202,16 @@ eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
201qrProp m = q <> r |~| m && unitary q && upperTriang r 202qrProp m = q <> r |~| m && unitary q && upperTriang r
202 where (q,r) = qr m 203 where (q,r) = qr m
203 204
205rqProp m = r <> q |~| m && unitary q && utr
206 where (r,q) = rq m
207 upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
208 where t = f-c
209 utr = upptr (rows r) (cols r) * r |~| r
210
211upperTriang' m = rows m == 1 || down |~| z
212 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
213 z = constant 0 (dim down)
214
204hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h 215hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
205 where (p,h) = hess m 216 where (p,h) = hess m
206 217
@@ -226,3 +237,8 @@ multProp1 (a,b) = a <> b |~| mulH a b
226multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a 237multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a
227 238
228linearSolveProp f m = f m m |~| ident (rows m) 239linearSolveProp f m = f m m |~| ident (rows m)
240
241linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
242 where q = min (rows a) (cols a)
243 b = a <> x
244 wc = rank a == q