Hessian-Free Self-Normalized Inference

• The paper demonstrates that perturbation-based sensitivity recovers predictive variance without explicitly computing, storing, or inverting the Hessian.
• Empirical results show competitive performance in uncertainty estimation, especially for out-of-distribution and in-between tasks.
• Variants such as amortized pre-training and finite-difference normalization offer trade-offs between computational efficiency and statistical accuracy.

Hessian-Free Self-Normalized Inference is not presented as a canonical method name in the cited literature. In the available arXiv record, it is better understood as a composite description for inferential procedures that avoid explicit Hessian formation, storage, inversion, or eigendecomposition while producing an uncertainty quantity through a data-dependent normalization or sensitivity scaling. Three nearby research programs define the technical landscape: Hessian-Free Laplace, which recovers linearized Laplace predictive variance from sensitivity to a small prediction-dependent regularization perturbation (McInerney et al., 2024); Hessian-free second-order variational inference based on reparameterized Hessian-vector products (Fan et al., 2015); and approximate Newton-based statistical inference, which estimates sandwich or HAC covariance matrices from stochastic approximate Newton replicates built only from gradients (Li et al., 2018). The “Hessian-free” component is therefore direct and well developed, whereas the “self-normalized” component ranges from a weak finite-difference scaling by $1/\lambda$ to a covariance-estimating, studentization-like normalization.

1. Terminological scope and conceptual boundary

The exact phrase “self-normalized” does not appear in the Hessian-Free Laplace paper, and that work does not present an importance-weight normalization or ratio estimator. The closest normalization is the scaling by $1/\lambda$, which turns a finite perturbation into a derivative estimate (McInerney et al., 2024). In the second-order variational-inference paper, there is likewise no self-normalized importance sampling, normalized weights, or SNIS-style estimator; Monte Carlo quantities are ordinary sample averages under a reparameterized Gaussian base distribution (Fan et al., 2015). The approximate Newton paper comes closest to a self-normalization interpretation, because it constructs a data-driven covariance estimate from stochastic gradient replicates and then uses that estimated covariance for confidence intervals and p-values, although it does not adopt the term “self-normalized” either (Li et al., 2018).

Paper	Core inferential object	Relation to the topic
(McInerney et al., 2024)	Linearized Laplace predictive variance	Hessian-free, weak normalization by $1/\lambda$
(Fan et al., 2015)	Curvature of variational objectives	Hessian-free, not self-normalized
(Li et al., 2018)	Sandwich/HAC covariance	Hessian-free, studentization-like normalization

This suggests that “Hessian-Free Self-Normalized Inference” is best treated as an umbrella label rather than a standardized doctrine. The common structure is indirect access to curvature-sensitive uncertainty without explicit second-order matrix algebra; the substantive differences lie in what is being inferred, how normalization is implemented, and which theoretical regime is assumed.

2. Hessian-Free Laplace as prediction-sensitivity inference

The most direct bridge between Hessian-free computation and uncertainty quantification appears in “Variation Due to Regularization Tractably Recovers Bayesian Deep Learning” (McInerney et al., 2024). The paper starts from MAP estimation for a neural network $f_\theta$,

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$

and recalls the Laplace approximation around $\hat\theta$, with precision

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$

For a locally linearized network, the predictive epistemic variance is

$$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$

The core construction replaces explicit computation of $P^{-1}$ by a perturb-and-reoptimize procedure. The training objective is modified to

$$\mathcal L^{(x,\lambda)} = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i)) + \log p(\theta) + \lambda f_\theta(x),$$

with optimum

$1/\lambda$0

The uncertainty estimate is then obtained from the sensitivity of the prediction at $1/\lambda$1 to the perturbation,

$1/\lambda$2

for small $1/\lambda$3. The method is deterministic and optimization-based: there is no posterior sampling during inference, and no ensemble of independently trained models is required. One computes one baseline optimum and one perturbed optimum, then uses their prediction difference.

By differentiating the perturbed first-order optimality condition at $1/\lambda$4, the paper proves that

$1/\lambda$5

and therefore

$1/\lambda$6

This is the precise sense in which regularization variation recovers Laplace: under standard Laplace assumptions, the directional derivative of the prediction induced by prediction-regularized re-optimization equals the local linearized Laplace predictive variance.

The assumptions stated are that $1/\lambda$7 exists, $1/\lambda$8 exists, and $1/\lambda$9 is continuously differentiable with respect to $1/\lambda$0. The method is therefore not Hessian-free in a philosophical sense, because the Hessian must exist and be nonsingular for the derivation, but it is Hessian-free in the computational sense that it does not explicitly form, invert, or eigendecompose the Hessian. In this formulation, epistemic uncertainty is read off from how much the optimum can move in prediction-relevant directions under a tiny perturbation of the objective.

3. Finite-difference normalization, amortization, and variants

The practical Hessian-Free Laplace estimator replaces the derivative by a finite difference,

$1/\lambda$1

leading to the approximation

$1/\lambda$2

In algorithmic practice, the paper returns

$1/\lambda$3

using an absolute value to ensure a nonnegative variance estimate in finite-$1/\lambda$4 settings. That absolute-value operation is presented as a practical stabilization rather than part of the infinitesimal derivation (McInerney et al., 2024).

The main computational tradeoff is point specificity. In the plain form, a new fine-tuning solve is required for each query $1/\lambda$5. The paper addresses this by proposing pre-trained HFL, which amortizes the perturbation over evaluation points $1/\lambda$6. It defines

$1/\lambda$7

and trains with

$1/\lambda$8

The paper also gives an absolute-valued regularizer,

$1/\lambda$9

an in-sample pre-training variant,

$f_\theta$0

a data-augmentation form under Gaussian likelihood,

$f_\theta$1

and a parameter-uncertainty variant,

$f_\theta$2

The experiments are regression only. The evaluation uses the synthetic tasks Quadratic-Uniform, Quadratic-Inbetween, Sin-Uniform, and Sin-Inbetween, with a feedforward network with 50 hidden units and $f_\theta$3 activations. The paper compares Exact Hessian Laplace, GGN Laplace, an Eigenvector approximation of GGN, and HFL, using PICP, CRPS, and NLL. Its headline empirical finding is that HFL performs comparably to exact or approximate Hessian-based methods and is often competitive or best, especially on OOD or in-between uncertainty tasks. Because HFL re-optimizes the entire parameter vector $f_\theta$4, the method in principle captures uncertainty contributions from all layers rather than only the last layer.

4. Reparameterized Hessian-free variational inference

“Fast Second-Order Stochastic Backpropagation for Variational Inference” develops a different Hessian-free program: second-order optimization for Gaussian variational inference via reparameterized derivatives and Hessian-vector products (Fan et al., 2015). The basic reparameterization is

$f_\theta$5

which rewrites

$f_\theta$6

The paper first recalls exact Gaussian second-order identities, then avoids the associated third- and fourth-order derivatives by differentiating through $f_\theta$7. Writing $f_\theta$8 and $f_\theta$9, it derives

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$0

and

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$1

The paper therefore derives the true Hessian of the reparameterized expectation, but practical optimization uses either Hessian-vector products or quasi-Newton approximations rather than explicit dense Hessians.

The Hessian-free mechanism is implemented with Pearlmutter’s $$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$2-operator,

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$3

This yields an analytical expression for Hessian-vector multiplication, avoids finite-difference instability, and does not require storing a dense Hessian. The practical optimizer HFSGVI uses conjugate gradient to solve

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$4

followed by

$$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$5

The paper states that CG is stopped early, typically at a small fixed number $$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$6 such as $$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$7, even for thousands of parameters; the appendix notes that preconditioned CG is used, but does not specify the preconditioner.

A central complexity claim is that the reparameterized second-order formulas have the same order in $$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$8 as first-order stochastic backpropagation, namely $$\hat{\theta}=\arg\max_\theta \mathcal L_\theta, \qquad \mathcal L_\theta = \sum_{i=1}^n \log p(y_i\mid f_\theta(x_i))+\log p(\theta),$$9 with respect to latent dimension. One stochastic gradient iteration costs $\hat\theta$0, and one HF step with at most $\hat\theta$1 CG iterations costs $\hat\theta$2. The paper also proves a dimension-free variance bound for Gaussian Monte Carlo: if $\hat\theta$3 is $\hat\theta$4-Lipschitz differentiable and $\hat\theta$5, then

$\hat\theta$6

and then notes that the proof can be sharpened to

$\hat\theta$7

Empirically, the paper studies Bayesian logistic regression on DukeBreast, Leukemia, and a9a, and VAE training on Frey Face, Olivetti Face, and MNIST. On DukeBreast and Leukemia, HFSGVI converges within 3 iterations; on a9a, L-BFGS-SGVI performs best overall in convergence rate and final lower bound. In the VAE experiments with $\hat\theta$8, HFSGVI and L-BFGS-SGVI converge faster than Ada in wall-clock time, and Ada takes “at least four times as long” to achieve similar lower bound. These results support curvature-based variational optimization, but they do not constitute self-normalized inference in the SNIS sense.

5. Approximate Newton covariance inference from gradients alone

“Approximate Newton-based statistical inference using only stochastic gradients” addresses a more classical inferential target: covariance estimation for empirical-risk minimizers without exact second-order information (Li et al., 2018). In low-dimensional convex $\hat\theta$9-estimation, the asymptotic distribution is

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$0

with

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$1

The paper’s key reduction is that

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$2

so inference can proceed if one can repeatedly approximate inverse-Hessian-conditioned gradients.

The algorithm does this through a stochastic Newton subproblem solved only with gradients. The outer loop forms

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$3

and the inner loop updates

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$4

with

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$5

The finite-difference approximation

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$6

is the paper’s Hessian-free device. Averaging the inner iterates,

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$7

the method uses $$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$8 for statistical inference, and estimates the covariance by

$$P = - \nabla\nabla_\theta \mathcal L_{\hat\theta}, \qquad q(\theta)=\mathcal N(\hat\theta, P^{-1}).$$9

A main theorem gives

$$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$0

The paper extends this framework in two directions. For high-dimensional sparse linear regression, it introduces a modified strongly convex objective based on a soft-thresholded covariance estimate $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$1, then defines the de-biased estimator

$$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$2

For dependent time-series data, it replaces iid outer minibatches with contiguous blocks, producing a Hessian-free HAC procedure whose covariance matches a Newey-West estimator with block-induced weights. The method is Hessian-free in derivative order and storage; in the high-dimensional implementation, the paper states $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$3 space rather than $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$4, though the arithmetic cost can still be large.

The inferential validation is explicit. In low-dimensional linear and logistic regression, the paper reports confidence interval coverage and lengths close to bootstrap and inverse-Fisher baselines, while averaged SGD intervals perform poorly. In high-dimensional sparse regression, two-sided $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$5-test p-values from the de-biased estimator are described as close to uniform. For time series with moving-average noise, average 95% interval coverage is reported around $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$6. These are covariance-estimation and studentization results, not Monte Carlo self-normalization in the importance-sampling sense.

6. Common misconceptions, limitations, and precise relevance of the label

A frequent misconception is that “Hessian-free” means the Hessian is absent from the theory. The available papers do not support that interpretation. In HFL, the derivation explicitly relies on the Hessian of the log posterior and on the existence of $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$7, even though computation avoids explicit matrix operations (McInerney et al., 2024). In second-order variational inference, the Hessian is derived exactly and then accessed through Hessian-vector products (Fan et al., 2015). In approximate Newton inference, the target covariance is still the sandwich form $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$8, with the Hessian action approximated by gradient differences (Li et al., 2018).

A second misconception is that these methods instantiate self-normalized importance sampling. They do not. The variational-inference paper uses ordinary sample averages under a Gaussian base law, not normalized weights (Fan et al., 2015). HFL uses a normalized perturbation size $$\operatorname{Var}\big(\tilde f_\theta(x)\big) = \nabla_\theta f_{\hat\theta}(x)^\top P^{-1}\nabla_\theta f_{\hat\theta}(x).$$9, but the paper does not present a formal self-normalization framework or classwise normalization scheme (McInerney et al., 2024). The approximate Newton paper is closest to self-normalization only in a studentization-like sense, because the covariance used for inference is estimated from the same stochastic gradient machinery that generates the inferential replicates (Li et al., 2018).

The limitations are equally distinct. HFL inherits the usual Laplace assumptions: locality around one optimum, neglect of multiple modes, and sensitivity to optimization and local geometry. Its plain form also requires a separate perturbed solve for each query, with pre-trained HFL introduced as an amortized approximation (McInerney et al., 2024). The second-order variational-inference theory is formulated for Gaussian variational families and smooth objectives, and the method is an optimizer for variational objectives rather than a direct uncertainty-calibration framework (Fan et al., 2015). The approximate Newton theory is strongest for convex, strongly convex, twice-differentiable objectives with Lipschitz Hessians; nonconvex applications are empirical rather than theoretically covered, and the high-dimensional theory is specialized to sparse Gaussian linear models with sparse, diagonally dominant covariance (Li et al., 2018).

Taken together, the cited works support a precise but narrow interpretation. The strongest direct reading of “Hessian-Free Self-Normalized Inference” is not a unified named method, but a family resemblance among procedures that recover uncertainty-relevant curvature effects without explicit Hessian algebra and that normalize those effects either by perturbation size or by an estimated covariance operator. Within that space, HFL is the clearest example on the Bayesian deep-learning side, because it shows that epistemic predictive variance can be recovered from the sensitivity of the optimum to a tiny prediction-dependent regularization perturbation, and in the infinitesimal limit that sensitivity equals

$P^{-1}$0

That identity is the sharpest available bridge between Hessian-free computation and normalized uncertainty estimation in the sources considered here (McInerney et al., 2024).