GLMoment() convergence issue for small objects

[reikonakajima]

The sigma size must be larger than e (2.718) for GLMoment to converge.
This must be a bug in the algorithm. GLMoment works in the mu=ln(sigma) space, so mu<1 must be causing some convergence issue. This is the issue to deal with this.

[reikonakajima]

A simple comparison between FDNT::GLMoment() and GalSim::AdaptiveMoment() show that, for the low S/N and small objects (note: no PSF involved!), the failure rate for GLMoment is way too high (>30%) compared to GalSim::AdaptiveMoment() (essentially 0%).

While the number 2.718 is a red herring, there definitely is an algorithmic issue. The GalSim::AdaptiveMoment also runs much faster (1ms per galaxy) than GLMoments (10ms per galaxy).

[reikonakajima]

Fixed a few things on branch #22:

1. dilating option was not implemented properly in GLSimple3.cpp, so that was fixed (does not change previous performance, where dilating was always set to true);
2. enforce a minimum stamp size of 32 pixels on the side (even though the predicted working space was very small for the small objects being tested---6 pixels!). This fixed unstable results due to the very limited number of pixels of small objects;
3. Having the correct initial size estimate was very important in making the dogleg iteration converge (fixed in megalut_tests/stampgrid.py).

However, after these fixes, I'm finding that

1. GalSim::AdaptiveMom() has a higher success rate in determining the ellipse than FDNT::GLMoments(). For a S/N = 18 object with HLR=0.9 pixel (sigma=0.9/1.17741=0.764 pixel) Gaussian galaxy (no PSF), GalSim::AdaptiveMom() fails 30+/-10% of the time, where FDNT::GLMoments() fail 60+/-10% of the time (+/-10% for a sample of 100).
2. In all cases, convergence fails by a run-away mu which becomes too small, causing the trustRadius of the dogleg to become too small (resulting in a DidNotConverge error).
3. The size sigma found in GLMoments() also tends to be biased high (where as AdaptiveMom() doesn't seem to have this problem).
   • Perhaps changing the step size of mu will improve the results....

GalSim::AdaptiveMom() iterates over the weighted moment to "center" the centroid, size and shear, while FDNT::GLMoments() applies an iterating least squares fitting algorithm. I believe this is where the advantage of GLMoments() is: it will allow multiple exposure fitting of the GL polynomials. [I do not know if a multi-exposure fitting of weighted moments can be calculated. Perhaps it is possible.]

[reikonakajima]

Perhaps changing the step size of mu will improve the results....

In GLSimple.cpp, I changed

const double MINIMUM_TRUST_RADIUS=0.03;
const double INITIAL_TRUST_RADIUS=0.4;

to

const double MINIMUM_TRUST_RADIUS=0.003;
const double INITIAL_TRUST_RADIUS=0.3;

which should decrease the step size. However, the GLMoment() failure rate for the similar set of sigma=0.76 Gaussians did not change (~60%).

The failure flag is always DidNotConverge (==128), due to trustRadius being too small.

[reikonakajima]

• Having the correct initial condition (i.e., calculating sigma from HLR, instead of feeding the HLR as the ic) doesn't seem to affect the bias on the size or the failure rate.
• Tinkering with dE (db349d3, 44f2d78) seems to help a bit, but within margin of error (55% failure rate instead of 60%).

[reikonakajima]

At this point, it's probably better just to start comparing precision and failure rates between GLMoments() and AdaptiveMom(). We might want to use AdaptiveMom() to measure our shapes.

[reikonakajima]

I am going to implement the GalSim::AdaptiveMom() algorithm of shape measurements as an option to the FDNT::GLMoment() fits. Then we will make the comparison in convergence and performance.

[reikonakajima]

Minor test on IC on convergence on small objects (tru_rad==0.9), on a set of images which gives AdaptiveMom() 27% failure. sigma = tru_rad / 1.3774 => 92% failure sigma = tru_rad / 1.2774 => 71% failure sigma = tru_rad / 1.1774 => 59% failure [this is the default IC, which should be the exact result] sigma = tru_rad / 1.0774 => 52% failure sigma = tru_rad / 0.9774 => 50% failure sigma = tru_rad / 0.8774 => 52% failure sigma = tru_rad / 0.7774 => 48% failure [this seems optimal... for intrinsically small radii] sigma = tru_rad / 0.6774 => 50% failure sigma = tru_rad / 0.5774 => 48% failure sigma = tru_rad / 0.4774 => 50% failure sigma = tru_rad / 0.3774 => 52% failure sigma = tru_rad / 0.2774 => 71% failure sigma = tru_rad / 0.1774 => 89% failure sigma = tru_rad / 0.0774 => 95% failure

[reikonakajima]

At sigma = tru_rad / 0.7774, GLMoments() has less than twice the failure rate than AdaptiveMom() for the small radii (tru_rad == 0.9). (e.g., 56% vs 29% failure rate, respectively)

[reikonakajima]

At large radii (tru_rad == 2.0), for a set of images with AdaptiveMom() having 4% failure: sigma = tru_rad / 1.4774 => 13% failure sigma = tru_rad / 1.3774 => 12% failure sigma = tru_rad / 1.2774 => 6% failure sigma = tru_rad / 1.1774 => 4% failure [default IC, which should be the exact result] sigma = tru_rad / 1.0774 => 4% failure sigma = tru_rad / 0.9774 => 2% failure sigma = tru_rad / 0.8774 => 1% failure [seems optimal... for normal size objects] sigma = tru_rad / 0.7774 => 3% failure sigma = tru_rad / 0.6774 => 2% failure sigma = tru_rad / 0.5774 => 6% failure sigma = tru_rad / 0.4774 => 15% failure sigma = tru_rad / 0.3774 => 34% failure

[reikonakajima]

It seems that GLMoments() prefers "slightly larger than truth" ICs.... at least for round objects.

[reikonakajima]

Now let's try tru_rad==2.0 and tru_g1,tru_g2 == 0.1,0.4. For a set of images with AdaptiveMom() giving 13% failure: sigma = tru_rad / 1.4774 => 21% failure sigma = tru_rad / 1.3774 => 13% failure sigma = tru_rad / 1.2774 => 8% failure sigma = tru_rad / 1.1774 => 6% failure [default IC, the exact result] sigma = tru_rad / 1.0774 => 5% failure sigma = tru_rad / 0.9774 => 4% failure sigma = tru_rad / 0.8774 => 5% failure sigma = tru_rad / 0.7774 => 3% failure [again, optimal around here] sigma = tru_rad / 0.6774 => 5% failure sigma = tru_rad / 0.5774 => 12% failure sigma = tru_rad / 0.4774 => 20% failure

[reikonakajima]

One last experiment: tru_rad==0.9 and tru_g1,tru_g2 == 0.1,0.4. For a set of images with AdaptiveMom() giving 40% failure: sigma = tru_rad / 1.4774 => 90% failure sigma = tru_rad / 1.3774 => 82% failure sigma = tru_rad / 1.2774 => 72% failure sigma = tru_rad / 1.1774 => 66% failure [default IC, the exact result] sigma = tru_rad / 1.0774 => 64% failure sigma = tru_rad / 0.9774 => 63% failure sigma = tru_rad / 0.8774 => 63% failure sigma = tru_rad / 0.7774 => 63% failure [seems optimal around here] sigma = tru_rad / 0.6774 => 63% failure sigma = tru_rad / 0.5774 => 63% failure sigma = tru_rad / 0.4774 => 62% failure sigma = tru_rad / 0.3774 => 63% failure sigma = tru_rad / 0.2774 => 73% failure sigma = tru_rad / 0.1774 => 89% failure sigma = tru_rad / 0.0774 => 86% failure (with long thinking time)

[reikonakajima]

Seems like ~50% larger-than-truth IC is optimal for fitting with GL moments.

[reikonakajima]

Probably already noted before, but: In general, GLMoment() failure happens when the trust radius shrinks too much (DidNotConverge or FitFlag==128), which seems to coincide with the object size sigma shrinking too much.

[reikonakajima]

Change mu updating scheme, when mu < 0, to some fraction of the original. 63% (original, 1.0) -> 58% (with 0.8) -> 57% (with 0.7 and 0.5) Compare with AdaptiveMom() failure rate of 40%.

Instead of fraction, try adding some number instead 63% (original, +0.0) -> 63% (+0.05) -> 68% (+0.15) -> 99% (+0.25)

This doesn't seem to help much.

[reikonakajima]

Since the bug (stamp size estimate on RunFDNT.cpp()) have been fixed, I will close this issue, and get back to testing. If any algorithmic improvement is needed, I will open a new issue.