Computing eigen values using python different from Eigen Library in Cpp in non-positive components

Question

I am using

scipy.linalg import eigh

to compute Eigenvalues. Everything works correctly, however, while computing a non-positive eigenvalue.The results provided from Scipy is different from Eigen library in Cpp.

I tried the same array data on 2 Cpp and Python. It seems to be correct with non-negative results. But not with negative one.

'''python

    from scipy.linalg import eigh
    eigenVal, eigenVec = eigh(XtX)
    idx = eigenVal.argsort()[::-1]
    eigenVal = eigenVal[idx] / n  # Descending Sort

'''

'''Cpp

   // Compute Eigendecomposition:
   Eigen::MatrixXf XtX (r,r);
   if (m <= n
        XtX.template triangularView<Eigen::Lower>() = X * X.transpose();
   else
        XtX.template triangularView<Eigen::Lower>() = X.transpose() * X;
   Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> eig (XtX);
   // eigenvalues provide squared singular values:
   Eigen::VectorXf s = eig.eigenvalues();

   s=s.reverse().eval()/n;

'''

output in python: '''

array([  1.97322578e+05,   8.04158203e+03,   7.24907227e+03,           
    6.59655957e+03,   5.78095068e+03,   5.29064844e+03,
    4.74032520e+03,   4.46454346e+03,   4.23509229e+03,
    4.09968652e+03,   3.52356104e+03,   3.17308252e+03,
    3.14735107e+03,   2.88025342e+03,   2.80748145e+03,
    2.62388794e+03,   2.50723022e+03,   2.42145996e+03,
    2.11814282e+03,   2.03532385e+03,   1.79338770e+03,
    1.75364771e+03,   1.69037793e+03,   1.57554895e+03,
    1.45615894e+03,   1.33216003e+03,   1.26516211e+03,
    1.18433081e+03,   1.14839172e+03,   9.90758057e+02,
    8.59174011e+02,   8.20709900e+02,   7.86408569e+02,
    7.27681091e+02,   6.52041077e+02,   5.90487732e+02,
    5.26059814e+02,   4.99097687e+02,   4.49730896e+02,
    4.16578094e+02,   3.22004700e+02,   3.04409729e+02,
    2.79609161e+02,   2.15276245e+02,   1.84194916e+02,
    1.57032639e+02,   1.26647469e+02,   9.45600433e+01,
    7.59200439e+01,   5.26064911e+01,   4.35399475e+01,
    3.20723610e+01,   2.21866665e+01,   1.30586596e+01,
    3.81547445e-03,   1.16530398e-03,   1.07570493e-03,
    9.38778510e-04,   7.38311443e-04,   5.67958807e-04,
    5.35471656e-04,   5.30579535e-04,   5.17265638e-04,
    5.16666041e-04,   5.10115584e-04,   5.03800577e-04,
    5.00108406e-04,   4.70151805e-04,   4.33832814e-04,
    4.32742119e-04,   4.15480201e-04,   3.78245371e-04,
    3.75256117e-04,   3.71288595e-04,   3.53171228e-04,
    3.44453088e-04,   3.06655653e-04,   3.04728514e-04,
    2.60292058e-04,   2.57157430e-04,   2.02826763e-04,
    1.92918960e-04,   1.73933528e-04,  -9.30694770e-03], 
dtype=float32)$

output in C:

    197323     8041.61    7249.11    6596.59    5780.98    5290.66    
    4740.34    4464.56    4235.1    4099.7    3523.58    3173.09
    3147.36    2880.26    2807.48    2623.89    2507.23    2421.46
    2118.15    2035.33    1793.39    1753.65    1690.38    1575.55
    1456.16    1332.16    1265.16    1184.33    1148.39    990.759
    859.177    820.711    786.409    727.682    652.043    590.488
    526.061    499.099    449.731    416.578    322.005     304.41
    279.611    215.276    184.195    157.033    126.647    94.5601
    75.9213    52.6073    43.5399    32.0722    22.1875    13.0586
    0.0208003    0.00212419    0.00183017    0.00154398    
    0.00125332   0.00120334    0.00113175    0.000939068  
    0.000862687  0.000722841     0.000684395  0.000629895  
    0.000528611  0.000499585  0.0004829     0.000405769  0.000357388  
    0.000167705  8.70496e-05 -1.52459e-05     -0.000112141  
    -0.000141369  -0.000237786  -0.000504876  -0.000550637     
    -0.000704483  -0.000729621  -0.000882211  -0.00105727  -0.00215314

Well which one is correct? Can you try that? Also you might want to switch to doubles. Floats are usually not faster on CPUs anyway. — Quimby, May 22 '19 at 15:52
Certainly you are correct that there can only be one set of eigenvalues for a matrix. If I had to guess I would say numerical precision. The eigenvalues nearest 0 are the hardest to compute precisely, especially for a large matrix. — Owen, May 22 '19 at 15:52
Hi. @Quimbly, yes. Definitely will switch to double. To be honest, this is the first time I am working with Eigenvalues, and I am not sure which one is correct. Which way could I try? My matrix size is 84x84. — anhbeo, May 22 '19 at 19:07

score 0 · Answer 1 · answered May 22 '19 at 20:43

You are using single precision float. Their relative precision is about 1e-7, so all eigenvalues below 197323*1e-7 can be considered as zero. With that in mind, both versions can be considered numerically equivalent. For more accuracy, switch to double.

Computing eigen values using python different from Eigen Library in Cpp in non-positive components

1 Answers1