Variational Laplace in Bayesian Inference

• Variational Laplace is a Bayesian inference technique that couples variational optimization with Gaussian second-order approximations to derive local posterior estimates.
• It systematically applies the Laplace method within variational frameworks to handle conjugate and nonconjugate models, enhancing estimation in hierarchical, sparse, and neural-network settings.
• Empirical results show improved model evidence approximation, accurate uncertainty quantification, and efficient parameter recovery across various applications.

Searching arXiv for recent and foundational papers on “Variational Laplace” and closely related formulations. Variational Laplace denotes a family of approximate Bayesian inference schemes in which variational optimization is combined with a local Gaussian, second-order approximation of the posterior; in Daunizeau’s formulation, variational-Laplace or VL schemes rely on Gaussian approximations to posterior densities on model parameters (Daunizeau, 2017). Closely related usages include Laplace variational inference for nonconjugate mean-field updates (Wang et al., 2012), hybrid Laplace–variational Gaussian approximations that keep Hessian structure and variationally refine a restricted subset of parameters (Gianniotis, 2019, Niekerk et al., 2021), and modern neural formulations that use Laplace-based posteriors over network weights or latent states (Unlu et al., 2020, Park et al., 2022). The term also appears, with a different meaning, in variational problems for Laplace transforms, Laplace eigenvalues, and $p$-Laplace equations (Ustunel, 2014, Kokarev, 2011, Hu et al., 27 May 2026).

1. Variational and Laplace components

In the canonical Bayesian formulation, VL starts from the variational free energy

$$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$

with equivalent decompositions

$$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$

and

$$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$

Introducing the variational energy

$$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$

the unrestricted optimum satisfies $q(\theta)\propto \exp(I(\theta))$; VL imposes the Gaussian restriction

$q(\theta)\approx \mathcal N(\mu,\Sigma),$

with

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$

At convergence, the free energy reduces to the familiar Laplace form

$$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$

so the Gaussian posterior is simultaneously a local second-order approximation and a constrained variational optimum (Daunizeau, 2017).

A closely related but not identical construction appears in nonconjugate mean-field variational inference. There the intractable coordinate update for a nonconjugate factor is written

$q(\theta)\propto \exp\{f(\theta)\},$

and Laplace variational inference replaces $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$0 by a second-order Taylor expansion around its maximizer $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$1, yielding

$$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$2

This use of Laplace occurs inside each variational coordinate update rather than as a one-shot approximation to the full posterior, but the operational idea is the same: a local Gaussian is extracted from curvature of a log-density-like objective (Wang et al., 2012).

2. Continuous, categorical, and hierarchical formulations

For nonlinear models with conditionally Gaussian observations,

$$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$3

and Gaussian prior $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$4, the variational energy becomes the sum of a Gaussian data-fit term and a Gaussian prior term. If $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$5 and higher derivatives of $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$6 are neglected, the posterior covariance is the Gauss–Newton expression

$$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$7

The corresponding free energy separates into an accuracy term involving the residual $$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$8 and a complexity term involving both prior deviation and posterior contraction relative to the prior (Daunizeau, 2017).

The same calculus extends to categorical data. For Bernoulli observations,

$$F(q)=\log p(y\mid m)-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid y,m)\right),$$9

the gradient and Hessian of $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$0 can be written explicitly, and for logistic regression

$$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$1

the curvature simplifies to

$$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$2

Analogous formulas are derived for binomial and multinomial likelihoods, with the covariance again obtained from a negative inverse Hessian or its first-order approximation. This places logistic regression, multinomial regression, and nonlinear Gaussian models inside a single VL template (Daunizeau, 2017).

Hierarchical extensions introduce hyperparameters $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$3 controlling observation and prior precision through linear precision-component decompositions such as

$$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$4

With the mean-field factorization $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$5, the parameter block is updated by replacing fixed precisions with their variational expectations, while Gamma hyperpriors

$$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$6

yield closed-form Gamma updates for $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$7 and $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$8. The resulting free energy acquires an additional hyperparameter correction $$F(q)=\langle \log p(y\mid \theta,m)\rangle_{q(\theta)}-D_{\mathrm{KL}}\!\left(q(\theta);p(\theta\mid m)\right)$$9, so model evidence approximation remains available in the hierarchical case (Daunizeau, 2017).

3. Nonconjugate and hybrid Gaussian approximations

Laplace variational inference in nonconjugate models systematizes the Gaussian coordinate update

$$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$0

for models of the form $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$1, where $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$2 is nonconjugate and $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$3 is conjugate. The approach is generic, requires only differentiability of $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$4, and recovers the standard Laplace approximation in Bayesian logistic regression, where there is effectively no additional latent block $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$5 (Wang et al., 2012).

Several later methods use this same Gaussian-curvature core but restrict how much of the Gaussian is variationally adapted. Mixed Variational Inference begins from the Laplace posterior

$$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$6

and then optimizes a stochastic variational objective over reduced Gaussian families: mean only in $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$7, mean plus axis rescaling in $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$8, or mean plus a rank-1 covariance-root update in $$F(q)=\left\langle \log p(y\mid \theta,m)+\log p(\theta\mid m)\right\rangle_{q(\theta)}+S(q(\theta)).$$9. The goal is to retain posterior correlations from the Laplace covariance while avoiding the full $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$0 parameter burden of unrestricted Gaussian VI (Gianniotis, 2019).

Low-Rank Variational Bayes Correction uses a different restriction. Starting from the Laplace precision $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$1 at the posterior mode $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$2, it keeps the covariance $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$3 fixed and variationally corrects only the mean: $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$4 where $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$5 collects selected columns of the inverse precision and $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$6 is low-dimensional. The reduced objective is a KL/ELBO/Zellner variational criterion over this affine subspace, so the method is explicitly a variational refinement of a Laplace Gaussian approximation in a reduced subspace (Niekerk et al., 2021).

A more corrective use of Laplace appears in differential-equation models. In Laplace-aided variational inference for ODE systems, mean-field state-space variational Bayes gives fast point estimates but underestimates posterior covariance; the method therefore treats the SSVB mean as an approximate mode and uses the inverse Hessian of the relaxed or exact ODE posterior as a covariance repair. This is not a full Gaussian variational Laplace scheme, but rather mean-field VB plus a post hoc Laplace covariance correction (Yang et al., 2021).

4. Sparse and structured priors inside variational-Laplace

A distinctive structured extension is the use of parameter transformations to emulate sparse priors within standard VL/VBA machinery. For a generative model

$$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$7

Daunizeau proposes the signed quadratic remapping

$$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$8

If latent parameters $$I(\theta)=\log p(y\mid \theta,m)+\log p(\theta\mid m),$$9 have i.i.d. Gaussian prior $q(\theta)\propto \exp(I(\theta))$0 and the effective model parameter is $q(\theta)\propto \exp(I(\theta))$1, then $q(\theta)\propto \exp(I(\theta))$2, so the quadratic latent penalty becomes linear in $q(\theta)\propto \exp(I(\theta))$3: $q(\theta)\propto \exp(I(\theta))$4 This produces an $q(\theta)\propto \exp(I(\theta))$5-like sparse prior in the transformed parameterization while keeping Gaussian priors and Gaussian posteriors in latent space (Daunizeau, 2017).

The exact signed-square map is not differentiable at zero, so a smoothed version is introduced,

$q(\theta)\propto \exp(I(\theta))$6

with temperature $q(\theta)\propto \exp(I(\theta))$7. VL then proceeds unchanged except that the likelihood uses $q(\theta)\propto \exp(I(\theta))$8 instead of $q(\theta)\propto \exp(I(\theta))$9. The Gaussian posterior remains

$q(\theta)\approx \mathcal N(\mu,\Sigma),$0

and the free energy

$q(\theta)\approx \mathcal N(\mu,\Sigma),$1

continues to approximate $q(\theta)\approx \mathcal N(\mu,\Sigma),$2 (Daunizeau, 2017).

The empirical setting is an underdetermined Gaussian linear model with $q(\theta)\approx \mathcal N(\mu,\Sigma),$3 and $q(\theta)\approx \mathcal N(\mu,\Sigma),$4. In sparse simulations, the sparse-transformed model outperformed the ordinary Gaussian-prior model and the model comparison favored sparse priors by $q(\theta)\approx \mathcal N(\mu,\Sigma),$5; in Gaussian simulations, where sparsity was the wrong inductive bias, the ordinary Gaussian-prior model won with $q(\theta)\approx \mathcal N(\mu,\Sigma),$6. The systematic Monte Carlo study further showed that relative evidence $q(\theta)\approx \mathcal N(\mu,\Sigma),$7 and relative recovery accuracy $q(\theta)\approx \mathcal N(\mu,\Sigma),$8 both increase with true sparsity and data precision, while support recovery measured by TPR and TNR is generally above chance but conservative (Daunizeau, 2017).

5. Neural-network and latent-state formulations

In Bayesian neural networks, VL has been formulated as a deterministic approximation to the ELBO for a factorized Gaussian posterior

$q(\theta)\approx \mathcal N(\mu,\Sigma),$9

A local quadratic expansion of the expected log-likelihood around the posterior mean, followed by a Fisher approximation to the negative Hessian, yields the objective

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$0

which the paper summarizes as the log-likelihood, plus weight-decay, plus a squared-gradient regularizer. On CIFAR-10, CIFAR-100, SVHN, and Fashion-MNIST, this deterministic VL objective gave better test performance and expected calibration errors than MAP inference and standard sampling-based VI, despite using the same variational approximate posterior (Unlu et al., 2020).

In deep generative models, Variational Laplace Autoencoders replace the standard encoder posterior

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$1

by a Laplace Gaussian

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$2

where $$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$3 is a refined posterior mode. For ReLU decoders with Gaussian output, local piecewise linearity yields the PPCA-like update

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$4

optionally damped. The method addresses both limited posterior expressiveness and amortization error, and on MNIST, Omniglot, Fashion-MNIST, SVHN, and CIFAR-10 it outperformed VAE, SA-VAE, VAE+HF, and VAE+IAF in the reported settings (Park et al., 2022).

Laplacian Autoencoders move the Laplace approximation to weight space inside an autoencoder. The parameter posterior is approximated online by

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$5

the network is linearized in parameter space,

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$6

and the precision update is

$$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$7

A new Hessian approximation scales linearly with data size, and the reported outcome is well-calibrated uncertainties in both latent and output space together with improved downstream performance (Miani et al., 2022).

A sequential latent-state version appears in Bayesian meta-reinforcement learning with Laplace variational recurrent networks. There the RNN hidden state $$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$8 is treated as a local MAP estimate for a latent task variable $$\mu=\arg\max_\theta I(\theta),\qquad \Sigma^{-1}=-\left.\frac{\partial^2 I}{\partial \theta^2}\right|_{\theta=\mu}.$$9, and Laplace is applied in latent space: $$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$0 In the Gaussian special case, the precision reduces to a sum of Jacobian outer products,

$$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$1

This approximation can be attached before, during, or after learning without changing the base architecture, performs on par with variational baselines while using much fewer parameters, and exposes the overconfidence of point-estimate recurrent meta-RL agents (Vries et al., 24 May 2025).

6. Asymptotics, limitations, and broader meanings

Recent theory has separated two questions: the statistical asymptotics of VL approximations and the computational asymptotics of the underlying optimization problems. In a frequentist analysis of variational Laplace, point estimates generated by VL were shown to enjoy asymptotic consistency and efficiency in two toy examples, and sufficient conditions were derived for these properties in a general setting. The same work also studies convergence of the associated Gaussian approximations in total variation distance, thereby linking VL to both classical Bernstein–von Mises considerations and modern asymptotic analyses of variational inference (Keck, 23 Jul 2025). A complementary line of work analyzes Gaussian variational inference and the Laplace approximation as nonconvex optimization problems whose global optima become locally strongly convex in the asymptotic regime; this yields Consistent Laplace Approximation and Consistent Stochastic Variational Inference, both initialized by a smoothed posterior mode and both shown empirically to improve the likelihood of obtaining the global optimum of their respective objectives (Xu et al., 2021).

The limitations emphasized across the literature are similarly consistent. VL is a local Gaussian approximation and is therefore sensitive to posterior multimodality, skewness, and heavy tails. Hybrid methods that freeze the Laplace covariance or confine mean updates to low-rank subspaces are designed mainly to correct location bias rather than covariance pathology (Niekerk et al., 2021). Mean-field or blockwise Gaussian variants can remain underdispersed even when point estimates are asymptotically efficient (Keck, 23 Jul 2025). In the sparse-prior construction, the nonlinear remapping can make optimization harder and more local-mode prone than an ordinary Gaussian-prior GLM, and support recovery is conservative in poorly conditioned problems (Daunizeau, 2017). In Bayesian neural networks, the curvature-based objective requires second derivatives, so the paper explicitly uses smooth activations rather than ReLU because piecewise linear activations have pathological second derivatives (Unlu et al., 2020).

Outside approximate Bayesian inference, the phrase acquires broader meanings. In spectral geometry, the variational problem is to maximize $$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$2 within a conformal class, and Kokarev reformulates conformal metrics as Radon measures to apply the direct method of the calculus of variations to Laplace eigenvalues on surfaces (Kokarev, 2011). On Wiener space, Üstünel studies the variational representation

$$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$3

and relates equality in the corresponding adapted-shift problem to almost sure invertibility of adapted perturbations of identity (Ustunel, 2014). In nonlinear PDEs, the dual variational neural network for the $$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$4-Laplace problem decomposes the flux through an $$F(q)\approx I(\mu)+\frac12\log|\Sigma|+\frac{n_\theta}{2}\log(2\pi),$$5-Helmholtz decomposition into a linear Poisson problem and a convex minimization over divergence-free fields, so even there “variational Laplace” names a variational treatment of Laplace-type structure rather than a single canonical algorithm (Hu et al., 27 May 2026).