Cholesky decomposition

[EnricaBelfiore]

@luca-fiorito-11 @AitorBengoechea

I and @GrimFe noticed that the last modification introduced in the sampling procedure does not work for all nuclear data covariance matrices. More in the details, we tested it for covariances of Pu239, U235, U233 and H1.
First of all, was this method already tested with real test cases? In case it was, which nuclides did you use for testing @AitorBengoechea?

We think that this problem is related to the fact that the covarinace for nuclear data are not always positive defined. As far as we understand, to_positive method handles this by transforming such matrices into SEMI-positive defined matrices. LU decomposition implemented inside sparse_table_cholesky method does not properly work for singular matrix (det=0).

We propose to back up to the previous decomposition scheme (cov = UDU.T) and to return L=U*sqrt(D), so that cov=L*L.T, where cov is the semipositive approximation of the nuclear data covariance matrix.
Unfortunately, scipy.linalg.ldl (implementing UDU.T decomposition) does not handle sparse matrix. Therefore our suggestion is to directly implement scipy.linalg.ldl decomposition in get_L. Further possible memory optimization can and should be considered.

[GrimFe]

Correct implementation of the covariance matrix decomposition is also needed to test the convergence of mean/std for the sampling according to a lognormal distribution.

[luca-fiorito-11]

@EnricaBelfiore ,this is very helpful! Can you report here a simple test where the current method fails.
@AitorBengoechea , can you implement Enrica's test in a PR and if needed go back to the previous QR decomposition?

[EnricaBelfiore]

This is a simple test that we did with the covariance of the hydrogen:

endf6 = sandy.get_endf6_file('jeff_33','xs', 10010)
Cov = sandy.XsCov.from_endf6(endf6)
Cov = sandy.CategoryCov(Cov)
Cov.sampling(20)
RuntimeError: Factor is exactly singular

[luca-fiorito-11]

I see three methods already developed.

non-sparse QR
non-sparse LU
sparse LU

Is sparsity needed?

How do we handle negative eigenvalues?
Aitor developed to_positive to make them zero.
Other two functions _reduce_size and _restore_size where developed to reduce the size of the covariance matrix removing the zero eigenvalues and then to restore it. These methods (now used in method invert) were previously used for the covariance decomposition. Should they be put back?

[EnricaBelfiore]

Are we considering the UDU.T decomposition in the LU decomposition?

The method to_positive should rather be called to_semipositive, in my opinion, and it should be working.
One could consider implementing the real Cholesky decompositon to CategoryCov.to_positive()._reduce_size(), which is positive definite, and to use _restore_size to get back to the full matrix. Otherwise, we can use the square root Cholesky decomposition (UDU.T) to decompose the semi-positive approximation of the covariance matrix.

Note: the method to_positive() is here considered to be the one implemented, which returns a semi-positive definite matrix.

The definition of a proper measure for the magnitude of the error introduced by the matrix approximation would be needed to compare the two proposed methods.

[AitorBengoechea]

It will not work. If you reduce size, do cholesky decomposition and then restore size. If you do L*L^T with the restore size matrix, it will not return the original covariance matrix

[AitorBengoechea]

The problem is that there are rows full of zeros, so the determinant is automatically zero.
The qr decomposition method automatically removed them, because it worked with eigenvalues and eigenvectors. But, cholesky uses the entire matrix.

[luca-fiorito-11]

Anyways, none of the considered methods actually use the cholesky decomposition

[luca-fiorito-11]

Can we agree that we need methods:

to convert the covariance matrix into a semi-positive definite one
to convert the semi-positive covariance e matrix into a positive covariance matrix
to decompose the positive definite covariance matrix and extract L (for this step all decomposition methods should be equivalent, right?!)
to restore the covariance size if it changed

Please, provide feedback

[EnricaBelfiore]

Ok we should try also with the QR factorization but not for the covariance matrix.
To be clear: we should perform this decomposition for a matrix that will be equal to eigenvector * sqrt(eigenvalues) (to be consistent with Cholesky) of the covariance matrix. In this way the R.T obained will be equal to the L matrix that satisfies the requested constraints. If you have the previous code of sandy, can you please check if the QR factorization was performed in this way or in a similar way? Otherwise it cannot converge to the Cholesky factorization for me. Let me know if it is correct

[luca-fiorito-11]

The code is in the link above

[luca-fiorito-11]

Sorry, you can find the QR in line 1983 of https://github.com/luca-fiorito-11/sandy/blob/master/sandy/formats/utils.py

[AitorBengoechea]

Be careful with this method, it loses precision as the problem scales to a larger size. In my case, in 600 X 600 matrices I could already see numerical fluctuations.

[EnricaBelfiore]

Thank you for the link! It is consistent with what I wrote in my previous comment and it solves my doubts.
I had found the explanation in the following link: https://scicomp.stackexchange.com/questions/34645/how-can-i-get-cholesky-decomposition-from-eigenvalue-decomposition
I hope it will be useful for your implementation @AitorBengoechea

[AitorBengoechea]

qr decomposition is not working for a real covariance matrix (Pu-241 compute with njoy for 241 groups with a covariance matrix 11809 X 11809). There is a lot of numerical fluctuation.

[luca-fiorito-11]

@AitorBengoechea , can you fix the implementation? Only one method. If sparsity is too difficult to implement we neglect it for the moment.
This PR will follow the one on eigenvalues by @EnricaBelfiore .

To make this clear, the purpose of the eig method is ro return the eigenvalues and eigenvectors. In that sense, the existing implementation is ok as long as the parameter k is taken into account (for instance as a kwarg). This implementation should in principle be independent from the decomposition.

The decomposition could call the eig method if needed. The eig method could be adjusted if the decomposition cannot handle the existing implementation. In that case the justification should be reported in the eig docstring.

[AitorBengoechea]

Thanks to @EnricaBelfiore and @GrimFe for the help.

For the moment, with the methods available in sandy.cov and the qr decomposition (introducing a threshold and sparse calculations), the results of the decomposition and reconstruction of the covariance matrix is equal to the original one (size 637 X 637)

[luca-fiorito-11]

But with what matrix?

[AitorBengoechea]

Pu-241 covariance matrix generated with njoy for 12 energy groups

[GrimFe]

Not having NJOY installed, I have been working with covariance matrices from endf6 files. Is this fine or is one forced to process covariance matrices in NJOY?

[luca-fiorito-11]

@EnricaBelfiore and @GrimFe , is this issue fully associated to the eigenvalues problem? If we change the eig method, will this still be an issue.

Can we clearly state what the problem is here?

[AitorBengoechea]

A real covariance matrix, for 241 groups of energy generated with njoy has a shape (11809 X 11809). Q matrix or B matrix taken from the decay data have 4000 X 4000 ... and inside category Cov you have the GLS method for the covariance. In this case, only taking fission yields, you are working with a diagonal matrix of 1000 X 1000.
The numpy and pandas algorithm works well for matrices smaller than 4000 X 4000. Also you dont have to take care of sparse matrices, sparse_tables functions are already implemented for doing multiplication and inversion of very big matrices.

[GrimFe]

Ok, thank you, I see your point.
My question was related to the fact that not all decompositions are implemented in scipy.linalg.sparse module, if I am not mistaken. And I am not sure switching from sparse matrices to dense ones helps code performances and readability. For this reason, this should also be considered when it comes to the identification of the best decomposition technique.

[luca-fiorito-11]

Let's not forget our goal. We need one method to decompose the covariance matrix. This method was the QR decomposition for dense matrices that was implemented and which is still available on the master branch. Effiency is problem dependent and can be considered later. Sparsity is preferable but not needed at the moment. Other methods such as to_positive etc. mean more maintenance. They will be kept only if needed. If the QR decomposition is still the preferable one, then we'll revert to what we had.

[luca-fiorito-11]

About the symmetry check. This check is welcome but it cannot become a burden. There is no reason for checking 14 digits. We could check less or none, or we could weaken the error to a warning. The point is that in no way the check should become a problem. Maybe we should just remove since in the development of SANDY we start from the principle that "the users of the API know what they are doing"....then they will not load a no symmetrical matrix to CategoryCov

[GrimFe]

@AitorBengoechea FYI: this (https://github.com/yig/PySPQR) might be relevant for sparse matrices-based optimization in case we revert to QR.

[AitorBengoechea]

@EnricaBelfiore , @luca-fiorito-11 and @GrimFe :

Decomposition problem solved. We are going to select the qr decomposition. My results are the following:

Problem size and description: Pu-241 xs computed with njoy for 241 energy groups with a covariance matrix 11809 X 11809

Eigenvalues and eigenvectors obtained in 8-10 minutes

A = V.dot(E), then the transpose and finally the qr decomposition: between 20s to 2 min

Reconstruction of the covariance matrix: between 3 min to 4 min

Results: comparison of the relative difference of the original covariance matrix and the reconstructed one:

• The first values (more accurate because they perform less operations due to the triangular matrix structure): highest relative difference of 4.31880e-06 %
• The last values ( less accurate because there are more operation in the calculations): highest relative difference of 1.19921e-02 %

[GrimFe]

@AitorBengoechea
I am not quite sure I understand what you mean by

A = V.dot(E), then the transpose and finally the qr decomposition

Is A the matrix to QR decompose in order to get L == R.T, where R comes from QR decomposition and L is such that Cov == LL.T?

[luca-fiorito-11]

• the matrix size being 11809, and considering 241 energy groups, means 11809/241=49 MT reactions. This is a large number of reactions, is it true? 2nd note, even if we consider that pu241 has covariance data for 49 MTs, most of the associated reactions should have a threshold, meaning that many variance terms should be zero. This could in principle be largely simplified. To be checked if these zero terms affect the efficiency of the eigenvalue calculation.
• 10 minutes for the eig calculation is a lot. I don't recall having such problems with the master branch
• to set up a reference, can we assume ~20000 as the largest matrix size that we need to handle at the moment?

[AitorBengoechea]

If i do using reduce size and restore size, the precision doesnt change, but the times are much lower than with all the matrix:
Matrix size reduction: from 11809 to 6283

Eigenvalues and eigenvectors obtained in 1-2 minutes

A = V.dot(E), then the transpose and finally the qr decomposition: approx 30s

Reconstruction of the covariance matrix: between 2min

[luca-fiorito-11]

FYI, have a look at the Wiki page about cholesky and its proof by QR decomposition

[luca-fiorito-11]

@AitorBengoechea , CategoryCov.eig method was updated with PR #139.
You can push your changes for the covariance decomposition.