Regularisation

On noisy or under-constrained data the retrieved spectrum can overfit, picking up ripples in amplitude or oscillations in phase that fit the noise rather than the pulse. Every solver except COPRA (the L-BFGS family, the Levenberg–Marquardt family and the global CMAES) accepts mild penalties that suppress this, at the cost of a bias towards the chosen spectrum, smoothness or time window. The objective is defined in retrieval theory; the summary table lists which option each solver takes.

Two penalties are on by default: amplitude smoothness (reg_amp=0.03) and, when an independent spectrum is supplied, the spectral-match penalty (reg_spectrum=0.01). The LM family uses 3e-5 and 1e-5 instead, because it minimises squared residuals. Pass 0 to disable either. COPRA and copra-jax take none of these penalties; warm-lbfgs applies them in its L-BFGS stage.

Second-difference smoothness

Two penalties damp curvature via a discrete second difference (\(A_{k+1} - 2A_k + A_{k-1}\), the discrete analogue of \(\mathrm{d}^2/\mathrm{d}\omega^2\)):

• reg_amp — penalises curvature of the spectral amplitude;

• reg_phase — penalises curvature of the spectral phase.

res = croak.retrieve(trace, g.omega, delays, "pg", algorithm="lbfgs",
                    reg_amp=0.01, reg_phase=0.01, maxiters=300)

A second difference vanishes for a constant or linear amplitude and for a linear spectral phase (a time shift). Quadratic phase (pure GDD) has a constant, nonzero second difference and so pays a small, uniform penalty, while ripple and bin-to-bin oscillation have large second differences and are suppressed strongly. reg_phase therefore prefers smooth phase rather than zero chirp.

The penalties have closed-form gradients, accumulated into the solver gradient (L-BFGS) or appended as extra residual rows (LM).

Choosing the weights

croak’s defaults are reg_spectrum=0.01 and reg_amp=0.03. The validation retrievals in the papers use reg_spectrum=0.03 and reg_amp=0.03; use those to reproduce them. Three points matter when retuning:

• The natural diagnostics fail. The retrieved duration is nearly flat in the weight, and the agreement with the target spectrum improves monotonically with it by construction: a regulariser cannot be tuned on the quantity it pulls towards. On simulated data, tune against the complex-field error relative to the known truth; on measured data, prefer the defaults and check the spread over random starts.

• The two penalties work together. Too weak a spectral weight lets the amplitude collapse on chirped traces (a spiky spectrum a fraction of the true bandwidth); the amplitude-smoothness penalty suppresses exactly that collapse, which is why the pair works better than either alone.

• The measured spectrum has its own error. A spectral weight strong enough to impose it exactly imports that error, or a difference in how it was collected, into the retrieved pulse. Excessive smoothing likewise removes real structure.

L-BFGS vs LM weight scaling

LBFGS adds the penalty to the trace error \(R\); LM minimises the sum of squares, so its weight scales the penalty relative to \(R^2\). The penalty functions are identical, but a given numeric weight is not directly comparable between the two: at stationarity the correspondence is roughly \(\lambda_{\rm LM} \approx 2R\,\lambda\), which for the default reg_spectrum=0.01 at typical converged \(R\) puts the LM-side weight near 1e-5–3e-5. Retune when switching solver families.

reg_phase is most valuable in phase-only retrieval, where the amplitude is fixed and only the phase can absorb noise.

Spectral-match penalty

When you have an independently measured fundamental spectrum, you can pull the retrieved amplitude toward it (full amplitude-and-phase mode only) with reg_spectrum and spectrum_target — the measured spectral intensity on the retrieval grid, e.g. td.Iomega from preprocessing:

res = croak.retrieve(td.trace, td.omega, td.delays, td.interaction,
                    algorithm="lbfgs", reg_spectrum=0.05,
                    spectrum_target=td.Iomega, omega0=td.omega0_pulse)

This keeps the retrieval consistent with a trusted spectrometer measurement while still letting the trace determine the phase. retrieve_from_tracedata() wires td.Iomega in for you (and warns if reg_spectrum > 0 but no spectrum is available).

Temporal regularisation

Some artefacts live in the time domain rather than the spectrum: a noisy or under-constrained retrieval can grow spurious satellite pulses or a low pedestal far from the main feature. The temporal penalty reg_time suppresses the fraction of pulse energy \(|E(t)|^2\) that falls outside a time window time_window (a (t_lo, t_hi) tuple in seconds; default the measured delay range), with Planck-tapered edges so it does not introduce a hard cut:

res = croak.retrieve(trace, g.omega, delays, "pg", algorithm="lbfgs-ad",
                    reg_time=0.05, time_window=(-40e-15, 40e-15), maxiters=300)

Unlike the cosmetic post-retrieval temporal filter, this acts during retrieval and keeps the spectrum and field an exact Fourier pair, so it introduces no spectral artefacts; a window narrower than the pulse will, of course, suppress real satellites. It is available on LBFGSAD and CMAES (it needs the JAX field transform), and works in every parameterisation including the B-spline basis.

The B-spline phase basis

The penalties above are soft constraints added to the objective. The B-spline phase basis is a harder, structural one: instead of optimising the spectral phase per frequency bin, phase_basis="bspline" represents it as a cubic B-spline with a handful of control points (n_nodes, default 20) over the non-zero-spectrum support, holding the amplitude at the guess.

res = croak.retrieve(trace, g.omega, delays, "pg", algorithm="lbfgs-ad",
                    phase_basis="bspline", n_nodes=24, maxiters=300)

Because the phase then has only n_nodes degrees of freedom and is \(C^2\)-smooth by construction, it cannot ripple from bin to bin: a strong regulariser that also makes the problem low-dimensional enough for the derivative-free CMA-ES global search. It implies phase-only retrieval. A P-spline roughness term on the control points is folded into reg_phase, so you can still smooth the spline itself. Supported by lbfgs-ad, lbfgs-optx, lm, lm-optx and cma-es; see croak._jax_pulse.make_spline_phase_basis().

Per-frequency scaling

Not a penalty, but related: R_omega=True replaces the single intensity scale factor with a per-frequency vector, absorbing an uncalibrated frequency-dependent detection efficiency. It relaxes a constraint, so it pairs well with a spectral-match penalty that re-adds physical structure.

Rω lowers the reported error while degrading the spectrum

Do not use R_omega=True on its own. The per-frequency scale factors are exactly the information that pins the retrieved spectral amplitude: one free factor per row is the frequency marginal. With those extra parameters the amplitudes are constrained only through the delay structure of each row, which is weakest where the trace is weakest, at the edges of the band.

Because \(|E(\omega)|^2 \ge 0\) and the true value sits on that bound, the freed bins can only drift upward. The signature is the retrieved spectrum rising at both ends of the grid, and because \(I_\lambda = I_\omega\,2\pi c/\lambda^2\) the short-wavelength end looks far worse than its long-wavelength twin, by a factor \((\lambda_\mathrm{max}/\lambda_\mathrm{min})^2\), easily 5–10×. edge_energy_fraction() puts a number on it, and the plot_retrieval() spectrum panel reports it, in red above EDGE_ENERGY_WARN (1 %).

On a perfect, noiseless, exactly modelled PG trace whose true pulse has \(10^{-6}\) % of its energy in the outer 10 bins of each end:

	edge energy	trace error R
R_omega=False	0.02 %	7.9e-5
R_omega=True	5.4 %	4.0e-3

Rω usually reports the lower R, since it has one extra free parameter per frequency row, so the error alone cannot show that the spectrum has gone wrong. On our measured traces it was routine for Rω to halve R while tripling the edge energy.

The fix is to put back what Rω removed, from an independent measurement: phase_only=True (hold the amplitude at the measured spectrum) or a light reg_spectrum. Widening lam_min/lam_max helps too, because it moves the edge taper out onto empty trace. On one of our deep-UV TG-FROG datasets, widening the band from 180–500 nm to 140–750 nm, with a reference spectrum carrying 0.01 % edge energy:

	edge energy
Rω, tight band	15.6 %
Rω off, tight band	6.0 %
Rω off + reg_spectrum=1e-2, wide band	0.01 %

COPRA supports phase-only retrieval but not the spectral penalty.

Summary

Option	Solvers	Effect
reg_amp	LBFGS, LBFGSAD, LBFGSHand, OptxLBFGS, LM, OptxLM, CMAES	smooth spectral amplitude
reg_phase	LBFGS, LBFGSAD, LBFGSHand, OptxLBFGS, LM, OptxLM, CMAES	smooth spectral phase (great for phase-only)
reg_spectrum + spectrum_target	LBFGS, LBFGSAD, LBFGSHand, OptxLBFGS, LM, OptxLM, CMAES (full mode)	pull amplitude to a measured spectrum
reg_time + time_window	LBFGSAD, CMAES	suppress energy outside a time window
phase_basis="bspline"	LBFGSAD, OptxLBFGS, LM, OptxLM, CMAES	low-dimensional, intrinsically smooth phase
R_omega	all solvers	per-frequency intensity scaling — never alone, see above

(COPRA and its JAX twin take none of the penalties — their random-order local sweep already smooths implicitly — but they do support R_omega; see per-frequency scaling.)