# Gaussian process regression + deep neural network autoencoder
for probabilistic surrogate modeling in nonlinear mechanics of solids

Saurabh Deshpande  Hussein Rappel  Mark Hobbs  Stéphane P.A. Bordas  Jakub
Lengiewicz<EMAIL_ADDRESS>Department of Engineering; Faculty of Science,
Technology and Medicine; University of Luxembourg  Institute of Fundamental
Technological Research, Polish Academy of Sciences  Department of Engineering,
Faculty of Environment, Science and Economy, University of Exeter, UK  Future
Methods, Rolls-Royce plc, Derby, UK

###### Abstract

Many real-world applications demand accurate and fast predictions, as well as
reliable uncertainty estimates. However, quantifying uncertainty on high-
dimensional predictions is still a severely under-invested problem, especially
when input-output relationships are non-linear. To handle this problem, the
present work introduces an innovative approach that combines autoencoder deep
neural networks with the probabilistic regression capabilities of Gaussian
processes. The autoencoder provides a low-dimensional representation of the
solution space, while the Gaussian process is a Bayesian method that provides
a probabilistic mapping between the low-dimensional inputs and outputs. We
validate the proposed framework for its application to surrogate modeling of
non-linear finite element simulations. Our findings highlight that the
proposed framework is computationally efficient as well as accurate in
predicting non-linear deformations of solid bodies subjected to external
forces, all the while providing insightful uncertainty assessments.

###### keywords:

Surrogate Modeling, Deep Neural Networks, Gaussian Process, Autoencoders,
Uncertainty Quantification, Finite Element Method

††journal: -

##  1 Introduction

In computational mechanics, high-fidelity modeling methods such as the finite
element method (FEM) are standard for achieving detailed insights into the
mechanical behaviors of materials and structures. Despite their capability to
provide accurate predictions for complex non-linear problems, these methods
are often hindered by the substantial computational effort they demand,
thereby constraining their applicability in domains requiring real-time
solutions. As a remedy, various approximate methods, such as surrogate
modeling, have emerged as an indispensable tool to alleviate the computational
burden associated with traditional simulations. By providing the balance
between computational efficiency and the accuracy, these methods play a
crucial role in scenarios requiring fast- or real-time responses.

In conventional methods like FEM, numerical challenges primarily arise from
two main factors: the non-linearities present in the problem and the high
number of unknowns (degrees of freedom) required to attain the desired level
of prediction accuracy. The majority of approximate methods concentrate on
addressing the latter aspect by directly reducing the number of unknowns of
the problem. For instance, the conventional reduced ordered modelling (ROM)
methods accomplish this by linearly projecting higher-order information to
lower-dimensional space using dimensionality reduction techniques such as
Principal Component Analysis (PCA) [1, 2], and Proper Orthogonal Decomposition
(POD) [3, 4]. However, these approaches yield undesirable results when
simulating highly nonlinear phenomena [5, 6, 7]. In response to the linearity
limitation, numerous nonlinear dimensionality reduction techniques have been
proposed [8, 9]. However, obtaining real-time simulations with these methods,
in particular when considering high degree of non-linearity, remains a
challenge [10].

Rapid advancements in machine learning has naturally fueled developments to
approximate modelling. There has been a widespread adoption of machine
learning techniques to accelerate computationally intensive numerical methods
[11, 12, 13, 14, 15]. In particular, deep learning (DL) models have been
successfully used as accurate and computationally inexpensive surrogates to
solve non-linear problems in mechanics [16, 17, 18, 19, 20]. In context of
ROM, DL methods are being extensively used to find _non-linear_ reduced order
representations of the high dimensional data [21]. At the same time, they are
computationally inexpensive at the prediction stage, hence, they are highly
suitable in a wide area of applications requiring real time simulations [22,
23].

This work focuses on another important aspect, which is the capabilities of
surrogate models to quantify uncertainty. Providing reliable uncertainty
estimates alongside predictions is needed in various real-life applications,
such as surgical simulations [24] or autonomous driving [25], to reduce risks
of harmful consequences in these crucial tasks. There exist many machine
learning approaches allowing to quantify model- and data uncertainties. In
deep learning, the most common approach is to use stochastic deep neural
network architectures with relevant suitable training procedures. These are,
for instance, [26, 27, 28]. They are known to be capable of providing
uncertainty information, however, they can be computationally much more
intensive as compared to their deterministic versions. Another class of
approaches, employed specifically in this work, are probabilistic models that
are based on Gaussian process regression [29, 30] (GP regression or GPR). The
excellent predictive power and intrinsic capability of GPRs to quantify
uncertainties make them a perfect candidate for the application of this
contribution. Indeed, GPs have been employed also in mechanics in a wide range
of studies both as surrogate models [31, 32, 33, 34] and as an approach to
model spatially varying parameter fields [35, 36, 37]. Ding et al. [38]
employed a combination of GPR and PCA to model displacement fields in problems
with nonlinear materials (i.e., elastoplastic material), [39] and [40] used
GPR to assess uncertainties in the full-field measurement experimental
context.

Despite being a powerful technique for constructing probabilistic surrogate
models, GPR is known for its poor computational scaling with respect to number
data points [41]. The scalability issue is substantially more significant for
multioutput GPs (MOGP), where surrogate models gives multiple outputs
simultaneously (i.e., $\mathcal{O}(N^{3}M^{3})$, where $N$ and $M$ denote the
number of data points and the dimension of the multiple output, respectively
[42]. Numerous approximated techniques have been developed to significantly
reduce this very poor scaling. The primary effort was to reduce the scaling
with respect to $N$, making the GPR more suitable for big data, see [41] for a
more extensive review. For the multi-output GPs, the most common
simplification is to apply independent GPs for each dimension of the output
(Independent Multi-Output Models), which reduces the scaling to $N^{3}M$.

The framework proposed in this contribution focuses on alleviating part of the
scaling issue by directly reducing the dimension of outputs, $M$. This is
achieved by applying the compressed latent representation of outputs via, so
called, autoencoder deep neural networks—a strategy akin to approaches found
in related works [43, 44]. Autoencoder neural networks, as an unsupervised
learning technique, have been commonly adopted to efficiently learning
compressed representations across various domains [45, 46, 47], also in the
context of computational fluid dynamics [48, 49] and computational solid
mechanics [50, 51]. The primary advantage of autoencoder DNNs is their ability
to reduce the dimensionality of data in a non-linear way [52]. As such, they
perform much better than POD or PCA when high non-linearities are present in
data [53], and are thus used in the present work as an integral component of
the proposed probabilistic framework.

In summary, in this work we combine GPR with autoencoder DNNs to provide
surrogate model for force-displacement predictions in mechanics of solids. In
the proposed framework, we reduce the high-dimensional 2D/3D finite element
solutions (full-field displacement) into their reduced representation by
applying autoencoder neural networks, and then use GPR to provide
probabilistic mapping between the input forces and the reduced-order
displacements. During the prediction phase, for a given unseen force, GPR
predicts the reduced-order (latent space) displacement which is then projected
to the full field using the decoder part of autoencoder network. We
demonstrate the robustness and versatility of our framework by applying it to
synthetic datasets generated from the non-linear finite element simulations.

The remainder of this paper is organized as follows. In Section 2 we introduce
the surrogate modeling framework, and describe its constituent parts. In
Section 3 we present two benchmark examples and study the performance of the
framework. The conclusions and future research directions are outlined in
Section 4.

##  2 Methodology

As outlined in the introduction, in this work, we focus on the problem of
predicting non-linear deformation of solids under external loads. The goal is
to introduce surrogate machine learning framework that is capable of accurate
and fast predictions of full field displacements along with the associated
uncertainties. The core idea is to combine the autoencoder deep neural network
with the Gaussian Process Regression, see Section 2.1 and Section 2.2,
respectively. The training of both models is based on the numerically
generated force-displacement dataset

$\mathcal{D}^{f}=\\{({\bf f}_{i},{\bf u}_{i})\\}_{i=1,..,N}.$ (1)

which is obtained with the classical FEM, see Section 2.5.

Figure 1: Two-stage training. (a) First, the autoencoder neural network is
trained to compress full field displacement data to corresponding latent space
representations. (b) Afterwards, the GP model is trained to provide the
probabilistic force-displacement mapping in the latent space.

The training of the two constituent models is performed sequentially in two
stages, see Figure 1. In the first stage, Figure 1a, the autoencoder network
is trained only on full-field displacement data to get the corresponding
compressed representations, ${\bf u}^{l}$. In the second stage, Figure 1b, the
GP regression model is trained between the forces, ${\bf f}$, and the latent
displacements vector ${\bf u}^{l}$.

Figure 2: For an unseen input force $\mathbf{f}^{*}$, the GP is used to
predict the latent displacement distribution (i.e., our uncertainty about the
displacements is described by a probability distribution). Subsequently, these
latent displacements are projected to the full-field displacement distribution
by using the decoder component of the autoencoder network.

Following the training, the GPR model is stacked with the Decoder part of the
autoencoder network, which provides the probabilistic force-full-displacement
model, see also Section 2.3. Specifically, for a given unseen input force, GPR
model predicts the latent displacement distribution. These latent
displacements are then reconstructed to the full space using the decoder part
of the autoencoder network, see Figure 2.

###  2.1 Obtaining latent representations with autoencoder neural networks

Autoencoders are neural networks that are primarily used for unsupervised
learning and dimensionality reduction tasks [52]. They are designed to encode
high-dimensional input data into a lower-dimensional latent space
representation, and then decode it to the original input space, aiming to
reconstruct the input data as accurately as possible.

The architecture of an autoencoder network typically consists of two main
components: an encoder and a decoder. The encoder takes the input data and
maps it to a compressed representation in the latent space, while the decoder
takes this compressed representation and reconstructs the original input data.
The encoder and decoder are usually symmetrical, meaning that they have the
same number of layers and the layer sizes are mirrored. One of the key aspects
of autoencoders is that the latent space, also referred to as the bottleneck
layer, has a lower dimensionality than the input space. This forces the
network to learn a more efficient and meaningful representation of the data.
By compressing the data into a lower-dimensional space and then reconstructing
it, autoencoders can capture the most salient features by discarding less
relevant information.

The autoencoder network, $\mathcal{U}$, is constructed by subsequently
applying the encoder and decoder networks

$\mathcal{U}({\bf
u},\bm{\theta}_{\text{auto}})=\mathcal{U}_{\text{decoder}}(\mathcal{U}_{\text{encoder}}({\bf
u},\bm{\theta_{\text{auto}}}),\bm{\theta_{\text{auto}}}),$ (2)

where ${\bf u}$ stands for the full field displacements and ${\bf
u}^{l}=\mathcal{U}_{\text{encoder}}({\bf u},\bm{\theta_{\text{auto}}})$ is
their corresponding latent/compressed representation. For a given full field
displacement dataset, $\\{{\bf u}_{i}\\}_{i=1,..,N}$, the autoencoder network
is trained by minimizing the following mean squared error loss

$\mathcal{L}(\mathcal{D}^{f},\bm{\theta}_{\text{auto}})=\frac{1}{N}\sum_{i=1}^{N}\lVert\mathcal{U}({\bf
u}_{i},\bm{\theta}_{\text{auto}})-{\bf u}_{i}\rVert_{2}^{2},$ (3)

where $\bm{\theta}_{\text{auto}}$ are trainable parameters of the autoencoder
network. The optimal parameters, $\bm{\theta}_{\text{auto}}^{*}$, are
retrieved by minimizing the loss function:

$\bm{\theta_{\text{auto}}}^{*}=\text{arg
}\underset{\bm{\theta_{\text{auto}}}}{\text{min}}~{}\mathcal{L}(\mathcal{D}^{f},\bm{\theta_{\text{auto}}}).$
(4)

Therefore, once the autoencoder network is trained, for a given full field
displacement field ${\bf u}$, the corresponding latent representation ${\bf
u}^{l}$ is computed by applying the encoder part of the autoencoder network:

${\bf u}^{l}=\mathcal{U}_{\text{encoder}}({\bf
u},\bm{\theta_{\text{auto}}}^{*})$ (5)

Hence, we obtain latent state dataset $\mathcal{D}^{l}=\\{({\bf f}_{i},{\bf
u}^{l}_{i})\\}_{i=1,..,N}$, from the full field dataset
$\mathcal{D}^{f}=\\{({\bf f}_{i},{\bf u}_{i})\\}_{i=1,..,N}$. Note that forces
are identical for both datasets and they are always provided as smaller
dimensional arrays, which is discussed in detail in Section 3.1. In the next
section, we will use the latent dataset, $\mathcal{D}^{l}$, to elaborate
Gaussian process regression model, which will provide the transformation
between the input forces, ${\bf f}$, and the latent solutions, ${\bf u}^{l}$.

_Remark:_ In our study, we employ two distinct autoencoder architectures, see
Section 3.2.1. The first type utilizes fully connected networks, suitable for
handling arbitrary unstructured meshes. Meanwhile, the second type utilizes
convolutional neural networks (CNNs) specifically designed for structured mesh
scenarios [10, 54].

###  2.2 Gaussian Process Regression in the latent space

As discussed in the introduction, GPs face scalability challenges as the
number of their inputs (data points), $N$, and dimension of outputs, $M$,
increase. In this work we focus on reducing the bottleneck related to $M$. It
is common practice to implement independent GPs for each dimension of the
output, assuming that all output components are independent of each other.
This approach simplifies the complexity of a multi-output GPR problem from
$\mathcal{O}(N^{3}M^{3})$ to $M$ problems with $\mathcal{O}(N^{3})$ complexity
(i.e., $\mathcal{O}(N^{3}M)$), where $N$ is the number of data points and $M$
is the dimensionality of the output. However, for problems with a large number
of degrees of freedom (DOF), such as the examples in this contribution, even
$\mathcal{O}(N^{3}M)$ scaling proves to be computationally costly since $M$
independent GPRs must be implemented for each dimension of the outputs.
Therefore, in this study, we project the high-dimensional displacements (with
dimension $M$) to a significantly lower-dimensional space ($L\ll M$) and then
employ $L$ independent GPs to establish a probabilistic relationship between
input forces and the latent representation of displacements, which scales as
$\mathcal{O}(N^{3}L)$. Because our framework uses independent GPRs for each
latent output, we provide a brief and practical introduction to GPs in this
subsection. Interested readers are referred to [29] for more details.

GPs are an extension of multivariate normal distributions into an infinite-
dimensional Gaussian distribution [55]. Consider two data points $({\bf
f},w({\bf f})),({\bf f}^{\prime},w({\bf f}^{\prime}))$, a GP is completely
specified by its mean function and covariance function. We define mean
function $m({\bf f})$ and the covariance function $k({\bf f},{\bf
f}^{\prime})$ of a real process $w({\bf f})$ as:

$\begin{split}m({\bf f})&=\mathbb{E}(w({\bf f})),\\\ k({\bf f},{\bf
f}^{\prime})&=\mathbb{E}\left((w({\bf f})-m({\bf f}))(w({\bf
f}^{\prime})-m({\bf f}^{\prime}))\right).\end{split}$ (6)

As mentioned, a GP is used to describe a distribution over functions.
Therefore, a realization of a GP reads as:

$w({\bf f})\sim\mathcal{GP}(m({\bf f}),k({\bf f},{\bf f}^{\prime})),$ (7)

where ${\bf f}$ and ${\bf f}^{\prime}$ are two GP input vectors of dimension
$D$, which are force vectors in the scope of the present framework. Hence, in
general, for a finite collection of inputs ${\bf F}=[{\bf f}_{1},\cdots,{\bf
f}_{N}]$:

${\bf w}\sim\mathcal{N}({\bf m},{\bf K}),$ (8)

where ${\bf w}=\begin{bmatrix}w_{1}&\cdots&w_{N}\end{bmatrix}^{T}$, ${\bf
m}=\begin{bmatrix}m({{\bf f}}_{1})&\cdots&m({{\bf f}}_{N})\end{bmatrix}^{T}$
and ${\bf K}$ denotes the covariance matrix between two inputs, i.e., $({\bf
K})_{ij}=k({\bf f}_{i},{\bf f}_{j})$. Frequent choices for the covariance
function are given in [29].

In simple terms, in GPR, we set a GP (Eq.(7)) as a prior for the function that
maps inputs to outputs and then update this prior with available observations
using Bayes’ rule (see [56]). In our framework, the GP predicts displacements
corresponding to latent space representations, as shown in Figures1 and 2.
Therefore, we will use the superscript-$l$ notation, $(\cdot)^{l}$, for all GP
outputs.

Since the proposed framework uses independent GPRs for latent outputs, we will
focus on a single latent output to explain the GP setup. Let
$\tilde{u}^{l}_{i}$ be an observation for a particular latent output
contaminated by Gaussian noise, i.e., $\tilde{u}^{l}_{i}\sim\mathcal{N}(w({\bf
f}_{i}),\sigma^{2})$, with $i=1,\cdots,N$. Without loss of generality, we
assume that the function mapping inputs to the output is a realization of a GP
with zero mean and a covariance function $k({\bf f},{\bf f}^{\prime})$,
defined by its parameter set $\bm{\theta}_{\text{GP}}$. Furthermore, let
$\mbox{\boldmath{$y$}}=[\tilde{u}^{l}_{1},\cdots,\tilde{u}^{l}_{N}]$. The
predictions ${\tilde{u}}^{l*}$ of a GPR for a new data point input ${\bf
f}^{*}$ will be as follows:

${\tilde{u}}^{l*}|\bm{\theta}_{\text{GP}},\sigma^{2},\tilde{{\bf u}}^{l},{\bf
F},{\bf
f}^{*}\sim\mathcal{N}(\mathbb{E}({\tilde{u}}^{l*}),\mathbb{V}({\tilde{u}}^{l*})),$
(9) $\mathbb{E}({\tilde{u}}^{l*})={\bf k}^{*T}({\bf K}+\sigma^{2}{\bf
I}_{N})^{-1}\mbox{\boldmath{$y$}},$ (10)

and

$\mathbb{V}({\tilde{u}}^{l*})=k({\bf f}^{*},{\bf f}^{*})-{\bf k}^{*T}({\bf
K}+\sigma^{2}{\bf I}_{N})^{-1}{\bf k}^{*},$ (11)

where $\mathbb{E}({\tilde{u}}^{l*})$ denotes mean of the predicted
$\tilde{u}^{l*}$ at ${\bf f}^{*}$, $\mathbb{V}({\tilde{u}}^{l*})$ denotes the
prediction variance, ${\bf I}_{N}$ is an $N\times N$ identity matrix, ${\bf
K}$ denotes the $N\times N$ covariance matrix given in Eq. (8), ${\bf k}^{*}$
is an $N\times 1$ vector with $({\bf k}^{*})_{i}=k({\bf f}^{*},{\bf f}_{i})$.
Furthermore, to obtain the optimal parameters for GPR, we train the GP by
maximizing the following equation:

$\log p(\mbox{\boldmath{$y$}}|{\bf
F},\bm{\theta_{\text{GP}}},\sigma^{2})=-\frac{1}{2}\mbox{\boldmath{$y$}}^{T}({\bf
K}+\sigma^{2}{\bf
I}_{n})^{-1}\mbox{\boldmath{$y$}}-\frac{1}{2}\log[\text{det}({\bf
K}+\sigma^{2}{\bf I}_{n})]-\frac{n}{2}\log(2\pi).$ (12)

Note that, $p(\mbox{\boldmath{$y$}}|{\bf F},\bm{\theta},\sigma^{2})$ in Eq.
(12), is the conditional likelihood of observations given the parameters set
$\bm{\theta_{\text{GP}}}$, the input matrix ${\bf F}$, $\sigma^{2}$, hence one
also can use Bayesian inference for identification of parameters by setting
priors on unknowns (i.e., $\bm{\theta_{\text{GP}}}$, and $\sigma^{2}$, see
[56]). Moreover, GP and autoencoder network are trained separately.
Autoencoder network needs to be trained first in order to generate the latent
state dataset, which is then used to train the GP.

###  2.3 Projecting latent GP predictions to the full field space using
decoder

As depicted in the Figure 2, latent space predictions obtained using GP are
projected to full space using decoder part of the autoencoder network. At the
inference stage, for a new input force ${\bf f^{*}}$, GP first generates
latent displacement distribution $p({\tilde{u}}^{*l}|{\bf f^{*}})$ for each
latent output, which is a Gaussian distribution with the mean and variance
obtained using Eq.(10) and Eq.(11) respectively.

$p({\tilde{u}}^{*l}|{\bf
f^{*}})=\mathcal{N}({\bf\mathbb{E}}_{\text{GP}},{\bf\mathbb{V}}_{\text{GP}}),$
(13)

Subsequently, the distributions of full field displacements are obtained by
translating samples from the latent predictive distributions through the
decoder component. Let ${\bf u}_{s}^{l}$ denote the vector of samples drawn
from distributions corresponding to all latent outputs, obtained using
independent GPs for an input ${\bf f}^{*}$. The displacement corresponding to
the full field solution for each degree of freedom (DOF) is a Gaussian
distribution, whose mean and variance are computed as follows:

$\begin{split}(\tilde{u}^{*}_{\mu})_{i}&\approx\frac{1}{S}\sum_{s=1}^{S}(\mathcal{U}_{\text{decoder}}({\bf
u}_{s}^{l}))_{i}\\\
{(u_{\sigma}^{*})_{i}}^{2}&\approx\frac{1}{S-1}\sum_{s=1}^{S}((\mathcal{U}_{\text{decoder}}({\bf
u}_{s}^{l}))_{i}-(u^{*}_{\mu})_{i})^{2}\end{split}$ (14)

where $i$ represents index of DOF. In order to get the full field displacement
distributions, we generate $S$ samples in the latent space, and these samples
are then projected to the original space using the decoder. Now the final
Gaussian distributions of displacements of the full field solution is
represented through the mean ($\tilde{{\bf u}}^{*}_{\mu}$) and standard
deviation (${{\bf u}_{\sigma}^{*}}$) of these $S$ samples. We set $S=300$ for
all the implementations presented in this work.

###  2.4 Understanding of errors

As depicted in the Figure 2, for a new unseen case, the GP is first used to
make displacement prediction in the latent space, which is subsequently
projected to full space solution using the decoder part of the autoencoder
network. Hence the error in the framework arises from GP fitting as well as
the compression and reconstructions steps using the autoencoder network. For
the sake of simplicity, in Figure 3, we demonstrate the schematic for a single
DOF prediction.

Figure 3:  (a) As described in the Figure 2, for an input force ${\bf f}$, the
latent distribution is predicted using GP. $\tilde{u}^{l}_{\mu}$ and $u^{l}$
represent the mean of the latent distribution and true latent solution for a
particular latent DOF, respectively. (b) Latent distribution is then projected
to the full space as described in Section 2.3. $\tilde{u}_{\mu}$, $u_{r}$, $u$
represent the mean of the full field displacement prediction, reconstructed
solution and the true FEM solution for a particular full field DOF,
respectively. $e_{\text{f}}$ represents the framework error while $e_{r}$
represents the reconstruction error, for the particular full field DOF.

The mean prediction of the framework is considered as the displacement
solution obtained using the framework. Hence the absolute difference between
this mean prediction and the true FEM solution is considered as the error of
the framework, $e_{\text{f}}$, as follows:

$e_{\text{f}}=|\tilde{u}_{\mu}-u|$ (15)

The error of reconstruction is computed by taking absolute difference between
the reconstructed solution (represented by) and the FEM solution.
Reconstructed solution is obtained by projecting the true latent solution (as
obtained by Eq. (5)) to the full field, thus neglecting any contribution of
the GP. The reconstruction solution is computed as ${\bf
u}_{r}=\mathcal{U}_{\text{decoder}}({\bf u}^{l})$, and the corresponding
reconstruction error is (for a particular DOF in this case):

$e_{r}=|u_{r}-u|$ (16)

Finally, the contribution of GP fitting error corresponding to the full field
representation is computed by taking the difference between the mean
prediction and the reconstruction solution as follows:

$e_{\text{GP}}=|\tilde{u}_{\mu}-u_{r}|$ (17)

Note: All the errors are computed in the decoded space.

###  2.5 Finite element formulation for non-linear deformations of solid
bodies

We consider the problem of deformation of a body occupying the domain
$\Omega$, that is fixed at a part of its boundary $\Gamma_{\text{u}}$
(Dirichlet boundary conditions), and loaded over the surface
$\Gamma_{\text{t}}$ (Neumann boundary conditions) and/or by body forces
$\bar{b}$. This boundary value problem is expressed by the virtual work
principle

$\int_{\Omega}\mbox{\boldmath{$P$}}(\mbox{\boldmath{$F$}}(\mbox{\boldmath{$u$}}))\cdot\nabla\delta\mbox{\boldmath{$u$}}\,dV-\int_{\Omega}\rho\,\mbox{\boldmath{$\bar{b}$}}\cdot\delta\mbox{\boldmath{$u$}}\,dV-\int_{\Gamma_{t}}\mbox{\boldmath{$\bar{t}$}}\cdot\delta\mbox{\boldmath{$u$}}\,dS=0\quad\quad\forall\delta\mbox{\boldmath{$u$}},$
(18)

where $u$ and $\delta\mbox{\boldmath{$u$}}$ belong to appropriate functional
spaces, $\mbox{\boldmath{$u$}}=\bar{\mbox{\boldmath{$u$}}}$ and
$\delta\mbox{\boldmath{$u$}}={\bf 0}$ on $\Gamma_{\text{u}}$, and
$\mbox{\boldmath{$P$}}(\mbox{\boldmath{$F$}})$ is the first Piola-Kirchhoff
stress tensor. $\mbox{\boldmath{$P$}}(\mbox{\boldmath{$F$}})$ is given by the
constitutive relationship through the strain energy potential,
$W(\mbox{\boldmath{$F$}})$, as:

$\mbox{\boldmath{$P$}}(\mbox{\boldmath{$F$}})=\frac{\partial
W(\mbox{\boldmath{$F$}})}{\partial\mbox{\boldmath{$F$}}},$ (19)

where
$\mbox{\boldmath{$F$}}=\mbox{\boldmath{$I$}}+\nabla\mbox{\boldmath{$u$}}$ is
the deformation gradient tensor. In this work we use following Neo-hookean
(hyperelastic) strain energy density:

$W(\mbox{\boldmath{$F$}})=\frac{\mu}{2}(I_{c}-3-2\ln{J})+\frac{\lambda}{4}(J^{2}-1-2\ln{J}),$
(20)

where $J=\text{det}(\mbox{\boldmath{$F$}})$ and
$I_{\text{c}}=\text{tr}(\mbox{\boldmath{$F$}}^{T}\mbox{\boldmath{$F$}})$.

Following conventional finite element methodology, the problem represented by
Eq. (18) is discretized and posed as a system of nonlinear equations for the
unknown nodal displacements, which is then solved iteratively using Newton-
Raphson method. As the effect, for an input vector of prescribed external
forces, the vector of nodal displacements is computed. The solutions for
different input forces provide the necessary force-displacement dataset, in
which each element consist of the pair $({\bf f},{\bf u})$, where ${\bf f}$
encodes the information on external loading and ${\bf u}$ is the vector of
full-field displacement of FE mesh.

_Remark:_ The vector ${\bf f}$ provides a compressed information about the
external loading. The solid body can be subjected to either loads applied to a
single node, set of nodes, or body forces acting throughout the entire region.
When dealing with point loads, we describe the force using its magnitude
($f_{x}$, $f_{y}$) and the position at which the force is applied (relative to
the fixed boundary). In scenarios involving body forces, we provide the force
density components ($f_{x}$, $f_{y}$) or ($f_{x}$, $f_{y}$, $f_{z}$), for 2D
or 3D problems, respectively. See also more details in Section 3.1.

##  3 Results

###  3.1 Dataset generation information

In this research, we examine two distinct scenarios, as illustrated in Figure
4: a 2D beam and a 3D liver. These scenarios are subjected to varied loading
conditions and are simulated using a non-linear finite element method during
the data generation phase.

2D beam case: In this scenario, the 2D beam is applied with point forces with
random direction and magnitude. These forces are applied to the nodes located
along the top line, which is visually represented by red line in the Figure
4a. For the 2D beam case, displacement data is provided for all degrees of
freedom (DOFs) of the mesh, capturing the full-field displacements of the
beam. However, the forces used as inputs to the surrogate framework are
represented as 3-dimensional arrays. Specifically, the applied point load is
characterized by two components, namely $f_{x}$ and $f_{y}$, representing the
force magnitudes in the x and y directions, respectively. Additionally, the
distance, $d$, from the fixed boundary to the point of application is
considered. This information is visually depicted in Figure 4a.

3D liver case: In this case, the 3D liver model is subjected to body forces
with random directions and magnitudes distributed throughout its volume. The
liver is also constrained at specific nodes located at the right end. This
simulation is aimed at understanding the deformations of the liver under the
influence of these distributed body forces and boundary constraints. Similar
to the beam case, displacement data is provided for all DOFs of the liver
mesh, representing the full-field displacements of the organ. However, in
contrast to the beam case, the forces used as inputs to the surrogate
framework are the body force density components in the x, y, and z directions,
respectively ($b_{x},b_{y},b_{z}$).

Further details, such as the range of external forces applied, as well as
other information used to generate the datasets are given in the Table 1.

Problem |  N.of DOFs ($\mathcal{F})$ |  Range (External forces/ body force density) |  Young’s modulus $E$ [Pa], Poisson’s ratio $\nu$, density $\rho$ [kg/m3] |  Dataset size (train+test)
---|---|---|---|---
2D beam |  128 |  $f_{x},f_{y}$ = -2.5 to 2.5 N |  500, 0.4, - |  $5700+300$
3D liver |  9171 |  $b_{x},b_{y},b_{z}$= -0.02 to 0.02 $\text{N}/\text{kg}$ |  5000, 0.45, 1000 |  $7600+400$
Table 1: Details of FE datasets.

(a)

(b)

Figure 4: Schematics of examples considered in this work (a) 2D beam
discretised with quad elements is applied with point loads on the nodes lying
on red line, this example has been taken from [54]. (b) 3D liver is applied
with random body forces within a prescribed region.

###  3.2 Implementation details

####  3.2.1 Autoencoder networks: Architectures & training details

In this study, we employ two types of autoencoder networks. Firstly, we
utilize a CNN-based architecture, illustrated in Figure 5, tailored for the 2D
beam case. Secondly, we utilize a fully connected autoencoder network,
depicted in Figure 6, implemented for the liver case.

Both of these networks consist of an encoder and a decoder component. The
primary function of these autoencoder networks is to compress the displacement
solutions of the full field space into corresponding latent states. It must be
noted that while the autoencoders handle displacement solutions, the forces
are initially provided in a sparse format and are not processed by the
autoencoder networks.

Architecture for the 2D beam case
Encoder: The input to the CNN autoencoder is a mesh displacement tensor, which
is represented with 2 channels. These channels correspond to the displacements
along the x and y DOFs of the mesh. The input tensor undergoes application of
two convolutional layers, each with 256 channels and $3\times 3$ filters as
shown in Figure 5. Convolutional operations enable feature extraction and
dimensionality reduction. Subsequently, a max-pooling operation is applied,
which reduces spatial dimensions of the tensors while preserving the number of
channels. This process is repeated again with 128 channels, further capturing
important patterns in the data. At the latent level of the network, two
convolution layers with 64 channels are applied, which is followed by
flattening of the tensor. Now the flattened tensor is applied with a dense
layer with $l$ units, which stands for the dimension of the compressed/latent
state. This compressed representation contains crucial information about the
input mesh, capturing the most relevant features.

Decoder: The Decoder takes a flattened tensor of dimension $L$ as its input.
Initially, this input tensor undergoes processing via a dense layer. This
strategic utilization of the dense layer facilitates a transformation of the
arbitrary latent dimensions into a new size that can be conveniently reshaped
into a grid structure, a format that seamlessly aligns with the subsequent
convolutional layer operations. For example with using a dense layer with 640
units as shown in the Figure 5 allows us to reshape the flat tensor to a
$64\times 5\times 2$ tensor, on which series of upsampling and convolution
operations are performed until the original spatial dimensions is achieved. In
the end, a 2 channel convolution with $1\times 1$ filter is applied to get
back the original mesh shape.

Figure 5: CNN autoencoder architecture for used for the 2D beam case. Input is
provided to the network in the structured mesh format, the number at the top
of the convolution layer output indicates the number of channels in the
respective layer. Input is applied series of convolution and
pooling/upsampling layers until the original dimension is retrieved. In the
bottleneck level, fully connected layers (dense layer) are applied and the
number the top of their output represents the number of units present in it.

Convolutional neural networks are inherently designed to work with grid-like
inputs, such as structured meshes. As a result, this methodology isn’t
immediately adaptable to inputs that lack this grid-based nature. To address
this limitation, we introduce an alternative solution in the form of a fully
connected autoencoder network as illustrated in Figure 6. This network can
handle diverse and unstructured meshes.

Architecture for the 3D liver case
Encoder: This component receives input in the form of displacements linked to
all degrees of freedom (DOFs) of the mesh. This information is encoded within
a flattened tensor, which serves as the input layer. It is further applied
with two dense layers with 4096 units each. We also introduce skip connections
as illustrated in the Figure 6, which avoid vanishing gradient issues and
stabilize the training procedure. This procedure applying blocks of dense
layers and skip connections is repeated three more times with 2048, 1024, 512
units for dense layers respectively.

Decoder: Decoder takes the $L=16$ dimensional 1D tensor as the input
(compressed/latent tensor), a similar procedure of applying two dense layers
and skip connections is followed until the original size tensor is obtained.
The number of units used in each dense layer is detailed in the Figure 6.

Figure 6: Fully connected autoencoder architecture used for the liver case.
Skip connections are indicated with the (+) sign. The number of units for each
dense layer output is indicated at the bottom. The latent representation
dimension is denoted as $L$, set to $L=16$ for the 3D liver case.

Training of autoencoder networks
Autoencoder network is trained on the full field displacement dataset, using
the loss function as described by Eq.(3). Minimization is carried out using
the Adam optimizer with recommended parameter configurations as presented by
Kingma [57]. In both instances, a batch size of 16 is employed, with the 2D
beam case trained across 32000 epochs and the 3D liver case for 4000 epochs.
The initial learning rate is set at 1e-4 and linearly decays to 1e-6 during
the training process. Autoencoder networks are trained using Tensorflow [58]
library on a Tesla V100-SXM2 GPU, utilizing the high-performance computing
facilities at the University of Luxembourg [59].

####  3.2.2 Gaussian process details

All the Gaussian Process (GP) implementations presented in this work are
performed using Matérn covariance functions, as described below:

$k_{\text{Mat\'{e}rn}}(r)=\lambda^{2}\Big{(}1+\frac{\sqrt{5}r}{\psi}+\frac{5r^{2}}{3\psi^{2}}\Big{)}\text{exp}\bigg{(}-\frac{\sqrt{5}r}{\psi}\bigg{)},$
(21)

where $r=\lVert\mathbf{f}_{i}-\mathbf{f}_{j}\rVert$, $\mathbf{f}_{i}$ and
$\mathbf{f}_{j}$ denote the $i^{\text{th}}$ and $j^{\text{th}}$ input vectors,
$\lVert\cdot\rVert$ denotes the Euclidean norm, $\psi$ is the length scale,
and $\lambda^{2}$ denotes the variance. Note that we have chosen the
covariance function to avoid overfitting.

The kernel hyperparameters are obtained by training GPs using the scikit-learn
library [60] on the latent datasets, as described in Section 2.2. Optimization
is performed using the L-BFGS-B optimizer [61], with the number of optimizer
restarts chosen as 9 for both cases.

As mentioned earlier, the output of the autoencoder (i.e., the latent
representations) is a multioutput function. Therefore, the GP must map the
input space $\mathbf{f}$ into a vector function $\mathbf{u}^{l}$. Different
methods for such mapping are given in [62]. In this contribution, we model the
Gaussian Process Regression (GPR) without considering any correlations between
the components of the multidimensional output, implying that we can model the
mapping from the inputs to each output component as an independent GP. Note
that the outputs of the autoencoders ($\mathbf{u}^{l}$) do not have any
physical meaning; hence, we can not make any judgment about their inter-
correlation based on physics. However, the results in [37] demonstrate that
independent fields can capture the response of multioutput functions with a
reasonable accuracy.

###  3.3 Predictions of the GP+autoencoder framework

We analyse the test example of 2D beam case, with the maximum nodal
displacement of 1.16 m. Figure 7a shows the displacement field prediction
obtained using the proposed framework which is in great agreement with the FEM
solution as presented in the Figure 7b. Figure 7c provides nodal error of the
entire framework. While Figure 7d and Figure 7e provide reconstruction error
and GP error components as described by Eq.(16) and Eq.(17) respectively.
Error plots show that all the errors are fairly low when compared to the
magnitude of the displacement solution thereby displaying the capability of
the framework in making accurate predictions. Figure 7f provides nodal
uncertainty estimates predicted by the framework (2 standard deviations), as
obtained from Eq.(14). One can observe that the GP component errors are within
the uncertainty bounds for respective nodal predictions, i.e., in other words,
true solution of the framework is contained within the displacement
distribution for all the nodes.

(a) Framework solution (b) FEM solution

(c) Absolute error (d) Reconstruction error

(e) GP component error (f) Uncertainty prediction

Figure 7: Deformation of the 2D beam when applied with the point load on the
top right corner node. Nodal displacements obtained using the (a) proposed
framework (b) FEM respectively. (c) Nodal error and (d) reconstruction error
of the framework. (e) GP component error. (f) Uncertainty predicted by the
framework for the full field solution.

Similarly, we analyse the test example of 3D liver case with the maximum nodal
displacement of 5.74 m. Figure 8 shows that the framework is able to
accurately predict the deformation responses of the 3D liver geometry
subjected to external body forces. Figure 8a, 8b show that the nodal
displacements predicted by the framework are in a great agreement with the FEM
solution. Figure 8c, 8d present the framework and reconstruction errors which
are extremely low compared to the magnitude of displacements. For the maximum
displacement node, the prediction error is $0.2\%$ of its magnitude. Directing
attention to Figure 8e and Figure 8f, these figures present the error
contribution of the GP component and the predicted uncertainty respectively.
Again, the nodal GP errors are within the respective nodal uncertainty values
thereby indicating that the framework is able to make reliable predictions.

Overall in addition to the framework error, its breakdown into reconstruction
and GP error provides a clear and accessible understanding of how each
component contributes to the overall performance of the framework.

(a) Framework solution (b) FEM Solution

(c) Absolute error (d) Reconstruction error

(e) GP component error (f) Uncertainty prediction

Figure 8: Deformation of the 3D liver for the highest nodal displacement case,
subjected to external body force. Nodal displacements obtained using the (a)
proposed framework (b) FEM respectively. (c) Nodal error and (d)
reconstruction error of the framework. (e) GP component error. (f) Uncertainty
predicted by the framework for the full field solution.

###  3.4 Overall accuracy over the test sets

The performance of the framework over the entire test set is computed using
the mean absolute error metric. For the $m^{th}$ test example, the mean error
is given as follows:

$e_{m}=\frac{1}{\mathcal{F}}\sum_{i=1}^{\mathcal{F}}{|\tilde{{\bf
u}}^{*}_{\mu}({\bf f}_{m})^{i}-{\bf u}_{m}^{i}|}.$ (22)

$\mathcal{F}$ is the number of DOFs of the mesh representing the full field
space, $\tilde{{\bf u}}^{*}_{\mu}({\bf f}_{m})$ is the mean prediction of the
GP+Decoder framework (full field prediction) as described through Eq.(14) and
${\bf u}_{m}$ is the finite element solution for the full field space. To have
a single validation metric over the entire test set, we compute the average
mean norm $\bar{e}$ and the corrected sample standard deviation $\sigma(e)$ as
follows:

$\begin{split}\bar{e}=\frac{1}{M}\sum_{m=1}^{M}e_{m},\qquad\sigma(e)=\sqrt{\frac{1}{M-1}\sum_{m=1}^{M}\left(e_{m}-\bar{e}\right)^{2}}.\end{split}$
(23)

Finally, in addition to that, we also use the maximum error per degree of
freedom over the entire test set

$e_{\text{max}}=\max_{m,i}|\tilde{{\bf u}}^{*}_{\mu}({\bf f}_{m})^{i}-{\bf
u}_{m}^{i}|.$ (24)

The proposed framework predicts probabilistic displacement fields
corresponding to the full field mesh. First, we analyse the prediction
performance by comparing mean predictions of the framework to the true FEM
solutions. Table 2 presents error metrics for both benchmark cases.
Additionally, the maximum nodal displacements for both cases are also
provided, i.e., we compute displacements of all the nodes for every single
example, and provide the maximum nodal displacement in the table. Table 2
shows that both mean and maximum errors are fairly low when compared to the
displacement magnitudes. The proposed framework is capable of predicting
mechanical deformation responses with very low errors.

Example | $M$ |  $\bar{e}$ [m] | $\sigma(e)~{}[\text{m}]$ |  $e_{\text{max}}$ [m] | Max nodal disp. [m]
---|---|---|---|---|---
2D Beam | 300 | 0.1 E-3 | 0.8 E-4 | 7.5 E-3 | 1.16
3D liver | 400 | 0.9 E-3 | 0.9 E-3 | 7.7 E-2 | 5.74
Table 2: Error metrics for 2D and 3D test sets for predictions using the
proposed GP+autoencoder framework. $M$ stands for the number of test examples,
and $\bar{e}$, $\sigma(e)$ and $e_{\text{max}}$ are error metrics defined in
Section LABEL:sec:_validation_metric.

Further, we analyze test examples with the maximum nodal displacement for both
benchmark cases. In both cases, we plot nodal error contours representing the
absolute difference between the mean displacement predicted by the framework
and the true FEM solution. We also plot uncertainty estimates obtained using
the proposed framework (see Eq.(14)).

###  3.5 Latent vs. reconstructed uncertainties

GPs are utilized to predict outcomes in the latent space. The GP forecasts a
Gaussian probability distribution for each latent DOF. To ensure the
efficiency and accuracy of GP’s predictions, we provide additional accuracy
measures. First, we define prediction of a particular DOF as _healthy_ , if
the corresponding true solution lies within two standard deviations of the
mean value.

To give an example, for the liver case, we plot latent probability
distributions for a randomly chosen test example and check if the
corresponding true latent values lie within the distribution. These true
latent values are derived from the autoencoder compression of full field
Finite Element Method (FEM) solutions. The displacement variables for all the
latent distributions are transformed to an uniform scale. If the true latent
solution for that latent distribution lies within ($u_{\mu}\pm 2u_{\sigma}$),
it is indicated by the green line, otherwise with the red line as shown in the
Figure 9. For a particular test examples, if all the latent predictions are
_healthy_ , then that prediction is regarded as _correct_ , see, e.g., Figure
9a. Even if a single latent prediction is not _healthy_ , that prediction is
regarded as _incorrect_ , see Figure 9b. This is because a single non-healthy
latent prediction can affect the full field solution during the reconstruction
step.

(a) Correct prediction (b) Incorrect prediction

Figure 9: Latent distribution prediction for two randomly chosen test examples
from the liver test dataset. All the 16 true latent values are transformed to
an uniform scale, along with the displacement random variables. True latents
lying within the two standard deviations ($\pm 2\sigma$ region represented by
shaded gray region) are indicated with green lines, while the latents lying
outside are indicated with red colors. (a) All true latent displacements lie
within $\pm\sigma$, which we classify as a correct prediction. (b) Three true
latent values don’t lie within $\pm 2\sigma$, which we classify as an
incorrect prediction.

The percentage of _healthy_ and correct predictions together contribute to a
comprehensive assessment of the GP’s predictive capabilities in predicting
latent distributions. Table 3 provides these accuracy measures for both
training and test datasets of both cases. The study on test datasets reveals
that even though $\sim 95\%$ of predicted latent components are _healthy_ ,
the percentage of _correct_ GPs predictions can be significantly lower. We
expect this unwanted effect to be more pronounced with the increasing
dimension of the latent space.

Case | % correct predictions | % healthy predictions
---|---|---
Train | Test | Train | Test
2D beam | 95.1 | 77.7 | 99.1 | 94.1
3D liver | 95.6 | 83.7 | 99.2 | 95.9
Table 3: Performance of the GP in correctly predicting latent distributions of
both cases.

###  3.6 Uncertainty for missing data region

As presented in the Table 1, the original 2D beam dataset of 6000 examples was
generated by applying random forces of magnitude -2.5 N to 2.5 N in X and
Y-direction. In order to further test the capability of the framework, we
remove all the cases whose nodal input force magnitude is less than 1 N (the
region represented by white circle in Figure 10). This subset of original data
(5019 examples) is now used to train the proposed framework. Now the framework
trained on this masked dataset is used to make predictions on two new input
datasets.

In the first case, we apply a one dimensional Y-directional force whose
magnitude changes from -3 to 3 N, as indicated by red dot line in the Figure
10a. In Figure 11a, we plot the uncertainty of selected latent DOF predicted
by GP. For comparison, we also plot the error of latent DOF prediction as
well. The figure shows that the uncertainty of prediction for all the latent
DOFs is low in the region supported by the data and the prediction errors are
within uncertainty bounds. The uncertainty grows quickly in both interpolated
and extrapolated missing data region. Especially, the uncertainty grows
rapidly in the extrapolated region (for magnitude of the force greater than
2.5 N).

In the second case, we use a new test dataset in which two dimensional input
forces are randomly generated as denoted by red dots in Figure 10b. In Figure
11b we plot the uncertainty for the first latent DOF. It is evident that the
uncertainty shows an increasing trend in the region not supported by training
data. All the 20 latent DOFs show a very similar trend.

Elevated uncertainty levels are clearly elevated for areas of sparse training
data. In the presented examples, these uncertainties correlate with the
prediction errors, which can provide useful information. However, it’s crucial
to emphasize that the predictive capabilities in the data-sparse and
extrapolated regions can never be assured, with the tendency to deteriorate
with the distance from the training data.

(a)

(b)

Figure 10: Reduced training dataset. The subset of the original dataset is
marked with gray color. The distributions of two new test sets are illustrated
with red dots in (a) and (b).

(a)

(b)

Figure 11: Uncertainty prediction for 2D beam case, for the GP trained using
the missing dataset as indicated in Figure 10. In (a), the uncertainties (red
line) and true latent prediction errors (blue line) are shown for selected
four latent components (component number is indicated in the inset) for the
test data depicted in Figure 10a. In (b), the uncertainties are shown for the
first latent component for the test data depicted in Figure 10b. In (a) and
(b) we observe the trend of increased uncertainty in the region not supported
by the training data.

##  4 Conclusions

This work presents a novel approach for probabilistic surrogate modeling of
high-fidelity simulations. It combines Gaussian processes with reduced-order
modeling to efficiently simulate the mechanics of solids while accounting for
uncertainties. The reduced representations of high-dimensional displacement
data are achieved using autoencoder neural networks, allowing for drastic non-
linear compression. For instance, this method reduced approximately 9100
degrees of freedom (DOFs) to 16 latent-space components, as demonstrated with
the liver example. This strategy significantly reduces the computational cost
of probabilistic predictions using Gaussian processes by lowering the output
dimensionality, which is crucial because this cost can scale linearly to
cubically with the output dimension size. We demonstrate that the proposed
framework can accurately predict non-linear hyper-elastic deformations of
solid bodies, along with the associated model uncertainties.

In addition to demonstrating the capabilities of our probabilistic framework,
we also studied its possible limitations. We emphasized that while the
uncertainties in the missing data region provide a useful indicator that is
correlated with the distance from training data, they do not necessarily need
to represent the true solution errors. The development of methods to improve
the quantitative nature of those predictions remains an open topic for future
research. We also pointed out that even if a single component of the latent
solution is predicted incorrectly, it can impact the entire reconstructed
solution. This led us to propose more meaningful error metrics for quantifying
the performance of the proposed probabilistic framework. This issue can be
alleviated by developing a suited training strategy, which poses an
interesting future research direction.

This study opens up several other opportunities for future research. One
exciting direction is applying this approach to probabilistic simulations of
time and path dependent problems. One could use latent state to track how
solutions evolve in such cases, which would make the whole process much faster
and more efficient [63]. Another direction would be to exploit state of the
art graph autoencoder networks to find compressed representations for
arbitrary high-dimensional meshes [64]. This advancement would overcome the
limitations of fully connected networks, which are recognized for their
challenges in handling problems with a high number of dimensions. Finally, we
believe that the proposed approach can find direct application in research and
development, which is not restricted to problems in mechanics but covers a
broader spectrum of engineering and scientific domains.

## Contributions

Saurabh Deshpande: Conceptualization, Data curation, Formal analysis,
Methodology, GP+ANN implementation, Validation, Visualization, Writing –
original draft, Writing – review & editing.
Hussein Rappel: Formal analysis, Methodology, GP implementation, Writing –
original draft, Writing – review & editing.
Mark Hobbs: Reviewing the manuscript.
Stéphane P.A. Bordas: Funding acquisition, Reviewing the manuscript,
Supervision.
Jakub Lengiewicz: Conceptualization, Data curation, Investigation,
Methodology, Project administration, Supervision, Validation, Visualization,
Writing – original draft, Writing – review & editing.

## References

  * [1] Lorenzo Grassi, Enrico Schileo, Christelle Boichon, Marco Viceconti, and Fulvia Taddei.  Comprehensive evaluation of pca-based finite element modelling of the human femur.  Medical engineering & physics, 36, 08 2014.
  * [2] Gianmarco Aversano, Aurélie Bellemans, Zhiyi Li, Axel Coussement, Olivier Gicquel, and Alessandro Parente.  Application of reduced-order models based on pca & kriging for the development of digital twins of reacting flow applications.  Computers & Chemical Engineering, 121:422–441, 2019.
  * [3] O. Goury and C. Duriez.  Fast, generic and reliable control and simulation of soft robots using model order reduction.  IEEE Transactions on Robotics, 34(6):1565 – 1576, 2018.
  * [4] Paris Papavasileiou, Eleni D. Koronaki, Gabriele Pozzetti, Martin Kathrein, Christoph Czettl, Andreas G. Boudouvis, and Stéphane P.A. Bordas.  Equation-based and data-driven modeling strategies for industrial coating processes.  Computers in Industry, 149:103938, 2023.
  * [5] Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller.  Nonlinear component analysis as a kernel eigenvalue problem.  Neural computation, 10(5):1299–1319, 1998.
  * [6] Gaëtan Kerschen, J.-C Golinval, ALEXANDER VAKAKIS, and LAWRENCE BERGMAN.  The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview.  Nonlinear Dynamics, 41:147–169, 08 2005.
  * [7] P. Kerfriden, P. Gosselet, S. Adhikari, and S.P.A. Bordas.  Bridging proper orthogonal decomposition methods and augmented newton–krylov algorithms: An adaptive model order reduction for highly nonlinear mechanical problems.  Computer Methods in Applied Mechanics and Engineering, 200(5):850–866, 2011.
  * [8] Miguel A Mendez.  Linear and nonlinear dimensionality reduction from fluid mechanics to machine learning.  Measurement Science and Technology, 34(4):042001, jan 2023.
  * [9] Pedro Díez, Alba Muixí, Sergio Zlotnik, and Alberto García-González.  Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition.  International Journal for Numerical Methods in Engineering, 122(24):7306–7327, 2021.
  * [10] A. Mendizabal, P. Márquez-Neila, and S. Cotin.  Simulation of hyperelastic materials in real-time using deep learning.  Medical Image Analysis, 59:101569, 10 2019.
  * [11] Michael Frank, Dimitris Drikakis, and Vassilis Charissis.  Machine-learning methods for computational science and engineering.  Computation, 8(1), 2020.
  * [12] Steven L. Brunton, Bernd R. Noack, and Petros Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics, 52(Volume 52, 2020):477–508, 2020.
  * [13] Anina Šarkić Glumac, Onkar Jadhav, Vladimir Despotović, Bert Blocken, and Stephane P.A. Bordas.  A multi-fidelity wind surface pressure assessment via machine learning: A high-rise building case.  Building and Environment, 234:110135, 2023.
  * [14] Saurabh DESHPANDE.  Data Driven Surrogate Frameworks for Computational Mechanics: Bayesian and Geometric Deep Learning Approaches.  PhD thesis, Unilu - University of Luxembourg [Faculty of Science, Technology and Medicine], Luxembourg, 2023.
  * [15] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia.  Learning skillful medium-range global weather forecasting.  Science, 382(6677):1416–1421, 2023.
  * [16] Saurabh Deshpande, Raúl I. Sosa, Stéphane P. A. Bordas, and Jakub Lengiewicz.  Convolution, aggregation and attention based deep neural networks for accelerating simulations in mechanics.  Frontiers in Materials, 10, 2023.
  * [17] V. Krokos, V. Bui Xuan, S. P. A. Bordas, P. Young, and P. Kerfriden.  A bayesian multiscale cnn framework to predict local stress fields in structures with microscale features.  Computational Mechanics, 69:733–766, 05 2022.
  * [18] Vasilis Krokos, Stéphane P. A. Bordas, and Pierre Kerfriden.  A graph-based probabilistic geometric deep learning framework with online physics-based corrections to predict the criticality of defects in porous materials, 2022.
  * [19] Saurabh Deshpande, Stéphane P.A. Bordas, and Jakub Lengiewicz.  Magnet: A graph u-net architecture for mesh-based simulations.  Engineering Applications of Artificial Intelligence, 133:108055, 2024.
  * [20] Karol Frydrych, Maciej Tomczak, and Stefanos Papanikolaou.  Crystal plasticity parameter optimization in cyclically deformed electrodeposited copper—a machine learning approach.  Materials, 17(14), 2024.
  * [21] Yasi Wang, Hongxun Yao, and Sicheng Zhao.  Auto-encoder based dimensionality reduction.  Neurocomputing, 184:232–242, 2016.
  * [22] G Thippa Reddy, M Praveen Kumar Reddy, Kuruva Lakshmanna, Rajesh Kaluri, Dharmendra Singh Rajput, Gautam Srivastava, and Thar Baker.  Analysis of dimensionality reduction techniques on big data.  Ieee Access, 8:54776–54788, 2020.
  * [23] Stefania Fresca, Luca Dede’, and Andrea Manzoni.  A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes.  Journal of Scientific Computing, 87:1–36, 2021.
  * [24] Arnaud Mazier, Sophie Ribes, Benjamin Gilles, and Stéphane P.A. Bordas.  A rigged model of the breast for preoperative surgical planning.  Journal of Biomechanics, 128:110645, 2021.
  * [25] Antonio Loquercio, Mattia Segu, and Davide Scaramuzza.  A general framework for uncertainty estimation in deep learning.  IEEE Robotics and Automation Letters, 5(2):3153–3160, 2020.
  * [26] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra.  Weight uncertainty in neural network.  In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1613–1622, Lille, France, 07–09 Jul 2015. PMLR.
  * [27] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * [28] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes, 2022.
  * [29] C. E. Rasmussen and C. K. I. Williams.  Gaussian processes for machine learning, volume 1.  MIT press Cambridge, 2006.
  * [30] S. Rogers and M. Girolami.  A first course in machine learning, second edition.  Chapman & Hall/CRC, 2nd edition, 2016.
  * [31] M.C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.
  * [32] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C. H. Lin, and J. Tu.  A framework for validation of computer models.  Technometrics, 49(2):138–154, 2007.
  * [33] P.D. Arendt, D.W. Apley, and W. Chen.  Quantification of model uncertainty: Calibration, model discrepancy, and identifiability.  Journal of Mechanical Design, 134(10):100908, 2012.
  * [34] Arnaud Vadeboncoeur, Ieva Kazlauskaite, Yanni Papandreou, Fehmi Cirak, Mark Girolami, and Omer Deniz Akyildiz.  Random grid neural processes for parametric partial differential equations.  In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 34759–34778. PMLR, 23–29 Jul 2023.
  * [35] P. S. Koutsourelakis.  A multi-resolution, non-parametric, Bayesian framework for identification of spatially-varying model parameters.  Journal of Computational Physics, 228(17):6184–6211, 2009.
  * [36] H. Rappel, L. Wu, L. Noels, and L. A. A. Beex.  A Bayesian Framework to Identify Random Parameter Fields Based on the Copula Theorem and Gaussian Fields: Application to Polycrystalline Materials.  Journal of Applied Mechanics, 86(12):121009, 10 2019.
  * [37] Hussein Rappel, Mark Girolami, and Lars A.A. Beex.  Intercorrelated random fields with bounds and the bayesian identification of their parameters: Application to linear elastic struts and fibers.  International Journal for Numerical Methods in Engineering, 123(15):3418–3463, 2022.
  * [38] Chensen Ding, Hussein Rappel, and Tim Dodwell.  Full-field order-reduced gaussian process emulators for nonlinear probabilistic mechanics.  Computer Methods in Applied Mechanics and Engineering, 405:115855, 2023.
  * [39] Carl Jidling.  Strain field modelling using gaussian processes.  Report: Uppsala University, Division of Systems and Control, (17004):82, 2017.
  * [40] Dario Poloni, Daniele Oboe, Claudio Sbarufatti, and Marco Giglio.  Towards a stochastic inverse finite element method: A gaussian process strain extrapolation.  Mechanical Systems and Signal Processing, 189:110056, 2023.
  * [41] Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai.  When gaussian process meets big data: A review of scalable gps.  IEEE transactions on neural networks and learning systems, 31(11):4405–4423, 2020.
  * [42] Wessel Bruinsma, Eric Perim, William Tebbutt, Scott Hosking, Arno Solin, and Richard Turner.  Scalable exact inference in multi-output Gaussian processes.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1190–1201. PMLR, 13–18 Jul 2020.
  * [43] Romit Maulik, Themistoklis Botsas, Nesar Ramachandra, Lachlan R. Mason, and Indranil Pan.  Latent-space time evolution of non-intrusive reduced-order models using gaussian process emulation.  Physica D: Nonlinear Phenomena, 416:132797, 2021.
  * [44] James Donnelly, Alireza Daneshkhah, and Soroush Abolfathi.  Forecasting global climate drivers using gaussian processes and convolutional autoencoders.  Engineering Applications of Artificial Intelligence, 128:107536, 2024.
  * [45] Jonathan Masci, Ueli Meier, Dan Cireşan, and Jürgen Schmidhuber.  Stacked convolutional auto-encoders for hierarchical feature extraction.  In Artificial Neural Networks and Machine Learning–ICANN 2011: 21st International Conference on Artificial Neural Networks, Espoo, Finland, June 14-17, 2011, Proceedings, Part I 21, pages 52–59. Springer, 2011.
  * [46] Cheng-Yuan Liou, Wei-Chen Cheng, Jiun-Wei Liou, and Daw-Ran Liou.  Autoencoder for words.  Neurocomputing, 139:84–96, 2014.
  * [47] Zilong Wang and Young-Jin Cha.  Unsupervised deep learning approach using a deep auto-encoder with a one-class support vector machine to detect damage.  Structural Health Monitoring, 20(1):406–425, 2021.
  * [48] Pranshu Pant, Ruchit Doshi, Pranav Bahl, and Amir Barati Farimani.  Deep learning for reduced order modelling and efficient temporal evolution of fluid simulations.  Physics of Fluids, 33(10), 2021.
  * [49] Takaaki Murata, Kai Fukami, and Koji Fukagata.  Nonlinear mode decomposition with convolutional neural networks for fluid dynamics.  Journal of Fluid Mechanics, 882:A13, 2020.
  * [50] Stefania Fresca, Giorgio Gobat, Patrick Fedeli, Attilio Frangi, and Andrea Manzoni.  Deep learning-based reduced order models for the real-time simulation of the nonlinear dynamics of microstructures.  International Journal for Numerical Methods in Engineering, 123(20):4749–4777, 2022.
  * [51] Krushna Shinde, Vincent Itier, José Mennesson, Dmytro Vasiukov, and Modesar Shakoor.  Dimensionality reduction through convolutional autoencoders for fracture patterns prediction.  Applied Mathematical Modelling, 114:94–113, 2023.
  * [52] G.E. Hinton and R.R. Salakhutdinov.  Reducing the dimensionality of data with neural networks.  Science (New York, N.Y.), 313:504–7, 08 2006.
  * [53] Jasem Almotiri, Khaled Elleithy, and Abdelrahman Elleithy.  Comparison of autoencoder and principal component analysis followed by neural network for e-learning using handwritten recognition.  In 2017 IEEE Long Island Systems, Applications and Technology Conference (LISAT), pages 1–5, 2017.
  * [54] Saurabh Deshpande, Jakub Lengiewicz, and Stéphane P.A. Bordas.  Probabilistic deep learning for real-time large deformation simulations.  Computer Methods in Applied Mechanics and Engineering, 398:115307, 2022.
  * [55] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  Bayesian data analysis.  Chapman & Hall/CRC Texts in Statistical Science. Chapman & Hall/CRC, 2003.
  * [56] H. Rappel, L. A. A. Beex, J. S. Hale, L. Noels, and S. P. A Bordas.  A tutorial on bayesian inference to identify material parameters in solid mechanics.  Archives of Computational Methods in Engineering, 27, 2020.
  * [57] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2017.
  * [58] J. V. Dillon, I. Langmore, D. Tran, E. Brevdo, S. Vasudevan, D. Moore, B. Patton, A. Alemi, M. D. Hoffman, and R. A. Saurous.  Tensorflow distributions.  CoRR, abs/1711.10604, 2017.
  * [59] S. Varrette, P. Bouvry, H. Cartiaux, and F. Georgatos.  Management of an academic hpc cluster: The ul experience.  07 2014.
  * [60] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  Journal of machine learning research, 12(Oct):2825–2830, 2011.
  * [61] Roger Fletcher.  Practical methods of optimization.  John Wiley & Sons, 2000.
  * [62] M.A. Álvarez, L. Rosasco, and N.D. Lawrence.  Kernels for vector-valued functions: A review.  Foundations and Trends® in Machine Learning, 4(3):195–266, 2012.
  * [63] Stefanos Nikolopoulos, Ioannis Kalogeris, and Vissarion Papadopoulos.  Non-intrusive surrogate modeling for parametrized time-dependent partial differential equations using convolutional autoencoders.  Engineering Applications of Artificial Intelligence, 109:104652, March 2022.
  * [64] Shivam Barwey, Varun Shankar, Venkatasubramanian Viswanathan, and Romit Maulik.  Multiscale graph neural network autoencoders for interpretable scientific machine learning, 2023.