-
Notifications
You must be signed in to change notification settings - Fork 26
PAHFIT Speed
Python PAHFIT is somewhat slow, a bit slower (on my machine) than IDL PAHFIT was. Once typical spectra have 5-10x as many wavlelength samples, and the model is expanded to have 3-5x as many parameters, this will really start to hurt. More than several minute fit times on decent hardware are quite likely. Once the model has settled a bit more, striving for performance gains will be important.
Example: M101_Nucleus_irs.ipac (from /data/)
With default settings (xtol=ftol=gtol=1D-10), the fit converges in:
- 71 “iterations”, using 6754 function evals
- 3.9 seconds
- 1732 FEVs/second
With the default acc=eps=1D-10, the fit converges in:
- 6.6s
- 814 function evaluations
- 123 FEVs/second
So Python PAHFIT (including overhead) is about 14x slower per model function call, but makes up some of this by converging in fewer calls.
Fitting the Akari + Spitzer SL/LL template spectrum is quite slow at present (v2.1)…
- 58.7s
- 4566 FEVs
- 77.8 FeVs/sec (oy!)
Much slower with the new science pack. A profiler report indicates much of this time in spent
Anecdotal evidence is that these spectra take minutes to hours to converge.
Allowing all lines and features to vary in position and width is very expensive. And hopefully (for unresolved lines), unnecessary.
A good instrument and science pack model should nail the position and widths of lines precisely. Any mismatch in the residuals should be considered an “opportunity to improve” either the rest line wavelength in the sci-pack, or the fwhm estimation at each wavelength in the instrument pack. There is some concern about NIRSPEC wavelength calibration repeatability, but this should be corrected prior to PAHFIT’ing.
astropy.modeling includes bounds of significance or bounding boxes. The classic case is a 2D star PSF on a large image; evaluating one stars PSF across the entire image is needlessly wasteful. Some savings may occur for lines and features.
5 * max(fwhm) should do it (or 3x for Gaussians, 5x for Drudes; can investigate).
Should profile slow fits to see where the time is being spent. Mostly in the model function(s)?
Numba is a nice and very easy-to-use decorator-based framework for transpiling to C, with numpy support. It is not compatible with scipy functions (which themselves call to compiled C/FORTRAN code). But the bottommost objective function apparently can be @jit’d: Example of optimize.leastsq with Numba.
To use numba we would probably need to re-write the model as a set of simple functions call from the objective function, and give up on CompoundModel.
The fit_deriv callable for fittable models can be used to speed up derivative estimation.
This would mean giving up on the conveniences of the CompoundModel, which in practice might mean just switching to scipy.least_squares (newer than the scipy.leastsq that LevMarLSQFitter uses internally). It natively supports bounds.
A very common source of optimization not converging well is human error in the computation of the gradient. You can use scipy.optimize.check_grad() to check that your gradient is correct. It returns the norm of the different between the gradient given, and a gradient computed numerically.
Some consider analytic derivatives not worth the overhead. Probably this would be a last resort.