negative eigenvalues for weak lensing gaussian covariance

[gouyoub]

Hello,
I am trying to use pymaster to compute the partial-sky Gaussian covariance matrix of the decoupled bandpowers for 3x2pt. While I have no issue with the covariance of photometric galaxy clustering bandpowers, the one for weak lensing or galaxy-galaxy lensing is not positive definite.

Here is an example code I am using to compute the Gaussian covariance for weak lensing. I am using the latest version of pymaster (v2.2.2)

# Loading my theory Cl's and footprint
data_ref = '/home/hidra2/gouyou/euclid/nl_bias_flagship/data/'
ref = '/home/hidra2/gouyou/cosmosis-standard-library/working_dir/output/fs2_nlbias/simulation/fs2_firstchain_linbias_unbinned/'
nzbins = 3
cls_wl, keys_wl = load.get_cosmosis_Cl(ref, nzbins, 'WL', True)
mask = hp.read_map(data_ref+'mask/flagship2_mask_binary_NS512.fits')

# Initializing Nmt objects
nside = 512
fmask = nmt.NmtField(mask, None, spin=0)
b = nmt.NmtBin.from_nside_linear(nside, 20)
w = nmt.NmtWorkspace.from_fields(fmask, fmask, b)
cw = nmt.NmtCovarianceWorkspace.from_fields(fmask, fmask, fmask, fmask)

# Initializing the covariance array
nell = b.get_effective_ells().size
ncl = len(keys_wl)
print(ncl*nell)
covmat = np.zeros((ncl*nell, ncl*nell))

# Assigning symmetric keys to Cl dictionnary
for key in keys_wl:
    probeA, probeB = key.split('-')
    cls_wl['-'.join([probeB, probeA])] = cls_wl[key]

# Computation of all the blocks of the covariance matrix
for (idx1, key1), (idx2, key2) in itertools.combinations_with_replacement(enumerate(keys_wl), 2):
    print(key1, key2, flush=True)
    probeA, probeB = key1.split('-')
    probeC, probeD = key2.split('-')

    covmat[idx1*nell:(idx1+1)*nell, idx2*nell:(idx2+1)*nell] =\
        nmt.gaussian_covariance(cw, 0, 0, 0, 0,
                                [cls_wl['-'.join([probeA, probeC])]],
                                [cls_wl['-'.join([probeB, probeC])]],
                                [cls_wl['-'.join([probeA, probeD])]],
                                [cls_wl['-'.join([probeB, probeD])]],
                                w, wb=w)

covmat = covmat + covmat.T - np.diag(covmat.diagonal())

All the individual blocks of the covariance (corresponding to a given pair of Cl's) are positive definite but when considering the entire matrix it is not. To test this I use np.linalg.cholesky(covmat).

Note that I am only computing the E-mode covariance and considering it is a spin-0 field. I also tried with spin-2 fields with a zero contribution from B-modes and the result is the same.

Have you ever encountered this issue ?

I made many tests (varying probes, mask, resolution, binning, lmin, lmax and more) so I can give more details if needed.

Thank you in advance !

Sylvain GB.