Cholesky decomposition

[luca-fiorito-11]

Discussed in #133

^{Originally posted by EnricaBelfiore February 18, 2022}
@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.

[luca-fiorito-11]

Closed with PR #143