LAPACK dgeev difference from R

GregorySchwartz commented 5 years ago

The following matrix gives different results from their eigenvectors in hmatrix's eig and R's eigen(x, symmetric = FALSE).

x = (13 >< 13) [ 1.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 ,-0.5773502691896257 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 ,1.0 , -0.223606797749979 ,0.0 , 0.0 ,0.0 , -0.408248290463863 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 , -0. 
223606797749979 ,1.0 ,-0.1414213562373095 , 0.0 , -0.31622776601683794 , -0.18257418583505536 , -0.223606797749979 , -0.223606797749979 , -0.223606797749979 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 , 0.0 ,0.0 ,- 
0.1414213562373095 ,1.0 , -0.4472135954999579 ,0.0 ,0.0 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 ,0.0 ,-0.4472135954999579 , 1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 , -0. 
31622776601683794 ,0.0 , 0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , -0.5773502691896257 , -0.408248290463863 , -0.18257418583505536 ,0.0 , 0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 , -0.223606797749979 , -0.3162277660168379 
4 , 0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 
 ,0.0 ,0.0 , 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 , 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 , 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0. 
0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ] :: Matrix Double

They are close, but different. If they both use dgeev from LAPACK, why would there be any difference at all? While this difference may seem minor, if using eigSH and symmetric = TRUE for other matrices I can get flipped signs in some, but not all, elements of a vector, leading to downstream issues (I can understand if they are all flipped, but not a subset of the vector). I would expect identical results from identical LAPACK usage.

idontgetoutmuch commented 5 years ago

I'll try and take a look but my initial guess would be that R uses floats as default and hmatrix uses doubles - what size differences are you seeing - if it's about 10^-6 then I would feel more confident about my guess.

GregorySchwartz commented 5 years ago

Here is an example with inconsistent signs (also 13x13):

[ 0.13999163487187666, -0.5470825180776895, -0.116729147517097, -0.1650987696536079, -0.4012370645003707, 0.11443366282230742, -5.1247275079565946e-2, -0.24792033595207114, 0.13668153974474845, 0.16338403369120205, 0.32089983957652196, -0.4716352921999764, 0.17677669529663678, -6.1450695637302176e-2, -0.3868457584013466, -0.2246751539090412, 0.2334849191752439, 0.5674348983431925, -0.161833637955337, 7.247459145218677e-2, 0.3506123014915127, -0.1932968872330616, -0.2310599163213195, 2.7453604182822576e-2, -0.3334965133615023, , , 0.25, 0.5886721231798339, 6.590923086328981e-16, 0.41024659846700234, 3.0562234986292693e-17, 5.240220576747221e-17, 3.0757761033646486e-17, , , 0.0, -6.569361062107281e-17, -1.7021764939055836e-17, -6.189357470450551e-17, -0.41552720709623947, -1.9657491327809715e-16, 0.5590169943749473, 0.32678868630431923, 0.22594485526486643, -0.5940805855497716, 1.5375779902026008e-16, -2.3344392182498285e-16, -3.076977501618843e-16, -1.3877787807814457e-17, 4.1042451057582157e-16, 4.3796962362036487e-16, -2.7825760528471957e-16, 0.4082108622127935, 0.4078482888613904, 0.3952847075210477, -0.17812934066909622, -0.14290007364681268, 0.4116536373344582, -0.10861207338475154, -0.29970811010179715, -0.4113915283856835, -5.282389403377295e-2, 0.3865947523498372, -0.3207101215713721, 3.2583051578164844e-2, 0.39188987240798434, 0.2579459065208551, 0.17677669529663698, -0.22689604205941058, 5.113403605559701e-16, -0.2010095022693317, 0.19114669256674677, 8.318702434630204e-3, -0.3832474578106228, 9.113688168423065e-2, -0.19641117400815863, -0.13666807168493128, , 0.740439722187514, -0.2820745802051518, 1.2785076160344199e-16, 0.1767766952966369, -0.1989341133852227, 0.6700365082436466, 0.13048732910909594, 3.405564424780567e-17, -1.3652331978027063e-16, 3.890519714849777e-17, -6.938893903907228e-17, -4.2757640395395525e-16, -1.3181182814588419e-16, -4.767921531785965e-16, 0.2589200259181121, -0.5776329052891948, 0.30618621784789724, -0.2863961946800039, -0.1010456111077192, 0.1489478963217259, -1.975218980698884e-4, 0.12973198549782955, -0.19474536584512603, 0.2712422631875739, -3.798588578893797e-2, 0.8052747893670356, -7.44635477338513e-2, 7.765113779460535e-2, 0.18239529968020743, 0.2500000000000001, -0.28639619468000377, -0.10104561110771879, 0.14894789632172553, 9.227969471707231e-2, 7.663333403713142e-2, 0.14917163046120416, 0.25506170297389075, -0.6616167022550241, -0.38996043120524415, -0.30377808935812367, 7.76511377946055e-2, 0.18239529968020807, 0.25000000000000017, -0.2863961946800039, -0.10104561110771826, 0.14894789632172592, 6.151849439917506e-2, 0.21748595452420902, 0.6273692142723512, -0.45159969880154344, 0.1528750461285547, 3.823824535476732e-2, 0.33216224364663216, 7.76511377946054e-2, 0.18239529968020793, 0.2500000000000001, -0.22689604205941052, -8.508164251468192e-16, -0.20100950226933154, 0.21034616773994066, -0.3019179655680041, -0.25341246282630076, -0.6541112752338364, -0.10742645822606708, 9.576155101863087e-2, -0.38614398065679806, -0.28207458020515164, -1.622770545193915e-16, 0.17677669529663684, -0.22689604205941058, -4.439690316613866e-16, -0.2010095022693317, 0.20341474289691452, -0.5061703127908951, 0.35546692609855784, 0.4589032244362664, 0.38795095641429367, -4.9085470186553876e-2, -9.923266198428042e-2, -0.2820745802051518, -1.775687488265898e-16, 0.1767766952966369, -0.22689604205941039, -2.9741932511886096e-16, -0.2010095022693318, -0.8786184462419617, 9.88244013221014e-2, -1.576487102501008e-2, , , 0.0, 5.456109221769827e-2, -9.403659097376879e-2, -5.9095994277069575e-2, -0.2820745802051516, -4.906249560798318e-17, 0.1767766952966369 ]

Haskell's second smallest eigenvector is:

    -0.4716352921999764
    -0.3334965133615023
-1.9657491327809715e-16
     0.4078482888613904
     0.2579459065208551
 1.2785076160344199e-16
    -0.5776329052891948
    0.18239529968020743
    0.18239529968020807
    0.18239529968020793
 -1.622770545193915e-16
 -1.775687488265898e-16
 -4.906249560798318e-17

R's is:

  4.716353e-01 
  3.334965e-01 
 -4.440892e-16 
 -4.078483e-01 
 -2.579459e-01 
 -3.053113e-16 
  5.776329e-01 
 -1.823953e-01 
 -1.823953e-01 
 -1.823953e-01 
 -4.440892e-16 
 -4.440892e-16 
 -4.996004e-16

As you can see, the signs are switched inconsistently (which is a problem when separating by sign).

GregorySchwartz commented 5 years ago

I believe that was with eigSH.

idontgetoutmuch commented 5 years ago

I am only getting differences at roughly the double precision limit at the moment

{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main
    ( main
    ) where

import Data.List
import Numeric.LinearAlgebra
import qualified Language.R as R
import Language.R.QQ

x = (13 >< 13) [ 1.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 ,-0.5773502691896257 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,1.0 , -0.223606797749979 ,0.0 , 0.0 ,0.0 , -0.408248290463863 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 , -0.223606797749979 ,1.0 ,-0.1414213562373095 , 0.0 , -0.31622776601683794 , -0.18257418583505536 , -0.223606797749979 , -0.223606797749979 , -0.223606797749979 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,
                 0.0 ,0.0 ,-0.1414213562373095 ,1.0 , -0.4472135954999579 ,0.0 ,0.0 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 ,0.0 ,-0.4472135954999579 , 1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 -0.5773502691896257 , -0.408248290463863 , -0.18257418583505536 ,0.0 , 0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ] :: Matrix Double

y :: [Double]
y = concat $ toLists x

v :: R.R s [Double]
v = fmap R.dynSEXP $ do
  [r| eigen(matrix(y_hs, nrow = 13, ncol = 13), symmetric=FALSE)$values  |]

main = do
  u <- R.runRegion v
  print u
  let w = reverse $ sort $ toList $ cmap realPart $ fst $ eig x
  print w
  print $ maximum $ map abs $ zipWith (-) u w
  return ()

*Main> main
[1.8204394785379518,1.707106781186546,1.645399167126937,1.0000000000000004,1.0000000000000002,1.0,1.0,0.9999999999999999,0.9999999999999998,0.9999999999999996,0.5341613543351121,0.2928932188134526,-2.7755575615628914e-16]
[1.8204394785379523,1.7071067811865483,1.6453991671269381,1.0000000000000004,1.0000000000000002,1.0000000000000002,1.0000000000000002,1.0,0.9999999999999998,0.9999999999999994,0.5341613543351111,0.2928932188134528,5.551115123125783e-17]
2.4424906541753444e-15

idontgetoutmuch commented 5 years ago

{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main
    ( main
    ) where

import Data.Ord
import Data.List
import Numeric.LinearAlgebra
import qualified Language.R as R
import Language.R.QQ

bigA = (13 >< 13) [ 1.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 ,-0.5773502691896257 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,1.0 , -0.223606797749979 ,0.0 , 0.0 ,0.0 , -0.408248290463863 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 , -0.223606797749979 ,1.0 ,-0.1414213562373095 , 0.0 , -0.31622776601683794 , -0.18257418583505536 , -0.223606797749979 , -0.223606797749979 , -0.223606797749979 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,
                 0.0 ,0.0 ,-0.1414213562373095 ,1.0 , -0.4472135954999579 ,0.0 ,0.0 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 ,0.0 ,-0.4472135954999579 , 1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 -0.5773502691896257 , -0.408248290463863 , -0.18257418583505536 ,0.0 , 0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ] :: Matrix Double

y :: [Double]
y = concat $ toLists bigA

v :: R.R s [Double]
v = fmap R.dynSEXP $ do
  [r| eigen(matrix(y_hs, nrow = 13, ncol = 13), symmetric=FALSE)$values  |]

w'' :: R.R s [Double]
w'' = fmap R.dynSEXP $ do
  [r| eigen(matrix(y_hs, nrow = 13, ncol = 13), symmetric=FALSE)$vectors  |]

s :: Vector (Complex Double)
w :: Matrix (Complex Double)
(s, w) = eig bigA

eigValsVecs :: [(Complex Double, Vector (Complex Double))]
eigValsVecs = reverse $ sortBy (comparing (magnitude . fst)) (zip (toList s) (toColumns w))

s' :: Vector (Complex Double)
s' = fromList $ map fst eigValsVecs

w' :: Matrix (Complex Double)
w' = fromColumns $ map snd eigValsVecs

foo = (cmap (:+ 0.0) bigA) <> w
bar = w <> diag s

foo' = (cmap (:+ 0.0) bigA) <> w'
bar' = w' <> diag s'

main = do
  u <- R.runRegion v
  putStrLn "\nMaximum eigenvalue discrepancy"
  print $ maximum $ map abs $ zipWith (-) u (map realPart $ toList s')
  u' <- R.runRegion w''
  putStrLn "\nMaximum discrepancy per eigenvector"
  mapM_ putStrLn $ map show $ maximum $ map (map magnitude) $ map toList $ toColumns $
                   (cmap (:+ 0.0) $ (13 >< 13) u') - (tr w')
  return ()

*Main> main

Maximum eigenvalue discrepancy
2.4424906541753444e-15

Maximum discrepancy per eigenvector
1.1773442463596682
4.05351182757118e-15
0.8204931969340032
1.0151707304829045e-15
2.496342145796157e-16
2.1402506898478911e-16
1.5470182145299065e-16
2.7120135468264146e-17
2.5320434873859485e-16
2.3651256266174785e-16
0.8310544141924785
3.39491484346944e-16
4.440892098500626e-16
*Main>

idontgetoutmuch commented 5 years ago

I am going to leave my notes here but here is the problem as R tells us:

> solve(x)
Error in solve.default(x) : 
  system is computationally singular: reciprocal condition number = 2.14591e-18

Here are some notes:

{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main
    ( main
    ) where

import Data.Ord
import Data.List
import Numeric.LinearAlgebra
import qualified Language.R as R
import Language.R.QQ

bigA = (13 >< 13) [ 1.0 ,0.0 ,0.0 ,0.0 , 0.0 ,0.0 ,-0.5773502691896257 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,1.0 , -0.223606797749979 ,0.0 , 0.0 ,0.0 , -0.408248290463863 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 , -0.223606797749979 ,1.0 ,-0.1414213562373095 , 0.0 , -0.31622776601683794 , -0.18257418583505536 , -0.223606797749979 , -0.223606797749979 , -0.223606797749979 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,
                 0.0 ,0.0 ,-0.1414213562373095 ,1.0 , -0.4472135954999579 ,0.0 ,0.0 , -0.31622776601683794 , -0.31622776601683794 , -0.31622776601683794 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 ,0.0 ,-0.4472135954999579 , 1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 -0.5773502691896257 , -0.408248290463863 , -0.18257418583505536 ,0.0 , 0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.223606797749979 , -0.31622776601683794 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ,0.0 ,
                 0.0 ,0.0 , -0.31622776601683794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,1.0 ] :: Matrix Double

y :: [Double]
y = concat $ toLists bigA

v :: R.R s [Double]
v = fmap R.dynSEXP $ do
  [r| eigen(matrix(y_hs, nrow = 13, ncol = 13), symmetric=FALSE)$values  |]

w'' :: R.R s [Double]
w'' = fmap R.dynSEXP $ do
  [r| eigen(matrix(y_hs, nrow = 13, ncol = 13), symmetric=FALSE)$vectors  |]

s :: Vector (Complex Double)
w :: Matrix (Complex Double)
(s, w) = eig bigA

eigValsVecs :: [(Complex Double, Vector (Complex Double))]
eigValsVecs = reverse $ sortBy (comparing (magnitude . fst)) (zip (toList s) (toColumns w))

s' :: Vector (Complex Double)
s' = fromList $ map fst eigValsVecs

w' :: Matrix (Complex Double)
w' = fromColumns $ map snd eigValsVecs

foo = (cmap (:+ 0.0) bigA) <> w
bar = w <> diag s

foo' = (cmap (:+ 0.0) bigA) <> w'
bar' = w' <> diag s'

main = do
  u <- R.runRegion v
  putStrLn "\nMaximum eigenvalue discrepancy"
  print $ maximum $ map abs $ zipWith (-) u (map realPart $ toList s')
  u' <- R.runRegion w''
  putStrLn "\nMaximum discrepancy per eigenvector"
  mapM_ putStrLn $ map show $ maximum $ map (map magnitude) $ map toList $ toRows $
                   (tr $ cmap (:+ 0.0) $ (13 >< 13) u') - w'
  let foo'' = (tr $ (13 >< 13) u') <> diag (fromList u)
      bar'' = bigA <> (tr $ (13 >< 13) u')
  putStrLn "\nEigenvalue equation test"
  print $ maximum $ concat $ toLists $ cmap abs $ foo'' - bar''
  print $ maximum $ concat $ toLists $ cmap magnitude $ foo' - bar'
  print $ maximum $ concat $ toLists $ cmap magnitude $ (cmap (:+ 0.0) foo'') - bar'
  return ()

GregorySchwartz commented 5 years ago

Why would a singular matrix be the issue? It's still the same algorithm and library being used, so they should have identical results in both the eigenvalues and eigenvectors.

idontgetoutmuch commented 5 years ago

It's almost singular - so small changes to an input to a calculation using such a matrix will lead to large changes in the output. The R and Haskell libraries will marshal values slightly differently hence the results we see. I don't mean to be presumptuous but if it were me, I would investigate why I have an almost singular matrix. Once you have eliminated such matrices, the differences between R and Haskell should be at floating point precision (give or take). More info here: https://scicomp.stackexchange.com/questions/31016/condition-number-of-a-matrix

GregorySchwartz commented 5 years ago

They are mostly singular or near singular matrices because they are Laplacians. Finding the eigenvectors of Laplacians is a common procedure, so I was surprised to see the discrepancies.

idontgetoutmuch commented 5 years ago

It doesn't matter where the matrix comes from, if it is ill-conditioned then it is probably telling you there is something wrong with your model (something like non-identifiability, multi-collinearity?). Without knowing what you are trying to do, it's a bit hard to say more. But I don't think there is anything wrong with hmatrix. Maybe you should ask a question on https://scicomp.stackexchange.com? If you cross post it here, I can see if I can help.

haskell-numerics / hmatrix

LAPACK dgeev difference from R #288