# Multi-Fidelity Residual Neural Processes for Scalable Surrogate Modeling

Ruijia Niu 1
& Dongxia Wu 1 2
& Kai Kim 1
& Yi-An Ma 2 1
& Duncan Watson-Parris 2 3
& Rose Yu 1 2

###### Abstract

  1. Multi-fidelity surrogate modeling aims to learn an accurate surrogate at the highest fidelity level by combining data from multiple sources. Traditional methods relying on Gaussian processes can hardly scale to high-dimensional data. Deep learning approaches utilize neural network based encoders and decoders to improve scalability. These approaches share encoded representations across fidelities without including corresponding decoder parameters. At the highest fidelity, the representations are decoded with different parameters, making the shared information inherently inaccurate. This hinders inference performance, especially in out-of-distribution scenarios when the highest fidelity data has limited domain coverage. To address these limitations, we propose Multi-fidelity Residual Neural Processes (MFRNP), a novel multi-fidelity surrogate modeling framework. MFRNP optimizes lower fidelity decoders for accurate information sharing by aggregating lower fidelity surrogate outputs and models residual between the aggregation and ground truth on the highest fidelity. We show that MFRNP significantly outperforms current state-of-the-art in learning partial differential equations and a real-world climate modeling task.

11footnotetext: Department of Computer Science and Engineering, University of
California San Diego, La Jolla, California, USA22footnotetext: Halıcıoğlu Data
Science Institute, University of California San Diego, La Jolla, California,
USA33footnotetext: Scripps Institution of Oceanography, University of
California San Diego, La Jolla, California, USA††footnotetext: Preprint. Under
Review. Correspondence to: Rose Yu<EMAIL_ADDRESS>
##  1 Introduction

From engineering to climate science, a computational model, often realized by
simulation, is frequently used to characterize the input-output relationship
of a physical system. These computational models can be simulated at multiple
levels of sophistication. The high-fidelity simulators are more accurate but
resource demanding, whereas lower-fidelity models are less accurate but more
computationally efficient. In many applications, calibrating results from
computational models with real-world observations plays a crucial role
Hersbach et al. (2020); Karami and Kashef (2020). While the calibrated data
has refined details and accuracy, the domain coverage is very limited. Multi-
fidelity surrogate modeling Peherstorfer et al. (2018) aims to balance the
computation efficiency-accuracy trade-off by utilizing data across fidelities
to learn an accurate surrogate at the highest fidelity.

Gaussian Processes (GPs) Seeger (2004) are a popular tool for surrogate
modeling. Recent works have attempted to extend GP to multi-fidelity setting
(Perdikaris et al., 2016; Wang et al., 2021). However, they inherit the
limited scalability from GP due to the inversion of the kernel matrix
(Williams and Rasmussen, 1995; Rasmussen, 2003). To solve this issue, many
have proposed deep learning-based surrogate models (Damianou and Lawrence,
2013; Raissi and Karniadakis, 2016; Salimbeni and Deisenroth, 2017; Wilson et
al., 2016). Neural Processes (NPs) (Wang and Lin, 2020; Hebbal et al., 2021)
stand out as one of the most appealing approaches regarding inference
performance and scalability. With neural networks as encoders and decoders,
NPs are capable of encoding fidelity-specific data into low-dimensional latent
representations and use them to improve inference performance at the highest
fidelity, bypassing the scalability issue from high-dimensional data.

For example, Deep Multi-fidelity Active Learning (DMFAL) Li et al. (2022a, b)
proposed to pass information from lower fidelities to higher fidelities with
hidden representations via neural network. This method requires a hierarchical
structure passing information from low to high fidelity level, which can lead
to the error propagation issue. Disentangled Multi-fidelity Deep Bayesian
Active Learning (D-MFDAL) Wu et al. (2023) alleviates this issue by
redesigning the NP using local and global representations. Nevertheless, these
methods rely on latent representations from encoders for cross-fidelity
information sharing. However, the decoder parameters varies across fidelities
and are not shared. At the highest fidelity, shared representations are
decoded with different parameters, making the decoded information inherently
inaccurate.

This limits the inference performance, especially when the model needs to
extrapolate with shared information given input that is out-of-distribution
(OOD) w.r.t the highest fidelity training data.

In this work, we introduce a novel multi-fidelity surrogate modeling
framework, Multi-fidelity Residual Neural Process (MFRNP), to address the
aforementioned issues. MFRNP aggregates the predictions from surrogate models
across lower-fidelity levels and employs an NP surrogate to capture the
residual between the aggregated prediction and the ground truth at the highest
fidelity level. By directly utilizing the outputs from lower-fidelity
surrogates to share information, MFRNP includes decoders in cross-fidelity
information sharing, improving accuracy while maintaining scalability.
Moreover, we developed a tailored Evidence Lower Bound, named Residual-ELBO,
to serve as our loss function. This novel loss function ensures the highest
fidelity latent variable $z_{K}$ depends on all the other latent variables and
decoder parameters across fidelities. Thus, MFRNP ensures accurate information
aggregation from lower fidelities to promote residual modeling at the highest
fidelity.

To summarize, our contributions include:

  * •

A novel multi-fidelity surrogate model, Multi-fidelity Residual Neural Process
(MFRNP). Its architecture shares input-specific information from lower
fidelities, tackles the varying decoder problem with no error propagation yet
preserving scalability.

  * •

A novel Residual-ELBO to simultaneously promote learning across fidelities and
optimize lower fidelity decoders for residual modeling at the highest
fidelity.

  * •

Superior performance in global scale real-world climate modeling and numerous
benchmark tasks on partial differential equations. MFRNP outperform the state-
of-the-art baseline by $\sim$90% in average.

##  2 Background

#### Multi-Fidelity Modeling.

Multi-fidelity modeling aims to capture the complex mapping from low-
dimensional input variables $\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$ to high-
dimensional output $\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$ of the function
$f:\mathcal{X}\rightarrow\mathcal{Y}$. In $K$ fidelity modeling where $K>1$,
the cost $c_{k}$ of evaluating $f_{k}\in\\{f_{1},\cdots,f_{K}\\}$ increases
with the fidelity level $(c_{1}<\cdots<c_{K})$ as $f_{k}$ conveying more
detailed information in approximating $f$. Our goal is to learn a deep
surrogate model $\hat{f}_{K}$ to approximate $f_{K}$ by combining data from
$K$ fidelities and $N$ parameter settings
$\\{x_{k,n},y_{k,n}\\}_{k=1,n=1}^{K,N}$.

#### Neural Processes.

Combining Gaussian Processes (GPs) and neural networks (NNs), Neural Processes
(NPs) (Garnelo et al., 2018a) constitute a family of latent variable models
for implicit stochastic processes (Wang and Van Hoof, 2020). NPs represent
distributions over functions and estimate prediction uncertainties like GPs
while featuring scalability in high dimensions (Jha et al., 2022).

Formally, NP consists of latent variables $z\in\mathbb{R}^{d_{z}}$ and model
parameters $\theta$, trained on context set
$\mathcal{D}^{c}\equiv\\{x^{c}_{n},y^{c}_{n}\\}_{n=1}^{N}$ and target sets
$\mathcal{D}^{t}\equiv\\{x^{t}{m},{y^{t}_{m}}\\}_{m=1}^{M}$. $\mathcal{D}^{c}$
and $\mathcal{D}^{t}$ are randomly split from the training set $\mathcal{D}$.
Learning the posterior of $z$ and $\theta$ equals to maximize the posterior
likelihood below:

$p(y^{t}_{1:M}|x^{t}_{1:M},\mathcal{D}^{c},\theta)=\int
p(z|\mathcal{D}^{c},\theta)\prod^{M}_{m=1}p(y^{t}_{m}|z,x^{t}_{m},\theta)dz$
(1)

Due to the intractability of marginalizing over latent variables $z$, the NP
family utilize approximate inference, deriving the corresponding evidence
lower bound (ELBO):

$\displaystyle\log
p(y^{t}_{1:M}|x^{t}_{1:M},\mathcal{D}^{c},\theta)\geq\mathbb{E}_{q_{\phi}(z|\mathcal{D}^{c}\cup\mathcal{D}^{t})}\big{[}\sum_{m=1}^{M}\log
p(y^{t}_{m}|z,x^{t}_{m},\theta)+\log\frac{q_{\phi}(z|\mathcal{D}^{c})}{q_{\phi}(z|\mathcal{D}^{c}\cup\mathcal{D}^{t})}\big{]}$
(2)

This variational approach approximates the intractable true posterior
$p(z|\mathcal{D}^{c},\theta)$ with the approximate posterior
$q_{\phi}(z|\mathcal{D}^{c})$. Here, $\phi$ parameterize the encoder and
$\theta$ parameterize the decoder of the model.

##  3 Methodology

Figure 1:  Graphical model comparison of MFRNP against the state-of-the-art
D-MFD baseline model. Right: D-MFD disentangles the latent representations
$r_{k,n}$ into local and global representations $L_{k,n}$ and $G_{k,n}$ to
infer $z_{k}$. Each $G_{k,n}$ is from a different fidelity level, while no
information about the decoder parameters $\theta_{k}$ is shared, making
$G_{k,n}$ inaccurate regarding $\theta_{i}$ where $i\neq k$. Left: MFRNP
dynamically constructs $\mathcal{D}^{train}=\\{x_{K},R(x_{K})\\}$ in cross-
fidelity optimization step and learn fidelity-specific information with in-
fidelity learning step. The residual function $R$ makes $z_{K}$ dependant on
$z_{1:K-1}$ and $\theta_{1:K-1}$ without inflating any dimensions. Thus, MFRNP
optimizes lower fidelity decoders for better information sharing. We use bold
letters in our model graph to denote a set of variables.

The high level goal of MFRNP is to ensure the output aggregation from lower
fidelities is informative for modeling the residual between the aggregated
prediction and the ground truth at the highest fidelity level.

For each fidelity $k\in\\{1...K-1\\}$, MFRNP learns a NP surrogate
$\hat{f}_{k}$ to approximate $f_{k}$. For the highest fidelity $K$, MFRNP
learns a NP to approximate the residual function
$R(x)=f_{K}(x)-[\frac{1}{K-1}\sum_{m=1}^{K-1}\hat{f}_{m}(x)]$. The lower
fidelity predictions are aggregated by linearly interpolating to the
resolution at the highest fidelity and taking the average. We choose to
average over other aggregation methods because this stabilizes the aggregation
and ensures the equal contribution of lower fidelity surrogates. Unlike
previous work where all the information from a fidelity is shared via the
latent variable during encoding, MFRNP shares information via the decoded
outputs. This explicitly involves the decoder at each fidelity to improve
information-sharing. The input-specific information-sharing scheme further
tackles the error propagation issue and ensures the accuracy and full
exploitation of lower-fidelity information.

Approximating the Residual Function. The core idea of MFRNP is to make the
highest fidelity latent variable $z_{K}$ dependent on all other latent
variables as well as the decoder parameters at lower fidelities. We introduce
the dependency by approximating $R(x)$ at the highest fidelity, which models
the residual between ground truth and the proposed aggregation from lower
fidelities given input $x$.

Intuitively, $\hat{f}_{k\in\\{1,...,K-1\\}}$ is optimized with two objectives.
As shown in Figure 1, the in-fidelity learning step captures fidelity-specific
information while the cross-fidelity optimization step encourages lower
fidelity surrogates to propose accurate and informative aggregations w.r.t the
highest fidelity given input $x_{K}$.

During inference, given target input $x^{t}$, MFRNP approximates
$f_{K}(x^{t})$ by aggregating predictions from fidelities $1...K-1$ and adding
it with the correction term from fidelity
$K:\hat{f}_{K}(x^{t})=R(x^{t})+[\frac{1}{K-1}\sum_{m=1}^{K-1}\hat{f}_{m}(x^{t})]$.

MFRNP generates $\hat{y}^{t}$ by fully exploiting the rich information in
$\hat{y}_{1\cdots K-1}$. In the OOD setup where $\mathcal{D}^{train}_{K}$ and
$\mathcal{D}^{test}_{K}$ covers different input domains, MFRNP can effectively
explore the regions beyond the $\mathcal{D}^{train}_{K}$, thus enhancing its
input-domain extrapolation capabilities in modeling $f_{K}$.

Residual-ELBO. We design a Residual-ELBO (R-ELBO) for MFRNP. For each fidelity
$k\in\\{1\cdots K-1\\}$, we infer the latent variable $z_{k}$ with the NN
encoder $q_{\phi_{k}}(z_{k}|\mathcal{D}^{c})$ and decoder
$p_{\theta_{k}}(y^{t}_{k}|z_{k},x^{t}_{k})$. At fidelity $K$,
$\mathcal{D}_{K}$ is dynamically constructed with $\\{X_{K},R(X_{K})\\}$. Here
$X_{K}$ denotes a set of input variables $x_{K}$. We infer the latent variable
$z_{K}$ with the NN encoder
$q_{\phi_{K}}(z_{K}|z_{1..K-1},\theta_{1..K-1},\mathcal{D}_{K}^{c})$ and
decoder $p_{\theta_{k}}(R(X_{K})|z_{K},X_{K},\theta_{K})$, where $R(X_{K})$
depends on the proposed aggregations of lower fidelity predictions via
ancestral sampling Wang and Van Hoof (2020). We derive the R-ELBO for $K>1$
fidelities in two terms:

$\displaystyle\log
p(R(x_{K}^{t})|x^{t}_{K},\mathcal{D}_{1:K}^{c}\cup\mathcal{D}_{1:K-1}^{t},\theta_{K})$
$\displaystyle\geq\mathbb{E}_{q_{\phi_{K}}(z_{K}|\mathcal{D}_{1:K}^{c}\cup\mathcal{D}_{1:K}^{t})}\big{[}\log
p(R(x_{K}^{t})|z_{K},x^{t}_{K},\theta_{K})+\log\frac{q_{\phi_{K}}(z_{K}|\mathcal{D}_{1:K}^{c}\cup\mathcal{D}_{1:K-1}^{t})}{q_{\phi_{K}}(z_{K}|\mathcal{D}_{1:K}^{c}\cup\mathcal{D}_{1:K}^{t})}\big{]}$
$\displaystyle=\mathbb{E}_{q_{\phi_{K}}(z_{K}|z_{1:K-1},\theta_{1:K-1},\mathcal{D}_{K}^{c}\cup\mathcal{D}_{K}^{t}))}\big{[}\log
p(R(x_{K}^{t})|z_{K},x^{t}_{K},\theta_{K})+\log\frac{q_{\phi_{K}}(z_{K}|z_{1:K-1},\theta_{1:K-1},\mathcal{D}_{K}^{c})}{q_{\phi_{K}}(z_{K}|z_{1:K-1},\theta_{1:K-1},\mathcal{D}_{K}^{c}\cup\mathcal{D}_{K}^{t}))}\big{]}$
(3)

$\displaystyle\log
p(y^{t}_{1:K-1}|x^{t}_{1:K-1},\mathcal{D}_{1:K-1}^{c},\theta_{1:K-1})$
$\displaystyle\geq\mathbb{E}_{q_{\phi}(z_{1:K-1}|\mathcal{D}_{1:K-1}^{c}\cup\mathcal{D}_{1:K-1}^{t})}\big{[}\log
p(y^{t}_{1:K-1}|z_{1:K-1},x^{t}_{1:K-1},\theta_{1:K-1})+..+\log\frac{q_{\phi}(z_{1:K-1}|\mathcal{D}_{1:K-1}^{c})}{q_{\phi}(z_{1:K-1}|\mathcal{D}_{1:K-1}^{c}\cup\mathcal{D}_{1:K-1}^{t})}\big{]}$
$\displaystyle=\mathbb{E}_{q_{\phi_{1}}(z_{1}|\mathcal{D}_{1}^{c}\cup\mathcal{D}_{1}^{t})..q_{\phi_{K-1}}(z_{K-1}|\mathcal{D}_{K-1}^{c}\cup\mathcal{D}_{K-1}^{t})}\big{[}\log
p(y^{t}_{1}|z_{1},x_{1}^{t},\theta_{1})+..+\log
p(y_{K-1}^{t}|z_{K-1},x_{K-1}^{t},\theta_{K-1})$
$\displaystyle+\log\frac{q_{\phi_{1}}(z_{1}|\mathcal{D}_{1}^{c})}{q_{\phi_{1}}(z_{1}|\mathcal{D}_{1}^{c}\cup\mathcal{D}_{1}^{t})}+..+\frac{q_{\phi_{K-1}}(z_{K-1}|\mathcal{D}_{K-1}^{c})}{q_{\phi_{K-1}}(z_{K-1}|\mathcal{D}_{K-1}^{c}\cup\mathcal{D}_{K-1}^{t})}\big{]}$
(4)

Equation 4 is a unified ELBO accounting for learning from fidelity specific
datasets at lower fidelities. Equation 3 is the ELBO for the residual function
at the highest fidelity. This term introduces dependency to every fidelity and
optimizes the output aggregations to better approximate $R(x)$.

For training, we calculate the R-ELBO with Monte Carlo (MC) sampling and
ancestral sampling (AS) to optimize the objective function below:

$\displaystyle\mathcal{L}^{\hat{f}}_{MC}$
$\displaystyle=\sum^{K-1}_{k=1}\bigg{[}\frac{1}{S}\sum^{S}_{s=1}\log
p(y_{k}^{t}|x_{k}^{t},z_{k}^{(s)})-\text{KL}[q(z_{k}|\mathcal{D}_{k}^{c},\mathcal{D}_{k}^{t})\|p(z_{k}|\mathcal{D}_{k}^{c})]\bigg{]}$
(5) $\displaystyle\mathcal{L}^{R}_{MC}$
$\displaystyle=\frac{1}{S}\sum^{S}_{s=1}\log
p(R(x_{K}^{t})|x_{K}^{t},z_{K}^{(s)})-\text{KL}[q(z_{K}|\mathcal{D}_{K}^{c},\mathcal{D}_{K}^{t})\|p(z_{K}|\mathcal{D}_{K}^{c})]$
(6) $\displaystyle\mathcal{L}_{MC}$
$\displaystyle=\mathcal{L}^{R}_{MC}+\mathcal{L}^{\hat{f}}_{MC}$ (7)

where the time for $q_{\phi}(z|\mathcal{D}^{c})$ to sample $z_{k}^{(s)}$
scales linearly with the number of fidelity levels.

Training & Inference. As shown in Algorithm 1, MFRNP calculates
$\mathcal{L}^{\hat{f}}_{MC}$ by predicting the set $\hat{Y}_{k}^{t}$ given
input target set $X_{k}^{t}$ for every fidelity $k\in\\{1..K-1\\}$. For
$\mathcal{L}^{R}_{MC}$, the AS steps are highlighted in orange in Algorithm 1.

MFRNP introduce decoder to the information aggregation step by predicting with
the highest fidelity level input set $X_{K}$ at lower fidelity levels. We note
the prediction as set $\\{\hat{o}_{k,K}\\}$ where $k\in\\{1..K-1\\}$,
representing surrogate predictions given $X_{K}$ at fidelity $k$. We linearly
interpolate each $\\{\hat{o}_{k,K}\\}$ to math the dimension of $Y_{K}$.

Then we dynamically construct the training dataset at $K$ as
$\mathcal{D^{\prime}}_{K}=\\{X_{K},R(X_{K})\\}$. Finally, we randomly split
$\mathcal{D^{\prime}}_{K}$ into
$\mathcal{D^{\prime}}_{K}^{c},\mathcal{D^{\prime}}_{K}^{t}$ to obtain
$\mathcal{L}^{R}_{MC}$ and perform back propagation with $\mathcal{L}_{MC}$.
Here, $\mathcal{L}^{R}_{MC}$ encourages MFRNP to optimize lower fidelity
surrogates for residual modeling at the highest fidelity,

while $\mathcal{L}^{\hat{f}}_{MC}$ regulates lower fidelity levels to learn
from fidelity-specific datasets.

The inference process is demonstrated in Algorithm 2 where the AS steps are
highlighted in orange. Given input set $\\{x_{i,K}\\}$, MFRNP propose
aggregated predictions $\\{a_{i,K}\\}$ from lower fidelities and the predicted
residual $\hat{r}_{i}$ at fidelity $K$. The final prediction is given as
$\\{\hat{y}_{i,K}\\}=\\{a_{i,K}\\}+\\{\hat{r}_{i}\\}$.

Input: Dataset $\mathcal{D}_{1\cdots K}$, number of fidelities $K$.

for $k=1$ to $K-1$ do

Randomly split $\mathcal{D}_{k}$ into
$\\{\mathcal{D}_{k}^{c},\mathcal{D}_{k}^{t}\\}$

Sample $\\{z_{i,k}\\}$ where $z_{i,k}\sim q_{k}(.|\mathcal{D}_{k}^{c})$

Predict $\\{\hat{y}_{i,k}^{t}\\}$ where $\hat{y}_{i,k}^{t}\sim
p_{k}(.|z_{i,k},x^{t}_{i,k})$

Sample $\\{z^{\prime}_{i,k}\\}$ where $z^{\prime}_{i,k}\sim
q_{k}(.|\mathcal{D}_{k}^{c})$

Predict $\\{\hat{o}_{i,k}\\}$ where $\hat{o}_{i,k}\sim
p_{k}(.|z^{\prime}_{i,k},x_{i,K})$

Linearly interpolate $\\{\hat{o}_{i,k}\\}$ to the resolution at $K$.

end for

Get the residual set $\\{r_{i}\\}$ at K where
$r_{i}=y_{i,K}-\frac{\sum(\hat{o}_{i,1},..,\hat{o}_{i,K-1})}{K-1}$

Random-split $\mathcal{D^{\prime}}_{K}=\\{x_{i,K},r_{i}\\}$ into
$\\{\mathcal{D^{\prime}}_{K}^{c},\mathcal{D^{\prime}}_{K}^{t}\\}$

Sample $\\{z_{i,K}\\}$ where $z_{i,K}\sim q_{K}(.|\mathcal{D}_{K}^{c})$

Predict $\\{\hat{r}^{t}_{i}\\}$ where $\hat{r}^{t}_{i}\sim
p_{K}(.|z^{t}_{i,K},x^{t}_{i,K})$

Back propagate with
$\mathcal{L}_{MC}=\mathcal{L}^{R}_{MC}+\mathcal{L}^{\hat{f}}_{MC}$ in Eqn
LABEL:eqn:L_ii

Algorithm 1 MFRNP Training Process

Input: Latent variables $z_{1..K}$, input $\\{x_{i,K}\\}$

for $k=1$ to $K-1$ do

Sample $\\{z^{\prime}_{i,k}\\}$ from $z_{k}$.

Predict $\\{\hat{o}_{i,k}\\}$ where $\hat{o}_{i,k}\sim
p_{k}(.|z^{\prime}_{i,k},x_{i,K})$

Linearly interpolate $\\{\hat{o}_{i,k}\\}$ to the resolution at $K$.

end for

Obtain aggregation
$\\{a_{i,K}\\}=\frac{\sum(\hat{o}_{i,1},..,\hat{o}_{i,K-1})}{K-1}$

Sample $\\{z^{\prime}_{i,K}\\}$ from $z_{K}$.

Predict $\\{\hat{r}_{i}\\}$ where $\hat{r}_{i}\sim p_{K}(.|z_{i,K},x_{i,K})$

Return $\\{\hat{y}_{i,K}\\}=\\{a_{i,K}\\}+\\{\hat{r}_{i}\\}$

Algorithm 2 MFRNP Inference Process

##  4 Related Work

#### Multi-fidelity Modeling.

Multi-fidelity surrogate modeling is prevalent across scientific and
engineering domains, including applications in aerospace systems Brevault et
al. (2020) and climate science Hosking (2020); Valero et al. (2021) Valero et
al. (2021). The foundational work of Kennedy and O’Hagan (2000) employs GPs to
connect models of varying fidelity levels, introducing an autoregressive
model. Le Gratiet and Garnier (2014) introduces a recursive GP with a nested
structure in the input domain to facilitate rapid inference. Perdikaris et al.
(2015, 2016) addresses high-dimensional GP scenarios by leveraging the Fourier
transformation of the kernel function. (Perdikaris et al., 2017) puts forth
the concept of multi-fidelity Gaussian processes (NARGP), assuming a nested
structure in the input domain for sequential training at each fidelity level.

Wang et al. (2021) presents a Multi-Fidelity High-Order GP model for
accelerating physical simulations. They extend the Linear Model of
Coregionalization (LMC) to the nonlinear case, incorporating a matrix GP prior
on the weight functions. Deep Gaussian processes (DGPs) (Cutajar et al., 2019)
formulate a unified objective to optimize kernel parameters jointly across
fidelity levels. DGPs face scalability challenges with high-dimensional data.
Infinite-Fidelity Coregionalization (IFC) from Li et al. (2022c) models the
output space as a continuous function of fidelity and input based on neural
ODEs, allowing the model to extrapolate to higher fidelities. However, IFC
faces scalability issue due to the computation-demanding ODE solver.

Multi-fidelity modeling has been integrated with deep learning. Guo et al.
(2022) employs deep neural networks to merge parameter-dependent output
quantities. Meng and Karniadakis (2020) proposes a composite neural network
for multi-fidelity data in inverse PDE problems, while Meng et al. (2021)
introduces Bayesian neural nets for multi-fidelity modeling. De et al. (2020)
employs transfer learning to fine-tune high-fidelity surrogate models using
deep neural networks trained on low-fidelity data. Despite advancements,
existing deep GP models (Cutajar et al., 2019; Hebbal et al., 2021) struggle
with cases where different fidelities involve data with varying dimensions.
Additionally, multi-fidelity methods have been investigated in Bayesian
optimization, active learning, and bandit problems (Li et al., 2020, 2022b;
Perry et al., 2019; Kandasamy et al., 2017).

#### Neural Processes.

Neural Processes (NPs) Garnelo et al. (2018b); Kim et al. (2018); Louizos et
al. (2019); Singh et al. (2019) emerge as scalable and expressive alternatives
to GPs for modeling stochastic processes. Previous work by Raissi and
Karniadakis (2016) combines multi-fidelity GPs with deep learning, placing a
GP prior on features learned by deep neural networks. Wang and Lin (2020)
proposes a multi-fidelity neural process with physics constraints (MFPC-Net),
leveraging NPs to capture correlations between multi-fidelity data.
Nonetheless, this model requires paired data and cannot utilize unpaired data
at the low-fidelity level.

The recent work of Wu et al. (2023) proposes to disentangle the latent
variable at each fidelity into local and global representations and share the
global part across fidelities. However, the fidelity-specific decoder
parameters are not included in information sharing. Thus, the highest fidelity
decoder expresses shared representations differently, hindering the inference
performance, especially in OOD scenarios.

#### Climate Modeling.

Climate modeling is a central component of modern climate science and the
primary tool for predicting future climate states Flato et al. (2014). Various
modeling centers around the world have developed distinct climate models. The
semi-independent development process has led to many plausible, but
disagreeing, climate models representing the same earth system Knutti et al.
(2010); Flato et al. (2014). Averaging these models often leads to improved
results compared to using individual models Lambert and Boer (2001); Gleckler
et al. (2008); Knutti et al. (2010). Proper aggregation of different climate
models for a consensus estimate is therefore an important topic Tebaldi and
Knutti (2007). Averaging all models with equal or varying weights has been the
most common approach Giorgi and Mearns (2002); Abramowitz et al. (2019), known
as (weighted) ensemble averaging (EA). However, EA techniques often do not
retain much spatial information and can cause severe blurring, corrupting
regional signals.

To address this problem, various alternative approaches have been developed
based on Bayesian hierarchical models, regression and machine learning that
all use observational data to improve model aggregation. DeepESD Baño-Medina
et al. (2022) utilize CNNs to learn the mapping from the EA of $8$ regional
climate models to the observation-calibrated data Dee et al. (2011) at finer
resolution. Similarly, NNGPR Harris et al. (2023) utilizes deep kernel GPs to
model the residual between the EA of $16$ low-resolution global climate models
and the up-to-date observation-calibrated ERA5-reanalysis dataset Hersbach et
al. (2020). A lot of progress have been made in downscaling the climate
simulators with observational data Kotamarthi et al. (2021), yet lacking the
method to accurately infer the long-term climate scenarios directly from
climate drivers.

##  5 Experiments

Task (Full) | DMF | NARGP | MFHNP | D-MFD | SF-NP | MFRNP
---|---|---|---|---|---|---
Heat 2 | 0.138 ± 4.0e-8 | 0.31±2.12e-6 | 0.026±4.01e-5 | 0.015±1.42e-5 | 0.308±6.38e-05 | 0.005±3.27e-4
Heat 3 | 0.137±1.23e-7 | 0.309±3.46e-6 | 0.111±4.82e-6 | 0.108±4.85e-8 | 0.307±8.62e-05 | 0.0039±2.94e-4
Heat 5 | 0.135±2.55e-4 | 0.306±1.14e-4 | 0.115±1.22e-2 | 0.106±1e-4 | 0.306±6.90e-05 | 0.0045±2.94e-4
Poisson 2 | 0.107 ± 6.58e-5 | 0.585±9.84e-5 | 0.093±2.55e-4 | 0.07 ±2.99e-4 | 0.575±1.39e-04 | 0.0076±7.49e-4
Poisson 3 | 0.121±1.47e-5 | 0.58 ±1.02e-4 | 0.335±2.37e-5 | 0.101±1.81e-4 | 0.572±1.95e-04 | 0.0073±5.25e-4
Poisson 5 | 0.101±2.24e-3 | 0.571±1.29e-4 | 0.299±8.84e-3 | 0.279±3.35e-3 | 0.571±1.82e-04 | 0.0046±1.2e-4
Fluid | 0.275±4.59e-7 | 0.353±9.28e-4 | 0.234±4.82e-6 | 0.207±1.31e-5 | 0.383±5.53e-05 | 0.129±8.19e-4
Table 1: Performance (nRMSE) comparison of 6 different models applied to the Heat and Poisson simulators with two, three, five fidelities and fluid simulation with Navier-Stokes equation with two fidelities. The full setting means same domain coverage of $\mathcal{D}_{K}^{train}$ and $\mathcal{D}_{K}^{test}$. Task (OOD) | DMF | NARGP | MFHNP | D-MFD | SF-NP | MFRNP
---|---|---|---|---|---|---
Heat 2 | 0.168±1.36e-4 | 0.313±1.30e-4 | 0.033±1.37e-2 | 0.213±1.65e-3 | 0.312±1.14e-04 | 0.005±1.33e-4
Heat 3 | 0.163±6.04e-4 | 0.309±5.00e-4 | 0.143±5.88e-3 | 0.141±4.94e-3 | 0.310±7.03e-05 | 0.004±3.82e-4
Heat 5 | 0.187±8.60e-4 | 0.308±1.05e-4 | 0.15±2.72e-3 | 0.145±2.87e-3 | 0.308±5.25e-05 | 0.012±1.05e-2
Poisson 2 | 0.183±6.85e-4 | 0.749±1.47e-5 | 0.103±2.03e-2 | 0.214±4.28e-2 | 0.749±1.79e-03 | 0.017±2.72e-3
Poisson 3 | 0.186±8.06e-4 | 0.744±8.24e-5 | 0.189±1.14e-2 | 0.2±1e-2 | 0.745±7.62e-04 | 0.018±1.27e-3
Poisson 5 | 0.16±4.47e-4 | 0.743±3.07e-4 | 0.399±1.07e-2 | 0.375±8.56e-3 | 0.744±1.85e-03 | 0.013±2.76e-4
Table 2: Performance (nRMSE) comparison of 6 different models applied to the
Heat and Poisson simulators with two, three, five fidelities. The OOD setting
here indicates $\mathcal{D}_{K}^{test}$ is OOD w.r.t the training domain at
fidelity $K$.

###  5.1 Datasets

We include 6 Partial Differential Equation benchmark tasks, a more complicated
fluid simulation task, as well as a real-world earth surface temperature
prediction task.

#### Partial Differential Equations (PDEs).

We include Heat and Poisson’s equations Olsen-Kettle (2011) from computational
physics. We test MFRNP for predicting the spatial solution fields of these
equations. We use numerical solvers to generate the ground-truth data with
dense and coarse meshes for different fidelity levels. For both Heat and
Poisson’s equations with $2$ fidelity setting, we use $16\times 16$ and
$32\times 32$ meshes as low and high fidelities. For three fidelity scenarios,
we run the solvers with $64\times 64$ meshes at the highest fidelity level.
For five fidelity scenarios, we run the solvers with $96\times 96$ and
$128\times 128$ meshes as the two additional fidelity levels. The heat
equation has an input dimension of 3, corresponding to the thermal diffusivity
coefficient and boundary conditions at the two edges. The Poisson’s equation
has an input dimension of 5, corresponding to the $4$ boundary conditions and
the magnitude of flow at the centered point source.

#### Fluid Simulation.

This task simulates the dynamics of a circular smoke cloud propelled by an
inflow force within a $50\times 50$ grid. The spatial-temporal solution field
is derived through the application of the incompressible Navier-Stokes
equations and the Boussinesq approximation Holl et al. (2020); Chorin (1968).
Initiated with the introduction of an incompressible static smoke cloud of
radius $5$ at the lower center of the grid, a persistent inflow force is
subsequently applied at the center of the cloud. The input parameters
encompass two variables governing the magnitude of the inflow force along the
$x$ and $y$ directions. The output is the initial component of the velocity
field following a temporal evolution of $30$ steps with the inflow force. The
simulation generates low-fidelity ground-truth with $32\times 32$ mesh, and
high-fidelity with $64\times 64$ mesh.

#### Climate Modeling: Earth Surface Temperature.

In this task, we take one step further from previous works in Section 4 with
multi-fidelity surrogate modeling to directly learn the mapping from climate
drivers Watson-Parris et al. (2022) to the observation-calibrated
ERA5-reanalysis data Hersbach et al. (2020), together with $13$ low-resolution
computational climate model predictions from the ScenarioMIP projectO’Neill et
al. (2016) at global scale. We provide details about the computational models
used in Appendix A.2. We group these climate data into $9$ fidelities based on
their original resolutions. The fidelity dimensions are: $144\times 192$,
$160\times 320$, $192\times 288$, $180\times 288$, $120\times 180$, $143\times
144$, $80\times 96$, $192\times 384$ and $721\times 1440$ at the highest
fidelity level, representing air temperature at 2 meters above the earth
surface. The climate drivers as inputs here consists of $12$ variables
representing the total emission of greenhouse gases (CO2, CH4) and aerosol
gases (BC, SO2) from year $1850$ to $2015$. Each of the aerosol gas is split
into $5$ signals via principal component analysis Wold et al. (1987); Watson-
Parris et al. (2022). For years after $2015$, we test $4$ hypothetical global
gas emission scenarios: ssp126, ssp245, ssp370, ssp585. Larger numbers
correspond to more total gas emissions, leading to more severe climate change.
Ssp126 represents the condition in which future gas emissions are well
controlled and gradually decrease. This scenario makes the gas emission data
as input to have similar domain coverage as of training data. Whereas in the
other three scenarios, future gas emissions would keep increasing, introducing
out-of-distribution inputs at the highest fidelity.

Figure 2: Perfect model test performance for climate modeling task. Measured
in latitude-weighted nRMSE. MFRNP outperforms other models in 3 out of the 4
scenarios and maintains consistent performance across all the scenarios from
year $2015\sim 2100$.

###  5.2 Experiment Setup

We implement MFRNP with PyTorch Paszke et al. (2019) and compare the average
nRMSE in three random runs with following baselines. Details are provided in
Appendix A.4.

  * •

SF-NP. The naive single-fidelity neural process trained on the highest
fidelity data, as the lower performance bound for multi-fidelity modeling.

  * •

DMFLi et al. (2021) performs multi-fidelity learning based on a multi neural
network (NN) structure. Each NN corresponds to one fidelity and the NN on the
next fidelity adapts the latent output from the NN at current fidelity to
propagate information.

  * •

NARGPPerdikaris et al. (2017) utilizes multi-fidelity Gaussian Processes by
sequentially train Gaussian Processes from low to high fidelity levels under
the assumption of a nested input domain.

  * •

MFHNPWu et al. (2022) learns distribution over functions at each fidelity
level with a latent variable. They designed a hierarchical structure of the
latent variables from low to high fidelity levels to pass information.

  * •

D-MFDWu et al. (2023) tackles the error propagation problem in MFHNP and the
inter-fidelity NN overfitting problem in DMF by introducing local and global
latent representations at each fidelity.

#### PDEs Setup.

For PDE tasks, we consider two data composition scenarios.

  * •

In-distribution (Full) Scenario. In this setting, the training set
$\mathcal{D}^{train}$ and testing set $\mathcal{D}^{test}$ at each fidelity
level have the same coverage on dataset domain $\mathcal{D}$.

  * •

Out-of-distribution (OOD) Scenario. This scenario simulates real-world
conditions where the highest fidelity data has limited domain coverage. We
evaluate model trained using $\mathcal{D}^{train}\subset\mathcal{D}$ at the
highest fidelity level, and test the model performance on $\mathcal{D}^{test}$
consists of data covering the rest of $\mathcal{D}$ s.t.
$\mathcal{D}^{train}\cap\mathcal{D}^{test}=\emptyset$ and
$\mathcal{D}^{train}\cup\mathcal{D}^{test}=\mathcal{D}$.

We construct the OOD scenario with Heat and Poisson’s equations. We did not
include Fluid data because the low fidelity mesh is not detailed enough to
capture the differences between $\mathcal{D}^{train}$ and
$\mathcal{D}^{test}$. Under the OOD setting, information covering the input
domain of $\mathcal{D}^{test}$ is missing in all the fidelities, making it
impossible for models to extrapolate with lower fidelity information. For Heat
equation, we test the OOD scenario of the thermal diffusivity coefficient by
limiting the corresponding input parameter to cover $80\%$ of original domain.
We uniformly sample $X$ in the constrained domain and obtain corresponding $Y$
with numerical solver as $\mathcal{D}^{train}$. $\mathcal{D}^{test}$ is
constructed with the rest $20\%$ of original domain. Similarly, for Poisson’s
equation, we construct $\mathcal{D}^{train}$ and $\mathcal{D}^{test}$ with the
$80/20$ coverage split on $3$ boundary conditions. Details are demonstrated in
Appendix A.3.

Figure 3: Absolute residual plot between ground truth and model predictions on
ERA5-reanalysis data at year $2020$, with global resolution of $721\times
1440$. The first row presents global view, second row focuses on smaller
geographic regions. The third row shows the latitude weighted nRMSE to measure
prediction accuracy from year 2015 to 2021. MFRNP outperforms other models
staring from $2017$.

#### Climate Modeling Setup.

In this task, we also compare MFRNP with recent works in the climate science
community of downscaling global climate projections. We include deep kernel
GPs from NNGPR Harris et al. (2023) and Gaussian CNN from DeepSED Baño-Medina
et al. (2022). Since no previous work simultaneously downscales while modeling
the mapping between climate drivers and scenarios, to ensure fair comparison
under the same task, we set up a naive extension for these methods. We first
fit $13$ linear surrogate models with the climate drivers as input $X$, and
the computational climate model predictions as target $Y$. Then we train the
models in NNGPR and DeepSED to downscale the averaged predictions from linear
surrogates given $X$.

We use $13$ "perfect model" tests Knutti et al. (2017) as the performance
metric. In each test, one climate model is held out as the ground truth and we
only use its historical data ($1950\sim 2015$) for training at the highest
fidelity. The remaining $12$ climate models are grouped at lower fidelity
levels based on their resolutions. We compare the average latitude weighted
nRMSE across $13$ tests. Thus, we are able to empirically measure the
performance with physically based future projections. We then include the
ERA5-reanalysis data as the highest fidelity level to measure the performance
on a refined resolution. We use year $1950\sim 2014$ as the highest fidelity
training data and $2015\sim 2021$ for testing.

###  5.3 Results

#### PDEs Performance.

We test the performance of MFRNP and baselines across $7$ tasks under the full
setting and $6$ tasks under the OOD setting. Table 1 shows the full setting
performance and table 2 shows the OOD setup performance. Among all the tasks
and settings, our model consistently outperforms baselines by an average of
$\sim 90\%$ in nRMSE.

The results indicate that MFRNP is more efficient in utilizing information
across fidelities to generate good predictions. Unlike other baselines for
which the performance generally gets worse as the number of fidelities
increase, the performance of MFRNP in the $5$ fidelity tasks is even better
than the corresponding $2$ fidelity tasks in $3$ out of $4$ five-fidelity
experiments. This indicates MFRNP gathers lower fidelity information
efficiently enough to outrun the increasing task complexity in modeling Heat
and Poisson’s equations.

#### Climate Modeling Performance.

As shown in Figure 2, for the perfect model test, MFRNP performs consistently
across all emission scenarios and outperforms other baselines in 3 out of the
4 scenarios. The perfect model test empirically shows the capability of MFRNP
to efficiently aggregate information from lower fidelities and its versatility
for OOD scenarios. Although CNN and Gaussian CNN perform relatively well, they
have very limited scalability and cannot scale to the global ERA5-reanalysis
dataset. We also had memory issues running NARGP and DMF with 9 fidelities for
this task. Among the baseline models that scale to the ERA5-reanalysis dataset
(MFHNP, D-MFD and NNGPR), the performance of MFRNP is significantly better.

For the ERA5-reanalysis data, MFRNP also outperforms other baseline models, as
shown in Figure 3. Our model keeps a consistent performance from year
$2015\sim 2021$, while the performance of other baseline models worsens for
further future predictions. MFRNP makes a cleaner residual plot compared to
the other models. Although we can only evaluate performance up to year $2021$
due to observation limitations, we have empirically shown the long term
prediction consistency of MFRNP from the perfect model test and we expect it
to further outperform other baselines for years after $2021$ with the
ERA5-reanalysis data.

We have demonstrated that MFRNP is not only capable of fully exploiting lower
fidelity information for reliable predictions, but also able to provide great
scalability with consistent performance to solve real-world problems. We also
conduct ablation study to thoroughly evaluate the effectiveness of our
aggregation step by changing lower-fidelity model structures. The results are
included in Appendix A.1.

##  6 Conclusion

In this paper, we present Multi-fidelity Residual Neural Processes (MFRNP), a
novel Neural Process-based multi-fidelity surrogate model. MFRNP utilizes a
residual modeling framework with our tailored Residual-ELBO loss. This allows
MFRNP to leverage the rich input space coverage from lower fidelities while
preserving accuracy from the highest fidelity data. Experimental results on
partial differential equations and climate modeling demonstrate that MFRNP
outperforming state-of-the-art methods by more than $90\%$, highlighting
scalability, efficiency in information fusion, and versatility of MFRNP for
real-world surrogate modeling. For computationally expensive encoders and
decoders, such as neural ODEs, the inference time penalty is doubled for MFRNP
due to ancestral sampling. For future work, we plan to address this limitation
and incorporate physics-informed structure to better model expensive physics
simulations.

## Impact Statement

The work in this paper aims to advance the performance and scalability of
multi-fidelity surrogate modeling. The model proposed has a potential impact
on climate science, computational chemistry, structural engineering, and many
other communities, none of which we foresee potential negative impacts.

## Acknowledgement

This work was supported in part by the U.S. Army Research Office under Army-
ECASE award W911NF-07-R-0003-03, the U.S. Department Of Energy, Office of
Science, IARPA HAYSTAC Program, CDC-RFA-FT-23-0069, NSF Grants #2205093,
#2146343, #2134274, SCALE MoDL-2134209, CCF-2112665 (TILOS) and DARPA AIE .

## References

  * Hersbach et al. (2020) Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al.  The era5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.
  * Karami and Kashef (2020) Zahra Karami and Rasha Kashef.  Smart transportation planning: Data, models, and algorithms.  _Transportation Engineering_ , 2:100013, 2020.
  * Peherstorfer et al. (2018) Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  _Siam Review_ , 60(3):550–591, 2018.
  * Seeger (2004) Matthias Seeger.  Gaussian processes for machine learning.  _International journal of neural systems_ , 14(02):69–106, 2004.
  * Perdikaris et al. (2016) Paris Perdikaris, Daniele Venturi, and George Em Karniadakis.  Multifidelity information fusion algorithms for high-dimensional systems and massive data sets.  _SIAM Journal on Scientific Computing_ , 38(4):B521–B538, 2016.
  * Wang et al. (2021) Zheng Wang, Wei Xing, Robert Kirby, and Shandian Zhe.  Multi-fidelity high-order gaussian processes for physical simulation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 847–855. PMLR, 2021.
  * Williams and Rasmussen (1995) Christopher Williams and Carl Rasmussen.  Gaussian processes for regression.  _Advances in neural information processing systems_ , 8, 1995.
  * Rasmussen (2003) Carl Edward Rasmussen.  Gaussian processes in machine learning.  In _Summer school on machine learning_ , pages 63–71. Springer, 2003.
  * Damianou and Lawrence (2013) Andreas Damianou and Neil D Lawrence.  Deep gaussian processes.  In _Artificial intelligence and statistics_ , pages 207–215. PMLR, 2013.
  * Raissi and Karniadakis (2016) Maziar Raissi and George Karniadakis.  Deep multi-fidelity gaussian processes.  _arXiv preprint arXiv:1604.07484_ , 2016.
  * Salimbeni and Deisenroth (2017) Hugh Salimbeni and Marc Peter Deisenroth.  Doubly stochastic variational inference for deep gaussian processes.  In _NIPS_ , 2017.
  * Wilson et al. (2016) Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing.  Deep kernel learning.  In _Artificial intelligence and statistics_ , pages 370–378. PMLR, 2016.
  * Wang and Lin (2020) Yating Wang and Guang Lin.  Mfpc-net: Multi-fidelity physics-constrained neural process.  _arXiv preprint arXiv:2010.01378_ , 2020.
  * Hebbal et al. (2021) Ali Hebbal, Loic Brevault, Mathieu Balesdent, El-Ghazali Talbi, and Nouredine Melab.  Multi-fidelity modeling with different input domain definitions using deep gaussian processes.  _Structural and Multidisciplinary Optimization_ , 63(5):2267–2288, 2021.
  * Li et al. (2022a) Shibo Li, Zheng Wang, Robert Kirby, and Shandian Zhe.  Deep multi-fidelity active learning of high-dimensional outputs.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1694–1711. PMLR, 2022a.
  * Li et al. (2022b) Shibo Li, Jeff Phillips, Xin Yu, Robert Kirby, and Shandian Zhe.  Batch multi-fidelity active learning with budget constraints.  In _Advances in Neural Information Processing Systems_ , 2022b.
  * Wu et al. (2023) Dongxia Wu, Ruijia Niu, Matteo Chinazzi, Yian Ma, and Rose Yu.  Disentangled multi-fidelity deep bayesian active learning.  _arXiv preprint arXiv:2305.04392_ , 2023.
  * Garnelo et al. (2018a) Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo J Rezende, SM Eslami, and Yee Whye Teh.  Neural processes.  _arXiv preprint arXiv:1807.01622_ , 2018a.
  * Wang and Van Hoof (2020) Qi Wang and Herke Van Hoof.  Doubly stochastic variational inference for neural processes with hierarchical latent variables.  In _International Conference on Machine Learning_ , pages 10018–10028. PMLR, 2020.
  * Jha et al. (2022) Saurav Jha, Dong Gong, Xuesong Wang, Richard E Turner, and Lina Yao.  The neural process family: Survey, applications and perspectives.  _arXiv preprint arXiv:2209.00517_ , 2022.
  * Brevault et al. (2020) Loïc Brevault, Mathieu Balesdent, and Ali Hebbal.  Overview of gaussian process based multi-fidelity techniques with variable relationship between fidelities, application to aerospace systems.  _Aerospace Science and Technology_ , 107:106339, 2020.
  * Hosking (2020) S. Hosking.  Multifidelity climate modelling, github.  https://github.com/scotthosking/mf_modelling, 2020.
  * Valero et al. (2021) Mario Miguel Valero, Lluís Jofre, and Ricardo Torres.  Multifidelity prediction in wildfire spread simulation: Modeling, uncertainty quantification and sensitivity analysis.  _Environmental Modelling & Software_, 141:105050, 2021.
  * Kennedy and O’Hagan (2000) Marc C Kennedy and Anthony O’Hagan.  Predicting the output from a complex computer code when fast approximations are available.  _Biometrika_ , 87(1):1–13, 2000.
  * Le Gratiet and Garnier (2014) Loic Le Gratiet and Josselin Garnier.  Recursive co-kriging model for design of computer experiments with multiple levels of fidelity.  _International Journal for Uncertainty Quantification_ , 4(5), 2014.
  * Perdikaris et al. (2015) Paris Perdikaris, Daniele Venturi, Johannes O Royset, and George Em Karniadakis.  Multi-fidelity modelling via recursive co-kriging and gaussian–markov random fields.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 471(2179):20150018, 2015.
  * Perdikaris et al. (2017) Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D Lawrence, and George Em Karniadakis.  Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2198):20160751, 2017.
  * Cutajar et al. (2019) Kurt Cutajar, Mark Pullin, Andreas Damianou, Neil Lawrence, and Javier González.  Deep gaussian processes for multi-fidelity modeling.  _arXiv preprint arXiv:1903.07320_ , 2019.
  * Li et al. (2022c) Shibo Li, Zheng Wang, Robert Kirby, and Shandian Zhe.  Infinite-fidelity coregionalization for physical simulation.  _Advances in Neural Information Processing Systems_ , 35:25965–25978, 2022c.
  * Guo et al. (2022) Mengwu Guo, Andrea Manzoni, Maurice Amendt, Paolo Conti, and Jan S Hesthaven.  Multi-fidelity regression using artificial neural networks: efficient approximation of parameter-dependent output quantities.  _Computer methods in applied mechanics and engineering_ , 389:114378, 2022.
  * Meng and Karniadakis (2020) Xuhui Meng and George Em Karniadakis.  A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse pde problems.  _Journal of Computational Physics_ , 401:109020, 2020.
  * Meng et al. (2021) Xuhui Meng, Hessam Babaee, and George Em Karniadakis.  Multi-fidelity bayesian neural networks: Algorithms and applications.  _Journal of Computational Physics_ , 438:110361, 2021.
  * De et al. (2020) Subhayan De, Jolene Britton, Matthew Reynolds, Ryan Skinner, Kenneth Jansen, and Alireza Doostan.  On transfer learning of neural networks using bi-fidelity data for uncertainty propagation.  _International Journal for Uncertainty Quantification_ , 10(6), 2020.
  * Li et al. (2020) Shibo Li, Wei Xing, Robert Kirby, and Shandian Zhe.  Multi-fidelity bayesian optimization via deep neural networks.  _Advances in Neural Information Processing Systems_ , 33:8521–8531, 2020.
  * Perry et al. (2019) Daniel J Perry, Robert M Kirby, Akil Narayan, and Ross T Whitaker.  Allocation strategies for high fidelity models in the multifidelity regime.  _SIAM/ASA Journal on Uncertainty Quantification_ , 7(1):203–231, 2019.
  * Kandasamy et al. (2017) Kirthevasan Kandasamy, Gautam Dasarathy, Jeff Schneider, and Barnabás Póczos.  Multi-fidelity bayesian optimisation with continuous approximations.  In _International Conference on Machine Learning_ , pages 1799–1808. PMLR, 2017.
  * Garnelo et al. (2018b) Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and SM Ali Eslami.  Conditional neural processes.  In _International Conference on Machine Learning_ , pages 1704–1713. PMLR, 2018b.
  * Kim et al. (2018) Hyunjik Kim, Andriy Mnih, Jonathan Schwarz, Marta Garnelo, Ali Eslami, Dan Rosenbaum, Oriol Vinyals, and Yee Whye Teh.  Attentive neural processes.  In _International Conference on Learning Representations_ , 2018.
  * Louizos et al. (2019) Christos Louizos, Xiahan Shi, Klamer Schutte, and Max Welling.  The functional neural process.  _Advances in Neural Information Processing Systems_ , 2019.
  * Singh et al. (2019) Gautam Singh, Jaesik Yoon, Youngsung Son, and Sungjin Ahn.  Sequential neural processes.  _Advances in Neural Information Processing Systems_ , 32:10254–10264, 2019.
  * Flato et al. (2014) Gregory Flato, Jochem Marotzke, Babatunde Abiodun, Pascale Braconnot, Sin Chan Chou, William Collins, Peter Cox, Fatima Driouech, Seita Emori, Veronika Eyring, et al.  Evaluation of climate models.  In _Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change_ , pages 741–866. Cambridge University Press, 2014.
  * Knutti et al. (2010) Reto Knutti, Reinhard Furrer, Claudia Tebaldi, Jan Cermak, and Gerald A Meehl.  Challenges in combining projections from multiple climate models.  _Journal of Climate_ , 23(10):2739–2758, 2010.
  * Lambert and Boer (2001) Steven J Lambert and George J Boer.  Cmip1 evaluation and intercomparison of coupled climate models.  _Climate Dynamics_ , 17:83–106, 2001.
  * Gleckler et al. (2008) Peter J Gleckler, Karl E Taylor, and Charles Doutriaux.  Performance metrics for climate models.  _Journal of Geophysical Research: Atmospheres_ , 113(D6), 2008.
  * Tebaldi and Knutti (2007) Claudia Tebaldi and Reto Knutti.  The use of the multi-model ensemble in probabilistic climate projections.  _Philosophical transactions of the royal society A: mathematical, physical and engineering sciences_ , 365(1857):2053–2075, 2007.
  * Giorgi and Mearns (2002) Filippo Giorgi and Linda O Mearns.  Calculation of average, uncertainty range, and reliability of regional climate changes from aogcm simulations via the “reliability ensemble averaging”(rea) method.  _Journal of climate_ , 15(10):1141–1158, 2002.
  * Abramowitz et al. (2019) Gab Abramowitz, Nadja Herger, Ethan Gutmann, Dorit Hammerling, Reto Knutti, Martin Leduc, Ruth Lorenz, Robert Pincus, and Gavin A Schmidt.  Esd reviews: Model dependence in multi-model climate ensembles: Weighting, sub-selection and out-of-sample testing.  _Earth System Dynamics_ , 10(1):91–105, 2019.
  * Baño-Medina et al. (2022) Jorge Baño-Medina, Rodrigo Manzanas, Ezequiel Cimadevilla, Jesús Fernández, Jose González-Abad, Antonio S Cofiño, and José Manuel Gutiérrez.  Downscaling multi-model climate projection ensembles with deep learning (deepesd): contribution to cordex eur-44.  _Geoscientific Model Development_ , 15(17):6747–6758, 2022.
  * Dee et al. (2011) Dick P Dee, S M Uppala, Adrian J Simmons, Paul Berrisford, Paul Poli, Shinya Kobayashi, U Andrae, MA Balmaseda, G Balsamo, d P Bauer, et al.  The era-interim reanalysis: Configuration and performance of the data assimilation system.  _Quarterly Journal of the royal meteorological society_ , 137(656):553–597, 2011.
  * Harris et al. (2023) Trevor Harris, Bo Li, and Ryan Sriver.  Multimodel ensemble analysis with neural network gaussian processes.  _The Annals of Applied Statistics_ , 17(4):3403–3425, 2023.
  * Kotamarthi et al. (2021) Rao Kotamarthi, Katharine Hayhoe, Donald Wuebbles, Linda O Mearns, Jennifer Jacobs, and Jennifer Jurado.  _Downscaling techniques for high-resolution climate projections: From global change to local impacts_.  Cambridge University Press, 2021.
  * Olsen-Kettle (2011) Louise Olsen-Kettle.  Numerical solution of partial differential equations.  _Lecture notes at University of Queensland, Australia_ , 2011.
  * Holl et al. (2020) Philipp Holl, Vladlen Koltun, and Nils Thuerey.  Learning to control pdes with differentiable physics.  _arXiv preprint arXiv:2001.07457_ , 2020.
  * Chorin (1968) Alexandre Joel Chorin.  Numerical solution of the navier-stokes equations.  _Mathematics of computation_ , 22(104):745–762, 1968.
  * Watson-Parris et al. (2022) Duncan Watson-Parris, Yuhan Rao, Dirk Olivié, Øyvind Seland, Peer Nowack, Gustau Camps-Valls, Philip Stier, Shahine Bouabid, Maura Dewey, Emilie Fons, et al.  Climatebench v1. 0: A benchmark for data-driven climate projections.  _Journal of Advances in Modeling Earth Systems_ , 14(10):e2021MS002954, 2022.
  * O’Neill et al. (2016) Brian C O’Neill, Claudia Tebaldi, Detlef P Van Vuuren, Veronika Eyring, Pierre Friedlingstein, George Hurtt, Reto Knutti, Elmar Kriegler, Jean-Francois Lamarque, Jason Lowe, et al.  The scenario model intercomparison project (scenariomip) for cmip6.  _Geoscientific Model Development_ , 9(9):3461–3482, 2016.
  * Wold et al. (1987) Svante Wold, Kim Esbensen, and Paul Geladi.  Principal component analysis.  _Chemometrics and intelligent laboratory systems_ , 2(1-3):37–52, 1987.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Li et al. (2021) Shibo Li, Zheng Wang, Robert M. Kirby, and Shandian Zhe.  Deep multi-fidelity active learning of high-dimensional outputs, 2021.
  * Wu et al. (2022) Dongxia Wu, Matteo Chinazzi, Alessandro Vespignani, Yi-An Ma, and Rose Yu.  Multi-fidelity hierarchical neural processes.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , pages 2029–2038, 2022.
  * Knutti et al. (2017) Reto Knutti, Jan Sedláček, Benjamin M Sanderson, Ruth Lorenz, Erich M Fischer, and Veronika Eyring.  A climate model projection weighting scheme accounting for performance and interdependence.  _Geophysical Research Letters_ , 44(4):1909–1918, 2017.
  * Bi et al. (2020) Daohua Bi, Martin Dix, Simon Marsland, Siobhan O’farrell, Arnold Sullivan, Roger Bodman, Rachel Law, Ian Harman, Jhan Srbinovsky, Harun A Rashid, et al.  Configuration and spin-up of access-cm2, the new generation australian community climate and earth system simulator coupled model.  _Journal of Southern Hemisphere Earth Systems Science_ , 70(1):225–251, 2020.
  * Wu et al. (2021) Tongwen Wu, Rucong Yu, Yixiong Lu, Weihua Jie, Yongjie Fang, Jie Zhang, Li Zhang, Xiaoge Xin, Laurent Li, Zaizhi Wang, et al.  Bcc-csm2-hr: a high-resolution version of the beijing climate center climate system model.  _Geoscientific Model Development_ , 14(5):2977–3006, 2021.
  * Cherchi et al. (2019) Annalisa Cherchi, Pier Giuseppe Fogli, Tomas Lovato, Daniele Peano, Doroteaciro Iovino, Silvio Gualdi, Simona Masina, Enrico Scoccimarro, Stefano Materia, Alessio Bellucci, et al.  Global mean climate and main patterns of variability in the cmcc-cm2 coupled model.  _Journal of Advances in Modeling Earth Systems_ , 11(1):185–209, 2019.
  * Dunne et al. (2020) John P Dunne, LW Horowitz, AJ Adcroft, P Ginoux, IM Held, JG John, John P Krasting, Sergey Malyshev, V Naik, F Paulot, et al.  The gfdl earth system model version 4.1 (gfdl-esm 4.1): Overall coupled model description and simulation characteristics.  _Journal of Advances in Modeling Earth Systems_ , 12(11):e2019MS002015, 2020.
  * Volodin (2022) EM Volodin.  Possible climate change in russia in the 21st century based on the inm-cm5-0 climate model.  _Russian Meteorology and Hydrology_ , 47(5):327–333, 2022.
  * Boucher et al. (2020) Olivier Boucher, Jérôme Servonnat, Anna Lea Albright, Olivier Aumont, Yves Balkanski, Vladislav Bastrikov, Slimane Bekki, Rémy Bonnet, Sandrine Bony, Laurent Bopp, et al.  Presentation and evaluation of the ipsl-cm6a-lr climate model.  _Journal of Advances in Modeling Earth Systems_ , 12(7):e2019MS002010, 2020.
  * Byun et al. (2019) Young-Hwa Byun, Yoon-Jin Lim, Sungbo Shim, Hyun Min Sung, Minah Sun, Jisun Kim, Byeong-Hyeon Kim, Jae-Hee Lee, and Hyejin Moon.  Nims-kma kace1. 0-g model output prepared for cmip6 scenariomip.  2019\.
  * Stouffer (2019) Ronald Stouffer.  Ipcc ddc: Ua mcm-ua-1-0 model output prepared for cmip6 scenariomip.  2019\.
  * Yukimoto et al. (2019) Seiji Yukimoto, Hideaki Kawai, Tsuyoshi Koshiro, Naga Oshima, Kohei Yoshida, Shogo Urakawa, Hiroyuki Tsujino, Makoto Deushi, Taichu Tanaka, Masahiro Hosaka, et al.  The meteorological research institute earth system model version 2.0, mri-esm2. 0: Description and basic evaluation of the physical component.  _Journal of the Meteorological Society of Japan. Ser. II_ , 97(5):931–965, 2019.
  * Swart et al. (2019) Neil C Swart, Jason NS Cole, Viatcheslav V Kharin, Mike Lazare, John F Scinocca, Nathan P Gillett, James Anstey, Vivek Arora, James R Christian, Sarah Hanna, et al.  The canadian earth system model version 5 (canesm5. 0.3).  _Geoscientific Model Development_ , 12(11):4823–4873, 2019.
  * Gutjahr et al. (2019) Oliver Gutjahr, Dian Putrasahan, Katja Lohmann, Johann H Jungclaus, Jin-Song von Storch, Nils Brüggemann, Helmuth Haak, and Achim Stössel.  Max planck institute earth system model (mpi-esm1. 2) for the high-resolution model intercomparison project (highresmip).  _Geoscientific Model Development_ , 12(7):3241–3281, 2019.
  * Seland et al. (2020) Øyvind Seland, Mats Bentsen, Dirk Olivié, Thomas Toniazzo, Ada Gjermundsen, Lise Seland Graff, Jens Boldingh Debernard, Alok Kumar Gupta, Yan-Chun He, Alf Kirkevåg, et al.  Overview of the norwegian earth system model (noresm2) and key climate response of cmip6 deck, historical, and scenario simulations.  _Geoscientific Model Development_ , 13(12):6165–6200, 2020.
  * Sellar et al. (2019) Alistair A Sellar, Colin G Jones, Jane P Mulcahy, Yongming Tang, Andrew Yool, Andy Wiltshire, Fiona M O’connor, Marc Stringer, Richard Hill, Julien Palmieri, et al.  Ukesm1: Description and evaluation of the uk earth system model.  _Journal of Advances in Modeling Earth Systems_ , 11(12):4513–4558, 2019.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.

##  Appendix A Appendix

###  A.1 Ablation Study

Task (Full) | MFRNP-H | MFRNP | Task (OOD) | MFRNP-H | MFRNP
---|---|---|---|---|---
Heat 3 | 0.013±4.46e-4 | 0.0039±2.94e-4 | Heat 3 | 0.047±2.72e-4 | 0.004±3.82e-4
Heat 5 | 0.013±2.76e-4 | 0.0045±2.94e-4 | Heat 5 | 0.055±1.87e-4 | 0.012±1.05e-2
Poisson 3 | 0.138±6.1e-3 | 0.0073±5.25e-4 | Poisson 3 | 0.128±1.3e-2 | 0.018±1.27e-3
Poisson 5 | 0.184±1.77e-2 | 0.0046±1.2e-4 | Poisson 5 | 0.151±7.9e-3 | 0.013±2.76e-4
Table 3: Performance (nRMSE) comparison of MFRNP-H and MFRNP applied to the
Heat and Poisson simulators with three, five fidelities. The OOD setting here
indicates $\mathcal{D}_{K}^{test}$ is OOD w.r.t the training domain at
fidelity $K$.

In this study, we create a variant of MFRNP, namely MFRNP-H, for which the
fidelity levels below the highest are implemented with a hierarchical
structure that shares information via the latent space. We use the surrogate
output at fidelity $K-1$ as the aggregated information and model the residual
between this aggregation and the highest fidelity ground truth. This
modification aims to evaluate the impact of including decoders in the
information sharing step. The same experiment settings are maintained for
consistency in the comparison. Our results in Table 3 show that including
decoders in the information aggregation step significantly impacts the
performance.

###  A.2 Climate Modeling Task Dataset

We list the $13$ computational climate models used here, together with the
literatures. They are all included in the ScenarioMIP project O’Neill et al.
[2016]. ACCESS-CM2; BCC-CSM2-MR; CMCC-CM2-SR5; GFDL-ESM4; INM-CM5-0; IPSL-
CM6A-LR; KACE-1-0-G; MCM-UA-1-0; MRI-ESM2-0; CESM5; MPI-ESM; NorESM2; UKESM.Bi
et al. [2020], Wu et al. [2021], Cherchi et al. [2019], Dunne et al. [2020],
Volodin [2022], Boucher et al. [2020], Byun et al. [2019], Stouffer [2019],
Yukimoto et al. [2019], Swart et al. [2019], Gutjahr et al. [2019], Seland et
al. [2020], Sellar et al. [2019]

###  A.3 OOD Dataset Construction

Generally, we follow the experiment setup as the full scenario. The following
describes the difference in building the highest fidelity dataset
$\mathcal{D}_{K}$. For Heat equation, we uniformly sample $X_{train}$ in the
training scope: ((0,0.8),(-1,0),(0.01,0.1)) and use numerical solver to
generate the corresponding $Y_{train}$. We then uniformly sample $X_{test}$ in
the testing scope ((0.8, 1),(-1,0),(0.01,0.1)) and use numerical solver to
generate $Y_{test}$.

For the Poisson’s equation, we uniformly sample $X_{train}$ in the training
scope: ((0.1, 0.74),(0.1, 0.74),(0.1, 0.74),(0.1, 0.9),(0.1, 0.9)) and use
numerical solver to generate the corresponding $Y_{train}$. We then uniformly
sample $X_{test}$ in the testing scope ((0.74, 0.9),(0.74, 0.9),(0.74,
0.9),(0.1, 0.9),(0.1, 0.9)) and use numerical solver to generate $Y_{test}$.

Training of the models is conducted with a patience threshold of 10,000 steps
and a maximum of 40,000 epochs. We use $10\%$ of training data as validation.

$W\\_nRMSE=\frac{\sqrt{\sum_{i=1}^{N}\left(\frac{\cos({\text{lat}_{i}})}{\frac{\sum_{j=1}^{N}\cos({\text{lat}_{j}})}{N}}\right)*\left(\frac{(y_{i}-\hat{y}_{i})^{2}}{N}\right)}}{std(\\{y_{1:N}\\})}$
(8)

$nRMSE=\frac{\sqrt{\sum_{i=1}^{N}\left(\frac{(y_{i}-\hat{y}_{i})^{2}}{N}\right)}}{std(\\{y_{1:N}\\})}$
(9)

###  A.4 Experimental Setup

#### Metrics.

We use latitude weighted nRMSE as shown in Equation 8 to measure model
performance in the climate modeling task. For PDE tasks, we use nRMSE as shown
in Equation 9.

#### PDE task Training Configurations.

For training, we use Adam optimizer Kingma and Ba [2014] with base learning
rate of $1e-3$. We use $10\%$ of training data set validation set, and run our
model with maximum epoch of $40000$ and patience $10000$. We use learning rate
decay of $0.85$ and stepsize $10000$. For context-target split, we randomly
select $20\%\sim 25\%$ of training data as our context set, the rest as our
target set for each fidelity. For Heat and Poisson’s equation, we run our
model with latent dimension and encoder/decoder dimension of $32$. For fluid
simulation, we use latent dimension and encoder/decoder dimension $128$. We
normalize the data before training and measure nRMSE on the de-normalized
space.

For the baseline models, we follow the same setup as above.

#### Climate Modeling Training Configurations.

We follow the similar setup as above, but set the latent and hidden dimensions
to $512$ to incorporate more fidelities and data at higher dimensions. We do
not normalize the data before training. For climate methods, we use their
original model setup. The linear surrogates for low-resolution climate
simulators are built with $4$ layers with ReLU activation. Dimensions are
$12$, (lat*lon)//128, (lat*lon)//32, (lat*lon)//4. Here, (lat,lon) refers to
the climate simulator data dimensions. These models are run until convergence.