Develop capability to load inverse covariance instead of covariance

[zatkins2]

After doing some profiling, I think the solution to #109 will not run sufficiently quickly. For ASO, the covariance matrix is expected to be O(80k x 80k), or about 50GB. On one tiger node (112 cores), direct inversion of such an array takes ~8.5 minutes, but this job resource allocation is probably not what will be used for running a chain (eg, a single task for running a chain typically has 1-10 cores). 112 cores for such a task may waste a lot of resources, but waiting O(2 hours) on 10 cores to invert the matrix before starting a chain is prohibitive and also wasteful if done repeatedly.

The motivation for passing the covariance to mflike is that cuts can be applied dynamically. However, in the vast majority of cases, two things are usually true: the cuts are repeated over many instances of mflike (usually, they are decided on as a group and then set as a standard for many runs) and/or the number of elements cut is much smaller than the total data vector (if restricted to TT/TE/EE to begin with, ie, we aren't also cutting B crosses). This suggests two solutions:

1. mflike should be able to load an inverse covariance that has been "pre-cut" for some standard cuts, along with the info necessary to apply the same cuts to the forward-modeled data vector. This assumes that the inverse covariance has also been precomputed, which may be challenging for a public user. This doesn't support any dynamic cuts to be applied, the cuts are fixed to whatever went into the precomputed inverse-covariance.
2. mflike can load the entire, not-cut inverse covariance matrix, and assuming the number of elements to be cut is small, use the Schur's complement to compute the inverse of the cut inverse covariance on the fly. This requires inverting the subsection of the entire inverse covariance that corresponds to the cut elements. If that subsection is small, then this computation should be fast. Of course, if we are always cutting B crosses, then this subsection won't be small, so somehow there should be a flag to signal that we are restricting ourselves to a universe without B crosses in the first place. Overall, I think this solution is better in the long term.

This should be accompanied by an investigation into the numerical stability of large matrix inversion, and the additional computation required for solution (2).

[sgiardie]

Thanks @zatkins2. If we go for option 1, I think the sacc file accepts only one covariance per file (@damonge correct me if I'm wrong), so we should be extra careful to read and use just the pre-cut inverse covariance throughout mflike. So far, we have a couple of points in which we are using self.cov (here and here), but I think it should be easy to solve that (computing the self.logp_const is not vital and the GaussianData class from SOLikeT accepts also inv_cov).

In the long run, we should also think if we want to make the code flexible to both options. I.e., if there should be a flag in the yaml that allows us to use approach 1 or 2. And of course the user should know which sacc file to use in each case. But of course we should first investigate how feasible each option is.

[zatkins2]

Thanks @zatkins2. If we go for option 1, I think the sacc file accepts only one covariance per file (@damonge correct me if I'm wrong), so we should be extra careful to read and use just the pre-cut inverse covariance throughout mflike.

Yes, for both 1 and 2, I opened a corresponding issue in sacc to track this, but I think we want to decide what we want/need on the mflike side first, and that will dictate what we need to ask of sacc.

So far, we have a couple of points in which we are using self.cov (here and here), but I think it should be easy to solve that (computing the self.logp_const is not vital and the GaussianData class from SOLikeT accepts also inv_cov).

Good point. Right, does the logp_const do anything? Why is it there at all?

Haven't thought at all about SOLikeT, it looks like cov gets used more there, also with MultiGaussianData. I wonder if there will ever be a case when we need to use the full likelihood (not CMB-only) when doing cross-covs with other data. I kind of think we shouldn't worry too much about SOLikeT unless there really is a use-case for using SOLikeT on the full likelihood instead of CMB-only.

In the long run, we should also think if we want to make the code flexible to both options. I.e., if there should be a flag in the yaml that allows us to use approach 1 or 2. And of course the user should know which sacc file to use in each case. But of course we should first investigate how feasible each option is.

Yeah let's see what works first and then can decide on design stuff like this. I do kind of think we should push for 2, I only suggest 1 in case it's annoying to implement quickly.

[sgiardie]

Yes, I agree pretty much on everything.

About the logp_const, I think we realized that is was more of a nuisance than a useful thing when computing the chi2 for each run in ACT dr6. @HTJense am I wrong?

About the usefulness of self.cov vs self.inv_cov for SOLikeT: @itrharrison or @ggalloni, do you think we really need mflike to pass self.cov to SOLikeT (vs just passing self.inv_cov)?

[HTJense]

I think we could leave logp_const but hide it after a flag that disables it by default. As far as I know, there are people who would like to keep it (since it normalizes the likelihood), so I think this would be a good compromise. You could also use the inv_cov for computing the logp_const, it would just incur a minus sign in the log-determinant factor.

[zatkins2]

Ok nice! I like the minus sign thing, we should just do that.

For @itrharrison and @ggalloni, I think the other question is whether you are anticipating a need to use SOLikeT on the full (multi-tube, all fg and systematics) likelihood, vs. were planning on practically always using it for the CMB-only (no tubes, fg and systematics marginzalized out) likelihood.