Choosing a solver#

croak has ten solvers in four families: COPRA, L-BFGS, Levenberg–Marquardt and the global CMA-ES, several with on-device JAX or Optimistix twins. Choose first by the model and constraints a solver supports, then by its speed and final accuracy on your data. The local solvers minimise the same least-squares objective but differ in penalty scaling and stopping rules. This page is the short version; the algorithm tutorial runs an executed comparison.

croak’s default, warm-lbfgs with reg_spectrum=0.01 and reg_amp=0.03, is a general-purpose starting point. The validation retrievals in the papers instead use lbfgs-ad from a transform-limited start for 400 iterations with reg_spectrum=0.03 and reg_amp=0.03; use those settings to reproduce them.

At a glance#

Start here; these six cover almost every case:

Solver	Gradient	Best for	Notes
`WarmLBFGS`	JAX autodiff	the default; when you do not want to think about the initial guess	`lbfgs-ad` seeded from a COPRA local sweep: same model, far less sensitive to where it starts
`COPRA`	analytic adjoint	full amplitude + phase in the core model	fast initial descent, robust from poor guesses
`LBFGS`	analytic adjoint	phase-only; custom regularisation	steady descent; analytic gradients for the core model
`LBFGSAD`	JAX autodiff	anything the core model does not cover	the solver with the widest model coverage: smearing, fitted thickness/τ₀/smearing, B-spline basis, temporal penalty
`LM`	JAX Jacobian	polishing to the lowest error	superlinear; use after COPRA/L-BFGS; its penalty weights are scaled differently from L-BFGS
`CMAES`	none (population search)	escaping a poor local minimum with no good guess	population search; pair with `phase_basis="bspline"`; needs `croak[evo]`

Reach for lbfgs-ad (or lm/lm-optx) whenever you need something the analytic adjoint does not carry — see Gradients for exactly where the boundary lies.

The remaining four are on-device reimplementations of the same algorithms:

Solver	Twin of	Why reach for it
`COPRAJax`	`copra`	same COPRA, on-device JAX numerics
`LBFGSHand`	`lbfgs`	analytic gradient, jitted/vectorised in JAX (no AD); typically the fastest L-BFGS
`OptxLBFGS`	`lbfgs-ad`	the L-BFGS loop runs entirely on-device
`OptxLM`	`lm`	an on-device JAX-Jacobian LM

By task#

General amplitude-and-phase retrieval, clean to moderately noisy data → copra. It converges fast and tolerates a rough initial guess.
Phase-only retrieval (amplitude fixed to a measured spectrum) → lbfgs with phase_only=True, optionally reg_phase.
Noisy experimental data → lbfgs with mild smoothness regularisation (reg_amp/reg_phase).
Uncalibrated spectral response → R_omega=True fits a gain for each frequency row, but never on its own: add phase_only=True or reg_spectrum with a measured spectrum, or the retrieved spectrum grows spurious energy at the band edges while \(R\) falls. See per-frequency scaling.
Squeezing out the last factor in \(R\) → run copra or lbfgs, then a short lm polish (20–50 iterations).
Stuck in a poor local minimum / no usable initial guess → a cma-es global search, usually with phase_basis="bspline" (≈20–30 spline nodes), then refine the result with lbfgs-ad. See global retrieval.
Spurious satellite pulses or a temporal pedestal → lbfgs-ad (or cma-es) with the temporal penalty reg_time/time_window (see Regularisation).
Keeping the loop on a JAX device → the twins copra-jax, lbfgs-hand, lbfgs-optx and lm-optx. Their options and stopping rules follow their own optimisers; time the first (compiling) call separately from repeated runs.
Checking a hand-derived gradient or a new interaction → compare the derivative with finite differences and autodiff on complex test fields; see numerical validation.
lm-optx iterations feel slow on a large grid → it already defaults to the faster normal-equations solve (linear_solver="normal"); the remaining cost is the Jacobian and linear solve, which grow with the parameter count, so phase_basis="bspline" is the big lever. Fall back to linear_solver="qr" only if the solve misbehaves (see Retrieval algorithms).

Run from several random restarts#

The objective is non-convex, so a single run can land in a poor local minimum. The standard practice is to retrieve from 5–10 random initial guesses and keep the result with the lowest error:

import numpy as np, croak

rng = np.random.default_rng(0)
best = None
for _ in range(8):
    res = croak.retrieve(trace, g.omega, delays, "pg",
                        algorithm="copra", guess=None, rng=rng, maxiters=300)
    if best is None or res.error < best.error:
        best = res
print(best.error)

guess=None draws a random Gaussian-spectrum start. Reusing one seeded generator draws a different, reproducible start for each run; resetting it to the same seed every time would repeat the same start.

For a penalised fit, error is only the trace term, so also compare the penalised objective and the plausibility of the pulse. Start-to-start variation is a convergence diagnostic, not a measurement uncertainty.

A two-stage recipe#

For the lowest final error, converge cheaply and then polish:

# 1. converge with COPRA (or LBFGS)
res = croak.retrieve(trace, g.omega, delays, "pg", algorithm="copra", maxiters=300)

# 2. polish with Levenberg–Marquardt, seeded from the COPRA spectrum
res = croak.retrieve(trace, g.omega, delays, "pg", algorithm="lm",
                    guess=res.spectrum, maxiters=50)

Dispersive problems#

Solver choice is independent of whether you model a slab: pass material/thickness/npoints/omega0 to any solver to retrieve through a dispersive medium (PG/SD only — see Interactions). “D-COPRA” in the literature is exactly COPRA with a material set.

Geometric smearing is the exception: pass smearing to one of the autodiff solvers (warm-lbfgs, lbfgs-ad, lm, lm-optx, lbfgs-optx, cma-es) only. The incoherent sum over quadrature nodes breaks COPRA’s magnitude-replacement projection and is not covered by the analytic adjoints, so those solvers raise instead of silently ignoring the kernel. See Geometric smearing.

Convergence and stopping#

COPRA: maxiters, plus optional reltol (stop on a plateau of the best error) and abstol (stop at a target error). By default it runs all maxiters.
LBFGS / LBFGSAD / LM: maxiters caps the evaluations; reltol/abstol are the optimiser’s function/step tolerances (defaults 1e-4/1e-8). For LM these become SciPy’s ftol and xtol/gtol.

Tip

Polishing with LM: tighten abstol, not reltol. For the SciPy route, abstol sets both the step (xtol) and gradient (gtol) tolerances, and its 1e-8 default is a step-size floor: it is what stops an LM polish, usually while the error is still falling steeply. On the clean PG trace of the algorithm tutorial, tightening it to 1e-12 lowers the final error by about three orders of magnitude for about 1.4× the evaluations, while tightening reltol changes nothing.

By default croak.retrieve() shows a live progress bar and prints a short end-of-run summary (progress="auto"); pass progress=False to silence it or progress=True to force it (e.g. when redirecting output). For programmatic monitoring, supply your own callback(iteration, R, best_R) (which replaces the bar), or inspect res.errors after the fact with croak.plotting.plot_convergence(). See croak.progress.