gh-293: ensure delete-2 JK remains posdef #292
Conversation
Codecov ReportAll modified and coverable lines are covered by tests ✅
Flags with carried forward coverage won't be shown. Click here to find out more.
🚀 New features to boost your workflow:
|
| # Diagonalise the correction | ||
| for key in Q: | ||
| q = Q[key] | ||
| *_, length = q.shape | ||
| q_diag = np.diagonal(q, axis1=-2, axis2=-1) | ||
| q_diag_exp = np.zeros_like(q) | ||
| diag_indices = np.arange(length) # Indices for the diagonal | ||
| q_diag_exp[..., diag_indices, diag_indices] = q_diag | ||
| Q[key] = Result(q_diag_exp, axis=q.axis, ell=q.ell) |
There was a problem hiding this comment.
Why is this no longer necessary here? Doesn't the change also change the result of this function?
There was a problem hiding this comment.
Indeed it does but I think this way it is better. Rather than diagonalizing the correction and then simply adding it to the cov_jk, we can save the fully correction (in case somebody wishes to have a look) and then substract its diagonal to the diagonal of the cov_jk.
| ell=cov_jk[key].ell, | ||
| axis=cov_jk[key].axis, | ||
| ) | ||
| debiased_cov = deepcopy(cov_jk) |
There was a problem hiding this comment.
I don't think the deepcopy is useful here, if you set all entries to a new Result down below anyway.
There was a problem hiding this comment.
I think otherwise the entries of cov_jk get accidentally overwritten in the process.
| # abs is only needed when too few Jackknife regions are used | ||
| c[..., diag_indices, diag_indices] = abs(c_diag - q_diag) |
There was a problem hiding this comment.
Isn't this materially changing the output of the debiasing step? In the sense that the absolute value is something we are making up on the spot to fix a "problem" that is really intrinsic to the process?
There was a problem hiding this comment.
for large enough jk regions, q_diag should always be smaller than c_diag. The goal here is to cover for when only a few regions are used. We could something more kosher such as enfocing pos-def for each entry in the matrix stack but this would require more involved coding.
| _debiased_cov[key] = cov_jk[key].array - Q[key] | ||
| q = Q[key].array | ||
| c = _debiased_cov[key].array | ||
| *_, length = q.shape | ||
| q_diag = np.diagonal(q, axis1=-2, axis2=-1) | ||
| c_diag = np.diagonal(c, axis1=-2, axis2=-1) | ||
| # Indices for the diagonal | ||
| diag_indices = np.arange(length) | ||
| # abs is only needed when too few Jackknife regions are used | ||
| c[..., diag_indices, diag_indices] = abs(c_diag - q_diag) | ||
| _debiased_cov[key] = c |
There was a problem hiding this comment.
Isn't this test the exact same code used in the function? If so, can we have a test where you construct a known output for a known input by hand? Otherwise, we might just end up with the same bug in 2 places...
There was a problem hiding this comment.
So was before to be honest
Ensures that the delete2 covariance remains positive definite even when using a few jackknife regions.
Closes: #293