# Automatic selection of inducing points in sparse Gaussian process for approximations of finite element analyses

Heine Havneraas Røstum, Sebastien Gros, Ketil Aas-Jakobsen, Joseph Morlier

► To cite this version:

Heine Havneraas Røstum, Sebastien Gros, Ketil Aas-Jakobsen, Joseph Morlier. Automatic selection of inducing points in sparse Gaussian process for approximations of finite element analyses. Engineering Applications of Artificial Intelligence, 2025, 151, pp.110629. 10.1016/j.engappai.2025.110629. hal-05293629

HAL Id: hal-05293629

https://hal.science/hal-05293629v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Automatic selection of inducing points in sparse Gaussian process for approximations of finite element analyses

Heine Havneraas Røstum<sup>a,b,D*</sup>, Sebastien Gros<sup>a</sup>, Ketil Aas-Jakobsen<sup>b,D</sup>, Joseph Morlier<sup>c,D</sup>

<sup>a</sup> Engineering Cybernetics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

<sup>b</sup> Dr Ing A Aas-Jakobsen AS, Oslo, Norway

<sup>c</sup> ICA, Université de Toulouse, ISAE-SUPAERO, Toulouse, France

## ABSTRACT

Gaussian process regression (GPR) is a widely used regression model, but it has poor scalability. Sparse approximation methods improve scalability by using inducing points to approximate the GPR, but determining the optimal number and placement of these points is challenging. Increasing the number of inducing points generally improves the predictive accuracy, but it comes at a computational cost. This article presents a method to estimate the necessary number of inducing points for accurate predictions of finite element method (FEM) analyses using approximate GPR. The approach leverages the proper orthogonal decomposition (POD) technique, using its modes to determine the inducing points. Results demonstrate that the proposed method identifies a sufficient number of inducing points for approximate GPR to achieve predictive accuracy comparable to full GPR, but with half the training time. This approach ensures computational efficiency without significant loss in accuracy, making it a valuable tool for scalable regression in engineering applications. POD has previously been combined with GPR to provide computationally efficient predictions for the full solution field across unseen variable combinations, treating spatial components separately via reduced basis functions. However, this work treats the spatial component as a variable within the GPR approximation, allowing continuous spatial predictions. This ensures that the covariance in the spatial dimension is captured by a single GPR. The method is applied to simulations of a three-span, post-tensioned concrete girder bridge.

- Approximate Gaussian process regression

- Sparse approximations

- Proper orthogonal decomposition

- Machine learning

- Bridge engineering

## 1. Introduction

The last decades have seen a rapid increase in computational power, enabling more complex high-fidelity simulations to be performed. These simulations often involve an ever increasing size and complexity, such as a higher number of degrees of freedom (DOF) and more intricate material or loading behaviors. For simulations of engineering structures, the finite element method (FEM) is the most popular tool. Multi-query problems, such as in structural optimization and reliability analyses, require the system to be solved for multiple parameter values. The increase in computational power still renders multi-query problems infeasible in many cases due to constraints associated to the computational time and/or memory. To ease computational burden, reduced order modeling (ROM) has been widely used (Lucia et al., 2004; Peherstorfer and Willcox, 2015; Chiplunkar et al., 2018; Guo and Hesthaven, 2018; Drakoulas et al., 2024), which replaces the full order model with a reduced order model with a much smaller dimension, reducing memory and computational cost.

A widely used technique for ROM is Proper Orthogonal Decomposition (POD) (Liang et al., 2002; Amsallem and Farhat, 2012; Paul-Dubois-Taine and Amsallem, 2015; Guo and Hesthaven, 2018; Zhao, 2021; Rocha et al., 2023). In POD, singular value decomposition (SVD) is applied to a collection of “snapshots”, which are results from numerical simulations, to derive reduced basis (RB) functions. Guo and Hesthaven (2018) used Gaussian process regression (GPR) to establish mapping from parameter values to projection coefficients onto the RB space for non-linear FEM analyses, enabling predictions for unseen variable combinations of the full solution field. Guo and Hesthaven (2018) reported good performance in both accuracy and computational efficiency. Demo et al. (2021) and Ortali et al. (2022) applied a similar technique to computational fluid dynamic (CFD) simulations, for the purpose of shape optimization. The Gaussian process is one of many supervised learning techniques available (Zhao et al., 2023; Waqas et al., 2023). One of the benefits of GPR is that the entire dataset is typically used to infer patterns and correlations, resulting in a model that generalizes well (Williams and Rasmussen, 2006). The application

* Corresponding author at: Engineering Cybernetics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway.
E-mail address: hhr@aaj.no (H.H. Røstum).
of GPR to POD opens the possibility of removing human interference from the process, enabling a data-driven process.

The flexibility, interpretability, uncertainty estimation, and the possibility of fully automating the training procedure have made the GPR a popular regression model (Ober et al., 2024). Its drawback is that the exact GPR scales poorly, but its many benefits have led to extensive effort in the machine learning community in recent decades to derive sparse approximation methods to the exact GPR that have better scaling properties (Snelson and Ghahramani, 2005; Quinonero-Candela and Rasmussen, 2005; Titsias, 2009; Liu et al., 2020). Sparse approximations to GPR provide methods to reduce computational complexity by relying on a set of inducing points to approximate the posterior distribution (Snelson and Ghahramani, 2005; Titsias, 2009). The accuracy of sparse approximations depend on the location and number of inducing points that are not known a priori, which poses a challenge for its applicability. As a general rule, increasing the number of inducing points improves the accuracy of the predictions, but increasing the number of inducing points comes at a computational cost. It is therefore a trade-off in having a sufficient number of inducing points, which is high enough to ensure sufficient accuracy of the predictions, but low enough to ensure a computational gain.

Burt et al. (2020) provided bounds for the number of inducing points, but in many situations they may be impractically large. Often, the number of inducing points is determined based on a computational budget (Wen et al., 2022; Iong et al., 2024). In this article, a process for determining the necessary number of inducing points is presented for approximate GPR for FEM analyses. The process is also applicable to other problems that involve approximating observations with a fixed structure, such as simulations for Computational Fluid Dynamics (CFD). As the process for determining the number of points is based on already sampled points, the name active points will be used instead of inducing points. The process relies on the POD modes and provides a computationally efficient method for determining the necessary number of active points. The method may also be used to determine the location of the active points. Although POD with GPR approximates the full solution field, which could be the solution at each node in the FEM analyses, the spatial component is not treated as a variable in the GPR approximation, but is captured by the linear combination of the RB functions. For the approximate GPR applied in this article, the spatial component is treated as a variable when constructing the GPR. Thus, the approximate GPR has additional information, and is also enabling continues predictions in the spatial dimension. In this article, we consider various approximate GPR methods that could be suitable to capture the spatial component as a continuous variable, while retaining computational efficiency. The purpose of this article is twofold:

1. Introduce a method for determining the number of active points when applying approximate GPR for FEM analyses.

2. Comparing the predictive accuracy of regression-based reduced basis and approximate GPR, when applied for approximating the full solution field of FEM simulations.

The methods are applied to approximate the results of simulations of a three-span post-tensioned concrete girder bridge. This case study exemplifies a scenario where the data from each simulation is consistently structured across all variable combinations, a characteristic typical in approximations of FEM and CFD simulations. Such simulations often feature spatial fields with high correlation, making them well suited for POD. However, problems involving time series or spatial modeling, where the data lack a recurring structure and spatial fields are not highly correlated, are unsuitable for POD. Consequently, the proposed method for determining the number of active points is not suitable in these cases.

The article is structured as follows; Section 2 presents the methodology for regression-based RB, GPR, and approximate GPR. Section 3 describes the test-case used; the three-span, post-tensioned concrete bridge girder. In Section 4 various results from applying approximate GPR and POD with GPR are presented and discussed, followed by a conclusion and recommendations for further work in Section 5.

## 2. Methodology

The methodology for approximating FEM simulations using the regression-based reduced basis method (GPR-RB) and approximate Gaussian process regression is initially introduced. In this article, it is separated between full and approximate Gaussian processes by how the dataset that is available to the regressor is treated. The full Gaussian process applies the entire available dataset for the regression process, while the approximate methods approximate the full dataset either by applying only a subset of the available data or by applying a subset of the data more actively than other parts of the dataset.

To address the major challenge for approximate GPR concerning not knowing a priori the necessary number of inducing points, a method for determining the number of inducing (active) points based on the POD modes is proposed in this article. This is followed by two methods for determining the location of the inducing points. The first method proposed in this article is based on the POD mode shapes and is called Peak Singular Values (PSV). The other method was proposed by Burt et al. (2020) which greedily selects points based on the difference between the true and approximate covariance matrix. The process of constructing GPR-RB and approximate GPR, where the number of active points is determined using the proposed PSV method, is illustrated in Fig. 1. In the figure, the arrows represent the workflow and illustrate which subprocesses provide inputs for subsequent steps. Both approximate GPR and GPR-RB begin with sampling, indicated by the gray box. Predictions using GPR-RB follow the subprocesses outlined in the blue box, while predictions using approximate GPR proceed through the red box. As shown in the figure, the proposed PSV method, outlined in a green box, utilizes the same initial subprocesses required for constructing GPR-RB to determine the number of active points for approximate GPR.

In this paper, we use bold uppercase letters to denote matrices, bold lowercase letters to denote vectors, and plain lowercase letters to denote a scalar or variable.

### 2.1. Regression-based reduced basis method

This section provides a brief description of proper orthogonal decomposition (POD) and how Gaussian process regression (GPR) is used to determine the reduced basis (RB) for unseen variable combinations. The approach is abbreviated as GPR-RB. A more in-depth description of GPR-RB applied to FEM simulations is given by Guo and Hesthaven (2018).

#### 2.1.1. Proper orthogonal decomposition

Proper orthogonal decomposition (POD) is a data analysis method with the purpose of obtaining a low-dimensional approximation of a high-dimensional process (Chatterjee, 2000). The results from the FEM simulations, e.g. strains, stresses, displacement, etc., for different design variable combinations $X \in \mathbb{R}^{N_s \times N_d}$ are typically stored in a vector $y \in \mathbb{R}^{(N_h \times N_s)}$, given by;

$$ y = [y_{1,1}, \dots, y_{N_h,1}, \dots, y_{N_h,N_s}]^T, \quad (1) $$

where $N_h$ is the number of degrees of freedom (DOF), $N_d$ is the number of design variables, and $N_s$ is the number of simulations. The vector $y$ can be transformed into a matrix $Y \in \mathbb{R}^{N_h \times N_s}$:

$$ Y = \begin{bmatrix} y_{1,1} & \cdots & y_{1,N_s} \\ \vdots & \ddots & \vdots \\ y_{N_h,1} & \cdots & y_{N_h,N_s} \end{bmatrix}. \quad (2) $$

In the context of POD, $Y$ is often known as the snapshot matrix, and $N_s$ is the number of snapshots. POD is commonly achieved using the singular value decomposition (SVD) (Guo and Hesthaven, 2018; Demo et al., 2021; Ortali et al., 2022). Applying SVD to the matrix $Y$ yields;

$$ Y = U \Sigma V^T. \quad (3) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Process of constructing GPR-RB and approximate GPR.</figcaption></figure>

Here, $U \in \mathbb{R}^{N_h \times N_h}$ and $V \in \mathbb{R}^{N_s \times N_s}$ are orthogonal matrices, and $\Sigma = \text{diag}\{\sigma_1, \sigma_2, \dots, \sigma_{N_s}\}$ are the singular values (eigenvalues) in descending order. The columns of $U$ represent the spatial modes, the diagonal of $\Sigma$ the importance of each mode, and the columns of $V$ how the modes evolve over the simulations/snapshots. A low dimensional truncated POD basis, $L \in \mathbb{R}^{N_h \times N_{POD}}$, which is derived by preserving the first $N_{POD}$ number of modes from $U$ that corresponds to the largest singular values, can be used to calculate the projection coefficients $A \in \mathbb{R}^{N_{POD} \times N_s}$, also sometimes known as the amplitude matrix (Buljak and Maier, 2011; de Gooijer et al., 2021), using the relationship;

$$A = L^{\top} \cdot Y. \qquad (4)$$

Due to the orthogonality of $U$, an approximation of the snapshot matrix $Y$ can be made from the low dimensional POD basis $L$ and the projection coefficients $A$ as;

$$Y \approx L \cdot A. \qquad (5)$$

A common method to select the number of POD modes ($N_{POD}$) to include in $L$ is by evaluating the projection error $\epsilon_{POD}$ by including the

first $N_{POD}$ left singular values of $\Sigma$ (Eckart and Young, 1936; Guo and Hesthaven, 2018; Brunton and Kutz, 2022). Thus $N_{POD}$ can be defined as the minimum integer such that;

$$\frac{\sum_{i=1}^{N_{POD}} \sigma_i^2}{\sum_{i=1}^{N_s} \sigma_i^2} > 1 - \epsilon_{POD}^2. \qquad (6)$$

Other techniques includes selecting the desired rank $N_{POD}$ by identifying a transition in the amount of variance described, which could be “elbows” or “knees” when plotting the cumulative variance, which may signify a transition from singular values that describe important features to those that describe noise (Brunton and Kutz, 2022). Care should be taken when considering the number of modes to include. The excluded modes may still contain critical, albeit less dominant, information, introducing truncation errors. These errors can lead to a loss of fine-scale spatial details, which may be significant in some engineering applications. POD assumes that the spatial field can be linearly decomposed into orthogonal modes. This assumption may be inadequate for highly nonlinear systems, leading to suboptimal representations.
### 2.1.2. Regression-based reduced order solution

For the regression-based reduced order solution, when the required rank $N_{POD}$ of the reduced bases $L$ has been decided, a regression model is used to predict the projection coefficients $\hat{A}(x)$ for unseen variable combinations $x$. Instead of treating the full high-dimensional spatial field as an input, the regression-based reduced order solution operates on a compact set of POD coefficients, reducing the input space dimensionality significantly. Using the training set $D = \{(X, A_i)\}_{i=1}^{N_{POD}}$, one regression model $\hat{A}_i(x)$ is created to approximate each of the $N_{POD}$ projection coefficients in $A \in \mathbb{R}^{N_{POD} \times N_s}$, resulting in the predicted projection coefficients;

$$ \hat{A}(x) = [\hat{A}_1(x), \hat{A}_2(x), \dots, \hat{A}_{N_{POD}}(x)]^T. \quad (7) $$

Approximations of the full set solution $\hat{Y}$ can be retrieved for new variable combinations $x$ using the relationship from Eq. (5);

$$ \hat{Y}(x) = L \cdot \hat{A}(x). \quad (8) $$

In this article, Gaussian process regression (GPR) is applied as the regression model, which will be further described in the following section.

## 2.2. Full Gaussian process regression

### 2.2.1. Gaussian process

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution (Williams and Rasmussen, 2006). A real process $f(x)$ can be formulated completely as a Gaussian process by its mean function $\mu(x)$ and covariance function $k(x, x')$ as:

$$ \mu(x) = \mathbb{E}[f(x)], \quad (9a) $$

$$ k(x, x') = \mathbb{E}[(f(x) - \mu(x))(f(x') - \mu(x'))]. \quad (9b) $$

With the Gaussian process formulated as;

$$ f(x) \sim GP(\mu(x), k(x, x')). \quad (10) $$

The mean is often taken to be zero, which will also be used in this article for notational simplicity, which can be approximated by mean centering the distribution. The covariance function, $k(x, x')$, typically determines the GP's smoothness, and characteristic length-scale. In this article, the stationary Matérn 5/2 kernel is used. The choice of kernel must always be carefully considered, as truncation of higher-order modes could introduce localized distortions, potentially making the effective representation less stationary. The Matérn 5/2 kernel is formulated as (Williams and Rasmussen, 2006);

$$ k(x, x') = \sigma_f^2 \left(1 + \sqrt{5}r + \frac{5}{3}r^2\right) \exp\left(-\sqrt{5}r\right), \quad (11) $$

where $\sigma_f^2$ is the signal variance, and $r$ is a distance metric. The distance metric computes the similarity between points in a continuous space, thus continuous variables $x$ are generally assumed. It is possible to handle discrete variables in GPR, but this is considered to be outside the scope of this article. In this article, a Euclidean distance metric with an anisotropic length-scale $l_i$ is used;

$$ r = \sqrt{\sum_{i=1}^{N_{var}} \left( \frac{x_i - x_i'}{l_i} \right)^2}. \quad (12) $$

When results from FEM simulations are considered, the results are typically noise-free, but it may be beneficial to include a term that allows the regressor function to somewhat deviate from the observations. It is assumed that the results obtained from the FEM simulations, $y$, may be formulated as; $y = f(x) + \epsilon$, where $\epsilon$ is the deviation (sometimes referred to as noise when dealing with noisy observations) and $f(x)$ is the true underlying function at location $x$.

### 2.2.2. Gaussian process regression

When a set of $n$ observations are available, the observations $D = \{(x_i, y_i)\}_{i=1}^n$, can be used to in order to infer a joint Gaussian posterior distribution for unobserved data, $y_*(X_*)$;

$$ y_* | X_*, D \sim \mathcal{N}(\hat{\mu}(X_*), \hat{\sigma}^2(X_*)). \quad (13) $$

where the posterior mean, $\hat{\mu}(X_*)$, and variance, $\hat{\sigma}^2(X_*)$, are calculated as (Williams and Rasmussen, 2006);

$$ \hat{\mu}(X_*) = K(X_*, X)[K(X, X) + \sigma_n^2 I]^{-1} y, \quad (14a) $$

$$ \hat{\sigma}^2(X_*) = K(X_*, X_*) - K(X_*, X)[K(X, X) + \sigma_n^2 I]^{-1} K(X, X_*) \quad (14b) $$

For compactness, $K(X, X)$ will be abbreviated as $K_{ff}$, and $K(X_*, X_*)$ as $K_{**}$, and so forth. The $\sigma_n^2 I$ term is the deviation between the observations and the predictions. $\sigma_n^2$ is the deviation variance, which is often referred to as the noise variance, and $I \in \mathbb{R}^{n \times n}$ is an identity matrix. Although the results from the FEM simulations are noise free, $\sigma_n^2$ is included to allow the predictions some deviation from the observations.

The marginal likelihood is commonly used to determine the model's hyperparameters $\theta$, which includes the noise and signal variance $\sigma_n^2$ and $\sigma_f^2$, and parameters associated to the covariance function. In this article, the hyperparameters are $\theta = \{\sigma_f^2, l_1, \dots, l_{N_{var}}, \sigma_n^2\}$. The logarithmic marginal likelihood is formulated as (Williams and Rasmussen, 2006);

$$ \log p(y | X) = -\frac{1}{2} \log(|K_{ff} + \sigma_n^2 I|) - \frac{1}{2} y^T (K_{ff} + \sigma_n^2 I)^{-1} y - \frac{n}{2} \log(2\pi). \quad (15) $$

The optimal hyperparameters $\theta_{opt}$ are then typically obtained by the maximization problem;

$$ \theta_{opt} = \arg \max_{y} \log p(y | X). \quad (16) $$

The first term in Eq. (15) penalizes high complexity in the covariance function, the second term ensures a good data fit, and the last term is a normalization constant. Optimizing the hyperparameters using the marginal likelihood thus ensures a good fit to the observations and reduces the risk of overfitting by penalizing complex models. A limitation for the application of GPR is the time complexity $O(n^3)$ for the matrix inversion in Eq. (15), where $n$ is the number of datapoints. For large datasets, this computational complexity can result in impractical or unfeasible computational times.

In FEM simulations, the results are obtained at the $N_h$ number of spatial or nodal positions. To approximate the full set of solutions $y$ of Eq. (1) using a single GPR, the matrix containing the different variable combinations $X$ must be expanded to include the spatial position of the results. If the position of all nodes are stored in a position matrix $X^x \in \mathbb{R}^{N_h \times N_x}$ where each row represents the position of each node in the FEM analyses, where $N_x$ is the number of spatial dimensions, which is 3 if three spatial dimensions are considered, the matrix $X$ can be expanded to matrix $X_{exp} \in \mathbb{R}^{(N_h \times N_x) \times (N_x + N_d)}$, which includes the spatial dimension for each variable combination;

$$ X_{exp} = \begin{bmatrix} X_{1,1}^x & \cdots & X_{1,N_x}^x & \cdots & X_{1,1} & \cdots & X_{1,N_d} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ X_{N_h,1}^x & \cdots & X_{N_h,N_x}^x & \cdots & X_{1,1} & \cdots & X_{1,N_d} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ X_{N_h,1}^x & \cdots & X_{N_h,N_x}^x & \cdots & X_{N_s,1} & \cdots & X_{N_s,N_d} \end{bmatrix}. \quad (17) $$

The number of datapoints for approximating the results using a single GPR is thus $n = N_h \cdot N_s$, which may result in a prohibitively large computational time. When applying GPR on the reduced basis (GPR-RB), the correlation between the results is captured by the modal shapes $L$, and GPR is used only to approximate each of the $N_{POD}$ number of modes, where the size of each dataset is $n = N_s$. Fortunately, methods for approximating the full GPR exist, which enables better scaling properties, and will be presented in the following section.
## 2.3. Approximate Gaussian process regression

The time complexity $O(n^3)$ for the matrix inversion in Eq. (15) can be prohibitively large for approximating simulations that involve a high number of degrees of freedom. This section presents methods that may significantly reduce the complexity of the approximations. Starting with perhaps the most intuitive method, Subset of Datapoints (SD), and then presenting two methods where $m$ latent points are projected to the full $n$ points of the dataset. These two methods are known as Projected Process approximation (PP) and its extension Variational Approximation (VA) (Williams and Rasmussen, 2006; Titsias, 2009).

### 2.3.1. Subset of Datapoints

Subset of Datapoints (SD) is simply using a smaller subset of size $m$ of the original $n$ number of datapoints, but still using Eq. (14) for predictions, and Eq. (15) for training. This method is perhaps the most intuitive method to reduce the computational complexity of GPR. Though this may in many cases be considered wasting data, when the results are obtained from FEM simulations there is often little variance between spatially close observations. Thus, some datapoints can often be omitted without the loss of significant information.

### 2.3.2. Projected Process approximation

The Projected Process (PP) approximation is called a projected process because it represents $m < n$ latent function values, but computes a likelihood involving all $n$ datapoints by projecting the $m$ latent points onto $n$ dimensions (Csató and Opper, 2002; Williams and Rasmussen, 2006). This ensures that the information of all points $n$ is absorbed by the $m$ points in the active set. The predictive mean and variance of the PP approximation is given as follows (Williams and Rasmussen, 2006):

$$ \hat{\mu}_* = \mathbf{K}_{m*}^{\top} \Sigma^{-1} \mathbf{K}_{mf} \mathbf{y}, \quad (18a) $$

$$ \hat{\sigma}_*^2 = \mathbf{K}_{**} - \mathbf{k}_{m*}^{\top} \mathbf{K}_{mm}^{-1} \mathbf{K}_{m*} + \sigma_n^2 \mathbf{K}_{m*}^{\top} \Sigma^{-1} \mathbf{K}_{m*}, \quad (18b) $$

where $\Sigma = \sigma_n^2 \mathbf{K}_{mm} + \mathbf{K}_{mf} \mathbf{K}_{fm}^\top$.

The marginal likelihood of the PP is calculated similarly to that of the full GPR in Eq. (15) except that the full covariance matrix $\mathbf{K}_{ff}$, is replaced with that of the Nyström approximation $\mathbf{Q}_{ff} = \mathbf{K}_{mf}^{-1} \mathbf{K}_{mm} \mathbf{K}_{mf}$;

$$ \log p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2} \log(|\mathbf{Q}_{ff} + \sigma_n^2 \mathbf{I}|) - \frac{1}{2} \mathbf{y}^{\top} (\mathbf{Q}_{ff} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{y} - \frac{n}{2} \log(2\pi). \quad (19) $$

### 2.3.3. Variational Approximation

Titsias (2009) proposed a Variational Approximation (VA) method where the predictive distribution is equivalent to that of the PP in Eq. (18), but an extra regularization term is added to the log marginal likelihood. The log marginal likelihood of the variational approximation proposed by Titsias (2009), also known as the evidence lower bound (ELBO), is formulated as;

$$ \log p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2} \log \left( |\mathbf{Q}_{ff} + \sigma_n^2 \mathbf{I}| \right) - \frac{1}{2} \mathbf{y}^{\top} (\mathbf{Q}_{ff} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{y} - \frac{n}{2} \log(2\pi) \\ - \frac{1}{2\sigma_n^2} \text{tr}(\mathbf{K}_{ff} - \mathbf{Q}_{ff}). \quad (20) $$

The last term in Eq. (20) is the KL divergence regularization term, which seeks to minimize the error of predicting the full number of latent values $n$ from the active set, or inducing variables, $m$. The variational method is currently one of the most popular approximations to the full GPR.

## 2.4. Computational complexity

The space and time complexity of the GPR methods relative to the size of the full dataset $n$ and that of the active dataset $m$ is shown in Table 1. The values are based on the values presented by Williams and Rasmussen (2006). The training complexity is due to the time needed to perform the necessary matrix inversions. For large datasets, it is evident that the approximation methods can provide significant computational time and space savings if $m \ll n$.

<table><caption>**Table 1**
Computational space and time complexity of the different GPR methods.</caption><thead><tr><th rowspan="2">Method</th><th rowspan="2">Storage</th><th rowspan="2">Training</th><th colspan="2">Prediction</th></tr><tr><th>Mean</th><th>Variance</th></tr></thead><tbody><tr><td>Full</td><td>$\mathcal{O}(n^2)$</td><td>$\mathcal{O}(n^3)$</td><td>$\mathcal{O}(n)$</td><td>$\mathcal{O}(n^2)$</td></tr><tr><td>SD</td><td>$\mathcal{O}(m^2)$</td><td>$\mathcal{O}(m^3)$</td><td>$\mathcal{O}(m)$</td><td>$\mathcal{O}(m^2)$</td></tr><tr><td>PP</td><td>$\mathcal{O}(mn)$</td><td>$\mathcal{O}(m^2n)$</td><td>$\mathcal{O}(m)$</td><td>$\mathcal{O}(m^2)$</td></tr><tr><td>VA</td><td>$\mathcal{O}(mn)$</td><td>$\mathcal{O}(m^2n)$</td><td>$\mathcal{O}(m)$</td><td>$\mathcal{O}(m^2)$</td></tr></tbody></table>

## 2.5. Selection of active dataset

The variational and projected process approximations involve $m$ latent points, which are projected to the full $n$ points. The quality of the approximations is greatly affected by the number and location of the latent points (Burt et al., 2020). These $m$ points do not have to be part of the already sampled $n$ points, but may be considered as pseudo-data points, commonly called *inducing points*. In this article, the $m$ points are a subset of the $n$ points; therefore, it is considered more accurate to denote them as the *active dataset*. In this section, a novel method is presented for selecting the number and location of the points in the active dataset when using approximate GPR for FEM analyses. We call this method Peak Singular Values (PSV). This method utilizes the POD modes $L$ for selecting the number and location of the active points. The method can be used either to solely determine the number of points in the active dataset, or to determine both the number and the locations of the active points.

Burt et al. (2020) proposed a greedy method to select the location of active points, applicable for variational approximation, demonstrating excellent predictive performance. They called this method greedy variance selection. This method can be combined with the number of active points determined using PSV, and will be briefly presented at the end of this section.

### 2.5.1. Peak Singular Values (PSV)

The approximation methods described in Section 2.3 converges to the full GPR with a sufficient number of points in the active set. The problem is that one does not know a priori what sufficient is. Burt et al. (2020) provided bounds for some methods, but as pointed out by Ober et al. (2024), the methods “...require detailed knowledge of the data generating process, and involve impractically large constants to provide practical recommendations.”. The literature generally does not advise on the number of active points, but rather points out that the number of active points depends on the data and kernel (an appropriate kernel should be used for a dataset), and that accuracy improves with an increase in the number of active points. In Section 2.1.1 the concept of applying SVD to the FEM simulations structured as a snapshot matrix was described. After decomposing the snapshot matrix $\mathbf{Y}$ using SVD, the columns of the $U$ component represents the spatial modes, which decreases in importance from the left to right column. The scheme proposed in this article for determining the number of active points, which can also be used to determine their location, is to identify the most important nodes in each of the $N_{POD}$ number of modes, which are represented by their peaks and troughs. We call this method Peak Singular Values (PSV) and it is performed as follows:

1. **Generate snapshot matrix:** Perform FEM simulations with a limited set of parametric configurations $\mathbf{X}$. The results are stored in a matrix $\mathbf{Y}$, where the columns represent the results for the different parametric configurations, and the rows are the spatial locations, which are the nodal coordinates in FEM analyses.

2. **Decompose the result matrix:** The result matrix is decomposed using SVD. A limit is set to the number of modes $N_{POD}$ to include by setting a lower limit to the amount of variance captured by the most important modes $L$, using Eq. (6).
<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of important nodes found on POD mode.</figcaption></figure>

### 3. Select most important nodes:

For each of the $N_{POD}$ modes in L, the nodes associated with the peaks and troughs are identified. These $N_{PSV}$ number of nodes are the most influential in the representation of the data. Only these nodes from each FEM analysis are added to the active dataset.

### 4. Create active dataset:

The $N_{PSV}$ number of nodes identified have the same location in each analysis. The number of active points can thus be calculated as; $m = N_{PSV} \cdot N_s$. Here $N_s$ is the number of simulations, or FEM analyses, performed.

In Fig. 2 a POD mode is illustrated with a blue line, the nodal positions are illustrated with black dots, and the most important nodes identified using the PSV method is illustrated with red stars.

The proposed PSV method has good scaling properties. The time complexity of the peak selection is $O(N_h N_{POD})$, and the time complexity of performing SVD on the snapshot matrix is $O(N_h N_s^2)$.

### 2.5.2. Greedy maximum variance

Burt et al. (2020) proposed a method for selecting the active dataset using a method they called Greedy Variance selection (GV). This method determines the active points by greedily selecting the points that maximize the difference between the true and approximated covariance matrix; $arg \max (\text{diag} [K_{ff} - Q_{ff}])$. They showed that this method converged more efficiently, i.e. with fewer inducing points, to the log marginal likelihood of the full Gaussian process, Eq. (15), than by using k-means or uniform sub-sampling. They also proposed to alternately optimize the hyperparameters and then reinitialize the inducing points using greedy variance selection with the updated hyperparameters. Ober et al. (2024) extended this to alternately optimize the hyperparameters and then select the inducing points using GV selection until the ELBO, Eq. (20), no longer improved. This method showed excellent predictive capabilities. In this article, variational inference is combined with the PSV method to determine the number of active points $m$. GV selection is then used to select the location of the active points. Approximations are performed both with reinitialization, until the ELBO no longer improves, and without.

## 3. Case study

A case study of a post-tensioned concrete bridge girder is used to evaluate the applicability of the proposed PSV method for determining the size of the active dataset $m$, as well as to investigate the performance of the two approaches, GPR-RB and approximate GPR, in approximating the full solution field of FEM simulations. This section presents the bridge, applied loading, variables, and results that were used as the basis for constructing the regression models.

<table><caption>**Table 2**
Girdler section geometry values.</caption><thead><tr><th></th><th>(H)</th><th>Varies</th></tr></thead><tbody><tr><td>Total height</td><td>(W)</td><td>7000 mm</td></tr><tr><td>Total width</td><td>(T1)</td><td>300 mm</td></tr><tr><td>Minimum cantilever thickness</td><td>(T2)</td><td>350 mm</td></tr><tr><td>Maximum cantilever thickness</td><td>(W1)</td><td>3300 mm</td></tr><tr><td>Base width</td><td>(W2)</td><td>350 mm</td></tr></tbody></table>

<table><caption>**Table 3**
Materials applied in the bridge girder.</caption><thead><tr><th>Material</th><th>Quality</th></tr></thead><tbody><tr><td>Concrete</td><td>B45</td></tr><tr><td>Reinforcement bars</td><td>B500NC</td></tr><tr><td>Post-tension tendons</td><td>Y1860S7</td></tr></tbody></table>

<table><caption>**Table 4**
Applied loads.</caption><tbody><tr><th colspan="3">Permanent loads</th></tr><tr><td>Selfweight concrete</td><td>25.00</td><td>kN/m<sup>3</sup></td></tr><tr><td>Surfacing</td><td>2.50</td><td>kN/m<sup>2</sup></td></tr><tr><td>Edge beams</td><td>8.75</td><td>kN/m</td></tr><tr><td>Railing</td><td>0.50</td><td>kN/m</td></tr><tr><td>Creep</td><td colspan="2">Annex B.1 in Norge (2021)</td></tr><tr><td>Shrinkage</td><td colspan="2">Annex B.2 in Norge (2021)</td></tr><tr><td>Post-tension</td><td colspan="2">Section 5.10 in Norge (2021)</td></tr><tr><th colspan="3">Variable loads</th></tr><tr><td>Uniform traffic - Lane 1</td><td>5.40</td><td>kN/m<sup>2</sup></td></tr><tr><td>Uniform traffic - Lane 2</td><td>2.50</td><td>kN/m<sup>2</sup></td></tr><tr><td>Uniform traffic - Remaining area</td><td>2.50</td><td>kN/m<sup>2</sup></td></tr><tr><td>Tandem traffic - Lane 1</td><td>2 x 300</td><td>kN</td></tr><tr><td>Tandem traffic - Lane 2</td><td>2 x 200</td><td>kN</td></tr><tr><td>Max. air temperature</td><td>36.0</td><td>°C</td></tr><tr><td>Min. air temperature</td><td>-30.0</td><td>°C</td></tr></tbody></table>

### 3.1. Bridge layout

#### 3.1.1. Geometry

The bridge is a continuous three-span, post-tensioned, concrete road bridge, with span lengths of 20 m–30 m–20 m. The girder has a massive cross-section and is illustrated in Fig. 3. In the analyses, only the height of the girder is considered a variable, the remaining values are fixed as indicated in Table 2.

The geometry of the post-tensioned cables along the bridge is illustrated in Fig. 4. The cables start at the girder's center of gravity at each end, drop to their lowest point at each mid-span, and are at their highest at each of the central axes. The side spans, L1 and L3, are of equal length, 20 m. The center span, L2, has a span length of 30 m. The girder is considered to have a pin-support in axis 1, and is freely supported at axis 2 and axis 3.

#### 3.1.2. Materials

The girder consists of concrete, reinforcement bars and a post-tension system. The materials applied to the girder are shown in Table 3.

#### 3.1.3. Loading

The girder is subjected to permanent and variable loads. The permanent loads are comprised of self-weight, creep, shrinkage, super dead-loads and the post-tension force. The variable loads comprise of traffic and temperature loads. The total carriage width of the girder is 7.0 m, which enables two notional lanes. The edge beams, on which the railing is mounted, are assumed to be suspended on the girder cantilevers. Table 4 shows the applied loads. As the cross section and amount of post-tension varies, references to the applied regulations are given when relevant.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Geometry of the girder section.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Geometry of post-tensioned cables along the bridge.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Illustration of the FEM analyses for the bridge girder.</figcaption></figure>

## 3.2. Numerical analyses

The FEM model is a beam element model, see Fig. 5, that is modeled and analyzed using a software called NovaFrame, which is an in-house software at Aas-Jakobsen (Dr. Ing. A.Aas-Jakobsen A/S, 2020). NovaFrame was chosen because of the authors' experience using it and because of the possibility to perform multiple analyses. The analyses are performed using linear-elastic theory to derive sectional forces. The software also has a module which allows the stresses and strains in the concrete, reinforcement and tendons at the section to be calculated by determining the equilibrium between the external and internal forces of a monotone stress field using a discrete approximation to Green's theorem.

### 3.2.1. Level of discretization

Analyses were performed for two levels of discretization. This was done to assess how the methods behave with more data and to indicate the computational scaling when the size of the observed dataset, *n*, is doubled, but where the data not necessarily provide more information, assuming that the model is sufficiently discretized.

1. Each span is modeled using 20 elements (60 elements in total and 61 nodes)

2. Each span is modeled using 40 elements (120 elements in total and 121 nodes)

### 3.2.2. Variables

In the analyses four design variables are considered, as illustrated in Fig. 6. The variables are approximated as continuous, which many engineering design variables may be. For the post-tension system, it is assumed that 8 tendons are applied, with 19 strands in each tendon. The variable is the cross-sectional area of each strand $A_T$, thus the total tendon area is; $8 \times 19 \times A_T$. The variables for the reinforcement bars, $A_{s,U}$ and $A_{s,L}$, are the total reinforcement area applied in the upper and lower parts of the section.

The bounds of the variables, presented in Table 5, ensure that the initial sampling is conducted using reasonable variable values.

<table><caption>Table 5
Upper and lower variable bounds.</caption><thead><tr><th>Variable</th><th></th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>Girder height</td><td>(H)</td><td>1.10 m</td><td>2.00 m</td></tr><tr><td>Strand area</td><td>(A<sub>T</sub>)</td><td>120 mm<sup>2</sup></td><td>200 mm<sup>2</sup></td></tr><tr><td>Upper reinforcement area</td><td>(A<sub>s,U</sub>)</td><td>3000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Lower reinforcement area</td><td>(A<sub>s,L</sub>)</td><td>3000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr></tbody></table>

### 3.2.3. Sampling

A total of 30 analyses were performed, of which 20 were used for training and 10 were used to test/validate the surrogate models. The sampling was performed using a Latin Hypercube Sampling (LHS) scheme, which was first proposed by McKay et al. (1979). The LHS strategy implemented in this article is to randomly generate multiple Latin hypercubes with uniform distribution and then select the sample which has the highest sum of minimum distances, which is considered to have the best space-filling ability of the generated samples. The distribution of the design variables used for training and validation is illustrated in the pairwise plot of Fig. 7.

### 3.2.4. Response: Tendon utilization ratio

The response which is approximated using the surrogate models is the tendons utilization ratio for Ultimate Limit State (ULS). The ULS load factors for the relevant loads applicable to road bridges according to the Norwegian annex of Norge (2003) is presented in Table 6. As described in the beginning of Section 3.2, the tendon strains along the bridge are calculated on the basis of the sectional forces in the software, and the utilization ratio (UR) is the strain calculated by the software ($\epsilon_{Ed}$) divided by the maximum allowable tendon strain ($\epsilon_{Rd}$), which is 10.0 %, according to the Norwegian national annex of the Eurocodes (Norge, 2021).

The tendon utilization ratio along the bridge, where the station number indicates the location on the bridge, for the test and training samples is shown in Figs. 8 and 9, respectively for the model with 61 nodes and 121 nodes, $N_h = 61$ and $N_h = 121$. As the training data
<figure><img src="image_7.png" /><figcaption>Fig. 6. Variables considered in the FEM analyses.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Distribution of initial design samples.</figcaption></figure>

<table><caption>**Table 6**
ULS load factors applied in the analyses.</caption><thead><tr><th rowspan="2">Dominant load</th><th colspan="3">Permanent loads<sup>a</sup></th><th colspan="2">Variable loads</th></tr><tr><th>Selfweight</th><th>Deformation</th><th>Tendons</th><th>Traffic</th><th>Temperature</th></tr></thead><tbody><tr><td>Selfweight</td><td>1.35/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>0.95</td><td>0.84</td></tr><tr><td>Traffic</td><td>1.20/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>1.35</td><td>0.84</td></tr><tr><td>Temperature</td><td>1.20/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>0.95</td><td>1.20</td></tr></tbody></table>

consisted of 20 FEM simulations, $N_s = 20$, the number of datapoints is $n = 61 \cdot 20 = 1220$ and $n = 121 \cdot 20 = 2420$. The figures show that there is little difference in the results when the number of elements is doubled. As expected, the increased level of discretization results in rounder peaks and troughs, but the overall shape of the analyses appears to be the same.

### 3.3. Constructing the GPR models

The methods described in Section 2 are implemented in Python using TensorFlow (Abadi et al., 2015). The estimation of hyperparameters is performed using the gradient-based L-BFGS algorithm from SciPy, where the gradients are calculated using algorithmic differentiation in TensorFlow (Virtanen et al., 2020).

The results of the FEM simulations are approximated using the methods described in Section 2. For the GPR approximation methods,

<table><caption>**Table 7**
Abbreviations used in the results for the Gaussian process regressors.</caption><thead><tr><th>Gaussian process regressor</th><th>Selection of active points</th><th>Abbreviation</th></tr></thead><tbody><tr><td>Full gaussian process</td><td><i>Not applicable</i></td><td>Full</td></tr><tr><td>Subset of datapoints</td><td>Peak singular value</td><td>SD-PSV</td></tr><tr><td>Projected process</td><td>Peak singular value</td><td>PP-PSV</td></tr><tr><td>Variational approximation</td><td>Greedy maximum variance</td><td>VA-GV<sup>a</sup></td></tr></tbody></table>

the location of the active dataset for the subset of datapoints (SD) and projected process (PP) approximation are selected using the peak singular values method (PSV), described in Section 2.5.1. The PP method might be the most intuitive method and it is therefore interesting to see how it compares to the methods where the subset is projected to the full dataset. For the variational approximation (VA) method, the location of the active points is solely selected using greedy maximum variance, as the marginal likelihood of the VA includes a regularization term which penalizes divergence between the true and approximated covariance matrix, that is the $tr(K_{ff} - Q_{ff})$ term. The abbreviations used in the results for the different methods are presented in Table 7.

The regression-based reduced basis method, described in Section 2.1 is abbreviated GPR-RB. There is a significant difference between the size of the datasets that are approximated using the GPR-RB method and the approximate GPR. This is because the spatial variable is included for the training data for the approximate GPR, $X_{exp} \in$

<footnote>a The load factors are given as unfavorable/favorable.</footnote>

<footnote>a Is also performed with re-initialization, as described in Section 2.5.2 and abbreviated VA-GV (reinit).</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Model with 61 nodes — Tendon utilization ratio (UR) along bridge.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Model with 121 nodes — Tendon utilization ratio (UR) along bridge.</figcaption></figure>

$\mathbb{R}^{(N_h \cdot N_s) \times (N_x + N_d)}$. For the GPR-RB method, the spatial dimension correlation is captured by the POD basis $L$ and the regressor captures how the variables affect the projection coefficients $A$. For that purpose, a GPR is constructed for each of the $N_{\text{POD}}$ bases of $A$ using $X \in \mathbb{R}^{N_x \times N_d}$, where $N_d$ for the case-study is the four design variables.

As the marginal likelihood, which is used for selecting the hyperparameters, is multi-modal each regression model is constructed 20 times. Approximate GPRs are constructed for the SD and PP methods using active datasets determined by PSV using 1 to 20 modes, where $N_s = 20$. The GPR-RB method is constructed similarly using $N_{\text{POD}}$ ranging from 1 to 20. The variational approximation is constructed using the sizes of the active datasets $m$ ranging from 0.01$n$ to 0.80$n$ with a step size of $(0.80n - 0.01n)/50$.

# 4. Results

In the following section the results from approximating the FEM simulations are presented. In general, the median of the 20 approximations for each method is presented, unless stated otherwise. The results from applying SVD on the snapshot matrix is first presented in order to see how these results may explain the data. This is followed by the predictive capability of the different methods and how this is affected by the number of active points, before the computational time is presented.

## 4.1. Eigenvalue analyses

Fig. 10 illustrates the first four POD modes, i.e. the first four columns from the matrix $U$. Here, “Mode 1” is the leftmost column of $U$,

and then continues from left to right with increasing mode number. As described in Section 2.5.1, modal shapes are used to select the active data for the SD and PP methods by identifying the modal peaks and troughs. Comparing the first modal shape with that of the overall shape of the training data results, which are illustrated in Figs. 8 and 9, shows that the first mode captures the general behavior well.

The cumulative variance explained by the POD modes is illustrated in Fig. 10. As the first mode captures the shape of the training data well, it also captures a significant part of the cumulative variance. The figure shows that about 85% of the variance is captured by the first mode and that the first four modes capture about 99% of the variance for the FEM simulations that consisted of 61 nodes, while five modes where necessary to capture 99% of the cumulative variance when the number of nodes was doubled.

By using the PSV method described in Section 2.5.1, the total number of points added to the active dataset by including a number of POD modes is illustrated in Fig. 12. The total number of points for the $N_s = 20$ simulations is 1220 for the model with 61 nodes and 2420 for the model with 121 nodes. The figure shows that about 140 points (7·20) were identified to ensure a good approximation to the first mode, and that approximately 560 points (28·20) were identified to ensure a good fit to the first 4 modes for the training data that were simulated using a model with 61 nodes, which represented approximately 99% of the cumulative variance. For the training data that were based on the model with 121 nodes, 580 points (29·20) were identified for the first four modes, which corresponded to slightly less than 99% of the cumulative variance. Fig. 12 also illustrates that up to 5 modes, the number of points identified using PSV for the results of the model with
<figure><img src="image_11.png" /><figcaption>Fig. 10. First four POD modes.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Cumulative variance explained by number of POD modes.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Cumulative number of SVD peak points.</figcaption></figure>

61 and 121 nodes is about the same, but that it starts to deviate at that point. This corresponds well with the results presented in Fig. 11, where most of the variance is explained using the first 4 to 5 modes, and that many of the remaining modes may represent noise.

## 4.2. Predictive capability and effect of the active set's size

The metrics Root Mean Square Error (RMSE) and Negative Log Predictive Density (NLPD) are used to assess the predictive accuracy of the different methods for an increasing number of points in the active dataset. The predictions are made to simulate the test data, which consists of the simulation results for combinations of design variables that were not part of the training data. Figs. 13 and 14 show, respectively,

the RMSE and NLPD with an increasing size of the active set. RMSE evaluates how well the model's mean predictions align with observed values, while NLPD evaluates the entire predictive distribution by incorporating both the mean prediction and the associated uncertainty. The RMSE and NLPD of a Gaussian distribution are calculated as;

$$ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{\mu}_i)^2}, \qquad (21) $$

$$ \text{NLPD} = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{(y_i - \hat{\mu}_i)^2}{2\hat{\sigma}_i^2} + \frac{1}{2} \log(2\pi\hat{\sigma}_i^2) \right). \qquad (22) $$

Here n is the number of observations, y_i the actual observed value, and ̂μ_i and ̂σ_i^2 are, respectively, the mean and variance of the Gaussian process predictor.

As GPR-RB does not have an active dataset, the number of points determined by the PSV method using the associated N<sub>POD</sub> modes is used in the plots. The vertical dashed-dotted lines signify the number of points determined by the PSV method using the first and the first four leftmost columns of U, recalling from Fig. 11 that about 99% of the variance was explained by the first four POD modes. We note that there is little difference in the approximation results for simulations using N<sub>h</sub> = 61 and N<sub>h</sub> = 121, which means that the methods are relatively independent of the discretization level, as long as the discretization level is sufficient to capture the overall behavior. The exception is the VA methods, which yield somewhat poorer NLPD results for the N<sub>h</sub> = 121 simulations until the active set includes about 580 points.

The figures also show that the general trend is an increased predictive capability with an increase in the size of the active dataset, as expected, but we also note a trough for the VA methods when the size of the active set is about the total number of points identified using PSV for the first POD mode, represented by the leftmost vertical dashed-dotted line. The VA method with reinitialization and the GPR-RB method appear to converge to the predictive capabilities of the full GP when the size of the active dataset is equal to that found by the PSV method using the first four POD modes, represented by the rightmost vertical dashed-dotted line. Here, the PP method also shows a trough, but then start to diverge from the predictive capabilities of the full GP. The GPR-RB generally shows good predictive capabilities, even with a limited number of modes included in N<sub>POD</sub>. The SD method generally obtains the poorest predictive results of RMSE, which indicates that the predictive mean ̂μ of the SD method is the least capable of predicting unseen data.

The generally inferior predictive accuracy of the SD method is because the posterior mean and variance are calculated on the basis of significantly fewer points, involving only the subset of active points
<figure><img src="image_14.png" /><figcaption>Fig. 13. Median RMSE of the surrogate models to the test data with increasing number of points in the active set.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Median NLPD of the surrogate models to the test data with increasing number of points in the active set.</figcaption></figure>

m. The other methods benefit by using the information of the full set of available points *n* for the computations of the predictive mean and variance. For the NLPD results, where the probabilistic property of the methods is included, i.e. the predictive variance $\hat{\sigma}^2$, the SD method is shown to perform better, and the VA method without re-initializing performs the poorest. It is evident that the reinitialization of the VA-GV method greatly improves its predictive capabilities. The reinitialization ensures that the hyperparameters are better aligned with the location of the active points, which explains the enhanced predictive capabilities.

The predictive metrics give a good indication on how the different methods are able to reproduce unseen data on average. For uncertainty quantification, one may require the surrogate to be able to predict the response of the entire structure well, while for optimization purposes one may be more interested in the peaks and troughs of the response, which represent areas where the structure may have to be strengthened to ensure an adequate design, or it could be weakened, while still ensuring a safe design, but reducing costs.

Fig. 15 illustrates how the different approximation methods are able to predict the tendon utilization for one of the test datasets. The results are for the FEM simulation with 61 nodes. The size of the active dataset is 560 points, which corresponds to the first four POD modes, and similarly $N_{POD} = 4$ for the GPR-RB method. The full Gaussian predictor is shown with a solid black line and the dashed black line represents its $\pm 1.65\hat{\sigma}$, and is included for reference. For the approximate GPR methods, the shaded area represents their $\pm 1.65\hat{\sigma}$.

The results of the GPR-RB method, which is discrete, are shown as squares representing the mean value, and the height of the vertical line represents its $\pm 1.65\hat{\sigma}$. The top row, illustrating the SD and PP approximation methods, shows that the uncertainty grows rapidly in certain areas, which is due to the distance from the points in the active set. They are able to capture the peaks well, but the uncertainty grows rapidly away from the peaks and is less able to predict the response near the troughs. The second row, representing the VA methods, also shows good approximations, especially the VA with reinitialization. The VA methods are generally more confident, that is, have a lower predictive variance $\hat{\sigma}^2$, and especially the VA method without reinitialization is perhaps overconfident, which explains the poor results for the NLPD metric in Fig. 14. The GPR-RB method shows predictive capabilities similar to that of the VA method with reinitialization, but with less confidence in the areas close to the peaks.

Figs. 13 and 14 showed a trough in the VA methods' prediction at about the level where the number of active points was equal to the number of points identified using PSV with the first POD mode. The prediction of the VA method with reinitialization for the same test case as illustrated in Fig. 15, but using 140 and 560 points in the active set, which is the number of points identified using PSV for the first and four most important POD modes, is illustrated in Fig. 16. The figure shows that the approximation using only 140 points in the active set is able to provide a prediction accuracy that does not deviate much from the predictions of the same method with an active set that is four times greater.
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Predicting response of a simulation for the test-data using the different methods (shaded area is $\hat{\mu} \pm 1.65\hat{\sigma}$ of the approximate GPRs and black dashed lines are the $\hat{\mu} \pm 1.65\hat{\sigma}$ of the full GPR).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Predicting response using VA-GV (reinit) with an active set ($m$) of 140 and 560 points (shaded area is $\hat{\mu} \pm 1.65\hat{\sigma}$ of prediction with $m = 560$ and dashed lines are the $\hat{\mu} \pm 1.65\hat{\sigma}$ of prediction with $m = 140$).</figcaption></figure>

### 4.3. Computational time

The calculations were performed on a personal computer with a dual processor of 3.60 GHz × 3.60 GHz and 48 GB RAM. Although the

training time could be further optimized, they are included in order to illustrate how the different methods perform relative to each other during training. The clock wall time for training the different methods is illustrated in Fig. 17. The solid lines in the figure illustrate the median
<figure><img src="image_18.png" /><figcaption>Fig. 17. Median clock wall training time with increasing number of points in active set (shaded area is the 50% confidence interval).</figcaption></figure>

values and the shaded areas are the 50% confidence interval. The training time for the simulation approximations with $N_h = 61$ shows that there is little, if any, time saving by using approximate methods compared to the full GPR. For the GPR-RB method, the number of points determined by the PSV method using the associated $N_{POD}$ modes is applied as the values along the horizontal axes in the figure. The figure shows that the SD and GPR-RB methods are the most computationally effective methods. The construction of regression models for the GPR-RB method involves a significantly smaller dataset, since the spatial variable is not included. The lower computational cost of the SD method is due to the fact that it only uses the covariance matrix of the active dataset $K_{mm}$, while the remaining approximation methods involve the more complex covariance matrix $Q_{ff}$. Locating the position of the active dataset is a little more time consuming for the VA methods, and the reinitialization process also comes at a computational cost.

The NLPD and the RMSE of both the VA with reinitialization and the GPR-RB converged to the full GPR when the active dataset for the VA was about 580 points and the GPR-RB included 4 modes. For simulations where $N_h = 121$, both methods required considerably less computational time than the full GPR. Thus, both methods can provide a significant reduction in computational time, provided that the dataset is of a certain magnitude, while retaining accuracy.

## 5. Conclusion

This article has demonstrated that the peaks and troughs of modes derived from proper orthogonal decomposition (POD) can effectively be utilized to determine the required number of active points needed to ensure accurate predictions when applying approximate GPR to FEM analysis results. The method for selecting the number and/or location of the active dataset is referred to as Peak Singular Values (PSV) in this article. Efficiency and predictive accuracy were evaluated using results from FEM simulations of a three-span post-tensioned concrete bridge girder, where the tendon strain utilization ratio was approximated. The study also demonstrated that the method ensures good scaling properties, making it suitable for efficiently approximating large-scale FEM simulations. This article also demonstrated that both approximate Gaussian process regression (GPR) and the regression-based reduced basis method (GPR-RB) are effective in efficiently approximating the full solution set of FEM analyses, in terms of computational time, for unobserved variable combinations.

The location and number of points in the active set are critical to ensuring predictive accuracy. The PSV method, used to determine the number of active datapoints, combined with the variational approximation (VA) method with reinitialization, demonstrated excellent

predictive performance. In the case study, the predictive metrics of the VA method with reinitialization converged to those of the full GPR when the number of points in the active dataset matched those identified by the PSV method, utilizing the most significant modes, specifically, the modes accounting for 99% of the POD variance. The PSV method was also employed to select the location of active points for the approximate GPR methods, namely the subset of datapoints (SD) and the projected process (PP). While both SD and PP achieved reasonable accuracy, their predictive performance was somewhat inferior to that of the VA method with reinitialization when using the same number of active points.

The results showed that the GPR-RB and the variational approximation (VA) of the GPR with reinitialization were most data-efficiently, i.e. with the least amount of active usage of the high-dimensional data, able to converge to the predictive accuracy of the full GPR. Whereas the spatial variable for the GPR-RB is captured by the linear combination of the RB functions, enabling the full discrete solution field of the FEM simulations to be recovered, the spatial variable is included in the construction of the approximate GPR which allows predictions to be made in the continues spatial dimension. Including the spatial variable provides the GPR with additional information, and the VA methods were generally more confident in their predictions. That is, they showed a lower predictive variance than the GPR-RB, especially in the vicinity of the peak response. The additional variable for the approximate GPR methods results in a higher computational cost.

As the VA method was generally more confident in its predictions than the GPR-RB method, it would be interesting to study which approach is more efficient when applied in optimization. Adopting a deep Gaussian process to the SD where the subset is selected using the PSV method may increase the predictive accuracy, while retaining computational efficiency. The article focuses solely on continuous variables due to the inherent properties of the GPR kernel. However, it would be valuable to investigate how discrete variables, such as material properties, could be incorporated into the process.

## CRediT authorship contribution statement

**Heine Havneraas Røstum:** Writing – review & editing, Writing – original draft, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Sebastien Gros:** Supervision, Project administration. **Ketil Aas-Jakobsen:** Supervision, Project administration. **Joseph Morlier:** Visualization, Supervision.
## Replication of results

Results from the finite element analyses derived from the case study of Section 3 can be made available upon request by e-mailing the corresponding author. The methods used for the creation of surrogate models were performed as described in the article.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work is funded by the Research Council of Norway [project number 332126] and by Dr Ing A Aas-Jakobsen AS.

## Data availability

Data will be made available on request.

## References

Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G.S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., Zheng, X., 2015. TensorFlow: Large-scale machine learning on heterogeneous systems. URL https://www.tensorflow.org/, Software available from tensorflow.org.

Amsallem, D., Farhat, C., 2012. Stabilization of projection-based reduced-order models. Internat. J. Numer. Methods Engrg. 91 (4), 358-377.

Brunton, S.L., Kutz, J.N., 2022. Data-driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, Cambridge, England.

Buljak, V., Maier, G., 2011. Proper orthogonal decomposition and radial basis functions in material characterization based on instrumented indentation. Eng. Struct. 33 (2), 492-501.

Burt, D.R., Rasmussen, C.E., Van Der Wilk, M., 2020. Convergence of sparse variational inference in Gaussian processes regression. J. Mach. Learn. Res. 21 (131), 1-63.

Chatterjee, A., 2000. An introduction to the proper orthogonal decomposition. Current Sci. 808-817.

Chiplunkar, A., Bosco, E., Morlier, J., 2018. Gaussian process for aerodynamic pressures prediction in fast fluid structure interaction simulations. In: Advances in Structural and Multidisciplinary Optimization: Proceedings of the 12th World Congress of Structural and Multidisciplinary Optimization (WCSMO12) 12. Springer, pp. 221-233.

Csató, L., Opper, M., 2002. Sparse on-line Gaussian processes. Neural Comput. 14 (3), 641-668.

de Gooijer, B.M., Havinga, J., Geijselaers, H.J., van den Boogaard, A.H., 2021. Evaluation of POD based surrogate models of fields resulting from nonlinear FEM simulations. Adv. Model. Simul. Eng. Sci. 8 (1), 25.

Demo, N., Tezzele, M., Mola, A., Rozza, G., 2021. Hull shape design optimization with parameter space and model reductions, and self-learning mesh morphing. J. Mar. Sci. Eng. 9 (2), 185.

Dr. Ing. A.Aas-Jakobsen A/S, 2020. NovaFrame v.6.2.

Drakoulas, G., Gortsas, T.V., Polyzos, D., 2024. Physics-based reduced order modeling for uncertainty quantification of guided wave propagation using bayesian optimization. Eng. Appl. Artif. Intell. 133, 108531.

Eckart, C., Young, G., 1936. The approximation of one matrix by another of lower rank. Psychometrika 1 (3), 211-218.

Guo, M., Hesthaven, J.S., 2018. Reduced order modeling for nonlinear structural analysis using Gaussian process regression. Comput. Methods Appl. Mech. Engrg. 341, 807-826.

Iong, D., McAnear, M., Qu, Y., Zou, S., Chen, G.T.Y., 2024. Sparse variational contaminated noise Gaussian process regression for forecasting geomagnetic perturbations. arXiv preprint arXiv:2402.17570.

Liang, Y., Lee, H., Lim, S., Lin, W., Lee, K., Wu, C., 2002. Proper orthogonal decomposition and its applications—Part I: Theory. J. Sound Vib. 252 (3), 527-544.

Liu, H., Ong, Y.-S., Shen, X., Cai, J., 2020. When Gaussian process meets big data: A review of scalable GPs. IEEE Trans. Neural Netw. Learn. Syst. 31 (11), 4405-4423.

Lucia, D.J., Beran, P.S., Silva, W.A., 2004. Reduced-order modeling: new approaches for computational physics. Prog. Aerosp. Sci. 40 (1-2), 51-117.

McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239-245.

Norge, S., 2003. Eurocode: Basis of Structural Design. Standard Norge, Oslo, Norway, EN 1990:2002+A1+AC+NA.

Norge, S., 2021. Eurocode 2: Design of Concrete Structures, Part 1-1: General Rules and Rules for Buildings. Standard Norge, Oslo, Norway, EN 1992-1-1:2004+A1+NA.

Ober, S.W., Artemev, A., Wagenländer, M., Grobins, R., van der Wilk, M., 2024. Recommendations for baselines and benchmarking approximate Gaussian processes. arXiv preprint arXiv:2402.09849.

Ortali, G., Demo, N., Rozza, G., et al., 2022. A Gaussian process regression approach within a data-driven POD framework for engineering problems in fluid dynamics. Math. Eng. 4 (3), 1-16.

Paul-Dubois-Taine, A., Amsallem, D., 2015. An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models. Internat. J. Numer. Methods Engrg. 102 (5), 1262-1292.

Peherstorfer, B., Willcox, K., 2015. Dynamic data-driven reduced-order models. Comput. Methods Appl. Mech. Engrg. 291, 21-41.

Quinonero-Candela, J., Rasmussen, C.E., 2005. A unifying view of sparse approximate Gaussian process regression. J. Mach. Learn. Res. 6, 1939-1959.

Rocha, P.R.B., de Sousa Almeida, J.L., de Paula Gomes, M.S., Junior, A.C.N., 2023. Reduced-order modeling of the two-dimensional Rayleigh-Bénard convection flow through a non-intrusive operator inference. Eng. Appl. Artif. Intell. 126, 106923.

Snelson, E., Ghahramani, Z., 2005. Sparse Gaussian processes using pseudo-inputs. Adv. Neural Inf. Process. Syst. 18.

Titsias, M., 2009. Variational learning of inducing variables in sparse Gaussian processes. In: Artificial Intelligence and Statistics. PMLR, pp. 567-574.

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al., 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 17 (3), 261-272.

Waqas, M., Tahir, M.A., Qureshi, R., 2023. Deep Gaussian mixture model based instance relevance estimation for multiple instance learning applications. Appl. Intell. 53 (9), 10310-10325.

Wen, H., Ma, J., Gu, J., Yuan, L., Jin, Z., 2022. Sparse variational Gaussian process based day-ahead probabilistic wind power forecasting. IEEE Trans. Sustain. Energy 13 (2), 957-970.

Williams, C.K., Rasmussen, C.E., 2006. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA.

Zhao, H., 2021. A reduced order model based on machine learning for numerical analysis: An application to geomechanics. Eng. Appl. Artif. Intell. 100, 104194.

Zhao, H., Wu, Y., Deng, W., 2023. An interpretable dynamic inference system based on fuzzy broad learning. IEEE Trans. Instrum. Meas. ..