PaperlayerView exact revision on arXiv

Source-linked AI summary

Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming

Gabriel Riutort-Mayol, Paul-Christian Bürkner, Michael R. Andersen, Arno Solin, Aki Vehtari

arXiv:2004.11408v2stat.COstat.ME

TL;DR

Exact Gaussian processes become computationally costly for large datasets, particularly with full Bayesian sampling. This paper analyzes a Laplace-eigenfunction basis approximation and develops practical recommendations and diagnostics for choosing its basis count and boundary. The authors report improved sampling efficiency over exact GPs on simulated and real datasets, while noting that computation grows exponentially with input dimensionality.

  • Problem

    Exact GP computation is costly for large datasets, especially when full Bayesian inference repeatedly inverts the covariance matrix during sampling.

    Source: §1 · p. 1

  • Method

    The paper analyzes a low-rank approximation of stationary GPs using Laplace eigenfunctions and develops diagnostics and an iterative procedure for selecting basis count and prediction-domain boundary.

    Source: §1 · p. 2

  • Results

    The authors demonstrate practical applicability and improved sampling efficiency over exact GPs on simulated and real datasets.

    Source: §6 · p. 21

  • Takeaways & Limitations

    Recommendations and diagnostics help assess basis-count and domain-boundary choices for accurate approximations with minimal computational cost.

    Source: §1 · p. 2, §6 · p. 21

    The multivariate basis count grows exponentially with input dimensions, making the method more costly as dimensionality increases.

    Source: §3.2 · p. 6

Abstract

from arXiv · show

Gaussian processes are powerful non-parametric probabilistic models for stochastic functions. However, the direct implementation entails a complexity that is computationally intractable when the number of observations is large, especially when estimated with fully Bayesian methods such as Markov chain Monte Carlo. In this paper, we focus on a low-rank approximate Bayesian Gaussian processes, based on a basis function approximation via Laplace eigenfunctions for stationary covariance functions. The main contribution of this paper is a detailed analysis of the performance, and practical recommendations for how to select the number of basis functions and the boundary factor. Intuitive visualizations and recommendations, make it easier for users to improve approximation accuracy and computational performance. We also propose diagnostics for checking that the number of basis functions and the boundary factor are adequate given the data. The approach is simple and exhibits an attractive computational complexity due to its linear structure, and it is easy to implement in probabilistic programming frameworks. Several illustrative examples of the performance and applicability of the method in the probabilistic programming language Stan are presented together with the underlying Stan model code.

1 Introduction

Exact Gaussian-process inference becomes costly for large datasets, especially with full Bayesian sampling. This paper analyzes a Laplace-eigenfunction basis approximation and provides practical guidance for its use in probabilistic programming.

1 Introduction

  • Exact GP inference requires O(n^3) computation and O(n^2) memory, with each full-Bayesian sampling step requiring O(n^3) matrix inversion.These costs limit applications to datasets of at most a few tens of thousands of observations.Source: §1 · p. 1
  • The introduction situates the method among sparse and basis-function GP approximations, including spectral methods and random Fourier features.It also notes precision-matrix approaches that exploit sparsity to improve computational scaling.Source: §1 · pp. 1–2
  • The paper analyzes a stationary-GP approximation using Laplace eigenfunctions, which represents GPs as linear models and allows accuracy–complexity trade-offs.The eigenfunctions are computed analytically and do not depend on the covariance function or its hyperparameters.Source: §1 · p. 2
  • The basis-function approach reduces repeated covariance-parameter learning to O(mn + m) per optimizer or sampler step after O(m^2n) precomputation.Its parameterization also supports non-centered inference, reducing posterior dependency between function and covariance parameters.Source: §1 · p. 2
  • Practical recommendations and diagnostics assess the number of basis functions and prediction-domain boundary, with an iterative procedure to seek accurate approximations at minimal computational cost.The analysis relates these choices to the domain and properties of the function being learned.Source: §1 · p. 2
  • The authors implement the method in Stan and brms and demonstrate it with simulated and real datasets, providing corresponding Stan code.Their focus is modular use in probabilistic programming rather than the fastest inference for a restricted set of GP models.Source: §1 · p. 2

2 Gaussian process as a prior

The section establishes Gaussian processes as priors over functions and introduces covariance and spectral representations that underpin the later Hilbert-space approximation.

2 Gaussian process as a prior

  • A GP prior defines a distribution over functions, with finite sets of function values jointly Gaussian; the paper uses this prior for latent functions in Bayesian models.The goal is inference for function values at new inputs, including models with non-Gaussian observations or multiple unknown functions.Source: §2 · p. 3

2.1 Covariance functions and spectral density

  • Covariance functions encode prior assumptions and correlations; stationary covariances depend on input differences, and isotropic ones depend only on their norm.Stationary covariance functions have spectral-density representations, which are Fourier duals of the covariance functions.Source: §2.1 · pp. 3–4

3.1 Unidimensional GPs

  • The Hilbert-space construction treats a stationary covariance operator on a compact domain with boundary conditions, with the domain boundary affecting approximation accuracy.For univariate inputs, the domain is Ω = [−L, L]; the section introduces covariance expansions and spectral densities for common kernel cases.Source: §2.1 · p. 4, §3.1 · p. 4

3 Hilbert space approximate Gaussian process model

The model approximates stationary GPs with a truncated Laplace-eigenfunction expansion, yielding a lower-dimensional linear representation. Its multidimensional extension multiplies basis counts across dimensions, while periodic covariance functions require a separate low-rank representation.

3 Hilbert space approximate Gaussian process model

  • The univariate HSGP truncates a Laplace-eigenfunction expansion, scaling basis functions by square roots of spectral-density values to approximate the function.The eigenfunctions are independent of covariance hyperparameters, and log-posterior evaluation costs O(nm + m).Source: §3 · pp. 4–5
  • The expansion gives a non-centered model with independent standard-normal weights, reducing inference to an m-dimensional parameter space where m is much smaller than n.The covariance and its hyperparameters enter through the prior distribution of the regression weights.Source: §3 · p. 5

3.2 Generalization to multidimensional GPs

  • In D dimensions, the approximation combines univariate eigenfunctions across dimensions, so the total basis count is the product of the per-dimension counts and grows exponentially with D.Log-posterior evaluation scales as O(nm* + m*), where m* is the number of multivariate basis functions.Source: §3.2 · pp. 5–6

3.3 Linear representation of a periodic squared exponential covariance function

  • Periodic covariance functions do not fit the HSGP framework described here because they lack a spectral representation, but they have a separate low-rank representation.The paper presents and analyzes an approximate linear representation for the periodic squared-exponential covariance function.Source: §3.3 · p. 6

4 The accuracy of the approximation

The section analyzes how basis count, boundary choice, and function length-scale affect HSGP approximation accuracy, and develops recommendations and diagnostics for choosing settings. It also examines these relationships for selected covariance kernels.

4 The accuracy of the approximation

  • Approximation accuracy and speed depend on interacting choices of basis count and boundary condition, so the section develops selection recommendations and diagnostics for the obtained approximation.The recommendations relate m, boundary L, and length-scale ℓ; their relationships vary with kernel properties.Source: §4 · p. 6

4.1 Dependency on the number of basis functions and the boundary condition

  • With the boundary factor sufficiently large, too few basis functions make the estimated function overly smooth because high-frequency components are missing.The study illustrates this effect using data generated from a squared-exponential GP with length-scale 0.3 and marginal variance 1.Source: §4.1 · p. 6, Figure 1 · p. 7
  • With basis count held sufficiently large, changing c affects the predictive mean mainly near boundaries but affects standard deviation across the domain, with larger error for standard deviation.Figure 2 compares the exact GP and HSGP means, standard deviations, and covariance functions across boundary-factor values.Source: §4.1 · p. 7, Figure 2 · p. 7
  • For this example, dynamic HMC jointly estimated function values and hyperparameters; the optimal precision–computation-time choice was m = 15 basis functions with boundary factor c between 1.5 and 2.5.The comparison fitted both exact GP and HSGP models.Source: §4.1 · p. 7
  • As c increases, m must also increase, and below a minimum c an accurate approximation cannot be achieved by increasing m alone.The relationships depend on the kernel; the section studies squared-exponential and Matérn kernels with ν = 3/2 and ν = 5/2.Source: §4 · p. 6, §4.1 · p. 7

4.2 Near linear proportionality between m, c and ℓ

The required basis count depends on length-scale and boundary factor: for squared exponential and Matérn kernels, their relationship is near-linear above a recommended boundary-factor minimum. Empirical results and practical equations guide parameter selection, with multidimensional use remaining an approximation.

4.2 Near linear proportionality between m, c and ℓ

  • Cumulative spectral densities can estimate how many basis functions explain nearly all variation for given c and ℓ.The approach applies to any covariance function for which the spectral density is available.Source: §4.2 · p. 7
  • For squared exponential and Matérn kernels above the recommended minimum c, m grows near-linearly with c and 1/ℓ.Shorter length-scales require more high-frequency basis functions; larger c uses a smaller range of basis functions, also requiring more high-frequency terms.Source: §4.2 · p. 9

4.3 Empirical discovering of the functional form of the relationships between m, c and ℓ

  • As ℓ increases, fewer basis functions m are needed; below a minimum c for a given ℓ, accurate approximation cannot be achieved regardless of m.The empirical accuracy criterion is a relative total-variation difference below 1% of the exact covariance’s total area.Source: §4.3 · p. 10

4.3.1 Numerical equations

  • The practical equations require c ≥ 1.2; the Matérn ν=3/2 formula is a conservative choice for Matérn kernels with ν ≥ 3/2.The one-dimensional equations can be applied dimension by dimension as an approximation, but multidimensional accuracy may depend on ratios among the dimensions’ length-scales.Source: §4.3.1 · pp. 10–11

4.4 Relationships between m and ℓfor a periodic squared exponential covariance function

For periodic squared exponential kernels, the approximation uses a sinusoidal series whose basis-count relationship with length-scale does not depend on a boundary condition. The paper also gives diagnostics and an iterative procedure for selecting basis count and boundary factor, while warning that the length-scale threshold alone does not guarantee adequate accuracy.

4.4 Relationships between m and ℓfor a periodic squared exponential covariance function

  • The periodic squared exponential approximation uses a sinusoidal series, so its required basis count depends on length-scale but not on a boundary condition.The paper derives a minimum number of terms for a close approximation to the exact periodic kernel.Source: §4.4 · p. 11

4.5 Diagnostics of the approximation

  • The equations give minimum m and c values for computational efficiency and a minimum inferable length-scale for diagnosing approximation accuracy.If the estimated length-scale falls below the threshold for the chosen m and c, increase m and/or c; exceeding it does not guarantee accuracy, so also check predictive measures.Source: §4.5 · p. 11

4.5.1 A step-by-step user-guide to apply the diagnostics

  • The recommended procedure iteratively updates m and c using the estimated length-scale, then checks predictive accuracy and stability before stopping.It begins from a length-scale guess, fits the HSGP, and repeats adjustments until the length-scale threshold and predictive checks are satisfactory and stable.Source: §4.5.1 · pp. 11–12, §4.5 · p. 11

4.6 Performance analysis of the diagnostics

The diagnostics reached optimal m and c values with minimum computational requirements in the examined datasets. Accurate HSGP length-scale estimates were associated with small RMSE against the exact GP.

4.6 Performance analysis of the diagnostics

  • Accurate HSGP length-scale estimates implied small RMSE against the exact GP across noisy datasets with varying smoothness.The comparison used squared exponential kernels and set c to a valid, optimal value for each m.Source: §4.6 · p. 12
  • The iterative procedure was illustrated on some datasets also used to compare exact-GP and HSGP length-scale estimates.Concrete examples of two datasets are described in Appendix C.Source: §4.6 · p. 12
  • Two to four iterations were sufficient to reach optimal m and c, depending on function wiggliness and the distance between the initial and true length-scales.The reported cases achieved these values with minimum computational requirements by following the recommendations.Source: §4.6 · p. 12

4.7 Other covariance functions

For stationary covariance functions with an available spectral density, the basis-function approximation can be applied beyond the Matérn and squared exponential families. When family-specific diagnostic plots or equations are unavailable, the paper proposes covariance-distance and spectral-density checks.

4.7 Other covariance functions

  • The basis-function approximation can be implemented for any stationary covariance function with an available spectral density.Source: §4.7 · p. 12
  • Without family-specific diagnostic plots or equations, relative total variational distance can assess approximation quality for a relevant length-scale.A small bound on this distance implies a high-quality approximation.Source: §4.7 · p. 12
  • Cumulative spectral-density sums offer another way to assess accuracy by measuring how much total variation the included series terms explain.Source: §4.7 · p. 12

4.8 The computational cost in the multi-dimensional setting

HSGP can outperform exact GPs in low-dimensional settings, but its cost rises with dimension because the multivariate basis count grows rapidly. The speed advantage depends on dataset size, function smoothness, and implementation.

4.8 The computational cost in the multi-dimensional setting

  • In 1D and 2D, HSGP was computationally superior even for highly wiggly functions, except when n ≲300, where exact GPs were already reasonably fast.The exact GP’s cubic scaling in n can still make HSGP faster for larger datasets.Source: §4.8 · p. 12
  • For most nonlinear 2D problems and moderate-to-large 3D datasets (n ≳1000), HSGP was faster than exact GP, even for highly wiggly 3D functions.For small datasets, HSGP may instead be slower for highly to moderately wiggly 3D functions and smooth 4D functions.Source: §4.8 · p. 12
  • With an optimal boundary factor, the diagnostic can reduce the multivariate basis count to a sufficient minimum, supporting highly wiggly 3D and smooth 4D applications.For a squared exponential covariance, the recommended per-dimension basis count increases as length-scales decrease.Source: §4.8 · p. 14
  • Which method is faster depends on implementation details that can substantially change scaling constants, so finer speed recommendations are software-specific.Source: §4.8 · p. 14
  • Thin plate splines had roughly similar or worse scaling than HSGP, with serious computation-time challenges above 10 basis functions per dimension in 3D or 40 in 2D.These thresholds describe the reported spline inference experiments.Source: §4.8 · p. 14

5 Case studies

The case studies apply HSGP fitting and diagnostics to simulated and real datasets. They examine predictive accuracy, computation, and the use of HSGP components in structured models.

5.1 Simulated data for a univariate function

  • For this univariate dataset, HSGP was roughly 400 times faster than exact GP and 10 times faster than the spline model, while approximating GP predictions in interpolation and extrapolation.The spline model did not extrapolate properly; HSGP computation time increased slowly with basis count.Source: §5.1 · p. 15, Figure 9 · p. 15, Figure 10 · p. 15

5.2 Birthday data

  • For the birthday data, the model decomposes birth counts into long-term, yearly, weekly, and special-day effects, with figures displaying the trend and seasonal components.Periodic covariance functions model the year and week effects, while a Student’s t prior models special-day effects.Source: §5.2 · pp. 15–16

5.3 Case study: Simulated data for 2D and 3D functions

  • In the 2D and 3D simulations, four diagnostic iterations reached the minimum basis counts needed for accurate approximation; in 2D, HSGP and exact GP were accurate, while the spline model was significantly worse.The final counts were 242 multivariate basis functions for 2D and 2520 for 3D; the figure compares 2D function errors and RMSE across methods and settings.Source: §5.3 · p. 18, Table 4 · p. 18, Table 5 · p. 18

5.3.1 Computation requirements in 2D, 3D, and 4D input spaces

  • In the reported computation study, HSGP was faster for most nonlinear 2D functions and moderate-to-large 3D datasets, but small datasets can favor exact GP; above D > 5, HSGP becomes impractical even for smooth functions.Multivariate basis counts grow multiplicatively across dimensions; lower-dimensional additive components may still use the method.Source: §5.3.1 · p. 19

5.4 Leukemia data

  • In the leukemia analysis, conditional estimates showed nonlinear patterns and a WBC–TDI interaction; with 3 input variables and m = 15, HSGP was slower than exact GP.A larger boundary factor and few basis functions were needed for a good approximation of the smooth functions, and increasing the boundary factor raised computation time.Source: §5.4 · p. 21

6 Conclusion

The paper analyzes a Laplace-eigenfunction low-rank approximation for Bayesian Gaussian processes and provides diagnostics and an iterative procedure for choosing its key settings. It reports practical use in Stan and improved sampling efficiency over exact GPs, while noting computational and scope limits.

6 Conclusion

  • The approximation represents stationary GPs with Laplace eigenfunction basis functions, yielding a linear-model structure with per-evaluation cost O(nm + m).The design matrix is constructed once at cost O(nm), and the posterior is m-dimensional, usually with m much smaller than n.Source: §6 · p. 21
  • The paper analyzes accuracy against basis-function count, boundary conditions, and function nonlinearity, then provides diagnostics and an iterative procedure for accurate approximations with minimal computational cost.The diagnostics and procedure guide selection of the number of basis functions and boundary factor.Source: §6 · p. 21
  • The method demonstrated practical applicability and improved sampling efficiency compared with exact GPs across simulated and real datasets.Its computational complexity grows exponentially with input dimensionality; inputs above three dimensions can be demanding even for moderately wiggly functions and few basis functions per dimension.Source: §6 · p. 21
  • The factor relationships and diagnostics were studied primarily for univariate inputs, leaving more thorough multivariate investigation for future work.For high-dimensional cases, the method may still be used for low-dimensional additive components without very high-dimensional interactions.Source: §6 · p. 21

A Approximation of the covariance function using Hilbert space methods

The appendix summarizes how stationary covariance functions can be approximated through Laplace-operator eigenfunctions. The spectral density supplies approximate covariance-operator eigenvalues, while the domain and boundary condition determine the basis eigenfunctions.

A Approximation of the covariance function using Hilbert space methods

  • For a stationary covariance function, the covariance operator is translation invariant and has a Fourier representation whose transfer function is the GP spectral density.That spectral density gives approximate eigenvalues for the operator.Source: A Approximation of the covariance… · p. 22
  • In the isotropic case, a regular spectral density can be expressed as a polynomial in squared frequency, linking the covariance operator to powers of the negative Laplace operator.The operator is written as a series of Laplace-operator powers, with a corresponding covariance-kernel representation.Source: A Approximation of the covariance… · pp. 22–23
  • On a compact domain, Laplace eigenfunctions satisfy an eigenvalue problem under a chosen boundary condition and are orthonormal, with real positive eigenvalues.The appendix uses Dirichlet boundary conditions while noting that other boundary conditions could also be used.Source: A Approximation of the covariance… · p. 23
  • The resulting covariance approximation uses the spectral density evaluated at Laplace eigenvalues together with the corresponding eigenfunctions on the selected domain.The appendix identifies S(·) as the spectral density, λj as the jth eigenvalue, and φj(·) as its Laplace eigenfunction.Source: A Approximation of the covariance… · p. 23

B Low-rank Gaussian process with a periodic covariance function

The appendix presents a low-rank cosine-series approximation for the periodic squared-exponential covariance, which is outside the HSGP framework discussed in the study. It assesses how many terms are needed for a close approximation and reports that this count grows linearly with the kernel length-scale.

B Low-rank Gaussian process with a periodic covariance function

  • The periodic squared-exponential covariance is approximated by truncating a Taylor-derived cosine series, which converges to the exact covariance as the number of terms J increases.Solin and Särkkä (2014) derive the cosine-series representation.Source: B Low-rank Gaussian process with a… · p. 24
  • Unlike the HSGP approximation covered in the study, this periodic-covariance approximation does not depend on a boundary condition.Its series expansion uses sinusoidal functions.Source: B Low-rank Gaussian process with a… · pp. 23–24
  • The empirical assessment finds a linear proportionality between the minimum term count J and the length-scale ℓ.The minimum J is assessed using a close-approximation criterion with ε = 0.005, corresponding to 0.5% of the exact covariance's total area.Source: B Low-rank Gaussian process with a… · p. 24
  • The periodic covariance admits a low-rank linear representation with cosine and sine terms whose basis functions do not depend on the covariance length-scale ℓ.The hyperparameter dependence enters through the coefficients.Source: B Low-rank Gaussian process with a… · p. 24

C Description of the iterative steps applied to perform diagnostic on two of the data sets in Figure 1

The examples illustrate an iterative HSGP diagnostic procedure that updates length-scale, boundary factor, and basis count, then checks whether the approximation appears sufficiently accurate. The procedure ends when diagnostic and performance measures stabilize.

C Description of the iterative steps applied to perform diagnostic on two of the data sets in Figure 1

  • The other example starts with ℓ1 = 0.5 and computes the boundary factor and minimum basis count before fitting the HSGP model.The supplied text identifies this as a squared-exponential dataset with true length-scale ℓGP = 0.08.Source: C Description of the iterative steps… · p. 24
  • The procedure stops when RMSE, R2, and ELPD are stable relative to the previous iteration, indicating the approximation is likely sufficiently accurate.The examples show this stopping check after iterations 2 and 4.Source: C Description of the iterative steps… · p. 25
  • One example starts with ℓ1 = 1 for a squared-exponential dataset whose true length-scale is ℓGP = 1.4.The example's diagnostic is checked after fitting the initial HSGP model.Source: C Description of the iterative steps… · p. 25
  • Each iteration fits an HSGP model to estimate the length-scale, then uses that estimate to update the next iteration's settings.The steps compute an optimal boundary factor and a minimum basis count from the current length-scale; some iterations add five basis functions when the diagnostic is true.Source: C Description of the iterative steps… · p. 25
Loading 2004.11408v2…