Hierarchical Gaussian Processes

• Hierarchical Gaussian processes are probabilistic models that compose multiple GPs via hyperpriors to capture complex dependencies and non-stationarity.
• They enable scalable inference and improved uncertainty quantification, outperforming traditional methods in large-scale regression and multi-task settings.
• Key applications include Bayesian deep learning, surrogate modeling, and advanced optimization leveraging structured kernels and mixture-of-experts.

A hierarchical Gaussian process (HGP) is a family of probabilistic models in which multiple Gaussian processes are composed, stacked, or jointly structured via hyperpriors, enabling the modeling of complex dependencies, multi-level uncertainty, non-stationarity, multi-output/replicated structures, or scalable inference in large datasets. HGPs generalize classical GPs and arise naturally in Bayesian deep learning, scalable regression, multi-task settings, structure-exploiting surrogates, and advanced Bayesian optimization. Hierarchical structure enables sharing of statistical strength across tasks, input domains, outputs, or latent features, yielding flexible representation of uncertainty and complex data hierarchies.

1. Structural Taxonomy and Core Model Classes

Hierarchical Gaussian processes encompass several structurally distinct modeling paradigms, unified by multi-level or compositional structure:

• Deep/Stacked GPs: Compositions $f_L(\cdot) = f^L( f^{L-1}( ... f^1(\cdot) ))$, where each $f^l$ is a GP, yielding nested stochastic mappings and non-Gaussian marginals (Havasi et al., 2018, Zhao et al., 2020, Popescu et al., 2020, Popescu et al., 2022).
• Hierarchical priors (hyper-GPs): Covariance or kernel hyperparameters of one or multiple GPs are themselves random variables or latent functions drawn from separate GPs, enabling random, spatially-varying kernel structure (Karaletsos et al., 2020, Monterrubio-Gómez et al., 2018, Zhao et al., 2020).
• Hierarchical shrinkage: GP prior weights in a basis expansion are endowed with structured shrinkage hyperpriors to induce effect-sparsity, heredity, or hierarchy in main/interaction terms (Tang et al., 2023).
• Multi-output and tree-structured HGPs: Multiple related outputs or replicate groups are modeled jointly via hierarchical or tree-structured kernels encapsulating within- and between-output dependencies, often combined with latent variable parameterizations (Ma et al., 2023).
• Hierarchical mixtures-of-experts: Conditional independence between local GPs is assumed and aggregation is performed via product or mixture operators in a tree, scaling exact or approximate inference to massive datasets (Ng et al., 2014, Liu et al., 2023).
• Hierarchical surrogates for mixed/hierarchical inputs: GP kernels encode variable “roles” in a hierarchy (meta, neutral, decreed), with structure-aware covariance models (Saves et al., 2023).
• Hierarchical matrix methods and scalable algebra: Covariance matrices are decomposed recursively (e.g., HMAT) to enable O$(n\log n)$ computation for massive datasets without explicit prior stacking (Keshavarzzadeh et al., 2021).

This taxonomy is not exhaustive; many hybrid models combine features across these axes.

2. Mathematical Formulations and Hierarchical Constructions

• Deep (L-layer) GPs: Hierarchy is established by recursive functional composition:

$$f(\mathbf{x}) = f^L(f^{L-1}(\cdots f^1(\mathbf{x}) \cdots))$$

Each $f^l$ is endowed with an independent or conditional GP prior, frequently parameterized via inducing points $Z^l, U^l$ for variational inference (Havasi et al., 2018, Zhao et al., 2020, Popescu et al., 2020).

• Hierarchical kernel priors: Parameters such as length-scales or amplitudes receive structured priors, e.g.,

$$\ell(x) = g(u(x)), \quad u(x) \sim \mathcal{GP}(m_\ell, k_\ell)$$

leading to covariance $k(x, x')$ with spatially varying hyperparameters. The overall model becomes

$$f \sim \mathcal{GP}(0, k(x,x';\ell(\cdot),\sigma(\cdot)))$$

where $\ell(\cdot),\sigma(\cdot)$ are themselves latent GPs (Monterrubio-Gómez et al., 2018, Zhao et al., 2020).

• Weight-space hierarchical GPs: For BNNs, weight matrices are generated by hierarchical GPs, e.g., latent unit embeddings $f^l$0 carry GP priors, with weights $f^l$1 drawn from another GP (Karaletsos et al., 2020).
• Shrinkage priors in function-space GPs: Function expansion

$f^l$2

with $f^l$3 governed by hierarchical, stick-breaking (e.g., cumulative shrinkage) or global-local (horseshoe) priors, enforces structured sparsity and effect heredity (Tang et al., 2023).

• Multi-output hierarchical kernels: For repeated measurements/replicates,

$f^l$4

with latent-variable coregionalization $f^l$5 for inter-output associations (Ma et al., 2023).

• Hierarchical mixture-of-experts aggregation: Local GP predictors at leaves are recursively aggregated. If child experts output $f^l$6, internal nodes compute

$f^l$7

(Ng et al., 2014).

3. Inference Methods and Computational Frameworks

Hierarchical architectures present unique inference and scalability challenges. Several principal methodologies include:

• Sparse Variational Inference: Deep and hierarchical GPs utilize inducing points and variational posteriors, allowing stochastic ELBO optimization and scalable training. Decoupled inducing sets for the predictive mean and covariance break previous bottlenecks and facilitate deep architectures (Havasi et al., 2018).
• State-Space Methods for Time Series: Deep hierarchies of conditional GPs can be cast as non-linear SDEs, permitting O$f^l$8 filtering and smoothing algorithms (CKF, EKF, particle smoothers) for potentially thousands of time steps (Zhao et al., 2020).
• Hierarchical MCMC: Hierarchical shrinkage, non-stationary, and logistic models deploy block Gibbs, marginal elliptical slice sampling, and adaptive random-walk Metropolis steps, exploiting sparse or banded structure to address posterior coupling and tuning (Tang et al., 2023, Monterrubio-Gómez et al., 2018, Rau et al., 2019).
• Distributed Training and Aggregation: Hierarchical mixture-of-experts models are structured for parallelism; closed-form aggregation at internal nodes and independent, small local GP inferences permit scaling to $f^l$9 on modern clusters (Ng et al., 2014).
• Matrix Decomposition Algorithms: Hierarchical low-rank matrix decompositions (e.g., HMAT) achieve O$(n\log n)$0 per-iteration cost for kernel inversion and determinant evaluation by recursively compressing off-diagonal blocks; randomized SVD and interpolative decompositions are critical enablers (Keshavarzzadeh et al., 2021).
• Tailored Kernel Parameterization for Hierarchical Inputs: Surrogate modeling frameworks implement kernels that couple meta, decreed, and categorical variables via algebraic or arc distance computations, with analytic derivatives for hyperparameter learning (Saves et al., 2023).
• Random Feature Expansions in MoE Models: Tree-structured mixtures with GP-based nonlinear gating are tractably approximated by random Fourier features and trained with doubly-stochastic variational inference (Liu et al., 2023).

4. Applications and Empirical Behavior

Hierarchical GPs have demonstrated superiority over unstructured and shallow counterparts in a variety of substantive domains:

• Large-scale Regression and Surrogate Modeling: Hierarchical GP-MoEs outperform sparse and full GPs on million-point regression, with parallel/distributed aggregation and no explicit sparsification (Ng et al., 2014).
• Bayesian Neural Networks: HGP priors over weights in BNNs (MetaGP) yield improved uncertainty quantification, principled extrapolation (e.g., recovery of periodic structure), and calibrated out-of-distribution detection compared to mean-field and deep kernel learning methods (Karaletsos et al., 2020).
• Multi-task and Replicated Data: Latent-variable multi-output HGPs dominate flat MOGPs on hierarchically replicated genomics and motion-capture benchmarks, with untied levels for intra- and inter-replicate dependencies (Ma et al., 2023).
• Non-stationary and Structured Sparsity Modeling: HGPs with non-stationary priors or shrinkage assign model capacity adaptively, capturing abrupt changes, structured sparsity, and effect heredity not available to standard or purely additive GPs (Monterrubio-Gómez et al., 2018, Tang et al., 2023).
• Out-of-Distribution Detection and Image Segmentation: Hierarchical convolutional GPs employing Wasserstein-2 kernels propagate distributional uncertainty, allowing OOD detection in medical imaging and achieving near state-of-the-art segmentation performance, outperforming Bayesian CNNs and deep GPs without distributional layers (Popescu et al., 2022, Popescu et al., 2020).
• Surrogate Optimization with Hierarchical and Mixed Inputs: Design-space-aware hierarchical Kriging surrogates outperform imputation-based or mixed-integer methods in regression and global optimization with hierarchical search spaces (Saves et al., 2023).
• Bayesian Optimization across Heterogeneous Search Spaces: Universal pre-trained hierarchical GP priors (HyperBO+) enable transfer learning and optimization across varying input domains, outperforming baselines not capable of domain adaptation (Fan et al., 2022).

5. Comparison with Classical and Deep GP Frameworks

• Relation to Deep GPs: Deep GPs are a specific case of hierarchical GPs, where GPs are composed via functional mappings, as opposed to hyperparameter hierarchies or mixture-of-expert aggregations. Deep GPs provide universal function approximation, non-Gaussianity, input warping, and layered uncertainty, but suffer from inference complexity and potential variance collapse in hidden layers without appropriate kernel or mean design (Havasi et al., 2018, Popescu et al., 2020, Popescu et al., 2022).
• Hierarchical Priors vs. Shallow GPs: Hierarchical GPs enable adaptive, context-dependent uncertainty and structured correlation across outputs or latent groups, which classical GPs with fixed or shallow kernels cannot represent. They naturally implement information-pooling across tasks, outputs, or sparse signals (Ma et al., 2023, Tang et al., 2023).
• Scalability and Parallelism: Hierarchical MoE and matrix-decomposition HGPs scale exact or nearly-exact GP performance to millions of points, contrasting with the limitations of standard covariance inversion or finite-inducing representations (Ng et al., 2014, Keshavarzzadeh et al., 2021).
• Latent Variable and Structured Input Representations: HGPs leveraging latent-variable kernels or input-aware representations (e.g., for mixed/hierarchical variables) can encode complex, structured relationships in a manner not accessible to flat kernel designs or simple additive structures (Ma et al., 2023, Saves et al., 2023).

6. Limitations and Current Research Directions

Hierarchical GPs, while flexible and expressive, entail unique statistical and computational challenges:

• Model Selection and Hyperparameter Complexity: Deep or hierarchical models risk overfitting, non-identifiability, and hyperparameter non-robustness, especially in low-data or deep-architecture settings (Rau et al., 2019, Tang et al., 2023).
• Inference Scalability: While mixture and low-rank strategies scale well, fully Bayesian inference in deep or hyper-GP models may still be hindered by O$(n\log n)$1 or $(n\log n)$2-layered kernel complexity, requiring careful exploitation of Kronecker or banded structure (Havasi et al., 2018, Monterrubio-Gómez et al., 2018, Keshavarzzadeh et al., 2021).
• Extension to Non-Gaussian Likelihoods and Arbitrary Hierarchies: Many current hierarchical GP frameworks are limited to regression or Gaussian likelihoods. Work on deep/hierarchical GPs for non-Gaussian and structured output spaces is ongoing (Ma et al., 2023, Liu et al., 2023).
• Automatic Structure Discovery: Model architectures (depth, kernel choice, hierarchy) are often specified a priori; automatic discovery and selection mechanisms remain an active area.
• Interpretability: The increased flexibility of hierarchical structures may impair interpretability (e.g., deep GPs) compared to tree-structured or shrinkage models with explicit sparsity or gating (Liu et al., 2023, Tang et al., 2023).

Active research is exploring richer kernel classes (e.g., spectral mixtures), online or sequential hierarchical model adaptation, expressive approximations (e.g., interdomain, Kronecker, semiseparable factorizations), and meta-learning of hierarchical priors for transfer learning or automated Bayesian optimization across heterogeneous task families (Fan et al., 2022).

---

Hierarchical Gaussian processes represent a unifying and highly adaptive probabilistic modeling framework that enables principled, flexible, and scalable uncertainty quantification, multi-task and multi-modal learning, and interpretable compositional structure across diverse domains of scientific machine learning. Key advances have resolved previous bottlenecks in depth, scalability, and cross-domain generalization, while motivating further work in inference algorithms, hierarchical kernel theory, and broader application to complex structured data.