Here is a plot of ellipticity residuals for three cases: ellipticity calculated from 1) HsmSourceMoments and HsmPsfMoments (s-p) 2) HsmSourceMoments and HsmPsfDebiasedMoments (s-pdb) 3) HsmPsfDebiasedMoments and HsmPsfMoments (pdb-p). `deepCoadd_meas` for tract=9697, all patches and only objects with `calib_psf_used=True` are used. The uncertainty of (pdb-p) and (s-p) are the same, as they should be and the uncertainty of (s-pdb) should be sqrt(2) times higher, which they are approximately. The mean values of (s-p)-(s-pdb) = mean (pdb-p), but a small difference exists because of some NaNs. However, that difference is much smaller than the statistical uncertainties.

The point to note is, the mean bias in (pdb-p) must match that of (s-p) so that (s-pdb) is debiased. However, mean(pdb-p) is about 2\sigma whereas mean(s-p) is about 10\sigma. The biases are way off. The bias in (pdb-p) and (s-p) have the same sign at least which reduces the bias in (s-pdb), but not nearly enough.

Does this give any clue? Here's a similar plot but for T = Ixx + Iyy. Note that there is no normalization by the zeroth moments done here - so it is fully linear, except for some centroiding bias. HsmSourceMoments have quite a few outliers, so not giving the mean and stdev, but from the plot, it is clear that the debiased moments aren't as wide as the source moments at all. The number of objects are similar (~35000). So it looks like the widths look fine when normalized by the zeroth moments, implicitly in the case of ellipticities, but not in raw moments themselves. This is different from the results found in the unit tests for debiased moments. I'm inclined to think that there is some assumption about the variance plane for coadds in the code computing PSF debiased moments that is incorrect.

I wrote a Python implementation to compute the debiased PSF moments and get identically small scatter for both the implementations. So I doubt there's a bug in the C++ implementation. The idea behind calculatinig the debiased moments is straightforward (as seen from the Python implementation), yet it is far from achieving the desired results!

The difference in the scatter between the debiased PSF moments and SourceMoments for stars are better on single-frame images than they are for coadds, but the SourceMoments still have error values that are almost 1.5 times larger. I think a factor of 1.2 difference is already known, but this is larger than that. It also appears a bit skewed, but that's later.

The trend observed in calexp, wherein the debiased PSF moments are less noisy compared to moments of stars is reproducible in controlled simulations. By remaking the simulations where the instFlux of the star is set to the measured value as opposed to the true value reproduces this trend, albeit with a SNR-dependent manner. This indicates that debiased PSF moments will require SNR-dependent corrections to remove the noise bias.