Hilbert Space Gaussian Process Approximation

• Hilbert space Gaussian process approximation is a framework that represents infinite-dimensional Gaussian measures using finite basis expansions and conditional observations.
• It unifies methods such as Laplace-eigenfunction truncations, optimal low-rank posterior updates, and function-space variational inference to balance computational efficiency with statistical accuracy.
• These techniques enable scalable GP inference and uncertainty quantification by preserving critical geometric and probabilistic structures in high-dimensional and inverse problem settings.

Hilbert space Gaussian process approximation denotes a family of methods that represent a Gaussian process prior or posterior as a Gaussian measure on a separable Hilbert space and then approximate that infinite-dimensional object by finite-dimensional subspaces, finitely many linear observations, or finite-rank operator updates. In the literature, the term encompasses at least four closely related strands: Laplace-eigenfunction truncations for stationary kernels used in scalable computation, approximation of conditional Gaussian measures by finite-dimensional conditioning, optimal low-rank posterior updates for linear Gaussian inverse problems, and function-space variational or Langevin methods on reproducing kernel Hilbert spaces (RKHSs) (Riutort-Mayol et al., 2020, Steinwart, 2024, Carere et al., 2024, Wild et al., 25 Feb 2025).

1. Function-space formulation and scope

A Hilbert-space-valued Gaussian random variable is a measurable map $X:\Omega\to H$ such that every scalar projection $\langle X,h\rangle_H$ is Gaussian for all $h\in H$, where $H$ is a separable Hilbert space. Its covariance operator is self-adjoint, positive, and trace-class; in the Banach-space generalization used for conditioning theory, covariance operators are nuclear (Steinwart, 2024). For jointly Gaussian $X$ and $Y$, the pair $(X,Y)$ is Gaussian in the product space, and conditioning becomes an operator-theoretic problem rather than a purely matrix-valued one.

For continuous-path Gaussian processes on a compact metric index set $T$, one may identify the process with a $\mathcal C(T)$-valued Gaussian random variable. The mean function $m(t)=\mathbb E X_t$ and covariance kernel $\langle X,h\rangle_H$0 then correspond to a Gaussian measure on a function space, and the associated RKHS $\langle X,h\rangle_H$1 provides the canonical Cameron–Martin geometry (Steinwart, 2024).

A distinct but compatible RKHS-based formulation appears in function-space Bayesian inference. There one starts from a Gaussian random element $\langle X,h\rangle_H$2 on an RKHS $\langle X,h\rangle_H$3, with covariance operator

$\langle X,h\rangle_H$4

and the induced Gaussian process has kernel

$\langle X,h\rangle_H$5

An important caveat is that a Gaussian process with kernel $\langle X,h\rangle_H$6 does not define a measure on $\langle X,h\rangle_H$7 in infinite dimensions: Driscoll’s theorem says sample paths of a GP with kernel $\langle X,h\rangle_H$8 almost surely lie outside $\langle X,h\rangle_H$9 if $h\in H$0 is infinite-dimensional (Wild et al., 25 Feb 2025). This distinction is central in function-space treatments of Hilbert space GP approximation.

Strand	Core approximation object	Representative papers
Basis-expansion HSGP	Stationary kernel via Laplace eigenfunctions and spectral density	(Riutort-Mayol et al., 2020, Mukherjee et al., 22 May 2025, Zhu et al., 22 Apr 2026)
Conditioning approximation	Conditional Gaussian laws from finitely many linear observations	(Steinwart, 2024, Winkle et al., 19 Aug 2025)
Low-rank posterior update	Posterior mean/covariance via finite-rank operators in Cameron–Martin directions	(Carere et al., 2024, Carere et al., 31 Mar 2025)
Function-space variational inference	RKHS-valued Langevin diffusion projected onto KL coordinates	(Wild et al., 25 Feb 2025)

This multiplicity of meanings is not merely terminological. It reflects a common principle: an infinite-dimensional Gaussian prior or posterior is approximated through a carefully chosen finite-dimensional structure while preserving Gaussianity, support, or measure equivalence as far as possible.

2. Basis-expansion Hilbert space Gaussian processes

In the computational literature, a Hilbert space Gaussian process is typically a stationary GP approximated by a finite basis of Laplace eigenfunctions on a bounded domain. For $h\in H$1, the Dirichlet eigenproblem for $h\in H$2 yields

$h\in H$3

If $h\in H$4 is stationary with spectral density $h\in H$5, then the kernel is approximated by

$h\in H$6

and the GP admits the weight-space representation

$h\in H$7

(Riutort-Mayol et al., 2020, Mukherjee et al., 22 May 2025).

This converts GP inference into Bayesian linear regression with deterministic basis functions and Gaussian coefficients. In one dimension, the dominant per-iteration cost becomes $h\in H$8; in $h\in H$9 dimensions, with $H$0, it becomes $H$1, in contrast with the $H$2 Cholesky cost of exact GP regression (Riutort-Mayol et al., 2020). For multi-output models with $H$3 inputs, $H$4 outputs, and $H$5 basis functions, the exact cost is roughly $H$6, whereas the HSGP approximation yields $H$7 (Mukherjee et al., 22 May 2025).

A practical issue is the interaction between basis size and the artificial domain boundary. Writing $H$8 and $H$9, the boundary factor $X$0 must be large enough to avoid Dirichlet artefacts, but larger $X$1 requires more basis functions. For one-dimensional kernels, the practical recommendations in Stan-based HSGP implementations include

$X$2

$X$3

$X$4

subject to corresponding lower bounds on $X$5 (Riutort-Mayol et al., 2020). In the multi-output latent-variable extension, analogous heuristics are written as

$X$6

(Mukherjee et al., 22 May 2025).

These basis-expansion constructions have been generalized beyond scalar regression. One formulation starts from $X$7 independent scalar GPs $X$8, approximates each with its own HSGP expansion, and then couples outputs by a correlation matrix $X$9. The same framework also accommodates latent inputs $Y$0 through a measurement model $Y$1, with HMC inference implemented in Stan (Mukherjee et al., 22 May 2025). The cited simulations report that these approximate Gaussian processes were not only faster, but also provides similar or even better uncertainty calibration and accuracy of latent variable estimates compared to exact Gaussian processes (Mukherjee et al., 22 May 2025). A plausible implication is that the eigenfunction truncation is especially effective when the computational bottleneck lies in repeated dense-kernel factorizations rather than in low-dimensional output coupling.

3. Conditioning by finite information and convergence of conditional laws

A second major strand treats Hilbert space GP approximation as approximation of conditional Gaussian measures by conditioning on finitely many functionals. For jointly Gaussian $Y$2 and $Y$3, with $Y$4 a separable Banach space, the conditional law is Gaussian: $Y$5 with

$Y$6

$Y$7

where $Y$8 is the Moore–Penrose pseudoinverse (Steinwart, 2024).

For infinite-dimensional observations $Y$9, direct inversion of $(X,Y)$0 is generally impossible because the covariance operator is compact and its range need not be closed. The remedy is a filtering sequence $(X,Y)$1, with $(X,Y)$2, such that the $(X,Y)$3-algebras generated by the $(X,Y)$4 approximate $(X,Y)$5. Then the finite-dimensional conditionals

$(X,Y)$6

are Gaussian, their conditional means converge strongly in $(X,Y)$7, their covariance operators converge in nuclear norm, and the conditional measures converge weakly to the infinite-dimensional conditional law $(X,Y)$8 (Steinwart, 2024). In Hilbert spaces, a natural choice is

$(X,Y)$9

for an orthonormal basis $T$0, while in RKHS settings one may use evaluations on a countable dense subset.

For continuous Gaussian processes $T$1 on compact metric $T$2, conditioned on the restriction of the path to a closed subset $T$3, this approximation has a concrete GP form. If $T$4 is an increasing dense sequence, then conditioning on $T$5 yields posterior mean and covariance

$T$6

$T$7

with uniform convergence on $T$8 and $T$9, respectively. On the observed set $\mathcal C(T)$0, the conditional process interpolates the observation: $\mathcal C(T)$1 for all $\mathcal C(T)$2 and $\mathcal C(T)$3 (Steinwart, 2024). This yields a rigorous continuum-limit interpretation of classical GP regression formulas.

A complementary result links covariance approximation directly to approximation of realizations in Banach or Hilbert norm. If

$\mathcal C(T)$4

and

$\mathcal C(T)$5

then, under the stated regularity assumptions on $\mathcal C(T)$6,

$\mathcal C(T)$7

(Winkle et al., 19 Aug 2025). In $\mathcal C(T)$8, the covariance error has the exact representation

$\mathcal C(T)$9

This shows that posterior variance decay controls approximation of the entire random element, not only its finite-dimensional marginals.

4. Low-rank posterior approximation in Hilbert-space inverse problems

In linear Gaussian inverse problems on a separable Hilbert space $m(t)=\mathbb E X_t$0, one observes

$m(t)=\mathbb E X_t$1

with Gaussian prior

$m(t)=\mathbb E X_t$2

where $m(t)=\mathbb E X_t$3 is trace-class, self-adjoint, positive, and nondegenerate (Carere et al., 2024). The posterior is Gaussian,

$m(t)=\mathbb E X_t$4

with

$m(t)=\mathbb E X_t$5

$m(t)=\mathbb E X_t$6

and

$m(t)=\mathbb E X_t$7

Since $m(t)=\mathbb E X_t$8 has rank at most $m(t)=\mathbb E X_t$9, the prior-to-posterior update acts only on a finite-dimensional subspace (Carere et al., 2024).

The infinite-dimensional difficulty is that not every low-rank covariance approximation yields a Gaussian measure equivalent to the posterior or prior. The Feldman–Hájek theorem therefore becomes structural rather than auxiliary. Within the admissible family

$\langle X,h\rangle_H$00

and the corresponding precision family

$\langle X,h\rangle_H$01

the optimal rank-$\langle X,h\rangle_H$02 covariance and precision approximations are

$\langle X,h\rangle_H$03

$\langle X,h\rangle_H$04

where $\langle X,h\rangle_H$05 arise from the spectral analysis of the covariance ratio or prior-preconditioned Hessian (Carere et al., 2024). These approximations are simultaneously optimal for a large class of spectral losses, including the Kullback–Leibler divergence, the Rényi divergences, the Amari $\langle X,h\rangle_H$06-divergences for $\langle X,h\rangle_H$07, and the Hellinger metric (Carere et al., 2024).

Posterior mean approximation requires additional care because equivalence of Gaussian measures imposes Cameron–Martin constraints on mean shifts. Two families are considered: a structure-preserving class and a structure-ignoring class. The optimal structure-ignoring approximation is

$\langle X,h\rangle_H$08

with minimal average loss

$\langle X,h\rangle_H$09

The optimal structure-preserving approximation is

$\langle X,h\rangle_H$10

with minimal average loss

$\langle X,h\rangle_H$11

(Carere et al., 31 Mar 2025). Both are simultaneously optimal for averaged KL, Rényi, Amari, and Hellinger losses when covariance is held fixed.

For reverse KL, the optimal low-rank covariance $\langle X,h\rangle_H$12 and the optimal low-rank mean $\langle X,h\rangle_H$13 combine to give an optimal joint approximation of mean and covariance. Moreover, the joint approximation with $\langle X,h\rangle_H$14 admits a parameter-space projector interpretation: it is the exact posterior of a projected inverse problem with forward operator $\langle X,h\rangle_H$15, where $\langle X,h\rangle_H$16 projects onto the most data-informed directions in $\langle X,h\rangle_H$17 (Carere et al., 31 Mar 2025). This identifies a likelihood-informed reduced basis for GP posterior approximation in function space.

5. Function-space variational inference and projected Langevin dynamics

A different formulation starts from the full posterior on an RKHS $\langle X,h\rangle_H$18 as the minimizer of the variational functional

$\langle X,h\rangle_H$19

where $\langle X,h\rangle_H$20 is a Gaussian random element on $\langle X,h\rangle_H$21 and $\langle X,h\rangle_H$22 is the negative log-likelihood (Wild et al., 25 Feb 2025). The corresponding Wasserstein gradient flow leads to the $\langle X,h\rangle_H$23-valued Langevin diffusion

$\langle X,h\rangle_H$24

whose stationary distribution is the Bayes posterior under the stated assumptions (Wild et al., 25 Feb 2025).

To approximate this infinite-dimensional dynamics, one projects onto the first $\langle X,h\rangle_H$25 terms of the Kosambi–Karhunen–Loève expansion

$\langle X,h\rangle_H$26

and keeps only

$\langle X,h\rangle_H$27

The projected coefficient process is an $\langle X,h\rangle_H$28-valued Langevin diffusion with stationary density proportional to $\langle X,h\rangle_H$29, where

$\langle X,h\rangle_H$30

Here $\langle X,h\rangle_H$31 and $\langle X,h\rangle_H$32 is the evaluation of the truncated KL expansion at the data sites (Wild et al., 25 Feb 2025).

The induced approximate posterior on function space is built by the law of total probability and a sufficiency assumption. Define

$\langle X,h\rangle_H$33

Then the exact Bayes posterior can be written through the marginal posterior of $\langle X,h\rangle_H$34, and the approximation replaces that marginal by the stationary law of the projected Langevin system while assuming

$\langle X,h\rangle_H$35

(Wild et al., 25 Feb 2025). This yields a nonparametric variational family on $\langle X,h\rangle_H$36, not merely a Gaussian family.

The central theoretical statement is a near-optimality bound. If $\langle X,h\rangle_H$37 is convex and $\langle X,h\rangle_H$38-Lipschitz, then the projected Langevin posterior $\langle X,h\rangle_H$39 satisfies

$\langle X,h\rangle_H$40

and, for i.i.d. inputs,

$\langle X,h\rangle_H$41

(Wild et al., 25 Feb 2025). The tail $\langle X,h\rangle_H$42 is $\langle X,h\rangle_H$43 for general trace-class covariance, $\langle X,h\rangle_H$44 for squared exponential kernels with Gaussian or compactly supported $\langle X,h\rangle_H$45, and $\langle X,h\rangle_H$46 for Matérn kernels of smoothness $\langle X,h\rangle_H$47 under the stated assumptions (Wild et al., 25 Feb 2025).

This framework also clarifies the relationship to sparse variational Gaussian processes. The method recovers the posterior arising from the sparse variational Gaussian process as a special case, owed to the fact that the sufficiency assumption underlies both methods. However, whereas the SVGP is parametrically constrained to be a Gaussian process, this method is based on a non-parametric variational family $\langle X,h\rangle_H$48, consisting of all probability measures on $\langle X,h\rangle_H$49, and coincides with SVGP for the special case of a Gaussian error likelihood (Wild et al., 25 Feb 2025).

6. Error bounds, applications, and interpretive issues

Approximation-theoretic results show that the efficiency of Hilbert-space GP approximation depends strongly on kernel regularity and on the approximation space. For analytic kernels on $\langle X,h\rangle_H$50, including the Gaussian kernel

$\langle X,h\rangle_H$51

linear approximation in the RKHS by function evaluations has $\langle X,h\rangle_H$52-th minimal error bounded above and below by multiples of

$\langle X,h\rangle_H$53

up to sub-exponential factors (Karvonen et al., 2022). This places Gaussian-kernel approximation in a markedly different regime from Sobolev- or Matérn-type kernels, where decay is algebraic rather than super-algebraic.

In sequential design, the same Laplace-eigenfunction machinery has been used to approximate the IMSE acquisition function. For stationary kernels on $\langle X,h\rangle_H$54, a truncated eigenbasis representation yields

$\langle X,h\rangle_H$55

so that the acquisition can be evaluated in closed form without numerical integration (Zhu et al., 22 Apr 2026). The paper establishes sharp global non-asymptotic bounds for both kernel approximation and acquisition approximation: for Gaussian kernels, both aliasing and truncation errors decay exponentially; for Matérn kernels, the truncation error is polynomial in $\langle X,h\rangle_H$56, and balancing the parameters gives rates of order $\langle X,h\rangle_H$57 under the stated conditions (Zhu et al., 22 Apr 2026).

Related work studies Gaussian approximation in Hilbert spaces from a probabilistic rather than Bayesian perspective. Quantitative central limit theorems for $\langle X,h\rangle_H$58-valued random variables use Stein’s method, Gamma calculus, and chaos expansions to control distances between non-Gaussian Hilbert-valued functionals and Gaussian laws. This includes functional Breuer–Major limits in $\langle X,h\rangle_H$59 and Poisson-to-Brownian approximation in Besov–Liouville spaces (Bourguin et al., 2019, Bourguin et al., 2021). Such results are not GP regression methods, but they broaden the mathematical meaning of Hilbert-space Gaussian approximation by showing how Gaussian process limits arise for random elements in infinite-dimensional Hilbert spaces.

Several misconceptions recur in the literature. First, Hilbert space Gaussian process approximation is not a single algorithm: it may refer to Laplace-eigenfunction truncation, approximation of conditional measures, optimal low-rank posterior updates, or RKHS-valued variational inference. Second, a GP with kernel $\langle X,h\rangle_H$60 should not be identified with a Gaussian measure on the RKHS $\langle X,h\rangle_H$61; in infinite dimensions, sample paths almost surely lie outside $\langle X,h\rangle_H$62 (Wild et al., 25 Feb 2025). Third, finite-rank covariance approximation in infinite dimensions is constrained by measure equivalence: many seemingly natural low-rank updates produce Gaussian measures that are mutually singular with the target posterior, which is why Feldman–Hájek conditions play a central role in the inverse-problem literature (Carere et al., 2024). Fourth, the familiar finite-dimensional Gaussian conditioning formula does not transfer verbatim to infinite-dimensional observation spaces; it must be recovered as a limit through filtering sequences or related approximation devices (Steinwart, 2024).

Taken together, these results show that Hilbert space Gaussian process approximation is best understood as a function-space paradigm. The prior or posterior is treated as a Gaussian measure on an infinite-dimensional Hilbert or Banach space, and approximation is carried out by finite basis expansions, finitely many observations, or finite-rank operator modifications. What varies across the literature is the target of approximation—kernel, conditional law, posterior covariance, posterior mean, acquisition function, or entire posterior measure—but the underlying problem remains the same: to replace an infinite-dimensional Gaussian object by a finite-dimensional surrogate while retaining the geometry and probabilistic structure needed for inference, uncertainty quantification, and asymptotic control.