# Scalable Sparse Co-Kriging for Multi-Fidelity Data Fusion: An Application to Aerodynamics

Mauricio Castaño-Aguirre, Andrés F. López-Lopera, Nathalie Bartoli, Franck

Massa, Thierry Lefebvre

► To cite this version:

Mauricio Castaño-Aguirre, Andrés F. López-Lopera, Nathalie Bartoli, Franck Massa, Thierry Lefebvre. Scalable Sparse Co-Kriging for Multi-Fidelity Data Fusion: An Application to Aerodynamics. 2026. [hal-05530684](https://hal.science/hal-05530684v1)

HAL Id: hal-05530684

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

Preprint submitted on 27 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Scalable Sparse Co-Kriging for Multi-Fidelity Data Fusion: An Application to Aerodynamics

Mauricio Castaño-Aguirre$^{a, e, f, *}$, Andrés F. López-Lopera$^{b, c}$, Nathalie Bartoli$^{a, d}$, Franck Massa$^{e, f}$ and Thierry Lefebvre$^{a, d}$

**Keywords.** Surrogate modeling Multi-fidelity learning Sparse approximation Large-scale aerodynamic data fusion

**Abstract.** Multi-fidelity Kriging surrogate modeling combines data of varying accuracy, such as experimental tests and numerical simulations, to improve predictive performance. Autoregressive models are commonly used to capture correlations across fidelity levels, offering interpretability compared to purely data-driven approaches. However, existing frameworks often suffer from high computational costs or rely on restrictive assumptions, such as requiring nested designs of experiments (DoEs). In an engineering context, the latter assumption implies that experimental tests must be conducted at the same input conditions as those used in numerical simulations. This requirement cannot always be met due to safety constraints in tests, which limit flexibility in data acquisition. In this work, we revisit the multi-fidelity Kriging approach and introduce a generalized co-Kriging framework that jointly models multiple fidelity levels through a correlation structure based on the autoregressive model. To ensure scalability, we extend sparse approximation techniques to the multi-fidelity setting, maintaining the efficiency of classical sparse Kriging while improving accuracy. We evaluate our framework on both nested and non-nested DoEs. On a large-scale multi-fidelity aerodynamic dataset that combines wind tunnel experiments with CFD simulations, the percentage of standardized residuals falling outside the confidence intervals is reduced by approximately 75% compared to standard sparse Kriging, highlighting the benefit of incorporating multi-fidelity information in complex engineering applications.

## 1. Introduction

In modern aerospace design, aerodynamic data for aerospace vehicles are derived from various sources, including different types of computational fluid dynamics (CFD) simulations, wind tunnel (WT) experiments, and in-flight tests (Johnson, Tinoco and Yu, 2005; Stoliker, 2005; Malik and Bushnell, 2012; Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Mavripilis, 2014). Each of these sources provides information about the same aircraft but with distinct characteristics in terms of accuracy, cost, cycle time, availability, and flight envelope limits (Malik and Bushnell, 2012). Flight test data, for instance, are the most representative of real-world conditions but are constrained by high costs, safety considerations, and limited accessibility during most of the design phase. Conversely, numerical simulations, particularly CFD, are highly beneficial in the early design stages, given the substantial progress made in this field over recent decades (Johnson et al., 2005). However, their accuracy is still limited by factors such as turbulence modeling and discretization errors. Furthermore, while CFD allows for simulating complex flow conditions, the ability to explore thousands of flight scenarios is constrained by the significant computational resources required. Finally, WT testing lies somewhat in between.

Although the aforementioned data sources were sometimes seen as competitors in the past, the growing consensus now recognizes their complementary nature (Malik and Bushnell, 2012). Their differences in characteristics and the aerodynamic data they produce highlight this complementarity. For example, CFD simulations can provide detailed mappings of the local pressure distribution across an aircraft's surface but are typically limited to a small number of flow conditions due to the high computational cost of each simulation. WT experiments, on the other hand, offer limited pressure measurements, constrained by physical space in the model, but allow for rapid exploration of new flow conditions by adjusting the tunnel settings. Given these differing but complementary strengths, there is a pressing need for robust data fusion methods that unify these datasets into a single framework (Meliani, Bartoli, Lefebvre, Bouhlel, Martins and Morlier, 2019; Charayron, 2023). Such approaches aim to fully exploit the knowledge provided

<footnote>$^a$DTIS, ONERA, Université de Toulouse, Toulouse, 31000, France</footnote>

<footnote>$^b$IMAG, Université de Montpellier, CNRS, Montpellier, 34090, France</footnote>

<footnote>$^c$Inria, LEMON, Montpellier, 34095, France</footnote>

<footnote>$^d$Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, Toulouse, 31000, France</footnote>

<footnote>$^e$LAMIH UMR CNRS 8201, Université Polytechnique Hauts-de-France, Valenciennes, 59313, Hauts-de-France, France</footnote>

<footnote>$^f$INSA Hauts-de-France, Valenciennes, 59313, Hauts-de-France, France</footnote>

<footnote>*Corresponding author</footnote>
by the various data sources, maximizing their collective utility in aerospace design and analysis. This characteristic
of leveraging complementary data sources is frequently encountered in real-world scenarios (Miele, Karve and
Mahadevan, 2023).

As demonstrated in (Kennedy and O'Hagan, 2000; O'Hagan, 2006; Le Gratiet, 2013; Meliani et al., 2019), data fusion can be achieved using multi-fidelity Kriging (MFK) based on an autoregressive model, which serves to link the different data sources. Kriging, also known as Gaussian process (GP) interpolation in broader scientific literature beyond geostatistics<fnref n="1" />, is widely used and is a powerful tool for approximating functions based on a set of observations (Laslett, 1994; Stein, 1999; Wang, Qiu, Gao and Jiang, 2023; Marrel and Iooss, 2024). These models belong to a specific class of surrogate models, where the function response is assumed to be a realization of a GP. In the multi-fidelity context, loosely speaking, each data source is modeled as a GP, with the relationship between fidelity levels often established through an autoregressive model (Kennedy and O'Hagan, 2000; O'Hagan, 2006; Le Gratiet, 2013; Shang, Su, Fang, Zeng and Zhang, 2023; Lu, Li, Mi and Huang, 2024).

A particularly efficient recursive formulation for MFK has been extensively studied in (Le Gratiet and Garnier, 2014), seeking to promote scalability of the work in (Kennedy and O'Hagan, 2000). This approach has been widely adopted for aerodynamic data fusion (Meliani et al., 2019; Brevault, Balesdent and Hebbal, 2020; Arenzana, López-Lopera, Mouton, Bartoli and Lefebvre, 2021; Charayron, 2023) and other engineering fields (Keane, 2012; Fernández-Godino, 2023; Toal, 2023). Despite the benefits led by the recursive MFK formulation, the model still presents some limitations. First, regarding input conditions, the model requires nested-input structures, i.e., input data at high fidelity (HF) levels must be a subset of low fidelity (LF) input data. In our aerodynamic context, this implies that WT experiments must be conducted at the same input conditions as those used in CFD simulations. This requirement cannot always be fulfilled due to safety constraints in tests, which limit flexibility in data acquisition. Therefore, such an assumption, while simplifying comparisons between fidelity levels, imposes major constraints on design space exploration. Moreover, the nested structure constrains the acquisition of additional data in informative regions, as HF observations must conform to a predefined subset (Meliani et al., 2019), thereby reducing the model's flexibility and applicability.

Han and Görtz (2012) have adapted the work in (Kennedy and O'Hagan, 2000) to handle non-nested Design of Experiments (DoEs). The proposed framework, referred to as Hierarchical Kriging (HK), maps the trend of the LF function to the sampled HF data, instead of considering the full LF distribution. Thus, as only information from point estimations is transferred across fidelity levels, one limitation of HK relies on the lack of uncertainty propagation.

Both the recursive MFK and HK frameworks inherit the computational burden of standard Kriging-based models, which scales cubically with the size of the input data in each fidelity level. This well-known limitation renders Kriging computationally expensive when considering large datasets (Rasmussen and Williams, 2005). As a result, exact inference becomes computationally intractable for datasets exceeding a few thousand observations.

One approach to address this scalability issue is through the use of sparse approximations, which aim to retain the predictive power of exact Kriging while significantly reducing computational complexity. The key idea behind these sparse models is to approximate the full covariance structure using a smaller set of inducing points. Rather than directly modeling dependencies across all training points, sparse Kriging (SK), more generally known as sparse GPs in the literature, introduces a structured approximation that enables efficient computation without compromising accuracy significantly. Several methods have been proposed to promote sparsity. One of the most widely used approaches is the Fully Independent Training Conditional (FITC) method (Snelson and Ghahramani, 2005), which approximates the covariance structure by assuming independence between training points conditioned on a set of inducing points. This leads to a low-rank approximation of the covariance matrix that significantly reduces computational complexity while maintaining local accuracy. Another widely used sparse approach is the Variational Free Energy (VFE) method, developed by (Titsias, 2009a), which formulates inference as a variational optimization problem. The VFE approach selects inducing points by maximizing a lower bound on the marginal likelihood, leading to a more theoretically grounded approximation that minimizes information loss.

A few works have been proposed to handle the computational complexity of MFK-based approaches using sparse approximations. For a specific two-fidelity Kriging formulation, a Nyström approximation of the sample covariance matrices is used in (Zaytsev and Burnaev, 2017) to reduce computational cost. This approach, known as the Subset of Regressors (SoR) method (Quinonero-Candela and Rasmussen, 2005), offers computational advantages. However, it relies on a rank-deficient prior that causes predictive uncertainty to collapse away from the inducing points, resulting

<footnote n="1">¹In this work, we also use the term Kriging to refer to GP regression. While this term was introduced to refer to GP interpolation tasks (Krige, 1951; Matheron, 1963), it has also been adopted for regression settings that account for noisy data (Christianson, Pollyea and Gramacy, 2023).</footnote>
in overconfident and potentially biased predictions. In (Perdikaris, Venturi, Royset and Karniadakis, 2015), another computationally efficient framework is proposed based on sparse precision matrices of Gaussian–Markov random fields. Their work achieves computational efficiency by imposing strong locality assumptions, which often reduce their ability to model complex and global dependencies. The resulting methodology is limited to model random fields, i.e., $D = 2$. Other potential approaches rely on projections onto basis functions to control computational complexity through the number of basis components. For example, (Rumpfkeil, Bryson and Beran, 2022) propose a nested MFK framework in which responses are represented using sparse Polynomial Chaos Expansions (PCEs). This formulation enables accurate predictions at a relatively low overall computational cost by incorporating a hybrid bridge-function strategy that captures both local variations and global trends. Nevertheless, the performance of the model depends on the choice of polynomial terms for a given task, assumptions on model smoothness, and restrictive nested-input structures, which may require advanced domain expertise for fine tuning. Finally, if the multi-fidelity problem is interpreted as a multi-output one, then it is possible to consider sparse Multi-output GPs (SMoGPs) (Álvarez and Lawrence, 2008), where each fidelity level will be treated as a distinct output. In that case, SMoGPs can be seen as a purely data-driven approach as correlations between outputs will be learned directly from the data without accounting for the autoregressive driven model.

In this work, we revisit the MFK approach originally introduced in (Kennedy and O'Hagan, 2000) and extend it to address three key challenges: handling non-nested DoEs, allowing uncertainty propagation across fidelity levels, and enabling scalability to large-scale datasets. The resulting multi-fidelity co-Kriging (MFCK) framework correlates all fidelity levels through an autoregressive structure while accommodating non-nested data. To ensure scalability, we further adapt the FITC and VFE sparse approximations to a general multi-fidelity setting by addressing both limitations simultaneously. Therefore, MFCK provides a generalized framework integrated with more theoretically grounded approaches, like the VFE framework of (Titsias, 2009a) directly into a multi-fidelity setting, without altering the underlying correlation model. As a result, our proposed approach extends sparse variational inference to multi-fidelity in a principled and scalable manner.

The organization of this paper is as follows. Section 2 further reviews related work and discusses the limitations of existing MFK approaches in (Kennedy and O'Hagan, 2000; Le Gratiet, 2013; Han and Görtz, 2012). Section 3 presents our MFCK formulation while deriving the cross-covariance expressions between fidelity levels. The extension of sparse approximation techniques (i.e., FITC and VFE) to the proposed MFCK framework is detailed in Section 4 with numerical examples illustrating the models. Section 5 provides experimental results on both a synthetic benchmark and a real multi-fidelity aerodynamic dataset. Finally, Section 6 concludes the paper and outlines directions for future research.

## 2. Multi-fidelity Kriging based on the autoregressive model

As suggested by (Kennedy and O'Hagan, 2000), multi-fidelity modeling can be addressed by adopting an autoregressive structure across fidelity levels. For $\kappa \in \{0, ..., L\}$ with $L \ge 1$, let $y_\kappa : D \to \mathbb{R}$, where $D \subseteq \mathbb{R}^D$, denote the function corresponding to the $\kappa$-th fidelity level. By convention, $y_0(\cdot)$ and $y_L(\cdot)$ represent the lowest- and highest-fidelity functions, respectively. The autoregressive framework then establishes a relationship between two consecutive fidelity levels, given by:

$$y_{\kappa}(\mathbf{x}) = \rho_{\kappa} y_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}), \qquad (1)$$

where $\mathbf{x} \in D$. This model considers that $y_\kappa(\cdot)$ is a scaled version of the previous fidelity level $y_{\kappa-1}(\cdot)$ by a factor $\rho_\kappa \in \mathbb{R}$. Note that $\rho_\kappa$ is treated here as constant parameter to be estimated as in (Kennedy and O'Hagan, 2000), but it can be a functional, i.e., $\rho_\kappa : D \to \mathbb{R}$ (Le Gratiet, 2013; Meliani et al., 2019). Finally, the additive extra term $\gamma_\kappa : D \to \mathbb{R}$ is considered to capture the non-linear discrepancies between $y_\kappa(\cdot)$ and $y_{\kappa-1}(\cdot)$.

From a probabilistic perspective, one can place independent GP priors over the functions $y_0(\cdot), y_1(\cdot), ..., y_L(\cdot)$. Due to the linear structure of Eq. (1), it can be shown that, for all $\kappa \in \{1, ..., L\}$, $y_\kappa(\cdot)$ is also GP-distributed, with a kernel function that incorporates the autoregressive model. As demonstrated by (Kennedy and O'Hagan, 2000), Eq. (1) follows from the Markov assumption

$$\operatorname{cov}(y_{\kappa}(\mathbf{x}), y_{\kappa-1}(\mathbf{x}') | y_{\kappa-1}(\mathbf{x})) = 0, \qquad (2)$$

for all $\mathbf{x} \neq \mathbf{x}'$ and all $\kappa \in \{1, ..., L\}$. This property means that if $y_{\kappa-1}(\mathbf{x})$ is known, then nothing more can be learned about $y_\kappa(\mathbf{x})$ from any observation of the lower fidelity $y_{\kappa-1}(\mathbf{x}')$ for $\mathbf{x}' \neq \mathbf{x}$ (Le Gratiet, 2013). This is particularly
useful in scenarios where data from fidelity level $\kappa - 1$ provides supplementary information that is inaccessible at fidelity level $\kappa$, owing to real-world constraints such as cost and physical limitations that are typically absent in simulations, facilitating the incorporation of this additional knowledge to improve predictions derived from higher-fidelity observations.

The main limitation of the MFK framework proposed in (Kennedy and O'Hagan, 2000), lies in the covariance matrix inversion of the joint GP, which is required for both covariance estimation and prediction. Its matrix size scales with the total number of observations (i.e., training data). Let $\mathbf{X}_{\kappa} = (\mathbf{x}_{\kappa,i})_{1 \le i \le N_{\kappa}}$ denote the DoE at the $\kappa$-th fidelity level, containing $N_{\kappa} \in \mathbb{N}$ input points. The computational complexity of the resulting Kriging model is then $\mathcal{O}((N_0 + \cdots + N_L)^3)$.

Le Gratiet and Garnier (2014) have shown that the framework introduced in (Kennedy and O'Hagan, 2000) can be implemented more efficiently by considering nested DoEs, i.e., $\mathbf{X}_{\kappa} \subseteq \mathbf{X}_{\kappa-1}$. Under this assumption, a recursive formulation becomes possible. There, loosely speaking, the GP $y_{\kappa}(\cdot)$ is conditioned on the posterior of $y_{\kappa-1}(\cdot)$, which interpolates the dataset $(\mathbf{X}_{\kappa-1}, \mathbf{y}_{\kappa-1})$, where $\mathbf{y}_{\kappa-1} = y_{\kappa-1}(\mathbf{X}_{\kappa-1}) \in \mathbb{R}^{N_{\kappa-1}}$ is the observation vector of the $(\kappa - 1)$-th fidelity level. This recursive structure allows for a reduction in computational complexity to $\mathcal{O}(N_0^3 + \cdots + N_L^3)$.

Although the recursive MFK formulation offers improved computational efficiency, one of its main limitations lies in the requirement for nested DoEs. In practice, this condition is often difficult to satisfy. For instance, in aerodynamic applications, experimental and numerical setups cannot always be conducted under identical input conditions due to experimental safety constraints and numerical convergence issues, respectively.

As stated in the introduction, the HK framework proposed in (Han and Görtz, 2012) may be seen as an alternative for handling non-nested DoEs. However, it lacks uncertainty propagation from one fidelity to the other. This is because HK does not account for the posterior variance $\text{var}(y_{\kappa-1}(\cdot)|\mathbf{y}_{\kappa-1})$ in the Kriging-based prediction at the fidelity level $y_{\kappa}(\cdot)$. More precisely, the autoregressive model considered by HK is given by $y_{\kappa}(\mathbf{x}) = \rho_{\kappa} \hat{y}_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x})$ where $\hat{y}_{\kappa-1}(\cdot) = \mathbb{E}(y_{\kappa-1}(\cdot)|\mathbf{y}_{\kappa-1})$ is the point prediction of $y_{\kappa-1}(\cdot)$. Notably, the nested DoE constraint is no longer required in the HK construction, since $\hat{y}_{\kappa-1}(\cdot)$ can be evaluated at any input $\mathbf{x}$ once the model has been trained.

Both the recursive MFK and HK frameworks suffer from additional limitations, one of which pertains to the size of the DoE at each fidelity level. In particular, when large low-fidelity datasets are available from fast and inexpensive computer simulations, the scalability issues of standard Kriging-based frameworks are encountered, as the computational complexity of both models remains $\mathcal{O}(N_0^3)$ for the $\kappa$-th fidelity level due to the inversion of the corresponding covariance matrix. Another limitation is related to the unidirectional nature of knowledge transfer, which occurs only from $y_{\kappa-1}(\cdot)$ to $y_{\kappa}(\cdot)$, as the data $\mathbf{y}_{\kappa}$ are not incorporated when training $y_{\kappa-1}(\cdot)$.

In the next section, we revisit the multi-fidelity framework to address three key limitations: handling non-nested DoEs while allowing uncertainty propagation across fidelity levels and enabling scalability to large-scale datasets. To tackle the latter, in Section 4, we adopt sparse approximations for GPs (Quinonero-Candela and Rasmussen, 2005; Titisias, 2009a). We illustrate numerically that our model construction also supports bidirectional knowledge transfer, i.e., predictions at low-fidelity levels can be improved by incorporating information from high-fidelity data.

### 3. Multi-fidelity co-Kriging framework

As sparse approximations will be introduced in Section 4, it is convenient to incorporate a noise term in the autoregressive model, since training data can no longer be exactly interpolated. Including this noise term is both a standard and essential modeling assumption in Kriging, as in many practical scenarios, observations can be inherently noisy due to measurement errors, numerical approximations, or system uncertainties (Rasmussen and Williams, 2005).

Under this setting, the autoregressive scheme in Eq. (1) must be reformulated as follows:

$$f_{\kappa}(\mathbf{x}) = \rho_{\kappa} f_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}), \qquad (3)$$

$$y_{\kappa}(\mathbf{x}) = f_{\kappa}(\mathbf{x}) + \epsilon_{\kappa},$$

where $\epsilon_{\kappa} \sim \mathcal{N}(0, \tau_{\kappa}^2)$, with noise variance $\tau_{\kappa}^2 > 0$, represents the additive Gaussian noise term of the $\kappa$-th fidelity level. Notice that, in Eq. (3), it is assumed that the autoregressive framework is driven by latent functions $f_{\kappa}(\cdot)$, while only noisy observations $y_{\kappa}(\cdot)$ are available.
## 3.1. Co-Kriging framework

Before introducing the new formulation of the multi-fidelity framework, we first rewrite the autoregressive scheme in Eq. (3) as a single expression. This will allow simplification in the computation of cross-covariance functions across fidelity levels.

We first denote the lowest-fidelity latent function $f_0(\cdot)$ by $\gamma_0(\cdot)$. Using this notation, then the $\kappa$-th fidelity function $f_\kappa(\cdot)$ can be expressed recursively in terms of the latent functions $\gamma_0(\cdot), \dots, \gamma_\kappa(\cdot)$ as follows (see Appendix A.1):

$$f_{\kappa}(\mathbf{x}) = \sum_{j=0}^{\kappa} \eta_{j, \kappa} \gamma_j(\mathbf{x}), \qquad (4)$$

with

$$\eta_{j, \kappa} = \begin{cases} \prod_{i=j+1}^{\kappa} \rho_i, & \text{if } j < \kappa, \\ 1, & \text{if } j = \kappa. \end{cases} \qquad (5)$$

It can be easily verified from Eq. (4) that $f_0(\mathbf{x}) = \gamma_0(\mathbf{x})$ when $\kappa = 0$, and that $f_1(\mathbf{x}) = \rho_1 f_0(\mathbf{x}) + \gamma_1(\mathbf{x})$ when $\kappa = 1$. Both expressions correspond to the formulation used in a two-fidelity framework. By replacing Eq. (4) in the autoregressive scheme in Eq. (3), we have that

$$y_{\kappa}(\mathbf{x}) = \sum_{j=0}^{\kappa} \eta_{j, \kappa} \gamma_j(\mathbf{x}) + \epsilon_{\kappa}. \qquad (6)$$

Now, following (Kennedy and O'Hagan, 2000), we place GP priors over the functions $\gamma_0(\cdot), \dots, \gamma_L(\cdot)$. We assume that the GPs $\gamma_\kappa(\cdot)$ are independent, that is, $\gamma_\kappa(\cdot) \perp \gamma_{\kappa'}(\cdot)$ for all $\kappa \neq \kappa'$, with $\kappa, \kappa' \in \{0, \dots, L\}$. Each GP prior $\gamma_\kappa \sim GP(0, k_{\gamma\kappa})$ is assumed to be centered and defined by a kernel function $k_{\gamma\kappa}(\mathbf{x}, \mathbf{x}') := \text{cov}(\gamma_\kappa(\mathbf{x}), \gamma_\kappa(\mathbf{x}'))$. In addition, we assume that the independent and identically distributed (i.i.d.) Gaussian noise terms $\epsilon_0, \dots, \epsilon_L$ are independent of $\gamma_\kappa(\cdot)$ for all $\kappa \in \{0, \dots, L\}$.

From Eq. (6), it follows that the functions $y_\kappa(\cdot)$ are centered GPs. Therefore, to fully characterize the joint distribution of $y_0(\cdot), \dots, y_L(\cdot)$, it suffices to compute the cross-covariance functions between $y_\kappa(\cdot)$ and $y_{\kappa'}(\cdot)$ for all $\kappa, \kappa' \in \{0, \dots, L\}$. These covariances adhere to the following general structure:

$$c_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') := \text{cov}(y_{\kappa}(\mathbf{x}), y_{\kappa'}(\mathbf{x}')) = k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') + \tau_{\kappa}^2 \delta_{\kappa, \kappa'} \delta_{\mathbf{x}, \mathbf{x}'}, \qquad (7)$$

where

$$k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') := \text{cov}(f_{\kappa}(\mathbf{x}), f_{\kappa'}(\mathbf{x}')) = \sum_{j=0}^{\min(\kappa, \kappa')} \eta_{j, \kappa} \eta_{j, \kappa'} k_{\gamma_j}(\mathbf{x}, \mathbf{x}'), \qquad (8)$$

with $\mathbf{x}, \mathbf{x}' \in D$. Here, $\delta_{a,b}$ is the Kronecker delta function that is equal to 1 if $a = b$, and 0 otherwise. Note that when $\kappa = \kappa' = 0$, we obtain $\text{cov}(f_0(\mathbf{x}), f_0(\mathbf{x}')) = \eta_{0,0}^2 k_{\gamma_0}(\mathbf{x}, \mathbf{x}') = k_{\gamma_0}(\mathbf{x}, \mathbf{x}')$, which is consistent with our convention $f_0(\cdot) = \gamma_0(\cdot)$. The entire development to show Eq. (7) and Eq. (8) is detailed in Appendix A.2.

Note that this framework assumes that all fidelity levels are correlated, with the correlation structure defined by Eq. (7), which incorporates the autoregressive model of Eq. (3). As a result, our model establishes a multi-fidelity co-Kriging (MFCK) framework, which jointly models multiple variables. In particular, from Eq. (8), we observe that the resulting MFCK exhibits the structure of a multi-task GP (Bonilla, Chai and Williams, 2007), which corresponds to a particular class of MoGP commonly used in the machine learning literature. This similarity suggests that the MFCK formulation can leverage advanced GP-based methodologies for sparse approximation and variational inference (Álvarez and Lawrence, 2009; Van Der Wilk, 2019; López-Lopera, Idier, Rohmer and Bachoc, 2022), as we will show in Section 4. Here, unlike in (Bonilla et al., 2007), where output correlations are modeled through coregionalization parameters, the dependencies in our framework arise naturally from the autoregressive structure defined in Eq. (3).
## 3.2. MFCK predictions

For prediction purposes, we aim to condition the GPs $y_0(\cdot), \dots, y_L(\cdot)$ on observed data. For each fidelity level $\kappa$, let $X_\kappa = (x_{\kappa,i})_{1\le i \le N_\kappa}$ denote the set of input points, and define $y_\kappa := y_\kappa(X_\kappa) \in \mathbb{R}^{N_\kappa}$ as the observation vector.

Assuming the autoregressive structure in Section 3.1, we can now write the joint distribution of the concatenated vector $\mathbf{y} = [\mathbf{y}_0, \dots, \mathbf{y}_L]^\mathrm{T}$ as:

$$ \mathbf{y} \sim \mathcal{N}\left(\begin{bmatrix} \mathbf{0} \\ \vdots \\ \mathbf{0} \end{bmatrix}, \mathbf{K} + \begin{bmatrix} \tau_0^2 \mathbf{I}_{0,0} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \tau_L^2 \mathbf{I}_{L,L} \end{bmatrix}\right), \qquad (9) $$

where $\mathbf{K} = (K_{\kappa,\kappa'})_{0\le\kappa,\kappa'\le L}$ denotes the full block-wise covariance matrix of the training data, with each block $K_{\kappa,\kappa'} = k_{\kappa,\kappa'}(X_\kappa, X_{\kappa'})$ corresponding to the $N_\kappa \times N_{\kappa'}$ matrix containing the covariance between the input sets $X_\kappa = (x_{\kappa,i})_{1\le i \le N_\kappa}$ and $X_{\kappa'} = (x_{\kappa',i})_{1\le i \le N_{\kappa'}}$. The matrices $I_{\kappa,\kappa'}$ denote identity matrices of sizes $N_\kappa$.

Without loss of generality, consider the prediction of $y_\kappa(\cdot)$ at a new input point $\mathbf{x} \in \mathcal{D}$ and for any $\kappa \in \{0, \dots, L\}$. Following the standard GP conditioning procedure (Rasmussen and Williams, 2005), the conditional distribution $p(y_\kappa(\mathbf{x})|\mathbf{y})$ is Gaussian with mean and variance given by:

$$ m_{\kappa}(\mathbf{x}) = k_{\kappa}^{\mathrm{T}}(\mathbf{x})(\mathbf{K} + \mathbf{D})^{-1}\mathbf{y}, \qquad (10) $$

$$ v_{\kappa}(\mathbf{x}) = k_{\kappa,\kappa}(\mathbf{x}, \mathbf{x}) + \tau_{\kappa}^{2} - k_{\kappa}^{\mathrm{T}}(\mathbf{x})(\mathbf{K} + \mathbf{D})^{-1}k_{\kappa}(\mathbf{x}), $$

where $k_\kappa(\mathbf{x}) = [k_{\kappa,0}(\mathbf{x}, X_0), \dots, k_{\kappa,L}(\mathbf{x}, X_L)]^\mathrm{T}$ with $k_{\kappa,\kappa'}(\cdot, \cdot)$ as in Eq. (8), and $\mathbf{D} = \mathrm{diag}(\tau_0^2 \mathbf{I}_{0,0}, \dots, \tau_L^2 \mathbf{I}_{L,L})$. These expressions can be evaluated for each input point $\mathbf{x}_\kappa \in \mathcal{D}$ to obtain a prediction of $y_\kappa(\mathbf{x}_\kappa)$ together with the corresponding confidence intervals.

## 3.3. Parameter estimation via maximum likelihood

The parameters of the kernels $k_{\gamma_\kappa}(\cdot, \cdot)$ (denoted here by $\theta_\kappa$ for all $\kappa \in \{0, \dots, L\}$), those of the autoregressive model (i.e., $\rho_1, \dots, \rho_L$), and the noise variances (i.e., $\tau_0^2, \dots, \tau_L^2$) can be estimated via maximum likelihood. In this case, the complete set of parameters is given by

$$ \theta = ((\theta_\kappa)_{0\le\kappa\le L}, (\rho_\kappa)_{1\le\kappa\le L}, (\tau_\kappa^2)_{0\le\kappa\le L}). $$

As for standard GPs, the marginal log-likelihood $\mathcal{L}(\theta) := \log p_\theta(\mathbf{y})$ is defined as

$$ \mathcal{L}(\theta) = -\frac{1}{2} \log |\mathbf{K}_{\theta} + \mathbf{D}_{\theta}| - \frac{1}{2} \mathbf{y}^{\mathrm{T}} (\mathbf{K}_{\theta} + \mathbf{D}_{\theta})^{-1} \mathbf{y} + c, \qquad (11) $$

where $\mathbf{K}_{\theta} = (K_{\kappa,\kappa'})_{0\le\kappa,\kappa'\le L}$, $c \in \mathbb{R}$ a constant term independent of $\theta$, and $|\cdot|$ denotes the determinant operator. Note that the dependence on the parameters $\theta$ has been made explicit: while $\mathbf{K}_{\theta}$ depends on the kernel parameters and the autoregressive coefficients, $\mathbf{D}_{\theta}$ encodes the noise variances associated with each fidelity level.

Unlike the recursive MFK formulation in (Le Gratiet and Garnier, 2014), the parameters $\theta$ must now be estimated jointly. As a result, the main limitation of this model (similarly to the framework in (Kennedy and O'Hagan, 2000)) lies in the block-wise inversion matrix $\mathbf{K}_{\theta}+\mathbf{D}_{\theta} \in \mathbb{R}^{N \times N}$, with $N = \sum_{\kappa=0}^{L} N_\kappa$. This computational burden can be alleviated by exploiting inversion properties of partitioned matrices (see, e.g., Appendix A.3 in Rasmussen and Williams, 2005), which enables more efficient implementations. Similar properties can be used to compute the likelihood gradients. An alternative, and more scalable solution, discussed later in Section 4, is to consider sparse approximations (Snelson and Ghahramani, 2005; Titisias, 2009a). There, briefly speaking, the computational complexity is controlled by a set of inducing points $M \in \mathbb{N}$ rather than $N$, such that $M \ll N$.

## 3.4. Numerical illustration

To illustrate the benefits of the MFCK framework, we consider the one-dimensional Forrester function, a benchmark commonly used in the literature to evaluate multi-fidelity models (Forrester, Keane and O'Hagan, 2007). The two-fidelity Forrester function is defined as follows:

$$ f_{HF}(x) = (6x - 2)^2 \sin(12x - 4), \qquad (12) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison between (a) the recursive MFK framework introduced in (Le Gratiet and Garnier, 2014), and (b) MFCK (ours). For both models, the target functions and observations are represented by dashed lines and points, respectively. The high-fidelity (HF) function $y_{HF}(\cdot)$ is shown in blue, while the low-fidelity (LF) one $y_{LF}(\cdot)$ is in orange. Using the same color code, the predictions are displayed as solid lines, with the corresponding confidence intervals (95%) indicated by shaded regions. This example considers a non-nested DoE with 4 HF observations (blue dots) and 8 LF observations (orange dots).</figcaption></figure>

$$f_{LF}(x) = \frac{1}{2} f_{HF}(x) + 10(x - 1),$$

with $x \in [0, 1]$, $f_{HF}(\cdot)$ denoting the HF level $f_L(\cdot)$ with $L = 1$, and $f_{LF}(\cdot)$ denoting the LF level $f_0(\cdot)$. These functions exhibit a complex behavior with multiple local minima, making it particularly challenging for testing multi-fidelity models (Yıldız, Pehlivan-Solak, Diez, Gören and Nikbay, 2022; Mainini, Serani, Rumpfkeil, Minisci, Quagliarella, Pehlivan, Yıldız, Ficini, Pellegrini, Fiore, Bryson, Nikbay, Diez and Beran, 2022).

Here, we seek to compare the MFCK predictions with those obtained using the recursive MFK formulation under interpolation assumptions. Therefore, we have limited the analysis to noise-free observations, which implies that $y_{HF}(\cdot) = f_{HF}(\cdot)$ and $y_{LF}(\cdot) = f_{LF}(\cdot)$. Regardless of this assumption, the noise variances in the MFCK model are still accounted for and estimated via maximum likelihood. Following estimation, the noise variances for both fidelity levels were found to be small, on the order of $10^{-5}$. For MFK, we have used the implementation available in the Python toolbox SMT (Bouhlel, Hwang, Bartoli, Lafage, Morlier and Martins, 2019).

Figure 1 shows the comparison between MFK and MFCK for a random non-nested DoE with $N_{HF} = 4$ and $N_{LF} = 8$. Both models are trained via maximum likelihood. We observe from Figure 1(a) that, as expected, the recursive MFK framework propagates undesired uncertainty from $y_{LF}(\cdot)$ to $y_{HF}(\cdot)$ at locations where $y_{HF}(\cdot)$ is observed. This results from a violation of nested DoE assumptions, under which the recursive predictive formulas are no longer valid. In particular, this issue comes from the conditional variance, which recursively transfers uncertainty through the fidelity hierarchy (Le Gratiet and Garnier, 2014). Such an undesired behavior may prevent the model from accurately capturing the HF behavior as observed in our example. This limitation tends to become more sensitive in higher-dimensional settings.

In contrast, from Figure 1(b), the MFCK framework successfully infers both $y_{LF}(\cdot)$ and $y_{HF}(\cdot)$, demonstrating its ability to learn the discrepancy between both fidelity levels. Moreover, unlike the recursive MFK that allows knowledge transfer only from the LF to HF levels (Le Gratiet and Garnier, 2014) due to the nested DoE assumption, MFCK enables bidirectional knowledge transfer. This means that information from the HF data can also improve predictions of the LF function in regions where it has not been directly observed (e.g., near $x = 0.3$ and $x = 0.8$ in our example). This bidirectional transfer learning property may be less relevant in multi-fidelity problems, where the primary interest lies in predicting the HF target function. However, it may be of greater importance in more general transfer learning settings, where the objective is to improve predictions of all target functions by leveraging information from the others.
## 4. Sparse approximations of the MFCK

The computational complexity of MFCK is $O(N^3)$ with $N = \sum_{\kappa=0}^{L} N_{\kappa}$, making the framework intractable for large datasets. To mitigate this limitation, we suggest adopting sparse approximations that reduce the complexity to $O(NM^2)$, where $M \ll N$ denotes the number of a smaller set of points called inducing points. For this purpose, we first need to introduce the notion of inducing variables.

Let $Z_{\kappa} = (z_{\kappa,i})_{1 \le i \le M_{\kappa}}$ be the set of inducing points for the $\kappa$-th fidelity, where $M_{\kappa} \in \mathbb{N}$ such that $M_{\kappa} \ll N_{\kappa}$. We define the inducing variables as:

$$ \mathbf{u}_{\kappa} = f_{\kappa}(Z_{\kappa}), \quad \forall \kappa \in \{0, \dots, L\}. \qquad (13) $$

As $\mathbf{u}_{\kappa}$ is a vector composed of evaluations of $f_{\kappa}(\cdot)$ at inducing points $Z$, then it is Gaussian distributed. Denote $\mathbf{f} = [f_0, \dots, f_L]^{\top}$ and $\mathbf{u} = [u_0, \dots, u_L]^{\top}$. The joint distribution of the function values and the inducing variables is given by

$$ \begin{bmatrix} \mathbf{f} \\ \mathbf{u} \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \begin{bmatrix} \mathbf{K}_{f,f} & \mathbf{K}_{f,u} \\ \mathbf{K}_{f,u}^{\top} & \mathbf{K}_{u,u} \end{bmatrix} \right), \qquad (14) $$

where

$$ \begin{align} \mathbf{K}_{f,f} &= (k_{\kappa,\kappa'}(\mathbf{X}_{\kappa}, \mathbf{X}_{\kappa'}))_{0 \le \kappa, \kappa' \le L}, \notag \\ \mathbf{K}_{f,u} &= (k_{\kappa,\kappa'}(\mathbf{X}_{\kappa}, \mathbf{Z}_{\kappa'}))_{0 \le \kappa, \kappa' \le L}, \tag{15} \\ \mathbf{K}_{u,u} &= (k_{\kappa,\kappa'}(\mathbf{Z}_{\kappa}, \mathbf{Z}_{\kappa'}))_{0 \le \kappa, \kappa' \le L}, \notag \end{align} $$

are block-wise covariance matrices constructed from kernel evaluations $k_{\kappa,\kappa'}(\cdot, \cdot)$ (see definition in Eq. (8)). Note that here we write the indices in the covariance matrices to distinguish if the covariance is computed between input and inducing points.

### 4.1. Fully Independent Training Conditional (FITC) approximation

**Framework.** In FITC, we need to approximate the conditional prior $p(\mathbf{f}|\mathbf{u})$ by the factorized distribution:

$$ p(\mathbf{f}|\mathbf{u}) \approx \prod_{\kappa=0}^{L} \prod_{i=1}^{N_{\kappa}} p(f_{\kappa}(\mathbf{x}_i)|\mathbf{u}). \qquad (16) $$

That is, each $f_{\kappa}(\mathbf{x}_i)$ is assumed to depend exclusively on $\mathbf{u}$, with no direct dependence on other outputs $f_{\kappa'}(\mathbf{x}_j)$ for $j \ne i$ and $\kappa' \in \{0, \dots, L\}$ (Quinonero-Candela and Rasmussen, 2005). This allows us to construct a matrix $\mathbf{Q}$ that will approximate the true covariance matrix $\mathbf{K}_{f,f}$ using the Nyström approximation:

$$ \mathbf{Q} = \mathbf{K}_{f,f} \mathbf{K}_{u,u}^{-1} \mathbf{K}_{u,f}, \qquad (17) $$

with $\mathbf{K}_{u,f} = \mathbf{K}_{f,u}^{\top}$. In contrast to the FITC approximation used in sparse Kriging, the three covariance matrices in Eq. (17) exhibit block-wise structures that capture the cross-covariances among all fidelity levels (see Eq. (15)).

As suggested by (Quinonero-Candela and Rasmussen, 2005), beyond the Nyström approximation, we also apply a correction to the diagonal by incorporating the exact diagonal covariance matrix elements $\mathbf{K}_{f,f}$:

$$ \tilde{\mathbf{Q}} = \mathbf{Q} + \mathrm{diag}(\mathbf{K}_{f,f} - \mathbf{Q}). \qquad (18) $$

This matrix $\tilde{\mathbf{Q}}$ will serve to enjoy computationally tractable implementations as shown later in Eq. (??).

**Parameter estimation.** The marginal likelihood $p(\mathbf{y})$ is obtained by marginalizing the conditional distribution $p(\mathbf{y}|\mathbf{f})$ with respect to $\mathbf{f}$. As in standard Kriging, the posterior distribution is given by $p(\mathbf{y}|\mathbf{f}) = \mathcal{N}(\mathbf{f}, \mathbf{D})$ with $\mathbf{D} = \mathrm{diag}(\tau_0^2 \mathbf{I}_{0,0}, \dots, \tau_L^2 \mathbf{I}_{L,L})$. Note that, in contrast to standard SK, the matrix $\mathbf{D}$ in this context exhibits a block-diagonal
structure. Considering that $p(f) = \mathcal{N}(0, \tilde{Q})$ (FITC approximation), then marginalizing $p(y|f)$ over $f$ yields a Gaussian distribution. Consequently, the marginal likelihood $p(y)$ is also Gaussian-distributed:

$$p(y) = \mathcal{N}(0, \tilde{Q} + D). \qquad (19)$$

Equation (19) can be used for parameter estimation. Recall that MFCK parameters, denoted by $\theta$ in Section 3.3, are involved in both $\tilde{Q}$ and $D$. For notational simplicity, we omit the explicit dependence on $\theta$ in the following computations.

The maximization of $p(y)$ (or $\log p(y)$) involves the inversion of the matrix $\tilde{Q} + D$, which can be computed efficiently using the matrix inversion lemma (see, Appendix A.3 in Rasmussen and Williams, 2005). Define $\Lambda = \text{diag}(K_{f,f} - Q + D)$, then we have:

$$(\tilde{Q} + D)^{-1} = (K_{f,u} K_{u,u}^{-1} K_{u,f} + \Lambda)^{-1} = \Lambda^{-1} - \Lambda^{-1} K_{f,u} P^{-1} K_{u,f} \Lambda^{-1},$$

with $P = K_{u,u} + K_{u,f} \Lambda^{-1} K_{f,u}$. This reduces the complexity as inverting this new approximated covariance matrix scales as $O(NM^2)$ with $N = \sum_{\kappa=0}^{L} N_{\kappa}$ and $M = \sum_{\kappa=0}^{L} M_{\kappa}$.

Note that the new tractable objective function in Eq. (19) is derived by modifying the prior, replacing the exact model with a Nyström approximation. As a result, the optimized model only approximates the exact GP for $M \ll N$, which may lead to a degradation in predictive accuracy. However, as shown in (Snelson and Ghahramani, 2005), the FITC approximation still provides accurate performance in regression tasks.

**Prediction.** Recall that the prior considered in FITC replaces the real covariance matrix $K_{f,f}$ using the corrected Nyström approximation in Eq. (18). The same approximation is used to replace the cross-covariance between training points and test points. Following the notation in Section 3.2, we consider prediction at a new input point $x \in D$ for any fidelity level $\kappa \in \{0, \dots, L\}$. The joint distribution of $(f, f_{\kappa}(x))$ is then given by

$$p(f, f_{\kappa}(x)) \sim \mathcal{N} \left( \mathbf{0}, \begin{bmatrix} \tilde{Q} & q_{\kappa}(x) \\ q_{\kappa}^{\top}(x) & k_{\kappa, \kappa'}(x, x) \end{bmatrix} \right), \qquad (20)$$

where $q_{\kappa}(x) = [q_{\kappa,0}(x, X_0), \dots, q_{\kappa,L}(x, X_L)]^\top$ with $q_{\kappa,\kappa'}(x, X_{\kappa'}) = k_{\kappa,\kappa'}(x, Z_{\kappa'}) k_{\kappa',\kappa'}^{-1}(Z_{\kappa'}, Z_{\kappa'}) k_{\kappa,\kappa'}(Z_{\kappa'}, X_{\kappa'})$. Note that $q_{\kappa,\kappa'}(x, X_{\kappa'})$ serves as a proxy of $k_{\kappa,\kappa'}(x, X_{\kappa'})$.

When predicting $y_{\kappa}(\cdot)$, we compute the conditional distribution $p(y_{\kappa}(x)|y)$, which is Gaussian with the following mean and variance:

$$m_{\kappa}(\mathbf{x}) = q_{\kappa}^{\top}(\mathbf{x})(\tilde{Q} + D)^{-1}\mathbf{y}, \qquad (21)$$

$$v_{\kappa}(\mathbf{x}) = k_{\kappa, \kappa}(\mathbf{x}, \mathbf{x}) + \tau_{\kappa}^{2} - q_{\kappa}(\mathbf{x})^{\top}(\tilde{Q} + D)^{-1}q_{\kappa}(\mathbf{x}).$$

Observe that the predictive mean is obtained by projecting the observations through the inducing variables, while the predictive variance retains the exact prior variance $k_{\kappa,\kappa}(\mathbf{x}, \mathbf{x})$ but subtracts the explained component through the inducing space.

**Remark.** We note that other sparse methods based on low-rank approximations, e.g., the SoR approximation and the deterministic training conditional (DTC) (Quinonero-Candela and Rasmussen, 2005), can also be adapted to the multi-fidelity setting in a similar fashion. However, in this work, we focus exclusively on FITC, as it preserves the variances of the exact conditional while maintaining independence between training points, which is essential for computational tractability (Van Der Wilk, 2019).

## 4.2. Variational Free Energy (VFE) approximation

**Framework.** Consider now the computation of the posterior distribution $p(u|y) = \int p(u|f)p(f|y)d f$, where $p(u|f)$ denotes the conditional prior over any finite set of function points $u$. Instead of modifying the exact Kriging model and optimizing the resulting modified marginal likelihood as in FITC, the VFE approximation focuses on minimizing the KL divergence between the exact GP posterior and a variational approximation (Titsias, 2009a).
We need to consider a variational distribution $q(f, u)$ to approximate the true posterior $p(f, u|y)$. As the exact posterior can be factorized as $p(f, u|y) = p(f|u) p(u|y)$, it is natural to impose a similar factorization for the variational distribution:

$$q(f, u) = p(f|u) q(u). \qquad (22)$$

Since the joint distribution of $(f, u)$ is Gaussian, the conditional distribution $p(f|u)$ can be derived directly using standard properties of multivariate Gaussian vectors:

$$p(f|u) = \mathcal{N} \left( K_{f,u} K_{u,u}^{-1} u, K_{f,f} - Q \right),$$

with $Q$ as in Eq. (17).

Notice that the variational distribution $q(u)$ is a free-form multivariate Gaussian introduced to approximate the true posterior over the inducing variables. This construction preserves the exact GP prior conditional $p(f|u)$, while allowing flexibility in the marginal $q(u)$ to capture posterior uncertainty (Titsias, 2009b). The key assumption underlying this formulation is that the dependence between $f$ and $u$ is fully governed by the prior, such that $q(f|u) = p(f|u)$, and all variational flexibility is allocated to $q(u)$. The variational distribution $q(u)$ is then introduced to approximate $p(f, u|y)$ ensuring that $q(f, u) = p(f|u)q(u)$.

**Parameter estimation.** As in (Titsias, 2009a), the covariance parameter estimation use a new variational distribution $q(\cdot)$ to write the log marginal likelihood $p(y)$ by making appear the variational distribution $q(u)$ (see Appendix B).

This allows us to rewrite the log marginal likelihood as an expected value over the distribution $q(f, u)$.

Using Jensen's inequality, we can rewrite the bound:

$$\log p(y) \leq \mathbb{E}_{q(f,u)} \left( \log p(y|f) \frac{p(u)}{q(u)} \right) = \int q(u) \log \left( \mathcal{N}(K_{f,u} K_{u,u}^{-1} u, D) \frac{p(u)}{q(u)} \right) du - \frac{1}{2} \mathrm{tr} (D^{-1} [K_{f,f} - Q]). \qquad (23)$$

Equation (23) can be derived following a similar approach to that in (Titsias, 2009a), with the key difference that the noise variance cannot be factorized out from the trace, as the matrix $D$ is block-diagonal and incorporates a distinct noise variance parameter for each fidelity level. The entire development to obtain Eq. (23) is provided in Appendix B.

The right-hand side of Eq. (23), denoted here as $\mathcal{L}(\theta, q)$, is known as the variational lower bound on the true log marginal likelihood and can be optimized to infer model parameters. After applying Jensen's inequality and moving the logarithm outside the integral, we obtain the expression that serves as the objective function for optimization:

$$\mathcal{L}(\theta) = \log [\mathcal{N}(\mathbf{0}, \mathbf{Q} + \mathbf{D})] - \frac{1}{2} \mathrm{tr} (\mathbf{D}^{-1} [\mathbf{K}_{f,f} - \mathbf{Q}]). \qquad (24)$$

This formulation can be expressed to show its dependence on the inducing inputs $Z_k$, which allows for their optimization (Titsias, 2009a). However, in this work, we assume that the inducing points are fixed (see discussion in Section 4.4) and are not subject to estimation as their optimization may result in computational overhead. We recall that in the multi-fidelity setting, we have $M_k$ inducing points for each fidelity level.

Similar to the FITC method, we can apply the matrix inversion lemma to invert $(Q + D)$ following the same procedure as in Eq. (??). Finally, as shown in Eq. (24) the main difference between the likelihoods obtained using the FITC and VFE methods lies in how they handle the diagonal correction of the low-rank approximation: FITC adjusts only the diagonal entries to match the exact covariance, while VFE introduces a global trace term that regularizes the approximation and ensures a valid lower bound on the marginal likelihood.

**Predictions.** By applying the first-order optimality condition to the variational bound $\mathcal{L}(\theta, q)$ in Eq. (23), we derive closed-form expressions for the optimal variational posterior distribution, specifically its conditional mean and covariance. This derivation reveals that the dependence of the variational distribution $q(u)$ on the observed data is fully captured through the inducing inputs $Z$. Consequently, the optimal $q(u)$ can be obtained analytically and does not require explicit parameterization or additional optimization. This quantity takes the form (see Appendix B for the detailed derivation):

$$q^*(u) = \mathcal{N}(u | K_{u,u} \Sigma^{-1} K_{u,y} D^{-1} y, K_{u,u} \Sigma^{-1} K_{u,u}), \qquad (25)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Comparison of the SMFCK strategies: (a) FITC and (b) VFE. Non-nested DoE with $N_1 = 100$, $N_0 = 200$, inducing points estimated via k-means with $M_1 = 4$ and $M_0 = 8$. The panel description is the same as in Figure 1. Bottom-aligned rectangles represent the inducing point locations for both HF (blue) and LF (orange) levels.</figcaption></figure>

with $\Sigma = K_{u,u} + K_{u,y}D^{-1}K_{y,u}$. Equation (25) is then used to establish the expressions for the conditional mean and variance. We can follow a similar procedure to the one proposed in (Titsias, 2009b) to obtain the mean and covariance expressions for $p(y_k(\mathbf{x})|\mathbf{y})$:

$$m_k(\mathbf{x}) = k_k(\mathbf{x}, Z)K_{u,u}^{-1}\Sigma^{-1}K_{u,y}\mathbf{y}, \qquad (26)$$

$$v_k(\mathbf{x}) = k_{\kappa,\kappa}(\mathbf{x}, \mathbf{x}) - k_k(\mathbf{x}, Z)K_{u,u}^{-1}k_k^T(\mathbf{x}, Z) + k_k^T(\mathbf{x}, Z)\Sigma^{-1}k_k(\mathbf{x}, Z),$$

where $k_k(\mathbf{x}, Z) = [k_{\kappa,0}(\mathbf{x}, Z_0), \dots, k_{\kappa,L}(\mathbf{x}, Z_L)]^T$.

### 4.3. Numerical illustration

We now illustrate the sparse MFCK (SMFCK) frameworks through the one-dimensional Forrester example described in Section 3.4. This experiment is designed to assess the different sparse approximation methods introduced earlier. Here, we increase the DoE size to have $N_1 = 100$ and $N_0 = 200$ to illustrate the benefit of accounting for sparsity. We aim to approximate the exact MFCK with a set of inducing points with $M_1 = 4$ and $M_0 = 8$ with the same point allocation considered in Section 3.4. To generate noisy data, we consider adding noise in each fidelity level, defining the noise variance per fidelity level as $\tau_1^2 = 0.2$ and $\tau_0^2 = 0.4$. These parameters, together with the model parameters, are estimated later via maximum likelihood.

Figure 2 illustrates the predictive performance of FITC and VFE approximations, each achieving up to a 96% reduction in the number of active data points by using only 4 and 8 inducing points, respectively, to approximate exact models constructed from 100 and 200 real observations. Despite this substantial reduction in computational cost, both methods maintain high predictive accuracy, highlighting the capability of inducing point approximations to effectively capture the data underlying structure.

In addition, Figure 2 also shows that the confidence intervals provided by SMFCK-FITC are narrower than those obtained with SMFCK-VFE, even though both models are evaluated under the same experimental settings. These trends are in agreement with the comparative analysis of sparse approximations provided in (Bauer, van der Wilk and Rasmussen, 2016). This study has revealed that the FITC method often overestimates the marginal likelihood and tends to underestimate the noise variance. This behavior can lead to a mischaracterization of uncertainty and a poor posterior distribution approximation. In contrast, the VFE method tends to overestimate the noise variance, but generally succeeds in recovering the true posterior distribution more accurately.

### 4.4. Selection of inducing points

As discussed in (Quinonero-Candela and Rasmussen, 2005; Titsias, 2009a; Bauer et al., 2016), it is crucial to consider an adequate strategy to allocate the inducing points correctly. One might consider a random sampling strategy; however, such an approach does not guarantee coverage of the input domain nor account for data density. As a result,
a poorly chosen random set of inducing points may lead to suboptimal predictive performance. This issue becomes more sensitive in high dimensions.

An alternative approach is to consider data-driven clustering methods. As demonstrated numerically in (Hensman, Fusi and Lawrence, 2013; Van Der Wilk, 2019; Valayer, Bartoli, Castaño-Aguirre, Lafage, Lefebvre, López-Lopera and Mouton, 2024), k-means offers a data-adaptive strategy that balances spatial coverage with computational efficiency by concentrating inducing points in informative input space regions. This typically leads to more stable and accurate performance. For this reason, k-means is used in further experiments in Section 5. An alternative clustering method, such as DBSCAN, was also investigated in (Valayer et al., 2024), but its performance was still inferior to that of k-means.

To perform k-means for inducing point selection at each fidelity level, we must specify the number of clusters, which corresponds to the number of inducing points $M_k$ for each $\kappa \in \{0, ..., L\}$. These values can be chosen based on empirical rules. In (Burt, Rasmussen and van der Wilk, 2019), the convergence of different criteria is studied when increasing the number of inducing points. Here, we propose to follow an empirical rule based on the experiments of (Burt et al., 2019), which suggests that $M_k$ can scale with the input dimension $D$ and the number of training points as follows:

$$M_k = [cD \log(N_k)], \qquad (27)$$

where $[\cdot]$ denotes the ceiling operator, which returns the smallest integer greater than or equal to its argument, and $c \in [1, 10]$ is a constant that depends on the kernel smoothness and the function complexity. Recall that $N_k$ denotes here the number of training points at the $\kappa$-th fidelity level.

In all experiments discussed later in Section 5, we apply Algorithm 1 based on the empirical rule of Eq. (27) for selecting $M_k$:

$$M_k = \{[D \log(N_k)], ..., [10D \log(N_k)]\}.$$

For the sake of simplicity, the update of the inducing points $Z_0, ..., Z_L$ is not stated in Algorithm 1, but they need to be computed via k-means once $M_0, ..., M_L$ are updated.

The adaptive inducing-point selection in Algorithm 1 seeks to achieve finite termination and monotonic improvement of the objective function. This algorithm has a convergence support according to (Burt et al., 2019) for the VFE-based formulation. There, the convergence to a local optimum of the likelihood optimization is proved as more inducing points are added. While convergence to the exact GP is not guaranteed due to the imposed inducing-point budget, the method adaptively identifies a compact set of inducing points that captures the effective model complexity at each fidelity level.

## 5. Results

We now present the results of the models introduced in Sections 3 and 4, namely MFCK and SMFCK, respectively. For the latter, we consider both FITC and VFE sparse approximations. Implementations of both models have been integrated into the open-source Python toolbox SMT-release 2.10.1 (Saves, Lafage, Bartoli, Diouane, Bussemaker, Lefebvre, Hwang, Morlier and Martins, 2024; Bouhlel et al., 2019). The source code and example notebooks are available in the GitHub repository: https://github.com/SMTorg/smt. Additionally, the benchmark experiments presented in this section are available in https://github.com/SMTorg/smt/tree/master/tutorial/MultiFi. For the comparison with MFK, we use the implementation already available in SMT. The comparison with HK is conducted using the open-source python codes linked to (Yang, Kang and Yee, 2023).

For all Kriging implementations, including the numerical illustration in Section 4.3, we consider Squared Exponential (SE) kernels defined as

$$k_{\theta}=(\sigma^{2}, \theta_{1}, \ldots, \theta_{D})\left(\mathbf{x}, \mathbf{x}^{\prime}\right)=\sigma^{2} \prod_{i=1}^{D} \exp \left(-\theta_{i}\left(x_{i}-x_{i}^{\prime}\right)^{2}\right),$$

where $\sigma^2 > 0$ is the variance parameter and $\theta_i > 0$ denotes the scaling parameter of the $i$-th input dimension.

Our goal throughout this section is to assess the performance of these models across different dimensional settings. We first consider a 2D example based on the Rosenbrock function. This function is commonly used in multi-fidelity
<caption>**Algorithm 1 Selection of inducing points for SMFCK**</caption>

**Require:** Initial number of inducing points (IPs)
$M_{\kappa} = D \log(N_{\kappa})$ for all $\kappa \in \{0, ..., L\}$
**Require:** Maximum number of IPs $M_{\kappa}^{\max} = 10D \log(N_{\kappa})$
**Require:** Performance metric $\mathcal{L} := \mathcal{L}(\theta)$
Define $\mathcal{L}^{M_{\kappa}}$ as the likelihood of the $\sum_{k=0}^{L} M_{\kappa}$ IPs
**Require:** Tolerance $\epsilon = 10^{-6}$

▷ Inducing point selection per fidelity level

1: **for** $\kappa = 0, ..., L$ **do**

2: $\Delta = +\infty$

3: Train the model with likelihood $\mathcal{L}^{M_{\kappa}}$

4: **while** $\Delta > \epsilon$ and $M_{\kappa} < M_{\kappa}^{\max}$ **do**

5: Train model with inducing points $\mathcal{L}^{M_{\kappa}+1}$

6: $\Delta = |\mathcal{L}^{M_{\kappa}+1} - \mathcal{L}^{M_{\kappa}}|$

7: $M_{\kappa} = M_{\kappa} + 1$

8: $\mathcal{L}^{M_{\kappa}} = \mathcal{L}^{M_{\kappa}+1}$

9: **end while**

10: **end for**

11: **return** Final number of inducing sets $M_0, ..., M_L$

studies because of its nonlinearity, narrow curved valley, and scalability in dimensionality, which make it a challenging benchmark for assessing the accuracy and robustness of multi-fidelity surrogate models (Mainini et al., 2022). We then extend the analysis to a 3D Large Reference Model (LRM), which incorporates both WT and CFD data (see further details in Section 5.2).

Throughout this section, we adopt the convention of denoting the LF and HF levels by $y_{LF}(\cdot) := y_0(\cdot)$ and $y_{HF}(\cdot) := y_L(\cdot)$, respectively. Accordingly, $N_{LF}$ and $N_{HF}$ represent the number of observations available at the LF and HF levels, while $M_{LF}$ and $M_{HF}$ denote the corresponding number of inducing points. Both $M_{LF}$ and $M_{HF}$ are chosen by following the procedure detailed in Section 4.4.

## 5.1. 2D Rosenbrock example

The Rosenbrock function is a well-known D-dimensional function used in a multi-fidelity context (Yıldız et al., 2022; Mainini et al., 2022). For the two-fidelity case, the HF and LF functions are given by

$$f_{HF}(\mathbf{x}) = \sum_{i=1}^{D-1} 100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2,$$

$$f_{LF}(\mathbf{x}) = \sum_{i=1}^{D-1} 50(x_{i+1} - x_i^2)^2 + (2 + x_i)^2 - \sum_{i=1}^{D} 0.5x_i.$$

We generate data from these functions considering $D = 2$, and noise variances $\tau_{HF}^2 = 0.2$ and $\tau_{LF}^2 = 0.4$.

As demonstrated in the numerical examples using the Forrester functions, practical scenarios often involve data collected independently at different fidelity levels, typically without shared sampling locations. To reflect this realistic setting, we construct a non-nested DoE for the Rosenbrock function by generating independent Latin Hypercube Samples (LHS) for each fidelity level. Figure 3 illustrates this setup for the case where $(x_1, x_2) \in [-2, 2]^2$. Specifically, Figure 3(a) shows the DoE for the HF and LF functions, and Figure 3(b) shows the distribution of the corresponding inducing point locations via k-means clustering.

Figure 4 compares the original Rosenbrock HF function with predictions using SMFCK-FITC for two different inducing points configurations. Figure 4(b) shows the results for $M_{HF} = 8$ and $M_{LF} = 16$, showing a less accurate prediction. In contrast, Figure 4(c) demonstrates improved prediction performance with $M_{HF} = 20$ and $M_{LF} = 30$. This result underscores the importance of selecting appropriate inducing points when approximating the MFCK model, as discussed in Section 4.4.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Example of DoE for two-dimensional Rosenbrock functions with two fidelity levels for non-nested DoE. Panels show: (a) the original data distribution and (b) the estimated inducing points using *k*-means.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Comparison of predictions for the two-fidelity 2D Rosenbrock function using the SMFCK-FITC. The DoE includes $N_{HF} = 128$ and $N_{LF} = 256$. The panels show: (a) the ground truth HF function, (b) prediction with $M_{HF} = 8$ and $M_{LF} = 16$ (RMSE = 12.1), and (c) prediction with $M_{HF} = 20$ and $M_{LF} = 30$ (RMSE = 0.6).</figcaption></figure>

For a comprehensive comparison, we repeat the experiments using the Rosenbrock functions ten times for different values of $N_{HF}$ and $N_{LF}$ while fixing the inducing points based on Algorithm 1 applied over the DoE with more observations (i.e, $N_{LF} = 1024$, $N_{HF} = 512$), obtaining $M_{HF} = 20$ and $M_{LF} = 30$. In each replicate, we have independently generated LHS-based DoEs.

The mean and standard deviation of the RMSE are reported in Table 1 for HK, MoGP, MFCK, SMoGP, SMFCK-FITC, and SMFCK-VFE. The table shows that the exact MFCK model consistently achieves the lowest RMSE across different DoE sizes. Training times for models with $O(N^3)$ complexity, such as HK and MoGP, may differ due to the optimization strategies used for parameter tuning (Álvarez and Lawrence, 2008; Yang et al., 2023). On the other hand, the sparse approximations, SMoGP, SMFCK-FITC, and SMFCK-VFE, yield higher RMSE values but offer a significant reduction in computational cost, with SMFCK-FITC showing better balance between predictive accuracy and training time. In particular, sparse methods reduce the training time by approximately 99% compared to MFCK, making them attractive alternatives in scenarios where computational efficiency is critical. As shown in Table 1, the SMFCK model outperforms more general approaches such as SMoGP due to its autoregressive assumption. Unlike SMoGP, which searches the correlation between different fidelities in a more general way, SMFCK explicitly models the dependency between fidelity levels through an autoregressive structure. When the underlying function is compatible with this assumption (see Eq. (1)), this leads to superior predictive performance as well as reduced training time.
**Table 1**

Comparison of RMSE (mean ± standard deviation for 10 replicates) and training time for HK, MoGP, MFCK, SMoGP, and SMFCK (FITC and VFE) on the 2D Rosenbrock function using three different non-nested DoE configurations. For SMFCK and SMoGP, we fix $M_{HF} = 20$ and $M_{LF} = 30$ using Algorithm 1 applied over the DoE, $N_{LF} = 1024$, $N_{HF} = 512$. Inducing points selection is conducted via k-means. The smallest RMSE values and training times for each input configuration are highlighted in bold. Because HK implementations rely on a different optimization scheme with higher iteration requirements, a fair comparison of training times is not possible. Therefore, those results are not reported.

<table><thead><tr><th rowspan="2">N<sub>HF</sub></th><th rowspan="2">N<sub>LF</sub></th><th rowspan="2">Metric</th><th colspan="3">Exact Methods</th><th colspan="3">Sparse Methods</th></tr><tr><th>HK</th><th>MoGP</th><th>MFCK (ours)</th><th>SMoGP</th><th>SMFCK-FITC (ours)</th><th>SMFCK-VFE (ours)</th></tr></thead><tbody><tr><td rowspan="2">32</td><td rowspan="2">64</td><td>RMSE (HF)</td><td>1.87±0.89</td><td>4.19±1.25</td><td><strong>1.36 ± 0.78</strong></td><td>9.22±2.21</td><td>2.66 ± 0.79</td><td>5.80 ± 0.81</td></tr><tr><td>Training time (s)</td><td>–</td><td>0.55±0.12</td><td>1.19 ± 0.02</td><td>1.80±0.82</td><td><strong>0.16 ± 0.03</strong></td><td>0.18 ± 0.03</td></tr><tr><td rowspan="2">128</td><td rowspan="2">256</td><td>RMSE (HF)</td><td>0.34±0.10</td><td>0.22±0.14</td><td><strong>0.22 ± 0.11</strong></td><td>0.87±0.19</td><td>0.65 ± 0.28</td><td>0.70 ± 0.36</td></tr><tr><td>Training time (s)</td><td>–</td><td>1.99±1.24</td><td>8.08 ± 0.19</td><td>5.74±1.74</td><td><strong>0.20 ± 0.01</strong></td><td>0.25 ± 0.01</td></tr><tr><td rowspan="2">512</td><td rowspan="2">1024</td><td>RMSE (HF)</td><td>0.25±0.19</td><td>0.09±0.01</td><td><strong>0.08 ± 0.01</strong></td><td>0.42±0.03</td><td>0.51 ± 0.02</td><td>0.42 ± 0.02</td></tr><tr><td>Training time (s)</td><td>–</td><td>71.40±2.15</td><td>158.92 ± 5.32</td><td>8.01±2.02</td><td><strong>0.41 ± 0.01</strong></td><td>0.65 ± 0.01</td></tr></tbody></table>

**Table 2**

Comparison of RMSE values using the nested LRM database studied in (Arenzana et al., 2021) with $N_{LF} = 1949$ and $N_{HF} = 250$. Evaluation performed on 1000 test observations. SMFCK and SMoGP models use $M_{LF} = 50$ and $M_{HF} = 20$ inducing points obtained from Algorithm 1. The smallest RMSE values and training times for each output are highlighted in bold.

<table><thead><tr><th rowspan="2">Method</th><th colspan="3">RMSE</th><th colspan="3">Training time (s)</th></tr><tr><th>C<sub>x</sub></th><th>C<sub>z</sub></th><th>C<sub>M</sub></th><th>C<sub>x</sub></th><th>C<sub>z</sub></th><th>C<sub>M</sub></th></tr></thead><tbody><tr><td>HK</td><td>0.0552</td><td>0.6172</td><td>0.0504</td><td>–</td><td>–</td><td>–</td></tr><tr><td>MoGP</td><td>0.0555</td><td>0.6178</td><td>0.0515</td><td>82.8721</td><td>80.2545</td><td>79.7841</td></tr><tr><td>MFCK</td><td>0.0515</td><td>0.6172</td><td><strong>0.0504</strong></td><td>53.5781</td><td>57.0937</td><td>52.9531</td></tr><tr><td>MFCK (ours)</td><td>0.0555</td><td>0.6171</td><td>0.0530</td><td>109.4062</td><td>103.8593</td><td>102.4073</td></tr><tr><td>SMoGP</td><td>0.0555</td><td>0.6182</td><td>0.0519</td><td>42.1245</td><td>42.8742</td><td>41.7465</td></tr><tr><td>SMFCK-FITC (ours)</td><td>0.0513</td><td><strong>0.6058</strong></td><td><strong>0.0504</strong></td><td><strong>35.8221</strong></td><td><strong>31.4187</strong></td><td><strong>30.4583</strong></td></tr><tr><td>SMFCK-VFE (ours)</td><td><strong>0.0511</strong></td><td>0.6108</td><td><strong>0.0504</strong></td><td>39.2567</td><td>33.1232</td><td>33.6523</td></tr></tbody></table>

Additional experiments on the 2D Rosenbrock function with three fidelity levels are reported in Appendix C. The conclusions obtained there are consistent with those drawn from Table 1, showing improved performance and slightly increased training times as models are trained with progressively more data.

**Remark.** Theoretical results of (Quinonero-Candela and Rasmussen, 2005; Burt, Rasmussen and van der Wilk, 2020), showed that sparse models like FITC can recover the exact GP when the number of inducing points matches exactly the training points $M = N$. Likewise, Titsias (2009a) demonstrated that the VFE framework provides a lower bound on the marginal likelihood that becomes tight as the inducing set becomes dense in the input space. Theoretically, the sparse approximation recovers the exact GP in the limiting case. However, the idea behind sparsity is to introduce a set of inducing points $M \ll N$ to obtain similar predictive performance while reducing computational complexity.

## 5.2. Large Reference Model (LRM) applications

In this section, we present the results obtained by applying the sparse multi-fidelity methods described in Section 4 to a large-scale aerodynamic database based on the NASA Common Reference Model (CRM) (Vassberg, Dehaan, Rivers and Wahls, 2008).

The dataset used in this study consists of WT experiments conducted at ONERA's WT department, and numerical simulations generated using the elsA CFD solver (Cambier, Heib and Plot, 2013). Experimental testing for the WT experiments is performed in two different ONERA WT centers: the F1-WT, focusing on high Reynolds number values under low subsonic flow, and the S1-WT, which focuses on transonic speeds. The CRM shape is used to build two LRM: the first one for S1-WT with a wing span of 3.5m (scale 1/16.835) and the second one for F1-WT with a wing
<caption>**Table 3**</caption>

Ranges of the input parameters conducted at ONERA's WT department, and numerical simulations generated using the elsa CFD solver.

<table><thead><tr><th></th><th>WT</th><th>CFD</th></tr></thead><tbody><tr><td>Mach number <i>M</i><sub>a</sub></td><td>[0.01,0.95]</td><td>[0.01,0.97]</td></tr><tr><td>Reynolds number <i>Re</i> [10<sup>6</sup>]</td><td>[0.77,6.60]</td><td>[0.1,30.0]</td></tr><tr><td>Angle of attack <i>α</i> [deg]</td><td>[-9.3, 20.5]</td><td>[-10.0,30.0]</td></tr></tbody></table>

span of 3m (scale 1/19.5). The CFD data is acquired using Reynolds-averaged Navier-Stokes (RANS) equations. Here, the equation Spalart-Allmaras turbulence model is used. This CFD data uses a mesh named L1, composed of 63.9×10³ cells.

We must remark that CFD data have been taken in a parametric space larger than that of WT data. This is chosen to stretch out the study to the entire flight aircraft envelope. This is meant to illustrate the complementary dataset nature of the two datasets, one being more consistent with reality but covering only a limited domain of the parametric space, and the other one being used to expand that knowledge. We can see this characteristic in the reduced aerodynamic dataset used in (Arenzana et al., 2021) and later in Figure 5.

In this study, the multi-fidelity dataset consists of three input variables: the Mach number, *Ma*, defined as the ratio between the flow velocity and the local speed of sound, which characterizes compressibility effects in aerodynamics (Anderson, 2010); the Reynolds number, *Re*, a dimensionless quantity representing the ratio of inertial to viscous forces in the flow, thereby governing the transition between laminar and turbulent regimes (White, 2006); and the angle of attack, *α*, which denotes the inclination of the airfoil or body relative to the incoming flow direction and strongly influences lift and drag generation (Bertin and Smith, 2008). Table 3 shows the intervals of the input variables for each data source. In the implementation, the Reynolds number is transformed on a base-10 logarithmic scale to account for its wide range of values (see Figure 5).

The outputs considered are three aerodynamic coefficients: the drag coefficient, $C_x$, which quantifies the dimensionless resistance force acting opposite to the flow direction; the lift coefficient, $C_z$, which measures the dimensionless force component perpendicular to the flow; and the pitching moment coefficient, $C_M$, which represents the dimensionless moment about a reference point, typically the aerodynamic center, and characterizes the longitudinal stability of the configuration (Anderson, 2010). These input-output variables collectively define the aerodynamic response surface to be explored and modeled.

In our experiments, we analyze two configurations:

(i) The nested DoE setting studied in (Arenzana et al., 2021) with 250 WT observations and 1949 CFD simulations. The goal of this experiment is to enable comparison of the results with those of the baseline MFK model.

(ii) A larger non-nested DoE setting with 32595 WT observations and 1949 CFD simulations. This experiment allows us to assess the benefits of integrating sparse approximations into the MFCK model.

### 5.2.1. LRM experiment with a nested DoE

In this section, we consider the multi-fidelity dataset used in (Arenzana et al., 2021), which consists of a fully nested structure. Since this nested configuration is compatible with the recursive MFK model described in Section 2, we then propose comparing the performance of both MFCK and SMFCK using the results provided by MFK as a baseline. Model evaluation is performed using a WT validation set consisting of 1000 test samples constructed via LHS, with a new parametrization of the input space to guarantee hypercube domains, as explained in (Arenzana et al., 2021). In this experiment, the number of inducing points is determined by Algorithm 1 obtaining $M_{HF} = 20$ for the high-fidelity level and $M_{LF} = 50$ for the low-fidelity level.

Table 2 presents a comparison of the predictive performance and training time for HK, MoGP, MFK, MFCK, SMoGP, and SMFCK (both with FITC and VFE approximations). We observe that the predictive performance of the exact models (HK, MoGP, MFK, and MFCK) is comparable to the sparse approximations (SMoGP, SMFCK-FITC, and SMFCK-VFE). Among the exact models, MFK trains significantly faster than the other ones, thanks to its recursive formulation that avoids the block-wise structure used in MFCK. The nested design leads to smaller covariance matrices and simpler linear algebra operations, reducing both computational cost and memory usage.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Pairwise histograms of the aerodynamic input parameters: Mach number (Ma), Reynolds number (Re) on a base-10 logarithmic scale, and angle of attack (α). WT data are shown in orange and CFD simulations in blue. Diagonal panels display 1D marginal histograms, while off-diagonal panels show 2D histograms. Darker regions in the 2D plots indicate a higher density of samples. This non-nested dataset comprises 32595 WT experiments and 1949 CFD simulations.</figcaption></figure>

We can also note that the introduction of sparse approximations further reduces the computational complexity while preserving predictive accuracy, as evidenced by the RMSE values and the substantially lower training times achieved by SMFCK-FITC and SMFCK-VFE. In terms of computational efficiency, more general models such as SMoGP exhibit higher RMSE values and longer training times when compared to SMFCK-based approaches. Among them, SMFCK-FITC consistently achieves the lowest training times across all outputs, followed closely by SMFCK-VFE. Both sparse variants significantly outperform the corresponding exact models in terms of training time, with SMFCK-FITC being up to 40% faster than MFK.

These results demonstrate that sparse multi-fidelity approaches can maintain predictive performance while relying on a reduced number of inducing points. We should note that sparse methods can sometimes yield lower RMSE values than the exact GP-based models due to factors such as improved numerical stability, regularization effects, or favorable optimization dynamics, particularly in high-dimensional or ill-conditioned settings (Burt et al., 2020).

### 5.2.2. LRM experiment with a non-nested DoE

We now seek to leverage all available information from WT experiments. We recall that this database now contains 32595 WT observations, while preserving the same CFD data considered in the previous nested LRM experiment. Figure 5 shows the large multi-fidelity dataset, constructed following a non-nested DoE, where the WT experiments and CFD simulations were conducted independently without a nested DoE structure.

In this experiment, a fair comparison with the MFK model is not possible, as the nested DoE assumption is not satisfied. In addition, results for the HK, MoGP and MFCK models cannot be assessed due to the computational overhead of their training. No convergence has been observed after one day of computation. However, it is still possible to evaluate the benefit of combining both CFD and WT data sources when considering sparse implementations. To this end, we compare our approach to standard single-fidelity sparse Kriging (SK). We leverage two different SK models: one trained on WT data, referred to as HF-SK, and another trained on CFD simulations, referred to as LF-SK. The performance of both models is evaluated over WT test set. To assess the computational complexity of the sparse methods inducing points is set equally in all the methodologies, *M* = 250. In the context of the multi-fidelity approach, the points are distributed as follows: 200 points for the LF and 50 for the HF. All inducing point locations distributed using *k*-means.
**Table 4**

Comparison of RMSE using the non-nested multi-fidelity database with $N_{LF} = 1949$, $N_{HF} = 32595$. Test over 3260 observations taken from the WT experiments. Results are reported for two single-fidelity SK models (one trained on LF CFD simulations - LF-SK, and the other on HF WT data - HF-SK), as well as for SMoGP and SMFCK. The smallest RMSE values and training times for each output are highlighted in bold.

<table><thead><tr><th rowspan="2">Method</th><th rowspan="2">Approx.</th><th rowspan="2">Inducing points</th><th colspan="3">RMSE</th><th colspan="3">Training time (s)</th></tr><tr><th>C<sub>x</sub></th><th>C<sub>z</sub></th><th>C<sub>M</sub></th><th>C<sub>x</sub></th><th>C<sub>z</sub></th><th>C<sub>M</sub></th></tr></thead><tbody><tr><td rowspan="2">Single-fidelity LF-SK</td><td>FITC</td><td>M = 250</td><td>0.1514</td><td>0.4925</td><td>0.1782</td><td><b>9.2125</b></td><td><b>9.8715</b></td><td><b>9.7465</b></td></tr><tr><td>VFE</td><td>M = 250</td><td>0.1421</td><td>0.4918</td><td>0.1693</td><td>9.6036</td><td>10.2222</td><td>9.3622</td></tr><tr><td rowspan="2">Single-fidelity HF-SK</td><td>FITC</td><td>M = 250</td><td>0.0623</td><td>0.0215</td><td>0.0216</td><td>79.5911</td><td>79.9044</td><td>73.6507</td></tr><tr><td>VFE</td><td>M = 250</td><td>0.0611</td><td>0.0212</td><td>0.0215</td><td>80.5145</td><td>86.9760</td><td>87.1453</td></tr><tr><td>SMoGP</td><td>VFE</td><td>M<sub>LF</sub> = 200, M<sub>HF</sub> = 50</td><td>0.0650</td><td>0.0246</td><td>0.0255</td><td>226.3212</td><td>266.9339</td><td>142.3213</td></tr><tr><td rowspan="2">SMFCK (ours)</td><td>FITC</td><td>M<sub>LF</sub> = 200, M<sub>HF</sub> = 50</td><td>0.0511</td><td>0.0176</td><td>0.0108</td><td>122.1342</td><td>126.1241</td><td>121.114</td></tr><tr><td>VFE</td><td>M<sub>LF</sub> = 200, M<sub>HF</sub> = 50</td><td><b>0.0403</b></td><td><b>0.0157</b></td><td><b>0.0098</b></td><td>132.7643</td><td>138.8763</td><td>130.0121</td></tr></tbody></table>

Table 4 shows the RMSE performance and computational time for the LF-SK, HF-SK, SMoGP and SMFCK, with evaluation performed on 3260 WT test points excluded from the training dataset. These test points were randomly removed from the original dataset to serve as validation data. We observe that single-fidelity HF-SK models, trained with 250 inducing points, achieve moderate RMSE values across all aerodynamic coefficients. In contrast, as anticipated, the LF-SK models exhibit increased RMSE values, reflecting the lower reliability of the CFD data. On the other hand, we observe that models accounting for both CFD and WT test data overall outperform those relying on a single data source, highlighting the advantage of data fusion over the use of either source alone. Notably, SMFCK-VFE delivers the lowest errors: $C_x = 0.0403$, $C_z = 0.0157$, and $C_M = 0.0098$, representing improvements of approximately 51%, 31%, and 55%, respectively, over the single-fidelity HF-SK-VFE. Additionally, SMoGP presents higher RMSE values compared to SMFCK models. This difference reflects the advantage of the autoregressive assumption in SMFCK, which enhances model performance in scenarios where the fidelity functions can be effectively represented through an autoregressive formulation.

In terms of the training time, the computational cost of SK-based approaches is relatively low, reflecting the single-fidelity nature. Contrarily, since sparse methods incorporate both LF and HF data, they make more comprehensive use of the available information, leading to an increase in training time that remains reasonable. This difference in computational time can be better understood by examining their respective computational complexities. The standard SK model scales with a complexity of $O(NM^2)$, where N denotes the number of training points and M the number of inducing variables. In contrast, SMFCK, with its block-wise matrix formulation, exhibits a more involved complexity of $O(\sum_{\kappa=0}^{L} N_{\kappa} (\sum_{\kappa=0}^{L} M_{\kappa})^2)$, reflecting the contribution from each fidelity level $\kappa$. Although this formulation leads to increased computational cost, it maintains consistency with the structure of sparse approximations by preserving localized dependencies across fidelity levels. Empirically, even when the total number of inducing points is kept constant across both models, the SMFCK approach results in a computational overhead of approximately 30 seconds in the considered experiments. It is important to note, however, that the HF-SK model is trained on the full set of 32595 WT measurements alone, while SMFCK integrates this data with an additional 1949 CFD simulations. This multi-fidelity integration inherently modifies the structure and size of the matrices involved, influencing the observed runtime differences.

Finally, Figure 6 presents the standardized residuals at the validation points from the WT experiments for the outputs $C_x$, $C_z$, and $C_M$, which are computed as follows:

$$ \text{Standardized residual}(x_i) = \frac{y(x_i) - \hat{y}(x_i)}{s(x_i)}, $$

where $y(x_i)$ corresponds to the observed value at point $x_i$, $\hat{y}(x_i)$ to the Kriging prediction, and $s(x_i)$ to the estimated standard deviation. Residuals beyond ±3 standard deviations can be considered as outliers.

Figure 6 shows a clear reduction in the dispersion of the residuals, as well as a decrease in the number of outliers, when using the multi-fidelity approach. Specifically, the proportion of outliers decreases from 2.09%, 1.56%, and 0.98% in the single-fidelity HF-SK-VFE baseline to 0.28%, 0.77%, and 0.12% in SMFCK-VFE for the outputs $C_x$, $C_z$, and $C_M$, respectively, corresponding to an average reduction of 75%.
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Standardized residuals for $C_x$, $C_z$ and $C_M$ coefficients over 3260 test WT points. Results are shown for (top) the HF-SK-VFE model with $M = 250$, and (bottom) the SMFCK-VFE model with $M_{LF} = 200$, $M_{HF} = 50$. Residuals falling outside the [-3, 3] range are highlighted in red.</figcaption></figure>

# 6. Conclusions and perspectives

In this work, we revisit the MFK model originally introduced in (Kennedy and O'Hagan, 2000) and further explored in (Le Gratiet and Garnier, 2014), which has been widely adopted in the aerodynamics community and other related engineering fields. We propose a new formulation, referred to as MFCK, in which multiple fidelity levels are jointly modeled through a correlation structure that incorporates cross-covariances based on an autoregressive model. Reinforcing the framework with such a model enhances interpretability compared to purely data-driven approaches.

Two key limitations of the traditional MFK approach are addressed: (i) the requirement for nested DoEs, and (ii) the lack of scalability to large datasets. By enabling the use of non-nested DoEs, the proposed framework allows for independent data acquisition from experimental tests and numerical simulations, thus offering greater flexibility in real-world settings. Unlike other multi-fidelity frameworks capable of handling non-nested DoEs, such as HK (Han and Görtz, 2012), MFCK propagates uncertainty and also supports bidirectional knowledge transfer, meaning that information from the HF data can also improve predictions of the LF function. To improve scalability, we adapt the FITC and VFE sparse approximation techniques to the multi-fidelity context. Overcoming these limitations extends the applicability of the multi-fidelity modeling approach.

According to our experiments, we show that the proposed MFCK provides a similar prediction performance compared to MFK under nested DoE assumptions. When considering non-nested DoE, we observe that MFCK still leads to accurate predictions as for HK and MoGP, but unlike the recursive MFK which may propagate undesired uncertainty from the LF level to the HF level at locations where the HF function is observed (as illustrated in the experiment of Figure 1).

When using sparse approximations, MFCK-based models often require a small tractable number of inducing points to provide accurate predictions. While training independent single-fidelity SK models for each fidelity level results in comparable per-level computational complexity, this approach fails to exploit the correlations between fidelity levels. Consequently, although similar computational costs, SMFCK-based models provide a more efficient and informative alternative, particularly in large-scale settings where LF data can enhance HF predictions (as demonstrated in Table 4). In our experiments, we observed that SMFCK models have achieved improved accuracy while maintaining the computational efficiency of standard SK approaches.

In the presented results, we observed the behavior described in (Bauer et al., 2016), where the FITC approximation tends to underestimate the noise variance, while VFE tends to overestimate it. As shown in Figure 2, this difference is directly reflected in the resulting confidence intervals. Moreover, VFE provides a more principled variational bound
and exhibits increased robustness to outliers, as illustrated in Figure 5. From a practical perspective, our results suggest that FITC is well-suited for very fast, point-prediction-oriented tasks, whereas VFE should be preferred when reliable uncertainty quantification and robust performance are critical.

Convergence-rate analyses of sparse approximations (Burt et al., 2019) indicate that accurate and stable predictions can be achieved with several inducing points that grow only logarithmically with dataset size. This theoretical insight explains why, in practice, both FITC and VFE approximations achieve near-exact model performance using relatively few inducing points, while maintaining favorable computational scalability for large-scale problems. In our experiments, we show how selection criteria based on Algorithm 1 can provide reliable estimates of inducing points. Based on our experiments, including both synthetic and LRM-based cases, we observed that Algorithm 1 frequently converged to values of $M_k \in [[3D \log(N_k)], [5D \log(N_k)]]$. Consequently, thinner intervals with $c \in [3, 5]$ could be considered in practice, as they may offer a favorable trade-off between computational cost and predictive performance.

In this work, inducing point locations are estimated via $k$-means, providing a set of data-representative points. However, as shown in (Titsias, 2009a), their locations can be treated as variational parameters and optimized directly within the VFE framework. In that case, the inducing points will be incorporated into the marginal likelihood through the variational objective, removing the need for heuristic selection procedures. Then, a promising direction of our framework consists in constructing a fully variational multi-fidelity formulation. While this may result in increased computational cost, since the complexity of the optimization grows with the total number of inducing points $M = \sum_{k=0}^{L} M_k$, it could lead to improved predictive accuracy.

Although the MFCK and SMFCK frameworks are established for any input dimension $D \ge 1$, as they are based on Kriging, they suffer from the curse of dimensionality, whereby the amount of data required for accurate learning grows exponentially with the input dimension, leading to poor generalization and increased computational burden (Rasmussen and Williams, 2005). While sparse approximations mitigate the computational complexity with a set of inducing points $M \ll N$, they do not fully eliminate the curse of dimensionality, since an increasingly large number of inducing points may still be required to adequately cover high-dimensional domains. Future extensions of this work will therefore focus on integrating dimensionality-reduction techniques such as Partial Least Squares (Wold, 1975) to improve scalability and robustness in high-dimensional multi-fidelity settings. Regarding the complexity of parameter optimization for the MFCK model in high-dimensional settings or when considering a large number of fidelity levels, a cooperative component approach may also be adopted (Zhan, Wu, Xing and Li, 2024; Pretsch, Arsenyev, Bartoli and Duddeck, 2025), in which the full parameter set is divided into several subsets. Each subset is then optimized while keeping the other subsets fixed. Another potential future work is the application of our framework to multi-fidelity adaptive learning scenarios, such as global optimization (Meliani et al., 2019; Charayron, 2023) and adaptive DoE (Jakeman, Friedman, Eldred, Tamellini, Gorodetsky and Allaire, 2022) involving large-scale datasets.

**Acknowledgement.** This research was carried out as part of the ECOS Nord project C25M01. Part of the work was conducted while AFLL was affiliated with CERAMATHS, Univ. Polytechnique Hauts-de-France. We thank Rémi Lafage (ONERA) for his support in integrating the Python implementation of the MFCK and SMFCK models into the toolbox SMT. We also thank Sylvain Mouton (ONERA) for the insightful discussions and valuable advice regarding the LRM database.
# A. Mathematical developments related to the MFCK based on the autoregressive model

## A.1. Autoregressive model

Consider the autoregressive model given in Eq. (3), and let $f_0(\cdot) := \gamma_0(\cdot)$ by convention. Then, the $\kappa$-th fidelity function $f_\kappa(\cdot)$ can be expressed recursively in terms of the latent functions $\gamma_0(\cdot), \dots, \gamma_\kappa(\cdot)$ as follows:

$$ 
\begin{aligned}
f_{\kappa}(\mathbf{x}) &= \rho_{\kappa} f_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}) \\
&= \rho_{\kappa} \rho_{\kappa-1} f_{\kappa-2}(\mathbf{x}) + \rho_{\kappa} \gamma_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}) \\
&= \rho_{\kappa} \rho_{\kappa-1} f_{\kappa-2}(\mathbf{x}) + \rho_{\kappa} \gamma_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}) \\
&\vdots \\
&= (\rho_{\kappa} \cdots \rho_1) \gamma_0(\mathbf{x}) + (\rho_{\kappa} \cdots \rho_2) \gamma_1(\mathbf{x}) + \cdots + \rho_{\kappa} \rho_{\kappa-1} f_{\kappa-2}(\mathbf{x}) + \rho_{\kappa} \gamma_{\kappa-1}(\mathbf{x}) + \gamma_{\kappa}(\mathbf{x}) \\
&= \sum_{j=0}^{\kappa-1} \left( \left( \prod_{i=j+1}^{\kappa} \rho_i \right) \gamma_j(\mathbf{x}) \right) + \gamma_{\kappa}(\mathbf{x}) \\
&= \sum_{j=0}^{\kappa} \eta_{j,\kappa} \gamma_j(\mathbf{x}),
\end{aligned}
$$

with $\eta_{j,\kappa} = \begin{cases} \prod_{i=j+1}^{\kappa} \rho_i, & \text{if } j < \kappa, \\ 1, & \text{if } j = \kappa. \end{cases}$

## A.2. Computation of the cross-covariance function for the MFCK

Consider the autoregressive scheme given in Eq. (3). As shown in Appendix A.1, the latent function $f_\kappa(\cdot)$ can be expressed as $f_\kappa(\cdot) = \sum_{j=0}^\kappa \eta_{j,\kappa} \gamma_j(\cdot)$, where the coefficients $\eta_{j,\kappa}$ are defined in Eq. (5). Recall from Eq. (7) that the observation model is given by $y_\kappa(\cdot) = f_\kappa(\cdot) + \epsilon_\kappa$ with $\epsilon_\kappa \sim \mathcal{N}(0, \tau_\kappa^2)$. Assuming that the noise terms $\epsilon_0, \dots, \epsilon_L$ are i.i.d. and independent of the latent GPs $\gamma_0(\cdot), \dots, \gamma_L(\cdot)$, the covariance between $y_\kappa(\mathbf{x})$ and $y_{\kappa'}(\mathbf{x}')$ is then given by:

$$ 
\begin{aligned}
c_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') &= \underbrace{\operatorname{cov}(f_{\kappa}(\mathbf{x}), f_{\kappa'}(\mathbf{x}'))}_{k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}')} + \operatorname{cov}(\epsilon_{\kappa}, \epsilon_{\kappa'}) = \operatorname{cov}\left(\sum_{j=0}^{\kappa} \eta_{j, \kappa} \gamma_j(\mathbf{x}), \sum_{j=0}^{\kappa'} \eta_{j, \kappa'} \gamma_j(\mathbf{x}')\right) + \tau_{\kappa}^2 \delta_{\kappa, \kappa'} \delta_{\mathbf{x}, \mathbf{x}'} \\
&= \sum_{j=0}^{\kappa} \sum_{j'=0}^{\kappa'} \eta_{j, \kappa} \eta_{j', \kappa'} \operatorname{cov}(\gamma_j(\mathbf{x}), \gamma_{j'}(\mathbf{x}')) + \tau_{\kappa}^2 \delta_{\kappa, \kappa'} \delta_{\mathbf{x}, \mathbf{x}'},
\end{aligned}
$$

where $\mathbf{x}, \mathbf{x}' \in \mathcal{D}$, and $\delta_{a,b}$ is the Kronecker delta function that is equal to 1 if $a=b$, and 0 otherwise.

We now turn to the computation of the term $k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}')$ for all $\kappa, \kappa' \in \{0, \dots, L\}$. To this end, we analyze separately the cases where $\kappa = \kappa'$ and $\kappa > \kappa'$. The case where $\kappa < \kappa'$ follows immediately by symmetry, since $\operatorname{cov}(f_\kappa(\mathbf{x}), f_{\kappa'}(\mathbf{x}')) = \operatorname{cov}(f_{\kappa'}(\mathbf{x}'), f_\kappa(\mathbf{x}))$. In that case, it suffices simply to interchange the roles of $\kappa$ and $\kappa'$ for the case when $\kappa > \kappa'$. In both cases, we need to make use of the fact that the GP priors $\gamma_\kappa(\cdot)$ are independent, i.e., $\operatorname{cov}(\gamma_\kappa(\mathbf{x}), \gamma_{\kappa'}(\mathbf{x}')) = 0$ for all $\kappa \neq \kappa'$.

* Case when $\kappa = \kappa'$.

$$ k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') = \sum_{j=0}^{\kappa} \eta_{j, \kappa} \eta_{j, \kappa'} \operatorname{cov}(\gamma_j(\mathbf{x}), \gamma_j(\mathbf{x}')) = \sum_{j=0}^{\kappa} \eta_{j, \kappa} \eta_{j, \kappa'} k_{\gamma_j}(\mathbf{x}, \mathbf{x}'). $$

* Case when $\kappa > \kappa'$.

$$ k_{\kappa, \kappa'}(\mathbf{x}, \mathbf{x}') = \sum_{j=0}^{\kappa'} \sum_{j'=0}^{\kappa'} \eta_{j, \kappa} \eta_{j', \kappa'} \operatorname{cov}(\gamma_j(\mathbf{x}), \gamma_{j'}(\mathbf{x}')) = \sum_{j=0}^{\kappa'} \eta_{j, \kappa} \eta_{j, \kappa'} k_{\gamma_j}(\mathbf{x}, \mathbf{x}'). $$

We note that in both cases, the cross-covariance terms are equal to zero as $\gamma_\kappa(\cdot) \perp \gamma_{\kappa'}(\cdot) \forall \kappa \neq \kappa'$. The previous results can be unified into a single expression by summing up to $\min(\kappa, \kappa')$.
## B. Derivation of the VFE likelihood for the MFCK

As in Section 4, denote $\mathbf{y} = [y_0, \dots, y_L]^\top$ and $\mathbf{u} = [u_0, \dots, u_L]^\top$ where $y_k = y_k(\mathbf{X}_k)$ and $u_k = y_k(\mathbf{Z}_k)$ for all $k \in \{0, \dots, L\}$. The log-marginal likelihood $p(\mathbf{y})$ can be written as follows:

$$ 
\begin{aligned}
\log p(\mathbf{y}) &= \log \iint p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\mathbf{u})p(\mathbf{u})d\mathbf{f}d\mathbf{u} \\
&= \log \iint p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\mathbf{u})\frac{q(\mathbf{u})}{q(\mathbf{u})}p(\mathbf{u})d\mathbf{f}d\mathbf{u} \quad (\text{recall that } q(\mathbf{f}, \mathbf{u}) = p(\mathbf{f}|\mathbf{u})q(\mathbf{u})) \\
&= \log \iint p(\mathbf{y}|\mathbf{f})\frac{p(\mathbf{u})}{q(\mathbf{u})}q(\mathbf{f}, \mathbf{u})d\mathbf{f}d\mathbf{u} \\
&= \log \mathbb{E}_{q(\mathbf{f}, \mathbf{u})} \left( p(\mathbf{y}|\mathbf{f})\frac{p(\mathbf{u})}{q(\mathbf{u})} \right).
\end{aligned} 
\quad (28) $$

Using Jensen's inequality (Jensen, 1906), we can express the log of the expectation as a lower bound, the formulation of Jensen implies that $\varphi(\mathbb{E}(\mathbf{X})) \le \mathbb{E}(\varphi(\mathbf{X}))$. Applying this property, we have:<fnref n="2" />

$$ 
\begin{aligned}
\log p(\mathbf{y}) &\le \mathbb{E}_{q(\mathbf{f}, \mathbf{u})} \log \left( p(\mathbf{y}|\mathbf{f}) \frac{p(\mathbf{u})}{q(\mathbf{u})} \right) \\
&= \iint \left[ \log(p(\mathbf{y}|\mathbf{f}) + \log \frac{p(\mathbf{u})}{q(\mathbf{u})}) \right] q(\mathbf{f}, \mathbf{u}) d\mathbf{f}d\mathbf{u} \\
&= \iint \left[ q(\mathbf{u}) \int \log p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\mathbf{u})d\mathbf{f} + q(\mathbf{u}) \int \log \frac{p(\mathbf{u})}{q(\mathbf{u})}p(\mathbf{f}|\mathbf{u})d\mathbf{f} \right] d\mathbf{u} \\
&= \iint q(\mathbf{u}) \left[ \int \log p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\mathbf{u})d\mathbf{f} + \log \frac{p(\mathbf{u})}{q(\mathbf{u})} \right] d\mathbf{u}.
\end{aligned} 
$$

We first focus on the computation of the integral $G(\mathbf{f}, \mathbf{u}) := \int \log p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|\mathbf{u})d\mathbf{f}$. According to Eq. (3), we have that $p(\mathbf{y}|\mathbf{f}) = \mathcal{N}(\mathbf{f}, \mathbf{D})$ where $\mathbf{D} = \text{diag}(\tau_0^2 I_0, \dots, \tau_k^2 I_k)$. Then, we have:

$$ 
\begin{aligned}
G(\mathbf{u}, \mathbf{y}) &= \int \log p(\mathbf{y}|\mathbf{f}) p(\mathbf{f}|\mathbf{u}) d\mathbf{f} \\
&= \int \left[ -\frac{1}{2} \log((2\pi)^N |\mathbf{D}|) - \frac{1}{2} \mathbf{y}^\top \mathbf{D}^{-1} \mathbf{y} + \mathbf{y}^\top \mathbf{D}^{-1} \mathbf{f} - \frac{1}{2} \mathbf{f}^\top \mathbf{D}^{-1} \mathbf{f} \right] p(\mathbf{f}|\mathbf{u}) d\mathbf{f} \\
&= -\frac{1}{2} \log((2\pi)^N |\mathbf{D}|) - \frac{1}{2} \mathbf{y}^\top \mathbf{D}^{-1} \mathbf{y} + \mathbf{y}^\top \mathbf{D}^{-1} \int \mathbf{f} p(\mathbf{f}|\mathbf{u}) d\mathbf{f} - \frac{1}{2} \int \mathbf{f}^\top \mathbf{D}^{-1} \mathbf{f} p(\mathbf{f}|\mathbf{u}) d\mathbf{f}.
\end{aligned} 
$$

Considering the factorization $\mathbf{D} = \mathbf{D}^{\frac{1}{2}} \mathbf{D}^{\frac{1}{2}}$ with $\mathbf{D}^{\frac{1}{2}} = \text{diag}(\tau_0 I_0, \dots, \tau_k I_k)$, we obtain for the last term

$$ \int \mathbf{f}^\top \mathbf{D}^{-1} \mathbf{f} p(\mathbf{f}|\mathbf{u}) d\mathbf{f} = \int \text{tr} \left( \mathbf{D}^{-\frac{1}{2}} \mathbf{f} \mathbf{f}^\top \mathbf{D}^{-\frac{1}{2}} \right) p(\mathbf{f}|\mathbf{u}) d\mathbf{f} = \text{tr} \left( \mathbf{D}^{-1} \int \mathbf{f} \mathbf{f}^\top p(\mathbf{f}|\mathbf{u}) d\mathbf{f} \right). \quad (29) $$

Due to the conditioning properties of the Gaussian distribution, we know that $p(\mathbf{f}|\mathbf{u}) = \mathcal{N}(\alpha, \Gamma)$ with

$$ \alpha = K_{f,u} K_{u,u}^{-1} \mathbf{u}, \quad \Gamma = K_{f,f} - K_{f,u} K_{u,u}^{-1} K_{u,f}. $$

By definition, we have

$$ \int \mathbf{f} p(\mathbf{f}|\mathbf{u}) d\mathbf{f} := \mathbb{E}(\mathbf{f}|\mathbf{u}) = \alpha, \quad (30) $$

$$ \int \mathbf{f} \mathbf{f}^\top p(\mathbf{f}|\mathbf{u}) d\mathbf{f} := \mathbb{E}(\mathbf{f} \mathbf{f}^\top | \mathbf{u}) = \mathbb{E}(\mathbf{f}|\mathbf{u}) \mathbb{E}^\top(\mathbf{f}|\mathbf{u}) + \text{cov}(\mathbf{f}|\mathbf{u}) = \alpha \alpha^\top + \Gamma. $$

<footnote n="2">²This assumption is true since $\log_b(\cdot)$ is an increasing function for $b > 1$.</footnote>
Replacing the result in Eq. (29) and Eq. (30) in the computation of $G(u, y)$ we then have

$$ 
\begin{aligned}
G(u, y) &= -\frac{1}{2} \log((2\pi)^N |\mathbf{D}|) - \frac{1}{2} y^\top \mathbf{D}^{-1} y + y^\top \mathbf{D}^{-1} \boldsymbol{\alpha} - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} [\boldsymbol{\alpha}\boldsymbol{\alpha}^\top + \mathbf{\Gamma}]) \\
&= -\frac{1}{2} \log((2\pi)^N |\mathbf{D}|) - \frac{1}{2} y^\top \mathbf{D}^{-1} y + y^\top \mathbf{D}^{-1} \boldsymbol{\alpha} - \frac{1}{2} \boldsymbol{\alpha}^\top \mathbf{D}^{-1} \boldsymbol{\alpha} - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \\
&= \log \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}).
\end{aligned}
$$

Replacing the expression of $G(u, y)$ in the bound of the log-marginal likelihood, we find:

$$ 
\begin{aligned}
\log p(y) &\le \int q(u) \left[ \log \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) + \log \frac{p(u)}{q(u)} \right] du \quad (31) \\
&= \int q(u) \log \left( \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) \frac{p(u)}{q(u)} \right) du - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \\
&= \mathbb{E}_{q(u)} \left( \log \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) \frac{p(u)}{q(u)} \right) - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) =: \mathcal{L}(\theta, q).
\end{aligned}
$$

Note that we explicitly wrote the dependency of the bound $\mathcal{L}(\theta, q)$ with the variational distribution $q(\cdot)$ as the bound needs to be optimized w.r.t. to these parameters. This can be achieved by using the properties of Jensen's inequality. By inverting Jensen's inequality, taking the log outside the integral, we compute the marginal log likelihood $\mathcal{L}(\theta)$ given by:

$$ 
\begin{aligned}
\mathcal{L}(\theta) &= \log \mathbb{E}_{q(u)} \left( \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) \frac{p(u)}{q(u)} \right) - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \quad (32) \\
&= \log \int \mathcal{N}(y|\boldsymbol{\alpha}, \mathbf{D}) p(u) du - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \\
&= \log \left[ \mathcal{N} \left( y | \mathbf{0}, \mathbf{D} + \mathbf{K}_{f,u} \mathbf{K}_{u,u}^{-1} \mathbf{K}_{u,f} \right) \right] - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \\
&= \log \left[ \mathcal{N} \left( y | \mathbf{0}, \mathbf{D} + \mathcal{Q} \right) \right] - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}).
\end{aligned}
$$

Equation (32) defines the cost function used to estimate the model parameters $\theta$ via maximum likelihood estimation. Once this cost function is established for determining the optimal parameters under a given set of observations $(\mathbf{X}_\kappa)_{0\le\kappa\le L}$, the next step is to compute the optimal distribution $q^*(u)$ to obtain the conditional mean and variance. To do so, we use the variational objective $\mathcal{L}(\theta, q)$ given in Eq. (31):

$$ 
\begin{aligned}
\mathcal{L}(\theta, q) &= \int q(u) \log \left( \mathcal{N}(y | \boldsymbol{\alpha}, \mathbf{D}) \cdot \frac{p(u)}{q(u)} \right) du - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}) \\
&= \int q(u) \left[ \log \mathcal{N}(y | \boldsymbol{\alpha}, \mathbf{D}) + \log p(u) - \log q(u) \right] du - \frac{1}{2} \operatorname{tr}(\mathbf{D}^{-1} \mathbf{\Gamma}).
\end{aligned}
$$

To compute the optimal distribution $q^*(u)$ we apply the first optimality condition to $\mathcal{L}(\theta, q)$ with respect to $q(u)$:

$$ \frac{\partial \mathcal{L}(\theta, q)}{\partial q(u)} = \log \mathcal{N}(y | \boldsymbol{\alpha}, \mathbf{D}) + \log p(u) - \log q(u) = 0 \quad (33) $$

Where we obtain $\log q(u) = \log \mathcal{N}(y | \boldsymbol{\alpha}, \mathbf{D}) + \log p(u)$, and therefore, we can deduce that $q^*(u) \propto \mathcal{N}(y|\boldsymbol{\alpha})$. This optimal expression corresponds exactly to the true posterior distribution of the latent variables $u$ conditioned on the observed data $y$, $q^*(u) = p(u|y)$.

By applying Bayes' theorem, we can approximate $p(u|y)$ as we know that the marginal likelihood $p(y)$ acts as a normalization of the posterior:

$$ p(u|y) \propto p(y|u)p(u) = \mathcal{N}(y | \mathbf{K}_{y,u} \mathbf{K}_{u,u}^{-1} u, \mathbf{D}) \mathcal{N}(u | \mathbf{0}, \mathbf{K}_{u,u}). $$

This forms the basis for deriving the conditional mean $m_\kappa(\cdot)$ and conditional variance $v_\kappa(\cdot)$ presented in Eq. (26). For further details on the subsequent steps, we refer the reader to (Titsias, 2009a), as the resulting form of the posterior $p(u|y)$ is analogous.
**Table 5**

Comparison of RMSE (mean ± standard deviation for 10 replicates) and training time for SMoGP, and SMFCK (FITC and VFE) on the 2D Rosenbrock function using three different DoE configurations. For SMFCK and SMoGP, we fix $M_{HF} = 20$, $M_{MF} = 30$ and $M_{LF} = 40$ using the empirical rule described in Section 4.4. The inducing points selection is conducted via k-means. The smallest RMSE values and training times for each input configuration are highlighted in bold.

<table><thead><tr><th>N<sub>HF</sub></th><th>N<sub>MF</sub></th><th>N<sub>LF</sub></th><th>Metric</th><th>SMoGP</th><th>SMFCK-FITC</th><th>SMFCK-VFE</th></tr></thead><tbody><tr><td rowspan="2">32</td><td rowspan="2">64</td><td rowspan="2">128</td><td>RMSE (HF)</td><td>8.221±3.214</td><td>3.248± 0.878</td><td><strong>2.212 ± 0.0645</strong></td></tr><tr><td>Training time (s)</td><td>7.801±1.824</td><td><strong>0.735 ± 0.254</strong></td><td>0.878 ± 0.057</td></tr><tr><td rowspan="2">128</td><td rowspan="2">256</td><td rowspan="2">512</td><td>RMSE (HF)</td><td>0.754±0.089</td><td>0.584 ± 0.245</td><td><strong>0.512 ± 0.179</strong></td></tr><tr><td>Training time (s)</td><td>8.741±2.743</td><td><strong>2.075 ± 1.131</strong></td><td>2.850 ± 1.742</td></tr><tr><td rowspan="2">512</td><td rowspan="2">1024</td><td rowspan="2">2048</td><td>RMSE (HF)</td><td>0.425±0.035</td><td>0.507 ± 0.022</td><td><strong>0.416 ± 0.021</strong></td></tr><tr><td>Training time (s)</td><td>11.124±2.224</td><td><strong>2.566 ± 1.021</strong></td><td>2.798 ± 1.254</td></tr></tbody></table>

### C. 2D Rosenbrock example for 3 fidelity levels

Building on the experiment presented in Section 5.1, we consider a Rosenbrock function with three fidelity levels with the functions described as:

$$ 
\begin{aligned}
f_{HF}(\mathbf{x}) &= \sum_{i=1}^{D-1} 100(x_{i+1} - x_i^2)^2 + (1-x_i)^2, \\
f_{MF}(\mathbf{x}) &= \sum_{i=1}^{D-1} 50(x_{i+1} - x_i^2)^2 + (2+x_i)^2 - \sum_{i=1}^{D} 0.5x_i, \\
f_{LF}(\mathbf{x}) &= \frac{f_{HF}(\mathbf{x}) - (4 + 0.5 \sum_{i=1}^{D} x_i)}{10 + 0.25 \sum_{i=1}^{D} x_i}
\end{aligned}
$$

We generate data from these functions considering $D = 2$, and noise variances $\tau_{HF}^2 = 0.1$, $\tau_{MF}^2 = 0.2$, and $\tau_{LF}^2 = 0.4$.

Table 5 reports the results of an experiment on the 2D Rosenbrock function using three fidelity levels and three DoE configurations. For each configuration, we compare the predictive performance at the high-fidelity level in terms of RMSE, as well as the associated training time, for SMoGP and the SMFCK model using both FITC and VFE approximations. Across all configurations, the sparse SMFCK variants consistently achieve lower RMSEs than SMoGP, with the VFE formulation generally yielding the most accurate predictions, particularly in the low-data regime. In terms of computational cost, both SMFCK variants exhibit substantially reduced training times compared to SMoGP, with FITC being systematically the fastest.

Compared to Table 1, the results in Table 5 indicate a clear improvement in predictive performance when an additional intermediate fidelity level is incorporated. This highlights the benefit of leveraging richer hierarchical information to better inform the high-fidelity predictions. Nevertheless, the inclusion of more than three fidelity levels may lead to mixed outcomes, primarily due to the increased complexity of the resulting model and the associated optimization landscape, which becomes higher-dimensional and potentially more ill-conditioned. In practice, this added complexity can hinder stable parameter estimation and may offset the potential gains from additional fidelity information.

Moreover, to the best of our knowledge, real-world applications involving four or more fidelity levels remain scarce in the literature. For these reasons, we have limited our study of model stability and performance to a maximum of three fidelity levels. Finally, as the number of training samples increases at all fidelity levels, all considered methods exhibit systematic improvements in predictive accuracy, confirming the expected consistency of the models with respect to data enrichment.
## References

Anderson, J.D., 2010. Fundamentals of Aerodynamics. 5th ed., McGraw-Hill Education, New York.

Arenzana, R.C., López-Lopera, A.F., Mouton, S., Bartoli, N., Lefebvre, T., 2021. Multi-fidelity Gaussian process model for CFD and wind tunnel data fusion. *Aerobest 2021*.

Bauer, M., van der Wilk, M., Rasmussen, C.E., 2016. Understanding probabilistic sparse Gaussian process approximations, in: Lee, D., Sugiyama, M., von Luxburg, U., Guyon, I., Garnett, R. (Eds.), *Advances in Neural Information Processing Systems 29*, Curran Associates, Inc., pp. 1533–1541.

Bertin, J.J., Smith, R.M., 2008. Aerodynamics for Engineers. 5th ed., Pearson Prentice Hall, Upper Saddle River, NJ.

Bonilla, E.V., Chai, K.M.A., Williams, C.K.I., 2007. Multi-task Gaussian process prediction, in: International Conference on Neural Information Processing Systems, pp. 153–160.

Bouhlel, M.A., Hwang, J.T., Bartoli, N., Lafage, R., Morlier, J., Martins, J.R.R.A., 2019. A python surrogate modeling framework with derivatives. *Advances in Engineering Software*, 102662.

Brevault, L., Balesdent, M., Hebbal, A., 2020. Overview of Gaussian process based multi-fidelity techniques with variable relationship between fidelities, application to aerospace systems. *Aerospace Science and Technology* 107, 106339.

Burt, D.R., Rasmussen, C.E., van der Wilk, M., 2019. Rates of convergence for sparse variational gaussian process regression, in: *Proceedings of the 36th International Conference on Machine Learning*, PMLR. pp. 862–871.

Burt, D.R., Rasmussen, C.E., van der Wilk, M., 2020. Convergence of sparse variational inference in Gaussian processes regression, in: *Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS)*, PMLR. pp. 2700–2710.

Cambier, L., Heib, S., Plot, S., 2013. The onera elsA CFD software: Input from research and feedback from industry. *Mechanics & Industry* 14, 159–174.

Charayron, R., 2023. Advances in multi-fidelity Bayesian optimization for multi-disciplinary drone design. Ph.D. thesis, ISAE.

Christianson, R.B., Pollyea, R.M., Gramacy, R.B., 2023. Traditional Kriging versus modern Gaussian processes for large-scale mining data. *Statistical Analysis and Data Mining: The ASA Data Science Journal* 16, 488–506.

Fernández-Godino, M.G., 2023. Review of multi-fidelity models. *Advances in Computational Science and Engineering* 1, 351–400.

Forrester, A.I.J., Keane, A.J., O'Hagan, E.F.H., 2007. Multi-fidelity modeling for engineering design: A surrogate modeling approach, in: *Proceedings of the 2007 International Conference on Engineering Design*. springer, pp. 1–6.

Han, Z.H., Görtz, S., 2012. Hierarchical kriging model for variable-fidelity surrogate modeling. *AIAA journal* 50, 1885–1896.

Hensman, J., Fusi, N., Lawrence, N.D., 2013. Gaussian processes for big data. *Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence*, 282–290.

Jakeman, J.D., Friedman, S., Eldred, M.S., Tamellini, L., Gorodetsky, A.A., Allaire, D., 2022. Adaptive experimental design for multi-fidelity surrogate modeling of multi-disciplinary systems. *International Journal for Numerical Methods in Engineering* 123, 2760–2790.

Jensen, J.L.W.V., 1906. On convex functions and inequalities between mean values. *Acta Mathematica* 30, 175–193.

Johnson, F., Tinoco, E., Yu, N., 2005. 30 years of development and application of CFD at Boeing commercial airplanes, seattle. *Computers & Fluids* 34, 1115–1151.

Keane, A.J., 2012. Cokriging for robust design optimization. *AIAA Journal* 50, 2351–2364.

Kennedy, M., O'Hagan, A., 2000. Predicting the output from a complex computer code when fast approximations are available. *Biometrika* 87, 1–13.

Krige, D., 1951. A statistical approach to some basic mine valuation problems on the witwatersrand. *Journal of the Chemical, Metallurgical and Mining Society of South Africa* 52, 119–139.

Laslett, G.M., 1994. Kriging and splines: An empirical comparison of their predictive performance in some applications. *Journal of the American Statistical Association* 89, 391–400.

Le Gratiet, L., 2013. Multi-fidelity Gaussian process regression for computer experiments. Ph.D. thesis, Université Paris-Diderot-Paris VII.

Le Gratiet, L., Garnier, J., 2014. Recursive co-Kriging model for design of computer experiments with multiple levels of fidelity. *International Journal for Uncertainty Quantification* 4, 365–386.

López-Lopera, A.F., Idier, D., Rohmer, J., Bachoc, F., 2022. Multioutput gaussian processes with functional data: A study on coastal flood hazard assessment. *Reliability Engineering & System Safety* 218, 108139.

Lu, N., Li, Y., Mi, J., Huang, H., 2024. Amfgp: An active learning reliability analysis method based on multi-fidelity gaussian process surrogate model. *Reliability Engineering & System Safety* 246, 110020.

Mainini, L., Serani, A., Rumpfkeil, M.P., Minisci, E., Quagliarella, D., Pehlivan, H., Yildiz, S., Ficini, S., Pellegrini, R., Fiore, F.D., Bryson, D., Nikbay, M., Diez, M., Beran, P., 2022. Analytical benchmark problems for multifidelity optimization methods.

Malik, M., Bushnell, D., 2012. Role of computational fluid dynamics and wind tunnels in aeronautics. *Technical Report TP-2012-217602*. NASA.

Marrel, A., Iooss, B., 2024. Probabilistic surrogate modeling by gaussian process: A review on recent insights in estimation and validation. *Reliability Engineering & System Safety* 247, 110094.

Matheron, G., 1963. Principles of geostatistics. *Economic Geology* 58, 1246–1266.

Meliani, M., Bartoli, N., Lefebvre, T., Bouhlel, M.A., Martins, J., Morlier, J., 2019. Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design, in: *AIAA Aviation 2019 Forum*, Dallas, United States. pp. 1–18.

Miele, S., Karve, P., Mahadevan, S., 2023. Multi-fidelity physics-informed machine learning for probabilistic damage diagnosis. *Reliability Engineering & System Safety* 235, 109243.

O'Hagan, A., 2006. Bayesian analysis of computer code outputs: A tutorial. *Reliability Engineering & System Safety* 91, 1290–1300.

Perdikaris, P., Venturi, D., Royset, J.O., Karniadakis, G.E., 2015. Multi-fidelity modelling via recursive co-Kriging and Gaussian-Markov random fields. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 471, 20150018.

Pretsch, L., Arsenyev, I., Bartoli, N., Duddeck, F., 2025. Bayesian optimization of cooperative components for multi-stage aero-structural compressor blade design. *Structural and Multidisciplinary Optimization* 68, 84.
Quinonero-Candela, J., Rasmussen, C., 2005. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research 6, 1939–1959.

Rasmussen, C.E., Williams, C.K.I., 2005. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press.

Rumpfkeil, M.P., Bryson, D., Beran, P., 2022. Multi-fidelity sparse Polynomial Chaos and Kriging surrogate models applied to analytical benchmark problems. Algorithms 15.

Saves, P., Lafage, R., Bartoli, N., Diouane, Y., Bussemaker, J., Lefebvre, T., Hwang, J.T., Morlier, J., Martins, J.R.R.A., 2024. SMT 2.0: A surrogate modeling toolbox with a focus on hierarchical and mixed variables gaussian processes. Advances in Engineering Software 188, 103571.

Shang, X., Su, L., Fang, H., Zeng, B., Zhang, Z., 2023. An efficient multi-fidelity kriging surrogate model-based method for global sensitivity analysis. Reliability Engineering & System Safety 229, 108858.

Slotnick, J., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W., Lurie, E., Mavripilis, D., 2014. CFD vision 2030 study: A path to revolutionary computational aerosciences. Contract Report CR-2014-218178. NASA.

Snelson, E., Ghahramani, Z., 2005. Sparse Gaussian processes using pseudo-inputs. Advances in Neural Information Processing Systems 18, 1257–1264.

Stein, M.L., 1999. Interpolation of Spatial Data: Some Theory for Kriging. Springer-Verlag, New York.

Stoliker, F., 2005. Introduction to flight test engineering. AGARDograph AG-300, RTO/NATO.

Titsias, M., 2009a. Variational learning of inducing variables in sparse gaussian processes. Journal of Machine Learning Research - Proceedings Track 5, 567–574.

Titsias, M.K., 2009b. Variational model selection for sparse Gaussian process regression, in: Proceedings of the 26th International Conference on Machine Learning (ICML), ACM, Montreal, Canada. pp. 871–878.

Toal, D.J.J., 2023. Applications of multi-fidelity multi-output Kriging to engineering design optimization. Structural and Multidisciplinary Optimization.

Valayer, H., Bartoli, N., Castaño-Aguirre, M., Lafage, R., Lefebvre, T., López-Lopera, A.F., Mouton, S., 2024. A python toolbox for data-driven aerodynamic modeling using sparse gaussian processes. Aerospace 11.

Van Der Wilk, M., 2019. Sparse Gaussian Process Approximations and Applications. Ph.D. thesis. University of Cambridge.

Vassberg, J., Dehaan, M., Rivers, M., Wahls, R., 2008. Development of a common research model for applied CFD validation studies, in: American Institute of Aeronautics and Astronautics, pp. 1–22.

Wang, D., Qiu, H., Gao, L., Jiang, C., 2023. A subdomain uncertainty-guided kriging method with optimized feasibility metric for time-dependent reliability analysis. Reliability Engineering & System Safety 243, 109839.

White, F.M., 2006. Viscous Fluid Flow. 3rd ed., McGraw-Hill Education, New York.

Wold, H., 1975. Soft modeling by latent variables: The nonlinear iterative partial least squares approach, in: Perspectives in Probability and Statistics: Papers in Honour of M. S. Bartlett. Academic Press, pp. 520–540.

Yang, S., Kang, Y.E., Yee, K., 2023. Multi-fidelity modeling via regression-based hierarchical Kriging, in: The Proceedings of the 2021 Asia-Pacific International Symposium on Aerospace Technology (APISAT 2021), Volume 1, Springer Nature Singapore, Singapore. pp. 643–652.

Yıldız, S., Pehlivan-Solak, H., Diez, M., Gören, Ö., Nikbay, M., 2022. Advanced experiments on gaussian process-based multi-fidelity methods over diverse mathematical characteristics, in: 8th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2022), Oslo, Norway. pp. –.

Zaytsev, A., Burnaev, E., 2017. Large scale variable fidelity surrogate modeling. Annals of Mathematics and Artificial Intelligence 81, 167–186.

Zhan, D., Wu, J., Xing, H., Li, T., 2024. A cooperative approach to efficient global optimization. Journal of Global Optimization 88, 327–357.

Álvarez, M.A., Lawrence, N.D., 2008. Sparse convolved Gaussian processes for multi-output regression, in: Advances in Neural Information Processing Systems, pp. 57–64.

Álvarez, M.A., Lawrence, N.D., 2009. Computationally efficient convolved multiple output Gaussian processes. Journal of Machine Learning Research 12, 1459–1500.