PathFinder: Quasi-Newton Variational Inference

• PathFinder is a quasi-Newton variational inference method that converts optimization trajectories into local Gaussian approximations of the target posterior.
• It leverages curvature information via inverse-Hessian estimates and selects the best approximation by minimizing the Monte Carlo KL divergence (maximizing the ELBO).
• Multi-path aggregation with Pareto-smoothed importance resampling boosts diversity and efficiency, making it a robust initializer for MCMC.

PathFinder is a fast, parallelizable variational inference method that turns quasi-Newton optimization trajectories into a sequence of multivariate normal approximations to a target posterior, then selects the best approximation by minimizing Monte Carlo Kullback–Leibler divergence. It was developed for differentiable log densities on an unconstrained support, with constrained parameters handled by transformation to $\mathbb{R}^N$ and Jacobian adjustments in the log density. In the formulation reported in "Pathfinder: Parallel quasi-Newton variational inference" (Zhang et al., 2021), the method is intended both to provide approximate posterior draws with substantially fewer gradient and log-density evaluations than automatic differentiation variational inference (ADVI) and short dynamic Hamiltonian Monte Carlo (HMC) runs, and to serve as a robust initializer for MCMC.

1. Conceptual basis and operating assumptions

PathFinder’s core idea is to run a quasi-Newton optimizer, specifically BFGS or L-BFGS, on the log posterior and to reinterpret points along the optimization path as a sequence of local Gaussian approximations. At each iteration, the optimizer supplies not only a location but also an inverse-Hessian approximation, and that inverse-Hessian approximation is used as the covariance estimate of a local multivariate normal approximation. The method then evaluates the ELBO, equivalently the KL divergence, for those local approximations and chooses the approximation with the lowest estimated KL.

This design sharply distinguishes PathFinder from ADVI. ADVI directly optimizes the stochastic ELBO with small noisy steps and typically diagonal or dense covariance, whereas PathFinder optimizes the target log density itself with curvature-informed steps and only evaluates the ELBO at a modest number of points along the path. Those ELBO evaluations are embarrassingly parallel. Relative to short dynamic HMC chains, PathFinder avoids serial sampler adaptation and uses far fewer log-density and gradient calls, while returning a normal approximation whose draws are typically competitive with the last draws of short HMC warmup.

The method assumes a differentiable log density up to a normalizing constant on an unconstrained support. When parameters are constrained, they are transformed to $\mathbb{R}^N$ using, for example, logit, log, or Cholesky transforms, and the change-of-variables adjustment is applied to the log density. These assumptions make PathFinder a normal variational approximation rather than an asymptotically exact sampler, a point that is central to understanding both its speed and its limitations (Zhang et al., 2021).

2. Quasi-Newton geometry and variational selection

PathFinder ascends the log density with a quasi-Newton step

$\theta_{k+1} = \theta_k + H_k g_k,$

where $g_k = \nabla_\theta \log p(\theta_k)$ and $H_k$ is an inverse-Hessian approximation of $-\nabla^2 \log p$. In BFGS inverse form,

$$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$

with

$$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$

L-BFGS stores only a recent history of pairs $\{(s_j, y_j)\}$ of size $J$, so the inverse-Hessian estimate is represented in low-rank plus diagonal form rather than as a dense matrix.

At iteration $\mathbb{R}^N$0, the local approximation is

$\mathbb{R}^N$1

with

$\mathbb{R}^N$2

The corresponding Gaussian log density is

$\mathbb{R}^N$3

Selection among path points is based on Monte Carlo estimation of

$\mathbb{R}^N$4

Equivalently, PathFinder maximizes

$\mathbb{R}^N$5

using the same estimate, and chooses

$\mathbb{R}^N$6

The method therefore reframes variational inference as approximation selection along a curvature-informed optimization trajectory rather than direct stochastic optimization of a single variational family (Zhang et al., 2021).

3. Single-path algorithm and covariance factorization

The single-path procedure begins with random initialization on the unconstrained space, for example $\mathbb{R}^N$7, followed by computation of $\mathbb{R}^N$8. L-BFGS is then run up to $\mathbb{R}^N$9 iterations with history size $\theta_{k+1} = \theta_k + H_k g_k,$0 and convergence tolerance $\theta_{k+1} = \theta_k + H_k g_k,$1, while maintaining the recent update pairs $\theta_{k+1} = \theta_k + H_k g_k,$2 and $\theta_{k+1} = \theta_k + H_k g_k,$3.

For each path point, PathFinder reconstructs an inverse-Hessian approximation in low-rank plus diagonal form,

$\theta_{k+1} = \theta_k + H_k g_k,$4

with the diagonal preconditioner $\theta_{k+1} = \theta_k + H_k g_k,$5 recovered adaptively through “$\theta_{k+1} = \theta_k + H_k g_k,$6-recover,” and pairs implying extreme curvature filtered out by positive-definite safeguards. The local normal approximation is then centered at

$\theta_{k+1} = \theta_k + H_k g_k,$7

with covariance $\theta_{k+1} = \theta_k + H_k g_k,$8.

ELBO evaluation is performed in parallel over path points. For each $\theta_{k+1} = \theta_k + H_k g_k,$9, the method draws $g_k = \nabla_\theta \log p(\theta_k)$0 samples from $g_k = \nabla_\theta \log p(\theta_k)$1 and evaluates

$g_k = \nabla_\theta \log p(\theta_k)$2

The maximizing index $g_k = \nabla_\theta \log p(\theta_k)$3 determines the returned approximation, from which $g_k = \nabla_\theta \log p(\theta_k)$4 approximate posterior draws are produced.

A distinctive implementation feature is the Cholesky-free sampler for the factored L-BFGS covariance. Using a factorization $g_k = \nabla_\theta \log p(\theta_k)$5, a draw is obtained by sampling $g_k = \nabla_\theta \log p(\theta_k)$6, setting $g_k = \nabla_\theta \log p(\theta_k)$7, and then

$g_k = \nabla_\theta \log p(\theta_k)$8

The same factorization yields

$g_k = \nabla_\theta \log p(\theta_k)$9

and therefore efficient evaluation of $H_k$0. The inner sampling and $H_k$1 calculations scale as $H_k$2, with memory footprint $H_k$3 per trajectory point; the key point is that $H_k$4 is the history size, not the dimension (Zhang et al., 2021).

4. Multi-path aggregation, importance resampling, and parallel structure

A single L-BFGS trajectory may terminate in a minor mode, on a plateau, or near a saddle point. PathFinder addresses this with a multi-path variant that runs $H_k$5 independent single-path instances in parallel, each producing a selected $H_k$6 and associated draws. These selected normals are then treated as an equally weighted mixture proposal,

$H_k$7

against the augmented target

$H_k$8

The approximate posterior sample is refined by Pareto-smoothed importance resampling. For candidate draws $H_k$9, raw weights are computed as

$-\nabla^2 \log p$0

then stabilized with PSIS, normalized, and resampled with replacement. In the single-path or pooled setting, the same principle can be written as

$-\nabla^2 \log p$1

This multi-path design is important for two reasons. First, it improves diversity and robustness by down-weighting minor-mode trajectories and reducing failures from optimization pathologies. Second, it amplifies the method’s computational advantage because the expensive pieces are embarrassingly parallelizable: ELBO estimation across trajectory points, independent single-path runs, and target log-density evaluations in the PS-IR step. Only the L-BFGS trajectory generation inside one run remains serial, and that serial segment is typically fast because curvature-informed steps require few evaluations (Zhang et al., 2021).

5. Empirical behavior and comparative performance

The reported empirical study covers 20 posterior distributions drawn from models including generalized linear models, hierarchical meta-analysis, Gaussian processes, mixtures, ODE models, HMMs, and ARMA, AR, and GARCH time-series. Approximation quality is measured by empirical $-\nabla^2 \log p$2-Wasserstein distance to long HMC reference draws, with the Kantorovich–Rubinstein dual written as

$-\nabla^2 \log p$3

Across this benchmark, single-path PathFinder typically outperformed mean-field and dense ADVI, and in many models ADVI’s median $-\nabla^2 \log p$4 was at least twice that of PathFinder. PathFinder’s $-\nabla^2 \log p$5 was comparable to, and often slightly worse than, short dynamic HMC chains, while multi-path PathFinder narrowed or reversed that gap. The reported hidden Markov model is a notable counterexample in which mean-field ADVI performed very well because stochastic noise escaped minor modes; in that setting, multi-path PathFinder eliminated failures of single-path PathFinder by aggregating multiple paths.

The efficiency results are equally central. Averaged across the 20-model posteriordb set, Stan dynamic HMC phase I required approximately $-\nabla^2 \log p$6 more log-density evaluations and approximately $-\nabla^2 \log p$7 more gradient evaluations than PathFinder. Mean-field ADVI required approximately $-\nabla^2 \log p$8 more log-density and approximately $-\nabla^2 \log p$9 more gradient evaluations, while dense ADVI required approximately $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$0 more log-density and approximately $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$1 more gradient evaluations. The paper emphasizes that these differences grow for more challenging posteriors, including plateaus, funnels, and multimodality, where curvature information is especially valuable.

The case studies clarify the geometric behavior of the method. In Neal’s funnel and the centered eight-schools model, PathFinder often selects an approximation in the high-probability mass region before the optimizer dives down a pole of unbounded density. In Gaussian-process models with non-normal posterior shapes, the selected approximations can be conservative and underdispersed, yet still useful as initializations. In multimodal settings, multiple independent paths plus Pareto-smoothed importance resampling improve diversity and reduce $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$2 further (Zhang et al., 2021).

6. Tuning, diagnostics, applications, and limitations

The defaults reported in the experiments are $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$3, $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$4, $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$5, $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$6, and $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$7 for single-path PathFinder, and $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$8 with $$H_{k+1} = (I - \rho_k s_k y_k^{\top}) H_k (I - \rho_k y_k s_k^{\top}) + \rho_k s_k s_k^{\top},$$9 PS-IR draws for the multi-path variant. Sensitivity analyses are explicit. Reducing $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$0 or loosening $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$1 degraded approximation stability and fidelity. Increasing $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$2 reduced ELBO variance but produced only modest $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$3 improvements; $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$4 typically sufficed, while $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$5 could help when ELBO selection near the optimum required more precision. Increasing $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$6, for example from $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$7 to $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$8, could help capture stronger posterior dependencies at higher computational cost. Increasing $$s_k = \theta_{k+1} - \theta_k,\qquad y_k = g_{k+1} - g_k,\qquad \rho_k = \frac{1}{y_k^{\top} s_k}.$$9 beyond $\{(s_j, y_j)\}$0 yielded modest gains; decreasing $\{(s_j, y_j)\}$1 to $\{(s_j, y_j)\}$2 slightly worsened $\{(s_j, y_j)\}$3 and increased variability.

The method is especially useful when fast approximate posterior draws are needed under a tight computational budget, when ADVI is unstable, when the posterior benefits from curvature-informed exploration, or when embarrassingly parallel inference on multi-core hardware is desirable. It integrates naturally with probabilistic programming systems that expose differentiable log densities and automatic differentiation, and the implementation described in the paper used R’s stats::optim() with L-BFGS-B. A common practical use is as an MCMC initializer: PathFinder draws can eliminate much of transient bias, reduce the risk of chains getting stuck in minor modes, and skip most warmup. In a Gaussian-process “birthday problem” with 429 parameters, multi-path PathFinder initialization led to better traceplots and predictions and lower $\{(s_j, y_j)\}$4 than short adaptive HMC alone.

Its limitations are equally clear. Severe multimodality can defeat a single L-BFGS trajectory, and even multi-path PathFinder cannot recover if all runs find minor modes. Non-identifiability, plateaus, and saddle points can make curvature estimates unstable, yielding underdispersed or miscentered approximations. Because PathFinder is a normal variational method, far-from-normal posteriors generally lead to conservative, concentrated approximations, and mixture aggregation improves but does not eliminate that limitation. Diagnostics therefore matter: Pareto-$\{(s_j, y_j)\}$5 values from PSIS indicate whether importance weights are high-variance, and optimization failures are handled by returning the last point with infinite $\{(s_j, y_j)\}$6 so that it receives zero weight in multi-path resampling. These features underscore a recurring misconception: PathFinder is not a replacement for asymptotically exact sampling when exactness is required, but rather a fast curvature-aware approximation and initialization scheme whose strengths and failure modes are sharply defined (Zhang et al., 2021).