Acceleration and its Limits in Gaussian Variational Inference
Jaesung R. Park, Uijeong Jang, Yifan Chen
stat.ML, cs.LG, math.OC
2026-10-15
Exact Gaussian integrals let momentum solve GVI in √κ iterations. With only pointwise derivatives the query lower bound is κ^{1-o(1)}, unless Hessian variation is small.
On an α-strongly convex, β-smooth objective, gradient descent takes O(κ log(1/ε)) iterations, where κ=β/α is the condition number. Nesterov momentum reduces that to O(√κ log(1/ε)). Gaussian variational inference (GVI) is the same kind of optimization: fit q=N(m,Σ) by minimizing F(q)=Eq[V]-H(q). The excess F(q)-F is the reverse KL above the best Gaussian approximation of the target.
When αI⪯∇²V⪯βI, the Gaussian expectation of V inherits those curvature bounds, so GVI in the mean and a covariance factor is a convex composite problem. Exact expectations make accelerated proximal gradient immediate. Working algorithms usually see only pointwise derivatives and estimate the expectation by Monte Carlo. Accelerated stochastic gradients already pay a variance term σ²/(αε), and in GVI that variance depends on the target, the variational parameters, and the estimator. Curvature bounds alone do not settle whether the total number of target queries also scales as √κ.
Write Σ=LL^⊤ with a positive diagonal, or Σ=C² with C symmetric positive definite. The entropy -log det has an explicit prox: diagonals or eigenvalues go through ϕη(y)=(y+√(y²+4η))/2. Momentum can stay in Euclidean coordinates. There is no need to move onto Bures–Wasserstein geometry, which treats distance between Gaussians as a distance between distributions.
With exact expectations, Algorithm 1 takes step size 1/β and momentum (1-κ^{-1/2})/(1+κ^{-1/2}). Theorem 2.1 gives F(qk)-F≤(1-κ^{-1/2})^k Δ0, where Δ0 adds a distance-to-optimum term to the initial gap. The iteration count is ⌈√κ log(Δ0/ε)⌉, one exact energy gradient per step.
Pointwise derivatives do not keep this rate. Theorem 3.1 builds a family in dimension only d=Θ((log κ)²). Each potential cuts a small cap off a large ball along a hidden direction v. Outside the cap, the potential and every derivative of any fixed order match a reference that ignores v. Inside the cap, curvature is pushed up to κ. From N(0,I), every member of the family has objective gap 1+o(1). For a random v, a fixed query hits the cap with probability at most κ^{-1+o(1)}, a small patch on a high-dimensional ball. Miss the cap, and the transcript matches the reference, so the returned Gaussian carries no information about v. A Gaussian chosen without seeing v also cannot, for most directions, stay close to N(0,I) in KL and still thin the tail along v. Any randomized method, even one allowed derivatives of arbitrary fixed order, therefore needs at least κ exp(-C(log κ)^{3/4}) queries, that is κ^{1-o(1)}, to drive the expected gap down to a fixed ε0 uniformly on the family. Every bound of order poly(d)√κ is ruled out. These potentials have large Hessian variation, so they sit outside the hypothesis of the accelerated theorems.
The positive side adds a Lipschitz Hessian. The constant LH is normalized as λH=LH/α^{3/2}. λH=0 means V is quadratic and the target is already Gaussian. Four estimators confine the nonquadratic residual to the variance.
The covariance factor is kept inside a curvature box that contains the optimum in its interior, so the inverse factor stays bounded. All four methods use O(√κ log(2Δ0/ε)) iterations. The query total stays on that accelerated schedule when λH² d²≲ε√κ, which stops the variance term from dominating. Second-order estimators then need only a constant batch. First-order estimators need a batch of size O(d).
| Setting | Metric | Result |
| Exact expectations, accelerated prox | Iterations | ⌈√κ log(Δ0/ε)⌉ |
| Curvature only, any fixed-order pointwise derivatives | Queries to a fixed accuracy | ≥ κ exp(-C(log κ)^{3/4}); d=Θ((log κ)²), initial gap 1+o(1) |
| FD, STL | Gradient queries | O((d√κ+λH² d³/ε) log(2Δ0/ε)), no HVPs |
| Paired | Gradient queries | O((√κ+λH² d²/ε) log(2Δ0/ε)), HVPs cost d times that |
| Plain reparameterization | Variance term | The available bound can scale as κ² and swallow the √κ gain |
Table 1 places Price–STL on the same second-order query order as the paired estimator.
On Gaussian targets, LH=0. Projected STL with batch 64, started at N(0,I), is compared with the same update without momentum, for κ from 10 to 3000. Updates to reach mean KL 10^{-6} grow like √κ with momentum and like κ without. The κ=300 curves are means over five seeds, with 20–80 percentile bands.
The non-Gaussian targets are graph-probit posteriors on five OpenML binary-classification sets, binary labels with a Gaussian prior smoothed by a graph Laplacian. The figure includes Sonar, WDBC, Diabetes, and Banknote. Predictors form a 10-nearest-neighbor graph, and the prior precision is I+τLG with τ=49. The condition number κ=2+τλmax(LG) worsens as τ grows, while LH does not, which separates conditioning from Hessian variation. Seventeen methods, including the five baselines DGG-Prox, DGG-STL, Price-SPGD, SVRGVI, and FB-GVI, run 600 steps from N(0,I) over five seeds. Gaps are measured against a high-accuracy numerical reference. At batch 64, every accelerated variant except STL-proj reaches 10^{-5} of the initial mean gap sooner than its nonaccelerated twin, on all five datasets. STL's batch should grow with dimension, and at batch 64 momentum can be slower. At batch 256 it reaches the paper's 1% threshold sooner and ends at a smaller gap on all five.
Full-covariance GVI is still a standard Bayesian approximation when the output people use is a mean and a covariance. An exact expectation gives the square-root dependence on κ. With only pointwise derivatives, there is a family of strongly log-concave, log-smooth targets, in dimension (log κ) squared and with a constant initial gap, on which no randomized algorithm has a uniform √κ query bound.
Accelerated sampling has recently reached sublinear condition-number dependence for strongly log-concave targets. The same curvature assumption is not enough for GVI when the error is excess reverse KL and the oracle is a pointwise derivative. If the posterior is close to Gaussian, finite differences, paired gradients, or Price–STL, plus momentum, can keep the leading query term on the accelerated side. A clearly nonquadratic posterior has no such guarantee here.
Natural gradient faces the same query lower bound and still often works well in practice. It is not one of the seventeen experimental methods. Euclidean momentum is using composite convexity and variance reduction. Geometry-aware methods are using a different structure. Which one spends fewer queries on a real posterior is not settled here.
The dimension factors are loose. Variance terms for finite differences and STL carry d³, on top of a covariance that already has d² parameters. Existing non-accelerated analyses give iteration bounds linear in d, so in high dimension the polynomial factors can dominate the √κ improvement. Making the Price–STL sampling term linear in d takes extra structure and appears only in the appendix. The sufficient batch sizes also hide absolute constants with four digits.
Nothing in the paper tests whether λH² d²≲ε√κ before a run, and LH is usually unknown on a real posterior. The graph-probit suite is built so LH does not track τ, with τ=49, five datasets, 600 steps, and five seeds. That does not transfer to a Bayesian neural net. The lower bound is a uniform complexity over a worst-case family at one fixed accuracy ε0. An individual target can be much easier. The hard examples have large Hessian variation, so they do not contradict the accelerated upper bounds. The exact energy gradient is an ideal oracle that a nonquadratic V does not provide. Natural gradient, named in the conclusion as the comparison that matters, is missing from the experiments.
The paper discloses that GPT-5.6 Sol and GPT-6 Astra were used to explore proofs, check steps, draft experiment code, and edit the text. The authors say they checked the mathematics and the experimental results.