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.hs116
1 files changed, 77 insertions, 39 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index e0cc0a0..1544d7c 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -10,10 +10,9 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional 10Stability : provisional
11Portability : uses ffi 11Portability : uses ffi
12 12
13A generic interface for some common functions. Using it we can write higher level algorithms 13Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices.
14and testing properties both for real and complex matrices.
15 14
16In any case, the specific functions for particular base types can also be explicitly 15Specific functions for particular base types can also be explicitly
17imported from "Numeric.LinearAlgebra.LAPACK". 16imported from "Numeric.LinearAlgebra.LAPACK".
18 17
19-} 18-}
@@ -34,7 +33,11 @@ module Numeric.LinearAlgebra.Algorithms (
34-- * Matrix factorizations 33-- * Matrix factorizations
35-- ** Singular value decomposition 34-- ** Singular value decomposition
36 svd, 35 svd,
37 full, economy, --thin, 36 fullSVD,
37 thinSVD,
38 compactSVD,
39 sv,
40 leftSV, rightSV,
38-- ** Eigensystems 41-- ** Eigensystems
39 eig, eigSH, eigSH', 42 eig, eigSH, eigSH',
40-- ** QR 43-- ** QR
@@ -54,6 +57,7 @@ module Numeric.LinearAlgebra.Algorithms (
54-- * Nullspace 57-- * Nullspace
55 nullspacePrec, 58 nullspacePrec,
56 nullVector, 59 nullVector,
60 nullspaceSVD,
57-- * Norms 61-- * Norms
58 Normed(..), NormType(..), 62 Normed(..), NormType(..),
59-- * Misc 63-- * Misc
@@ -63,8 +67,8 @@ module Numeric.LinearAlgebra.Algorithms (
63 haussholder, 67 haussholder,
64 unpackQR, unpackHess, 68 unpackQR, unpackHess,
65 pinvTol, 69 pinvTol,
66 rankSVD, ranksv, 70 ranksv,
67 nullspaceSVD 71 full, economy
68) where 72) where
69 73
70 74
@@ -80,6 +84,8 @@ import Data.Array
80-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 84-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
81class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where 85class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
82 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 86 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
87 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
88 sv' :: Matrix t -> Vector Double
83 luPacked' :: Matrix t -> (Matrix t, [Int]) 89 luPacked' :: Matrix t -> (Matrix t, [Int])
84 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t 90 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
85 linearSolve' :: Matrix t -> Matrix t -> Matrix t 91 linearSolve' :: Matrix t -> Matrix t -> Matrix t
@@ -94,10 +100,18 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
94 multiply' :: Matrix t -> Matrix t -> Matrix t 100 multiply' :: Matrix t -> Matrix t -> Matrix t
95 101
96 102
97-- | Singular value decomposition using lapack's dgesvd or zgesvd. 103-- | Full singular value decomposition using lapack's dgesdd or zgesdd.
98svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 104svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
99svd = svd' 105svd = svd'
100 106
107-- | Thin singular value decomposition using lapack's dgesvd or zgesvd.
108thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
109thinSVD = thinSVD'
110
111-- | Singular values using lapack's dgesvd or zgesvd.
112sv :: Field t => Matrix t -> Vector Double
113sv = sv'
114
101-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'. 115-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
102luPacked :: Field t => Matrix t -> (Matrix t, [Int]) 116luPacked :: Field t => Matrix t -> (Matrix t, [Int])
103luPacked = luPacked' 117luPacked = luPacked'
@@ -163,7 +177,9 @@ multiply = multiply'
163 177
164 178
165instance Field Double where 179instance Field Double where
166 svd' = svdR 180 svd' = svdRd
181 thinSVD' = thinSVDRd
182 sv' = svR
167 luPacked' = luR 183 luPacked' = luR
168 luSolve' (l_u,perm) = lusR l_u perm 184 luSolve' (l_u,perm) = lusR l_u perm
169 linearSolve' = linearSolveR -- (luSolve . luPacked) ?? 185 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
@@ -178,7 +194,9 @@ instance Field Double where
178 multiply' = multiplyR 194 multiply' = multiplyR
179 195
180instance Field (Complex Double) where 196instance Field (Complex Double) where
181 svd' = svdC 197 svd' = svdCd
198 thinSVD' = thinSVDCd
199 sv' = svC
182 luPacked' = luC 200 luPacked' = luC
183 luSolve' (l_u,perm) = lusC l_u perm 201 luSolve' (l_u,perm) = lusC l_u perm
184 linearSolve' = linearSolveC 202 linearSolve' = linearSolveC
@@ -232,36 +250,61 @@ inv m | square m = m `linearSolve` ident (rows m)
232pinv :: Field t => Matrix t -> Matrix t 250pinv :: Field t => Matrix t -> Matrix t
233pinv m = linearSolveSVD m (ident (rows m)) 251pinv m = linearSolveSVD m (ident (rows m))
234 252
235-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. 253
236-- 254{-# DEPRECATED full "use fullSVD instead" #-}
237-- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@.
238full :: Element t
239 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
240full svdFun m = (u, d ,v) where 255full svdFun m = (u, d ,v) where
241 (u,s,v) = svdFun m 256 (u,s,v) = svdFun m
242 d = diagRect s r c 257 d = diagRect s r c
243 r = rows m 258 r = rows m
244 c = cols m 259 c = cols m
245 260
246-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations. 261-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
247-- 262--
248-- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. 263-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
249economy :: Element t 264fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
250 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) 265fullSVD m = (u,d,v) where
266 (u,s,v) = svd m
267 d = diagRect s r c
268 r = rows m
269 c = cols m
270
271
272{-# DEPRECATED economy "use compactSVD instead" #-}
251economy svdFun m = (u', subVector 0 d s, v') where 273economy svdFun m = (u', subVector 0 d s, v') where
252 r@(u,s,v) = svdFun m 274 (u,s,v) = svdFun m
253 d = rankSVD (1*eps) m r `max` 1 275 d = rankSVD (1*eps) m s `max` 1
276 u' = takeColumns d u
277 v' = takeColumns d v
278
279-- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
280--
281-- If @(u,s,v) = compactSVD m@ then @m == u \<> diag s \<> trans v@.
282compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
283compactSVD m = (u', subVector 0 d s, v') where
284 (u,s,v) = thinSVD m
285 d = rankSVD (1*eps) m s `max` 1
254 u' = takeColumns d u 286 u' = takeColumns d u
255 v' = takeColumns d v 287 v' = takeColumns d v
256 288
289vertical m = rows m >= cols m
290
291-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
292rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
293rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
294 | otherwise = let (_,s,v) = svd m in (s,v)
295
296-- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
297leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
298leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
299 | otherwise = let (u,s,_) = thinSVD m in (u,s)
257 300
258-- | Numeric rank of a matrix from the SVD decomposition. 301-- | Numeric rank of a matrix from the SVD decomposition.
259rankSVD :: Element t 302rankSVD :: Element t
260 => Double -- ^ numeric zero (e.g. 1*'eps') 303 => Double -- ^ numeric zero (e.g. 1*'eps')
261 -> Matrix t -- ^ input matrix m 304 -> Matrix t -- ^ input matrix m
262 -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m 305 -> Vector Double -- ^ 'sv' of m
263 -> Int -- ^ rank of m 306 -> Int -- ^ rank of m
264rankSVD teps m (_,s,_) = ranksv teps (max (rows m) (cols m)) (toList s) 307rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
265 308
266-- | Numeric rank of a matrix from its singular values. 309-- | Numeric rank of a matrix from its singular values.
267ranksv :: Double -- ^ numeric zero (e.g. 1*'eps') 310ranksv :: Double -- ^ numeric zero (e.g. 1*'eps')
@@ -316,12 +359,12 @@ pnormCV PNorm1 = norm1 . mapVector magnitude
316pnormCV Infinity = vectorMax . mapVector magnitude 359pnormCV Infinity = vectorMax . mapVector magnitude
317--pnormCV _ = error "pnormCV not yet defined" 360--pnormCV _ = error "pnormCV not yet defined"
318 361
319pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m 362pnormRM PNorm2 m = sv m @> 0
320pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m 363pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
321pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) 364pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
322--pnormRM _ _ = error "p norm not yet defined" 365--pnormRM _ _ = error "p norm not yet defined"
323 366
324pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC m 367pnormCM PNorm2 m = sv m @> 0
325pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m 368pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m
326pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) 369pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)
327--pnormCM _ _ = error "p norm not yet defined" 370--pnormCM _ _ = error "p norm not yet defined"
@@ -354,9 +397,9 @@ nullspaceSVD :: Field t
354 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'), 397 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
355 -- or Right \"theoretical\" matrix rank. 398 -- or Right \"theoretical\" matrix rank.
356 -> Matrix t -- ^ input matrix m 399 -> Matrix t -- ^ input matrix m
357 -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m 400 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
358 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 401 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
359nullspaceSVD hint a z@(_,_,v) = vs where 402nullspaceSVD hint a (s,v) = vs where
360 r = rows a 403 r = rows a
361 c = cols a 404 c = cols a
362 tol = case hint of 405 tol = case hint of
@@ -364,9 +407,8 @@ nullspaceSVD hint a z@(_,_,v) = vs where
364 _ -> eps 407 _ -> eps
365 k = case hint of 408 k = case hint of
366 Right t -> t 409 Right t -> t
367 _ -> rankSVD tol a z 410 _ -> rankSVD tol a s
368 nsol = c - k 411 vs = drop k $ toRows $ ctrans v
369 vs = drop k (toColumns v)
370 412
371 413
372-- | The nullspace of a matrix. See also 'nullspaceSVD'. 414-- | The nullspace of a matrix. See also 'nullspaceSVD'.
@@ -374,10 +416,7 @@ nullspacePrec :: Field t
374 => Double -- ^ relative tolerance in 'eps' units 416 => Double -- ^ relative tolerance in 'eps' units
375 -> Matrix t -- ^ input matrix 417 -> Matrix t -- ^ input matrix
376 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 418 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
377nullspacePrec t m = ns where 419nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
378 r@(_,_,v) = svd m
379 k = rankSVD (t*eps) m r
380 ns = drop k $ toRows $ ctrans v
381 420
382-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'. 421-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'.
383nullVector :: Field t => Matrix t -> Vector t 422nullVector :: Field t => Matrix t -> Vector t
@@ -405,7 +444,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
405-} 444-}
406--pinvTol :: Double -> Matrix Double -> Matrix Double 445--pinvTol :: Double -> Matrix Double -> Matrix Double
407pinvTol t m = v' `mXm` diag s' `mXm` trans u' where 446pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
408 (u,s,v) = svdR m 447 (u,s,v) = thinSVDRd m
409 sl@(g:_) = toList s 448 sl@(g:_) = toList s
410 s' = fromList . map rec $ sl 449 s' = fromList . map rec $ sl
411 rec x = if x < g*tol then 1 else 1/x 450 rec x = if x < g*tol then 1 else 1/x
@@ -460,15 +499,14 @@ uH (pq, tau) = (p,h)
460 499
461-------------------------------------------------------------------------- 500--------------------------------------------------------------------------
462 501
463-- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD. 502-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
464rcond :: Field t => Matrix t -> Double 503rcond :: Field t => Matrix t -> Double
465rcond m = last s / head s 504rcond m = last s / head s
466 where (_,s',_) = svd m 505 where s = toList (sv m)
467 s = toList s'
468 506
469-- | Number of linearly independent rows or columns. 507-- | Number of linearly independent rows or columns.
470rank :: Field t => Matrix t -> Int 508rank :: Field t => Matrix t -> Int
471rank m = rankSVD eps m (svd m) 509rank m = rankSVD eps m (sv m)
472 510
473{- 511{-
474expm' m = case diagonalize (complex m) of 512expm' m = case diagonalize (complex m) of