Hilbert Space Methods for Reduced-Rank Gaussian Process Regression

Arno Solin; Simo Särkkä

Source-linked AI summary

Hilbert Space Methods for Reduced-Rank Gaussian Process Regression

Arno Solin, Simo Särkkä

arXiv:1401.5508v3Jan 21, 2014stat.ML

TL;DR

Large Gaussian process datasets make exact inference computationally expensive because direct methods require cubic computation and quadratic memory. This paper uses Laplace-operator eigenfunctions and spectral-density evaluations to form a reduced-rank approximation, yielding O(nm^2) initial cost and O(m^3) marginal-likelihood evaluation while providing dimension-independent convergence rates under appropriate smoothness.

Problem
Direct Gaussian process regression requires O(n^3) computation and O(n^2) memory, limiting applicability as the number of training samples grows.
Method
The method approximates covariance eigendecompositions with Laplace-operator eigenfunctions on compact domains, assigning eigenvalues through the covariance spectral density.
Results
O(nm^2) initial cost and O(m^3) marginal-likelihood evaluation are achieved with basis functions independent of covariance hyperparameters.
Takeaways & Limitations
The convergence rate can remain independent of input dimensionality when covariance differentiability increases appropriately, and squared exponential covariance converges better than ~1/m.
Takeaways & Limitations
The hyperparameter-learning formulation still requires O(n^3) work per optimizer step for the full model because it inverts the n × n covariance matrix.

Abstract

from arXiv · show

This paper proposes a novel scheme for reduced-rank Gaussian process regression. The method is based on an approximate series expansion of the covariance function in terms of an eigenfunction expansion of the Laplace operator in a compact subset of $\mathbb{R}^d$. On this approximate eigenbasis the eigenvalues of the covariance function can be expressed as simple functions of the spectral density of the Gaussian process, which allows the GP inference to be solved under a computational cost scaling as $\mathcal{O}(nm^2)$ (initial) and $\mathcal{O}(m^3)$ (hyperparameter learning) with $m$ basis functions and $n$ data points. Furthermore, the basis functions are independent of the parameters of the covariance function, which allows for very fast hyperparameter learning. The approach also allows for rigorous error analysis with Hilbert space theory, and we show that the approximation becomes exact when the size of the compact subset and the number of eigenfunctions go to infinity. We also show that the convergence rate of the truncation error is independent of the input dimensionality provided that the differentiability order of the covariance function is increases appropriately, and for the squared exponential covariance function it is always bounded by ${\sim}1/m$ regardless of the input dimensionality. The expansion generalizes to Hilbert spaces with an inner product which is defined as an integral over a specified input density. The method is compared to previously proposed methods theoretically and through empirical tests with simulated and real data.

1 Introduction

The paper addresses the cubic computational cost of Gaussian process regression by proposing a reduced-rank covariance approximation based on Laplace-operator eigenfunctions, with basis functions independent of covariance hyperparameters.

1 Introduction: O(n^3) computation and O(n^2) memory make direct Gaussian process regression difficult when the training set is large.The bottleneck is inversion of the n × n Gram matrix K + σ^2_nI.
1 Introduction: Reduced-rank methods approximate the Gram matrix with rank m < n, typically reducing storage to O(nm) and computation to O(nm^2).Several existing approaches can be interpreted spectrally as modifications of the Nyström method.
1 Introduction: The paper proposes approximate covariance eigendecompositions using Laplace-operator eigenfunction expansions on compact subsets of R^d.The covariance function is interpreted as the kernel of a pseudo-differential operator and approximated with Hilbert-space methods.
1 Introduction: The resulting basis functions are independent of covariance functions and their parameters, while the approximation is shown to converge under well-defined conditions.The paper also provides theoretical and experimental comparisons with existing methods.
1 Introduction: The paper develops the approximation, applies it to GP regression, analyzes convergence, compares methods theoretically, and evaluates them on datasets.These topics are covered across Sections 2–6 before the final discussion.

2 Approximating the Covariance Function

The covariance operator is represented through the spectral density and approximated on a Laplace eigenbasis over a compact domain or a density-weighted Hilbert space. Truncation yields a practical series whose accuracy depends on domain, boundary, smoothness, and basis size.

2 Approximating the Covariance Function: For homogeneous covariance functions, the covariance operator can be represented through the spectral density, and isotropy makes that density depend only on frequency magnitude.The paper assumes stationary and isotropic covariance functions for the main development.
2 Approximating the Covariance Function: The covariance operator is interpreted as a pseudo-differential operator formed as a series of powers of the negative Laplace operator.An analytic spectral transfer function yields coefficients for the Laplace-operator series.
2 Approximating the Covariance Function: On a compact domain with boundary conditions, covariance eigenvalues are obtained by evaluating the spectral density at square roots of Laplacian eigenvalues.The corresponding basis functions are Laplacian eigenfunctions, producing simple closed-form expressions in many cases.
2 Approximating the Covariance Function: Retaining finitely many terms is justified by increasing Laplacian eigenvalues and rapidly decaying spectral density at higher frequencies.The approximation is expected to be better away from the boundary of the compact domain.
2 Approximating the Covariance Function: For a one-dimensional example, the squared exponential approximation is indistinguishable from the exact covariance already with m = 12 basis functions.The example uses Ω = [−L, L], L = 1, and ℓ = 0.1; less smooth Matérn functions require more basis functions.
2 Approximating the Covariance Function: The same construction can use an inner product weighted by an input density, with accuracy when the weight is nearly constant where input points lie.The weighted formulation is related to Sturm–Liouville theory and gives an alternative to a compact bounded domain.

3 Application of the Method to GP Regression

The method replaces the full GP covariance matrix with a truncated Laplace-eigenfunction approximation, enabling reduced-rank inference and faster hyperparameter learning. It also extends to operator-observation and latent-force models.

3.1 Gaussian Process Regression: The GP regression covariance matrix is approximated as K ≈ ΦΛΦT using m Laplacian eigenfunctions evaluated at the training inputs.The approximate eigenvalues are spectral-density evaluations at Laplacian eigenvalues, and Φ contains the corresponding eigenfunctions.
3.1 Gaussian Process Regression: When n > m, the matrix inversion lemma makes the reduced-rank approximation advantageous for GP prediction.The predictive mean and variance are obtained after rewriting the full n × n inversion in the reduced basis.
3.2 Learning the Hyperparameters: Hyperparameters are learned by maximizing the marginal likelihood, replacing the full covariance matrix with the approximate matrix in the likelihood and derivative expressions.The full-model optimizer requires repeated inversion of an n × n matrix.
3.3 Discussion on the Computational Complexity: O(nm^2) initial computation and O(m^3) marginal-likelihood evaluation follow because the basis is hyperparameter-independent and Cholesky factorization is performed on an m × m matrix.The product ΦTΦ can be computed once and cached with O(m^2) memory, or formed in blocks when n is large.
3.3 Discussion on the Computational Complexity: Conventional sparse GP approximations typically require O(nm^2) time at every marginal-likelihood step and O(nm) storage.The paper attributes this to re-evaluating basis-related quantities on each step.
3.3 Discussion on the Computational Complexity: On regular grids, ΦTΦ becomes diagonal and projections onto the basis can be evaluated with FFT-based methods in O(n log n) time.This connects the approximation to FFT-based GP regression methods.
3.4 Inverse Problems and Latent Force Models: The approximation extends to observations transformed by linear operators and to latent force models through transformed basis functions or operator projections.For latent force models, the relevant inverse operator can be approximated by projection or point collocation.

4 Convergence Analysis

The analysis establishes convergence of the Laplace-eigenfunction covariance approximation on bounded domains, including multivariate extensions and truncation-rate results. Posterior mean and covariance converge uniformly in the joint limit under stated growth conditions, while broader domain and input-density analyses remain open.

Univariate convergence: The approximation converges uniformly to the true covariance in the joint limit m, L →∞ when m/L →∞.The limits cannot be exchanged arbitrarily; the series must grow sufficiently fast relative to the domain size.
Posterior convergence: Uniform convergence of the prior covariance implies uniform convergence of the posterior mean and covariance as m, L →∞.This follows because GP regression equations use pointwise kernel evaluations.
Multivariate convergence: In d dimensions, uniform posterior convergence holds as m, L1, …, Ld →∞ under the condition m̂/Lk →∞.The result assumes stationary spectral densities satisfying the corresponding one-dimensional conditions in each variable.
Truncation-rate analysis: Selecting a sufficiently high covariance differentiability order can make the truncation convergence rate arbitrarily large; a = d yields approximately 1/m.With a = d/2, the rate reaches the Monte Carlo rate.
Squared exponential covariance: For squared exponential covariance, convergence is faster than approximately 1/m as a function of m, independently of dimensionality, although constants depend on d.The bound has a dimension-dependent term that is linear in d and a term controlling scaling with the number of series terms.
Scope and limitations: The convergence analysis is limited to bounded rectangular domains and is left for future work for spherical domains and input-density-based bases.Spherical domains require Bessel-function eigenbases, while general input densities lack explicit eigenfunctions and eigenvalues.

5 Relationship to Other Methods

The paper contrasts its Laplace-eigenfunction reduced-rank approximation with Nyström-family and direct spectral methods, emphasizing basis construction, probabilistic structure, and boundary effects. Its method uses a fuller spectrum at fixed rank and can achieve accurate approximations when the domain is chosen appropriately.

5.1 Methods from the Nystr¨om Family: Nyström-family methods approximate covariance eigendecompositions using inducing inputs and the eigendecomposition of their covariance submatrix.The inducing-input locations and number determine the approximation.
5.1 Methods from the Nystr¨om Family: SOR, DTC, and VAR differ in which covariance terms they approximate or retain, with DTC lacking a well-formed Gaussian-process interpretation and VAR adding a variational trace term.FIC instead adds a diagonal correction using the true covariance values.
5.1 Methods from the Nystr¨om Family: With 16 inducing points, FIC does not capture the exact squared-exponential correlation structure for the illustrated locations.The passage attributes the mismatch to the number or locations of inducing points being insufficient.
5.2 Direct Spectral Methods: SSGP samples spectral points to span a trigonometric basis, but random spectral coverage can match some contour directions while missing others.The method is interpreted as a Monte Carlo approximation and may require many spectral points, especially in higher dimensions.
5.2 Direct Spectral Methods: The proposed method instead makes the spectrum as full as possible at fixed rank through a quadrature-like Laplace-eigenfunction expansion whose basis is independent of covariance parameters.This provides a contrast with randomly sampled spectral points and inducing-input constructions.
5.2 Direct Spectral Methods: Extending the rectangular domain boundaries by 25% while retaining 16 basis functions gives good contour results, whereas boundary effects distort the process near the edges.The distortion drives the process toward zero at the boundary under Dirichlet conditions.

6 Experiments

Experiments evaluate convergence, accuracy, computational cost, and comparisons across simulated, spatial, spherical, and large-scale real datasets. The reduced-rank method approaches full-GP performance while substantially reducing evaluation time in several settings.

Experimental design: The experiments test convergence, practical use, and comparisons with other methods on simulated and real-world data, including spherical inputs and a very large dataset.Evaluation uses simulated one-dimensional data, real-world data, a spherical input domain, and a large dataset.
Toy data: For the toy dataset, m = 32 produces likelihood curves that are almost exact, while mean and uncertainty estimates differ from the full GP mainly near boundaries.The likelihood differs for small length-scales, roughly below 2L/m; the 95% confidence region differs mainly near the boundaries.
Toy data: m = 20 basis functions are enough for the proposed approximation to converge to the exact full GP solution on the toy data.The convergence is reported for both SMSE and MSLL once enough eigenfunctions account for the short length-scales.
Precipitation data: On the precipitation data, higher ranks converge to full-GP performance in SMSE and MSLL while retaining computational times comparable to conventional methods.With two input dimensions, the method requires more basis functions; spatial medium-size GP regression can be solved in seconds.
Precipitation data: m = 1728 yields precipitation interpolations practically identical to the full GP, using slightly less than 30 seconds instead of approximately 18 minutes.The reduced-rank timing includes initialization and hyperparameter learning.
Large-scale data: On airline arrival delays, Hilbert-GP performs practically equally to VFF on million-scale datasets while completing the full computation in 41 ± 2 seconds.The reported computation includes hyperparameter training and prediction; VFF takes 265 ± 6 seconds in the comparison.

7 Conclusion and Discussion

The paper introduces a Laplace-eigenbasis approximation whose basis functions are independent of covariance hyperparameters, enabling efficient reduced-rank GP inference and hyperparameter learning. It provides truncation and boundary-error analysis, with convergence unaffected by input dimensionality under suitable smoothness conditions.

The method approximates covariance eigendecompositions using Laplace-operator eigenfunctions and the covariance spectral density.
The basis functions are independent of covariance hyperparameters, while inference costs O(nm^2) initially and O(m^3) for marginal-likelihood evaluation.
The approximation supports theoretical analysis of truncation and boundary errors and enables approximation of eigenvalue-based learning-curve bounds.
Boundary conditions introduce edge effects, but the paper states these are known, can be accounted for, and are minor in practical use.
Increasing covariance differentiability appropriately keeps the convergence rate at ∼1/m^a regardless of input dimensionality.
For the squared exponential covariance function, convergence is always better than ∼1/m independently of input dimensionality.

A.1 Auxiliary Lemmas

The appendix proves convergence of the covariance expansion by bounding domain and finite-truncation errors under regularity assumptions on the spectral density. The multidimensional argument extends the one-dimensional bounds through coordinatewise decomposition and triangle inequalities.

The appendix establishes auxiliary lemmas for convergence proofs, relying on classical Riemannian-sum convergence results that are proved for completeness.
Under boundedness, integrability, and derivative assumptions on the spectral density, lemmas bound errors from replacing integrals with infinite spectral sums.
Additional lemmas bound terms involving cosine factors and show constants can be chosen independently of the truncation index under stated assumptions.
Finite truncation is analyzed separately, producing a bound on the residual between infinite- and m-term covariance approximations.
The one-dimensional convergence theorem combines infinite-sum and finite-truncation bounds using the triangle inequality, then takes L →∞ for the limiting result.
For d dimensions, the proof separates coordinatewise sums, repeats the one-dimensional bounds, and combines the resulting terms with telescoping sums and triangle inequalities.