📊 statistics

Multivariate Gaussian process emulation for multifidelity computer models with high-dimensional spatial outputs

This paper introduces a multivariate Gaussian process emulation framework for high-dimensional spatial outputs in multifidelity computer models, utilizing two novel cross-covariance structures based on autoregressive cokriging to enable computationally efficient and accurate storm surge risk assessment.

Original authors: Cyrus S. McCrimmon, Pulong Ma

Published 2026-03-31

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Cyrus S. McCrimmon, Pulong Ma

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict exactly how high the water will rise during a hurricane in a specific coastal town. You have two tools to help you:

The "Fast & Rough" Tool: A simplified computer model that runs in minutes. It's quick, but it ignores things like giant waves crashing against the shore, so it's not perfectly accurate.
The "Slow & Perfect" Tool: A super-detailed computer model that includes every wave, wind shift, and sandbar. It's incredibly accurate, but running it once takes days of supercomputer time.

To make good safety plans (like where to build sea walls), you need to run the "Slow & Perfect" tool thousands of times to test different storm scenarios. But you can't afford the time or money to do that.

The Solution: The "Smart Translator" (The Emulator)
This paper introduces a clever statistical trick called an emulator. Think of it as a "Smart Translator" that learns from the "Fast & Rough" tool to predict what the "Slow & Perfect" tool would have said, without actually running the slow one.

The authors, Cyrus McCrimmon and Pulong Ma, realized that predicting the water level for one house is easy, but predicting it for 9,000 houses at once is a nightmare. The water levels at neighboring houses aren't random; if the water is high at House A, it's likely high at House B next door. This is called spatial dependence.

The paper solves two big problems:

Speed: How do we handle 9,000 locations without the math taking forever?
Accuracy: How do we make sure our "uncertainty" (our guess of how wrong we might be) is correct when we look at the whole neighborhood, not just one house?

They propose two different "Smart Translators" to solve this:

1. The "Neighborhood Watch" Approach (SEP Model)

Imagine you are trying to guess the temperature of every house on a long street.

The Idea: You assume the temperature at House 10 depends mostly on House 9, which depends on House 8, and so on. You don't need to know the temperature of House 1 to guess House 10; you just need the immediate neighbors.
The Metaphor: This is like a domino effect. If you know the first few dominoes fall, you can predict the rest of the line without checking every single one.
How it works: The authors use a mathematical trick (called a "sparse Cholesky prior") that forces the model to only look at the "nearest neighbors" in the data. This makes the math super fast because it ignores the distant, irrelevant connections.
The Catch: It assumes the relationship between neighbors is the same everywhere (like a perfectly straight street). This works well, but sometimes the "street" is curvy or bumpy.

2. The "Orchestra Conductor" Approach (NONSEP Model)

Now imagine the storm surge isn't a straight line, but a complex symphony.

The Idea: Instead of looking at house-by-house, this model looks for the "main themes" of the storm. Maybe the first "theme" is the overall height of the water, and the second "theme" is how the water tilts from left to right.
The Metaphor: Think of the storm surge as a song. The "Fast & Rough" tool sings a simple version of the song. The "Slow & Perfect" tool sings a complex version. The authors use Principal Component Analysis (PCA) to break the complex song down into its main notes (the "basis functions").
How it works: The model learns how the "Fast" version of the song relates to the "Slow" version for just those main notes. It's much easier to translate a 5-note melody than a 9,000-note symphony. Once it translates the notes, it rebuilds the full 9,000-note song.
The Benefit: This is more flexible. It can handle "bumpy streets" and complex patterns where neighbors don't just depend on their immediate next-door neighbor.

Why Does This Matter?

The paper tested these two methods on a real-world problem: predicting storm surges in Cape Coral, Florida.

The Result: Both methods were fast and accurate. However, the "Orchestra Conductor" (NONSEP) approach was slightly better at capturing the complex, real-world patterns of the storm.
The Big Win: The most important finding was about Risk Assessment. If you ignore how houses are connected (spatial dependence), your safety estimates are wrong.
- If you treat 9,000 houses as independent, you might think the risk is low.
- But because the water rises together, the total risk to the whole town is actually much higher.
- These new models correctly calculate that "group risk," helping FEMA and engineers build better defenses.

In a Nutshell

The authors built a mathematical shortcut that lets us use cheap, fast computer models to predict the behavior of expensive, slow ones. They did this by realizing that storm surges are connected (like dominoes or a song) and creating two different ways to map those connections. This allows scientists to run thousands of "what-if" storm scenarios quickly, giving us better data to protect coastal cities from hurricanes.

1. Problem Statement

The paper addresses the challenge of Uncertainty Quantification (UQ) for deterministic computer models that produce high-dimensional spatial outputs.

Context: The motivating application is hurricane-driven storm surge modeling (using the ADCIRC and ADCIRC+SWAN simulators). These models simulate physical processes over large spatial domains (e.g., thousands of grid points) and are computationally expensive to run at high fidelity.
Multifidelity Setting: Simulators exist at varying levels of fidelity (accuracy vs. cost). Low-fidelity models are fast but less accurate; high-fidelity models are accurate but prohibitively expensive for extensive UQ.
The Core Challenge:
1. High Dimensionality: Outputs are vectors of size $N$ (e.g., $N \approx 9,000$ spatial locations), making standard multivariate Gaussian Process (GP) emulation computationally infeasible due to the $O(N^3)$ cost of inverting covariance matrices.
2. Spatial Dependence: Outputs at nearby spatial locations are correlated. Ignoring this cross-covariance leads to incorrect uncertainty quantification, particularly when aggregating predictions over a region (e.g., calculating average surge for a flood zone).
3. Multifidelity Integration: The emulator must leverage abundant low-fidelity data to predict high-fidelity outputs efficiently while maintaining the natural ordering of accuracy.

2. Methodology

The authors propose a Bayesian Multivariate Gaussian Process framework based on the Autoregressive Cokriging model (Kennedy & O'Hagan, 2000). They develop two distinct cross-covariance structures to handle the high-dimensional output space: SEP (Separable) and NONSEP (Nonseparable).

2.1 General Framework

Let $y_t(x)$ be the output vector at fidelity level $t$ and input $x$ . The model assumes an autoregressive structure:
$y_t(x) = \gamma_{t-1} y_{t-1}(x) + \delta_t(x)$
where $\gamma_{t-1}$ is a scale discrepancy parameter, and $\delta_t(x)$ is a residual process modeled as a Matrix-Variate Gaussian Process (MVGP).

2.2 SEP ARCokrig (Separable Approach)

Structure: Assumes the covariance between inputs and outputs is separable: $\text{Cov}(y) = R \otimes \Sigma$ , where $R$ is the input correlation matrix and $\Sigma$ is the $N \times N$ output cross-covariance matrix.
Innovation: To handle high $N$ , the authors do not estimate $\Sigma$ directly. Instead, they model the precision matrix $\Omega = \Sigma^{-1}$ using a sparse Cholesky decomposition: $\Omega = (I - A)^\top D^{-1} (I - A)$ .
Sparsity Induction: By ordering spatial locations (e.g., via maximum-minimum distance), they assume that the $j$ -th output depends only on a small set of its $p$ nearest predecessors ( $j-1, \dots, j-p$ ). This induces sparsity in matrix $A$ .
Inference: This transforms the estimation of the cross-covariance into a series of $N$ sparse autoregressive regressions. Priors are placed on the Cholesky factors (sparse normal and inverse-gamma), allowing for computationally efficient Bayesian estimation even when $N \gg n$ (number of model runs).

2.3 NONSEP ARCokrig (Nonseparable Approach)

Structure: Based on the Linear Model of Coregionalization (LMC) and basis representation (Higdon et al., 2008).
Basis Representation: Outputs are approximated via Principal Component Analysis (PCA):
$y_t(x) \approx \sum_{\ell=1}^{p_t} k_{t,\ell} w_{t,\ell}(x)$
where $k_{t,\ell}$ are basis vectors (principal components) and $w_{t,\ell}(x)$ are random weights.
Autoregressive Weights: The random weights $w_{t,\ell}(x)$ are modeled using an autoregressive cokriging structure across fidelity levels. This preserves the rank ordering of accuracy while allowing nonseparable dependence between input and output spaces.
Inference: The model fits independent autoregressive GPs for each principal component weight. This reduces the dimensionality of the problem significantly.

2.3 Posterior Inference

Both models allow for analytic integration of regression coefficients ( $\beta$ ) and variance parameters ( $\sigma^2$ ) due to conjugacy.
Posterior sampling is required only for the input correlation parameters ( $\theta$ ).
Crucially, the autoregressive structure allows parameters at each fidelity level to be inferred separately, avoiding the need to sample a massive joint posterior, thus ensuring computational scalability.

3. Key Contributions

Scalable Multivariate Emulation: The paper bridges the gap between multifidelity emulation and high-dimensional spatial outputs by introducing sparse Cholesky priors (SEP) and PCA-based dimension reduction (NONSEP) within an autoregressive framework.
Proper Uncertainty Quantification (UQ): The authors demonstrate that ignoring cross-covariance (as in independent models) leads to severe under-coverage or over-coverage when aggregating predictions over spatial regions. Their methods correctly quantify uncertainty for both individual locations and aggregated spatial resolutions.
Computational Efficiency: By exploiting conditional independence and sparse structures, the methods avoid the $O(N^3)$ bottleneck, making them feasible for outputs with thousands of spatial locations.
Unified Framework: The approach generalizes the standard univariate autoregressive cokriging to the multivariate setting, accommodating both separable and nonseparable dependence structures.

4. Results

The methods were evaluated on two datasets: a synthetic chemical spill testbed and the real-world storm surge application.

Testbed Example (Chemical Spill):
- Metrics: Root Mean Square Prediction Error (RMSPE), Coverage Probability (CVG), and Average Length of Credible Intervals (ALCI).
- Findings: The SEP model achieved the lowest RMSPE. Both SEP and NONSEP provided well-calibrated uncertainty (CVG $\approx$ 95%) for aggregated outputs, whereas models ignoring cross-covariance (PP, PP-ARCokrig) failed significantly in aggregated UQ (CVG dropped to ~20-23%).
Storm Surge Application (Cape Coral, FL):
- Data: 9,284 spatial locations, 2 fidelity levels (ADCIRC vs. ADCIRC+SWAN).
- Findings:
  - NONSEP outperformed SEP in terms of UQ (CVG and ALCI) for both individual and aggregated outputs, likely due to the complex, non-stationary nature of storm surge data.
  - SEP had slightly better marginal RMSPE but exhibited slight under-coverage in aggregated regions.
  - Models ignoring cross-covariance (PP, PP-ARCokrig) showed catastrophic failure in aggregated UQ (CVG < 5%), confirming the necessity of modeling spatial dependence.
  - The NONSEP model with 7 principal components captured >99% of variance and reduced prediction error by >50% compared to a standard kriging model (using only high-fidelity data).

5. Significance

Practical Impact: The proposed emulators enable rapid and accurate risk assessment for coastal flood hazards. Decision-makers can now obtain reliable uncertainty estimates for aggregated metrics (e.g., total flood risk for a county) without running thousands of expensive high-fidelity simulations.
Theoretical Advancement: The paper provides a robust solution to the "curse of dimensionality" in multivariate GP emulation. It demonstrates that sparse precision matrices and basis representations are viable strategies for high-dimensional outputs in multifidelity settings.
Future Directions: The framework is extensible to more than two fidelity levels and can be adapted for outputs at different spatial resolutions (change of support), offering a versatile tool for complex physical simulations.

In conclusion, McCrimmon and Ma successfully develop a computationally feasible, statistically rigorous framework for emulating high-dimensional multifidelity computer models, ensuring that uncertainty quantification remains valid even when predictions are aggregated over large spatial domains.