[table]capposition=top

# Probabilistic Neural Data Fusion for Learning from an Arbitrary Number of
Multi-fidelity Data Sets

Carlos Mora† Department of Mechanical and Aerospace Engineering, University of
California, Irvine  Jonathan Tammer Eweis-Labolle Equal Contribution.
Department of Mechanical and Aerospace Engineering, University of California,
Irvine  Tyler Johnson  Department of Mechanical and Aerospace Engineering,
University of California, Irvine  Likith Gadde  Northwood High School, Irvine
Ramin Bostanabad Corresponding Author<EMAIL_ADDRESS>
GitLab repository: https://gitlab.com/TammerUCI/pro-ndf Department of
Mechanical and Aerospace Engineering, University of California, Irvine

Abstract

In many applications in engineering and sciences analysts have simultaneous
access to multiple data sources. In such cases, the overall cost of acquiring
information can be reduced via data fusion or multi-fidelity (MF) modeling
where one leverages inexpensive low-fidelity (LF) sources to reduce the
reliance on expensive high-fidelity (HF) data. In this paper, we employ neural
networks (NNs) for data fusion in scenarios where data is very scarce and
obtained from an arbitrary number of sources with varying levels of fidelity
and cost. We introduce a unique NN architecture that converts MF modeling into
a nonlinear manifold learning problem. Our NN architecture inversely learns
non-trivial (e.g., non-additive and non-hierarchical) biases of the LF sources
in an interpretable and visualizable manifold where each data source is
encoded via a low-dimensional distribution. This probabilistic manifold
quantifies model form uncertainties such that LF sources with small bias are
encoded close to the HF source. Additionally, we endow the output of our NN
with a parametric distribution not only to quantify aleatoric uncertainties,
but also to reformulate the network’s loss function based on strictly proper
scoring rules which improve robustness and accuracy on unseen HF data. Through
a set of analytic and engineering examples, we demonstrate that our approach
provides a high predictive power while quantifying various sources
uncertainties. Our codes and examples can be accessed via GitLab.

Keywords: Multi-fidelity Modeling; Uncertainty Quantification; Bayesian Neural
Networks; Inverse Problems; Manifold Learning; Data Fusion.

##  1 Introduction

In an increasing number of applications in engineering and sciences analysts
have simultaneous access to multiple sources of information. For instance,
materials’ properties can be estimated via multiple techniques such as (in
decreasing order of cost and accuracy/fidelity) experiments, direct numerical
simulations (DNS), a host of physics-based reduced order models (ROMs), or
analytical methods [1, 2, 3]. In such applications, the overall cost of
gathering information about the system of interest can be reduced via multi-
fidelity (MF) modeling or data fusion where one leverages inexpensive low-
fidelity (LF) sources to reduce the reliance on expensive high-fidelity (HF)
data sources. In this paper, we employ neural networks (NNs) for MF modeling
in scenarios where data is scarce and obtained from multiple sources with
varying levels of fidelity and cost (i.e., data is unbalanced since more
samples are available from cheaper sources). In particular, our contributions
are as follows $(1)$ we introduce a unique NN architecture that not only
facilitates data fusion, but also quantifies and visualizes the
discrepancies/similarities between all data sources, and $(2)$ we illustrate
that a Bayesian treatment, besides alleviating overfitting and providing a
probabilistic surrogate (i.e., an emulator), provides the means to develop a
novel loss function (based on proper scoring rules) that improves the
performance and robustness of the resulting MF NN emulator.

Over the past few decades, many techniques have been developed for building MF
surrogates which are used in outer-loop applications such as design
optimization [4, 5], calibration of computer models [6], or Bayesian
optimization [7]. The main motivation behind these techniques is to leverage
the correlations between LF and HF data sources (and the fact that sampling
from the former is typically cheaper) to improve the predictive performance of
the surrogate while reducing the overall data acquisition costs. Early works
in this field focused primarily on hierarchically linking bi-fidelity data.
For instance, in space mapping [8, 9, 10] or multi-level [11, 12, 13]
techniques the inputs of the LF data are mapped following formulations such as
$\boldsymbol{x}^{l}=F(\boldsymbol{x}^{h}{})$ where $\boldsymbol{x}^{l}$ and
$\boldsymbol{x}^{h}{}$ are the inputs of LF and HF sources, respectively. In
this equation, $F(\cdot)$ is a transformation function whose predefined
functional form is calibrated such that $y^{l}(F(\boldsymbol{x}^{h}{}))$
approximates $y^{h}(\boldsymbol{x}^{h}{})$ as closely as possible. These
techniques are useful in applications where higher fidelity data are obtained
by successively refining the discretization in simulations [11, 12], e.g., by
refining the mesh when modeling the flow around an airfoil or estimating the
fracture toughness of a microstructure. The main disadvantages of space
mapping techniques are that $(1)$ they rely on iterative and time-consuming
analysis for choosing a near-optimal functional form for $F(\cdot)$, $(2)$
they cannot jointly fuse more than two data sources at a time, $(3)$ they
quantify similarity/discrepancy between the sources based on pre-defined
functions whose space may not include the true discrepancy, and $(4)$ they do
not quantify some uncertainty sources (such as lack of data) and are rarely
formulated within a Bayesian setting that leverages prior information.

A well-known hierarchical bi-fidelity modeling framework is that of Kennedy
and O’Hagan (KOH) [14] who assume that the discrepancy between the LF and HF
sources is additive (multiplicative terms have also been explored [15]) and
that both sources as well as the discrepancy between them can be modeled via
Gaussian processes (GPs). Upon this modeling assumption, KOH find the joint
posterior of GPs’ hyperparameters via either fully [16, 17] or modular
Bayesian inference [18, 19, 20, 21]. While KOH’s approach considers multiple
uncertainties and has been successfully applied to a broad range of
applications [22, 23, 24], it has three main limitations: $(1)$ it only
accommodates two data sources at a time, $(2)$ it places a priori independence
assumption between the GPs, and $(3)$ it does not provide a low-dimensional,
visualizable, and interpretable metric that quantifies the correlations
between the data sources.

Recent works have acknowledged the limitations of hierarchical methods and
devised new methodologies to address them. For instance, MF modeling can be
achieved via a recursive scheme [25] where a bi-fidelity method is repeatedly
applied from the lowest to the highest fidelities. However, such recursive
schemes inherit the limitations of bi-fidelity methods, cannot jointly fuse
multi-source data sets, and are sensitive to the ordering (i.e., the relative
accuracy of all sources must be known a priori).

As another example, [26] presents MF networks (MFNets): an approach based on
directed acyclic graphs that builds a MF surrogate using an arbitrary number
of data sources. MFNets accommodate noisy data and are trained via gradient-
based minimization of a nonlinear least squares objective. While MFNets can
learn non-hierarchical relations between data sources, they: $(1)$ rely on
having prior knowledge on a set of latent variables that explain the relations
between the sources, $(2)$ assume each source can be surrogated via a linear
subspace model, $(3)$ are not probabilistic and also require regularization,
$(4)$ impose independence assumption among the data sources to derive the
likelihood (i.e., the objective) function, and $(5)$ rely on iterative
approaches for finding the optimal graph structure.

Other notable works that have studied the limitations of hierarchical
techniques include [27, 28, 29] which are focused on identifying (and
correcting) non-additive discrepancies between LF and HF sources. However, the
proposed solution in these works is intrusive and relies on some rather strong
modeling assumptions that largely limit the applications. These limitations
arise because the formulation of the discrepancy is learned via an embedded
operator whose functional form and interaction with the LF source are
constructed a priori.

We have recently developed a GP-based approach [30] that addresses the above
issues by converting MF modeling into a manifold learning problem where the
relations between the sources are automatically quantified via an
appropriately learnt distance measure. The conversion is achieved via latent
map Gaussian processes [30] (LMGPs, see Section 2.1) which enable GPs to
handle categorical variables and, correspondingly, data fusion: by augmenting
the inputs via a categorical variable (which indicates the source of a data
point) and then concatenating all the data sets, LMGPs can simultaneously
learn from an arbitrary number of information sources. We have shown [30] that
LMGP-based MF modeling consistently outperforms KOH’ approach and can also
handle calibration problems.

Following the success of LMGPs in data fusion, in this work we examine the
potentials of NNs in matching (and, hopefully, improving) LMGPs’ efficiency in
MF modeling. Our current studies are motivated by the facts that $(1)$ when
viewed as (probabilistic or deterministic) graphical models [31], NNs provide
unique opportunities to use MF data sets to uncover complex hidden relations
between the corresponding sources, $(2)$ the recent hardware and software
advancements have dramatically accelerated architecture design and training of
NNs, and $(3)$ NNs scale to higher dimensions and big data significantly
better than GPs.

Over the past few years some NN-based approaches have been developed for MF
modeling [32, 26, 33, 34]. However, most of these works design the network
architecture primarily based on hierarchical methods and consequently inherit
their limitations. For instance, [32] builds two sequentially connected
deterministic networks based on KOH’s method where the first and second NNs
are tasked to emulate the LF and HF sources, respectively. In addition to
sharing the limitations of KOH’s method, such a sequential bi-fidelity NN
requires the LF and HF training data to be available at the same inputs
(unless the two parts of the network are trained separately) and also relies
on manual tuning of the architecture and loss function. It has been argued
[33] that such sequentially trained NNs bridge MF modeling with transfer
learning where the knowledge gained from the LF data is used in building the
NN module that surrogates the HF source.

Non-sequential NNs are rarely used for MF modeling (esp. with $>2$ sources)
due to the fact that searching for the optimum architecture (and effectively
training it with small data) is a difficult task. We address this challenge by
drawing inspiration from LMGPs where we design the architecture such that any
number of MF data sets can be simultaneously fused and the overall
discrepancies between sources are quantified with visualizable metrics. We
also illustrate that making specific parts of the network probabilistic, in
addition to being superior to both deterministic and all-probabilistic NNs,
enables us to infuse a proper scoring rule [35] into the loss function and, in
turn, improve the performance of the MF emulator. The particular rule that we
adopt is interval score which is frequently used in testing the quality of
probabilistic predictions but, to the best of our knowledge, has never been
used in the training stage of a probabilistic NN. In summary, our major
contributions are as follows:

  * •

We introduce a unique NN architecture for MF modeling that can fuse an
arbitrary number of data sets and quantify both epistemic and aleatoric
uncertainties.

  * •

We inversely learn the accuracy of the LF sources (with respect to the HF
source) and visualize the learned relations in an interpretable manifold.

  * •

We show that a probabilistic setting allows us to develop a novel loss
function (based on proper scoring rules) that improves the performance of the
emulator.

  * •

We validate the performance of our approach on analytical and real-world
examples and show that it performs on par with the state of the art while
providing improved scalability to high dimensions and big data.

The rest of this paper is organized as follows. We review the relevant
technical background in Section 2 and then introduce our approach in Section
3. We test the performance of our approach on a host of analytical problems
and real-world data sets in Section 4 and conclude the paper in Section 5.

##  2 Technical Preliminaries

In this section we first review LMGPs which are extensions of GPs that handle
categorical inputs and, thus, can readily fuse any number of data sets. Then,
we provide some background on Bayesian neural networks (BNNs) which form the
foundation of our neural data fusion framework.

###  2.1 Latent Map Gaussian Processes (LMGPs)

Let us denote the output and inputs in the training data by
$y\in\mathcal{Y}\equiv\mathbb{R}$ and
$\boldsymbol{x}=[x_{1},x_{2},\dots,x_{dx}]^{T}\in\mathbb{R}^{dx}$,
respectively, with an individual training point $i=1,\dots,n$ denoted by the
pair $(y^{(i)},\boldsymbol{x}^{(i)})$. Assume the training data is a
realization from a constant-mean111GPs (and LMGPs) can also be formulated by
using a linear combination of basis functions in place of the constant mean.
This formulation relies on prior knowledge of the functional form of the
output and can improve performance in extrapolation, see [30]. GP and that the
following relation holds:

$y(\boldsymbol{x})=m+\xi(\boldsymbol{x})\vspace{-0.2cm}$ (1)

where $m$ is the unknown constant mean and $\xi(\boldsymbol{x})$ is a zero-
mean GP whose covariance function or kernel is:

$\text{cov}(\xi(\boldsymbol{x}),\xi(\boldsymbol{x}^{{}^{\prime}}))=c(\boldsymbol{x},\boldsymbol{x}^{{}^{\prime}})=s^{2}r(\boldsymbol{x},\boldsymbol{x}^{{}^{\prime}})\vspace{-0.2cm}$
(2)

where $s^{2}$ is the variance of the process and $r(\cdot,\cdot)$ is a
parametric correlation function such as the Gaussian:

$r(\boldsymbol{x},\boldsymbol{x}^{{}^{\prime}})=\exp{\left(-\sum^{dx}_{i=1}10^{\omega_{i}}(x_{i}-x_{i}^{{}^{\prime}})^{2}\right)}=\exp{\left(-\left(\boldsymbol{x}-\boldsymbol{x}^{{}^{\prime}}\right)^{T}10^{\boldsymbol{\Omega}}\left(\boldsymbol{x}-\boldsymbol{x}^{{}^{\prime}}\right)\right)}\vspace{-0.2cm}$
(3)

where $\boldsymbol{\omega}=[\omega_{1},\dots,\omega_{dx}]^{T}$ are the
roughness or scale parameters and
$\boldsymbol{\Omega}=diag(\boldsymbol{\omega})$.

The training process and prediction formulas for a GP depend on the choice of
the correlation function, which relies on a weighted Cartesian distance metric
between any two inputs, see Equation 3. As we recently motivated in [36], to
directly use GPs for mixed-variable modeling we reformulate $r(\cdot,\cdot)$
as detailed below such that it can handle categorical (qualitative) inputs.

Let us denote the categorical inputs by
$\boldsymbol{t}=[t_{1},\dots,t_{dt}]^{T}$ where the total number of distinct
levels for qualitative variable $t_{i}$ is $\tau_{i}$. To handle mixed inputs,
LMGP learns a parametric function that maps categorical variables to some
points in a quantitative manifold or latent space222Multiple mapping functions
can also be used to build multiple manifolds. We leverage this in Section 4
where we build two manifolds for data fusion problems with categorical or
mixed inputs.. These points (and hence the mapping function) can be
incorporated into any standard correlation function, such as the Gaussian,
which is reformulated as follows for mixed inputs:

$r\left((\boldsymbol{x},\boldsymbol{t}),(\boldsymbol{x}^{{}^{\prime}},\boldsymbol{t}^{{}^{\prime}})\right)=\exp{\left(-\left\lVert\boldsymbol{z}(\boldsymbol{t})-\boldsymbol{z}(\boldsymbol{t}^{{}^{\prime}})\right\rVert^{2}_{2}-\left(\boldsymbol{x}-\boldsymbol{x}^{{}^{\prime}}\right)^{T}10^{\boldsymbol{\Omega}}\left(\boldsymbol{x}-\boldsymbol{x}^{{}^{\prime}}\right)\right)}\vspace{-0.2cm}$
(4)

or, equivalently,

$r\left((\boldsymbol{x},\boldsymbol{t}),(\boldsymbol{x}^{{}^{\prime}},\boldsymbol{t}^{{}^{\prime}})\right)=\exp{\left(-\sum^{dx}_{i=1}10^{\omega_{i}}(x_{i}-x_{i}^{{}^{\prime}})^{2}\right)}\times\exp{\left(-\sum^{dz}_{i=1}(z_{i}(\boldsymbol{t})-z_{i}(\boldsymbol{t}^{{}^{\prime}}))^{2}\right)}\vspace{-0.2cm}$
(5)

where $\left\lVert\cdot\right\rVert_{2}$ denotes the Euclidean 2-norm and
$\boldsymbol{z}(\boldsymbol{t})=[z_{1}(\boldsymbol{t}),\dots,z_{dz}(\boldsymbol{t})]_{1\times
dz}$ is the to-be-learned latent space point corresponding to the particular
combination of categorical variables denoted by $\boldsymbol{t}$. To find
these points in the latent space, LMGP assigns a unique vector (i.e., a prior
representation) to each combination of categorical variables. Then, it uses
matrix multiplication333More complex transformations based on, e.g., NNs, may
also be used, although we do not do so in this paper. to map each of these
vectors to a point in a latent space of dimension $dz$:

$\boldsymbol{z}(\boldsymbol{t})=\boldsymbol{\zeta}(\boldsymbol{t})\boldsymbol{A}\vspace{-0.2cm}$
(6)

where $\boldsymbol{\zeta}(\boldsymbol{t})$ is the
$1\times\sum^{dt}_{i=1}\tau_{i}$ unique prior vector representation of
$\boldsymbol{t}$ and $\boldsymbol{A}$ is a $\sum^{dt}_{i=1}\tau_{i}\times dz$
matrix that maps $\boldsymbol{\zeta}(\boldsymbol{t})$ to
$\boldsymbol{z}(\boldsymbol{t})$. In this paper, we use $d{z}=2$ since it
simplifies visualization and has been shown to provide sufficient flexibility
for learning the latent relations [36]. We construct $\boldsymbol{\zeta}$ via
a form of one-hot encoding where we first construct the $1\times\tau_{i}$
vector
$\boldsymbol{v}^{i}=\left[v^{i}_{1},v^{i}_{2},\dots,v^{i}_{\tau_{i}}\right]$
for each categorical variable $t_{i}$ such that $v^{i}_{j}=1$ when $t_{i}$ is
at level $k=j$ and $v^{i}_{j}=0$ when $t_{i}$ is at level $k\neq j$ for $k\in
1,2,\dots,\tau_{i}$. Then, we set
$\boldsymbol{\zeta}(\boldsymbol{t})=[\boldsymbol{v}^{1},\boldsymbol{v}^{2},\dots,\boldsymbol{v}^{d_{t}}]$.
For example, for the two categorical variables $t_{1}$ and $t_{2}$ with $2$
and $3$ levels, $\boldsymbol{\zeta}(\boldsymbol{t})=[0,1,0,1,0]$ encodes the
combination where both variables are at level $2$.

To train an LMGP, we use maximum likelihood estimation (MLE) to jointly
estimate all of its parameters:

$\left[\hat{m},\hat{s}^{2},\hat{\boldsymbol{\omega}},\hat{\boldsymbol{A}}\right]=\operatorname*{argmax}_{m,s^{2},\boldsymbol{\omega},\boldsymbol{A}}\hskip
14.22636pt\left|2\pi
s^{2}\boldsymbol{R}\right|^{-\frac{1}{2}}\times\exp{\left(-\frac{1}{2}(\boldsymbol{y}-\boldsymbol{1}m)^{T}(s^{2}\boldsymbol{R})^{-1}(\boldsymbol{y}-\boldsymbol{1}m)\right)}\vspace{-0.2cm}$
(7)

where $\left|\cdot\right|$ denotes the determinant operator,
$\boldsymbol{y}=\left[y_{1},\dots,y_{n}\right]^{T}$ is the $n\times 1$ vector
of outputs in the training data, $\boldsymbol{R}$ is the $n\times n$
correlation matrix with the $(i,j)^{\text{th}}$ element
$R_{ij}=r\left((\boldsymbol{x}^{(i)},\boldsymbol{t}^{(i)}),(\boldsymbol{x}^{(j)},\boldsymbol{t}^{(j)})\right)$
for $i,j=1,\dots,n$, and $\boldsymbol{1}$ is a $n\times 1$ vector of ones.

After estimating the hyperparameters, we use the conditional distribution
formulas to predict the response distribution at the arbitrary point
$\boldsymbol{p}^{*}=(\boldsymbol{x}^{*},\boldsymbol{t}^{*})$. The mean and
variance of this normal distribution are:

$\mathop{\mathbb{E}}\left[y\left(\boldsymbol{p}^{*}\right)\right]=\hat{m}+\boldsymbol{r}^{T}\left(\boldsymbol{p}^{*}\right)\boldsymbol{R}^{-1}\left(\boldsymbol{y}-\boldsymbol{1}\hat{m}\right)\vspace{-0.2cm}$
(8)

$\text{cov}\left(y(\boldsymbol{p}^{*}),y(\boldsymbol{p}^{{}^{\prime}})\right)=\hat{s}^{2}r(\boldsymbol{p}^{*},\boldsymbol{p}^{{}^{\prime}})=\hat{s}^{2}\left\\{1-\boldsymbol{r}^{T}(\boldsymbol{p}^{*})\boldsymbol{R}^{-1}\boldsymbol{r}(\boldsymbol{p}^{{}^{\prime}})+g(\boldsymbol{p}^{*})(\boldsymbol{1}^{T}\boldsymbol{R}^{-1}\boldsymbol{1})^{-1}g(\boldsymbol{p}^{{}^{\prime}})\right\\}\vspace{-0.2cm}$
(9)

where $\mathop{\mathbb{E}}$ denotes expectation,
$\boldsymbol{r}\left(\boldsymbol{p}^{*}\right)$ is an $(n\times 1)$ vector
with the $i^{\text{th}}$ element
$r\left(\boldsymbol{p}^{(i)},\boldsymbol{p}^{*}\right)$, and
$g\left(\boldsymbol{p}^{*}\right)=1-\boldsymbol{1}^{T}\boldsymbol{R}^{-1}\boldsymbol{r}\left(\boldsymbol{p}^{*}\right)$.

To perform data fusion via LMGP, we re-frame multi-fidelity modeling as a
manifold learning problem. Assume that we have $ds$ data sources whose inputs
and outputs are denoted by $\boldsymbol{x}^{s_{i}},y^{s_{i}}$, respectively,
with $i=1,\dots,ds$. We first pre-process the data by appending the inputs
with a single categorical variable $t^{s}$ with $ds$ levels (hereafter
referred to as the source index variable) that distinguishes the data sources.
Specifically, we add $t^{s}$ at level $i$ for source $s_{i}$, i.e.,
$\boldsymbol{x}^{s_{i}}\rightarrow\left[\boldsymbol{x}^{s_{i}},\boldsymbol{i}_{n^{s_{i}}\times
1}\right]$, where $\boldsymbol{i}_{n^{s_{i}}\times 1}$ is an $n^{s_{i}}\times
1$ vector of $i$’s and $n^{s_{i}}$ is the number of data points for source
$s_{i}$. We then combine the data for all sources into one unified data set
and fit an LMGP directly it, i.e., we fit LMGP to all of the data from all
sources at once.

The fitted LMGP can provide predictions for any desired data source based on
the level used for $t^{s}$ and as such is an emulator for all of the data
sources. Additionally, since the data sources are distinguished via a
categorical variable, LMGP learns the correlations between them via a
visualizable latent representation and uses these correlations to improve its
predictions [30]. In the case that the raw inputs contain categorical
variables $\boldsymbol{t}^{c}$, we use separate mappings for $t^{s}$ and
$\boldsymbol{t}^{c}$, i.e., we assign unique priors
$\boldsymbol{\zeta}(t^{s})$ and $\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ which
LMGP uses to find mapping matrices $\boldsymbol{A}^{s}$ and
$\boldsymbol{A}^{c}$. The latent points corresponding to each mapping are then
$\boldsymbol{z}^{s}$ and $\boldsymbol{z}^{c}$, respectively.

Note that the correlation function in Equation 5 depends directly on the
euclidean distance between a pair of latent points. This means that relative
distances in the latent space directly correspond to correlations, e.g., if a
pair of data sources $y^{s_{1}}$ and $y^{s_{2}}$ have corresponding latent
points with a distance $\Delta$ in the latent space then this directly implies
by Equation 5 that LMGP has found those two sources to have a correlation of
$\exp{(-\Delta^{2})}$.

###  2.2 Bayesian Neural Networks

Feedforward neural networks (FFNNs) are one of the most common models used in
deep learning and their main goal is to learn the underlying function
$f(\boldsymbol{x})$ that maps the inputs $\boldsymbol{x}$ to the target $y$
[37]. To this end, an FFNN defines the mapping
$\hat{f}(\boldsymbol{x};\boldsymbol{\theta})$ whose parameters
$\boldsymbol{\theta}$ are estimated such that
$\hat{y}=\hat{f}(\boldsymbol{x};\boldsymbol{\theta})$ best approximates
$f(\boldsymbol{x})$. NN-based approaches for MF emulation can provide
attractive advantages since they are universal function approximators [38] and
can handle high-dimensional inputs and large data sets. In this subsection, we
first describe the working principle of FFNNs and motivate the use of BNNs and
Bayes by backprop [39].

FFNNs propagate information from the inputs $\boldsymbol{x}$ to the output $y$
through intermediate computations that define $\hat{f}$. They are
traditionally built via a succession of $L$ layers where $L-2$ hidden layers
are placed between the input and output layers. The output of layer $k$ is
denoted by $\boldsymbol{z}_{k}$ and is obtained as follows:

$\displaystyle\boldsymbol{z}_{1}$ $\displaystyle=\boldsymbol{x},$ (10)
$\displaystyle\boldsymbol{z}_{k}$
$\displaystyle=\phi_{k}\left(\boldsymbol{W}_{k}\boldsymbol{z}_{k-1}+\boldsymbol{b}_{k}\right)\quad\forall\
k\in[2,L-1],$ (11) $\displaystyle\hat{y}$
$\displaystyle=\phi_{L}\left(\boldsymbol{W}_{L}\boldsymbol{z}_{L-1}+{b}_{L}\right)\vspace{-0.2cm}$
(12)

where $\phi$ is the (typically non-linear) activation function. The parameters
$\boldsymbol{\theta}_{k}=(\boldsymbol{W}_{k},\boldsymbol{b}_{k})$, where
$\boldsymbol{W}_{k}$ and $\boldsymbol{b}_{k}$ are the weight matrices and bias
vectors, respectively, correspond to the connections between the $(k-1)^{th}$
and $k^{th}$ layer. For brevity, we denote the parameters of the entire
network by $\boldsymbol{\theta}$.

From a statistical perspective, an FFNN aims to learn the conditional
distribution $P(y|\boldsymbol{x};\boldsymbol{\theta})$ given the noisy data
set $\mathcal{D}$ with independent and identically distributed samples:

$y^{(i)}=f(\boldsymbol{x}^{(i)})+\epsilon\approx\hat{f}(\boldsymbol{x}^{(i)};\boldsymbol{\theta})+\epsilon$
(13)

where $\epsilon\sim\mathcal{N}(0,\sigma^{2})$ represents noise. Equation 13
indicates that
$P(y^{(i)}|\boldsymbol{x}^{(i)})\sim\mathcal{N}(f(\boldsymbol{x}^{(i)}),\sigma^{2})$
and hence the conditional probability
$P\left(\mathcal{D}|\boldsymbol{\theta}\right)$ can be written as:

$\displaystyle P(\mathcal{D}|\boldsymbol{\theta})$
$\displaystyle=\prod_{i=1}^{n}P(y^{(i)}|\boldsymbol{x}^{(i)},\boldsymbol{\theta})P(\boldsymbol{x}^{(i)})=\prod_{i=1}^{n}\mathcal{N}(y^{(i)};\hat{y}^{(i)},\sigma^{2})P(\boldsymbol{x}^{(i)})$
$\displaystyle=\prod_{i=1}^{n}\frac{1}{\sigma\sqrt{2\pi}}\exp{\left(-\frac{1}{2\sigma^{2}}\left(y^{(i)}-\hat{y}^{(i)}\right)^{2}\right)}P(\boldsymbol{x}^{(i)})\vspace{-0.2cm}$
(14)

Since the likelihood function $L(\boldsymbol{\theta})\equiv
P(\mathcal{D}|\boldsymbol{\theta})$, the parameters $\boldsymbol{\theta}$ can
be estimated by maximizing $L(\boldsymbol{\theta})$ (the dependence on
$\boldsymbol{x}$ is dropped for brevity):

$\boldsymbol{\theta}_{MLE}=\arg\max_{\boldsymbol{\theta}}L(\boldsymbol{\theta})\equiv\arg\min_{\boldsymbol{\theta}}\hskip
5.69054pt-\log
P(\mathcal{D}|\boldsymbol{\theta})=\arg\min_{\boldsymbol{\theta}}\frac{1}{n}\sum_{i=1}^{n}\left(y^{(i)}-\hat{f}(\boldsymbol{x}^{(i)};\boldsymbol{\theta})\right)^{2}$
(15)

which is equivalent to minimizing the mean squared error (MSE) of the
predictions $\hat{y}=\hat{f}(\boldsymbol{x};\boldsymbol{\theta})$ with respect
to the targets $y$. Equation 15 can be updated via Bayes rule to consider
prior knowledge on $\boldsymbol{\theta}$ in the optimization. These maximum a
posteriori (MAP) estimates are obtained via:

$\displaystyle\boldsymbol{\theta}_{MAP}$
$\displaystyle=\arg\max_{\boldsymbol{\theta}}P\left(\boldsymbol{\theta}|\mathcal{D}\right)=\arg\max_{\boldsymbol{\theta}}\log
P\left(\boldsymbol{\theta}|\mathcal{D}\right)=\arg\max_{\boldsymbol{\theta}}\log
P\left(\mathcal{D}|\boldsymbol{\theta}\right)+\log
P(\boldsymbol{\theta})\vspace{-0.2cm}$ (16)

where the first term recovers MSE as in Equation 15 and the second term
depends on the prior distribution assigned to the parameters. Equation 16
illustrates that Gaussian and Laplacian priors are equivalent to $L2$ and $L1$
regularization, respectively [37, 39].

FFNNs are likely to overfit in scenarios where data is scarce. Additionally,
they cannot directly quantify prediction uncertainty and are often
overconfident in extrapolation [40]. BNNs are developed to address these
issues [41, 42]. In BNNs, the weights are endowed with probability
distributions (rather than single point estimates) which naturally results in
probabilistic predictions and can dramatically reduce overfitting via
parameter regularization and model averaging.

Predictions via a BNN requires sampling from the posterior distribution of the
parameters, i.e., $P\left(\boldsymbol{\theta}|\mathcal{D}\right)$, which does
not have a closed form and is highly complex. Over the past few years, various
techniques have been developed to obtain samples from
$P\left(\boldsymbol{\theta}|\mathcal{D}\right)$ (or an approximation thereof).
The most popular techniques are based on either Markov Chain Monte Carlo
(MCMC) [43] or variational inference (VI) [44] which, unlike MCMC, learns an
approximation of the posterior distribution.

Although MCMC methods are arguably the best techniques for sampling from the
exact posterior, their lack of scalability makes them inefficient for BNNs of
any practical size [45]. Hence, we employ Bayes by backprop [39] which is a
variational method that approximates
$P\left(\boldsymbol{\theta}|\mathcal{D}\right)$ with the parameterized
distribution $q(\boldsymbol{\theta}|\boldsymbol{\varphi})$ . The parameters
$\boldsymbol{\varphi}$ are learned by minimizing the Kullback–Leibler (KL)
divergence between the true and approximated posteriors:

$\displaystyle\operatorname{KL}[q(\boldsymbol{\theta}|\boldsymbol{\varphi})||P(\boldsymbol{\theta}|\mathcal{D})]=\int
q(\boldsymbol{\theta}|\boldsymbol{\varphi})\log\left(\frac{q(\boldsymbol{\theta}|\boldsymbol{\varphi})}{P(\boldsymbol{\theta}|\mathcal{D})}\right)d\boldsymbol{\theta}=\int
q(\boldsymbol{\theta}|\boldsymbol{\varphi})\log\left(\frac{q(\boldsymbol{\theta}|\boldsymbol{\varphi})P(\mathcal{D})}{P(\mathcal{D}|\boldsymbol{\theta})P(\boldsymbol{\theta})}\right)d\boldsymbol{\theta}$
$\displaystyle=\int q(\boldsymbol{\theta}|\boldsymbol{\varphi})\log
P(\mathcal{D})d\boldsymbol{\theta}+\int
q(\boldsymbol{\theta}|\boldsymbol{\varphi})\log\left(\frac{q(\boldsymbol{\theta}|\boldsymbol{\varphi})}{P(\boldsymbol{\theta})}\right)d\boldsymbol{\theta}-\int
q(\boldsymbol{\theta}|\boldsymbol{\varphi})\log
P(\mathcal{D}|\boldsymbol{\theta})d\boldsymbol{\theta}$ $\displaystyle=\log
P(\mathcal{D})+\operatorname{KL}[q(\boldsymbol{\theta}||\boldsymbol{\varphi})|P(\boldsymbol{\theta})]-\mathbb{E}_{q(\boldsymbol{\theta}|\boldsymbol{\varphi})}[\log
P(\mathcal{D}|\boldsymbol{\theta})]\vspace{-0.2cm}$ (17)

where Bayes rule is applied to $P(\boldsymbol{\theta}|\mathcal{D})$ in the
first line. Then, the parameters $\boldsymbol{\varphi}$ are estimated by
minimizing Equation 17:

$\displaystyle\boldsymbol{\varphi}^{*}=\operatorname*{argmin}_{\boldsymbol{\varphi}}\operatorname{KL}[q(\boldsymbol{\theta}|\boldsymbol{\varphi})||P(\boldsymbol{\theta}|\mathcal{D})]=\operatorname*{argmin}_{\boldsymbol{\varphi}}\operatorname{KL}[q(\boldsymbol{\theta}|\boldsymbol{\varphi})||P(\boldsymbol{\theta})]-\mathbb{E}_{q(\boldsymbol{\theta}|\boldsymbol{\varphi})}[\log
P(\mathcal{D}|\boldsymbol{\theta})]\vspace{-0.2cm}$ (18)

where the term $\log P(\mathcal{D})$ is excluded as it is constant. Equation
18 aims to minimize the sum of two terms. The second term corresponds to the
expectation of the negative log-likelihood while the first term acts as a
regularizer and corresponds to the KL divergence between the approximated
posterior and the prior.

##  3 Probabilistic Neural Data Fusion

Designing a multi-fidelity NN that leverages an ensemble of LF data sets to
better learn an HF source is a very challenging task because of the following
major reasons:

  1. 1.

The relations among the data sources can be unknown. For instance, in the
Rational example (see Table 4 in Appendix A) there are three LF sources whose
biases are not additive. Additionally, these LF sources are not hierarchically
ordered in the sense that the second LF source is more accurate than the first
one.

  2. 2.

There are typically (but not always) more LF data available since LF sources
are generally cheaper compared to the HF source. Learning from such an
unbalanced MF data is quite difficult especially in the presence of scarce HF
data (as an example, see the sample sizes for the engineering applications
described in Appendix B).

  3. 3.

NNs can be built in many ways and, as shown in Section 4, their performance
heavily depends on their architecture and training mechanism. Building an
optimum444We measure optimality in terms of NN’s error in predicting unseen
data from the HF source. NN with small, unbalanced, and MF data is even more
difficult since the sensitivity to the architecture and training mechanism
considerably increases.

We propose to address the above challenges by converting MF modeling to a
manifold555A manifold or a latent-space is a compact representation of a high-
dimensional object such as an image. learning problem which is then solved via
an NN. We design the architecture, loss function, and training mechanism of
this NN with a particular focus on uncertainty sources that include data
scarcity (especially HF samples), noise with unknown variance (which can
affect any of the data sources), non-trivial biases of LF sources, and data
imbalances.

As schematically demonstrated in Figure 1, we convert MF modeling to manifold
learning by augmenting the input space with the categorical variable $t^{s}$
whose levels (e.g., $\\{^{\prime}1^{\prime},^{\prime}2^{\prime},\cdots\\}$ or
$\\{a,b,\cdots\\}$) indicate the source that generates a sample. We then map
this source indicator variable to a low-dimensional manifold via a BNN (see
Block 1 in Figure 1). If the original input space has the categorical
variables $\boldsymbol{t}^{c}$, we similarly map them to a manifold (but this
time we use a deterministic NN, see Block 2 in Figure 1). Afterwards, we
combine the latent variables of these two manifolds with the quantitative
inputs $\boldsymbol{x}$ via a deterministic NN, see Block 3 in Figure 1. As
opposed to the other two blocks, we require Block 3 to produce a normal
probability distribution in order to capture aleatoric uncertainties. Finally,
we train the entire network on the entire666By entire, we mean the combined
data sets from all sources. data using our custom loss function that
noticeably improves the prediction intervals.

Figure 1: Pro-NDF): The proposed architecture allows to combine an arbitrary
number of sources by appending a source indicator variable to the data sets
and then concatenating them. Pro-NDF consists of three blocks that perform
separate tasks related to MF modeling: (1) Block 1 is a BNN that maps a
quantitative prior representation of the source indicator
$\boldsymbol{\zeta}(t^{s})$ to a continuous manifold, (2) Block 2 is an FFNN
that maps a quantitative prior representation of the categorical inputs
$\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ to a continuous manifold, and (3)
Block 3 is an FFNN with a probabilistic output that maps the numerical inputs
and the latent variables to a parametric distribution.

In the following subsections, we elaborate on our rationale for designing a
multi-block architecture and a custom loss function in Section 3.1 and Section
3.2, respectively. Then, we provide some details on the training and inference
stages in Section 3.3.

###  3.1 Multi-Block Architecture

Each block of our network is designed to address particular challenges
associated with MF modeling. Specifically, the BNN of Block 1 maps a
quantitative prior representation $\boldsymbol{\zeta}(t^{s})$ of the source
indicator variable $t^{s}$ to a continuous manifold $\boldsymbol{z}^{s}$. We
design $\boldsymbol{\zeta}(t^{s})$ by one-hot encoding $t^{s}$ to merely
inform the network about the source that generates a sample777If there is some
prior knowledge about the relation among the sources,
$\boldsymbol{\zeta}(t^{s})$ can be designed to reflect it. We do not pursue
designing such informative priors in this work.. We build $\boldsymbol{z}^{s}$
based on a categorical variable because it forces the manifold to uncover the
relations between sources (i.e., the levels of $t^{s}$). These relations are
represented as distances in $\boldsymbol{z}^{s}$ where sources that produce
similar data are encoded with close-by points (see Section 4 for multiple
examples). This distance learning is in sharp contrast to existing approaches
since $(1)$ it does not assume there is any hierarchy between the data
sources, $(2)$ it is scalable to an arbitrary number of data sets, $(3)$ it
enables training the entire network via all available samples, $(4)$ it is
visualizable and interpretable which helps in identifying anomalous data
sources, and $(5)$ it does not assume any specific form (e.g., additive,
multiplicative, etc.) for the biases of LF sources.

Block 1 is the only part of our network where the weights and biases are
endowed with probability distributions. We make this choice to better learn
model form errors and more accurately quantify the epistemic uncertainties due
to lack of data and source-wise discrepancies. We note that, while the outputs
of Block 1 do not parameterize a probability distribution, they are
probabilistic by nature since they are obtained by propagating the
deterministic vector $\boldsymbol{\zeta}(t^{s})$ through some probabilistic
hidden layers.

Block 2 is an FFNN that maps the quantitative prior representation
$\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ of the categorical inputs
$\boldsymbol{t}^{c}$ to the manifold $\boldsymbol{z}^{c}$ (Block 2 is omitted
if the original inputs are purely quantitative). Similar to Block 1, we design
$\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ via one-hot encoding and use
deterministic outputs. However, unlike Block 1 we use a deterministic FFNN in
Block 2 to map $\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ into
$\boldsymbol{z}^{c}$. We make this decision to reduce the number of parameters
and also because the meaning (and hence effects) of categorical inputs across
different sources is typically the same888Due to severe discrepancies such as
large model form errors, the effects of a categorical variable on the response
may be quite different across the sources..

We set the manifold dimension to $2$ for both Block 1 and Block 2, i.e.,
$dz^{s}=dz^{c}=2$. While higher dimensions provide more learning capacities,
our results in Section 4 and those reported elsewhere [46, 47, 48, 49, 50, 51]
indicate that low-dimensional manifolds are quite powerful in learning highly
complex relations. For instance, [52] shows that a single latent variable can
encode smiling in images of human faces which is a high-dimensional and
complex feature in the original data space. Additionally, our choice
simplifies the visualization of the manifolds and reduces the chances of
overfitting since we are primarily interested in scarce data applications.

Block 3 is also an FFNN that maps the numerical inputs and the latent
variables in both manifolds to a parametric distribution which represents the
output. Block 3 has deterministic weights and biases since source-wise
uncertainties are propagated to it via Block 1. However, we equip Block 3 with
a probabilistic output because it: $(1)$ quantifies aleatoric uncertainties
that are inherent to the data sets999The predicted variance also includes
epistemic uncertainties that are propagated from Block 1, see Section 3.3, and
$(2)$ enables designing a multi-task loss that considers the quality of the
prediction intervals (detailed in Section 3.2). Additionally, Block 3 is
responsible for learning the behavior for all data sources simultaneously,
which allows it to leverage correlations between sources to augment
predictions through a process akin to weight sharing.

###  3.2 Uncertainty-Focused Loss

NNs typically provide overconfident predictions especially when they are
trained on small and unbalanced data. As explained in Section 3.1, we aim to
address this issue by making Block 1 and the network’s final output
probabilistic. However, for these measures to work, we must develop an
effective optimization101010Recall that we use Bayes by backprop which takes a
variational approach towards finding the posteriors, see Section 2.2. scheme
where the loss function appropriately rewards prediction intervals (PIs) that
are sufficiently wide (but not too wide) to cover unseen data (especially HF
data). To design such a loss function, we draw inspiration from strictly
proper scoring rules [35] and augment Equation 18 with the negatively oriented
interval score. Our loss is defined as:

$\displaystyle\mathcal{L}=\mathcal{L}_{NLL}+\alpha_{1}\mathcal{L}_{KL}+\alpha_{2}\mathcal{L}_{IS}+\alpha_{3}\mathcal{L}_{2}$
(19)

where $\mathcal{L}_{NLL}$ refers to the negative log-likelihood,
$\mathcal{L}_{KL}$ is the KL divergence between the prior and the variational
posterior distributions on the parameters (only applicable for the BNN from
Block 1), $\mathcal{L}_{IS}$ denotes the interval score term, and
$\mathcal{L}_{2}$ is $L2$ regularization (only applicable for deterministic
NNs, i.e., Block 2 and 3). $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are
hyperparameters that, respectively, determine the relative strengths of
$\mathcal{L}_{KL}$, $\mathcal{L}_{IS}$ and $\mathcal{L}_{2}$ compared to
$\mathcal{L}_{NLL}$. The four terms in Equation 19 are calculated as:

$\displaystyle\mathcal{L}_{NLL}=-\frac{1}{N}\sum^{N}_{i=1}\log\mathcal{N}(y^{(i)};\hat{\mu}^{(i)},\left(\hat{\sigma}^{(i)}\right)^{2})$
(20)
$\displaystyle\mathcal{L}_{KL}=\operatorname{KL}[q(\boldsymbol{\theta}|\boldsymbol{\varphi})||P(\boldsymbol{\theta})]$
(21)
$\displaystyle\mathcal{L}_{IS}=\frac{1}{N}\sum^{N}_{i=1}[(\hat{{{u}}}^{(i)}-\hat{{l}}^{(i)})+\frac{2}{\gamma}(\hat{{l}}^{(i)}-y^{(i)})\mathbbm{1}\\{y^{(i)}<\hat{{l}}^{(i)}\\}+\frac{2}{\gamma}(y^{(i)}-\hat{{{u}}}^{(i)})\mathbbm{1}\\{y^{(i)}>\hat{{{u}}}^{(i)}\\}]$
(22) $\displaystyle\mathcal{L}_{2}=|\boldsymbol{\theta}|^{2}\vspace{-0.1cm}$
(23)

where $\mathcal{L}_{KL}$ is computed via a Monte Carlo approximation, $N$ is
the batch size, and $\mathbbm{1}\\{\cdot\\}$ denotes the indicator function
that returns $1$ if the event in brackets is true and $0$ otherwise. The three
terms of Equation 19 compose a multi-task loss where: $(1)$ the likelihood
term $\mathcal{L}_{NLL}$ penalizes the model if the predicted distribution
does not match the target distribution, $(2)$ the KL divergence term
$\mathcal{L}_{KL}$ favors variational posteriors that are similar to the
assumed prior as per Equation 18, and $(3)$ the interval score term
$\mathcal{L}_{IS}$ rewards narrow PIs while penalizing the model for each
observation $y^{(i)}$ that lies outside the $(1-\gamma)\times 100\%$
prediction interval that spans the range $[\hat{{l}}^{(i)},\hat{{{u}}}^{(i)}]$
where $\hat{{l}}^{(i)}=\hat{\mu}^{(i)}-1.96\hat{\sigma}^{(i)}$ and
$\hat{{{u}}}^{(i)}=\hat{\mu}^{(i)}+1.96\hat{\sigma}^{(i)}$. In this paper, we
use $\gamma=5\%$, thus implying that $\mathcal{L}_{IS}$ is minimized by a
distribution whose $95\%$ PI is as tight as possible while containing all the
training data.

###  3.3 Training and Prediction

In BNNs, the variational posterior of $\boldsymbol{\theta}$ is typically
defined layer-wise as a multivariate Gaussian with mean
$\boldsymbol{\mu}\in\mathbb{R}^{c_{k}}$ and covariance matrix
$\boldsymbol{\Sigma}\in\mathbb{R}^{c_{k}\times c_{k}}$, i.e.,
$\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, where $c_{k}$ is the
total number of connections between two consecutive layers. Estimating the
full covariance matrix requires learning $\mathcal{O}(c_{k}^{2})$ parameters
and is thus computationally prohibitive in most applications [45]. To reduce
the costs, some simplifications have been adopted in the literature, such as
learning diagonal or block diagonal [53] covariance matrices. However, our
approach does not suffer from this computational issue since the only Bayesian
part of our network is Block 1 (see Figure 1) whose size is typically very
small (we use one hidden layer with $5$ neurons for all the studies in Section
4). Hence, we estimate a dense covariance matrix between any two layers of
Block 1 to improve its uncertainty quantification capacity. As for the prior,
we use a zero mean Gaussian distribution with diagonal covariance matrix which
makes the KL term equivalent to $L2$ regularization with a rate defined by the
standard deviation of the prior distribution [54]. Thus, the standard
deviation is a hyperparameter that needs to be tuned specifically to each
problem.

BNNs represent their weights and biases by parameterized distributions which
in our case are multivariate normal with dense covariance matrices. In a
forward pass during either training or prediction, we take individual samples
from these distributions and assign them to the weights and biases. In this
way, instead of explicitly obtaining the true posterior distribution of the
output of Block 1 (i.e., $\boldsymbol{z}^{s}$, see Figure 1), we obtain an
empirical distribution in the $\boldsymbol{z}^{s}$ manifold by taking a number
of forward passes, see Figure 2. We refer to these forward passes as
realizations and as explained below we use different number of passes in
training versus prediction.

Figure 2: Outputs of Pro-NDF : We visualize the outputs of Pro-NDF after it is
trained on the MF data of the HOIP data set, which does not have any numerical
inputs (see Appendix B for more details). To provide probabilistic predictions
that quantify both epistemic as well as aleatoric uncertainties, Pro-NDF
learns a probabilistic fidelity manifold (where sources with similar behavior
are encoded with close-by distributions) and a deterministic manifold for the
categorical inputs.

To obtain the response (in training or testing) at the input $\boldsymbol{u}$
using Pro-NDF, which contains both a BNN component and a probabilistic output,
we use ensemble prediction formulas [34]:

$\hat{\mu}(\boldsymbol{u})=\frac{1}{M}\sum_{j=1}^{M}\hat{\mu}_{\boldsymbol{\theta}_{j}}(\boldsymbol{u})$
(24)
$\hat{\sigma}(\boldsymbol{u})=\frac{1}{M}\sum_{j=1}^{M}\left(\hat{\sigma}^{2}_{\boldsymbol{\theta}_{j}}(\boldsymbol{u})+\hat{\mu}^{2}_{\boldsymbol{\theta}_{j}}(\boldsymbol{u})\right)-\hat{\mu}^{2}(\boldsymbol{u})$
(25)

where $\hat{\mu}_{\boldsymbol{\theta}_{j}}(\boldsymbol{u})$ and
$\hat{\sigma}_{\boldsymbol{\theta}_{j}}(\boldsymbol{u})$ are, respectively,
the mean and standard deviation of the output distribution in the $j^{th}$
realization and $\boldsymbol{\theta}_{j}$ are the associated network
parameters. For predictions with a fitted NN, we use $M=1000$ since it
provides a higher accuracy in quantifying the uncertainty associated with
learning the fidelity manifold (i.e., $\boldsymbol{z}^{s}$). While training
the network, we use $M=200$ to reduce the computational costs.

The performance of an NN is highly sensitive to its architecture and
hyperparameters if the training data is small, unbalanced, and multi-fidelity.
To reduce this sensitivity and leverage the low costs of training a single NN
on small data, we perform automated hyperparameter tuning 111111We use this
approach for all NN-based data fusion approaches (including ours) in Section
4.. To this end, we use RayTune [55] and Hyperopt to find the optimum
hyperparameters and architecture by minimizing the five-fold cross-validation
errors on predicting the high-fidelity data.

For our approach specifically, we apply the above tuning strategy to the
architecture of Block 3, the learning rate of the Adam optimizer,
$\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ in Equation 19, the prior standard
deviation of weight matrices in Block 1, and the batch size. We fix the
architectures of Block 1 and Block 2 to one hidden layer with $5$ neurons and
the dimension of both manifolds to $2$. The activation function for all the
neurons of Block 1 and 3 is hyperbolic tangent, whereas for Block 2 it is the
sigmoid function. For more information and full details on implementation,
please see our GitLab repository.

##  4 Results and Discussions

In this section, we validate our approach on three analytic and two real-world
MF problems (detailed in Appendices A and B) and compare its performance
against LMGP and two other existing NN-based approaches which are based on
simple feedforward networks or sequential multi-fidelity (SMF) networks which
are described in Appendix C. The hyperparameters of all the NN-based
approaches are tuned as described in Section 3.3. We refer the reader to our
GitLab repository for specific details on implementation, estimated
hyperparameters, and training/test data. For LMGP, none of its architectural
parameters (such as the kernel type, mean function, latent map, etc.) are
tuned.

We first conduct an ablation study in Section 4.1 to quantify the impacts of
our designed architecture, loss function, and probabilistic elements. Then, we
test the performance of the four MF approaches on the analytic and real-world
problems in Section 4.2 and Section 4.3, respectively. In each problem, the
goal is to model the HF source as accurately as possible, i.e., to obtain the
lowest mean prediction error while maximizing the number of training/test
samples that fall in the $95\%$ PI. To this end, we use mean squared-error
(MSE) and mean negatively oriented interval score (IS). Note that the FFNN and
SMF approaches are not probabilistic, i.e., they provide point estimates
rather than PIs and therefore they are only evaluated based on MSE.

###  4.1 Ablation Study

To evaluate the impact of the key components of Pro-NDF, we perform an
ablation study on the Rational and DNS-ROM problems which are detailed in
Appendices A and B, respectively. Namely, we analyze the impact of:

  1. 1.

Using a BNN rather than a deterministic FFNN in Block 1 for probabilistically
learning the relations between the data sources.

  2. 2.

Considering $\mathcal{L}_{IS}$ in the loss function of Equation 19.

  3. 3.

Fitting the model to the parameters of a distribution instead of a scalar,
i.e., using a probabilistic output.

  4. 4.

Leveraging the fidelity map to detect the least accurate LF source and, in
turn, assessing whether this source helps emulating the HF source.

Regarding the third item above we note that we no longer use IS in the loss
once the probabilistic output is removed. However, we still calculate the IS
after training based on the empirical distribution of the fidelity manifold
which is produced by the multiple realizations of the BNN component.

Problem |  Model Version | Input data | Components | MSE | IS
---|---|---|---|---|---
HF | LF1 | LF2 | LF3 | $\mathcal{L}_{IS}$ | PB1 | PO
Rational | Base | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | $1.65\times 10^{-3}$ | $0.22$
V1 | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | $3.31\times 10^{-3}$ | $0.26$
V2 | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | $2.83\times 10^{-3}$ | $0.24$
V3 | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✗ | $2.01\times 10^{-3}$ | $0.40$
V4 | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | $5.68\times 10^{-3}$ | $0.98$
DNS-ROM | Base | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | $8.99\times 10^{9}$ | $4.84\times 10^{5}$
V1 | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | $1.23\times 10^{10}$ | $5.89\times 10^{5}$
V2 | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | $2.09\times 10^{10}$ | $1.11\times 10^{6}$
V3 | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✗ | $1.70\times 10^{10}$ | $4.07\times 10^{6}$
V4 | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | $\bm{8.17\times 10^{9}}$ | $5.06\times 10^{5}$
Table 1: Results of the ablation study: We evaluate the effect of removing
individual components of Pro-NDF from it by reporting the MSE and IS on unseen
HF data. All models are trained as discussed in Section 3 (e.g., all models
benefit from automatic hyperparameter tuning). For both MSE and IS, lower
numbers indicate better performance. The ticks indicate whether a component is
used. The acronyms and symbols are defined as: HF: high-fidelity, LF1: low-
fidelity 1, LF2: low-fidelity 2, LF3: low-fidelity 3, $\mathcal{L}_{IS}$:
negatively oriented interval score term in the loss function of Equation 19,
PB1: probabilistic Block 1, PO: probabilistic output.

We summarize the results of the ablation study on the two examples in Table 1.
For both problems, we observe that using all components minimizes the test MSE
and IS. Notably, both of our model’s probabilistic components significantly
increase the performance: the probabilistic output enables Pro-NDF to not only
capture aleatoric uncertainty, but also leverage IS in its loss function.
Additionally, using a BNN improves Pro-NDF’s HF emulation capabilities by
preventing overfitting in scarce data regions (since Block 1 is regularized)
and by partially disentangling epistemic and aleatoric uncertainties which
yields better PIs.

We observe that without a probabilistic output, the IS (and hence the
uncertainty quantification accuracy) drops quite significantly (compare V3 to
V1 and the base in either of the problems) since the model can no longer
account for aleatoric uncertainties. By comparing V1 to the base model in
either of the problems in Table 1 we see that for a model with a probabilistic
output the optimal performance is obtained when $\mathcal{L}_{IS}$ is used in
the loss. That is, leveraging the IS in training improves both mean prediction
and uncertainty quantification (measured via MSE and IS, respectively).

In both problems, evaluating V1 through V3 against one another indicates that
there is a trade-off between MSE and IS. That is, versions that perform well
in terms of MSE, do not generally provide the smallest IS. However, when all
of these components are included in Pro-NDF (see the base model in Table 1 for
either of the problems), both MSE and IS are reduced. This improvement is due
to the fact that the priors and $\mathcal{L}_{IS}$ effectively regularize the
model whose learning capacity is substantially increased by the probabilistic
natures of Block 1 and the output.

The probabilistic fidelity manifold (i.e., output of Block 1) provides an
intuitive and visualizable tools to learn the similarity/discrepancy among the
sources. Hence, once we fit the base model in each problem, we analyze the
learned fidelity manifold to determine the LF source that has the least
similarity to the HF source, see Figure 4(a) and Figure 5(a). Based on the
distances in the fidelity manifold of each problem, we conclude that the third
LF source is the least correlated one with the HF source in both cases. We
exclude this source and its data from MF modeling and refit the base model to
the rest of the data, see version V4 for both problems.

One of the major outputs of Pro-NDF is the learned fidelity manifold which
indicates which LF source has the highest discrepancy compared to the HF
source. Hence, after training a Pro-NDF and inversely identifying the least
accurate LF source, we can build another Pro-NDF while excluding the data from
this source. In the Rational problem, omitting the lowest-fidelity source
results in much worse MSE and IS. We explain this observation by noting that
this problem has an extremely small number of HF samples and therefore it is
important to judiciously use all available data in training. However, in the
DNS-ROM problem version V4 achieves the best MSE while Pro-NDF with all
components achieves the best IS and second best MSE (compare base to V4 in
Table 1). We explain this trend by noting that the size of the training data
in the DNS-ROM problem is significantly higher than that in the Rational
problem. Therefore, omitting a highly inaccurate data source improves mean
prediction accuracy for the HF source in the DNS-ROM problem since the input-
output relationships learned by Block 3 for the different sources are more
similar. Omitting data from this source also increases the ratio of HF data
available in the unified data set which helps in learning the HF behavior.
However, using all data sources provides Pro-NDF with more information which
improves the uncertainty quantification capability and hence a smaller error
on IS.

###  4.2 Analytic Problems

In this section, we validate our approach against LMGP and existing NN-based
technologies for the Rational, Wing-weight, and Borehole examples detailed in
Appendix A. These examples cover a wider range of input dimensionality, number
of sources, and model form errors (e.g., additive and nonlinear biases).
Similar to the previous section, we use MSE and IS on HF test data as the
performance metrics. The input space of these three examples does not have
categorical features and hence both Pro-NDF and LMGP learn a single manifold.
We visualize the fidelity manifolds learned by Pro-NDF and LMGP to examine
these models’ ability in inversely learning the relationships among the data
sources (note that the LF sources are not ordered based on their accuracy). We
highlight that, unlike Pro-NDF, the fidelity manifold of LMGP is not
probabilistic and hence each data source is encoded with a single point in the
manifold.

| Rational | Wing-weight | Borehole
---|---|---|---
Model | MSE | IS | MSE | IS | MSE | IS
Pro-NDF | $\bm{1.65\times 10^{-3}}$ | $0.22$ | $59.14$ | $37.09$ | $14.94$ | $19.24$
LMGP | $1.70\times 10^{-3}$ | $\bm{0.20}$ | $\bm{37.97}$ | $\bm{29.70}$ | $\bm{12.69}$ | $\bm{17.57}$
FFNN | $2.95\times 10^{-3}$ | - | $64.23$ | - | $21.16$ | -
SMF | $8.08\times 10^{-3}$ | - | $542.74$ | - | $172.87$ | -
Table 2: Results on the analytic examples for different models: We test the
performance of Pro-NDF against LMGP and existing NN-based technologies for the
Rational, Wing-weight and Borehole examples detailed in Table 4. The training
procedure for Pro-NDF, LMGP, FFNN and SMF is discussed in Section 3, Section
2.1, Section C.1 and, Section C.2 respectively. We report the MSE and IS on
unseen HF data.

The results for each approach on each problem are summarized in Table 2 and
demonstrate that the probabilistic approaches, i.e., LMGP and Pro-NDF,
significantly outperform the deterministic approaches in all problems. The
FFNN approach performs significantly worse than LMGP and sometimes approaches
the performance of Pro-NDF in MSE, while the SMF approach shows poor
performance for all problems. We explain SMF’s poor performance by noting
that, as explained in Appendix C, hierarchical MF techniques such as SMF
heavily rely on the knowledge of fidelity levels to process the data sources
sequentially in the order of increasing accuracy. Since we assume in the
problem setup that we only know which source has the highest fidelity and do
not know the relative fidelity levels of the LF sources, the LF sources are
ordered sub-optimally in the SMF approach which leads to a very poor
prediction accuracy. The FFNN approach, by contrast, does not rely on the
knowledge of fidelity levels and as such performs better than SMF. However,
its performance lags behind that of LMGP and Pro-NDF because the architecture
is not designed with MF problems in mind.

LMGP, which is considered as our gold standard for MF problems with small
data, outperforms Pro-NDF in both MSE and IS for the Wing-weight and Borehole
problems, and in IS for the Rational problem. The Rational problem is
simultaneously the most data deficient and least complex of the problems
examined in this paper: as shown in Table 4, there are 4 data sources with
only one being especially inaccurate, the input and output are both $1D$, and
there are only $5$ training samples provided for the HF source. Pro-NDF and
LMGP are well suited to tackle this problem as they both perform well for low-
dimensional problems with simple underlying functional forms and well-
correlated sources, and as such they have similar performance. Figure 3(a,b)
reveals that LMGP captures all of the training points in a narrower $95\%$ PI
compared to Pro-NDFwhich explains LMGP’s lower IS in Table 2. However, Pro-NDF
shows a better performance for this problem in terms of mean prediction
accuracy and it also has a higher degree of agreement with the true function
in extrapolation while LMGP reverts to its mean. We therefore conclude that
both methods perform on par on the Rational problem.

Figure 3: High-fidelity emulation on the Rational problem: Pro-NDF and LMGP
approaches produce similar results in terms of mean prediction in
interpolation. However, LMGP has a narrower $95\%$ PI and reverts to its mean
in extrapolation.

The learned fidelity manifold of Pro-NDF for the Rational problem is shown in
Figure 4(a) which indicates that the network has inversely learned the true
relationship between the data sources as $y^{l_{1}}$ and $y^{l_{2}}$ are
encoded close to $y^{h}$ while $y^{l_{3}}$ is quite far from $y^{h}$. These
relative distances are proportional to the accuracy of the LF sources with
respect to the HF source which are reported in Table 4. The fidelity manifold
also shows a high spread in the distributions of the realizations for
individual sources which indicates either a poor fit to the data or a lack of
training samples. In this case, we attribute this spread to lack of data since
the performance in IS and MSE is quite good.

Figure 4: Pro-NDF and LMGP fidelity manifolds for the analytic problems: (a)
Pro-NDF for Rational problem, (b) LMGP for Rational problem, (c) Pro-NDF for
Wing-weight problem, (d) LMGP for Wing-weight problem, (e) Pro-NDF for
Borehole problem, and (f) LMGP for Borehole problem.

The Wing-weight and Borehole problems are both high-dimensional problems with
relatively complex underlying functional forms and small amounts of data. LMGP
is very well suited to tackle this type of problem[30] because the number of
its hyperparameters scales much better than NN-based approaches such as Pro-
NDF . Accordingly, we observe that LMGP achieves lower MSE and IS for both
examples.

Comparing the performance of Pro-NDF across the two high-dimensional problems,
we observe that it performs much better on the Borehole problem. Examining the
fidelity manifold learned by Pro-NDF and LMGP for the Wing-weight problem, see
Figure 4(c) and Figure 4(d), respectively, we see that both approaches
accurately determine the relationship between the sources as they agree with
the RRMSEs reported in Table 4. Specifically, $y^{l_{1}}$ is closer to $y^{h}$
than $y^{l_{2}}$, which in turn is closer than $y^{l_{3}}$. Notably, both LMGP
and Pro-NDF have the same relative ordering and positioning of the sources,
i.e., $(1)$ the mean position of all sources lies on an axis, and $(2)$
$y^{l_{1}}$ is in the opposite direction relative to $y^{h}$ from $y^{l_{2}}$
and $y^{l_{3}}$. This reinforces our earlier assertion in Section 3: the
positions of the sources in the fidelity manifold learned by Pro-NDF reflect
correlations between the data sources. However, the relative distances between
the LF sources in the latent space found by LMGP more accurately represents
the true relationships between the sources because the position for
$y^{l_{3}}$ is much more distant from $y^{h}$ than encoded positions of the
other sources.

In Figure 4(a) we observe a large spread in the realizations (i.e., the
posterior distributions in the fidelity manifold are quite wide) which
partially explains the poor121212Poor with respect to LMGP. The performance of
Pro-NDF is still much better than the other two NN-based approaches.
performance in this problem. We attribute this performance level to the
relative accuracy of the data sources since only one source, $y^{l_{1}}$, is
at all accurate with respect to $y^{h}$ while the other LF sources are quite
inaccurate. LMGP’s performance is not inherently hampered by including poorly
correlated sources in the data fusion problem [30] since its performance, upon
successful optimization, is at worse on par with fitting separate GPs to each
source. By contrast, since Pro-NDF’s Block 3 is responsible for learning the
relations between all sources and uses weight sharing, including especially
inaccurate sources leads to relatively poor performance as shown in 4.1.

The Borehole problem has five total sources where two LF sources ($y^{l_{3}}$
and $y^{l_{4}}$) are accurate while two LF sources ($y^{l_{1}}$ and
$y^{l_{2}}$) are quite inaccurate with respect to $y^{h}$. Since there are
more total data available compared to the Wing-weight problem due to the
additional data source, and since there are more high-accuracy LF sources, we
observe that the spread of the realizations in the fidelity manifold of Pro-
NDF is much smaller than in the Wing-weight, see Figure 4(e). This narrow
spread indicates that a good fit has been achieved. We again observe that the
relative directions and distances of the LF sources from $y^{h}$ are nearly
identical in the manifolds of Pro-NDF and LMGP, see Figure 4(f), and that both
methods have correctly identified the relationships between the sources, see
Table 4. Based on these observations, it is no surprise that Pro-NDF achieves
very good performance and nearly matches LMGP in terms of MSE and IS.

###  4.3 Real-World Problems

In this section, we validate our approach against LMGP and existing NN-based
technologies on two engineering applications which are detailed in Appendix B.
We again use MSE and IS on HF test data as our performance metrics and examine
the manifolds learned by Pro-NDF and LMGP. In both of these applications, the
input space has categorical features (so Pro-NDF and LMGP each build two
manifolds) and we do not know the underlying relationships between the data
sources.

| DNS-ROM | HOIP
---|---|---
Model | MSE | IS | MSE | IS
Pro-NDF | $\bm{8.99\times 10^{9}}$ | $\bm{4.84\times 10^{5}}$ | $\bm{14.16}$ | $\bm{14.84}$
LMGP | $9.66\times 10^{9}$ | $6.60\times 10^{5}$ | $14.33$ | $20.08$
FFNN | $1.12\times 10^{10}$ | - | $21.55$ | -
SMF | $1.81\times 10^{10}$ | - | $28.09$ | -
Table 3: Results on the real-world examples for different models: We test the
performance of Pro-NDF against LMGP and existing NN-based technologies for the
DNS-ROM and HOIP data sets detailed in Appendix B. We report the MSE and IS on
unseen HF data. Figure 5: Pro-NDF and LMGP fidelity manifolds for the real-
world problems: (a) Pro-NDF for DNS-ROM data set, (b) LMGP for DNS-ROM data
set, (c) Pro-NDF for HOIP data set, (d) LMGP for HOIP data set.

The results for each approach on each problem are summarized in Table 3 and
demonstrate that the probabilistic approaches again significantly outperform
the deterministic ones. The FFNN approach performs nearly as well as LMGP and
Pro-NDF in the DNS-ROM problem, but lags behing Pro-NDF and LMGP in the HOIP
problem in terms of MSE. The SMF approach shows poor performance for both
problems for the same reasons provided in Section 4.2. Notably, Pro-NDF
outperforms LMGP for both problems in terms of both metrics which we partially
explain by noting that there are much more data are available in these real-
world problems compared to the analytical examples of Appendix A. Being an NN-
based approach, Pro-NDF scales very well with additional data while the
performance of LMGP has diminishing returns and eventually plateaus (recall
that the latent map and kernel of LMGP are not tuned which contribute to this
plateauing performance).

Figure 6: Pro-NDF categorical manifold for the HOIP problem: The combination
of the categorical variables’ levels are color-coded based on: (a) the levels
of $t^{c}_{1}$, (b) the levels of $t^{c}_{2}$, (c) the levels of $t^{c}_{3}$,
(d) the average output value. Figure 7: LMGP categorical manifold for the HOIP
problem: The combination of the categorical variables’ levels are color-coded
based on: (a) the levels of $t^{c}_{1}$, (b) the levels of $t^{c}_{2}$, (c)
the levels of $t^{c}_{3}$, (d) the average output value.

As shown in Figure 5(a-b), the fidelity manifolds learned by Pro-NDF and LMGP
for the DNS-ROM problem are nearly analogous as the relative distances are
quite similar. However, LMGP finds all sources to be on the diagonal axis
while Pro-NDF learns a more nuanced relationship between the sources, which
may contribute to its superior performance. We also observe that the spreads
in the individual realizations for each point are fairly tight, which
indicates that Pro-NDF is able to learn the relations between the sources
reasonably well and, accordingly, provide good performance in terms of MSE and
IS.

The HOIP problem has three categorical inputs with $10$, $3$, and $16$ levels
and as such Pro-NDF uses two separate latent transformations (one for the data
source and the other for the three categorical variables) that correspond to
Blocks 1 and 2 in Figure 1. The learned categorical manifolds for Pro-NDF and
LMGP are shown in Figures 6 and 7 where the $\boldsymbol{z}^{c}$ is visualized
four times as the combinations of the categorical variables are color-coded
based on the levels of each of the three categorical variables and based on
the average value of the output 131313This average is obtained using the
entire data set including both the training and test data. Since there are no
numerical features, the combined inputs are
$\boldsymbol{u}=\left[\boldsymbol{\zeta}(t^{s}),\boldsymbol{\zeta}(\boldsymbol{t}^{c})\right]$
and $\boldsymbol{\nu}=\left[\boldsymbol{z}^{s},\boldsymbol{z}^{c}\right]$ are
the inputs to Block 3 of Pro-NDF . Recall that we only use a BNN in Block 1
and as such we show only one realization for the manifold for Pro-NDF that
encodes the categorical variables.

Pro-NDF outperforms LMGP in terms of MSE by a small margin and IS by a
significant margin for this problem which we attribute to the size of the data
sets. Pro-NDF is able to leverage these additional data much more readily than
LMGP which only uses simple mapping functions to handle categorical variables
$t^{s}$ and $\boldsymbol{t}^{c}$. Pro-NDF also finds fairly tight spreads in
the probabilistic fidelity manifold, see Figure 5(c); indicating that it has
high certainty in its outputs and that we should expect good performance. We
note that all sources are found to be roughly on one axis and roughly spaced
evenly from each other, which may indicate that Pro-NDF has failed to learn
the more nuanced relationships between the sources. Equally likely, however,
is that the relationships between the sources are simple enough to be
represented in this way; since we do not know the underlying functional forms
for this problem, we cannot give a definitive answer.

We can also glean some information about the relationships between the
categorical variable levels and their impact on the output by examining the
corresponding manifolds in Figures 6 and 7. Figure 6(c) shows that Pro-NDF
finds distinct clusters for all $16$ levels of $t^{c}_{3}$ which indicates
that distinguishing between the levels of $t^{c}_{3}$ is important to learning
the output. Similarly, the levels of $t^{c}_{2}$ are distinguishable in Figure
6(b) as $t^{c}_{2}$ affects the response value. By contrast, Figure 6(a) shows
no apparent trend between the $10$ levels of $t^{c}_{1}$ which implies that
$t^{c}_{1}$ has little effect on the output as Pro-NDF does not learn to
distinguish the levels from each other. By contrast, the manifold found by
LMGP, shown in Figure 7 shows much less distinct clustering for each of the
three categorical variables, which may help explain why it achieves a lower IS
than Pro-NDF . Finally, we examine whether the latent positions for the
categorical combinations are influenced by the average output value in Figure
6(d) and Figure 7(d). The manifold for Pro-NDF shows a clear trend of the
average output value increasing as the latent points move from the bottom-left
of the space to the top-right, while for LMGP there is no obvious trend. Based
on these manifolds, Pro-NDF shows superior ability to discern relationships
between the categorical combinations and between levels of categorical
variables.

##  5 Conclusion

In this paper, we introduce Pro-NDF for data fusion under uncertainty. Pro-NDF
is based on a multi-block NN where each block is designed to take on specific
tasks for MF modeling that arise in typical engineering applications. One of
these blocks is probabilistic whose visualizable output can be used to detect
LF sources with large model form errors. The final output of Pro-NDF is also
probabilistic which enables to not only quantify aleatoric uncertainties, but
also leverage strictly proper scoring rules during training.

We validate each of the key components of Pro-NDF by performing an ablation
study on an analytic and a real-world example. We also demonstrate that Pro-
NDF outperforms other NN-based data fusion approaches by a large margin.
Moreover, Pro-NDF performs on par to LMGP in low-dimensional cases with small
data sets and slightly lags behind LMGP (a competing GP-based approach) in
high-dimensional examples with very small data sets. However, as the size of
the training data increases, Pro-NDF scales better than LMGP and provides
smaller errors. In these studies, we test the performance on unseen HF data
but note that Pro-NDF builds an MF emulator that probabilistically surrogates
all the data sources simultaneously.

A particularly useful output of Pro-NDF is its learnt fidelity manifold which
encodes source-wise similarities/discrepancies. While the learnt distances in
this manifold do not directly link correlation between the sources, we observe
that the fidelity manifold of Pro-NDF and LMGP look quite similar in our
studies. Since the fidelity manifold of LMGP is embedded in its kernel and
hence indicates the correlations, we believe the fidelity of Pro-NDF also
estimates a scaled version of correlation. An added benefit of Pro-NDF’s
fidelity manifold is that it is probabilistic where wide distributions can
indicate if Pro-NDF is able to learn the relation between the data sources.
Reducing this uncertainty via domain knowledge (especially qualitative
information in engineering applications) is a future direction that we plan to
investigate.

The performance of any data fusion approach (including ours) can drop if there
are one or more very inaccurate LF sources. With Pro-NDF , the learned
fidelity manifold can be used to identify-discard such sources and then
retrain Pro-NDF anew. This process can be repeated until all LF sources are
encoded close to the HF source in the fidelity manifold. This iterative
approach is, however, quite inefficient so we plan to develop an automated
mechanism that perhaps leverages the fidelity manifold to adjust the loss
function and, in turn, prevent Pro-NDF from learning the highly inaccurate LF
sources.

##  6 Acknowledgement

We appreciate the support from National Science Foundation (award numbers
OAC-2211908 and OAC-2103708) and the Early Career Faculty grant from NASA’s
Space Technology Research Grants Program (award number 80NSSC21K1809).

## Appendices

We provide the formulations of the analytic problems in Appendix A, the
background and details of the real-world problems in Appendix B, and the
methodology and details of the FFNN and SMF methods in Appendix C.

##  Appendix A Table of Analytic Examples

Name | Source ID | Formulation | $n$ | $\sigma^{2}$ | RRMSE
---|---|---|---|---|---
Rational | $y^{h}(x)$ | $\frac{1}{0.1x^{3}+x^{2}+x+1}$ | $5$ | $0.001$ | \-
$y^{l_{1}}(x)$ | $\frac{1}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.2}x^{3}+x^{2}+x+1}$ | $30$ | $0.001$ |  $0.23$
$y^{l_{2}}(x)$ | $\frac{1}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times x^{3}+x^{2}+x+1}$ | $30$ | $0.001$ |  $0.15$
$y^{l_{3}}(x)$ | $\frac{1}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times x^{3}+x^{2}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times x+1}$ | $30$ | $0.001$ |  $0.73$
Wing Weight | $y^{h}(\boldsymbol{x})$ | $0.036S_{\omega}^{0.758}W_{f\omega}^{0.0035}\left(\frac{A}{\cos^{2}(\Lambda)}\right)^{0.6}q^{0.006}\times$ | $15$ | $25$ |  -
$\lambda^{0.04}\left(\frac{100t_{c}}{\cos(\Lambda)}\right)^{-0.3}+\left(N_{z}W_{dg}\right)^{0.49}+S_{\omega}W_{p}$ |
$y^{l_{1}}(\boldsymbol{x})$ | $0.036S_{\omega}^{0.758}W_{f\omega}^{0.0035}\left(\frac{A}{\cos^{2}(\Lambda)}\right)^{0.6}q^{0.006}\times$ | $50$ | $25$ |  $0.20$
$\lambda^{0.04}\left(\frac{100t_{c}}{\cos(\Lambda)}\right)^{-0.3}+\left(N_{z}W_{dg}\right)^{0.49}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\times W_{p}$ |
$y^{l_{2}}(\boldsymbol{x})$ | $0.036S_{\omega}^{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.8}}W_{f\omega}^{0.0035}\left(\frac{A}{\cos^{2}(\Lambda)}\right)^{0.6}q^{0.006}\times$ | $50$ | $25$ |  $1.14$
$\lambda^{0.04}\left(\frac{100t_{c}}{\cos(\Lambda)}\right)^{-0.3}+\left(N_{z}W_{dg}\right)^{0.49}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\times W_{p}$ |
$y^{l_{3}}(\boldsymbol{x})$ | $0.036S_{\omega}^{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9}}W_{f\omega}^{0.0035}\left(\frac{A}{\cos^{2}(\Lambda)}\right)^{0.6}q^{0.006}\times$ | $50$ | $25$ |  $5.75$
$\lambda^{0.04}\left(\frac{100t_{c}}{\cos(\Lambda)}\right)^{-0.3}+\left(N_{z}W_{dg}\right)^{0.49}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times W_{p}$ |
Borehole | $y^{h}(\boldsymbol{x})$ | $\frac{2\pi T_{u}\left(H_{u}-H_{l}\right)}{\ln{\left(\frac{r}{r_{w}}\right)}\left(1+\frac{2LT_{u}}{\ln{\left(\frac{r}{r_{w}}\right)}r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}}\right)}$ | $15$ | $6.25$ | \-
$y^{l_{1}}(\boldsymbol{x})$ | $\frac{2\pi T_{u}\left(H_{u}-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.8}H_{l}\right)}{\ln{\left(\frac{r}{r_{w}}\right)}\left(1+\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}LT_{u}}{\ln{\left(\frac{r}{r_{w}}\right)}r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}}\right)}$ | $50$ | $6.25$ |  $3.67$
$y^{l_{2}}(\boldsymbol{x})$ | $\frac{2\pi T_{u}\left(H_{u}-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3}H_{l}\right)}{\ln{\left(\frac{r}{r_{w}}\right)}\left(1+\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}8}LT_{u}}{\ln{\left(\frac{r}{r_{w}}\right)}r_{w}^{2}k_{w}}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.75}\frac{T_{u}}{T_{l}}\right)}$ | $50$ | $6.25$ |  $3.73$
$y^{l_{3}}(\boldsymbol{x})$ | $\frac{2\pi T_{u}\left({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1.1}H_{u}-H_{l}\right)}{\ln{\left(\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}4}r}{r_{w}}\right)}\left(1+\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}3}LT_{u}}{\ln{\left(\frac{r}{r_{w}}\right)}r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}}\right)}$ | $50$ | $6.25$ |  $0.38$
$y^{l_{4}}(\boldsymbol{x})$ | $\frac{2\pi T_{u}\left({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1.05}H_{u}-H_{l}\right)}{\ln{\left(\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}r}{r_{w}}\right)}\left(1+\frac{2LT_{u}}{\ln{\left(\frac{r}{r_{w}}\right)}r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}}\right)}$ | $50$ | $6.25$ |  $0.19$
Table 4: Table of analytic functions: The analytic examples have different
input dimensionality, number of sources, and forms of model error. $n$ denotes
the number of samples, $\sigma^{2}$ is the variance of the noise, and RRMSE is
the relative root mean squared error of an LF source with respect to an HF
source, see Equation A-1.

Table 4 details the analytic functions used for the examples covered in
Section 4. For each multi-fidelity problem, we calculate the accuracy of each
LF source with respect to the HF source via relative root mean squared error
(RRMSE):

$\text{RRMSE}=\sqrt{\frac{(\boldsymbol{y}^{l}-\boldsymbol{y}^{h})^{T}(\boldsymbol{y}^{l}-\boldsymbol{y}^{h})}{10000\times\text{var}(\boldsymbol{y}^{h})}}$
(A-1)

where $\boldsymbol{y}^{l}$ and $\boldsymbol{y}^{h}$ are $10000\times 1$ arrays
of outputs sampled randomly via Sobol sequence from the LF and HF sources,
respectively. We use the same sample locations and outputs as our test data
when evaluating MSE and IS in Sections 4.1 and 4.2.

##  Appendix B Background on Real-World Examples

In the DNS-ROM problem, the goal is to predict the toughness of a multiscale
metallic component with spatially varying porosity by combining four sources
of data: $(1)$ high-fidelity: direct numerical simulations (DNS) and $(2,3,4)$
low-fidelity: a reduced-order model (ROM) with three different number of
clusters $(800,1600,3200)$ which balance accuracy against computational costs.
The data sets have six numerical inputs that include pore volume fraction,
number of pores, pore aspect ratio, average nearest neighbor distance among
the pores, evolutionary rate parameter, and critical effective plastic strain
(the last two inputs govern the damage response of the material under load).
The more clusters are used in the ROM, the more similar are its results
compared to those of DNS at the expense of a higher computational burden. The
data set contains $n^{h}=70,n^{l_{1}}=110,n^{l_{2}}=170,n^{l_{3}}=250$
samples. We use $80\%$ of the available samples for each source for training
and $20\%$ for testing. For further details on this data set, we refer the
reader to [2].

In the HOIP problem, the goal is to predict the inter-molecular binding energy
in hybrid organic-inorganic perovskite (HOIP) crystals. The data set has three
categorical inputs with $l_{1}=10$, $l_{2}=3$, and $l_{3}=16$ levels which
correspond to the elements present in each crystal. There are one HF and three
LF data sets with unknown levels of fidelity and
$n_{h}=480,n^{l_{1}}=480,n^{l_{2}}=179,n^{l_{3}}=240$. We use $90\%$ of the
available samples for each source for training and $10\%$ for testing.

##  Appendix C Other Multi-Fidelity NN-Based Approaches

###  C.1 Feedforward Neural Networks

Figure 8: FFNN for multi-fidelity modeling: As Pro-NDF , this approach allows
to use an arbitrary number of sources by appending a source indicator variable
to each data set and concatenating them. The FFNN maps the numerical inputs
$\boldsymbol{x}$, a priori representation of the source indicator
$\boldsymbol{\zeta}(t^{s})$, and categorical inputs
$\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ to the output.

As depicted in Figure 8, for MF modeling via an FFNN we simply feed the
numerical inputs $\boldsymbol{x}$, the prior representation of the source
indicator $\boldsymbol{\zeta}(t^{s})$ and the prior representation of the
categorical inputs $\boldsymbol{\zeta}(\boldsymbol{t}^{c})$ into the FFNN to
produce the output. This approach has two clear disadvantages with respect to
Pro-NDF: $(1)$ it does not provide a tool such as the fidelity manifold of
Pro-NDF that provides a direct visualization of the correlation between the
data sources, and $(2)$ it has a fully deterministic setting which does not
enable uncertainty quantification and thus using a loss function based on
proper scoring rules. In particular, we use the following loss function for
training the FFNN:

$\displaystyle\mathcal{L}=\mathcal{L}_{MSE}+\beta\mathcal{L}_{2}$ (C-2)

where $\mathcal{L}_{MSE}$ is the mean squared error of the predictions and
$\mathcal{L}_{2}$ is $L2$ regularization:

$\displaystyle\mathcal{L}_{MSE}$
$\displaystyle=\frac{1}{N}\sum^{N}_{i=1}(y^{(i)}-\hat{y}^{(i)})^{2}$ (C-3)
$\displaystyle\mathcal{L}_{2}$
$\displaystyle=|\boldsymbol{\theta}|^{2}\vspace{-0.1cm}$ (C-4)

We employ Adam as the optimizer and use RayTune [55] and Hyperopt with five-
fold cross-validation to find the optimum architecture and hyperparameters
which include the learning rate, regularization parameter $\beta$, and batch
size $N$. For further details on implementation, please see our GitLab
repository.

###  C.2 Sequential Multi-Fidelity Networks

Unlike the other methods presented in this paper, multi-fidelity modeling via
SMF requires training a separate sorrgate for each data source. As depicted in
Figure 9, individual FFNNs are trained for each source in the sequence that
ends with the HF source. After a sorrugate is trained for a data source, its
outputs are used to augment the inputs of the next model in the sequence and
hence the resulting input-output relationships are:

Figure 9: Sequential Multi-fidelity (SMF) Networks: SMF is a hierarchical
approach that relies on sequentially training a model (e.g., an FFNN) for each
data source in an ascending order based on the fidelities. The inputs of a
model are augmented with the outputs of the previous one until reaching the
model of the HF source.

$\displaystyle\hat{y}^{s_{1}}$
$\displaystyle=\hat{f}^{s_{1}}\left(\boldsymbol{u}^{s_{1}}\right)$
$\displaystyle\hat{y}^{s_{2}}$
$\displaystyle=\hat{f}^{s_{2}}\left(\boldsymbol{u}^{s_{2}},\hat{y}^{s_{1}}(\boldsymbol{u}^{s_{2}})\right)$
$\displaystyle\cdots$ $\displaystyle\hat{y}^{s_{ds}}$
$\displaystyle=\hat{f}^{s_{ds}}\left(\boldsymbol{u}^{s_{ds}},\hat{y}^{s_{ds-1}}(\boldsymbol{u}^{s_{ds}})\right)$
(C-5)

where $\hat{y}^{s_{i}}$ is the output of the FFNN , $\hat{f}^{s_{i}}$ is the
mapping defined by the FFNN, $\boldsymbol{u}^{s_{i}}$ is the combined numeric
and categorical input
$\boldsymbol{u}=\left[\boldsymbol{x},\boldsymbol{\zeta}(\boldsymbol{t}^{c})\right]$,
and $i$ denotes the data source with $i=ds$ being the HF source. Each
individual FFNN employs the same loss function and optimizer as in the FFNN
method presented in Section C.1.

Unlike the other three MF methods we study in this paper, the SMF approach is
highly sensitive to the ordering of the data sources in the sequence. In the
case that the fidelity levels are known, they are assigned in the order of
increasing fidelity, i.e., source 1 is the least accurate LF source while
source $ds-1$ is the most accurate. With this ordering, the SMF approach
leverages the entire data set to achieve good HF prediction accuracy by
minimizing the complexity of the mapping learned by each successive FFNN.
However, in the case that the fidelities are not known, the order of the LF
sources is assigned randomly. In this case, the mappings of the successive
FFNNs no longer monotonically approaches that of the HF function, and the SMF
approach is unable to properly leverage the additional LF data. In this paper,
we assume that the fidelity levels are unknown and therefore assign the data
source ordering randomly when using SMF.

Similar to the FFNN approach, the SMF approach does not provide a latent
mapping and is entirely deterministic. Like all hierarchical approaches, it
also requires knowledge of fidelity levels for good performance. These factors
lead to a marked disadvantage in the context of the problems examined in this
paper, and we therefore expect the SMF method to perform poorly.

We use RayTune and Hyperopt with five-fold cross-validation to find the
optimum architecture and hyperparameters for each FFNN in the SMF method.
Namely, we tune the learning rate, regularization parameter $\beta$, and batch
size N. We also tune an additional parameter that determines whether to use
the numeric and categorical inputs $\boldsymbol{u}$ in the final FFNN, since
the mapping may be simple enough to learn from just the previous FFNN outputs
in the case that the last LF source is highly accurate. For further details on
implementation, please see our GitLab repository.

## References

  * [1] Ghanshyam Pilania, James E Gubernatis and Turab Lookman  “Multi-fidelity machine learning models for accurate bandgap predictions of solids”  In _Computational Materials Science_ 129 Elsevier, 2017, pp. 156–163
  * [2] Shiguang Deng, Carlos Mora, Diran Apelian and Ramin Bostanabad  “Data-Driven Calibration of Multi-Fidelity Multiscale Fracture Models”  In _arXiv preprint arXiv:2205.12157_ , 2022
  * [3] Xiaotong Liu, Pierre-Paul De Breuck, Linghui Wang and Gian-Marco Rignanese  “A simple denoising approach to exploit multi-fidelity data for machine learning materials properties”  In _arXiv preprint arXiv:2204.10430_ , 2022
  * [4] Souvik Chakraborty, Tanmoy Chatterjee, Rajib Chowdhury and Sondipon Adhikari  “A surrogate based multi-fidelity approach for robust design optimization”  In _Applied Mathematical Modelling_ 47 Elsevier, 2017, pp. 726–744
  * [5] Péter Zénó Korondi, Mariapia Marchi, Lucia Parussini and Carlo Poloni  “Multi-fidelity design optimisation strategy under uncertainty with limited computational budget”  In _Optimization and Engineering_ 22.2 Springer, 2021, pp. 1039–1064
  * [6] Ghina N Absi and Sankaran Mahadevan  “Multi-fidelity approach to dynamics model calibration”  In _Mechanical Systems and Signal Processing_ 68 Elsevier, 2016, pp. 189–206
  * [7] Sanaz Zanjani Foumani, Mehdi Shishehbor, Amin Yousefpour and Ramin Bostanabad  “Multi-Fidelity Cost-Aware Bayesian Optimization”  In _Available at SSRN 4268166_ , 2022
  * [8] Siyu Tao et al.  “Input mapping for model calibration with application to wing aerodynamics”  In _AIAA journal_ 57.7 American Institute of AeronauticsAstronautics, 2019, pp. 2734–2745
  * [9] Slawomir Koziel, Qingsha S Cheng and John W Bandler  “Space mapping”  In _IEEE Microwave Magazine_ 9.6 IEEE, 2008, pp. 105–122
  * [10] John W Bandler et al.  “Space mapping technique for electromagnetic optimization”  In _IEEE Transactions on microwave theory and techniques_ 42.12 IEEE, 1994, pp. 2536–2544
  * [11] Anand Amrit, Leifur Leifsson and Slawomir Koziel  “Fast multi-objective aerodynamic optimization using sequential domain patching and multifidelity models”  In _Journal of Aircraft_ 57.3 American Institute of AeronauticsAstronautics, 2020, pp. 388–398
  * [12] Slawomir Koziel and Leifur Leifsson  “Multi-level CFD-based airfoil shape optimization with automated low-fidelity model selection”  In _Procedia Computer Science_ 18 Elsevier, 2013, pp. 889–898
  * [13] Leifur Leifsson and Slawomir Koziel  “Aerodynamic shape optimization by variable-fidelity computational fluid dynamics models: a review of recent progress”  In _Journal of Computational Science_ 10 Elsevier, 2015, pp. 45–54
  * [14] Marc C Kennedy and Anthony O’Hagan  “Bayesian calibration of computer models”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 63.3 Wiley Online Library, 2001, pp. 425–464
  * [15] John McFarland and Sankaran Mahadevan  “Multivariate significance testing and model calibration under uncertainty”  In _Computer methods in applied mechanics and engineering_ 197.29-32 Elsevier, 2008, pp. 2467–2479
  * [16] Matthew Plumlee  “Bayesian calibration of inexact computer models”  In _Journal of the American Statistical Association_ 112.519 Taylor & Francis, 2017, pp. 1274–1285
  * [17] Dave Higdon et al.  “Combining field data and computer simulations for calibration and prediction”  In _SIAM Journal on Scientific Computing_ 26.2 SIAM, 2004, pp. 448–466
  * [18] Daniel W Apley, Jun Liu and Wei Chen  “Understanding the effects of model uncertainty in robust design with computer experiments”, 2006
  * [19] Maria J Bayarri et al.  “A framework for validation of computer models”  In _Technometrics_ 49.2 Taylor & Francis, 2007, pp. 138–154
  * [20] Paul D Arendt et al.  “Improving identifiability in model calibration using multiple responses”, 2012
  * [21] Paul D Arendt, Daniel W Apley and Wei Chen  “Quantification of model uncertainty: Calibration, model discrepancy, and identifiability”, 2012
  * [22] David A Stainforth et al.  “Uncertainty in predictions of the climate response to rising levels of greenhouse gases”  In _Nature_ 433.7024 Nature Publishing Group, 2005, pp. 403–406
  * [23] Weizhao Zhang et al.  “A numerical Bayesian-calibrated characterization method for multiscale prepreg preforming simulations with tension-shear coupling”  In _Composites Science and Technology_ 170 Elsevier, 2019, pp. 15–24
  * [24] Robert B Gramacy et al.  “Calibrating a large computer experiment simulating radiative shock hydrodynamics”  In _The Annals of Applied Statistics_ 9.3 Institute of Mathematical Statistics, 2015, pp. 1141–1168
  * [25] Lluis Jofre et al.  “Multi-fidelity uncertainty quantification of irradiated particle-laden turbulence”  In _arXiv preprint arXiv:1801.06062_ , 2018
  * [26] Alex A Gorodetsky, John D Jakeman, Gianluca Geraci and Michael S Eldred  “MFNets: multi-fidelity data-driven networks for Bayesian learning and prediction”  In _International Journal for Uncertainty Quantification_ 10.6 Begel House Inc., 2020
  * [27] Rebecca E Morrison, Todd A Oliver and Robert D Moser  “Representing model inadequacy: A stochastic operator approach”  In _SIAM/ASA Journal on Uncertainty Quantification_ 6.2 SIAM, 2018, pp. 457–496
  * [28] Rebecca E Morrison  “Embedded discrepancy operators in reduced models of interacting species”  In _arXiv preprint arXiv:1910.08191_ , 2019
  * [29] Teresa Portone, Damon McDougall and Robert D Moser  “A stochastic operator approach to model inadequacy with applications to contaminant transport”  In _arXiv preprint arXiv:1702.07779_ , 2017
  * [30] Jonathan Tammer Eweis-Labolle, Nicholas Oune and Ramin Bostanabad  “Data Fusion With Latent Map Gaussian Processes”  In _Journal of Mechanical Design_ 144.9 American Society of Mechanical Engineers, 2022, pp. 091703
  * [31] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep learning”  MIT press, 2016
  * [32] Xuhui Meng and George Em Karniadakis  “A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems”  In _Journal of Computational Physics_ 401 Elsevier, 2020, pp. 109020
  * [33] Subhayan De et al.  “On transfer learning of neural networks using bi-fidelity data for uncertainty propagation”  In _International Journal for Uncertainty Quantification_ 10.6 Begel House Inc., 2020
  * [34] Suraj Pawar et al.  “Multi-fidelity information fusion with concatenated neural networks”  In _Scientific Reports_ 12.1, 2022, pp. 5900  DOI: 10.1038/s41598-022-09938-8
  * [35] Tilmann Gneiting and Adrian E Raftery  “Strictly proper scoring rules, prediction, and estimation”  In _Journal of the American statistical Association_ 102.477 Taylor & Francis, 2007, pp. 359–378
  * [36] Nicholas Oune and Ramin Bostanabad  “Latent map Gaussian processes for mixed variable metamodeling”  In _Computer Methods in Applied Mechanics and Engineering_ 387 Elsevier, 2021, pp. 114128
  * [37] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _nature_ 521.7553 Nature Publishing Group, 2015, pp. 436–444
  * [38] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Multilayer feedforward networks are universal approximators”  In _Neural networks_ 2.5 Elsevier, 1989, pp. 359–366
  * [39] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu and Daan Wierstra  “Weight uncertainty in neural network”  In _International conference on machine learning_ , 2015, pp. 1613–1622  PMLR
  * [40] Chuan Guo, Geoff Pleiss, Yu Sun and Kilian Q Weinberger  “On calibration of modern neural networks”  In _International conference on machine learning_ , 2017, pp. 1321–1330  PMLR
  * [41] John Mitros and Brian Mac Namee  “On the validity of Bayesian neural networks for uncertainty estimation”  In _arXiv preprint arXiv:1912.01530_ , 2019
  * [42] Agustinus Kristiadi, Matthias Hein and Philipp Hennig  “Being bayesian, even just a bit, fixes overconfidence in relu networks”  In _International conference on machine learning_ , 2020, pp. 5436–5446  PMLR
  * [43] W Keith Hastings  “Monte Carlo sampling methods using Markov chains and their applications”  Oxford University Press, 1970
  * [44] David M Blei, Alp Kucukelbir and Jon D McAuliffe  “Variational inference: A review for statisticians”  In _Journal of the American statistical Association_ 112.518 Taylor & Francis, 2017, pp. 859–877
  * [45] Laurent Valentin Jospin et al.  “Hands-on Bayesian neural networks—A tutorial for deep learning users”  In _IEEE Computational Intelligence Magazine_ 17.2 IEEE, 2022, pp. 29–48
  * [46] S.. Roweis and L.. Saul  “Nonlinear dimensionality reduction by locally linear embedding”  In _Science_ 290.5500, 2000, pp. 2323–6  DOI: 10.1126/science.290.5500.2323
  * [47] D.. Donoho and C. Grimes  “Hessian eigenmaps: locally linear embedding techniques for high-dimensional data”  In _Proc Natl Acad Sci U S A_ 100.10, 2003, pp. 5591–6  DOI: 10.1073/pnas.1031596100
  * [48] J.. Tenenbaum, V. Silva and J.. Langford  “A global geometric framework for nonlinear dimensionality reduction”  In _Science_ 290.5500, 2000, pp. 2319–23  DOI: 10.1126/science.290.5500.2319
  * [49] Ashutosh Saxena, Abhinav Gupta and Amitabha Mukerjee  “Non-linear dimensionality reduction by locally linear isomaps”  In _Neural Information Processing_ Springer, pp. 1038–1043
  * [50] Ronald R. Coifman and Stéphane Lafon  “Diffusion maps”  In _Applied and Computational Harmonic Analysis_ 21.1, 2006, pp. 5–30  DOI: 10.1016/j.acha.2006.04.006
  * [51] N. Lawrence  “Probabilistic non-linear principal component analysis with Gaussian process latent variable models”  In _Journal of Machine Learning Research_ 6.Nov, 2005, pp. 1783–1816  URL: %3CGo%20to%20ISI%3E://WOS:000236330700002
  * [52] Francois Chollet  “Deep learning with python”  Manning Publications Co., 2017
  * [53] Hippolyt Ritter, Aleksandar Botev and David Barber  “A scalable laplace approximation for neural networks”  In _6th International Conference on Learning Representations, ICLR 2018-Conference Track Proceedings_ 6, 2018  International Conference on Representation Learning
  * [54] Meire Fortunato, Charles Blundell and Oriol Vinyals  “Revisiting Bayes by Backprop”, 2018
  * [55] Richard Liaw et al.  “Tune: A Research Platform for Distributed Model Selection and Training”  In _arXiv preprint arXiv:1807.05118_ , 2018