diff options
author | Maxim Koltsov <kolmax94@gmail.com> | 2018-11-13 15:39:55 +0300 |
---|---|---|
committer | Maxim Koltsov <kolmax94@gmail.com> | 2018-11-13 15:39:55 +0300 |
commit | cc737710b4878a387ea8090f9c2330b1e8d73f38 (patch) | |
tree | 02fc9ff241d06ac168e6a6cc19ecfd8e6ce8f812 /packages | |
parent | acb57ec8dd3011534813bdbee84563f1a830f499 (diff) |
Fix geig for complex eigenvalues, add tests
eigG was incorrectly returning eigenvectors corresponding to complex
eigenvalues. This was discovered with tests for geig, which are also
added to the repo.
Diffstat (limited to 'packages')
4 files changed, 40 insertions, 1 deletions
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs index b834b20..c9f5d0f 100644 --- a/packages/base/src/Internal/LAPACK.hs +++ b/packages/base/src/Internal/LAPACK.hs | |||
@@ -302,6 +302,14 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | |||
302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) | 302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) |
303 | fixeig _ _ = error "fixeig with impossible inputs" | 303 | fixeig _ _ = error "fixeig with impossible inputs" |
304 | 304 | ||
305 | -- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0. | ||
306 | -- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1) | ||
307 | -- can be different. Therefore old 'fixeig' would fail for 'eigG'. | ||
308 | fixeigG [] _ = [] | ||
309 | fixeigG [_] [v] = [comp' v] | ||
310 | fixeigG ((ar1:+ai1) : an : as) (v1:v2:vs) | ||
311 | | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs | ||
312 | | otherwise = comp' v1 : fixeigG (an:as) (v2:vs) | ||
305 | 313 | ||
306 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | 314 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. |
307 | -- The eigenvalues are not sorted. | 315 | -- The eigenvalues are not sorted. |
@@ -316,7 +324,7 @@ eigG a b = (alpha', beta, v'') | |||
316 | (alpha, beta, v) = eigGaux dggev a b "eigG" | 324 | (alpha, beta, v) = eigGaux dggev a b "eigG" |
317 | alpha' = fixeig1 alpha | 325 | alpha' = fixeig1 alpha |
318 | v' = toRows $ trans v | 326 | v' = toRows $ trans v |
319 | v'' = fromColumns $ fixeig (toList alpha') v' | 327 | v'' = fromColumns $ fixeigG (toList alpha') v' |
320 | 328 | ||
321 | eigGaux f ma mb st = unsafePerformIO $ do | 329 | eigGaux f ma mb st = unsafePerformIO $ do |
322 | a <- copy ColumnMajor ma | 330 | a <- copy ColumnMajor ma |
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 4ed1462..c0c151a 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs | |||
@@ -756,6 +756,9 @@ runTests n = do | |||
756 | test (eigSHProp2 . cHer) | 756 | test (eigSHProp2 . cHer) |
757 | test (eigProp2 . rSq) | 757 | test (eigProp2 . rSq) |
758 | test (eigProp2 . cSq) | 758 | test (eigProp2 . cSq) |
759 | putStrLn "------ geig" | ||
760 | test (uncurry geigProp . rSq2WC) | ||
761 | test (uncurry geigProp . cSq2WC) | ||
759 | putStrLn "------ nullSpace" | 762 | putStrLn "------ nullSpace" |
760 | test (nullspaceProp . rM) | 763 | test (nullspaceProp . rM) |
761 | test (nullspaceProp . cM) | 764 | test (nullspaceProp . cM) |
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 59230e0..97cfd01 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs | |||
@@ -17,6 +17,7 @@ Arbitrary instances for vectors, matrices. | |||
17 | 17 | ||
18 | module Numeric.LinearAlgebra.Tests.Instances( | 18 | module Numeric.LinearAlgebra.Tests.Instances( |
19 | Sq(..), rSq,cSq, | 19 | Sq(..), rSq,cSq, |
20 | Sq2WC(..), rSq2WC,cSq2WC, | ||
20 | Rot(..), rRot,cRot, | 21 | Rot(..), rRot,cRot, |
21 | rHer,cHer, | 22 | rHer,cHer, |
22 | WC(..), rWC,cWC, | 23 | WC(..), rWC,cWC, |
@@ -105,6 +106,23 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where | |||
105 | 106 | ||
106 | shrink (Sq a) = [ Sq b | b <- shrink a ] | 107 | shrink (Sq a) = [ Sq b | b <- shrink a ] |
107 | 108 | ||
109 | -- a pair of square matrices | ||
110 | newtype (Sq2WC a) = Sq2WC (Matrix a, Matrix a) deriving Show | ||
111 | instance (ArbitraryField a, Numeric a) => Arbitrary (Sq2WC a) where | ||
112 | arbitrary = do | ||
113 | n <- chooseDim | ||
114 | l <- vector (n*n) | ||
115 | r <- vector (n*n) | ||
116 | l' <- makeWC $ (n><n) l | ||
117 | r' <- makeWC $ (n><n) r | ||
118 | return $ Sq2WC (l', r') | ||
119 | where | ||
120 | makeWC m = do | ||
121 | let (u,_,v) = svd m | ||
122 | n = rows m | ||
123 | sv' <- replicateM n (choose (1,100)) | ||
124 | let s = diag (fromList sv') | ||
125 | return $ u <> real s <> tr v | ||
108 | 126 | ||
109 | -- a unitary matrix | 127 | -- a unitary matrix |
110 | newtype (Rot a) = Rot (Matrix a) deriving Show | 128 | newtype (Rot a) = Rot (Matrix a) deriving Show |
@@ -204,6 +222,9 @@ cRot (Rot m) = m :: CM | |||
204 | rSq (Sq m) = m :: RM | 222 | rSq (Sq m) = m :: RM |
205 | cSq (Sq m) = m :: CM | 223 | cSq (Sq m) = m :: CM |
206 | 224 | ||
225 | rSq2WC (Sq2WC (a, b)) = (a, b) :: (RM, RM) | ||
226 | cSq2WC (Sq2WC (a, b)) = (a, b) :: (CM, CM) | ||
227 | |||
207 | rWC (WC m) = m :: RM | 228 | rWC (WC m) = m :: RM |
208 | cWC (WC m) = m :: CM | 229 | cWC (WC m) = m :: CM |
209 | 230 | ||
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs index 6cd3a9e..38aa977 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs | |||
@@ -36,6 +36,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( | |||
36 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, | 36 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, |
37 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, | 37 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, |
38 | eigProp, eigSHProp, eigProp2, eigSHProp2, | 38 | eigProp, eigSHProp, eigProp2, eigSHProp2, |
39 | geigProp, | ||
39 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, | 40 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, |
40 | hessProp, | 41 | hessProp, |
41 | schurProp1, schurProp2, | 42 | schurProp1, schurProp2, |
@@ -237,6 +238,12 @@ eigProp2 m = fst (eig m) |~| eigenvalues m | |||
237 | 238 | ||
238 | eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m | 239 | eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m |
239 | 240 | ||
241 | geigProp a b = a' <> v <> diag betas' |~| b' <> v <> diag alphas | ||
242 | where (alphas, betas, v) = geig a b | ||
243 | betas' = complex betas | ||
244 | a' = complex a | ||
245 | b' = complex b | ||
246 | |||
240 | ------------------------------------------------------------------ | 247 | ------------------------------------------------------------------ |
241 | 248 | ||
242 | qrProp m = q <> r |~| m && unitary q && upperTriang r | 249 | qrProp m = q <> r |~| m && unitary q && upperTriang r |