On Markov chain Monte Carlo methods for tall data

Rémi Bardenet; Arnaud Doucet; Chris Holmes

Source-linked AI summary

On Markov chain Monte Carlo methods for tall data

Rémi Bardenet, Arnaud Doucet, Chris Holmes

arXiv:1505.02827v1May 11, 2015stat.MEstat.COstat.ML

TL;DR

Tall datasets make standard MCMC computationally burdensome, despite cases where full Bayesian uncertainty quantification is valuable. This paper reviews scalable Metropolis-Hastings methods and proposes a subsampling approach that can reduce per-iteration cost in favourable scenarios, while performance remains strongest when the Bernstein-von Mises approximation is excellent.

Problem
MCMC for tall datasets often requires sweeping over all data points at each iteration, limiting practical Bayesian inference despite the need for precise uncertainty quantification.
Method
The paper reviews divide-and-conquer and subsampling approaches, then proposes an original subsampling method targeting a controlled approximation of the posterior under strong ergodicity assumptions.
Results
In favourable scenarios, the proposed methodology lowers the natural O(n) subsample cost per iteration to as low as O(1), and improves on other subsampling approaches when the Bernstein-von Mises approximation is excellent.
Takeaways & Limitations
The approach can make subsampling-based MCMC substantially cheaper for certain models and favourable posterior regimes.
Takeaways & Limitations
Vanilla confidence samplers still require O(n) samples at each iteration at equilibrium, while robustness to non-Gaussianity can also require O(n) subsamples per iteration.

Abstract

from arXiv · show

Markov chain Monte Carlo methods are often deemed too computationally intensive to be of any practical use for big data applications, and in particular for inference on datasets containing a large number $n$ of individual data points, also known as tall datasets. In scenarios where data are assumed independent, various approaches to scale up the Metropolis-Hastings algorithm in a Bayesian inference context have been recently proposed in machine learning and computational statistics. These approaches can be grouped into two categories: divide-and-conquer approaches and, subsampling-based algorithms. The aims of this article are as follows. First, we present a comprehensive review of the existing literature, commenting on the underlying assumptions and theoretical guarantees of each method. Second, by leveraging our understanding of these limitations, we propose an original subsampling-based approach which samples from a distribution provably close to the posterior distribution of interest, yet can require less than $O(n)$ data point likelihood evaluations at each iteration for certain statistical models in favourable scenarios. Finally, we have only been able so far to propose subsampling-based methods which display good performance in scenarios where the Bernstein-von Mises approximation of the target posterior distribution is excellent. It remains an open challenge to develop such methods in scenarios where the Bernstein-von Mises approximation is poor.

1 Introduction

MCMC inference on tall datasets is often too computationally intensive because Metropolis-Hastings requires a full data sweep at every iteration. The paper reviews scalable approaches and introduces an improved subsampling method that can reduce per-iteration likelihood evaluations below O(n) in favourable cases.

Tall datasets make fully Bayesian MCMC inference computationally difficult because Metropolis-Hastings sweeps over all n data points at each iteration.
Frequentist and variational Bayes methods are often preferred computationally, although complex likelihoods and the need for precise uncertainty quantification can favour full Bayesian analysis.
The paper reviews divide-and-conquer and subsampling approaches for scaling Metropolis-Hastings in Bayesian inference with independent data.
The improved confidence sampler targets a controlled approximation of the posterior and can break the O(n) barrier for likelihood evaluations per iteration in favourable cases.Its main limitation is requiring cheap-to-evaluate log-likelihood proxies with known error.

2 Bayesian inference, MCMC, and tall data

The paper formulates Bayesian inference with independent observations and uses Metropolis-Hastings to approximate an intractable posterior. For tall datasets, evaluating the full log-likelihood ratio becomes too costly, motivating batching or partial-data approximations.

The inference setup assumes conditionally independent data with likelihood contributions p(x_i|θ), a prior p(θ), and posterior target π.
Metropolis-Hastings simulates a Markov chain with invariant distribution π and uses ratios of an unnormalized target γ.
When n ≫ 1, evaluating the full log-likelihood ratio is too costly, motivating divide-and-conquer or partial-data approximations.
The running examples fit a one-dimensional Gaussian model to 10^5 Gaussian and lognormal observations, including a misspecified-model case.The algorithms are compared with vanilla MH, which requires n likelihood evaluations per iteration.
Vanilla MH matches a longer reference run in both examples, while the Bernstein-von Mises Gaussian approximation is very good for both posteriors.

3 Divide-and-conquer approaches

Divide-and-conquer methods run MCMC on data batches and combine the resulting posterior approximations. Their main challenges are support mismatch, poor scaling with the number of batches, limited theoretical guarantees, and difficult-to-interpret combinations.

Divide-and-conquer methods split the dataset into batches, run separate MCMC chains, and combine the batch posterior results.
Batch-posterior combinations include Gaussian or importance-sampling approximations, averaged samples, kernel-density products, and expectation-propagation-like procedures.
These methods are theoretically justified mainly when batch posteriors are Gaussian or batch sizes tend to infinity.
For kernel-density products, the available mean-squared-error bound is exponential in the number of batches and therefore uninformative for tall data.
Nearly disjoint batch-posterior supports can make multiplicative approximations poor, motivating Weierstrass transforms and related extended-model constructions.
Overall, efficient batch combination remains difficult because posterior multiplication scales poorly, supports may be disjoint, guarantees are often asymptotic, and cheap combinations are hard to interpret.

4 Exact subsampling approaches: Pseudo-marginal MH

Pseudo-marginal MH uses nonnegative unbiased likelihood estimators to retain the target posterior, but estimator variance and tall-data scaling can make proposed constructions impractical. Auxiliary-variable Firefly MH can reduce evaluations in favourable cases, yet its performance may deteriorate sharply for difficult targets.

Pseudo-marginal MH: Pseudo-marginal MH targets the correct posterior by substituting an almost-surely nonnegative unbiased estimator of the unnormalized target.The estimator is reused for the current state and newly computed for proposals.
Pseudo-marginal MH: Estimator variance increases asymptotic variance and can degrade the underlying chain’s convergence properties.Large estimator variance can produce misleadingly high likelihood estimates and poor mixing.
Pseudo-marginal MH: A log-likelihood-estimator variance around 1.5 is recommended when ideal-chain autocorrelation times are unknown; values below 0.2 or above 10 incur large penalties.Around 1.0 is preferred for quasi-i.i.d. ideal-MH samples and around 3.0 for very large integrated autocorrelation times.
Unbiased likelihood estimation: The Rhee–Glynn construction yields a nonnegative unbiased likelihood estimator, but controlling its relative variance requires batchsize O(n^2) and expected truncation length O(n).These requirements defeat the intended reduction in data usage.
Unbiased likelihood estimation: The resulting pseudo-marginal approach produced poorly mixing chains even on the Gaussian example and was judged unsuitable for tall-data MH.The authors tested multiple truncation parameters and batch sizes without satisfactory results.
Auxiliary variables: Firefly MH can use a subset of data under bounds ℓ_i(θ) ≥ b_i(θ), with likelihood evaluations controlled through attempted flips or the fraction flipped to 1.The expected equilibrium cost of a conditional parameter update remains O(n), with the constant tied to bound tightness.
Auxiliary variables: A tight Taylor bound yielded likelihood evaluations close to 10% in both examples, but bright points caused erratic lognormal behaviour and very slow mixing.The first 1,000 iterations were discarded as burn-in because poor initial auxiliary variables could move the chain out of the posterior mode.

5 Other exact approaches

Other exact approaches include SGLD and delayed-acceptance schemes that reduce per-iteration computation while preserving or approximating target-invariance under stated conditions. Their benefits are bounded by convergence, bias, noise, or full-data evaluation costs.

Stochastic approximation: SGLD replaces exact gradients with unbiased subsample-based score estimates and removes the MH accept–reject step while decreasing stepsizes.Under suitable stepsize conditions, the resulting approximation provides a consistent estimate of the target.
Stochastic approximation: SGLD’s theory leaves unclear whether using fewer data points per iteration offsets the additional iterations needed to reach the same accuracy.Its convergence rate is slower than the traditional Monte Carlo rate N^-1/2.
Stochastic approximation: In the lognormal example, SGLD remained far from convergence because subsampling moved the chain away from posterior support.Smaller stepsizes reduce overconfident moves but increase the variance of the final estimate.
Subsampling HMC: Subsampling HMC modifications include heuristic methods that rely on Gaussian subsampling noise, an assumption that can produce arbitrarily poor performance when violated.The cited modification is inspired by decreasing-stepsize SGLD.
Delayed acceptance: Delayed-acceptance MH can reject proposals using batches and parallel calculations without adding bias, but accepted points still require full-data evaluation.Its average data fraction is therefore bounded below by the acceptance rate, often around 25%.
Delayed acceptance: Delayed acceptance can have larger asymptotic variance, need not inherit the original chain’s ergodicity, and may reduce convergence speed.The paper characterizes it as easy to implement but inherently limited in benefit.

6 Approximate subsampling approaches

Approximate subsampling methods trade exact posterior targeting for reduced data use, but their accuracy and computational gains depend strongly on controlling subsampling noise and modeling assumptions. The reviewed approaches include naive subsampling, Gaussian or CLT-based corrections, adaptive tests, and confidence samplers with theoretical guarantees.

6.1 Naive subsampling: Naive subsampling targets an approximation to π rather than sampling exactly from the posterior, and its target can be difficult to interpret.Using exponentiated unbiased average log-likelihood estimates produces a nontrivial target rather than the desired posterior.
6.1 Naive subsampling: Naive subsampling broadens the target because likelihood contributions from small subsamples remain influential, with additional broadening when bias correction is absent.The largest substantially contributing subsample size is roughly nλ, while small subsample likelihoods broaden the resulting distribution.
6.2 Relying on the CLT: Gaussian-noise methods can be exact when the log-likelihood-ratio noise has known variance and the correct mean, but estimated-variance variants are inexact and rely crucially on Gaussianity.Heavy-tailed log-likelihood ratios can violate this assumption, causing tail contributions to be missed and acceptance decisions to become erroneous.
6.2.2 Adaptive subsampling with T-tests: Austerity MH uses sequential T-tests to decide acceptance with subsamples, but its performance can depart substantially from the target when the CLT-based Student approximation is poor.In the Gaussian example the posterior is centered but slightly too wide; in the lognormal example the chain departs more markedly from the actual posterior.
6.3 Exchanging acceptance noise for subsampling noise: Methods avoiding CLT-based approximations can provide weaker-assumption guarantees, yet the confidence sampler still requires O(n) samples per iteration at equilibrium.The proportionality constant is the variance of the log-likelihood ratio under subsampling.
6.3 Exchanging acceptance noise for subsampling noise: Robustness to non-Gaussianity can force O(n) subsamples per iteration even when the first three moments of subsampling noise are controlled.For logistic regression, proposal scaling implies that maintaining temperature β of order 1 requires t of order n.
6.4 Confidence samplers: A confidence sampler with a second-order Taylor proxy can leverage local target concentration to reduce likelihood evaluations while retaining the confidence sampler’s theoretical guarantee.The proxy acts as a control variate in a concentration inequality and replaces many likelihood evaluations with cheap Taylor-expansion evaluations.

7 An improved confidence sampler

The paper improves the confidence sampler by using likelihood proxies as control variates, preserving controlled approximation guarantees while potentially reducing subsampling cost. Taylor proxies can yield o(n) or, under especially favourable conditions, near-constant likelihood evaluations per iteration when the posterior is approximately Gaussian.

Proxy construction: Likelihood proxies act as control variates for individual likelihoods, reducing the variance of log likelihood ratios used by the confidence sampler.The method requires proxies whose approximation error or range can be bounded cheaply.
Guarantees: The proxy sampler retains the confidence sampler’s controlled approximation to the posterior and convergence guarantees under the underlying sampler’s ergodicity assumptions.Proposition 7.1 states uniform ergodicity inheritance, an O(δ) target discrepancy, and convergence-speed control by δ.
Proxy construction: Good proxies produce more accurate confidence bounds, which reduces the number of likelihood samples required for acceptance decisions.The improvement depends on the proxy being accurate for the proposed state transition.
Taylor proxies: Taylor expansions provide usable proxies for models with up to third derivatives, although finding a good problem-specific proxy remains a strong assumption.A second-order Taylor proxy is defined using gradients and Hessians at a reference point, with Taylor-Lagrange bounds controlling the residual.
Taylor proxies: When posterior mass concentrates near the maximum likelihood estimate, a single Taylor proxy can represent the target accurately; poor concentration or nonlocal proposals can reduce gains.Multiple proxies or regularly refreshed proxies are suggested when one reference point is inadequate.
Subsampling gain: Under an approximately Gaussian target and suitable derivative-growth conditions, proxy-based confidence sampling can require o(n) likelihood evaluations per iteration, with Gaussian data achieving O(1) up to polylogarithmic factors.The original confidence sampler instead requires O(n) samples per iteration at equilibrium.

8 Experiments

Experiments evaluate Taylor-proxy confidence samplers on synthetic and covtype logistic and gamma regression problems. The methods use substantially fewer likelihood evaluations than full MH in the reported settings, with gains associated with periodically recomputed proxies.

Experimental setup: The experiments use disk-based SQLite storage to avoid loading the dataset and proxy quantities into memory.This implementation serves as a proof of concept for handling large datasets.
8.1 Logistic regression: On synthetic logistic regression, confidence MH is tested across datasets with log10 n from 3 to 7 using a single proxy centered at the MAP.The experiment runs 10 000 iterations with δ = 0.1 and reports likelihood evaluations per iteration versus n.
8.1 Logistic regression: 27 to 42% of n likelihood evaluations per iteration are used on covtype logistic regression, with 50% of iterations requiring less than 5% of the dataset.Proxy recomputation every 10 iterations accounts for 20% of n, and the gains are two- or threefold in this experiment.
8.2 Gamma regression: Taylor proxies can be applied to gamma regression when the response model’s derivatives and covariate norms satisfy the required bounds.The experiment assumes known κ and uses the same covtype data with a nonnegative response feature.
8.2 Gamma regression: 33 to 54% of n likelihood evaluations per iteration are used on covtype gamma regression, with 50% of iterations requiring less than 10% of the dataset.Proxy recomputation every 10 iterations accounts for 20% of n, and subsampling converges two or three times faster in this example.

9 Discussion

The discussion concludes that the proposed subsampling approach can approximate the posterior with lower computational cost under strong ergodicity assumptions. Its largest observed gains occur when the Bernstein-von Mises approximation is excellent, leaving poorer-approximation regimes as an open challenge.

Discussion: Divide-and-conquer methods still face the problem of stably recombining distributions from chains run on many smaller datasets.Subsampling methods instead must control target approximation and likelihood-evaluation budgets.
Discussion: The proposed subsampling approach has a controlled posterior approximation under strong ergodicity assumptions, while experiments suggest applicability beyond those assumptions.The authors note that these assumptions are rarely satisfied in practice.
Scaling: The methodology can reduce the natural O(n) subsample cost per iteration to as low as O(1) in favourable scenarios.The observed improvement over other subsampling approaches is reported specifically where the Bernstein-von Mises approximation is excellent.
Limitations: The methods have not yet demonstrated comparable gains when Bernstein-von Mises is difficult to compute or provides a poor approximation.Further work is directed toward settings where n is not large enough for the approximation to be accurate.

Appendix A: proof of Proposition 4.1

The appendix derives the second moment of the nonnegative unbiased likelihood estimator used in the pseudo-marginal analysis. The derivation uses the conventions for the estimator’s truncation variables and the definition of N.

Proof: The appendix obtains the second moment of Y from the Rhee–Glynn result.This supports analysis of the estimator’s variance-related behaviour.
Proof: The derivation applies the conventions S−1 = 0 and S0 = ena(θ) before using the definition of N.These conventions determine the boundary terms in the calculation.