Possible issue in getting covariance matrix eigenvalues

[GrimFe]

@luca-fiorito-11 @AitorBengoechea

When trying to identify the relevant eigenvalues of the Pu239 covariance matrix, @EnricaBelfiore and I noticed a possible issue with the latest implementation.
Only six eigenvalues (the larger ones) are now returned by eig method of CategoryCov.
This is of course not coherent with the matrix shape and should be clarified.
Even if in the reported case further eigenvalues are low (despite them not being zero), this truncation to the sixth eigenvalue happens also when considering U233 covariance matrix. There the seventh eigenvalue is of the same order of magnitude as the sixth, suggesting the truncation is not order of magnitude-based. Moreover, in our opinion all the eigenvalues should be returned to get the best possible matrix representation.

It is worth noticing that np.linalg.eig(cov.data) reports all the eigenvalues, meaning there is no lack of information, but rather an error in eigenvalue calculation. Browsing for this issue, we found this scipy/scipy#11198 reporting a possible error in scipy.sparse.linalg.eigs(), which might be a good starting point for debugging.

Example of what reported in the following - considering U233:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
tape = sandy.get_endf6_file("jeff_33", "xs", 922330)
Cov = sandy.XsCov.from_endf6(tape)
cov = sandy.CategoryCov(Cov)
len(cov.eig()[0])
len(pd.Series(np.linalg.eig(cov.data)[0]).sort_values(ascending=False))
len(pd.Series(scipy.linalg.eig(cov.data)[0]).sort_values(ascending=False))

EE = pd.Series(np.linalg.eig(cov.data)[0]).sort_values(ascending=False)
plt.plot(np.array(range(len(EE))), EE, 'o-')
plt.grid()

[AitorBengoechea]

in my case, even if the the eig algorith is different, i obtain the same results doing qr decomposition and scipy method consume less memory, because they save in the memory only the important information

[EnricaBelfiore]

Did you a test a matrix with more than six eigenvalues "important"? Is it defined a kind of threshold to truncate the eigenvalues?

[GrimFe]

@AitorBengoechea What do you mean by important information? Loosing part of the covariance matrix before even working with it seems to possibly pose some problems. If scipy.sparse keeps the whole matrix information, then the eigenvalues should be the same as the ones of the original covariance matrix. In my opinion, no matter the implementation scheme, CategoryCov.eig() should return the full set of eigenvalues, even if the convergence will most likely be checked for the first one only.

Moreover, this is not strictly related to covariance matrix decomposition: as you can see in the example I reported in my previous comment, no decomposition is involved there. Could you please run it and check you get the same differences I highlighted?

[luca-fiorito-11]

Please, focus the conversation on the eigenvalues problem, not the decomposition

[AitorBengoechea]

E, V = np.linalg.eig(cov_mat.data.values)
M = V.dot(np.diag(np.sqrt(E)))
Q, R = np.linalg.qr(M.T)
E_2, V_2 = np.linalg.eig(cov_mat.data.values)
M_2 = V_2.dot(np.diag(np.sqrt(E_2)))
Q_2, R_2 = scipy.linalg.qr(M_2.T)
(Q.real == Q_2.real).all(), (R_2.real == R.real).all() # (True, True)

[GrimFe]

@AitorBengoechea Could you please provide an example one can run?
Not all the variables you are using are defined in the script you shared.

[luca-fiorito-11]

This is not a good example nor a test, let's focus on the real eigenvalue problem.

[luca-fiorito-11]

The scipy docs indicate that the sparse eigs solver only returns k eigenvalues (with k=6) by default.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html

This is no bug to me

[GrimFe]

Good to hear, I missed it.
But still, why are we cutting to k=6 eigenvalues returned by CategoryCov.eig()? Shouldn't we return all of them?

[luca-fiorito-11]

Yes, this was not the case before and it came in when we recently introduced the sparse algorithm. This is why we do tests and apparently the implemented tests were not the right ones.

@EnricaBelfiore , can you submit a quick PR and take care of this. Look whether we could keep working with the sparse algorithm (what k do we use), or not.

[luca-fiorito-11]

PR in #139

[luca-fiorito-11]

Recap of in-person discussion with @EnricaBelfiore and @GrimFe .

The eig method of CategoryCov shall:

• return eigenvalues and eigenvector of the covariance matrix
• allow for a eigenvalue truncation proportional to the magnitude of the largest positive eigenvalue

The latter will be implemented with a kwarg called tolerance

[luca-fiorito-11]

Just googling it is possible to find proof that the eigs of a covariance/correlation matrix must be >=0.
Let's close the discussion on eigenvalues once and for all. We shall remove all negative eigenvalues and replace them with zero.

[luca-fiorito-11]

See here the analytical solutions of the eigenvalues of a 3x3 correlation matrix as a function of the correlation coefficients

[luca-fiorito-11]

Issue #136 closed with PR #139