PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficiently Sampling Functions from Gaussian Process Posteriors

James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, Marc Peter Deisenroth

arXiv:2002.09309v4stat.MLcs.LGstat.CO

TL;DR

GP posterior quantities often require Monte Carlo, but accurate draws remain expensive at many test locations. The paper separates the prior from the data and combines bases suited to each component. Decoupled sample paths retain posterior fidelity while enabling linear-time evaluation and compatibility with sparse approximations.

  • Problem

    Accurate GP posterior sampling is costly at test time: exact and sparse methods scale cubically in test locations, while Fourier features can misrepresent posteriors.

  • Method

    The method decomposes the posterior into a prior and a data-driven update, representing them separately with sampling-suited bases and combining this with approximate GPs.

  • Results

    Decoupled sample paths evaluate in O(∗) time while retaining statistical properties across dense and sparse data regimes.

  • Takeaways & Limitations

    Decoupled sampling provides a general-purpose plugin for sample-based algorithms using sparse GPs, combining scalable evaluation with posterior-faithful draws.

  • Takeaways & Limitations

    Using an approximate prior introduces additional test-time sampling error, which is often smaller than the error from inducing-point approximations.

Abstract

from arXiv · show

Gaussian processes are the gold standard for many real-world modeling problems, especially in cases where a model's success hinges upon its ability to faithfully represent predictive uncertainty. These problems typically exist as parts of larger frameworks, wherein quantities of interest are ultimately defined by integrating over posterior distributions. These quantities are frequently intractable, motivating the use of Monte Carlo methods. Despite substantial progress in scaling up Gaussian processes to large training sets, methods for accurately generating draws from their posterior distributions still scale cubically in the number of test locations. We identify a decomposition of Gaussian processes that naturally lends itself to scalable sampling by separating out the prior from the data. Building off of this factorization, we propose an easy-to-use and general-purpose approach for fast posterior sampling, which seamlessly pairs with sparse approximations to afford scalability both during training and at test time. In a series of experiments designed to test competing sampling schemes' statistical properties and practical ramifications, we demonstrate how decoupled sample paths accurately represent Gaussian process posteriors at a fraction of the usual cost.

1. Introduction

The paper targets accurate, scalable sampling from GP posteriors for Monte Carlo estimation when posterior quantities are analytically intractable. It proposes separating the prior from the data so different bases can address sampling efficiency and posterior fidelity.

  • Predictive uncertainty from GPs supports decision-making tradeoffs, while nonlinear posterior functionals often require Monte Carlo estimation.
  • Existing GP posterior samplers trade statistical behavior against computational cost: naive methods scale poorly, while Fourier-feature approximations can misrepresent posteriors.
  • The proposed approach decomposes the GP posterior into a prior and an update, representing each term with a sampling-suited basis.
  • Decoupled sampling combines approximate-GP techniques with this decomposition to sample GP posteriors accurately in linear time.

2. Review of Gaussian processes

The review contrasts exact, sparse, and weight-space GP sampling. Exact and sparse approaches preserve posterior structure but can scale poorly at test time, whereas Fourier-feature methods sample efficiently yet may deteriorate as data increases.

  • Gaussian-process setup: A GP defines jointly Gaussian function values at any finite set of locations, with covariance determined by a kernel; the posterior is Gaussian after conditioning on observations.
  • Exact posterior sampling: Exact sampling applies a location-scale transform using a matrix square root, providing statistically exact draws up to numerical error.
  • Exact posterior sampling: O(∗^3) cost for computing the posterior covariance square root makes standard sampling poorly scalable in the number of test locations.
  • Function-space approximations: Sparse GPs represent functions through inducing values and can train in O(˜nm^2) time instead of exact O(n^3) scaling when m ≪ n.
  • Function-space approximations: Sparse-GP posterior moments are cheaper to compute, but standard sampling still costs O(∗^3), making deployment slow despite efficient training.
  • Weight-space approximations: Random Fourier features sample weight vectors defining functions at arbitrary locations, with linear-algebra cost O(min{ℓ,n}^3) enabled by the Woodbury identity.
  • Weight-space approximations: Finite-dimensional feature maps lose expressivity; Fourier-feature approximations can suffer variance starvation and increasingly deteriorate away from data as n increases.
  • Motivation: Function-space sampling scales poorly with test locations, whereas weight-space sampling scales gracefully there but incurs increasing approximation error with more training data.

3. Sampling with Matheron’s rule

Matheron’s rule decomposes GP posterior draws into a prior sample plus a corrective update, enabling separate bases for efficient sampling. Decoupled sampling combines an approximate weight-space prior with a function-space update to obtain scalable posterior paths while preserving statistical properties.

  • Matheron’s rule: Matheron’s rule represents a conditional Gaussian draw as a prior draw corrected using residuals from the observed conditioning value.The GP analogue updates realized prior sample paths rather than forming posterior covariance matrices or their square roots.
  • Pathwise updates: The pathwise GP construction draws a prior function, forms an update from residuals at inducing locations, and adds the two functions to obtain a posterior draw.The update uses an independently sampled inducing target and can be evaluated throughout the input space.
  • Decoupled bases: The decomposition diagnoses function-space sampling’s O(|X∗|^3) prior cost and weight-space sampling’s data-representation problem, motivating different bases for the two terms.The approach separates prior and data effects so each component can be represented using a basis suited to its computational or statistical weakness.
  • Decoupled bases: Decoupled sampling uses a weight-space approximate prior and canonical basis updates tied one-to-one to inducing locations, preserving data representation with fewer total basis functions.The paper reports that its statistical properties mirror the gold-standard sampler despite using half as many basis functions as the weight-space method.
  • Practical properties: Decoupled sampling scales linearly in the number of test locations and produces pathwise differentiable functions, supporting efficient evaluation and gradient-based analysis of GP extrema.These properties complement the strengths of function draws while avoiding the cubic test-location cost of standard sampling.
  • Error analysis: The method’s error separates into sparse-GP and approximate-prior terms, while its mean matches the sparse-GP posterior and its RFF covariance error decays dimension-free as ℓ increases.The prior-to-posterior error transfer is controlled by inducing locations Z, not by the inducing distribution q(u).

4. Experiments

The experiments evaluate decoupled sampling across uncertainty calibration, Thompson sampling, and dynamical-system simulation. Decoupled paths retain posterior-relevant statistical properties while enabling linear-time evaluation and simulation.

  • Decoupled sampling is evaluated for uncertainty calibration, Thompson sampling, and dynamical-system simulation.
  • Uncertainty calibration: Weight-space sampling deteriorates as inducing locations increase relative to the basis budget because variance starvation misrepresents extrapolatory posterior behavior.
  • Uncertainty calibration: Decoupled sampling retains its performance, and may even improve, as the number of inducing locations increases relative to the basis functions.
  • Thompson Sampling: Decoupled sample paths have linear evaluation complexity O(∗), support pathwise differentiation, and retain statistical properties when b is comparable to n or b ≪n.
  • Dynamical-system simulation: For dynamical-system simulation, decoupled sampling avoids iterative conditioning and enables trajectories to be simulated in linear time O(t).
  • Dynamical-system simulation: Neuron simulations requiring 10 hours with iterative sampling ran in 20 seconds with decoupled sampling while achieving similar accuracy.

5. Conclusion

The conclusion presents decoupled sampling as a general, extensible approximation strategy that combines separate representations of the prior and data. Its speed–accuracy balance remains task-dependent.

  • The best treatment of prior and update components depends on the task, especially when structured covariance matrices permit exact prior sampling.
  • Decoupled sampling separates the prior and data into bases better suited for sampling, combining strengths of previous approaches.
  • The method can serve as a plugin extension to existing sample-based algorithms driven by sparse GPs.
  • Decoupled sample paths can be evaluated in linear time without fear of misrepresenting GP posteriors.

A. Alternative decompositions

The paper considers direct weight-space and orthogonal decompositions as alternatives to representing GP posteriors as a weight-space prior plus a function-space update. Each alternative has sampling limitations.

  • A direct weight-space sparse posterior uses a Bayesian linear model f(·) = φ(·)⊤w with inducing information q(u).
  • The direct weight-space representation forces f to pass exactly through u because it lacks measurement noise σ2.
  • An orthogonal decomposition writes f(·) = f∥(·) + f⊥(·), with statistically independent components satisfying Cov(f∥, f⊥) = 0.
  • In the orthogonal decomposition, the residual f⊥ is independent of inducing values and vanishes at inducing locations Z.
  • Sampling the orthogonal residual is difficult, and random Fourier feature variants inherit the exact-interpolation issues of the direct representation.

B. Error analysis

The error analysis defines the GP approximations and derives Wasserstein-distance bounds using pathwise inequalities, operator norms, and continuity assumptions.

  • Setup: The analysis compares the true posterior with weight-space, inducing-point, and decoupled posterior approximations through their respective kernels.The setup assumes a GP on a compact domain with Gaussian likelihood and identifies each approximation explicitly.
  • Bounds: The resulting bounds relate function-space discrepancies to 2-Wasserstein distances over L2(X) and continuous functions equipped with the supremum norm.The relevant Wasserstein metrics and operator norm are defined over the corresponding function spaces.
  • Proof strategy: The proofs apply the triangle inequality, Hölder’s inequality, operator-norm bounds, and Matheron’s rule to control approximation errors.The argument lifts pathwise bounds to Wasserstein distances by integrating over couplings.
  • Bounds: Continuity of sample paths and separability of C(X) ensure that the Wasserstein distance on continuous functions is a proper metric.The proof uses this property to complete the bound.
  • Bounds: The proof concludes by combining the intermediate inequalities with a constant determined by kernel regularity and spectral properties.The constant depends on the kernel’s Lipschitz behavior, spectral-density decay, and dimension.

C. Additional experiments

Additional experiments evaluate posterior-sampling quality with Wasserstein distances and compare decoupled sampling against LOVE under structured experimental conditions.

  • Experimental setup: Experiments use zero-mean GP priors with Matérn-5/2 kernels, while dynamical-system hyperparameters are learned and other hyperparameters are specified.The specified settings include lengthscales l = d/100, noise variance σ2 = 10^-3, and kernel amplitude α = 1.
  • C.1. 2-Wasserstein sample tests: The 2-Wasserstein sample tests estimate true posterior moments from 100,000 draws generated across randomly sampled training and test locations.Distances are computed between estimated and true posterior moments.
  • C.1. 2-Wasserstein sample tests: Figure 6 reports medians and interquartile ranges of empirically estimated 2-Wasserstein distances over 32 independent trials.Weight-space and decoupled methods use b = m + ℓ basis functions, with ℓ ∈ {1024, 4096, 16384}.
  • C.1. 2-Wasserstein sample tests: LOVE improves as regularly spaced training grids become denser but still lags behind decoupled sampling for matching m.The comparison uses exact GPs with known hyperparameters and gives LOVE favorable structured-grid conditions.
  • C.1. 2-Wasserstein sample tests: LOVE has O(m + ∗) sample-generation complexity, whereas decoupled sampling has O(m × ∗) complexity.The methods also differ in whether they generate values at test locations or function draws evaluable throughout the domain.

C.2. Thompson sampling

The Thompson-sampling experiments compare function-space, weight-space, and decoupled sampling within parallel optimization procedures, using matched experimental runs and runtime comparisons.

  • Benchmark: The benchmark compares Thompson-sampling variants against Random Search and sequential Dividing Rectangles on minimization tasks drawn from a known GP prior.Global minima are estimated with gradient descent from many starting locations to measure regret.
  • Function-space TS: The method selects the s = 2048 smallest mesh values, jointly samples their posterior values, and chooses batch elements by minimizing those samples.The active set is formed before joint posterior sampling.
  • Procedure: Each Thompson-sampling iteration shares a newly generated mesh across batch elements while running the sampling and optimization steps independently.L-BFGS-B performs the multi-start gradient-based optimization in the function-draw procedure.
  • Evaluation: Figure 7 reports parallel Thompson-sampling results as quartiles over 32 independent runs with matched seeds.Weight-space and decoupled variants use a similar procedure with differences arising from their function draws.

C.3. Dynamical systems

The dynamical-systems experiment studies sequential Monte Carlo runtimes using sparse GPs to simulate the FitzHugh–Nagumo model and compares iterative posterior unrolling costs.

  • Experimental setup: The experiment uses the FitzHugh–Nagumo neuron model with diffusion coefficient Σ = 0.01 · I and step size Δt = 0.25.Training and simulation both use the same step size.
  • Training: Independent sparse GPs with m = 32 shared inducing locations model 3-dimensional state-control inputs during training.Training uses 10,000 uniformly random points and stochastic gradient descent.
  • Baseline unrolling: The baseline iteratively draws drift vectors, adds each current input to the inducing set, and augments the inducing distribution with the sampled drift.The covariance update uses the preceding inducing locations and maintains prior covariance factors online.
  • Runtime analysis: The Cholesky-factor update requires O(m^2) time per iteration because only the leading principal submatrix is modified.The remaining terms are zero because the sampled drift is directly observed.
  • Runtime analysis: Because the matrices become dense, the overall iterative unrolling procedure scales cubically in the number of steps.Posterior-marginal computation and online covariance maintenance incur O((m + t)^2) cost per step.
Loading 2002.09309v4…