I was using the sersic fits in the object table for DP2 and noticed the fits looked off when projected over an image. I took the sersic parameters and fed them into AstroPhot (another parametric fitting code) and the residuals looked pretty bad. After some back and forth with claude we pinned the issue down to the sersic_reff_major being too small. As far as it could tell the reff was too small by a factor of 2ln(2) which it traced to the knot tables and could be explained by an accidental Gaussian sigma/fwhm switch somewhere in the knot table generation. Here is a plot showing for a few random objects, what I get if I naively take the values from the object table, vs using the rescaled reff, and then another column for reproducing the multi-Gaussian expansion.
I didn’t realize that the Rubin Sersic models were actually Gaussian mixture models until finding this issue. The mixture seems remarkably close to an actual sersic, however it breaks down at high radii and for high sersic index. I guess this isn’t a big issue for most sources, but it did result in some surprises for me. I was using the fraction of light outside a certain radius as a metric and if I used an actual sersic model with the fitted sersic index I got unrealistically high values when the sersic index was high. I should also note that even with just 6 random galaxies, there is a visible difference in a real sersic model vs the gaussian mixture, perhaps Rubin should consider fitting with true sersic models?
I wasn’t really able to trace through the pipeline myself since its all in C, but claude was able to follow it and it highlighted the lines of code that are relevant to what it found. If you’re interested, here is the write up it made to help point you in the right direction for the bug.
Yes, unfortunately I made an error transcribing the tables for the Sersic mixture some time ago and multiplied by sqrt(2ln(2)) instead of dividing. I’m not sure why I didn’t notice this earlier (besides the unit test not checking the size values), because I have done comparisons between e.g. exponential_reff_major and i_cModel_exp_reff_major before and you can clearly see the difference there.
I’ve filed two tickets (DM-56611 and DM-56612) to fix the issue. Unfortunately, I don’t think we’ll be able to correct the DP2 tables in place so in the meantime please multiply [exponential/sersic]_reff_[major/minor/x/y] by 2ln(2) to get the correct values.
This is partly due to the nature of Gaussian mixtures and partly by design. See DMTN-312 and the notebook where I originally derived the weights for more information.
In brief, the weights were fit to reproduce a 1D Sersic profile weighted by surface area (so essentially the 2D surface brightness) out to 12R_eff, and excluding an increasing large fraction of the core for n>1. For now we are using the 4-Gaussian approximation rather than the 8-Gaussian one; we’ll probably switch to the latter for DR1 since this choice was mainly cautionary to avoid uncommon performance issues that have since been resolved. The consequence is that the profiles deviate most from a true Sersic at 4-5R_eff. They also rapidly drop to zero at 8R_eff, but this is unavoidable - adding more Gaussians just pushes this radius further out.
This is precisely because the Gaussian mixture doesn’t and can’t reproduce the outer profile for large Sersic n. If you fit the Sersic index with AstroPhot you’ll probably get a lower value, and the fraction of light outside a given radius (assuming you’re looking at a value larger than 2-3Reff) should agree better. Having said that, the fact that the Sersic profile has such a large fraction of light at large radii for n>6 is practically problematic (if not pathological) and part of the reason we didn’t fit the mixture beyond 12Reff or use the more accurate 8-component mixture by default.
Also, the effective radius for the Gaussian mixture Sersic profile isn’t and can’t be precisely the half-light radius except for n=0.5. That’s part of the reason I named the columns reff rather than rhalf or r50. If you want to compute a radius for a given integrated flux percentile you’ll need to use the individual Gaussian weights. I should probably provide a convenience function for that.
I’m not sure what you mean - the third and fourth columns in your plot look virtually identical except for the one n>4 object. If you run AstroPhot with similar methodology - a simultaneous fit on all bands - you should get more or less the same answer for most n<2 galaxies, modulo choices in fitting area, PSF model (are you using the cell PSF directly or the double Gaussian mixture approximation thereof?) and priors. I’m happy to discuss this on Slack or in another thread if you have a particular project in mind.
Hi @dtaranu , thanks for the very detailed answer! Indeed no matter how much unit testing you do, you can only limit towards zero bugs rather than eliminate them entirely.
This could cause a lot of silent bugs for downstream users. It if can’t be fixed in the tables, then maybe a note could be included in the schemas?
I agree, but this is something of an expectation for people working with a sersic profile. Breaking that expectation without a warning can cause silent errors (like what happened to me before I started tracking things down). I think that somewhere prominent (like the schema table) it should be made clear that the sersic parameters are actually a (very close) approximation of sersic parameters. And if I’m being bold I would say that you might want to do away with sersic all together and fully embrace the gaussian mixture model. A 4 component gaussian mixture model (with shared center, ellip, PA) doesn’t have many more parameters than a sersic and is actually a really great model for faint galaxies. But decisions like that are well above my pay grade.
True, they are all near identical except the n>4 one. But these were 6 random galaxies, so in my head I was thinking if 1/6 of the fits are off by several chi^2/pix then that’s a big difference. Also the n=1.13 one moves by 0.1 chi^2/pix which isn’t visually a big change, but internally represents a lot (100s) of pixels moving by ~1 sigma. And yes I’m using the PSF with the image, not a 2 gaussian approximation (which I assume is for compute speed reasons).
The project I’m working on here is transient host association, and I’m using the object table sersic parameters to do the association. Now that I know about these factors, I’ve updated my code and I’m quite happy with how it works now. The issues I am raising at this point are for future people also getting surprised by broken expectations (such as the equality of Reff and Rhalf).
Anyway, there are no issues now on my end, I’ll mark your response as the solution since it seems this has moved up the chain to Jira. Thanks @dtaranu for your fast response and hard work building these data products for us!
Yes, that’s what the second ticket is for. I’ll add a note in the schemas and the Known Issues page for DP2 and DP1. We just don’t have an approved process for patching or revising entire columns of already-released datasets even if it’s technically possible (and in this case, a trivial multiplication).
I did have some nomenclature in mind for the Gaussian mixture approximations since they are different models in a sense (even if we’re not using the 8-component version yet), but I suppose we can just expand “Sersic” to “Gaussian mixture Sersic approximation” everywhere in the schemas.
Completely free Gaussian mixtures have too many free parameters and don’t converge well. You can do a linear fit of the Gaussian weights, and I tried this along with a variety of more flexible models (including just adding a central point source) and none of them actually improved the photometry statistically when testing on injected sources. There is certainly room for improvement for bright galaxies but unfortunately the (inevitable!) improvements in chi^2 from adding more free parameters doesn’t necessarily translate into smaller bias or scatter in practice.
And yes, the Gaussian mixture models (including the older CModel code) use Gaussian mixture PSF models because the performance benefit comes from having analytic convolutions.
