summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Algorithms.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Algorithms.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs284
1 files changed, 157 insertions, 127 deletions
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