Skip to content

virgil.fitting

fit finds maximum a posteriori parameters, taking the same arguments as numpyro_model, plus optional regularizers. gauss_newton_mass turns a fit into a mass matrix for numpyro's NUTS, for sampling any posterior with one peak: images, and orbits and binaries with fitted error scales and gains. gauss_newton_diag gives the per-coordinate scales that fit(..., precondition="diag") divides by, which conditions a fit of an image or other badly scaled parameters. Angles with an AngleVector prior are fitted as 2-D vectors, with no wrap boundary at 0°/360°.

Maximum a posteriori fits of models, including images.

fit takes the same arguments as numpyro_model: a model (a template, or a function of the parameters), a dict of numpyro priors whose keys are the free parameters, and the data, plus optional regularizers (see virgil.imaging). It finds the maximum a posteriori parameters with Levenberg–Marquardt, L-BFGS or Adam, optimizing each parameter in unconstrained coordinates through the bijection to its prior's support, in float64 by default. A parameter whose prior is uniform in some coordinate (a log-uniform scale, an isotropic inclination) is fitted in that flat coordinate, where its prior adds nothing to the loss, so that Levenberg–Marquardt works with the Jeffreys priors. To sample the same posterior, pass the same arguments to numpyro_model; it samples in numpyro's coordinates unless told flat_coordinates=True.

FitResult dataclass

The result of fit.

Attributes:

Name Type Description
model SourceModel or list

The fitted model, or models (one per dataset) if the model function returned a list.
values dict

The fitted parameter values, keyed by path, in the model's own parameters (not the flat coordinates some are fitted in). An angle with an AngleVector prior is in degrees at its path and its vector at "<path>_vec", so that values can start numpyro (init_to_value).
info dict

method; converged (Adam: the ftol_abs window stop fired or grad_norm <= gtol); steps; loss_history (Adam: the unscaled loss before each step); loss (the unscaled negative log posterior); chi2 and ndata, per dataset then per likelihood term (a Toeplitz term's re-evaluated directly at the result; with fitted error terms, on the inflated errors); likelihood_forms ("toeplitz", "direct" or "direct (auto): <unmet conditions>"); chi2_red, the total χ² per data point; grad_norm, the largest absolute component of the per-point loss gradient in the unconstrained coordinates of the convergence test (gtol in fit); stop, why an unconverged fit ended (None if it converged): "limit" (max_steps), "time" (time_limit), "stalled" (an L-BFGS step no longer moved the parameters), "non-finite" (a NaN or infinite loss or gradient; any method ending on a non-finite loss is unconverged) or "failed" (another LM failure); at_bound, a tuple of the parameters and error terms that ended within 1/1000 of their prior's range of a finite end, in the coordinate they are fitted in (a fit at a bound, say an eccentricity at 0.9, is often not at a mode); and derived, quantities that likelihood terms compute and that are not parameters, such as a radio term's profiled error scales ŝ ("cal.sigma"). With fitted error terms values holds them too.

model instance-attribute

values instance-attribute

info instance-attribute

__init__(model, values, info)

GaussNewtonMass

Bases: dict

NUTS keyword arguments from gauss_newton_mass.

A dict that unpacks into NUTS(...) as usual, and records the flat_coordinates of its matrix as an attribute for chain_init_params. condition is the condition number of the matrix's correlation form (None for a diagonal mass).

flat_coordinates = False class-attribute instance-attribute

condition = None class-attribute instance-attribute

NoLeastSquaresForm

Bases: TypeError

The objective has no residual vector (use L-BFGS or Adam).

A TypeError for backward compatibility; fit picks L-BFGS on it and variational falls back to a diagonal start, while other errors propagate.

NoStartingValue

Bases: ValueError

A free parameter has no starting value (a function model needs one).

SingularCurvature

Bases: ValueError

The Gauss–Newton curvature is singular at the given values: some direction is constrained by neither the data nor a proper prior (an ImproperUniform or flat-in-coordinate parameter the data do not fix), or the matrix is not positive definite in the working precision.

fit(model, priors, data, regularizers=(), *, noise=None, init=None, method=None, max_steps=None, gtol=0.0001, max_step_size=2.0, lbfgs_memory=None, learning_rate=0.01, cg_steps=50, dtype='float64', likelihoods=(), time_limit=None, progress=None, precondition=None, adam_eps=1e-08, ftol_abs=None, ftol_window=100, post_step=None)

Find the maximum a posteriori parameters of a model given data.

Arguments are shared with numpyro_model (common arguments). Parameters under a prior uniform in some coordinate are optimized in it; values are in the model's own.

Parameters:

Name Type Description Default
model SourceModel or callable

A template whose leaves at the paths in priors are fitted, or a function of the parameters returning a model (or one per dataset).

 required
priors dict[str, Distribution]

A prior for each free parameter, keyed by its path.

 required
data OIData or sequence of OIData

The data, fitted jointly; () if likelihoods has it all.

 required
regularizers sequence and dict

Loss penalties (e.g. virgil.imaging) and priors on error terms; see common arguments.

 ()
noise sequence and dict

Loss penalties (e.g. virgil.imaging) and priors on error terms; see common arguments.

 ()
init dict and str

Starting values by path (required for a function model); "lm" (the default for a least-squares objective), "lbfgs" (for penalties such as maximum entropy, which LM also fits in float64) or "adam".

 None
method dict and str

Starting values by path (required for a function model); "lm" (the default for a least-squares objective), "lbfgs" (for penalties such as maximum entropy, which LM also fits in float64) or "adam".

 None
max_steps int and float

Step limit (default 1000 LM, 20000 L-BFGS, 2000 Adam), gradient tolerance and L-BFGS step cap; see fit stopping.

 None
gtol int and float

Step limit (default 1000 LM, 20000 L-BFGS, 2000 Adam), gradient tolerance and L-BFGS step cap; see fit stopping.

 None
max_step_size int and float

Step limit (default 1000 LM, 20000 L-BFGS, 2000 Adam), gradient tolerance and L-BFGS step cap; see fit stopping.

 None
lbfgs_memory (int, int and float)

L-BFGS history (default 100, cut to keep it under 1 GiB); CG steps per LM step (above 200 coordinates); Adam's rate: a float or an optax schedule (a callable of the step).

 None
cg_steps (int, int and float)

L-BFGS history (default 100, cut to keep it under 1 GiB); CG steps per LM step (above 200 coordinates); Adam's rate: a float or an optax schedule (a callable of the step).

 None
learning_rate (int, int and float)

L-BFGS history (default 100, cut to keep it under 1 GiB); CG steps per LM step (above 200 coordinates); Adam's rate: a float or an optax schedule (a callable of the step).

 None
dtype (float64, float32)

Precision; see common arguments.

 "float64"
likelihoods (sequence, float and callable)

Further terms, e.g. PositionData.term; wall-clock budget (s); callback (or True) on steps, loss, grad_norm, elapsed (grad_norm is NaN for Adam). A time_limit that runs out on Adam's or LM's last chunk reports "limit", not "time".

 ()
time_limit (sequence, float and callable)

Further terms, e.g. PositionData.term; wall-clock budget (s); callback (or True) on steps, loss, grad_norm, elapsed (grad_norm is NaN for Adam). A time_limit that runs out on Adam's or LM's last chunk reports "limit", not "time".

 ()
progress (sequence, float and callable)

Further terms, e.g. PositionData.term; wall-clock budget (s); callback (or True) on steps, loss, grad_norm, elapsed (grad_norm is NaN for Adam). A time_limit that runs out on Adam's or LM's last chunk reports "limit", not "time".

 ()
adam_eps (float, float and int)

Adam only: optax.adam's eps; and a plateau stop after 2 * ftol_window steps, at the first step where the last window's loss range is below ftol_abs and its mean is not lower than the previous window's by ftol_abs or more (unscaled; off if None).

 1e-08
ftol_abs (float, float and int)

Adam only: optax.adam's eps; and a plateau stop after 2 * ftol_window steps, at the first step where the last window's loss range is below ftol_abs and its mean is not lower than the previous window's by ftol_abs or more (unscaled; off if None).

 1e-08
ftol_window (float, float and int)

Adam only: optax.adam's eps; and a plateau stop after 2 * ftol_window steps, at the first step where the last window's loss range is below ftol_abs and its mean is not lower than the previous window's by ftol_abs or more (unscaled; off if None).

 1e-08
post_step callable

Adam only, not with precondition: a jax-traceable post_step(z) -> z on the unconstrained coordinates z (a dict by path, not values), applied after each update. Source: Charles et al. (2026).

 None
precondition (None, diag)

"diag" rescales each coordinate by the Gauss–Newton curvature diagonal for L-BFGS and Adam (provisional); see conventions. info["precondition_rounds"] records the spread per round. With a maximum-entropy regularizer at weak weights float32 can stall: use float64. Source: Skilling & Bryan (1984), Hutchinson (1989).

 None

Returns:

Type Description
FitResult

Model, values and diagnostics; warns if LM, L-BFGS or Adam did not converge (info["stop"] gives the reason, info["grad_norm"] the gradient).

gauss_newton_mass(model, priors, data, values, *, likelihoods=(), noise=None, flat_coordinates=False)

A NUTS mass matrix from the Gauss–Newton curvature at a fit.

Near the maximum a posteriori the posterior is close to a Gaussian with precision F = JᵀJ + diag(P) in numpyro's unconstrained coordinates; its inverse is the mass matrix. See the Gauss–Newton mass matrix for the construction and when not to use it. Use it at fixed field hyperparameters, with every sampled parameter in priors. Above 4096 coordinates the mass is diagonal, with a UserWarning.

Parameters:

Name Type Description Default
model

As for fit. Priors must be elementwise (a MultivariateNormal or Dirichlet raises NoLeastSquaresForm).

 required
priors

As for fit. Priors must be elementwise (a MultivariateNormal or Dirichlet raises NoLeastSquaresForm).

 required
data

As for fit. Priors must be elementwise (a MultivariateNormal or Dirichlet raises NoLeastSquaresForm).

 required
likelihoods

As for fit. Priors must be elementwise (a MultivariateNormal or Dirichlet raises NoLeastSquaresForm).

 required
noise

As for fit. Priors must be elementwise (a MultivariateNormal or Dirichlet raises NoLeastSquaresForm).

 required
values dict

Where to take the curvature, normally fit(...).values (model's own parameters).

 required
flat_coordinates bool

Pass the value given to numpyro_model; see common arguments.

 False

Returns:

Type Description
GaussNewtonMass

Keyword arguments for numpyro.infer.NUTS; condition is the correlation matrix's condition number (None if diagonal; above 1e10 the matrix is unreliable).

Raises:

Type Description
SingularCurvature

If a direction is fixed by neither the data nor a proper prior, or the matrix is not positive definite.

Examples:

>>> result = fit(scene, priors, data)
>>> kernel = NUTS(numpyro_model(result.model, priors, data),
...               init_strategy=init_to_value(values=result.values),
...               **gauss_newton_mass(scene, priors, data, result.values))

gauss_newton_diag(model, priors, data, values, *, regularizers=(), likelihoods=(), noise=None, n_probes=None, dtype='float64')

The diagonal of the Gauss–Newton curvature of a fit's loss.

Provisional

May change in any 0.4.x release without a deprecation cycle.

The curvature is F = JᵀJ + diag(P) + diag(R) in the unconstrained coordinates fit optimizes: J is the Jacobian of the whitened residuals, P the curvature of the priors and R the Hessian diagonal of the regularizers, floored at zero. On the pixels of a regularizer with metric_rows (maximum entropy) it is that metric alone, weight * b * (1 - b) (the softmax's, unsmoothed pixels). The data JᵀJ and other terms are dropped there, so it is a scale, not the curvature, at weak weights. It is the scale fit(..., precondition="diag") divides each coordinate by; see gauss_newton_mass for the dense matrix.

Source: probes, Bekas et al. (2007) and Hutchinson (1989); curvature, Skilling & Bryan (1984), Cornwell & Evans (1985), Martens (2010), Schraudolph (2002).

Parameters:

Name Type Description Default
model

As for fit.

 required
priors

As for fit.

 required
data

As for fit.

 required
regularizers

As for fit.

 required
noise

As for fit.

 required
likelihoods

As for fit.

 required
dtype

As for fit.

 required
values dict

Parameter values at which to take the curvature (e.g. fit(...).values).

 required
n_probes int

Up to 1024 coordinates the diagonal is exact, one Jacobian-vector product per coordinate. Above, it is estimated: JᵀJ from the mean of (Jᵀu)² over 16 Rademacher vectors u, R from 16 Rademacher Hessian-vector products, each good to a relative √(2 / probes). Give a number of probes to estimate at any size, or 0 for the exact diagonal at any size.

 None

Returns:

Type Description
dict

The diagonal by coordinate name, with each coordinate's shape.