Totally, yeah... In that case I only have one more idea. Estimate mean/var for each feature separately using all their observations, then standardize them before estimating their correlation with masking. Then multiplying stuff to get the covariance may have a slightly larger effective sample size?

One option is to take masked averages of products of pairs of features for cases where both are observed. But the resulting matrix is not positive definite, although you can project it... A more expensive way is EM in a multivariate normal, if you’re ok with normals. If not, maybe another model?