# Optimal ANOVA-Based Emulators of Models With(out) Derivatives

Matieyendou Lamboni

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

## To cite this version:

Matieyendou Lamboni. Optimal ANOVA-Based Emulators of Models With(out) Derivatives. Stats, 2025, 8 (1), pp.24. 10.3390/stats8010024. [https://doi.org/10.3390/stats8010024](https://doi.org/10.3390/stats8010024)

HAL Id: hal-05083387

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

Submitted on 24 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Optimal ANOVA-Based Emulators of Models With(out) Derivatives

Matieyendou Lamboni <fnref n="1" />,<fnref n="2" />

**Abstract:** This paper proposes new ANOVA-based approximations of functions and emulators of high-dimensional models using either available derivatives or local stochastic evaluations of such models. Our approach makes use of sensitivity indices to design adequate structures of emulators. For high-dimensional models with available derivatives, our derivative-based emulators reach dimension-free mean squared errors (MSEs) and a parametric rate of convergence (i.e., O(N⁻<fnref n="1" />)). This approach is extended to cope with every model (without available derivatives) by deriving global emulators that account for the local properties of models or simulators. Such generic emulators enjoy dimension-free biases, parametric rates of convergence, and MSEs that depend on the dimensionality. Dimension-free MSEs are obtained for high-dimensional models with particular distributions from the input. Our emulators are also competitive in dealing with different distributions of the input variables and selecting inputs and interactions. Simulations show the efficiency of our approach.

**Keywords:** derivative-based ANOVA; emulators; high-dimensional models; independent input variables; optimal estimators

**PACS:** 62J10; 62L20; 62Fxx; 49Q12; 26D10

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

check for
updates

Received: 31 January 2025

Revised: 1 March 2025

Accepted: 12 March 2025

Published: 17 March 2025

**Citation:** Lamboni, M. Optimal ANOVA-Based Emulators of Models With(out) Derivatives. *Stats* 2025, 8, 24. https://doi.org/10.3390/stats8010024

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Derivatives are sometimes available in modeling either according to the nature of observations of the phenomena of interest ([1,2] and the references therein) or low-cost evaluations of exact derivatives for some classes of PDE/ODE-based models, thanks to adjoint methods [3–8]. Models are defined via their rates of change with respect to their inputs; implicit functions (defined via their derivatives) are instances. Additionally, and particularly for complex models or simulators, efficient estimators of gradients and second-order derivatives using stochastic approximations are provided in refs. [9–13]. Being able to reconstruct functions using such derivatives is worth investigating [14], as having a practical, fast-evaluated model that links stochastic parameters and/or stochastic initial conditions to the output of interest of PDE/ODE models remains a challenge due to numerous uncertain inputs [15].

Moreover, the first-order derivatives of models are used to quickly select non-relevant input variables of simulators or functions, leading to effective screening measures. Efficient variance-based screening methods rely on the upper bounds of generalized sensitivity indices, including Sobol's indices (see refs. [16–22] for independent inputs and [23] for non-independent variables).

<footnote n="1">¹ Department DFR-ST, University of Guyane, Cayenne 97346, French Guiana; matieyendou.lamboni@gmail.com or matieyendou.lamboni@univ-guyane.fr</footnote>

<footnote n="2">² 228-UMR Espace-Dev, University of Guyane, University of Réunion, IRD, University of Montpellier, 34090 Montpellier, France</footnote>
For high-dimensional simulators, dimension reductions via screening input variables are often performed before building their emulators using Gaussian processes (kriging models) [24–28] polynomial chaos expansions [29–31], SS-ANOVA [32,33], or other machine learning approaches [34,35]. Indeed, such emulators rely on nonparametric or semi-parametric regressions and struggle to reconstruct simulators for a moderate to large number of inputs (e.g., [35]). Often, nonparametric rates of convergence are achieved by such emulators (see ref. [33] for SS-ANOVA and [36–38] for polynomial chaos expansions). Regarding the stability and accuracy of polynomial chaos expansions, the number of model runs needed is first estimated at the square of the dimension of the basis used [36] and then reduced at that dimension up to a logarithm factor [37,38]. Note that for *d* inputs, such a dimension is about (*w* + 1)⁴ for the tensor-product basis, including, at most, the monomial of degree *w*.

For models with a moderate to large number of relevant inputs, Bayesian additive regression trees have been used for building emulators of such models using only the input and output observations (see ref. [35] and the references therein). Such approaches rely somehow on rule-ensemble learning approaches by constructing homogenous or local base learners (see ref. [34] and references therein). Combining the model outputs and model derivatives can help to build emulators that account for both local and global properties of simulators. For instance, including derivatives in the Gaussian processes (considered in ref. [2]) allows for improving such emulators. Emulators based on Taylor series (see ref. [14]) combine both the model outputs and derivative outputs with interesting theoretical results, such as dimension-free rates of convergence. However, concrete constructions of such emulators are not provided in that paper.

Note that the aforementioned emulators are part of the class of global approximations of functions. While global emulators can be used to approximate models at any point in entire domain, local or point-based emulators require building different emulators for the same model. Conceptually, the main issues related to such practical emulators are the truncation errors and the biases due to epistemic uncertainty. Indeed, none of the above emulators rely on exact and finite expansions of functions in general. Thus, additional assumptions about the order of such errors are necessary to derive the rates of convergence of emulators. For instance, decreasing eigenvalues of kernels are assumed in kernel-based learning approaches (see, e.g., ref. [39]).

So far, the recent derivative-based (Db) ANOVA provides exact expansions of smooth functions with different types of distributions from inputs using the model outputs, such as first-order and second-order derivatives up to *d*th-order cross partial derivatives [22]. It was used in ref. [13] to derive the plug-in and crude emulators of models by replacing the unknown derivatives with their estimates. However, convergence analysis of such emulators is not known, and derivative-free methods are much more convenient for applications in which the computations of cross-partial derivatives are too expansive or impossible [9–11].

Therefore, for high-dimensional simulators for which all the inputs are important or not, it is worth investigating the development of their emulators based directly on available derivatives or derivative-free methods. The contribution of this paper is threefold:

* We designed adequate structures of emulators based on information gathered from global and derivative-based sensitivity analysis, such as unbiased orders of truncations and the selection of relevant ANOVA components (inputs and interactions);

* We constructed derivative-based or derivative-free global emulators that are easy to fit and compute and can cope with every distribution of continuous input variables;

* We examined the convergence analysis of our emulators, with a particular focus on the (i) dimension-free upper bounds of the biases and MSEs; (ii) the parametric rates of
convergence (i.e., O(N⁻¹)); and (iii) the number of model runs needed to obtain the stability and accuracy of our estimations.

In this paper, flexible emulators of complex models or simulators and approximations of functions are proposed using exact and finite expansions of functions involving cross-partial derivatives, known as Db-ANOVA. Section 2 first deals with the general formulation of Db-ANOVA using different continuous distribution functions and emulators of models with available derivatives. Such emulators reach the parametric rates of convergence, and their MSEs do not depend on dimensionality (i.e., d). Second, adequate structures of emulators are investigated using Sobol's indices and their upper bounds, as the components of such emulators are interpretable as the main effects and interactions of a given order. The orders of unbiased truncations have been derived, leading to the possibility of selecting the ANOVA components that were included in our emulators.

For non-highly smooth models and for high-dimensional simulators for which the computations of derivatives are too expansive or impossible, new, efficient simulator emulators have been considered and are shown in Section 3, including their statistical properties. First, we provide such emulators under the assumption of quasi-uniform distributions of inputs so as to (i) obtain practical conditions for using such emulators and (ii) derive emulators that enjoy dimension-free MSEs for particular distributions of inputs. Second, such an assumption is removed to cope with every distribution of inputs. The proposed emulators have dimension-free biases and reach the parametric rate of convergence as well. Numerical illustrations (see Section 4) and an application to a heat diffusion PDE model (see Section 5) have been considered to show the efficiency of our estimators, and we conclude this work in Section 6.

## General Notation

For an integer, $d > 0$, denote with $\mathbf{X} := (X_1, \dots, X_d)$ a random vector of $d$ independent and continuous variables with marginal cumulative distribution functions (CDFs), $F_j$, and probability density functions (PDFs), $\rho_j$, $j = 1, \dots, d$.

For a non-empty subset, $u \subseteq \{1, \dots, d\}$, $|u|$ stands for its cardinality, and $(\sim u) := \{1, \dots, d\} \setminus u$. Additionally, $\mathbf{X}_u := (X_u, \mathbf{X}_{\sim u})$ denotes a subset of such variables, and the partition $\mathbf{X} = (\mathbf{X}_u, \mathbf{X}_{\sim u})$ holds. Finally, we use $||\cdot||_2$ for the Euclidean norm, $\mathbb{E}[\cdot]$ for the expectation operator, and $\mathbb{V}[\cdot]$ for the variance operator.

## 2. New Insight into Derivative-Based ANOVA and Emulators

Given an integer, $n > 0$, and an open set, $\Omega \subseteq \mathbb{R}^d$, consider a weak partial differentiable function, $f: \Omega \to \mathbb{R}^n$ [40,41], and a subset, $v \subseteq \{1, \dots, d\}$, with $|v| > 0$. Denote with $\mathcal{D}^{|v|}f := (\prod_{k \in v} \frac{\partial}{\partial x_k})f$ the $|v|$th weak cross-partial derivatives of each component of $f$ with respect to each $x_k$ with $k \in v$ and $\mathcal{L}^2(\Omega) := \{f: \Omega \to \mathbb{R}^n: \mathbb{E}[||f(\mathbf{X})||_2^2] < +\infty\}$, e.g., the Hilbert space of functions. Consider the following Hilbert-Sobolev space:

$$W^{d,2} := \{f \in \mathcal{L}^2(\Omega) : \mathcal{D}^{|v|}f \in \mathcal{L}^2(\Omega); \quad \forall |v| \le d\}.$$

In what follows, assume the following:

**Assumption 1.** $\mathbf{X}$ is a random vector of independent and continuous variables, supported on an open domain, $\Omega$.

**Assumption 2.** $f(\cdot)$ is a deterministic function with $f(\cdot) \in W^{d,2}$.
## 2.1. Full Derivative-Based Emulators

Under Assumption 2, every sufficiently smooth function, $f(\cdot)$, admits the derivative-based ANOVA (Db-ANOVA) expansion (see refs. [13,22]), that is, $\forall \mathbf{x} \in \Omega$,

$$f(\mathbf{x}) = \mathbb{E}_{\mathbf{X}'}[f(\mathbf{X}')] + \sum_{v, v \subseteq \{1, \dots, d\} \atop |v| > 0} \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|v|} f(\mathbf{X}') \prod_{k \in v} \frac{G_k(X_k') - \mathbb{I}_{[X_k' \ge x_k]}}{g_k(X_k')} \right], \quad (1)$$

where $\mathbf{X}' := (X'_1, \dots, X'_d)$ is a random vector of independent variables, having the CDFs $X'_j \sim G_j$ and the PDFs $g_j := \frac{dG_j}{dx'_j}$. By evaluating $f(\cdot)$ at the random vector, $\mathbf{X}$, and taking $G_j = F_j$, yields the unique and orthogonal Db-ANOVA decomposition of $f(\mathbf{X})$, that is,

$$f(\mathbf{X}) = \mathbb{E}[f(\mathbf{X}')] + \sum_{v \subseteq \{1, \dots, d\} \atop |v| > 0} f_v(\mathbf{X}_v); \quad f_v(\mathbf{X}_v) := \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|v|} f(\mathbf{X}') \prod_{k \in v} \frac{F_k(X_k') - \mathbb{I}_{[X_k' \ge x_k]}}{\rho_k(X_k')} \right]. \quad (2)$$

When analytical cross-partial derivatives are available, or the derivative datasets are observed (see refs. [1,2] and the references therein), we are able to derive emulators of complex models that are time-demanding, bearing in mind the method of moments. Indeed, given a sample, $\{\mathbf{X}'_i\}_{i=1}^N := \{(X'_{i,1}, \dots, X'_{i,d})\}_{i=1}^N$, from $\mathbf{X}'$, and the associated sample of the (analytical or observed) derivatives' outputs, that is,

$$\left( \left\{ \mathcal{D}^{|v|} f(\mathbf{X}'_i) \right\}_{i=1}^N, \forall v \subseteq \{1, \dots, d\} \right),$$

the consistent (full) emulator or estimator of $f$ at any sample point, $\mathbf{x}$, of $\mathbf{X}$ is

$$\widehat{f}_N(\mathbf{x}) := \frac{1}{N} \sum_{i=1}^{N} \sum_{v, v \subseteq \{1, \dots, d\} \atop |v| \ge 0} \mathcal{D}^{|v|} f(\mathbf{X}'_i) \prod_{k \in v} \frac{F_k(X'_{i,k}) - \mathbb{I}_{[X'_{i,k} \ge x_k]}}{\rho_k(X'_{i,k})}, \quad (3)$$

with $\mathcal{D}^{|v|} f = f$ and $\prod_{k \in \mathcal{O}} c_k := 1$ for every real $c_k$ by convention. We can check that $\widehat{f}_N(\mathbf{x})$ is an unbiased estimator and that it reaches the parametric mean squared error (MSE) rate of convergence, that is, $\mathbb{E} \left[ \left( \widehat{f}_N(\mathbf{x}) - f(\mathbf{x}) \right)^2 \right] = \mathcal{O}(N^{-1})$. This rate is dimension-free, provided that $\mathbb{V} \left[ \sum_{v, v \subseteq \{1, \dots, d\} \atop |v| \ge 0} \mathcal{D}^{|v|} f(\mathbf{X}'_1) \prod_{k \in v} \frac{F_k(X'_{1,k}) - \mathbb{I}_{[X'_{1,k} \ge x_k]}}{\rho_k(X'_{1,k})} \right] < +\infty$.

For complex models without cross-partial derivatives, optimal estimators of such derivatives (i.e., $\widehat{\mathcal{D}^{|v|} f}$) have been used to construct the plug-in consistent emulator of $f$ [13]. Such an emulator is given by

$$\widehat{f}_{N,p}(\mathbf{x}) := \frac{1}{N} \sum_{i=1}^{N} \sum_{v, v \subseteq \{1, \dots, d\} \atop |v| \ge 0} \widehat{\mathcal{D}^{|v|} f}(\mathbf{X}'_i) \prod_{k \in v} \frac{F_k(X'_{i,k}) - \mathbb{I}_{[X'_{i,k} \ge x_k]}}{\rho_k(X'_{i,k})}.$$

While the estimator, $\widehat{\mathcal{D}^{|v|} f}$, provided in ref. [13] has a dimension-free upper bound for the bias and reaches the parametric rate of convergence, its MSE increases with $d^{|v|}$, showing the necessity of using the number of model runs, $N \propto d^{|v|}$, to expect a significant reduction in the MSE for higher-order cross-partial derivatives.
## 2.2. Adequate Structures of Emulators and Truncations

In high-dimensional settings, it is common to expect a reduction in the dimensionality before building emulators. The use of truncations is common practice in the polynomial approximation of functions [29,32,33,37] and in ANOVA-practicing communities [17,42–44], leading to truncated errors. When using Db-ANOVA, controlling such errors is made possible according to the information gathered from global sensitivity analysis [13,17]. Indeed, the variances in the terms in Db-ANOVA expansions of functions are exactly the main part, and interactions with Sobol’s indices up to a normalized constant occur when $G_j = F_j$. Thus, we are able to avoid non-relevant terms in our emulators according to the values of Sobol’s indices, suggesting that adequate truncations will not have any impact on the MSEs and the above parametric rate of convergence. For the sake of simplicity, $n = 1$ is considered in what follows.

**Definition 1.** Consider an integer, $d_0 \in \{1, \dots, d\}$, and the full Db-ANOVA given by (2). The truncated Db-ANOVA of $f$ (in the superpose sense [42–45]) of order $d_0$ is given by

$$f_{T,d_0}(\mathbf{X}) := \mathbb{E}[f(\mathbf{X}')] + \sum_{\substack{v \subseteq \{1, \dots, d\} \\ 0 < |v| \le d_0}} f_v(\mathbf{X}_v).$$

While $f_{T,d_0}$ is an approximation of $f$ in general, the equality holds for some functions. Given an integer, $\alpha \ge 0$, consider the space of functions:

$$\mathcal{L}_{\alpha,0} := \left\{ f : \mathbb{R}^d \to \mathbb{R}^n : \left| f(\mathbf{x}) - \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le \alpha}} \mathbb{E}_{\mathbf{X}'_w} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} \frac{G_k(\mathbf{X}'_k) - \mathbb{I}_{\mathbf{X}'_k \ge \mathbf{x}_k}}{g_k(\mathbf{X}'_k)} \right] \right| = 0 \right\}.$$

We can see that $\mathcal{L}_{\alpha,0}$ is a class of functions having, at most, $\alpha$-order interactions. Moreover, $\mathcal{L}_{\alpha,0}$ contains the class of functions given by

$$\mathcal{D}^{|v|,j_0}|f = 0, \forall v \subseteq \{1, \dots, d\}, j_0 \in (\sim v), |v| = \alpha.$$

It is clear, then, that we have $f(\mathbf{X}) = f_{T,d_0}(\mathbf{X})$ when $f \in \mathcal{L}_{\alpha=d_0,0}$.

**Definition 2.** A truncation of $f$ given by $f_{T,d_0}$ is said to be an unbiased one whenever $f(\mathbf{X}) = f_{T,d_0}(\mathbf{X})$.

Thus, $\mathcal{L}_{\alpha,0}$ is the class of unbiased truncations of $f$ up to the order of $d_0 = \alpha$.

Formally, based on available derivatives, we are able to derive unbiased truncations for some classes of functions. Consider the following Db expressions of Sobol’s indices of the input variables $X_j$, $j = 1, \dots, d$ and their upper bounds (see ref. [22]):

$$S_j := \frac{1}{\mathbb{V}[f(\mathbf{X})]} \mathbb{E} \left[ \frac{\partial f(\mathbf{X})}{\partial x_j} \frac{\partial f(\mathbf{X}')}{\partial x_j} \frac{F_j(\min(X_j, X'_j)) - F_j(X_j)F_j(X'_j)}{\rho_j(X_j)\rho_j(X'_j)} \right],$$

$$S_{T_j} := \frac{1}{\mathbb{V}[f(\mathbf{X})]} \mathbb{E} \left[ \frac{\partial f(\mathbf{X})}{\partial x_j} \frac{\partial f(\mathbf{X}')}{\partial x_j} (\mathbf{X}'_j, \mathbf{X}_{\sim j}) \frac{F_j(\min(\mathbf{X}_j, \mathbf{X}'_j)) - F_j(\mathbf{X}_j)F_j(\mathbf{X}'_j)}{\rho_j(\mathbf{X}_j)\rho_j(\mathbf{X}'_j)} \right],$$
and

$$S_j \le S_{T_j} \le UB_j := \frac{1}{2\mathbb{V}[f(\mathbf{X})]} \mathbb{E} \left[ \left( \frac{\partial f}{\partial x_j}(\mathbf{X}) \right)^2 \frac{F_j(X_j)(1-F_j(X_j))}{[\rho_j(X_j)]^2} \right].$$

Note that the computations of $S_j$ and $UB_j$ are straightforward, using given first-order derivatives (see ref. [13]), whereas those of $S_{T_j}$ require using i) nonparametric methods for given derivatives or ii) derivatives of specific input values. Using such indices, adequate structures of $f(\cdot)$ can be constructed. For instance, it is known that when $\sum_{j=1}^d S_j = 1$, $f$ is an additive function of the form $f(\mathbf{x}) = \sum_{j=1}^d h_j(\mathbf{x}_j)$, with $h_j$ being a real-valued function. Thus, a truncation (in the superpose sense) of order $d_0 = 1$ is an unbiased truncation. Other values of $d_0$ are given in Proposition 1.

**Proposition 1.** Consider the main and total indices given by $S_j$, $S_{T_j}$ with $j = 1, \dots, d$. Then,

* $\sum_{j=1}^d (S_j + S_{T_j}) = 2$ implies that using $d_0 = 2$ leads to unbiased truncations;

* $\sum_{j=1}^d (2S_j + S_{T_j}) \in \{2, 3\}$ implies that using $d_0 = 3$ leads to unbiased truncations;

* If there exists an integer, $\alpha$, such that $(\alpha - 1) \sum_{j=1}^d S_{T_j} + \sum_{j=1}^d S_j \ge \alpha$ and $\sum_{j=1}^d S_{T_j} + (\alpha - 1) \sum_{j=1}^d S_j \le \alpha$, then $d_0 = \alpha$ leads to unbiased truncations.

**Proof.** See Appendix ??.

In general, if $\mathcal{D}^{|\{v, j_0\}|} f = 0$, $\forall v \subseteq \{1, \dots, d\}$, $j_0 \in (\sim v)$, $|v| = \alpha$, then $d_0 = \alpha$ leads to unbiased truncations. Often, low-order derivatives are available or can be efficiently computed using fewer model runs, leading to truncated emulators in the superpose sense. It is worth noting that our emulators still enjoy the above parametric rate of convergence for unbiased truncations. In the presence of truncated errors, it is usually difficult to derive the rates of convergence without additional assumptions about the order of such errors.

In addition to such truncations in the superpose sense, screening measures allow for quickly identifying non-relevant inputs (i.e., $UB_j \approx 0$), leading to possible dimension reductions. For instance, we can see the following:

* $S_j = S_{T_j}$ or $S_j \approx UB_j$ implies removing all cross-partial derivatives or interactions involving $X_j$;

* $S_j = 0$ and $UB_j \ge S_{T_j} \gg 0$, $\forall j \in \{1, \dots, d\}$ suggest removing the first-order terms, corresponding to $d_0 = 1$;

* $1 - (S_i + S_{T_j}) \le S_{T_{\sim\{i,j\}}} = 0$ or, equivalently, $S_i + S_{T_j} = 1$ and $S_{T_k} = 0, \forall k \notin \{i, j\}$ implies keeping only $X_i$ and $X_j$.

For non-highly smooth functions and the models for which the computations of derivatives are too expansive or impossible, derivative-free methods combined with unbiased truncations remain an interesting framework.

# 3. Derivative-Free Emulators of Models

This section covers the development of the model emulators, even when all the inputs are important according to screening measures.

## 3.1. Stochastic Surrogates of Functions Using Db-ANOVA

Consider integers $L > 0$, $q > 0$; $\beta_l \in \mathbb{R}$ with $l = 1, \dots, L$; $h := (h_1, \dots, h_d) \in \mathbb{R}_+^d$, and denote with $V := (V_1, \dots, V_d)$ a $d$-dimensional random vector of independent variables satisfying: $\forall j \in \{1, \dots, d\}$,

$$\mathbb{E}[V_j] = 0; \quad \mathbb{E}[(V_j)^2] = \sigma^2; \quad \mathbb{E}[(V_j)^{2q+1}] = 0; \quad \mathbb{E}[(V_j)^{2q}] < +\infty.$$
Any random variable that is symmetrically about zero is an instance of $V_j$s. Additionally, denote $\beta_\ell h \mathbf{V} := (\beta_\ell h_1 V_1; \dots, \beta_\ell h_d V_d)$. For concise reporting of the results, elementary symmetric polynomials (ESPs) were used (see, e.g., [46,47]).

**Definition 3.** Given $u \subseteq \{1, \dots, d\}$ with $|u| > 0$ and $\mathbf{r}_u := (r_k, \forall k \in u) \in \mathbb{R}^{|u|}$, the $p^{th}$ ESP of $\mathbf{r}_u$ is defined as follows:

$$ \mathbf{e}_p^{(u)}(\mathbf{r}_u) := \begin{cases} 0 & \text{if } p > |u| \text{ or } p < 0 \\ 1 & \text{if } p = 0 \\ \sum_{w \subseteq u} \prod_{k \in w} r_k & \text{if } p = 1, \dots, |u| \end{cases}. $$

Note that $\mathbf{r} := \mathbf{r}_{\{1, \dots, d\}} = (r_1, \dots, r_d) \in \mathbb{R}^d$, $\mathbf{e}_p^{(1:d)}(\mathbf{r}) := \mathbf{e}_p^{(\{1, \dots, d\})}(\mathbf{r})$. In addition, given $X'_j \sim G_j, \forall j \in \{1, \dots, d\}$, define

$$ R_k(x_k, X'_k, V_k) := \frac{G_k(X'_k) - \mathbb{I}[X'_k \ge x_k]}{g_k(X'_k) h_k \sigma^2} V_k, \quad k = 1, \dots, d; $$

$$ \mathbf{R}_u(\mathbf{x}_u, \mathbf{X}'_u, \mathbf{V}_u) := (R_k(x_k, X'_k, V_k), \forall k \in u); \quad \mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V}) := \mathbf{R}_{\{1, \dots, d\}}(\mathbf{x}, \mathbf{X}', \mathbf{V}). $$

Without loss of generality, we are going to focus on the modified output, that is, $f^c(\mathbf{x}) := f(\mathbf{x}) - \mathbb{E}[f(\mathbf{X}')]$. Based on the above framework, Theorem 1 provides a new approximation of every function or surrogate of a deterministic simulator.

**Theorem 1.** Consider distinct $\beta_\ell$s. If $f$ is smooth enough and Assumption 1 holds, then there exists $\alpha_d \in \{1, \dots, L\}$ and real coefficients $C_1^{(p)}, \dots, C_L^{(p)}, p = 1, \dots, d$ such that

$$ f^c(\mathbf{x}) = \sum_{\ell=1}^{L} \sum_{p=1}^{d} C_{\ell}^{(p)} \mathbb{E} [f(\mathbf{X}' + \beta_{\ell} h \mathbf{V}) \mathbf{e}_{p}^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V}))] + \mathcal{O} (\|\mathbf{h}\|_2^{2\alpha_d}). \quad (4) $$

**Proof.** See Appendix B. □

The setting $L = 1, \beta_1 = 1, C_1^{(p)} = 1$ with $p = 1, \dots, d$ provides an approximation of order $\mathcal{O}(\|\mathbf{h}\|_2^2)$. Equivalently, the same order is obtained using the constraints:

$$ \left\{ \begin{array}{ll} \sum_{\ell=1}^{L} C_{\ell}^{(p)} \beta_{\ell}^{r} = \delta_{p,r}; & r = 0, \dots, L-1, \quad \text{if } p \le L-1 \\ \sum_{\ell=1}^{L} C_{\ell}^{(p)} \beta_{\ell}^{r} = \delta_{p,r}; & r = 0, \dots, L-2, p, \quad \text{otherwise} \end{array} \right. $$

with $p = 1, \dots, d$ and $L \le d+1$. In the same sense, taking $\sum_{\ell=1}^{L} C_{\ell}^{(p)} \beta_{\ell}^{r} = \delta_{p,r}$; $r = p, \dots, L+p-1$ with $p = 1, \dots, d$ yields an approximation of order $\mathcal{O}(\|\mathbf{h}\|_2^{2L})$. Such constraints implicitly define the coefficients $C_1^{(p)}, \dots, C_L^{(p)}, p = 1, \dots, d$, and they rely on the (generalized) Vandermonde matrices. Distinct values of $\beta_\ell$s (i.e., $\beta_{\ell_1} \ne \beta_{\ell_2}$) ensure the existence and uniqueness of such coefficients, as such matrices are invertible (see refs. [47,48]).

To improve the approximations of lower-order terms (i.e., lower values of $p$) in Equation (4), we are given an integer $r^* \in \{0, \dots, d-1\}$ with $r^* \le L-2$ and consider the following constraints:

$$ \left\{ \begin{array}{ll} \sum_{\ell=1}^{L} C_{\ell}^{(p)} \beta_{\ell}^{r} = \delta_{p,r}; & r = 0, \dots, r^*, p + 2\lambda_p + 2, \dots, p + 2\lambda_p + 2(L-1-r^*), \quad \text{if } p \le r^* \\ \sum_{\ell=1}^{L} C_{\ell}^{(p)} \beta_{\ell}^{r} = \delta_{p,r}; & r = 0, \dots, r^*, p, p + 2, \dots, p + 2(L-r^*-2), \quad \text{otherwise} \end{array} \right. \quad (5) $$

where $\lambda_p := \lfloor \frac{r^*-p}{2} \rfloor$ stands for the largest integer that is less than $\frac{r^*-p}{2}$. The above choice of coefficients requires $L \ge r^* + 2$, and $L^* := r^* + 2$ is the minimum number of model runs used for deriving surrogates of functions. Such coefficients are more suitable for truncated
surrogates. For instance, the truncated surrogate of order $d_0 \le d$ (in the superposition sense) is given by

$$ \tilde{f}_{T,d_0}^{c}(\mathbf{x}) := \sum_{\ell=1}^{L} \sum_{p=1}^{d_0} C_{\ell}^{(p)} \mathbb{E} [f(\mathbf{x}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \mathbf{e}_{p}^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V}))]. $$

Note that in this case, one must require $r^* \in \{0, \dots, d_0 - 1\}$, and we will see that $r^* = d_0 - 1$ is the best choice to improve the MSEs.

Likewise, when $\mathbf{X}_{uI}$ with $uI \subset \{1, \dots, d\}$ is the vector of the most influential input variables according to variance-based sensitivity analysis (see Section 2.2), the following truncated surrogate should be considered:

$$ \tilde{f}_{uI}^{c}(\mathbf{x}_{uI}) := \sum_{\ell=1}^{L} \sum_{p=1}^{|uI|} C_{\ell}^{(p)} \mathbb{E} [f(\mathbf{x}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \mathbf{e}_{p}^{(uI)}(\mathbf{R}_{uI}(\mathbf{x}_{uI}, \mathbf{X}'_{uI}, \mathbf{V}_{uI}))]. $$

Based on Equation (4), the method of moments allows for deriving the emulator of any simulator or the estimator of any function. To that end, we are given two independent samples of size $N$, that is, $\{\mathbf{X}'_i\}_{i=1}^N := \{(\mathbf{X}'_{i,1}, \dots, \mathbf{X}'_{i,d})\}_{i=1}^N$ from $\mathbf{X}'$ and $\{\mathbf{V}_i\}_{i=1}^N := \{(V_{i,1}, \dots, V_{i,d})\}_{i=1}^N$ from $\mathbf{V}$. The full and consistent emulator is given by

$$ \widehat{f}_{N}^{c}(\mathbf{x}) := \frac{1}{N} \sum_{i=1}^{N} \sum_{\ell=1}^{L} \sum_{p=1}^{d} C_{\ell}^{(p)} f(\mathbf{X}'_{i} + \beta_{\ell} \mathbf{h} \mathbf{V}_{i}) \mathbf{e}_{p}^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}'_{i}, \mathbf{V}_{i})). \quad (6) $$

The derivations of the truncated emulators (i.e., $\tilde{f}_{T,d_0}^{c}$ and $\tilde{f}_{uI}^{c}$) are straightforward. All these emulators rely on NL model runs with the possibility $L \ll d$. This property is useful for high-dimensional simulators.

## 3.2. Statistical Properties of Our Emulators

While the emulator $\widehat{f}_{N}^{c}(\mathbf{x})$ does not rely on the model derivatives, structural and technical assumptions about $f$ are needed to derive the biases of this emulator, such as the Hölder space of functions. Given $\tilde{\mathbf{r}} := (i_1, \dots, i_d) \in \mathbb{N}^d$, denote $\mathcal{D}^{(\tilde{\mathbf{r}})} f := (\prod_{k=1}^d \frac{\partial^{i_k}}{\partial x_k}) f$, $(\mathbf{x})^{\tilde{\mathbf{r}}} := \prod_{k=1}^d x_k^{i_k}, \tilde{\mathbf{r}}! = i_1! \dots i_d!$ and $||\tilde{\mathbf{r}}||_1 = i_1 + \dots + i_d$. Given $\alpha \ge 0$, the Hölder space of $\alpha$-smooth functions is given by $\forall \mathbf{x}, \mathbf{y} \in \mathbb{R}^d$,

$$ \mathcal{H}_{\alpha} := \left\{ f : \mathbb{R}^{d} \to \mathbb{R} : \left| f(\mathbf{x}) - \sum_{\|\tilde{\mathbf{r}}\|=1}^{\alpha-1} \frac{\mathcal{D}^{(\tilde{\mathbf{r}})} f(\mathbf{y})}{\tilde{\mathbf{r}}!} (\mathbf{x} - \mathbf{y})^{\tilde{\mathbf{r}}} \right| \le M_{\alpha} ||\mathbf{x} - \mathbf{y}||_{2}^{\alpha} \right\}, $$

with $M_{\alpha} > 0$, and $\mathcal{D}^{(\tilde{\mathbf{r}})} f(\mathbf{y})$ is a (weak) cross-partial derivative.

Moreover, given $B_{\alpha} \ge 0$ and CDFs $G_{js}$, define the following space of functions:

$$ \mathcal{L}_{\alpha, B_{\alpha}} := \left\{ f : \mathbb{R}^d \to \mathbb{R} : \left| f(\mathbf{x}) - \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le \alpha}} \mathbb{E}_{\mathbf{X}'_w} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} \frac{G_k(\mathbf{X}'_k) - \mathbb{I}_{\mathbf{X}'_k \ge \mathbf{X}_k}}{g_k(\mathbf{X}'_k)} \right] \right| \le B_{\alpha} \right\}. $$

We can see that $\mathcal{L}_{\alpha, B_{\alpha}}$ contains constants; $\mathcal{L}_{\alpha, 0}$ is a class of functions having, at most, $\alpha$-order of interactions, and $\mathcal{L}_{d, 0}$ is a class of all smooth functions, as $B_d = 0$. Lemma 1 provides the links between both spaces. To that end, consider $M'_{|w|} := ||\mathcal{D}^{|w|} f||_{\infty}$ for all $w \subseteq \{1, \dots, d\}$ with $||\cdot||_{\infty}$ being the infinity norm.

**Assumption 3.** $g_j \ge \rho_{\min} > 0$ for any $j \in \{1, \dots, d\}$.
Assumption 3 aims to cover the class of quasi-uniform distributions and other distributions for which the event $g_j \ge \rho_{\min}$ occurs with a high probability. It is the case for most unbounded distributions.

**Lemma 1.** Consider $0 < d_0 \le d$, and assume that $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, d\}$ and Assumptions 1 and 3 hold. Then, there exists $\gamma_0 > 0$ such that $f \in \mathcal{L}_{d_0, D_{d_0, \rho_{\min}}}$, with

$$D_{d_0, \rho_{\min}} := 2\gamma_0 M_0' \left[ \frac{1}{2\rho_{\min}} \left( \frac{M_0'}{M_0'} \right)^{1/d} + 1 \right]^{d_0} \left\{ \left[ \frac{1}{2\rho_{\min}} \left( \frac{M_0'}{M_0'} \right)^{1/d} + 1 \right]^{d-d_0} - 1 \right\}.$$

**Proof.** See Appendix C. □

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

Note that Lemma 1 also provides the upper bound of the remaining terms when approximating $f(\mathbf{x})$ using the truncated function

$$f_{d_0}(\mathbf{x}) := \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le d_0}} \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} \frac{G_k(X_k') - \mathbb{I}_{X_k' \ge x_k}}{g_k(X_k')} \right].$$

When $\frac{1}{2\rho_{\min}} \left(\frac{M_0'}{M_0'}\right)^{1/d} \to 0$, $D_{d_0, \rho_{\min}}$ is equivalent to

$$D_{d_0, \rho_{\min}} = \frac{(d-d_0)\gamma_0 (M_0')^{(d-1)/d} (M_0')^{1/d}}{\rho_{\min}} \left[ \frac{d_0}{2\rho_{\min}} \left( \frac{M_0'}{M_0'} \right)^{1/d} + 1 \right].$$

### 3.2.1. Biases of the Proposed Emulators

To derive the bias of $\widehat{f}_{N,d_0}^c$ (i.e., an estimator of $f_{d_0}^c$) in Theorem 2 using the aforementioned spaces of functions, denote with $\mathbf{Z} := (Z_1, \dots, Z_d)$ a $d$-dimensional random vector of independent variables that are centered about zero and standardized (i.e., $\mathbb{E}[Z_k^2] = 1, k = 1, \dots, d$), and $\mathcal{R}_c$ denotes the set of such random vectors. For any $r \in \mathbb{N}$ and $w \subseteq \{1, \dots, d\}$, define

$$\Gamma_r := \sum_{\ell=1}^{L} |C_{\ell}^{(|w|)} \beta_{\ell}^r|; \qquad K_{w,L} := \inf_{\mathbf{Z} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| \mathbf{Z}^2 \right\|_2^L \prod_{k \in w} Z_k^2 \right] \Gamma_{|w|+2L}; \\ L_w' := \left( \left[ \frac{r^* - |w|}{2} \right] + L - r^* \right) \mathbb{I}_{|w| \le r^*} + (L - r^* - 1) \mathbb{I}_{|w| > r^*}.$$

**Theorem 2.** Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2(L - r^* - 1))\}$ and Assumptions 1 and 3 hold. Then, we have

$$\left| \mathbb{E} \left[ \widehat{f}_{N, d_0}^c (\mathbf{x}) \right] - f^c (\mathbf{x}) \right| \le \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \sigma^{2L_w'} \left\| \mathbf{h}^2 \right\|_2^{L_w'} M_{|w|+2L_w'} K_{w, L_w'} \left( \frac{1}{2\rho_{\min}} \right)^{|w|} + D_{d_0, \rho_{\min}}. \quad (7)$$

Moreover, if $V_k \sim \mathcal{U}(-\zeta, \zeta)$ with $\zeta > 0$ and $k = 1, \dots, d$, then

$$\left| \mathbb{E} \left[ \widehat{f}_{N, d_0}^c (\mathbf{x}) \right] - f^c (\mathbf{x}) \right| \le \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \xi^{2L_w'} \left\| \mathbf{h}^2 \right\|_1^{L_w'} M_{|w|+2L_w'} \Gamma_{|w|+2L_w'} \left( \frac{1}{2\rho_{\min}} \right)^{|w|} + D_{d_0, \rho_{\min}}. \quad (8)$$

**Proof.** See Appendix D. □

<figure><img src="image_5.png" /></figure>
Using the fact $h_k \to 0$, the results provided in Theorem 2 have simple upper bounds (see Corollary 1). To provide such results, consider

$$K_{1,r^*,d_0}^{\max} := \max_{\substack{w \subseteq \{1, \dots, d\} \\ r^* < |w| \le d_0}} \{K_{w,(L-r^*-1)} M_{|w|+2(L-r^*-1)}\};$$

$$K_{2,r^*,d_0}^{\max} := \max_{\substack{w \subseteq \{1, \dots, d\} \\ r^* < |w| \le d_0}} \{M_{|w|+2(L-r^*-1)} \Gamma_{|w|+2(L-r^*-1)}\};$$

$$K_{1,\rho_{\min},r^*} := \left[ 2\rho_{\min} \left( \frac{d}{2\rho_{\min}} \right)^{r^*+1} \frac{\left( \frac{d}{2\rho_{\min}} \right)^{d_0-r^*} - 1}{d - 2\rho_{\min}} \right] \mathbb{I}_{r^* < d_0 - 1} + \left( \frac{d}{d_0} \right) \left( \frac{1}{2\rho_{\min}} \right)^{d_0} \mathbb{I}_{r^* = d_0 - 1}.$$

**Corollary 1.** Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2(L - r^* - 1))\}$ and Assumptions 1 and 3 hold. If $h_k \to 0$, then

$$\left| \mathbb{E} \left[ \widehat{f}_{N,d_0}^c(\mathbf{x}) \right] - f^c(\mathbf{x}) \right| \le \left\| \mathbf{h}^2 \right\|_2^{L-r^*-1} \sigma^{2(L-r^*-1)} K_{1,r^*,d_0}^{\max} K_{1,\rho_{\min},r^*} + D_{d_0,\rho_{\min}} + \mathcal{O}\left( \left\| \mathbf{h}^2 \right\|_2^{L-r^*-1} \right). \quad (9)$$

Moreover, if $V_k \sim \mathcal{U}(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, then

$$\left| \mathbb{E} \left[ \widehat{f}_{N,d_0}^c(\mathbf{x}) \right] - f^c(\mathbf{x}) \right| \le \left\| \mathbf{h}^2 \right\|_1^{L-r^*-1} \xi^{2(L-r^*-1)} K_{2,r^*,d_0}^{\max} K_{1,\rho_{\min},r^*} + D_{d_0,\rho_{\min}} + \mathcal{O}\left( \left\| \mathbf{h}^2 \right\|_1^{L-r^*-1} \right). \quad (10)$$

**Proof.** See Appendix E. $\square$

Using the above results, the bias of the full emulator of $f^c$ is straightforward when taking $d_0 = d$ and knowing that $D_{d_0,\rho_{\min}} = 0$. Moreover, Corollary 2 provides the bias of such an emulator under different structural assumptions about $f$ so as to cope with many functions. To that end, define

$$K_{1:d} := \inf_{\mathbf{Z} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| \mathbf{Z} \right\|_2 \prod_{k=1}^d Z_k^2 \right].$$

**Corollary 2.** Let $d_0 = d$; $r^* \le d - 1$ and $L = r^* + 2$. Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, d + 1\}$ and Assumptions 1 and 3 hold. If $h_k \to 0$, then

$$\left| \mathbb{E} \left[ \widehat{f}_{N}^{c}(\mathbf{x}) \right] - f^{c}(\mathbf{x}) \right| \le \sigma \left\| \mathbf{h} \right\|_{2} M_{d+1} K_{1:d} \Gamma_{d+1} \prod_{k=1}^{d} \mathbb{E}[|E_k|] + \mathcal{O}\left( \left\| \mathbf{h} \right\|_{2} \right). \quad (11)$$

Moreover, if $V_k \sim \mathcal{U}(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, then

$$\left| \mathbb{E} \left[ \widehat{f}_{N}^{c}(\mathbf{x}) \right] - f^{c}(\mathbf{x}) \right| \le \xi \left\| \mathbf{h} \right\|_{1} M_{d+1} \Gamma_{d+1} \prod_{k=1}^{d} \mathbb{E}[|E_k|] + \mathcal{O}\left( \left\| \mathbf{h} \right\|_{1} \right). \quad (12)$$

**Proof.** See Appendix F. $\square$

In view of Corollaries 1 and 2, Equation (12) can lead to a dimension-free upper bound of the bias. Indeed, using the uniform bandwidth $h_k = h$ and

$$\xi \le \frac{1}{d M_{d+1} \Gamma_{d+1} \prod_{k=1}^{d} \mathbb{E}[|E_k|]}$$

we can see that the upper bound of the bias is $\left| \mathbb{E} \left[ \widehat{f}_{N}^{c}(\mathbf{x}) \right] - f^{c}(\mathbf{x}) \right| \le h$. Furthermore, it is worth noting that $\left| \mathbb{E} \left[ \widehat{f}_{N,d_0}^{c}(\mathbf{x}) \right] - f^{c}(\mathbf{x}) \right| \le h^{L-r^*-1}$ for any function $f \in \mathcal{L}_{d_0,0}$.
### 3.2.2. Mean Squared Errors

We start this section with the variance of the proposed emulators, followed by their mean squared errors and different rates of convergence. For the variance, define

$$F_{r^*,d_0}^{\max} := \max_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le d_0}} \{M_{\min(r^*,|w|-1)+1} \Gamma_{\min(r^*,|w|-1)+1}\}; \quad F_{d_0}^{\max} := \max_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le d_0}} \{M_{|w|} \Gamma_{|w|}\}.$$

**Theorem 3.** Consider the coefficients given by Equation (5), and assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2(L - r^* - 1))\}$ and Assumptions 1–3 hold. Then,

$$\mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] \le \frac{(F_{r^*,d_0}^{\max})^2}{N} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \prod_{k \in w} \left( \frac{\mathbb{E}[E_k^2]}{h_k^2 \sigma^4} \mathbb{E}[V_1^2 || h \mathbf{V} ||_2^2 \frac{\min(r^*,|w|-1)+1}{|w|}] \right).$$

Moreover, if $r^* = d_0 - 1$, $h_k = h$ and $Z_k = V_k/\sigma$, then

$$\mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] \le \frac{d(F_{d_0}^{\max})^2 \mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{N(d\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] - 3\rho_{\min}^2)} \left[ \left( \frac{d\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} - 1 \right].$$

**Proof.** See Appendix G. □

It turns out that the upper bounds of the variance in Theorem 3 do not depend on the uniform bandwidths when $r^* = d_0 - 1$, leading to the derivations of the parametric MSEs of $\widehat{f}_{N,d_0}^c(\mathbf{x})$ and $\widehat{f}_N^c(\mathbf{x})$. To that end, consider the upper bound of the above variance, that is,

$$\bar{\Sigma}_{d_0} := \frac{d(F_{d_0}^{\max})^2 \mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{N(d\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] - 3\rho_{\min}^2)} \left[ \left( \frac{d\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} - 1 \right].$$

**Remark 1.** Based on the expression of $\bar{\Sigma}_{d_0}$, the random variable $V_j$ or $Z_j = V_j/\sigma$ having the smallest value of fourth moments or kurtosis should be used. Under the additional condition $\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] \ge 3\rho_{\min}^2$, we can check that (see Appendix G)

$$\bar{\Sigma}_{d_0} \le \frac{(F_{d_0}^{\max})^2}{N} \left[ 2 \left( d \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} \mathbb{I}_{d_0 \le d_0^*} + 2^d \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} \mathbb{I}_{d_0 > d_0^*} \right],$$

with $d_0^* := \frac{(d-1)\ln(2)}{\ln(d)}$.

In what follows, we are going to use $\mathbb{E}[(\widehat{f}_{N,d_0}^c(\mathbf{x}) - f^c(\mathbf{x}))^2]$ for the MSE of $\widehat{f}_{N,d_0}^c$ and

$$\mathbb{E}[(\widehat{f}_{N,d_0}^c - f^c)^2] := \mathbb{E}_{\mathbf{x}} \{ \mathbb{E}[(\widehat{f}_{N,d_0}^c(\mathbf{x}) - f^c(\mathbf{x}))^2] \},$$

for the expected or integrated MSE (IMSEs) of $\widehat{f}_{N,d_0}^c$.

**Corollary 3.** Given (5), $r^* = d_0 - 1$, assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2(L - d_0))\}$; $h_k = h \to 0$ and Assumptions 1–3 hold. Then, the MSE and IMSE share the same upper bound, given as follows:
$$ \mathbb{E} \left[ \left( \widehat{f}_{N,d_0}^c - f^c \right)^2 \right] \le 2D_{d_0, \rho_{\min}} \left\| \mathbf{h}^2 \right\|_2^{L-d_0} \sigma^{2(L-d_0)} K_{1,d_0-1,d_0}^{\max} \binom{d}{d_0} \left( \frac{1}{2\rho_{\min}} \right)^{d_0} + D_{d_0, \rho_{\min}}^2 \\ + \bar{\Sigma}_{d_0} + \left\| \mathbf{h}^2 \right\|_2^{2(L-d_0)} \sigma^{4(L-d_0)} \left( K_{1,d_0-1,d_0}^{\max} \right)^2 \binom{d}{d_0}^2 \left( \frac{1}{2\rho_{\min}} \right)^{2d_0}. $$

Moreover, if $G_j = F_j$ with $j = 1, \dots, d$, then the IMSE is given by

$$ \mathbb{E} \left\{ \left[ \left( \widehat{f}_{N,d_0}^c - f^c \right)^2 \right] \right\} \le \left\| \mathbf{h}^2 \right\|_2^{2(L-d_0)} \sigma^{4(L-d_0)} \left( K_{1,d_0-1,d_0}^{\max} \right)^2 \binom{d}{d_0}^2 \left( \frac{1}{2\rho_{\min}} \right)^{2d_0} \\ + D_{d_0, \rho_{\min}}^2 + \bar{\Sigma}_{d_0}. $$

**Proof.** See Appendix H. $\square$

The presence of $D_{d_0, \rho_{\min}}$ in Corollary 3 is going to decrease the rates of convergence of our estimators without additional assumptions about $D_{d_0, \rho_{\min}}^2$. Corollary 4 starts providing such conditions and the associated MSEs and IMSEs.

**Corollary 4.** Under the conditions of Corollary 3, assume that $f \in \mathcal{L}_{d_0,0}$. Then, the upper bound of the IMSE and MSE is given by

$$ \left\| \mathbf{h}^2 \right\|_2^{2(L-d_0)} \sigma^{4(L-d_0)} \left( K_{1,d_0-1,d_0}^{\max} \right)^2 \binom{d}{d_0}^2 \left( \frac{1}{2\rho_{\min}} \right)^{2d_0} + \bar{\Sigma}_{d_0}. $$

Moreover, if $V_k \sim U(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, then this bound becomes

$$ \left\| \mathbf{h}^2 \right\|_1^{2(L-d_0)} \xi^{4(L-d_0)} \left( M_{d_0+2(L-d_0)} \Gamma_{d_0+2(L-d_0)} \right)^2 \binom{d}{d_0}^2 \left( \frac{1}{2\rho_{\min}} \right)^{2d_0} + \bar{\Sigma}_{d_0}. $$

**Proof.** Using Corollary 1, the results are straightforward. $\square$

Based on the upper bounds of Corollary 4, interesting choices of $\sigma$ or $\xi$ on one hand and $\rho_{\min}$ and $h$ on the other hand help in obtaining the parametric rates of convergence due to the fact that $\bar{\Sigma}_{d_0}$ does not depend on $h$.

**Corollary 5.** Let $r^* = d_0 - 1$, $L = d_0 + 1$. Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2)\}$; $f \in \mathcal{L}_{d_0,0}$ and Assumptions 1–3 hold. If $h_k = h \propto N^{-\eta}$ with $\eta \in ]\frac{1}{4}, 1[$ and $\xi^2 \le \left( dM_{d_0+2}\Gamma_{d_0+2}\binom{d}{d_0}\left(\frac{1}{2\rho_{\min}}\right)^{d_0} \right)^{-1}$, then the upper bound of MSE and IMSE is

$$ \mathbb{E} \left[ \left( \widehat{f}_{N,d_0}^c - f^c \right)^2 \right] \le \frac{\left( F_{d_0}^{\max} \right)^2}{N} \left[ 2 \left( \frac{d(d+0.8)}{3\rho_{\min}^2} \right)^{d_0} \mathbb{I}_{d_0 \le d_0^*} + 2^d \left( \frac{d+0.8}{3\rho_{\min}^2} \right)^{d_0} \mathbb{I}_{d_0 > d_0^*} \right] + \mathcal{O}(N^{-1}), $$

provided that $d + 0.8 \ge 3\rho_{\min}^2$.

Moreover, if $\rho_{\min} = \sqrt{\frac{c_0 d(d+0.8)}{3}}$ with the real $1 < c_0 \le 2$, then

$$ \mathbb{E} \left[ \left( \widehat{f}_{N,d_0}^c - f^c \right)^2 \right] \le \frac{\left( F_{d_0}^{\max} \right)^2}{N(c_0 - 1)} + \mathcal{O}(N^{-1}). $$

**Proof.** See Appendix I. $\square$

It is worth noting that the parametric rates of convergence are reached for any function that belongs to $\mathcal{L}_{d_0,0}$ (see Corollary 5). Moreover, taking $c_0 \in ]1, 2]$ leads to the dimension-free MSEs, which hold for particular distributions of the inputs given by
$$X_j \sim F_j^*, \quad \rho_j^* \ge \rho_{\min}^* := \sqrt{\frac{c_0 d(d+0.8)}{3}}, \quad \forall j \in \{1, \dots, d\}.$$

For uniform distributions, $X_j \sim U(a_j, b_j)$, we must have $b_j - a_j \le \frac{1}{\rho_{\min}^*}$. Obviously, such conditions are a bit strong, as a few distributions are covered. In the same sense, using $\widehat{f}_{N,d_0}^c$ as an emulator of $f^c$ for any sample point, $\mathbf{x}$, of $\mathbf{X} \sim F$ requires choosing $X_j^j$ such that its support contains that of $X_j$, $j = 1, \dots, d$. Thus, Assumption 3, given by $g_j > \rho_{\min}$, implicitly depends on the distribution of $X_j$. For instance, given the bounded support of $X_j^j$, that is, $(a_j^j, b_j^j)$, we must have $\frac{1}{\rho_{\min}^*} > b_j^j - a_j^j \ge b_j - a_j$, with $(a_j, b_j)$ being the support of $X_j$, limiting our ability to deploy $\widehat{f}_{N,d_0}^c$ as a dimension-free, global emulator for some distributions of inputs.

Nevertheless, the assumption $g_j \ge \rho_{\min}^*$ is always satisfied for the finite dimensionality, $d$, if we are only interested in estimating $f(\mathbf{x}_0)$ for a given point $\mathbf{x}_0$, leading to local emulators. Indeed, taking $X_j^j \sim G_j$ to be depended on the point $\mathbf{x}^0$ at which $f$ must be evaluated allows for enjoying the parametric rate of convergence and dimension-free MSEs for sample points falling in a neighborhood of $\mathbf{x}_0$. An example of such a choice is

$$X_j^j \sim U\left(x_j^0 - \frac{1}{2\rho_{\min}^*}, x_j^0 + \frac{1}{2\rho_{\min}^*}\right).$$

However, different emulators are going to be built in order to estimate $f(\mathbf{x})$ for any value, $\mathbf{x}$, of $\mathbf{X}$. Constructions of balls of given nodes and the radius $1/\rho_{\min}^*$ are an interesting perspective.

**Remark 2.** When $f \notin L_{d_0,0}$, in-depth structural assumptions about $f$ that should allow for enjoying the above MSEs concern the truncation error, resulting from keeping only all the $|v|^h$ interactions or cross-partial derivatives with $|v| \le d_0$. One way to handle it consists of choosing $d_0 = d_{0,N}$ such that the residual bias is less than $1/N$ (i.e., $D_{d_0,N}^2/\rho_{\min} < 1/N$), thanks to sensitivity indices.

While truncations are sometimes necessary in higher dimensions, it is interesting to have the rates of convergence without any truncation to cover lower or moderate dimensional functions, for instance.

**Corollary 6.** Let $d_0 = d$; $r^* = d - 1$, $L = d + 1$ and $h_k = h$. If $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, d+1\}$, $h \to 0$, and Assumptions 1–3 hold, then the upper bound of MSE and IMSE is

$$\mathbb{E}\left[(\widehat{f}_N^c - f^c)^2\right] \le \sigma^2 ||\mathbf{h}||_2^2 M_{d+1}^2 K_{1:d}^2 \Gamma_{d+1}^2 \left(\frac{1}{2\rho_{\min}}\right)^{2d} + \frac{(F_d^{\max})^2}{N} \left[\left(\frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} + 1\right)^d - 1\right].$$

Moreover, if $V_k \sim U(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, then

$$\mathbb{E}\left[(\widehat{f}_N^c - f^c)^2\right] \le \xi^2 ||\mathbf{h}||_1^2 M_{d+1}^2 \Gamma_{d+1}^2 \left(\frac{1}{2\rho_{\min}}\right)^{2d} + \frac{(F_d^{\max})^2}{N} \left[\left(\frac{d+0.8}{3\rho_{\min}^2} + 1\right)^d - 1\right].$$

**Proof.** See Appendix J. $\square$

In the case of the full emulator of $f^c$, remark that

$$\frac{(F_d^{\max})^2}{N} \left[ \left( \frac{d+0.8}{3\rho_{\min}^2} + 1 \right)^d - 1 \right] \le \frac{d(F_d^{\max})^2(d+0.8)}{N(d(d+0.8) - 3\rho_{\min}^2)} \left[ \left( d \frac{d+0.8}{3\rho_{\min}^2} \right)^d - 1 \right].$$
Based on Corollary 6, different rates of convergence can be obtained depending on the support of the input variables $X_j \sim F_j$ via the choice of $\rho_{\min}$.

**Corollary 7.** Let $r^* = d-1$ and $L = d+1$. Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, d+1\}$; $\zeta \le [dM_{d+1}\Gamma_{d+1}(\frac{1}{2\rho_{\min}})^d]^{-1}$; $h_k = h \propto N^{-\eta}$ with $\eta \in ]\frac{1}{2}, 1[$ and Assumptions 1–3 hold. Then, the (MSE and IMSE) rates of convergence are given as follows:

$$ \mathbb{E}\left[(\hat{f}_N^c - f^c)^2\right] \le \frac{(F_d^{\max})^2}{N} \left[ \left(\frac{d+0.8}{3\rho_{\min}^2} + 1\right)^d - 1 \right] + \mathcal{O}(N^{-1}). $$

$$ (i) \quad \text{If } \rho_{\min} = \sqrt{\frac{d+0.8}{3}}, \text{ then } \mathbb{E}\left[(\hat{f}_N^c - f^c)^2\right] \le \frac{(F_d^{\max})^2}{N}(2^d - 1) + \mathcal{O}(N^{-1}). $$

$$ (ii) \quad \text{If } \rho_{\min} = \sqrt{\frac{c_0 d(d+0.8)}{3}} \text{ with } 1 \le c_0 \le 2, \text{ then} \\ \mathbb{E}\left[(\hat{f}_N^c - f^c)^2\right] \le \frac{(F_d^{\max})^2}{N} \min\left\{\frac{1}{c_0 - 1}, \left(\frac{1}{c_0 d} + 1\right)^d - 1\right\} + \mathcal{O}(N^{-1}), \\ \text{with } \left(\frac{1}{c_0 d} + 1\right)^d - 1 \equiv \frac{1}{c_0} \text{ when } c_0 d \to \infty. $$

**Proof.** It is straightforward bearing in mind Corollaries 5 and 6. $\square$

Again, the assumptions in Points (i)–(ii) are satisfied for fewer distributions, but they are always satisfied if we are only interested in estimating $f(x_0)$ for a given point $x_0$, leading to building different local emulators.

### 3.2.3. Mean Squared Errors for Every Distribution of Inputs

In this section, we are going to remove the assumption about the quasi-uniform distribution of inputs (Assumption 3) so as to cover any probability distribution of inputs. Note that Assumption 3 is used to derive $\mathbb{E}[E_k^2] \le \frac{1}{3\rho_{\min}^2}$ and $\mathbb{E}[|E_k|] \le \frac{1}{2\rho_{\min}^2}$. Such an assumption can be avoided by using the following inequalities:

$$ \mathbb{E}[|E_k|] = \mathbb{E}\left[\frac{[1 - G(X'_k)] \mathbb{I}_{X'_k \ge X_k} + G(X'_k) \mathbb{I}_{X'_k < X_k}}{g_k(X'_k)}\right] \le \kappa_1, $$

with

$$ \kappa_1 := \sup_{k \in \{1, \dots, d\}} \sup_{x'_k \in \Omega_k} \mathbb{E} \left[ \frac{F_k(X'_k) + G(X'_k) - 2G(X'_k)F_k(X'_k)}{g_k(X'_k)} \right]; $$

and

$$ \mathbb{E}[E_k^2] = \mathbb{E}\left[\frac{G^2(X'_k) + F_k(X'_k) - 2G(X'_k)F_k(X'_k)}{g_k^2(X'_k)}\right] \le \kappa_2, $$

with

$$ \kappa_2 := \sup_{k \in \{1, \dots, d\}} \sup_{x'_k \in \Omega_k} \mathbb{E} \left[ \frac{G^2(x'_k) + F_k(x'_k) - 2G(x'_k)F_k(x'_k)}{g_k^2(x'_k)} \right]. $$

Using such inequalities, the following results are straightforward keeping in mind Corollaries 5–7.

**Corollary 8.** Let $r^* = d_0 - 1$ and $L = d_0 + 1$. Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, \max(d, d_0 + 2)\}$; $f \in \mathcal{L}_{d_0, 0}$ and Assumption 1 and 2 hold. If $h_k = h \propto N^{-\eta}$ with $\eta \in ]\frac{1}{4}, 1[$ and $\zeta^2 \le (dM_{d_0+2}\Gamma_{d_0+2}(\frac{d}{d_0})\kappa_1^{d_0})^{-1}$, then we have
$$ \mathbb{E}\left[\left(\widehat{f}_{N,d_0}^c - f^c\right)^2\right] \le \frac{\left(F_{d_0}^{\max}\right)^2}{N} \left[2(\kappa_2 d(d+0.8))^{d_0} \mathbb{I}_{d_0 \le d_0^*} + 2^d (\kappa_2(d+0.8))^{d_0} \mathbb{I}_{d_0 > d_0^*}\right], \\ \text{provided that } \kappa_2(d+0.8) \ge 1. $$

**Corollary 9.** Let $r^* = d-1$ and $L = d+1$. Assume $f \in \mathcal{H}_\alpha$ with $\alpha \in \{0, d+1\}$; $\zeta \le [dM_{d+1}\Gamma_{d+1}\kappa_1^d]^{-1}$; $h_k = h \propto N^{-\eta}$ with $\eta \in ]\frac{1}{2}, 1[$; and Assumption 1 and 2 hold. Then,

$$ \mathbb{E}\left[\left(\widehat{f}_{N}^{c}-f^{c}\right)^{2}\right] \leq \frac{\left(F_{d}^{\max}\right)^{2}}{N}\left[(\kappa_{2}(d+0.8)+1)^{d}-1\right]. $$

Regarding the choice of $G_k$, from the expression

$$ \mathbb{E}[E_k^2] = \mathbb{E}\left[\frac{F_k(X'_k) + G(X'_k)[G(X'_k) - 2F_k(X'_k)]}{g_k^2(X'_k)}\right], $$

an interesting choice of $G_k$ must satisfy $G_k \le 2F_k$ in order to reduce the value of $\mathbb{E}[E_k^2]$. The following proposition gives interesting choices of $G_k$s. Recall that $X_k \sim F_k$ and $X'_k \sim G_k$ are supported on $\Omega_k, \Omega'_k$, respectively, with $\Omega_k \subseteq \Omega'_k$.

**Proposition 2.** Consider a PDF, $\rho'_k$, supported on $\Omega'_k \setminus \Omega_k$ and $\tau \in ]0, 1]$. The distribution, $G_k$, defined as a mixture of $\rho_k$ and $\rho'_k$, that is,

$$ g_k = \tau \rho_k \mathbb{I}_{\Omega_k} + (1 - \tau) \rho'_k \mathbb{I}_{\Omega'_k \setminus \Omega_k}, $$

allows for reducing $\mathbb{E}[E_k^2]$, provided that $\sup\{x: x \in \Omega\} \le y_0, \forall y_0 \in \Omega'_k \setminus \Omega_k$.

**Proof.** We can check that $G_k \le F_k$. $\square$

In what follows, we are going to consider $\tau = 1$ (i.e., $G_k = F_k$) and $\tau < 1$ with $\rho'_k$ being the uniform distribution.

## 4. Illustrations

In this section, we deploy our derivative-free emulators to approximate analytical functions. For the setting of different parameters needed, we rely on the results of Corollary 5. Indeed, we use $V_k \sim U(-\zeta, \zeta)$ with $\zeta = \left(d\binom{d}{d_0}\left(\frac{1}{2\rho_{\min}}\right)^{d_0}\right)^{-1/2}$; $L = d_0 + 1$ for the identified $d_0$. For each function, we use NL runs of the model to construct our emulators, corresponding to $f(X_i + \beta_l h V_i)$ with the following:

* $i = 1, \dots, N = 500$;

* $\beta_l \in \{0, \pm 2^{k-1}: k = 1, \dots, \left\lfloor \frac{L-1}{2} \right\rfloor\}$ if $L$ is odd and $\beta_l \in \{\pm 2^k: k = 1, \dots, \left\lfloor \frac{L}{2} \right\rfloor\}$; otherwise, $h = N^{-1}$.

Then, we predict the output for $N = 500$ sample points, that is, $f(X_i)$. Finally, sample values are generated using Sobol's sequence, and we compare the predictions associated with the initial distributions, that is, $G_k = F_k$, with a mixture distribution of $F_k$ (i.e., $\tau = 0.9$) (see Proposition 2).

### 4.1. Test Functions

#### 4.1.1. Ishigami's Function ($d = 3$)

The Ishigami function is given by

$$ f(\mathbf{x}) = \sin(x_1) + 7 \sin^2(x_2) + 0.1 x_3^4 \sin(x_1), $$
with $X_j \sim U(-\pi, \pi)$, $j = 1, 2, 3$. The sensitivity indices are $S_1 = 0.3139$, $S_2 = 0.4424$, $S_3 = 0.0$, $S_{T_1} = 0.567$, $S_{T_2} = 0.442$, and $S_{T_3} = 0.243$. Thus, we have $d_0 = 2$ because $\sum_{j=1}^3 (S_j + S_{T_j}) = 2.01$ (see Proposition 1). Moreover, we are going to remove (in our emulator) the main effect term corresponding to $X_3$, as $S_3 = 0$. Figure 1 depicts the predictions versus observations (i.e., simulated outputs) for 500 sample points. We can see that our predictions are in line with the simulated values of the output for both distributions used.

<figure><img src="image_6.png" /><figcaption>Figure 1. Predictions versus simulated outputs (observations) for Ishigami's function.</figcaption></figure>

### 4.1.2. Sobol's g-Function ($d = 10$)

The g-function [49] is defined as follows:

$$f(\mathbf{x}) = \prod_{j=1}^{d=10} \frac{|4x_j - 2| + a_j}{1 + a_j}.$$

with $X_j \sim U(0, 1)$, $j = 1, \dots, d = 10$. This function is differentiable almost everywhere and has different properties according to the values of $\mathbf{a} = (a_j, j = 1, 2, \dots, d)$ [17]. Indeed, the following applies:

* If $\mathbf{a} = [0, 0, 6.52, 6.52, 6.52, 6.52, 6.52, 6.52]^T$ (i.e., type A), we have $S_1 = S_2 = 0.39$, $S_j = 0.0069, \forall j > 2$, $S_{T_1} = S_{T_2} = 0.54$, and $S_{T_j} = 0.013, \forall j > 2$. We have $d_0 = 2$, as $\sum_{j=1}^d (S_j + S_{T_j}) = 2.01$. Moreover, we have $S_1 + S_{T_2} \approx 1$, suggesting that we should only include $X_1$ and $X_2$ in our emulator;

* If $\mathbf{a} = [50, 50, 50, 50, 50, 50, 50, 50]^T$ (i.e., type B), we have $S_j = S_{T_j} = 0.1$ and $\forall j \in \{1, 2, \dots, d\}$, leading to $d_0 = 1$;

* If $\mathbf{a} = [0, 0, 0, 0, 0, 0, 0, 0]^T$ (i.e., type C), we have $S_j = 0.02$ and $S_{T_j} = 0.27, \forall j \in \{1, 2, \dots, d\}$, and we can check that $d_0 = 4$ and that all the inputs are important. Thus, we have to include a lot of ANOVA components in our emulator with small effective effects since the variance of that function is fixed. More information is needed to design the structure of this function better.

Figures 2 and 3 depict the predictions versus the simulated outputs (i.e., observations) of the g-function of type A and type B, respectively, for 500 sample points. Even though we obtain quasi-perfect predictions in the case of the g-function of type B, those of type A face some difficulties in predicting small and null values.
<figure><img src="image_7.png" /><figcaption>Figure 2. Predictions versus simulated outputs (observations) for the g-function of type A.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3. Predictions versus simulated outputs (observations) for the g-function of type B.</figcaption></figure>

## 5. Application: Heat Diffusion Models with Stochastic Initial Conditions

We deployed our emulators to approximate a high-dimensional model defined by the one-dimensional (1-D) diffusion PDE with stochastic initial conditions, that is,

$$ \left\{ \begin{array}{ll} \frac{\partial f}{\partial t} - D \frac{\partial^2 f}{\partial x^2} = 0, & x \in [0, 1], t \in [0, T] \\ f(x, t = 0) = Z(x), & x \in [0, 1], \\ f(x = 0, t) = 0, \quad f(x = 1, t) = 1, & t \in [0, T] \end{array} \right. $$

where $D = 0.0011$ represents the diffusion coefficient. The quantity of interest (QoI) is given by $J(Z(x)) := \frac{1}{2} \int_{0}^{T} \int_{0}^{10} (f(x, t))^2 dxdt$. The spatial discretization consists of subdividing the spatial domain, $[0, 1]$, into $d$ equally-sized cells, leading to $d$ initial conditions or inputs, that is, $Z(x_j)$ with $j = 1, \dots, d$. We assume that the $d = 50$ inputs are independent, and $X_j := Z(x_j) \sim U(\sin(2\pi x_j) - 1.96, \sin(2\pi x_j) + 1.96)$. A time step of 0.025 was considered, starting from 0 up to $T = 5$.
It is known [12] that the exact gradient can be computed as follows: $\nabla_Z J(Z(x)) = f^a(x, 0)$, where $f^a(x, 0)$ stands for the adjoint model of $f$ evaluated at $(x, t = 0)$. Note that only one evaluation of such a function is needed to obtain the gradient of the QOI. The adjoint model is given by (see ref. [12])

$$ \left\{ \begin{array}{ll} -\frac{\partial f^a}{\partial t} - D \frac{\partial^2 f^a}{\partial^2 x} = f, & x \in ]0, 1[, t \in [0, T] \\ f^a(x = 0, t) = f^a(x = 1, t) = 0, & t \in [0, T] \\ f^a(x, T) = 0, & x \in [0, 1] \end{array} \right. $$

The values of the hyper-parameters derived in this paper (considered at the beginning of Section 4) were used to compute the results below. Using the exact values of the gradient (i.e., $f^a(x, 0)$), we computed the main indices ($S_{js}$) and the upper bounds of the total indices ($UB_{js}$) (see Figure 4, top-left panel). It appears that the upper bounds are almost equal to the main indices, showing the absence of interactions. This information is confirmed by the fact that $\sum_{j=1}^{d=50} S_j = 1.09$, leading to $d_0 = 1$. Based on this information, Figure 4 (top-right panel) depicts the predictions versus the observations (simulated outputs) using the derivative-based emulator with all the first-order partial derivatives. In the same sense, Figure 4 (bottom-left panel) depicts the observations versus predictions by including only the ANOVA components for which $UB_j > 0.01$, that is, 37 components. Both results are close together and are in line with the observations. Finally, Figure 4 (bottom-right panel) shows the observations versus predictions for derivative-free emulators using only the components for which $UB_j > 0.01$. It turns out that our emulators provide reliable estimations. As expected (see MSEs), the derivative-based emulator using exact values of derivatives performs better.

<figure><img src="image_9.png" /><figcaption>Figure 4. Main indices (◦) and upper bounds of total indices (+) of $d = 50$ inputs (top-left panel); observations versus predictions using either derivative-based emulators (see the top-right panel when including all components and the bottom-left panel for all other cases) or derivative-free emulators (bottom-right panel).</figcaption></figure>
## 6. Conclusions

In this paper, we have proposed simple, practical, and sound emulators of simulators or estimators of functions using either the available derivatives and distribution functions of inputs or derivative-free methods, such as stochastic approximations. Since our emulators or estimators rely exactly on Db-ANOVA, Sobol's indices and their upper bounds were used to derive the appropriate structures of such emulators so as to reduce epistemic uncertainty. The derivative-based and derivative-free emulators reach the parametric rate of convergence (i.e., O(N⁻¹)) and have dimension-free biases. Moreover, the former emulators enjoy dimension-free MSEs when all cross-partial derivatives are available and, therefore, can cope with higher-dimensional models. However, the MSEs of the derivative-free estimators depend on dimensionality, and we have shown that the stability and accuracy of such emulators require about N ≈ (d + 1)d model runs for full emulators and about N ≈ min(d²d₀, 2d²d₀) runs for unbiased, truncated emulators of order d₀.

To be able to deploy our emulators in practice, we have provided the best known values for the hyper-parameters needed. The numerical results have revealed that our emulators provide efficient predictions of models once the adequate structures of such models are used. While such results are promising, further improvements are going to be investigated in our next works by i) considering distributions of Vjs that may help in reducing the dimensionality in MSEs, ii) taking into account the discrepancy errors by using the output observations (rather than their mean only), and iii) considering local emulators. It is also worth investigating adaptations of such methods in the presence of empirical data.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No empirical datasets are used in this paper. All simulated data are already in the paper.

**Acknowledgments:** We would like to thank the three reviewers for their suggestions and comments that have helped improve our manuscript.

**Conflicts of Interest:** The author has no conflicts of interest to declare regarding the content of this paper.

## Appendix A. Derivations of Unbiased Truncations (Proposition 1)

Keeping in mind the Sobol indices, it is known that $\sum_{v \subset \{1, \dots, d\}} S_v = 1$, which comes down to $\sum_{\substack{v \subset \{1, \dots, d\} \\ |v| \le d_0}} S_v = 1$ for functions of the form $f(\mathbf{X}) = \sum_{\substack{v \subset \{1, \dots, d\} \\ |v| \le d_0}} f_v(\mathbf{X}_v)$. Thus, we have

$$d_0 = \sum_{\substack{v \subset \{1, \dots, d\} \\ |v| \le d_0}} d_0 S_v, \qquad \sum_{j=1}^{d} S_{T_j} = \sum_{\substack{v \subset \{1, \dots, d\} \\ |v| \le d_0}} |v| S_v.$$

Taking the difference yields

$$d_0 - \sum_{j=1}^{d} S_{T_j} = \sum_{\substack{v \subset \{1, \dots, d\} \\ |v| \le d_0 - 1}} (d_0 - |v|)S_v \implies d_0 - \sum_{j=1}^{d} S_{T_j} - (d_0 - 1) \sum_{j=1}^{d} S_j = \sum_{\substack{v \subset \{1, \dots, d\} \\ 2 \le |v| \le d_0 - 1}} (d_0 - |v|)S_v \ge 0$$

which implies that $\sum_{j=1}^{d} S_{T_j} + (d_0 - 1) \sum_{j=1}^{d} S_j \le d_0$.
Additionally, taking $d_0 - (d_0 - 1) \sum_{j=1}^{d} S_{T_j} = \sum_{v \subset \{1, \dots, d\}} [d_0 - (d_0 - 1)|v|]S_v$ yields
$|v| \le d_0$

$$d_0 - (d_0 - 1) \sum_{j=1}^{d} S_{T_j} - \sum_{j=1}^{d} S_j = \sum_{\substack{v \subset \{1, \dots, d\} \\ 2 \le |v| \le d_0}} [d_0 - (d_0 - 1)|v|]S_v \le 0,$$

which implies $(d_0 - 1) \sum_{j=1}^{d} S_{T_j} + \sum_{j=1}^{d} S_j \ge d_0$.

## Appendix B. Proof of Theorem 1

Denote $\vec{i} := (i_1, \dots, i_d) \in \mathbb{N}^d$, $||\vec{i}||_1 = i_1 + \dots + i_d$, $\vec{i}! := i_1! \dots i_d!$, $\vec{x}^i := x_1^{i_1} \dots x_d^{i_d}$, and $D^{(\vec{i})}f := \partial_1^{i_1} \dots \partial_d^{i_d} f$. The Taylor expansion of $f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V})$ about $\mathbf{x}$ of order $\alpha$ is given by

$$f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{r=0}^{\alpha} \sum_{||\vec{i}||_1=r} \frac{D^{(\vec{i})} f(\mathbf{x})}{\vec{i}!} \beta_{\ell}^{r} (\mathbf{h} \mathbf{V})^{\vec{i}} + O(||\beta_{\ell} \mathbf{h} \mathbf{V}||_1^{\alpha+1}).$$

For any $w \subseteq \{1, \dots, d\}$ with the cardinality $|w|$; using $\mathbb{I}_w(\cdot)$ for the indicator function and $\vec{w} := (\mathbb{I}_w(1), \dots, \mathbb{I}_w(d))$ lead to $D^{(\vec{w})}f = D^{|w|}f$. Additionally, using $E_k := \frac{G_k(X_k') - \mathbb{I}_{X_k' \ge x_k}}{G_k(X_k')}$ implies that $R_k = \frac{V_k}{h_k \sigma^2} E_k, k = 1, \dots, d$.

First, by evaluating the above expansion at $X'$ and taking the expectation with respect to $\mathbf{V}$, $A := \sum_{\ell=1}^L C_{\ell}^{(|u|)} \mathbb{E}_{\mathbf{V}} [f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) e_u^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V}))]$ can be written as

$$A = \sum_{r \ge 0} \sum_{||\vec{i}||_1=r} \frac{D^{(\vec{i})} f(\mathbf{X}')}{\vec{i}!} \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^r \right) \mathbb{E}_{\mathbf{V}} \left[ (\mathbf{h} \mathbf{V})^{\vec{i}} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w|=|u|}} \left( \frac{\mathbf{V}}{\mathbf{h} \sigma^2} \mathbf{E} \right)^{\vec{w}} \right] \\ = \sum_{r \ge 0} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ ||\vec{i}||_1=r}} \frac{D^{(\vec{i})} f(\mathbf{X}') \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^r \right)}{\vec{i}! \sigma^{2|u|}} \mathbb{E}_{\mathbf{V}} \left[ (\mathbf{V})^{\vec{i}+\vec{w}} (\mathbf{h})^{\vec{i}-\vec{w}} (\mathbf{E})^{\vec{w}} \right]$$

We can see that $\mathbb{E}_{\mathbf{V}} [(\mathbf{V})^{\vec{i}+\vec{w}} (\mathbf{h})^{\vec{i}-\vec{w}}] \ne 0$ iff $\vec{i} + \vec{w} = 2\vec{q}$, $\forall \vec{q} \in \mathbb{N}^d$. Equation $\vec{i} + \vec{w} = 2\vec{q}$ implies $i_k = 2q_k \ge 0$ if $k \notin w$; otherwise, $i_k = 2q_k - 1 \ge 0$. The last quantity is equivalent to $i_k = 2q_k + 1$ when $k \in w$ with $q_k \in \mathbb{N}$, and it leads to $\vec{i} = 2\vec{q} + \vec{w}$, $\forall \vec{q} \in \mathbb{N}^d$, which also implies that $||\vec{i}||_1 \ge ||\vec{w}||_1 = |u|$.

When $||\vec{q}||_1 = 0$, we have $\vec{i} = \vec{w}$; $D^{(\vec{i})}f = D^w f$ and $\mathbb{E}[(\mathbf{V})^{\vec{i}+\vec{w}}(\mathbf{h})^{\vec{i}-\vec{w}}] = \mathbb{E}[(\mathbf{V})^{2\vec{w}}] = \sigma^{2|u|}$, by independence. Thus, we have

$$A = \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w|=|u|}} D^{|w|} f(\mathbf{X}') \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^{|u|} \right) \prod_{k \in w} E_k + \sum_{\substack{s \ge 1 \\ ||\vec{q}||_1=s}} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w|=|u|}} \frac{D^{(2\vec{q}+\vec{w})} f(\mathbf{X}') \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^{|u|+2s} \right)}{(2\vec{q}+\vec{w})! \sigma^{2|u|}} \mathbb{E} \left[ (\mathbf{V})^{2(\vec{q}+\vec{w})} (\mathbf{h})^{2\vec{q}} (\mathbf{E})^{\vec{w}} \right],$$

using the change of variable $2s = r - |u|$. At this point, the setting $L = 1$, $\beta_{\ell} = 1$ and $C_{\ell}^{(|u|)} = 1$ gives the approximation of $\sum_{w \subseteq \{1, \dots, d\}} D^{|w|} f(\mathbf{X}') \prod_{k \in w} \frac{G_k(X_k') - \mathbb{I}_{X_k' \ge x_k}}{G_k(X_k')}$ of order $\mathcal{O}(|h|^2)$ for all $u \subseteq \{1, \dots, d\}$.
Second, taking the expectation with respect to $\mathbf{X}'$ and the sum over $|u| = 1, \dots, d$, we obtain the result, that is,

$$ \sum_{|u|=1}^{d} \mathbb{E}_{\mathbf{X}'}[A] = \sum_{|u|=1}^{d} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w|=|u|}} \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} \frac{G_k(X'_k) - \mathbb{I}_{X'_k \ge x_k}}{g_k(X'_k)} \right] + \mathcal{O}\left(\|h\|_2^2\right), $$

bearing in mind Equation (1).

Third, we can increase this order up to $\mathcal{O}(\|h\|_2^{2L})$ by using the constraints $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{2s+|u|} = \delta_{0,s}$ $s = 0, 1, \dots, L$ to eliminate some higher-order terms. Thus, the order $\mathcal{O}(\|h\|_2^{2L})$ is reached. Other constraints can lead to the order $\mathcal{O}(\|h\|_2^{2\alpha|u|})$ with $\alpha|u| = 1, \dots, L$.

## Appendix C. Proof of Lemma 1

Recall that $|E_k| = \left| \frac{G_k(X'_k) - \mathbb{I}_{X'_k \ge x_k}}{g_k(X'_k)} \right|$. By using the definition of the absolute value and the fact that $0 \le G(x) \le 1$, we can check that $\mathbb{E}|E_k| \le \frac{1}{2\rho_{\min}}$. Additionally, using the following inequality (see Lemma 1 in ref. [10])

$$ M'_{|w|} \le 2\gamma_0(M'_0)^{1-|w|/d} (M'_d)^{|w|/d} = 2\gamma_0 M'_0 \left(\frac{M'_d}{M'_0}\right)^{|w|/d}, $$

for a given $\gamma_0$, we can write

$$
\begin{align*}
D_{d_0, \rho_{\min}} &:= \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} M'_{|w|} \left(\frac{1}{2\rho_{\min}}\right)^{|w|} \le 2\gamma_0 M'_0 \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d}\right]^{|w|} \\
&= 2\gamma_0 M'_0 \left\{ \sum_{w \subseteq \{1, \dots, d\}} \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d}\right]^{|w|} - \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| \le d_0}} \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d}\right]^{|w|} \right\} \\
&\le 2\gamma_0 M'_0 \left\{ \sum_{\ell=0}^{d} \binom{d}{\ell} \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d}\right]^\ell - \sum_{\ell=0}^{d_0} \binom{d_0}{\ell} \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d}\right]^\ell \right\} \\
&= 2\gamma_0 M'_0 \left\{ \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d} + 1\right]^d - \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d} + 1\right]^{d_0} \right\} \\
&= 2\gamma_0 M'_0 \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d} + 1\right]^{d_0} \left\{ \left[\frac{1}{2\rho_{\min}} \left(\frac{M'_d}{M'_0}\right)^{1/d} + 1\right]^{d-d_0} - 1 \right\},
\end{align*}
$$

and the result holds.

## Appendix D. Proof of Theorem 2

Recall that $R_k = \frac{V_k}{h_k \sigma^2} E_k$ with $E_k := \frac{G_k(X'_k) - \mathbb{I}_{X'_k \ge x_k}}{g_k(X'_k)}$, $k = 1, \dots, d$. Using

$$ A_1 := \sum_{p=1}^{d_0} \sum_{\ell=1}^{L} C_{\ell}^{(p)} \mathbb{E} \left\{ \left[ f(\mathbf{X}' + \beta_{\ell} h \mathbf{V}) \right] e_p^{(1:d)} (\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V})) \right\}, $$
the bias $B := A_1 - f^c(\mathbf{x})$ becomes

$$
\begin{align*}
B &= \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \left\{ \sum_{\ell=1}^{L} C_{\ell}^{(|w|)} \mathbb{E} \left[ f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in w} \frac{V_k}{h_k \sigma^2} E_k \right] - \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} E_k \right] \right\} \\
&\quad - \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} E_k \right] \\
&= \sum_{\substack{w \subseteq \{1, \dots, d_0\} \\ |w| > 0}} \mathbb{E}_{\mathbf{X}'} \left\{ \left( \prod_{k \in w} E_k \right) \left( \sum_{\ell=1}^{L} C_{\ell}^{(|w|)} \mathbb{E}_V \left[ f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in w} \frac{V_k}{h_k \sigma^2} \right] - \mathcal{D}^{|w|} f(\mathbf{X}') \right) \right\} \\
&\quad - \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} \mathbb{E}_{\mathbf{X}'} \left[ \mathcal{D}^{|w|} f(\mathbf{X}') \prod_{k \in w} E_k \right].
\end{align*}
$$

Note that the quantity $\sum_{\ell=1}^L C_{\ell}^{(|w|)} \mathbb{E}_V [f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in w} \frac{V_k}{h_k \sigma^2}] - \mathcal{D}^{|w|} f(\mathbf{X}')$ has been investigated in ref. [13]. To make use of such results in our context given by Equation (5), let $L_w' := \left( \left[ \frac{r^* - |w|}{2} \right] + L - r^* \right) \mathbb{I}_{|w| \le r^*} + (L - r^* - 1) \mathbb{I}_{|w| > r^*}$. Thus, we have

$$
\left| \mathcal{D}^{|w|} f(\mathbf{x}) - \sum_{\ell=1}^{L} C_{\ell}^{(|w|)} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in w} \frac{V_k}{h_k \sigma^2} \right] \right| \le \sigma^{2L_w'} M_{|w|+2L_w'} K_{1,L_w'} \|\mathbf{h}^2\|_2^{L_w'} . \quad (A1)
$$

When $V_k \sim \mathcal{U}(-\zeta, \zeta)$ with $\zeta > 0$ and $k = 1, \dots, d$, then

$$
\left| \mathcal{D}^{|w|} f(\mathbf{x}) - \sum_{\ell=1}^{L} C_{\ell}^{(|w|)} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in w} \frac{V_k}{h_k \sigma^2} \right] \right| \leq M_{|w|+2L'_w} \zeta^{2L'_w} \|\mathbf{h}^2\|_1^{L'_w} \Gamma_{|w|+2L'_w} . \quad (A2)
$$

Using Equation (A1); $g_k > \rho_{\min}$ and the fact that $\mathbb{E}|E_k| \le \frac{1}{2\rho_{\min}}$, we can write

$$
\begin{align*}
|B| &\le \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \sigma^{2L'_w} M_{|w|+2L'_w} K_{1,L'_w} \|\mathbf{h}^2\|_2^{L'_w} \prod_{k \in w} \mathbb{E}[|E_k|] + \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} M'_{|w|} \prod_{k \in w} \mathbb{E}[|E_k|] \\
&\le \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \sigma^{2L'_w} M_{|w|+2L'_w} K_{1,L'_w} \|\mathbf{h}^2\|_2^{L'_w} \prod_{k \in w} \frac{1}{2\rho_{\min}} + \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w| > d_0}} M'_{|w|} \prod_{k \in w} \frac{1}{2\rho_{\min}},
\end{align*}
$$

where $M'_{|w|} = \left\| \mathcal{D}^{|w|} f \right\|_{\infty}$. The results hold by using Lemma 1.

Appendix E. Proof of Corollary 1

First, keeping in mind Theorem 2, we can see that the smallest value of $L_w'$ is $L - r^* - 1$, which is reached when $|w| > r^*$. Thus, the bias verifies

$$
\begin{align*}
|B| &\le \left\| \mathbf{h}^2 \right\|_2^{L-r^*-1} \sigma^{2(L-r^*-1)} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ r^* < |w| \le d_0}} K_{w,(L-r^*-1)} M_{|w|+2(L-r^*-1)} \left( \frac{1}{2\rho_{\min}} \right)^{|w|} \\
&\quad + D_{d_0,\rho_{\min}} + \mathcal{O}\left( \left\| \mathbf{h}^2 \right\|_2^{L-r^*-1} \right).
\end{align*}
$$

Second, using $K_{1,r^*,d_0}^{\max}$, we can write
$$A_3 := \sum_{\substack{w \subseteq \{1, \dots, d\} \\ r^* < |w| \le d_0}} K_{w,(L-r^*-1)} M_{|w|+2(L-r^*-1)} \left(\frac{1}{2\rho_{\min}}\right)^{|w|} \le K_{1,r^*,d_0}^{\max} \sum_{t=r^*+1}^{d_0} \binom{d}{t} \left(\frac{1}{2\rho_{\min}}\right)^t \\ = 2K_{1,d_0}^{\max} \rho_{\min} \left(\frac{d}{2\rho_{\min}}\right)^{r^*+1} \frac{\left(\frac{d}{2\rho_{\min}}\right)^{d_0-r^*} - 1}{d - 2\rho_{\min}},$$

$$\text{because } \binom{d}{t} \le d^t \text{ and } \sum_{l=r^*+1}^{d_0} \binom{d}{l} \left(\frac{1}{2\rho_{\min}}\right)^l \le \sum_{l=r^*+1}^{d_0} \binom{d}{l} \left(\frac{1}{2\rho_{\min}}\right)^l = \left(\frac{d}{2\rho_{\min}}\right)^{r^*+1} \frac{\left(\frac{d}{2\rho_{\min}}\right)^{d_0-r^*} - 1}{\frac{d}{2\rho_{\min}} - 1}.$$

Finally, if $V_k \sim U(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, the following bias is used to derive the result:

$$|B| \le ||h^2||_1^{L-r^*-1}\xi^{2(L-r^*-1)} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ r^* < |w| \le d_0}} M_{|w|+2(L-r^*-1)} \Gamma_{|w|+2(L-r^*-1)} \left(\frac{1}{2\rho_{\min}}\right)^{|w|} \\ + D_{d_0, \rho_{\min}} + \mathcal{O}(|h^2|_1^{L-r^*-1}).$$

## Appendix F. Proof of Corollary 2

For $r^* < |w| \le d - 1$, we can see that $L'_w = 1$, and the order of approximation in Corollary 1 becomes $\mathcal{O}(|h^2|_2)$ or $\mathcal{O}(|h^2|_1)$ because $|w| + 2 \le d + 1$ and $L = r^* + 2$. When $|w| = d$, the smallest order is obtained thanks to ref. [13], Corollary 2.

## Appendix G. Proof of Theorem 3

For the variance of our emulator, we can write

$$\begin{aligned} \mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] &= \frac{1}{N} \mathbb{V} \left[ \sum_{p=1}^{d_0} \sum_{\ell=1}^{L} C_{\ell}^{(p)} f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) e_p^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V})) \right] \\ &\le \frac{1}{N} \mathbb{E} \left[ \left\{ \sum_{p=1}^{d_0} \sum_{\ell=1}^{L} C_{\ell}^{(p)} f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) e_p^{(1:d)}(\mathbf{R}(\mathbf{x}, \mathbf{X}', \mathbf{V})) \right\}^2 \right] \\ &\le \frac{1}{N} \mathbb{E} \left[ \left\{ \sum_{p=1}^{d_0} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ |w|=p}} \left( \prod_{k \in w} \frac{V_k}{h_k \sigma^2} E_k \right) \left( \sum_{\ell=1}^{L} C_{\ell}^{(p=|w|)} f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) \right) \right\}^2 \right]. \end{aligned}$$

Using $A_2 := \sum_{\ell=1}^L C_{\ell}^{(|w|)} f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V})$, and keeping in mind Equation (5), we can write

$$|A_2| = \left| \sum_{\ell=1}^{L} C_{\ell}^{(|w|)} \left[ f(\mathbf{X}' + \beta_{\ell} \mathbf{h} \mathbf{V}) - \sum_{r=0}^{\min(r^*, |w|-1)} \sum_{||\tilde{\mathbf{t}}||_1=\tilde{r}} \frac{\mathcal{D}(\tilde{\mathbf{t}}) f(\mathbf{X}')}{\tilde{\mathbf{t}}!} \beta_{\ell}^r (\mathbf{h} \mathbf{V})^{\tilde{\mathbf{t}}} \right] \right| \\ \le M_{\min(r^*, |w|-1)+1} ||\mathbf{h} \mathbf{V}||_2^{\min(r^*, |w|-1)+1} \sum_{\ell=1}^{L} |C_{\ell}^{(|w|)} \beta_{\ell}^{\min(r^*, |w|-1)+1}|,$$

because $f \in \mathcal{H}_{|w|+1}$. Keeping in mind that $\Gamma_{\min(r^*,|w|-1)+1} = \sum_{\ell=1}^L |C_{\ell}^{(|w|)} \beta_{\ell}^{\min(r^*,|w|-1)+1}|$ and using $\Gamma_{r^*,|w|} := M_{\min(r^*,|w|-1)+1} \Gamma_{\min(r^*,|w|-1)+1}$, the variance becomes
$$
\begin{aligned}
\mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] &\le \frac{1}{N} \mathbb{E} \left[ \left\{ \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \mathcal{F}_{r^*,|w|} ||\mathbf{h}\mathbf{V}||_2^{\min(r^*,|w|-1)+1} \left( \prod_{k \in w} \frac{V_k}{h_k \sigma^2} E_k \right) \right\}^2 \right] \\
&\le \frac{(\mathcal{F}_{r^*,d_0}^{\max})^2}{N} \mathbb{E} \left[ \left\{ \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \prod_{k \in w} \left( \frac{V_k E_k}{h_k \sigma^2} ||\mathbf{h}\mathbf{V}||_2^{\min(r^*,|w|-1)+1} \frac{|w|}{|w|} \right) \right\}^2 \right] \\
&= \frac{(\mathcal{F}_{r^*,d_0}^{\max})^2}{N} \mathbb{E} \left[ \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \prod_{k \in w} \left( \frac{V_k E_k}{h_k \sigma^2} ||\mathbf{h}\mathbf{V}||_2^{\frac{\min(r^*,|w|-1)+1}{|w|}} \right)^2 \right],
\end{aligned}
$$

$$
\begin{aligned}
&\text{because } \mathbb{E}[V_k] = 0, \mathbb{E}\left[V_k ||\mathbf{h}\mathbf{V}||_2^{\frac{\min(r^*,|w|-1)+1}{|w|}}\right] = 0 \text{ and when } w \neq w', \\
&\mathbb{E}\left[\prod_{k \in w} \left(\frac{V_k E_k}{h_k \sigma^2} ||\mathbf{h}\mathbf{V}||_2^{\frac{\min(r^*,|w|-1)+1}{|w|}}\right) \prod_{\ell \in w'} \left(\frac{V_\ell E_\ell}{h_\ell \sigma^2} ||\mathbf{h}\mathbf{V}||_2^{\frac{\min(r^*,|w|-1)+1}{|w|}}\right)\right] = 0.
\end{aligned}
$$

For the second result, by expanding $E_k^2$ and knowing that $0 \le G(x) \le 1$, we can see that $\mathbb{E}[E_k^2] \le \frac{1}{3\rho_{\min}^2}$. Additionally, using $Z_k = V_k/\sigma$ and the fact that $h_k = h$, $r^* = d_0 - 1$ and $|w| \le d_0$, we have

$$
\begin{aligned}
\mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] &\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \prod_{k \in w} \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{V}||_2^2]}{\sigma^4} \mathbb{E}[E_k^2] \right) \\
&\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \prod_{k \in w} \left( \mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] \mathbb{E}[E_k^2] \right) \tag{A3} \\
&\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \sum_{\substack{w \subseteq \{1, \dots, d\} \\ 0 < |w| \le d_0}} \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{|w|} = \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \sum_{p=1}^{d_0} \binom{d}{p} \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^p \\
&\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \frac{d \mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] - 3\rho_{\min}^2} \left[ \left( \frac{d \mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} - 1 \right].
\end{aligned}
$$

Finally, when $\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2] \ge 3\rho_{\min}^2$ and knowing that $\sum_{p=1}^{d_0} \binom{d}{p} \le \sum_{p=1}^{d_0} d^p = \frac{d^{d_0-1}}{d-1} d \le 2d^{d_0}$, we can write

$$
\begin{aligned}
\mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] &\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \sum_{p=1}^{d_0} \binom{d}{p} \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^p \\
&\le \frac{(\mathcal{F}_{d_0}^{\max})^2}{N} \left( \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0} \sum_{p=1}^{d_0} \binom{d}{p} = \frac{2(\mathcal{F}_{d_0}^{\max})^2}{N} \left( d \frac{\mathbb{E}[Z_1^2 ||\mathbf{Z}||_2^2]}{3\rho_{\min}^2} \right)^{d_0}.
\end{aligned}
$$

Additionally, note that $\sum_{p=1}^{d_0} \binom{d}{p} \le 2^d$.
## Appendix H. Proof of Corollary 3

Since $r^* = d_0 - 1$ and $|w| \le d_0$, we have $K_{1,\rho_{\min},d_0-1} = \binom{d}{d_0} \left(\frac{1}{2\rho_{\min}}\right)^{d_0}$. The first result is obvious using the variance of the emulator provided in Theorem 3 and the bias from Corollary 1.

The second result is due to the fact that when $G_j = F_j$, the terms in the Db expansion of $f$ are $L_2$-orthogonal.

## Appendix I. Proof of Corollary 5

As $V_k \sim U(-\zeta, \bar{\zeta})$, then $Z_k \sim U(-\sqrt{3}, \sqrt{3})$ and $\mathbb{E}[Z_1^2||\mathbf{Z}||_2^2] = d + 4/5$. Thus, the first result holds.

For the second result, taking $\rho_{\min} = \sqrt{\frac{c_0d(d+0.8)}{3}}$ with $1 < c_0 \le 2$ yields

$$ \mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] \le \bar{\Sigma}_{d_0} = \frac{(F_{d_0}^{\max})^2}{N} \left( \frac{1 - (1/c_0)^{d_0}}{c_0 - 1} \right) \le \frac{(F_{d_0}^{\max})^2}{N(c_0 - 1)}. $$

## Appendix J. Proof of Corollary 6

Using Equation (A3) and the fact that $r^* = d - 1$, $d_0 = d$, we can check that

$$ \mathbb{V}[\widehat{f}_{N,d_0}^c(\mathbf{x})] \le \frac{(F_{d_0}^{\max})^2}{N} \left[ \left( \frac{\mathbb{E}[Z_1^2 || \mathbf{Z} ||_2^2]}{3\rho_{\min}^2} + 1 \right)^d - 1 \right] \le \bar{\Sigma}_d, $$

with

$$ \bar{\Sigma}_d := \frac{d(F_{d_0}^{\max})^2 \mathbb{E}[Z_1^2 || \mathbf{Z} ||_2^2]}{N(d\mathbb{E}[Z_1^2 || \mathbf{Z} ||_2^2] - 3\rho_{\min}^2)} \left[ \left( \frac{d\mathbb{E}[Z_1^2 || \mathbf{Z} ||_2^2]}{3\rho_{\min}^2} \right)^d - 1 \right]. $$

Thus, the results hold using Corollaries 2 and 3.

## References

1. Max D. Morris, T.J.M.; Ylvisaker, D. Bayesian Design and Analysis of Computer Experiments: Use of Derivatives in Surface Prediction. Technometrics 1993, 35, 243–255.

2. Solak, E.; Murray-Smith, R.; Leithead, W.; Leith, D.; Rasmussen, C. Derivative observations in Gaussian process models of dynamic systems. Adv. Neural Inf. Process. Syst. 2002, 15, 1–8.

3. Le Dimet, F.X.; Talagrand, O. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus A Dyn. Meteorol. Oceanogr. 1986, 38, 97–110. [CrossRef]

4. Le Dimet, F.X.; Ngodock, H.E.; Luong, B.; Verron, J. Sensitivity analysis in variational data assimilation. J. Meteorol. Soc. Jpn. 1997, 75, 245–255. [CrossRef]

5. Cacuci, D.G. Sensitivity and Uncertainty Analysis—Theory; Chapman & Hall/CRC: Boca Raton, FL, USA, 2005.

6. Gunzburger, M.D. Perspectives in Flow Control and Optimization; SIAM: Philadelphia, PA, USA, 2003.

7. Borzi, A.; Schulz, V. Computational Optimization of Systems Governed by Partial Differential Equations; SIAM: Philadelphia, PA, USA, 2012.

8. Ghanem, R.; Higdon, D.; Owhadi, H. Handbook of Uncertainty Quantification; Springer International Publishing: Berlin/Heidelberg, Germany, 2017.

9. Agarwal, A.; Dekel, O.; Xiao, L. Optimal Algorithms for Online Convex Optimization with Multi-Point Bandit Feedback. In Proceedings of the COLT, 27–29 June 2010, Haifa, Israel; Citeseer: Princeton, NJ, USA, 2010; pp. 28–40.

10. Bach, F.; Perchet, V. Highly-Smooth Zero-th Order Online Optimization. In Proceedings of the 29th Annual Conference on Learning Theory, New York, NY, USA, 23–26 June 2016; Volume 49, pp. 257–283.

11. Akhavan, A.; Pontil, M.; Tsybakov, A.B. Exploiting higher order smoothness in derivative-free optimization and continuous bandits. In Proceedings of the NIPS '20: Proceedings of the 34th International Conference on Neural Information Processing Systems, Vancouver, BC, Canada, 6–12 December 2020.
12. Lamboni, M. Optimal and Efficient Approximations of Gradients of Functions with Nonindependent Variables. *Axioms* **2024**, *13*, 426. [CrossRef]

13. Lamboni, M. Optimal Estimators of Cross-Partial Derivatives and Surrogates of Functions. *Stats* **2024**, *7*, 1–22. [CrossRef]

14. Chkifa, A.; Cohen, A.; DeVore, R.; Schwab, C. Sparse adaptive Taylor approximation algorithms for parametric and stochastic elliptic PDEs. *ESAIM Math. Model. Numer. Anal.* **2013**, *47*, 253–280. [CrossRef]

15. Patil, P.; Babaee, H. Reduced-Order Modeling with Time-Dependent Bases for PDEs with Stochastic Boundary Conditions. *SIAM/ASA J. Uncertain. Quantif.* **2023**, *11*, 727–756. [CrossRef]

16. Sobol, I.M.; Kucherenko, S. Derivative based global sensitivity measures and the link with global sensitivity indices. *Math. Comput. Simul.* **2009**, *79*, 3009–3017. [CrossRef]

17. Kucherenko, S.; Rodriguez-Fernandez, M.; Pantelides, C.; Shah, N. Monte Carlo evaluation of derivative-based global sensitivity measures. *Reliab. Eng. Syst. Saf.* **2009**, *94*, 1135–1148. [CrossRef]

18. Lamboni, M.; Iooss, B.; Popelin, A.L.; Gamboa, F. Derivative-based global sensitivity measures: General links with Sobol' indices and numerical tests. *Math. Comput. Simul.* **2013**, *87*, 45–54. [CrossRef]

19. Roustant, O.; Fruth, J.; Iooss, B.; Kuhnt, S. Crossed-derivative based sensitivity measures for interaction screening. *Math. Comput. Simul.* **2014**, *105*, 105–118. [CrossRef]

20. Roustant, O.; Barthe, F.; Iooss, B. Poincaré inequalities on intervals—Application to sensitivity analysis. *Electron. J. Statist.* **2017**, *11*, 3081–3119. [CrossRef]

21. Lamboni, M. Derivative-based integral equalities and inequality: A proxy-measure for sensitivity analysis. *Math. Comput. Simul.* **2021**, *179*, 137–161. [CrossRef]

22. Lamboni, M. Weak derivative-based expansion of functions: ANOVA and some inequalities. *Math. Comput. Simul.* **2022**, *194*, 691–718. [CrossRef]

23. Lamboni, M.; Kucherenko, S. Multivariate sensitivity analysis and derivative-based global sensitivity measures with dependent variables. *Reliab. Eng. Syst. Saf.* **2021**, *212*, 107519. [CrossRef]

24. Krige, D.G. A Statistical Approaches to Some Basic Mine Valuation Problems on the Witwatersrand. *J. Chem. Metall. Soc. S. Afr.* **1951**, *52*, 119–139.

25. Currin, C.; Mitchell, T.; Morris, M.; Ylvisaker, D. Bayesian Prediction of Deterministic Functions, with Applications to the Design and Analysis of Computer Experiments. *J. Am. Stat. Assoc.* **1991**, *86*, 953–963. [CrossRef]

26. Oakley, J.E.; O'Hagan, A. Probabilistic sensitivity analysis of complex models: A Bayesian approach. *J. R. Stat. Soc. Ser. B (Stat. Methodol.)* **2004**, *66*, 751–769. [CrossRef]

27. Conti, S.; O'Hagan, A. Bayesian emulation of complex multi-output and dynamic computer models. *J. Stat. Plan. Inference* **2010**, *140*, 640–651. [CrossRef]

28. Kennedy, M.C.; O'Hagan, A. Bayesian calibration of computer models. *J. R. Stat. Soc. Ser. B (Stat. Methodol.)* **2001**, *63*, 425–464. [CrossRef]

29. Xiu, D.; Karniadakis, G. The Wiener-Askey polynomial chaos for stochastic differential equations. *Siam J. Sci. Comput.* **2002**, *24*, [CrossRef]

30. Ghanem, R.G.; Spanos, P.D. *Stochastic Finite Elements: A Spectral Approach*; Springer: New York, NY, USA, 1991; pp. 1–214.

31. Sudret, B. Global sensitivity analysis using polynomial chaos expansions. *Reliab. Eng. Syst. Saf.* **2008**, *93*, 964–979. [CrossRef]

32. Wahba, G. *Spline Models for Observational Data*; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 1990. [CrossRef]

33. Wong, R.K.W.; Storlie, C.B.; Lee, T.C.M. A Frequentist Approach to Computer Model Calibration. *J. R. Stat. Soc. Ser. B Stat. Methodol.* **2016**, *79*, 635–648. [CrossRef]

34. Friedman, J.H.; Popescu, B.E. Predictive learning via rule ensembles. *Ann. Appl. Stat.* **2008**, *2*, 916–954. [CrossRef]

35. Horiguchi, A.; Pratola, M.T. Estimating Shapley Effects in Big-Data Emulation and Regression Settings using Bayesian Additive Regression Trees. *arXiv* **2024**, arXiv:2304.03809.

36. Migliorati, G.; Nobile, F.; Schwerin, E.; Tempone, R. Analysis of Discrete L2 Projection on Polynomial Spaces with Random Evaluations. *Found. Comput. Math.* **2014**, *14*, 419–456.

37. Hampton, J.; Doostan, A. Coherence motivated sampling and convergence analysis of least squares polynomial Chaos regression. *Comput. Methods Appl. Mech. Eng.* **2015**, *290*, 73–97. [CrossRef]

38. Cohen, A.; Davenport, M.A.; Leviatan, D. On the stability and accuracy of least squares approximations. *arXiv* **2018**, arXiv:math.NA/1111.4422. [CrossRef]

39. Tsybakov, A. *Introduction to Nonparametric Estimation*; Springer: New York, NY, USA, 2009.

40. Zemanian, A. *Distribution Theory and Transform Analysis: An Introduction to Generalized Functions, with Applications*; Dover Books on Advanced Mathematics; Dover Publications: Garden City, NY, USA, 1987.

41. Strichartz, R. *A Guide to Distribution Theory and Fourier Transforms*; Studies in Advanced Mathematics; CRC Press: Boca Raton, FL, USA, 1994.
42. Caflisch, R.E.; Morokoff, W.J.; Owen, A.B. Valuation of mortgage-backed securities using Brownian bridges to reduce effective dimension. *J. Comput. Financ.* **1997**, *1*, 27–46. [CrossRef]

43. Owen, A. The dimension distribution and quadrature test functions. *Stat. Sin.* **2003**, *13*, 1–17.

44. Rabitz, H. General foundations of high dimensional model representations. *J. Math. Chem.* **1999**, *25*, 197–233. [CrossRef]

45. Kuo, F.; Sloan, I.; Wasilkowski, G.; Woźniakowski, H. On decompositions of multivariate functions. *Math. Comput.* **2010**, *79*, 953–966. [CrossRef]

46. Alatawi, M.S.; Martinucci, B. On the Elementary Symmetric Polynomials and the Zeros of Legendre Polynomials. *J. Math. 2022*, 1–9. [CrossRef]

47. Arafat, A.; El-Mikkawy, M. A Fast Novel Recursive Algorithm for Computing the Inverse of a Generalized Vandermonde Matrix. *Axioms* **2023**, *12*, 27. [CrossRef]

48. Rawashdeh, E. A simple method for finding the inverse matrix of Vandermonde matrix. *Mathematiqki Vesnik* **2019**, *71*, 207–213.

49. Homma, T.; Saltelli, A. Importance measures in global sensitivity analysis of nonlinear models. *Reliab. Eng. Syst. Saf.* **1996**, *52*, 1–17. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.