Functional Gaussian Process Regression

• FGPR is a framework where the primary unknown is a functional object or operator rather than a finite-dimensional vector.
• It leverages covariance kernels, basis expansions, and structured relevance measures to encode smoothness and geometric properties.
• The method supports various regression geometries and scalable computational strategies for applications in physics, remote sensing, and PDE modeling.

Searching arXiv for recent and relevant papers on Functional Gaussian Process Regression and closely related formulations. Functional Gaussian Process Regression (FGPR) denotes a family of Gaussian-process-based regression formulations in which the primary unknown is a function, a collection of spatially indexed functions, or an operator acting on function spaces, rather than merely a scalar-valued map on finite-dimensional Euclidean inputs. Across the literature, FGPR appears in several mathematically equivalent or closely related forms: Gaussian processes viewed as priors over latent functions in regression; Gaussian measures on Hilbert spaces represented through basis expansions; Gaussian-process priors over coefficient functions in functional-input or functional-output models; and function-space corrections to mechanistic models such as linear PDEs. The common structure is that inference is carried out over an infinite- or high-dimensional functional object, with covariance kernels or covariance operators encoding smoothness, correlation length, stationarity, geometric invariance, or other structural assumptions (Gandrakota et al., 2022, Farooq et al., 2024, Nguyen et al., 2014).

1. Conceptual scope and formal definitions

A canonical FGPR formulation begins with a latent function $z(\cdot)$ and observations

$y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$

with $z(\cdot)$ assigned a Gaussian-process prior

$z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$

In this sense, FGPR is “Gaussian processes viewed as priors over functions” in a regression setting, where inference concerns the latent function itself rather than a finite vector of coefficients (Gandrakota et al., 2022). This functional interpretation is equally explicit in Hilbert-space approximations, where one works in $L^2(\Omega)$ or a related separable Hilbert space and writes

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$

so that the GP prior becomes a Gaussian distribution over coefficients in a basis of eigenfunctions (Farooq et al., 2024).

The term also covers regression problems in which at least one covariate is itself a function $X(t)$ defined over a continuum $\mathcal T$, so that the regression map acts on a function space rather than $\mathbb R^p$. In that setting, one may write $f:\mathcal X\to\mathcal Y$, where $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$0 is a space of curves and $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$1 is scalar- or function-valued, and define kernels directly on functional arguments through weighted $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$2-type distances or related constructions (Damiano et al., 2022). A different but compatible usage places the GP on functional responses: for replicated curves $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$3, one decomposes

$y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$4

with $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$5 modeled as a GP and $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$6 supplied by a functional regression mean structure (Konzen et al., 2021, Wang et al., 2014).

A further extension treats the GP as a Gaussian random functional in a weak formulation of a PDE. There the unknown object is not a scalar-valued regression function on Euclidean inputs, but a linear functional $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$7 on a Hilbert space $y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$8, inserted into a linear PDE model as a model-discrepancy term and inferred from data through adjoint states (Nguyen et al., 2014). This suggests that FGPR is best understood as a unifying perspective on Gaussian priors over function-space objects, rather than a single model class.

2. Functional priors, kernels, and covariance operators

In FGPR, the kernel or covariance operator determines the admissible regularity class of the latent functional object. In ordinary input-space formulations, standard kernels such as the radial basis function kernel,

$y_k = z(x_k) + \epsilon_k,\quad k=1,\dots,N,$9

the Matérn kernel with $z(\cdot)$0, and polynomial kernels define different smoothness and complexity regimes for the latent background or regression surface (Gandrakota et al., 2022). In the Hilbert-space view, these same priors are represented spectrally through eigenfunctions $z(\cdot)$1 and spectral densities $z(\cdot)$2,

$z(\cdot)$3

making the covariance operator explicit and finite-rank after truncation (Farooq et al., 2024).

For functional inputs, kernels are often constructed from weighted norms over the index domain. A central example is

$z(\cdot)$4

embedded inside a squared exponential covariance

$z(\cdot)$5

Here the weight function $z(\cdot)$6 encodes predictive relevance over the index space, so the kernel is not merely a measure of global curve similarity but a structured relevance-weighted functional distance (Damiano et al., 2022). The asymmetric Laplace functional weight (ALF) further parameterizes $z(\cdot)$7 through a location of peak relevance and asymmetric left/right decay rates, enforcing smooth unimodal relevance profiles with three unknowns per input variable (Damiano et al., 2022).

For functional outcomes with predictors at different scales, covariance modeling becomes additive across components. One example writes

$z(\cdot)$8

placing GP priors on the spatially varying coefficient surfaces $z(\cdot)$9 and on coefficient functions $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$0 inside the expansion

$z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$1

with Matérn covariance kernels on the spatial domain $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$2 (Andros et al., 10 Feb 2026). This produces a covariance that is a sum of separable terms in spatial and predictor-space components rather than a single tensor product (Andros et al., 10 Feb 2026).

On manifolds, covariance construction uses geometry. For wrapped Gaussian process functional regression, the response takes values on a Riemannian manifold $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$3, and Gaussianity is defined in tangent spaces via the exponential and logarithm maps; covariance kernels are specified for tangent-space vector fields and transported back to the manifold through $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$4 and $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$5 (Liu et al., 2024). For spatiotemporal fields on compact homogeneous manifolds, covariance operators are diagonalized in Laplace–Beltrami eigenfunctions, and isotropic kernels admit angular spectra $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$6 in the expansion

$z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$7

which supports time-adaptive truncation and Empirical Bayes estimation in function space (Ruiz-Medina et al., 22 Mar 2026).

3. Regression structures: scalar responses, functional inputs, and functional outputs

FGPR encompasses several distinct regression geometries. The simplest is scalar observation of a latent function, as in background estimation or signal extraction from binned counts. There the latent background is modeled nonparametrically by a GP, while a localized signal can be included through a parametric mean function such as

$z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$8

yielding a signal-plus-background decomposition in which the GP carries the background uncertainty and the mean function carries the structured signal (Gandrakota et al., 2022).

A second class treats inputs as functions. Rather than discretizing $z(x) \sim \mathcal{GP}\big(m(x),k(x,x')\big).$9 into a high-dimensional vector with one relevance parameter per grid point, one may define a kernel directly on curves through $L^2(\Omega)$0 and learn a smooth relevance function over $L^2(\Omega)$1. This replaces high-dimensional automatic relevance determination by automatic dynamic relevance determination, in which posterior inference targets $L^2(\Omega)$2 itself (Damiano et al., 2022). The paper reporting this formulation emphasizes that FPCA-based reduction targets variance in the input rather than predictive relevance for the response, and therefore need not align with the truly predictive regions of the functional domain (Damiano et al., 2022).

A third class treats responses as functions and often decomposes them into mean and residual GP components. In the Gaussian process functional regression framework implemented in GPFDA, one writes

$L^2(\Omega)$3

where $L^2(\Omega)$4 may depend on scalar or functional covariates through standard functional regression terms, and $L^2(\Omega)$5 models residual dependence across the functional domain (Konzen et al., 2021). The generalized extension to non-Gaussian responses replaces the Gaussian likelihood by an exponential-family observation model with inverse link $L^2(\Omega)$6,

$L^2(\Omega)$7

allowing binomial, Poisson, ordinal, and related functional responses while retaining GP structure in the latent layer (Wang et al., 2014).

Function-on-function regression constitutes a more explicitly operator-valued regime. The deep Gaussian process formulation for functional maps takes observed pairs $L^2(\Omega)$8, with $L^2(\Omega)$9 and $$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$0, and models the operator $$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$1 through a sequence of GP-based linear integral transforms and GP-sampled nonlinear activations (Lowery et al., 24 Oct 2025). A linear layer is defined by

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$2

while nonlinear layers use activation functions sampled from GPs. This explicitly treats the regression target as a map between function spaces, rather than a scalar regression surface (Lowery et al., 24 Oct 2025).

Mechanistic FGPR offers another regression structure: instead of regressing directly on function-valued inputs or outputs, it infers a Gaussian random functional $$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$3 augmenting a weak PDE formulation,

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$4

and uses observations of linear functionals of the field to learn $$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$5 through adjoint equations (Nguyen et al., 2014). This formulation embeds regression into operator equations and treats model inadequacy as a stochastic object in the dual space.

4. Likelihoods, posterior inference, and model selection

With Gaussian observations, FGPR retains the standard GP marginalization formulas. In the scalar regression setting with covariance matrix $$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$6, the marginal likelihood is

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$7

and the predictive mean and variance at a new input are given by the standard GP formulas (Gandrakota et al., 2022). These expressions support type-II maximum likelihood, Bayesian Information Criterion, and Akaike Information Criterion, with effective degrees of freedom

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$8

used as a complexity measure for GP regression (Gandrakota et al., 2022).

When responses are non-Gaussian, conjugacy is lost and approximation becomes necessary. The generalized functional GP model for exponential-family responses defines the marginal likelihood

$$f(\mathbf x) \approx \sum_{j=1}^M w_j \phi_j(\mathbf x),\qquad w\sim \mathcal N(0,\Lambda),$$9

and then uses Laplace-type approximations and INLA-like Gaussian approximations to integrate out the latent GP $X(t)$0 (Wang et al., 2014). This preserves the hierarchical GP structure while allowing Bernoulli, Poisson, and ordinal functional data.

Bayesian estimation also appears in functional-input kernels with learned relevance profiles. The ALF-based ADRD model places priors on $X(t)$1, evaluates the log marginal likelihood of the GP covariance matrix built from $X(t)$2, and uses a fully Bayesian workflow consisting of random initialization, local optimization, and NUTS sampling to obtain posterior relevance profiles $X(t)$3 with credible intervals (Damiano et al., 2022). The paper also introduces Permutation Feature Dynamic Importance as a screening and diagnostic tool for assessing whether the unimodal ALF assumption is broadly consistent with the data (Damiano et al., 2022).

Posterior inference in hierarchical functional-output models remains Gaussian when the observation model is Gaussian. In the multi-scale functional-outcome model, integrating out all latent coefficient functions yields

$X(t)$4

with

$X(t)$5

so that prediction proceeds by standard conditioning in a joint multivariate normal distribution (Andros et al., 10 Feb 2026). In the manifold-valued response setting, the wrapped GP predictive distribution takes the form

$X(t)$6

combining Euclidean GP inference in tangent coordinates with geometric reconstruction on $X(t)$7 (Liu et al., 2024).

5. Computational strategies and scalability

The principal computational obstacle in FGPR is covariance manipulation at the scale induced by discretized functional data. Exact GP regression typically requires $X(t)$8 factorization of dense covariance matrices, which becomes prohibitive when each observation is itself a finely sampled function (Farooq et al., 2024, Hoffmann et al., 2024). Several distinct strategies have emerged.

A first strategy is basis truncation in Hilbert space. The Hilbert-space Gaussian process approximation replaces the full $X(t)$9 kernel matrix by an $\mathcal T$0 system in basis coefficients,

$\mathcal T$1

with complexity scaling as $\mathcal T$2 or $\mathcal T$3 when $\mathcal T$4 (Farooq et al., 2024). The quantum-assisted variant preserves this same functional formulation and applies qPCA, conditional rotations, and Hadamard and Swap tests to accelerate evaluation of the posterior mean and variance in coefficient space, reporting total complexity

$\mathcal T$5

with data loading included (Farooq et al., 2024).

A second strategy exploits sampling design. In the multi-level functional GP model for repeated functions, completely regular designs induce a covariance of the form

$\mathcal T$6

which allows exact analytic expressions for the log-likelihood and posterior in terms of two $\mathcal T$7 matrices rather than one $\mathcal T$8 matrix (Hoffmann et al., 2024). This reduces log-likelihood evaluation from $\mathcal T$9 to essentially $\mathbb R^p$0, and empirical benchmarks report speedups of $\mathbb R^p$1–$\mathbb R^p$2 for log-likelihood computation and $\mathbb R^p$3–$\mathbb R^p$4 for posterior simulation under regular sampling designs (Hoffmann et al., 2024).

A third strategy is low-rank approximation by predictive processes. In function-on-function regression with functional responses, a predictive-process approximation replaces the full GP by its conditional expectation given values at a set of knot covariates and knot times, reducing the effective covariance rank from $\mathbb R^p$5 to $\mathbb R^p$6 (Lian, 2010). The paper stresses that unmodified predictive processes can severely underestimate predictive variance in this setting and therefore introduces diagonal corrections that restore valid uncertainty quantification (Lian, 2010).

A fourth strategy appears in deep functional GP maps. There, a key observation is that when projection and quadrature nodes coincide, discrete approximations to kernel integral transforms collapse into direct linear transforms on discretized function values, eliminating the need to explicitly propagate induced kernels through all layers (Lowery et al., 24 Oct 2025). Combined with inducing-point variational inference and whitening, this yields tractable deep probabilistic operator learning on irregular functional data (Lowery et al., 24 Oct 2025).

Time-adaptive manifold FGPR uses yet another form of reduction: covariance operators are diagonal in Laplace–Beltrami eigenfunctions, and truncation of the angular spectrum to $\mathbb R^p$7 terms converts infinite-dimensional posterior updating into a finite set of scalar spectral updates (Ruiz-Medina et al., 22 Mar 2026). This suggests a general principle: when the geometry or sampling design induces a natural diagonalization or Kronecker structure, exact or near-exact FGPR can be made computationally feasible without abandoning the underlying probabilistic model.

6. Applications, extensions, and recurrent points of comparison

FGPR has been applied to weak-signal extraction in high-energy physics, where a GP background model separated Higgs-resonance structure from a complex smooth background and yielded a signal estimate $\mathbb R^p$8 events at $\mathbb R^p$9 GeV and $f:\mathcal X\to\mathcal Y$0 GeV (Gandrakota et al., 2022). In remote sensing, ALF-based functional-input GPs were used with atmospheric profile covariates from NASA’s Microwave Limb Sounder, where the resulting functional-input GP models were competitive with full ARD and clearly superior to FPCA-based models while using far fewer relevance parameters (Damiano et al., 2022). In computer experiments and coastal flooding, a functional GP prior on spatially indexed nonlinear effects of global predictors was demonstrated on synthetic datasets and outputs from the Sea, Lake, and Overland Surges from Hurricanes model (Andros et al., 10 Feb 2026).

Manifold-valued FGPR extends the framework to non-Euclidean responses such as trajectories on the sphere or Kendall shape space. The wrapped Gaussian process functional regression model combines a Fréchet mean curve, tangent-space functional regression on scalar batch covariates, and GP covariance structure driven by functional covariates, and reports smaller RMSE than functional linear regression on manifolds and a wrapped GP with only a constant Fréchet mean in interpolation and extrapolation tasks (Liu et al., 2024). On functional maps between spaces, deep GP architectures have been evaluated on Burgers, Darcy, Navier–Stokes, Beijing-Air, SLC-Precipitation, and Quasar reverberation mapping, with reported advantages in both predictive performance and uncertainty calibration (Lowery et al., 24 Oct 2025).

Several recurrent comparison points appear across the literature. Against polynomial or finite-basis parametric fits, FGPR is repeatedly described as more flexible because uncertainty is placed over an infinite-dimensional function space rather than a small coefficient vector (Gandrakota et al., 2022, Nguyen et al., 2014). Against FPCA or basis truncation used only for dimension reduction, FGPR variants that model relevance directly over the original index space are presented as more interpretable when the scientific question concerns which regions of a function matter for prediction (Damiano et al., 2022). Against neural operators or deterministic function-on-function learners, deep functional GP models are positioned as offering calibrated uncertainty under noisy, sparse, and irregular sampling (Lowery et al., 24 Oct 2025). Against classical data assimilation or model calibration, the PDE-constrained formulation differs by placing a GP prior on a correction functional in weak form and learning it from observations through adjoints (Nguyen et al., 2014).

A common misconception is that FGPR refers only to one of these subfamilies. The surveyed literature suggests instead that the term spans a spectrum of constructions unified by Gaussian priors over function-space objects. A plausible implication is that kernel design, basis choice, or covariance-operator structure should be viewed as the central modeling decision, while the distinction between “input-space GP,” “Hilbert-space GP,” “functional-input GP,” “functional-output GP,” and “operator-valued GP” is often a matter of representation rather than of underlying probabilistic principle (Farooq et al., 2024, Konzen et al., 2021). Another recurrent point is that computational tractability is not intrinsic to or incompatible with FGPR: exact inference is feasible in some structured settings, whereas in others low-rank, spectral, inducing-point, or geometry-aware approximations are the decisive ingredient (Hoffmann et al., 2024, Farooq et al., 2024, Lowery et al., 24 Oct 2025).