Forward Generative Inversion (\(y \to x\))¶

← LatentNN | Regression Calibration → SIMEX

ForwardInversionRegressor estimates regression in the forward direction: a decoder \(g_\theta(y) \to x\) is trained to reproduce the observed features from the labels under an explicit measurement-error model, and the labels of a new object are then recovered by solving a weighted least-squares problem at test time. It addresses one specific failure mode of discriminative regression — attenuation when the features are much noisier than the labels — and it has known costs and failure modes of its own.

This page documents the estimator's assumptions, the statistical literature it belongs to, what the unit tests in this repository do and do not establish, and the practical hazards of test-time optimization.

Scope and evidence status¶

Why inverse regression can differ from direct regression¶

Consider the univariate noisy-input problem \(w = x + \epsilon\), \(\epsilon \sim \mathcal N(0, \sigma_x^2)\), with a linear relation \(y = \alpha + \beta x\). Least squares of \(y\) on the observed \(w\) estimates

so the fitted slope is diluted by a factor \(\lambda < 1\) and predictions of extreme labels are compressed toward the mean. This attenuation persists with training-set size, label precision and label-distribution shape. In the multivariate case it becomes a matrix shrinkage of the regression coefficients.

Fitting the other direction — a model of \(x\) given \(y\), then solving it for \(y\) — is the inverse regression approach to statistical calibration (Krutchkoff 1967; Osborne 1991). Berkson (1950) pointed out that the two regressions estimate different things and need not agree. The classical inverse estimator is not unbiased; its justification is a mean-squared-error argument, and the bias/MSE balance depends on the design and the noise level (Krutchkoff 1967; Osborne 1991; Fuller 1987). The estimator here is the likelihood version of that idea:

• under the assumed observation model and an identifiable, correctly specified decoder, the inversion is the maximum-likelihood estimate of \(y\); it is consistent as the information in the observables grows, but it carries the usual finite-sample bias of nonlinear MLEs;
• if the decoder or the noise model is misspecified, that bias propagates into \(\hat y\) — unlike the direct regression, which optimizes the prediction objective directly.

So "no attenuation" should be read as the estimator does not inherit the specific dilution factor of inverse regression, not as the estimator is unbiased in general.

Relation to auto-decoders and energy-based models¶

Structurally the estimator is not a new architecture:

• auto-decoder (Park et al. 2019): a decoder maps a low-dimensional code (the labels) to the high-dimensional observation; the code for a new sample is obtained by optimizing the decoder objective at test time instead of by a separate encoder. Here the codes are observed during training (supervised), which makes training easier than unsupervised auto-decoding, but inference is still an optimization problem.
• energy-based prediction (LeCun et al. 2006): prediction is \(\hat{y} = \arg\min_y E_\theta(x, y)\) with \(E_\theta(x, y) = \lVert \Sigma_x^{-1/2}(x - g_\theta(y)) \rVert_2^2\).
• non-amortised counterpart of an auto-encoder: an auto-encoder learns an explicit encoder (amortised inference); here the encoder is implicit. If inference cost matters, an encoder can be trained on the training set or on simulated inversions to amortise the fast path, with the decoder optimization used only for refinement. ForwardInversionRegressor implements only the decoder + inversion path.

The formulation is also a MAP inverse problem (Tarantola 2005), and closely related to likelihood-based event reconstruction in high-energy physics (matrix element method, Gainer et al. 2013) and to simulation-based inference more broadly (Cranmer et al. 2020).

Training: an explicit observation model¶

The default objective is the full heteroscedastic Gaussian likelihood with known per-observable measurement error \(\sigma_{x,i,d}\) and an optional learnable per-observable model-mismatch variance \(s_d^2\):

The log-determinant is not a constant: it lets the fit learn \(s_d\) and down-weight observables that the surrogate cannot reproduce. Two alternatives are available:

Plain unweighted MSE is not the default; it is obtained only by passing a constant sigma_x with learn_noise=False, and is the least informative choice when per-pixel errors are known. The full likelihood is what allows different observables to be weighted by their stated reliability.

Because \(y\) is the network input, the training label distribution \(p_{\text{train}}(y)\) shapes the design of the surrogate and, through it, the regions where the decoder is constrained. It does not enter the per-object estimator as a prior factor, unlike a discriminative model whose learned inverse map absorbs the marginal of \(y\).

Inference: weighted least squares¶

For a new object the estimate is

which is the Gaussian MLE (and the fixed-scale Student-t MLE). The training likelihood changes the fitted decoder and noise floor, not the form of the inversion. Two local optimizers are implemented:

The local covariance is the inverse Fisher information at the solution,

returned by predict_with_covariance. It is a linearised, local uncertainty: it propagates the measurement noise under the assumption that the objective is well approximated by a quadratic around the solution. It does not include decoder model error (except through the empirical \(s_d^2\)), and it can be badly misleading when the problem is weakly identified or strongly nonlinear.

Test-time optimization is nonlinear and often non-convex¶

This is a central practical caveat: at inference time the estimator solves a nonlinear least-squares problem, and neither implemented optimizer has a global-convergence guarantee.

• For a nonlinear decoder, \(E(y) = \lVert \Sigma_x^{-1/2}(x_{\text{obs}} - g_\theta(y)) \rVert_2^2\) can be multimodal, can have flat valleys, and can be locally ill-conditioned. Different initializations can converge to different stationary points that fit the same observation comparably well.
• Gauss–Newton has local quadratic convergence only near a solution with small residual and a well-conditioned Jacobian; far from it, the undamped step can overshoot or diverge. Levenberg–Marquardt damping and backtracking improve robustness but are heuristics, not guarantees (Bates & Watts 1988; Nocedal & Wright 2006).
• L-BFGS is a local descent method: it converges to a stationary point that depends on initialization and line-search behavior, and it can stall on non-smooth or badly scaled objectives. It should be treated as a cross-check, not as a trust-region guarantee.
• Identifiability can fail structurally (label directions that the decoder cannot distinguish) or practically (label directions with very little information at the observed noise level). In either case the optimizer may still return a well-converged but wrong estimate; the Fisher condition number and profile likelihood are the appropriate diagnostics (Raue et al. 2009).
• The reported covariance is local. Strong nonlinearity, bounded/periodic label directions, or multiple minima make the linearised covariance unreliable. For interval estimates use profile likelihood or parametric bootstrap; for multimodal posteriors use a multi-start or sampling approach.

Recommended protocol¶

1. Initialise sensibly. The default is the training-label mean; for extreme or structured objects pass init= from a physical prior, a nearest-neighbour label, or an amortised encoder prediction. Run several initializations and retain the solution with the lowest weighted SSE.
2. Inspect the diagnostics. predict(..., return_info=True) reports per-object final_weighted_sse, gradient_norm, relative_gradient_norm, and fisher_condition (condition number of \(J^\top\Sigma_x^{-1}J\)), plus the max_iter cap. There are no universal thresholds, but a large relative_gradient_norm, a very large fisher_condition, or a large final SSE relative to the stated noise indicates a problem with that object.
3. Cross-check the optimizers. Run Gauss–Newton and L-BFGS (and two initializations) and compare the weighted SSE. Agreement on the estimate is weak but useful evidence; disagreement is a red flag.
4. Check the model, not only the optimizer. Hold out objects with known labels, inject realistic systematics, and compare the recovered bias and scatter against a discriminative baseline. Optimizer convergence cannot repair a wrong decoder.
5. Use float64 for convergence studies. In float32, gradient norms and covariance entries can be dominated by round-off long before the tolerance is meaningful; the internal tolerances are heuristics. A relative_gradient_norm of order \(10^{-4}\)–\(10^{-3}\) in float32 is often the round-off floor rather than a statement about the optimum; rerun in float64 when stationarity is part of the reported result.
6. Mask known-bad observables. Defective pixels are not down-weighted by the Gaussian inversion; pass mask= for them. The Student-t likelihood robustifies training only, not the test-time objective.
7. Report the optimization settings. For scientific use, record the optimizer, initialization strategy, iteration budget, convergence statistics, masking and noise model, so the estimates are reproducible and auditable.

What changes relative to direct regression¶

When a discriminative model is the better choice¶

There is no general advantage of the generative direction; the following are reasons to expect a discriminative model to win at scale:

• Unmodelled structure. A forward likelihood must assign probability to every observed channel. Unmodelled systematics (bad columns, calibration drifts, missing physics, correlated detector effects) enter the energy unless they are masked or modelled, whereas a discriminative network can learn to down-weight or ignore defective regions. On messy data this robustness can dominate the attenuation benefit.
• Misspecification. The generative estimator inherits the decoder's bias; the discriminative estimator at least optimizes the prediction objective. The classical generative-versus-discriminative trade-off (Ng & Jordan 2002) applies: generative models can be statistically efficient when the model is correct, discriminative models are often more robust when it is not.
• Amortisation and throughput. A surrogate fit plus iterative Jacobian passes per object is far more expensive than one amortised forward pass. For \(10^8\)–\(10^9\) objects, a discriminative encoder (possibly with robust or EIV-corrected losses, or with a likelihood-based refinement on a subset) is usually the only practical option.
• Amortising changes the bias. An explicit encoder trained on the training set is again a discriminative estimator for the fast path, and inherits the same attenuation issue; the likelihood inversion only recovers the dynamic range for the objects on which it is run.

A defensible strategy is to compare both on matched splits and budgets, report tail-slice metrics and recovered dynamic range, and use the generative path where its assumptions hold and its diagnostics are clean — not as a default.

Comparison with other EIV methods¶

Example¶

import torch
from torchregress.algorithms import ForwardInversionRegressor

torch.manual_seed(0)

# Synthetic problem with a known linear generative model.
n_labels, n_observables, n = 2, 8, 512
y = torch.rand(n, n_labels) * 4 - 2
A = torch.randn(n_observables, n_labels)
X = y @ A.T + 0.05 * torch.randn(n, n_observables)  # known pixel noise = 0.05

reg = ForwardInversionRegressor(
    n_labels=n_labels,
    n_observables=n_observables,
    hidden_dims=(64, 64),
    sigma_x=0.05,
    learn_noise=True,
    likelihood="gaussian",     # full heteroscedastic likelihood (default)
    epochs=150,
    batch_size=128,
    seed=0,
)
reg.fit(X, y)

# Point estimates, local covariance, and inversion diagnostics.
y_hat, cov = reg.predict_with_covariance(X[:8])
y_hat, info = reg.predict(X[:8], return_info=True)
print("predicted labels:", y_hat[:3])
print("marginal std:", torch.diagonal(cov, dim1=1, dim2=2).sqrt()[:3])
print("final weighted SSE:", info["final_weighted_sse"])
print("relative gradient norm:", info["relative_gradient_norm"])
print("Fisher condition:", info["fisher_condition"])

A frozen, already-trained surrogate can be inverted without retraining:

from torchregress.algorithms import ForwardInversionRegressor

frozen = ForwardInversionRegressor(
    decoder=reg.decoder_.eval(),
    sigma_x=0.05,
    learn_noise=False,
    epochs=1,
)
frozen.fit(X, y, train_decoder=False)
y_hat = frozen.predict(X, method="lbfgs")

For data with defective pixels, train with the robust likelihood and mask the defects at inversion time (reusing X and y from above):

X_with_defects = X.clone()
X_with_defects[::10] += 20.0   # e.g. cosmic-ray hits / bad columns
defect_mask = torch.ones_like(X_with_defects, dtype=torch.bool)
defect_mask[::10] = False

robust = ForwardInversionRegressor(
    n_labels=n_labels,
    n_observables=n_observables,
    hidden_dims=(64, 64),
    sigma_x=0.05,
    likelihood="student_t",
    degrees_of_freedom=4.0,
    epochs=150,
    seed=0,
)
robust.fit(X_with_defects, y)
y_hat = robust.predict(X_with_defects, mask=defect_mask)

Diagnostics¶

Caveats¶

References¶

See also¶

• LatentNN — joint latent-input EIV training in the discriminative direction.
• Regression Calibration — moment-based correction of noisy inputs.
• SIMEX — simulation-extrapolation correction.
• Split-measurement & corrected-score losses — attenuation-free EIV objectives.
• EIV losses — losses that consume noisy inputs directly.
• Noisy features / measurement error — method-selection matrix.