Natural-Gradient Variational-Gaussian Inference

• Natural-gradient variational-Gaussian inference is a Bayesian method that approximates posteriors using Gaussian distributions optimized via Fisher information-based natural gradients.
• It employs geometric update rules, such as mirror descent and quasi-second-order steps, to improve convergence and efficiency in complex, high-dimensional models.
• Techniques like Metric Gaussian Variational Inference extend this framework to correlated posteriors and large-scale problems, balancing accuracy with computational speed.

Natural-gradient variational-Gaussian inference denotes a family of variational Bayesian methods in which the posterior is approximated by a Gaussian distribution and optimized in the information geometry induced by the Fisher information rather than in an ambient Euclidean parameter space. In this setting, the same optimization problem can be written as Fisher-preconditioned descent in natural parameters, mirror descent in expectation parameters, or structured updates of means and covariances or precisions; standardized-coordinate variants such as Metric Gaussian Variational Inference extend the framework to correlated posteriors in very high dimensions by replacing explicit covariance parametrization with implicit Fisher-metric operators (Knollmüller et al., 2019, Sun et al., 22 Oct 2025).

1. Information-geometric formulation

For exponential-family variational distributions, natural-gradient variational inference minimizes the negative ELBO

$$\ell(q) = - \mathbb{E}_{q(z)}[\log p(D \mid z)] + \mathrm{KL}(q(z)\|p(z)),$$

or equivalently maximizes the ELBO, using the Fisher information of the variational family as the metric tensor. In natural parameters $\eta$, the canonical update is

$$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$

with $A$ the log-partition function. The same step admits the dual expectation-parameter form

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$

which identifies NGVI with stochastic mirror descent under the mirror map $A^*$ (Sun et al., 22 Oct 2025).

This equivalence is central because, for exponential families, the Bregman divergence induced by $A^*$ coincides with a KL divergence between variational distributions: $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$ Consequently, the natural gradient is the steepest descent direction measured directly in distribution space rather than in a coordinate chart. A further consequence is reparameterization invariance: if $\phi=\phi(\theta)$ is a smooth change of coordinates, then the natural-gradient direction transforms covariantly, whereas the Euclidean gradient does not (Barfoot, 2020).

For Gaussian variational families, this geometric interpretation becomes especially concrete. The local Hessian of the KL between nearby Gaussian members coincides with the Fisher information matrix, so natural-gradient descent acts as a second-order method on the manifold of Gaussians. This is why NGVI is repeatedly described as a quasi-second-order method in the cited literature, but its curvature object is the Fisher geometry of the approximating family or, in MGVI-like constructions, an approximate posterior Fisher metric rather than the raw Hessian of the log posterior (Barfoot, 2020, Knollmüller et al., 2019).

2. Gaussian parametrizations and exact natural-gradient maps

The multivariate Gaussian family admits several equivalent parametrizations. In mean-covariance form,

$$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$

while in mean-precision form one writes $\eta$0. The exponential-family natural parameters are

$\eta$1

In mean-field Gaussian analyses, expectation parameters are often written as $\eta$2 with $\eta$3 and $\eta$4; in full Gaussian treatments they take the analogous form $\eta$5, $\eta$6 (Sun et al., 22 Oct 2025, Yu et al., 2023).

For these parametrizations, the Fisher blocks are available in closed form. In mean-covariance coordinates, the inverse Fisher blocks are $\eta$7 for the mean and $\eta$8 for the covariance block; in mean-precision coordinates, the corresponding block is $\eta$9. The resulting natural-gradient updates reduce to compact matrix formulas: $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$0 for an objective $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$1 such as the negative ELBO (Barfoot, 2020).

A notable consequence is that the $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$2 parametrization is often computationally attractive. One explicit recommendation is the hybrid update

$$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$3

derived for the KL-like objective studied in the Gaussian natural-gradient note. This form automatically respects symmetry and preserves sparsity when the mean Hessian is sparse (Barfoot, 2020).

Positive-definite constraints motivate Cholesky and square-root parametrizations. If $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$4 or $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$5 with lower-triangular $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$6 or $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$7, analytic natural-gradient directions can be written in factor coordinates without explicitly inverting the full Fisher matrix. For covariance factors, one has

$$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$8

with an analogous form for precision factors. These updates preserve triangular structure and, with positive diagonal parametrization, maintain SPD constraints while also encoding sparsity patterns such as block-diagonal or sparse-precision structures (Tan, 2021).

A related development uses the square-root covariance parameter $$\Delta \eta = -\gamma\, I(\eta)^{-1}\nabla_\eta \ell(\eta), \qquad I(\eta)=\nabla^2 A(\eta),$$9 directly. Under strong convexity and smoothness of the penalized loss, the expected loss is convex jointly in $A$0, whereas it is generally not convex jointly in $A$1 or in natural or expectation parameters. This leads to the square-root natural-gradient flow

$A$2

and its discrete-time SR-VN update

$A$3

which preserves positive semidefiniteness by construction (Kumar et al., 10 Jul 2025).

3. Correlated large-scale inference and Metric Gaussian Variational Inference

Metric Gaussian Variational Inference is a full-covariance NGVI method designed to go beyond mean-field while retaining linear scaling in computational time and memory. Its key construction is a standardized-coordinate representation of a hierarchical model: original parameters $A$4 are written as $A$5 with $A$6, obtained by reparameterizing priors through a multivariate distributional transform and then mapping to a standard normal. In these coordinates,

$A$7

and the information Hamiltonian is

$A$8

This makes the prior metric trivial and removes support restrictions from the Gaussian approximation (Knollmüller et al., 2019).

MGVI approximates the posterior covariance by the inverse Fisher information metric evaluated at the current mean. If $A$9 and $$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$0 is the likelihood Fisher information, then the MGVI metric in standardized coordinates is

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$1

with covariance

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$2

Because the Fisher information is positive semidefinite by construction, the posterior metric plus prior identity yields a positive definite precision everywhere; this distinguishes it from Hessian-based Laplace approximations, whose log-posterior Hessian may be indefinite away from modes (Knollmüller et al., 2019).

Optimization alternates between covariance estimation and mean optimization. For a Gaussian $$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$3 in standardized variables, the KL objective with fixed $$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$4 reduces to the cross-entropy, and the mean gradient is

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$5

MGVI then performs a natural-gradient step

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$6

with line search or damping. In practice, the metric is averaged over sample locations, and the natural-gradient direction is obtained by solving an implicit linear system rather than forming matrices explicitly (Knollmüller et al., 2019).

The scaling properties come from an implicit covariance representation. MGVI never stores $$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$7 explicitly; instead it samples from $$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$8 by exploiting the structure

$$\nabla A^*(\omega_{t+1}) = \nabla A^*(\omega_t) - \gamma_t \nabla \ell(\omega_t),$$9

using operator-vector products, draws in data space and parameter space, and conjugate gradients to solve $A^*$0. For independent data, $A^*$1 is diagonal in data space. The paper reports CG iteration counts of $A^*$2–$A^*$3, antithetic sampling for variance reduction, and linear scaling up to problems with about $A^*$4 or “up to a million” model parameters (Knollmüller et al., 2019).

Empirically, MGVI is positioned between mean-field VI and fully explicit full-covariance VI. In Poisson log-normal GP regression, RMS mean error relative to HMC was $A^*$5 for MGVI, versus $A^*$6 for fc-ADVI, $A^*$7 for mf-ADVI, and $A^*$8 for Laplace; RMS standard-deviation error was $A^*$9, $A^*$0, $A^*$1, and $A^*$2, respectively. In hierarchical logistic regression, RMS mean errors were $A^*$3 for MGVI, $A^*$4 for fc-ADVI, $A^*$5 for mf-ADVI, and $A^*$6 for Laplace; RMS standard-deviation errors were $A^*$7, $A^*$8, $A^*$9, and $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$0. The same study reports an order-of-magnitude speed advantage over mf-ADVI in a large binary GP classification problem and marked gains over fc-ADVI in wall-clock cost (Knollmüller et al., 2019).

4. Convergence theory, projections, and inversion-free natural gradients

A major line of recent work studies when NGVI is provably convergent for non-conjugate models. In the mean-field Gaussian case, the variational objective is analyzed in expectation parameters $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$1 with mirror map $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$2, and the algorithm is written as projected stochastic mirror descent: $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$3 The projection is non-Euclidean: standard parameters $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$4 are clipped into a bounded set $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$5 with $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$6 and $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$7, then mapped back to expectation parameters. Under relative smoothness on the bounded domain and a variance condition in mirror geometry, projected NGVI converges globally to a stationary point at $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$8 in the Bregman Forward–Backward Envelope, and under additional log-concavity and hidden-convexity assumptions it achieves $$D_{A^*}(\omega,\omega') = \mathrm{KL}(q_\omega \| q_{\omega'}).$$9 convergence with stochastic gradients or linear convergence in the noiseless case (Sun et al., 22 Oct 2025).

The square-root covariance analysis establishes a complementary result. Under strong convexity $\phi=\phi(\theta)$0, smoothness $\phi=\phi(\theta)$1, and bounded iterates, the KL satisfies a Riemannian Polyak–Łojasiewicz inequality in the Fisher geometry. The continuous-time natural-gradient flow then obeys

$\phi=\phi(\theta)$2

and the discrete-time SR-VN scheme converges at an exponential rate under explicit step-size conditions (Kumar et al., 10 Jul 2025).

Projected stochastic NGVI for general exponential families extends these analyses to different step-size and sample-size schedules. With fixed step size and fixed batch size, the method converges geometrically to a neighborhood of the optimum. With decreasing step sizes or increasing sample or batch sizes, it converges to the optimum at rates of the form $\phi=\phi(\theta)$3, possibly with $\phi=\phi(\theta)$4. Under the Linearly Extended Recoverability condition, the target posterior is “close” to the chosen exponential family in the sense that expectations of a larger sufficient statistic are linear in the smaller one; in that case the constrained optimum is a Bregman projection of the ideal parameter, and the objective is both $\phi=\phi(\theta)$5-strongly convex and $\phi=\phi(\theta)$6-smooth relative to $\phi=\phi(\theta)$7 (Guilmeau et al., 1 Apr 2026).

A different theoretical route avoids explicit Fisher construction entirely. The inversion-free natural-gradient method builds an online approximation of the inverse Fisher information through a Sherman–Morrison recursion applied to score outer products, optionally with a stabilizing spectral dither term. The resulting operator converges almost surely to the inverse Fisher, supports matrix-free application in linear memory and time in the number of variational parameters, and yields averaged-iterate convergence of order $\phi=\phi(\theta)$8 together with a central limit theorem (Godichon-Baggioni et al., 2023).

Taken together, these results narrow a persistent gap between empirical practice and theory. They also show that the sharpness of the guarantees depends strongly on parametrization and geometry: mean-field analyses need compactness or projection arguments, square-root analyses require strong convexity and smoothness, and exponential-family analyses rely on Bregman rather than Euclidean geometry (Sun et al., 22 Oct 2025, Kumar et al., 10 Jul 2025, Guilmeau et al., 1 Apr 2026).

5. Structured covariances, mixture families, and black-box or manifold variants

One branch of NGVI replaces direct covariance updates with optimization on matrix manifolds. In precision form, Manifold Gaussian Variational Bayes on the Precision matrix uses the identities

$\phi=\phi(\theta)$9

and updates the SPD precision matrix by a manifold retraction

$$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$0

Because the update is performed on the SPD manifold, positive definiteness is preserved without ad hoc projection. The same framework introduces vector transport for momentum and highlights a computational advantage of the precision parametrization, since it avoids repeated $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$1-type multiplications (Magris et al., 2022).

Another black-box line exploits the exponential-family identity between natural gradients in natural parameters and Euclidean gradients in expectation parameters. In the Quasi Black-box Variational Inference framework, Gaussian updates are written directly in terms of score-function expectations and do not require model gradients with respect to latent parameters or explicit Fisher prescriptions. For full-covariance Gaussians, the updates specialize to

$$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$2

$$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$3

with $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$4. The same paper also gives diagonal specializations and analytic control-variate denominators for variance reduction (Magris et al., 2022).

Gaussian mixture families extend NGVI beyond unimodal Gaussian approximations. A component-wise NGVI scheme for Gaussian mixtures derives natural-parameter updates

$$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$5

and natural-gradient updates for logits of the mixture weights. A later unifying analysis shows that the component-wise natural-gradient updates of VIPS and iBayes-GMM are equivalent, and that practical performance depends instead on sample selection, natural-gradient estimation, step-size adaptation, trust-region enforcement, and whether components are adapted during optimization. The same study reports that per-component sampling and first-order natural-gradient estimation with trust regions outperform mixture-sampling and zero-order estimation in higher dimensions (Mahdisoltani, 2021, Arenz et al., 2022).

A derivative-free extension targets Bayesian inverse problems with inaccessible forward-model gradients. Derivative Free Gaussian Mixture Variational Inference combines Fisher–Rao natural gradients with specialized quadrature rules that are exact for linear maps, affine invariant in continuous time, and covariance-positive for discrete steps $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$6. The covariance update is performed in inverse form, and the paper proves a covariance-positivity proposition for the resulting Gaussian-mixture flow (Che et al., 8 Jan 2025).

In Bayesian neural networks, noisy natural-gradient methods interpret adaptive weight noise as variational Gaussian inference. Natural-gradient ascent with covariance-matched weight noise yields Gaussian variational posteriors whose structure is determined by the curvature approximation: noisy Adam corresponds to fully factorized Gaussian posteriors, noisy K-FAC to matrix-variate Gaussian posteriors, and eigenvalue-corrected noisy K-FAC augments the K-FAC eigenbasis with a full diagonal rescaling $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$7 measured from projected gradient second moments. This corrects the diagonal variances that a pure Kronecker posterior cannot represent (Zhang et al., 2017, Bae et al., 2018).

6. Applications, empirical scope, and limitations

NGVI has become a methodological bridge across several research areas. In non-conjugate sparse Gaussian-process models, natural gradients for Gaussian variational parameters substantially improve wall-clock performance and robustness, especially for ill-conditioned posteriors; the cited study gives a practical example in which ordinary gradients were unusable, and reports an implementation integrated into GPflow (Salimbeni et al., 2018).

In stochastic motion planning, Gaussian Variational Inference Motion Planning formulates planning under uncertainty as KL minimization for a Gaussian trajectory posterior on a sparse factor graph. The method uses natural-gradient updates in $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$8 coordinates, while a dual control-theoretic algorithm, PGCS-MP, solves a proximal covariance-steering problem with closed-form Riccati updates. On 7-DOF WAM tasks with $$q(z\mid \mu,\Sigma) = (2\pi)^{-d/2} |\Sigma|^{-1/2}\exp\!\left\{-\tfrac12 (z-\mu)^\top \Sigma^{-1}(z-\mu)\right\},$$9 support states, reported runtimes are about $\eta$00–$\eta$01 s for GPMP2, $\eta$02–$\eta$03 s for PGCS-MP, and $\eta$04–$\eta$05 s for GVI-MP, with the variational methods returning trajectory distributions rather than single paths (Yu et al., 2023).

For latent diffusion-process inference, a site-based exponential-family parametrization of a Gaussian variational process replaces slow fixed-point procedures with a fast convex optimization scheme “akin to natural gradient descent.” The posterior process is represented by a block-tridiagonal precision over a discretized path, and Gaussian smoothing supplies the required moments. This formulation also yields a more favorable objective for learning model parameters (Verma et al., 2023).

In heterogeneous multi-output Gaussian processes, a fully natural-gradient scheme augments the usual inducing-variable Gaussian posterior with an exploratory Gaussian distribution over hyperparameters and inducing locations. The method was developed for both linear model of coregionalisation and convolution-process constructions, and the reported outcome is better local optima and higher test performance than adaptive gradient methods on both families of models (Giraldo et al., 2019).

The main limitations are consistent across the literature. A single Gaussian approximation is not well suited to multimodal, heavily skewed, or heavy-tailed posteriors, and MGVI explicitly excludes discrete parameters from its standardized-coordinate construction (Knollmüller et al., 2019). Some theoretical guarantees apply only to mean-field Gaussian families or to bounded domains enforced by projection, while extension to full-covariance Gaussian families remains an open direction in the non-conjugate theory (Sun et al., 22 Oct 2025). The strongest convergence results for square-root parametrizations require strong convexity, smoothness, and bounded iterates (Kumar et al., 10 Jul 2025). In MGVI specifically, the interaction of mini-batching with the Fisher-based covariance is identified as an open question, and strong nonlinearity in the generative map can degrade accuracy because curvature terms in the reparameterized Cramér–Rao bound are neglected (Knollmüller et al., 2019).

A persistent misconception is that natural gradients eliminate approximation error. The cited work does not support that claim. Natural gradients improve optimization geometry, conditioning, and invariance properties, but they do not by themselves make a Gaussian family expressive enough for non-Gaussian posteriors. Another common confusion is to equate Fisher-based NGVI with Laplace or Hessian-based Gaussianization. The opposite distinction is repeatedly emphasized: Fisher metrics remain positive semidefinite by construction and can be valid away from modes, whereas Hessians of the log posterior need not define valid covariances except in special regimes (Knollmüller et al., 2019).