Gaussian Hierarchies: Scalable High-Dimensional Inference

• Gaussian hierarchies are multilevel frameworks that decompose Gaussian random fields into coarse and fine components for scalable inference.
• They integrate Karhunen–Loève expansions at coarse levels with stochastic PDE sampling at fine levels to preserve full covariance and ergodicity.
• This approach significantly reduces computational cost while maintaining accurate uncertainty quantification in complex inverse problems.

A Gaussian hierarchy is a multilevel structure for representing, sampling, or inferring Gaussian random fields (GRFs) or Gaussian latent variables in high-dimensional Bayesian inference. Such hierarchical decompositions enable both scalable computational methods and effective dimensionality reduction, making them well-suited for uncertainty quantification and inverse problems in computational science and engineering. Recent advances have established hybrid approaches that integrate coarse-level Karhunen–Loève (KL) expansions with finer-level stochastic partial differential equation (SPDE) sampling, producing efficient and ergodic algorithms for sampling from GRFs in multilevel Markov Chain Monte Carlo (MLMCMC).

1. Hierarchical Decomposition of Gaussian Random Fields

A central concept in high-dimensional Gaussian hierarchies is the decomposition of a fine-level GRF sample as the sum of a projected coarse-level component and a fine-level complement:

$$\widetilde{\eta}_L = Q_L \eta_\ell^L + (I - Q_L)\eta_L$$

where:

• $\widetilde{\eta}_L$ is the full fine-level realization (e.g., on the finest mesh in a multigrid hierarchy).
• $Q_L$ is a projection operator mapping from the coarse space to the fine space.
• $\eta_\ell^L$ is the coarse-level random field, projected to the fine level.
• $\eta_L$ is an independent sample from the fine-level GRF.

This decomposition ensures that the resulting field $\widetilde{\eta}_L$ has the target distribution $\mathcal{N}(0, K_L)$, where $K_L$ is the fine-level covariance. The spectral equivalence assumption for covariances across levels is required to maintain the statistical properties of the hierarchy.

The hierarchical decomposition aligns closely with geometric multigrid hierarchies constructed from finite element spaces: $\Theta_0 \subset \Theta_1 \ldots \subset \Theta_L$, where prolongation, restriction, and projection operators enable multilevel transfer of information. At each level, samples are constructed recursively: the coarse-level sample is prolonged, and an independent complement from the orthogonal subspace is added, preserving Gaussianity and covariance.

2. Integration of Karhunen–Loève Expansions and SPDE Sampling

The hierarchical method integrates:

• Truncated Karhunen–Loève (KL) expansion at the coarsest level, providing a low-dimensional representation that captures the dominant modes of the GRF.
• SPDE-based sampling at finer levels, efficiently generating the remaining high-frequency structure by solving stochastic PDEs.

The hierarchical noise decomposition used for sampling at the finest level is:

$$\widetilde{W}_L = P M_\ell^{-1/2} V \xi_\ell + (I - P V^T \Pi) W_L$$

Here,

• $P$ is the prolongation operator,
• $M_\ell$ is the mass matrix at level $\ell$,
• $V$ is the matrix of KL eigenvectors (basis),
• $\xi_\ell$ are the KL coefficients, sampled from independent standard Gaussians,
• $W_L$ is fine-level white noise,
• $\Pi$ denotes the projection used for orthogonalization against the KL space.

The first term reconstructs the coarse KL modes at the fine level, while the second term (in the orthogonal complement) is sampled via SPDE, ensuring the overall field spans the full desired function space.

This methodology capitalizes on the computational efficiency of sparse-PDE solvers at high dimension and the low-rank efficiency of KL for dominant modes, sidestepping the pitfalls of pure KL truncation (which leads to non-ergodic, low-rank sampling) and expensive full SPDE solves at every level.

3. Dimensionality Reduction and Computational Efficiency

The hierarchical approach achieves dimensionality reduction and algorithmic scalability by:

• Restricting high-dimensional sampling (expensive in MCMC and data generation) to finer levels and operating only in the KL subspace at the coarse level.
• Conferring all statistical information for coarse (long-range) field structure to the low-dimensional KL modes.
• Sampling the fine-scale structure exclusively in the orthogonal subspace, with SPDE methods tailored for sparse, parallel solution of elliptic operators, and bypassing the high cost of high-rank eigendecomposition in KL across the entire domain.

Efficiency gains over conventional approaches are concretely demonstrated in numerical benchmarks. For instance, using 50 KL modes in the coarsest level, the total computational time in subsurface flow likelihood estimation is significantly reduced (e.g., 33,000 seconds for hybrid KL-SPDE vs. 186,000 seconds for pure SPDE). Integrated autocorrelation times (IACT) and effective sample sizes for posterior expectations are improved, substantiating the algorithmic and statistical superiority of the multilevel hierarchical design.

4. Multilevel Markov Chain Monte Carlo and Statistical Accuracy

In the context of MLMCMC, Gaussian hierarchies enable structured variances and covariance-consistent proposals at each level. The multilevel estimator for a quantity of interest (QoI) $Q$ is, for hierarchy levels $0,\ldots,L$:

$$\mathbb{E}_{\pi_L}[Q_L] = \mathbb{E}_{\pi_0}[Q_0] + \sum_{\ell=1}^L \left( \mathbb{E}_{\pi_{\ell}}[Q_{\ell}] - \mathbb{E}_{\pi_{\ell-1}}[Q_{\ell-1}] \right)$$

where $\pi_{\ell}$ is the posterior at level $\ell$. The variance and computational cost per effective sample are systematically reduced due to the hierarchical KL-SPDE design, offering superior convergence and mixing properties. Acceptance rates for fine-level proposals and variance decay across levels are empirically validated in subsurface flow inversion, with multilevel variance decaying monotonically and well-controlled across levels.

Posterior means and variances of the underlying GRF and associated quantities (e.g., pressure in Darcy law) are accurately estimated, matching ground truth with both visual and quantitative fidelity.

5. Representative Mathematical Formulations

Several key mathematical operations are central to Gaussian hierarchies in this context:

• KL Expansion of the Covariance Operator:

$$C(x,y) = \sum_{k=1}^\infty \lambda_k \psi_k(x) \psi_k(y)$$

$$\theta(x) = \sum_{i=1}^\infty \xi_i \sqrt{\lambda_i} \psi_i(x),\quad \xi_i \sim \mathcal{N}(0,1)$$

• SPDE Sampling (Whittle–Matérn Covariance Case):

$$\mathcal{A} \theta(x) = W(x),\quad \mathcal{A} = \frac{1}{g}(\Delta - \kappa^2)^{\alpha/2}$$

$\theta_h = -g A_h^{-1} M_h^{1/2} \xi_h$

• Hierarchical Decomposition for GRF Sampling:

$$\widetilde{\eta}_L = Q_L \eta^L_\ell + (I - Q_L)\eta_L$$

• Hybrid KL-SPDE Noise Decomposition:

$$\widetilde{W}_L = P M_\ell^{-1/2} V \boldsymbol{\xi} + (I - P V^T \Pi) W_L$$

• Variance and IACT:

$\tau = 1 + 2 \sum_{\chi=1}^{M} \rho_Q(\chi)$

where $\rho_Q(\chi)$ is the normalized autocorrelation function.

These formulations realize structure-preserving, scalable sampling algorithms, ensuring that at every hierarchical level, the marginals and conditionals remain consistent with the target GRF.

6. Empirical Performance and Superiority over Previous Methods

Benchmark studies on Bayesian inference for subsurface flow systems—characterized by high-dimensional GRFs and data-limited pressure observations—show that the hierarchical hybrid method outperforms pure KL and pure SPDE strategies. For a fixed estimator variance:

• Total cost (in computational time) is reduced, especially at moderate to large sample sizes.
• Variance decay and chain mixing properties, expressed through rapid reduction in variance of QoI differences and lower IACT, are improved.
• Structure-preserving ergodicity is maintained; the MCMC chain can stochastically explore the entire admissible parameter space.

These results are not contingent on ad hoc parameter tuning: efficiency is the product of the statistical and computational structure of the method.

7. Summary Table: Contrasting Approaches in Gaussian Hierarchies

Method	Dimensionality Reduction	Covariance Preservation	Ergodicity	Computational Cost (per sample)
Pure KL (truncated)	Yes (strong)	No (truncates modes)	No	High at high dim
SPDE-based	Partial (hierarchy)	Yes	Yes	High (repeated PDE solves)
Hybrid KL–SPDE (this approach)	Yes (coarse only)	Yes (full structure)	Yes	Lower overall

Summary and Outlook

The hierarchical integration of KL expansions and SPDE-based sampling for Gaussian random fields defines a robust computational framework for high-dimensional Bayesian UQ. By leveraging projections from coarse (KL) to fine (SPDE) levels within a consistent multigrid hierarchy, these Gaussian hierarchies achieve effective dimensionality reduction, full covariance and ergodicity preservation, and multilevel acceleration—all validated in realistic benchmark applications. This methodology offers a scalable, accurate, and theoretically well-founded solution to large-scale inverse problems in scientific computing and beyond.