Skip to content

replace Monte Carlo uncertainty estimation with Jacobian-based covariance #250

Description

@moustakas

@jdbuhler, below is a potential speed-up suggested by claude-code. It's probably a high-risk, high-reward change, but I welcome your thoughts before I proceed in this direction.


The emission-line fitter computes parameter uncertainties by running nmonte (default 50) full optimizer passes on noise-perturbed flux realizations (emlines.emline_specfit.get_results). Each pass calls EMFitTools.optimize, which runs scipy.optimize.least_squares with tr_solver='lsmr'. On a typical object these 50 MC passes account for ~65–70% of emline_specfit wall-clock time.

Because optimize already supplies an analytical Jacobian (jac=obj.jacobian), least_squares computes and stores fit_info.jac at the best-fit solution at no extra cost. The parameter covariance matrix follows directly:

  Cov(θ) ≈ (J^T W J)^{-1}

where W is the diagonal inverse-variance weight matrix. Parameter uncertainties are sqrt(diag(Cov)). This replaces all MC passes with a single matrix inversion.

Tradeoffs to evaluate before adopting:

  • Linear approximation: valid when the model is well-approximated as linear around the best fit, but may underestimate uncertainties near parameter bounds or for low-S/N detections.
  • Bound constraints: parameters at their bounds have non-Gaussian posteriors; the Jacobian covariance will underestimate those uncertainties.
  • Should be validated against the current MC results on a representative sample spanning a range of S/N and spectral types before replacing it in production.

Metadata

Metadata

Assignees

No one assigned

    Labels

    might-doMight be addressed.

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions