virgil.fitting

A template model whose leaves at the paths in priors are fitted (their values are the starting point), or a function called with the parameters as keyword arguments, as for numpyro_model. The function may return a list of models, one per dataset, sharing parameters: for example a scene and a Rotated copy of it, for two epochs between which it turns. Regularizers then act on the first model only: an Image that appears only in a later model is not regularized.

A prior for each free parameter, keyed by its path (e.g. "comp.flux" or "env.log_brightness"; see image_priors). Priors on fluxes must have non-negative support. An angle (degrees) with an AngleVector prior is fitted as a 2-D vector, with no wrap boundary; its von Mises prior, if any, has a least-squares form.

A prior that is uniform in some coordinate of its parameter has a flat coordinate, and fit optimizes the parameter in it: LogUniform(a, b) is uniform in log x on [log a, log b], an isotropic inclination in cos i and an isotropic latitude in sin(lat) (any prior with a flat_coordinate() method), and Uniform in the parameter itself. The prior is then constant and adds nothing to the loss, so the fit is the maximum of the likelihood (times the other priors) inside the prior's range, and Levenberg–Marquardt applies. This is the maximum a posteriori in the coordinates in which the invariant (Jeffreys) prior is uniform, the choice consistent with virgil's prior rule; a mode in x itself would depend on the parametrization (for LogUniform, the density 1/x would pull every scale towards a). Other priors, such as Normal, Beta or HalfNormal, are evaluated in the model's own parameters, without the Jacobian of the bijection, as before. Either way the returned values are in the model's own parameters.

The data, fitted jointly. May be empty (()) when likelihoods holds all the data.

Penalties added to the loss, e.g. from virgil.imaging.

Priors on error-inflation terms to fit with the parameters: vis_scale and phi_scale multiply the uncertainties, and vis_error_rel (a fraction of the model visibility) and phi_error (radians) are added in quadrature (see inflated_errors), and the widths of gains correlated across channels, vis_gain_<group> (see OIData.with_gains), and of closure-phase offsets, phi_offset_<group> (see OIData.with_closure_offsets), and the wavelength scale, wavel_scale and wavel_offset (see OIData.with_wavelength_scale), and the North angle, north_angle in degrees (see OIData.with_north_angle), whose priors may be of either sign (e.g. Normal(1, 2e-4)). A dict applies to every dataset (values "noise.<term>"); a list gives each dataset its own ("noise[i].<term>"). An entry may be a function of the fitted values instead of a prior, tying the term to parameters in priors (see hierarchical_scales). The loss is then the full Gaussian negative log likelihood, including Σ log σ, so the default method is L-BFGS. So it is for data with gains, whose covariance depends on the model. Fitting error terms with an image is degenerate (a smoother image with larger errors fits as well): estimate them with a parametric model first. Priors. These terms are scale parameters, so their default (Jeffreys) prior is log-uniform on stated bounds; a Uniform(0, ...) favours large values. The bounds must contain the plausible values: for the factors vis_scale and phi_scale, whose neutral value is 1, e.g. dist.LogUniform(0.1, 10.0); for the added errors and widths (vis_error_rel, phi_error, vis_gain_<group>, phi_offset_<group>), e.g. dist.LogUniform(1e-4, 0.3). wavel_scale is a scale too: log-uniform about 1 unless a calibration gives a Gaussian (Normal(1, 2e-4) for GRAVITY is such information).

Starting values by path (or noise site), overriding the template's (required for a function model).

"lm": Levenberg–Marquardt (optimistix) on the residuals. Its inner solve is a dense QR of the Jacobian for up to 200 unconstrained coordinates, and otherwise matrix-free (cg_steps conjugate-gradient steps on the normal equations), so that the Jacobian of a larger image is never formed. The default when the whole objective has a least-squares form: Normal priors, and priors that are flat in their fitted coordinate (Uniform, ImproperUniform, LogUniform and the other priors with a flat coordinate). "lbfgs": L-BFGS (optax) on the loss, for penalties such as maximum entropy and total variation; the default otherwise. No unconstrained coordinate moves by more than max_step_size per step, so that log-brightness pixels cannot be switched off by one long step. "adam": Adam with learning_rate, run for max_steps.

Step limit (defaults: 1000 for LM, 20000 for L-BFGS, 2000 for Adam). It sets the length of LM's and Adam's loops, so a new value recompiles them (not L-BFGS); other numbers do not.

LM and L-BFGS stop when no component of the gradient of the loss per data point exceeds gtol, nor 1/1000 of its largest starting value (so that a fit started near a solution, e.g. along an L-curve, still converges). The tolerance is never below √eps of the dtype times that starting value (3.5e-4 times it in float32), which rounding errors in the gradient would not let the fit reach, nor below 1e-6 * gtol, so that a fit started exactly at a zero-residual optimum, whose gradient is rounding noise, is converged at once.

L-BFGS moves no unconstrained coordinate by more than this per step (for a log-brightness pixel, a factor exp(max_step_size)). Uncapped, its line search accepts steps that drive pixels so dark that their gradient vanishes and they never recover, leaving the fit at a spurious stationary point far above the minimum. Coordinates with real support (e.g. a position under a Normal prior) are in their own units, so raise this if they must move far.

Number of past steps L-BFGS keeps to model the curvature (default 50; optax's own default is 10). On regularized images, 10 left the fits short of their optimum at many weights, and weakly regularized ones running to the step limit; 50 found lower losses and converged in fewer steps, at a higher cost per step.

Adam's learning rate, in unconstrained coordinates.

Conjugate-gradient steps per LM step, for more than 200 coordinates. The inner solve runs for exactly this many steps, or as many as there are coordinates if fewer: its tolerances are zero, because an inner solve that stops at a step limit would abort the outer one.

Precision of the fit. The default runs in float64 inside a local jax.enable_x64 context. The returned model and values are cast back to JAX's precision outside the fit (float32, unless x64 is enabled), so that they work with the rest of your code.

Further Gaussian likelihood terms that are not visibilities: each is a callable of the fitted values (a dict, by path or keyword) that returns whitened residuals, such as PositionData.term or RVData.term for an orbit's positions and radial velocities. Their χ² follow the datasets' in info.

Wall-clock budget (seconds) for the optimizer, for LM and L-BFGS. The optimizer then runs in chunks of steps, and stops after the first chunk that ends past the limit, unconverged, with info["stop"] == "time" and a warning; the time to compile the first chunk counts. L-BFGS carries its whole state from chunk to chunk, so its path is the same as without a limit. optimistix's LM restarts from the last point in each chunk (with its damping reset), so its path can differ slightly. Adam, whose cost is set by max_steps, takes neither this nor progress.

Report progress between chunks (as for time_limit): a function called with a dict of steps, loss (the negative log posterior), grad_norm (as in info) and elapsed (seconds), or True to print them.

The fitted model, parameter values and diagnostics. A warning is raised if LM or L-BFGS did not converge; for L-BFGS it says whether the fit reached max_steps or time_limit, met a non-finite loss or gradient, or stopped earlier because a step no longer changed the parameters (its precision ran out, as can happen in float32). info["stop"] holds the reason.

info["grad_norm"] is the infinity norm (largest absolute component) of the gradient of the loss per data point at the returned point, in the unconstrained coordinates the optimizer works in: the quantity LM and L-BFGS compare with their tolerance (see gtol). It is computed the same way for every method, so a fit that stopped early (at max_steps, or after Adam's fixed number of steps) shows a larger value. L-BFGS tests the gradient at the point before its last step, so a converged L-BFGS fit's grad_norm can be slightly above the tolerance.

As for fit. The priors must have a least-squares form, as for fit's Levenberg–Marquardt: Normal priors, AngleVector priors, and flat ones (Uniform, ImproperUniform, LogUniform and the other priors with a flat coordinate). Flat priors add no curvature, as for Uniform. An angle vector's block is keyed by its site, "<path>_vec". The matrix is in numpyro's unconstrained coordinates (biject_to of each prior's support), which NUTS samples, not in fit's flat coordinates: for LogUniform(a, b) that is the logit of (x - a) / (b - a), not of log x.

As for fit. The priors must have a least-squares form, as for fit's Levenberg–Marquardt: Normal priors, AngleVector priors, and flat ones (Uniform, ImproperUniform, LogUniform and the other priors with a flat coordinate). Flat priors add no curvature, as for Uniform. An angle vector's block is keyed by its site, "<path>_vec". The matrix is in numpyro's unconstrained coordinates (biject_to of each prior's support), which NUTS samples, not in fit's flat coordinates: for LogUniform(a, b) that is the logit of (x - a) / (b - a), not of log x.

As for fit. The priors must have a least-squares form, as for fit's Levenberg–Marquardt: Normal priors, AngleVector priors, and flat ones (Uniform, ImproperUniform, LogUniform and the other priors with a flat coordinate). Flat priors add no curvature, as for Uniform. An angle vector's block is keyed by its site, "<path>_vec". The matrix is in numpyro's unconstrained coordinates (biject_to of each prior's support), which NUTS samples, not in fit's flat coordinates: for LogUniform(a, b) that is the logit of (x - a) / (b - a), not of log x.

The parameter values at which to take the curvature, normally fit(model, priors, data).values.

Further likelihood terms, as for fit (e.g. PositionData.term); data may then be (). Their residuals join the data's in J. A term's log_norm (a fitted RV jitter's, or the Jacobian of position_angle_prior) has no residuals, so its curvature is left out: the matrix is a preconditioner, so that costs efficiency, not correctness.

Keyword arguments for numpyro.infer.NUTS: inverse_mass_matrix, dense_mass, and adapt_mass_matrix=False. Warmup adaptation is switched off because a dense covariance estimated from a few hundred draws in thousands of dimensions is far worse than this matrix.