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

# Optimal Estimators of Cross-Partial Derivatives and Surrogates of Functions

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

**Abstract:** Computing cross-partial derivatives using fewer model runs is relevant in modeling, such as stochastic approximation, derivative-based ANOVA, exploring complex models, and active subspaces. This paper introduces surrogates of all the cross-partial derivatives of functions by evaluating such functions at *N* randomized points and using a set of *L* constraints. Randomized points rely on independent, central, and symmetric variables. The associated estimators, based on *NL* model runs, reach the optimal rates of convergence (i.e., *O*(N⁻<fnref n="1" />)), and the biases of our approximations do not suffer from the curse of dimensionality for a wide class of functions. Such results are used for (i) computing the main and upper bounds of sensitivity indices, and (ii) deriving emulators of simulators or surrogates of functions thanks to the derivative-based ANOVA. Simulations are presented to show the accuracy of our emulators and estimators of sensitivity indices. The plug-in estimates of indices using the U-statistics of one sample are numerically much stable.

**Keywords:** derivative-based ANOVA; high-dimensional models; independent input variables; optimal estimators of derivatives; sensitivity analysis

**MSC:** 62Fxx; 62J10; 49-XX; 26D10

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

### ✓ check for updates

**Citation:** Lamboni, M. Optimal Estimators of Cross-Partial Derivatives and Surrogates of Functions. *Stats* **2024**, *7*, 697–718. https://doi.org/10.3390/stats7030042

Academic Editor: Wei Zhu

Received: 31 May 2024

Revised: 4 July 2024

Accepted: 12 July 2024

Published: 14 July 2024

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

**Copyright:** © 2024 by the author. 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 relevant in modeling, such as inverse problems, first-order and second-order stochastic approximation methods [1–4], exploring complex mathematical models or simulators, derivative-based ANOVA (Db-ANOVA) or exact expansions of functions, and active subspaces. First-order derivatives or gradients are sometime available in modeling. Instances are (i) models defined via their rates of change with respect to their inputs; (ii) implicit functions defined via their derivatives [5,6]; (iii) cases listed in [7,8] and the references therein.

In FANOVA and sensitivity analysis user and developer communities (see, e.g., [9–14]), screening of input variables and interactions of high-dimensional simulators is often performed before building emulators of such models using Gaussian processes [15–19], polynomial chaos expansions and SS-ANOVA [20,21], or other machine learning approaches. Emulators are fast-evaluated models that better approximate complex and/or too-expansive simulators. Efficient variance-based screening methods rely on the upper bounds of generalized sensitivity indices, including Sobol' indices (see [14,22–26] for independent inputs and [27–29] for non-independent variables). Such upper bounds require the computations of cross-partial derivatives, even for simulators for which these computations are time-demanding or impossible. Also, active subspaces rely on the first-order derivatives for performing dimension reduction and then for approximating complex models [30–33].

For functions with full interactions, all the cross-partial derivatives are used in the integral representations of the infinitesimal increment of functions [34], and in the unanchored decompositions of functions in the Sobolev space [35]. Recently, such derivatives have become crucial in the Db-ANOVA representation of every smooth function, such as

<footnote n="1">¹ Department DFR-ST, University of Guyane, 97346 Cayenne, France; matieyendou.lamboni@gmail.com</footnote>

<footnote n="2">² 228-UMR Espace-Dev, University of Guyane, University of Réunion, IRD, University of Montpellier, 34090 Montpellier, France</footnote>
high-dimensional PDE models. Indeed, it is known, in [14], that every smooth function $f$ admits an exact Db-ANOVA decomposition, that is, $\forall \mathbf{x} \in \Omega \subseteq \mathbb{R}^d$,

$$f(\mathbf{x}) = \mathbb{E}[f(\mathbf{X}')] + \sum_{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 $\mathcal{D}^{|v|} f$ stands for the cross-partial derivative with respect to $x_k$ for any $k \in v$; $\mathbf{X}' := (X'_1, \dots, X'_d)$ is a random vector of independent variables, supported on an open $\Omega$ with margins $G_j$s and densities $g_j$s (i.e., $X'_j \sim G_j$, $j = 1, \dots, d$).

Computing all the cross-partial derivatives using a few model runs or evaluations of functions is challenging. Direct computations of accurate cross-partial derivatives were considered in [36,37] using the generalization of Richardson's extrapolation. Such approximations of all the cross-partial derivatives require a number of model runs that strongly depends on the dimensionality (i.e., $d$). While adjoint-based methods can provide the gradients for some POE/PDE-based models using only one simulation of the adjoint models [38-43], note that computing the Hessian matrix requires running $d$ second-order adjoint models in general, provided that such models are available (see [38,44]).

Stochastic perturbations methods or Monte Carlo approaches have been used in stochastic approximations (see, e.g., [1-4]), including derivative-free optimization (see [2-4,45-48] and the references therein), for computing the gradients and Hessian of functions. Such approaches lead to the estimates of gradients using a number of model runs that can be less than the dimensionality [48,49]. While gradients computations and the convergence analysis are considered in the first-order stochastic approximations, estimators of the Hessian matrices are investigated in the second-order stochastic approximations [2,50-53]. Most of such approaches rely on the Taylor expansions of functions and randomized kernels and/or a random vector that is uniformly distributed on the unit sphere. Nevertheless, independent variables are used in [48,50,53], and the approaches considered in [53,54] rely on the Stein identity [55]. Note that the upper bounds of the biases of such approximations depend on the dimensionality, except in the work [48] for the gradients only. Moreover, the convergence analysis for more than the second-order cross-partial derivatives are not available according to our knowledge.

Given a smooth function defined on $\mathbb{R}^d$, the motivation of this paper consists of proposing new approaches for deriving surrogates of cross-partial derivatives and derivative-based emulators of functions that

* Are simple to use and generic by making use of $d$ independent variables that are symmetrically distributed about zero and a set of constraints;

* Lead to dimension-free upper bounds of the biases related to the approximations of cross-partial derivatives for a wide class of functions;

* Provide estimators of cross-partial derivatives that reach the optimal and parametric rates of convergence;

* Can be used for computing all the cross-partial derivatives and emulators of functions at given points using a small number of model runs.

In this paper, new expressions of cross-partial derivatives of any order are derived in Section 3 by combining the properties of (i) the generalized Richardson extrapolation approach so as to increase the approximations accuracy, and (ii) the Monte Carlo approaches based only on independent random variables that are symmetrically distributed about zero. Such expressions are followed by their order of approximations and biases. We also derive the estimators of such new expressions and their associated mean squared errors, including the rates of convergence for some classes of functions (see Section 3.3). Section 3.4 provides the derivative-based emulators of functions, depending on the strength of the interactions or (equivalently) the cross-partial derivatives, thanks to Equation (1). The strength of the interactions can be assessed using sensitivity analysis. Thus, Section 4 deals with the derivation of new expressions of sensitivity indices and their estimators by making use of
the proposed surrogates of cross-partial derivatives. Simulations based on test functions are considered in Section 5 to show the accuracy of our approach, and we conclude this work in Section 6.

## 2. Preliminary

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

For a non-empty subset $u \subseteq \{1, \dots, d\}$, we use $|u|$ for its cardinality (i.e., the number of elements in $u$) and $(\sim u) := \{1, \dots, d\} \setminus u$. Also, we use $\mathbf{X}_u := (X_j, \forall j \in u)$ for a subset of inputs, and we have the partition $\mathbf{X} = (\mathbf{X}_u, \mathbf{X}_{\sim u})$. Assume that:

**Assumption 1 (A1).** $\mathbf{X}$ is a random vector of independent variables, supported on $\Omega$.

Working with partial derivatives requires a specific mathematical space. Given an integer $n \in \mathbb{N} \setminus \{0\}$ and an open set $\Omega \subseteq \mathbb{R}^d$, consider a weak partial differentiable function $f: \Omega \to \mathbb{R}^n$ [56,57] and a subset $v \subseteq \{1, \dots, d\}$ with $|v| > 0$. Namely, we use $\mathcal{D}^{|v|}f := (\prod_{k \in v} \frac{\partial^{j_k}}{\partial x_k})f$ for the $|v|$-th weak cross-partial derivatives of each component of $f$ with respect to each $x_k$ with $k \in v$.

Likewise, given $\vec{\tau} := (i_1, \dots, i_d) \in \mathbb{N}^d$, denote $\mathcal{D}^{(\vec{\tau})}f := (\prod_{k=1}^d \frac{\partial^{j_k}}{\partial x_k})f$ and $\mathbb{I}_v(j) = 1$ if $j \in v$ and zero otherwise. Thus, taking $\vec{v} := (\mathbb{I}_v(1), \dots, \mathbb{I}_v(d))$ yields $\mathcal{D}^{|v|}f = \mathcal{D}^{(\vec{v})}f$. Moreover, denote $(\mathbf{x})^{\vec{\tau}} = \mathbf{x}^{\vec{\tau}} := \prod_{k=1}^d x_k^{j_k}, \vec{\tau}! := i_1! \dots i_d!$, and consider the Hölder space of $\alpha$-smooth functions 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_{0 \le i_1 + \dots + i_d \le \alpha - 1} \frac{\mathcal{D}^{(\vec{\tau})} f(\mathbf{y})}{\vec{\tau}!} (\mathbf{x} - \mathbf{y})^{\vec{\tau}} \right| \le M_{\alpha} ||\mathbf{x} - \mathbf{y}||_2^{\alpha} \right\} $$

with $\alpha \ge 1$, $M_{\alpha} > 0$, and $\mathcal{D}^{(\vec{\tau})}f(\mathbf{y})$ as weak cross-partial derivatives. We use $||\cdot||_2$ for the Euclidean norm, $||\cdot||_1$ for the $L_1$-norm, $\mathbb{E}(\cdot)$ for the expectation, and $\mathbb{V}(\cdot)$ for the variance.

## 3. Surrogates of Cross-Partial Derivatives and New Emulators of Functions

### 3.1. New Expressions of Cross-Partial Derivatives

This section aims at providing expressions of cross-partial derivatives using the model of interest, and new independent random vectors. We provide approximated expressions of $\mathcal{D}^{|u|}f(\mathbf{x})$ for all $u \subseteq \{1, \dots, d\}$ and the associated orders of approximations.

Given $L, q \in \mathbb{N} \setminus \{0\}$, consider $\beta_{\ell} \in \mathbb{R}$ with $\ell = 1, \dots, L$, $h := (h_1, \dots, h_d) \in \mathbb{R}_+^d$, and denote with $\mathbf{V} := (V_1, \dots, V_d)$ as $d$-dimensional random vectors 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. $$

Random vectors of $d$ independent variables that are symmetrically distributed about zero are instances of $\mathbf{V}$, including the standard Gaussian random vector and symmetric uniform distributions about zero.

Denote $\beta_{\ell}h\mathbf{V} := (\beta_{\ell}h_1V_1, \dots, \beta_{\ell}h_dV_d)$. The reals $\beta_{\ell}$s are used for controlling the order of derivatives (i.e., $|u|$) we are interested in, while $V_j$s help in selecting one particular derivative of order $|u|$. Finally, $h_j$s aim at defining a neighborhood of a sample point $\mathbf{x}$ of $\mathbf{X}$ that will be used. Thus, using $\beta_{\max} := \max(|\beta_1|, \dots, |\beta_L|)$ and keeping in mind the variance of $\beta_{\ell}h_jV_j$, we assume that $\forall j \in \{1, \dots, d\}$

**Assumption 2 (A2).** $\beta_{\max}h_j\sigma \le 1/2$ or, equivalently, $0 < \beta_{\max}h_j|V_j| \le 1$ for bounded $V_j$s.
Based on the above framework, Theorem 1 provides a new expression of the cross-partial derivatives of $f$. Recall that $|u|$ is the cardinality of $u$ and $D^{|u|}f = D^{(\vec{u})}f$.

**Theorem 1.** Consider distinct $\beta_{\ell}s$, and assume that $f \in \mathcal{H}_{\alpha}$ with $\alpha \ge |u| + 2L$ and (A2) holds. Then, for any $u \subseteq \{1, \dots, d\}$ with $|u| > 0$, there exists $\alpha_{|u|} \in \{1, \dots, L\}$ and coefficients $C_1^{(|u|)}, \dots, C_L^{(|u|)}$ such that

$$ \mathcal{D}^{|u|}f(\mathbf{x}) = \sum_{\ell=1}^{L} C_{\ell}^{(|u|)} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in u} \frac{V_k}{(h_k \sigma^2)} \right] + \mathcal{O} \left( \|\mathbf{h}\|_2^{2\alpha_{|u|}} \right). \quad (2) $$

**Proof.** The detailed proof is provided in Appendix A. □

In view of Theorem 1, we are able to compute all the cross-partial derivatives using the same evaluations of functions with the same or different order of approximations, depending on the constraints imposed to determine the coefficients $C_1^{(|u|)}, \dots, C_L^{(|u|)}$ (see Appendix A). While the setting $L = 1, \beta_1 = 1, C_1^{(|u|)} = 1$ or the constraints $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$; $r = 0, \dots, L-2, |u|, L-1 \le |u|$ lead to the order $\mathcal{O}(\|\mathbf{h}\|_2^2)$, one can increase that order up to $\mathcal{O}(\|\mathbf{h}\|_2^{2L})$ by using either $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{r+|u|} = \delta_{0,r}$; $r = 0, 2, \dots, 2(L-1)$ or the full constraints given by $\sum_{\ell=1}^{L+|u|} C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$; $r = 0, \dots, 2L + |u| - 1$. The last setting is going to improve the approximations and numerical computations of derivatives. Since increasing the number of constraints requires more evaluations of simulators, and in ANOVA-like decomposition of $f(\mathbf{X})$, it is common to neglect the higher-order components or, equivalently, the higher-order cross-partial derivatives thanks to Equation (1), the following parsimony number of constraints may be considered. Given an integer $r^* > 0$, controlling the partial derivatives of order up to $r^* \le L-2$ can be performed using the constraints

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

Equation (3) gives approximations of all the cross-partial derivatives of $\mathcal{O}(\|\mathbf{h}\|_2^{\alpha_{|u|}})$ where $o := \left\{ \begin{array}{ll} L - |u| & \text{if } 1 \le |u| \le r^* \\ L - r^* - 1 & \text{otherwise} \end{array} \right.$ and $\alpha_{|u|} = o$ if $o$ is even and $o+1$ otherwise. This equation relies on the Vandermonde matrices and the generalized Vandermonde matrices, which ensure the existence and uniqueness of the coefficients for distinct values of $\beta_{\ell}s$ (i.e., $\beta_{\ell_1} \ne \beta_{\ell_2}$) because the determinant is of the form $\prod_{1 \le \ell_1 < \ell_2 \le L'} (\beta_{\ell_1} - \beta_{\ell_2})$ (see [58,59] for more details and the inverse of such matrices).

**Remark 1.** When $L = 1$, we must have $\beta_1 = 1, C_1^{(|u|)} = 1, \forall u \subseteq \{1, \dots, d\}$. Thus, the coefficient $C_i^{(|u|)}$ does not necessarily depend on $|u|$. Taking $L$ for an even integer, the following nodes may be considered: $\{\beta_1, \dots, \beta_L\} = \{\pm 2^k, k = 0, \dots, \frac{L-2}{2}\}$. When $L$ is odd, one may add 0 to the above set. Of course, other possibilities can be considered provided that $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{|u|} = 1$.

**Remark 2.** For a given $u_0 \subseteq \{1, \dots, d\}$, if we are only interested in all the $|v|^{th}$ cross-partial derivatives with $v \subseteq u_0$, it is better to set $h_{\sim u} = 0$ in Equation (2).

**Remark 3.** Links to other works.

Consider $\beta_1 = -1, \beta_2 = 1, \beta_3 = 0$, and $V_j \sim \mathcal{N}(0, 1)$ with $j = 1, \dots, d$. Using $L = 1$ or $L = 2$ or $L = 3$, our estimators of the first-order and the second-order cross-partial derivatives are very similar to the results obtained in [53].

Using the uniform perturbations and $L = 2$ and $L = 3$, our estimators of the first-order and the second-order cross-partial derivatives are similar to those provided in [50]. However, we will
see later that specific uniform distributions allow for obtaining dimension-free upper bounds of the biases.

## 3.2. Upper Bounds of Biases

To derive precise biases of our approximations provided in Theorem 1, different structural assumptions on the deterministic functions $f$ and $V$ are considered. Assume that $f \in \mathcal{H}_\alpha$ with $\alpha \ge |u|$ is sufficient to define $\mathcal{D}^{|u|}f(x)$ for any $u \subseteq \{1, \dots, d\}$. Note that such an assumption does not depend on the dimensionality $d$. For the sequel of generality, we provide the upper bounds of the biases for any value of $L$ by considering two sets of constraints.

Denote with $\mathbf{R} := (R_1, \dots, R_d)$ a $d$-dimensional random vector of independent variables that are centered about zero and standardized (i.e., $\mathbb{E}[R_k^2] = 1, k = 1, \dots, d$), and $\mathcal{R}_c$ the set of such random vectors. For any $r \in \mathbb{N}$, define

$$ \Gamma_r := \sum_{\ell=1}^{L} |C_{\ell}^{(|u|)} \beta_{\ell}^{r}|; \qquad K_{1,L} := \inf_{\mathbf{R} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| \mathbf{R}^2 \right\|_2^L \prod_{k \in u} R_k^2 \right] \Gamma_{|u|+2L}. $$

**Corollary 1.** Consider distinct $\beta_\ell$s and the constraints $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{|u|+r} = \delta_{0,r}$ with $r = 0, 2, 4, \dots, 2(L-1)$. If $f \in \mathcal{H}_{|u|+2L}$ and (A2) hold, then there is $M_{|u|+2L} > 0$ such that

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

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

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

**Proof.** See Appendix B for the detailed proof. □

In view of Corollary 1, one obtains the upper bounds that do not depend on the dimensionality $d$ by choosing $\zeta = \sigma = \frac{1}{\sqrt{d}}$, $h_k = h$ for instance. When $\Gamma_{|u|+2L} > 1$, the choice $\zeta = \sigma = d^{-\frac{1}{2}} (\Gamma_{|u|+2L})^{-\frac{1}{2L}}$ is more appropriate. Corollary 1 provides the results for highly smooth functions. To be able to derive the optimal rates of convergence for a wide class of functions (i.e., $\mathcal{H}_{|u|+1}$), Corollary 2 starts providing the biases for this class of functions under a specific set of constraints. To that end, define

$$ K_1 := \inf_{\mathbf{R} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| \mathbf{R} \right\|_2 \prod_{k \in u} R_k^2 \right]. $$

**Corollary 2.** For distinct $\beta_\ell$s, consider $r^* \in \{0, \dots, |u| - 1\}$ and the constraints $\sum_{\ell=1}^{L=r^*+2} C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$ with $r = 0, 1, \dots, r^*, |u|$. If $f \in \mathcal{H}_{|u|+1}$ and (A2) hold, then there is $M_{|u|+1} > 0$ such that

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

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

$$ \left| \mathcal{D}^{|u|} f(\mathbf{x}) - \sum_{\ell=1}^{L=2} C_{\ell}^{(|u|)} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in u} \frac{V_k}{h_k \sigma^2} \right] \right| \le \zeta M_{|u|+1} \left\| \mathbf{h} \right\|_1 \Gamma_{|u|+1}. \quad (7) $$
**Proof.** See Appendix C. □

Note that the upper bounds derived in Corollary 2 depend on $r^* \in \{0, \dots, |u| - 1\}$ through $L = r^* + 2$ and $\Gamma_{|u|+1} = \sum_{\ell=1}^{L=r^*+2} |C_{\ell}^{(|u|)} \beta_{\ell}^{(|u|+1)}|$. Thus, taking $h_k = h$ and $\zeta = 1/(d\Gamma_{|u|+1})$ will give a dimension-free upper bound that does not increase with $r^*$. The crucial role and importance of $r^*$ is highlighted in Section 3.3.

**Remark 4.** When $r^* = 0$, which implies that $L = 2$, consider $\beta_1 = 1, \beta_2 = -1$; $C_1^{(|u|)} = 1/2$; $C_2^{(|u|)} = 1/2$ when $|u|$ is even and $C_2^{(|u|)} = -1/2$ otherwise. With the above choices, the upper bounds given by Equations (6) and (7) become, respectively,

$$ \sigma M_{|u|+1} K_1 ||\mathbf{h}||_2, \qquad \zeta M_{|u|+1} ||\mathbf{h}||_1. $$

We can check that the same results hold when using $L = 1, \beta_1 = 1$ and $C_1^{(|u|)} = 1$.

**Remark 5.** It is worth noting that we obtain exact approximations of $\mathcal{D}^{|u|}f(\mathbf{x})$ in Corollary 1 for the class of functions described by

$$ \mathcal{B}_0 := \{h \in \mathcal{H}_{\infty} : \mathcal{D}^{(\bar{t})} h = 0, \quad \forall \bar{t} \in \mathbb{N}^d \text{ and } ||\bar{t}||_1 \ge |u| + 2L \}. $$

In general, exact approximations of $\mathcal{D}^{|u|}f(\mathbf{x})$ are obtained when $L \to \infty$ for highly smooth functions.

### 3.3. Convergence Analysis

Given a sample of $\mathbf{V}$, that is, $\{\mathbf{V}_i := (V_{i,1}, \dots, V_{i,d})\}_{i=1}^N$ and using Equation (2), the method of moments implies that the estimator of $\mathcal{D}^{|u|}f(\mathbf{x})$ is given by

$$ \widehat{\mathcal{D}^{|u|}f}(\mathbf{x}) := \frac{1}{N} \sum_{i=1}^{N} \sum_{\ell=1}^{L} C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}_i) \prod_{k \in u} \frac{V_{i,k}}{(h_k \sigma^2)}. $$

Statistically, it is common to measure the quality of an estimator using the mean squared error (MSE), including the rates of convergence. The MSEs can also help in determining the optimal value of $\mathbf{h}$. Theorem 2 provides such quantities under different assumptions. To that end, define

$$ K_{2,r} := \inf_{\mathbf{R} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| \mathbf{R}^2 \right\|_2^r \prod_{k \in u} R_k^2 \right]. $$

**Theorem 2.** For distinct $\beta_{\ell}$s, consider $r^* \in \mathbb{N}$ and $\sum_{\ell=1}^{L=r^*+2} C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$ with $r = 0, 1, \dots, r^*, |u|$ and $r^* \le |u| - 1$. If $f \in \mathcal{H}_{|u|+1}$ and (A2) hold, then

$$ \mathbb{E} \left[ \left( \widehat{\mathcal{D}^{|u|}f}(\mathbf{x}) - \mathcal{D}^{|u|}f(\mathbf{x}) \right)^2 \right] \le \sigma^2 M_{|u|+1}^2 K_1^2 \Gamma_{|u|+1}^2 ||\mathbf{h}||_2^2 + \frac{M_{r^*+1}^2 \Gamma_{r^*+1}^2 K_{2,r^*+1}}{N \sigma^2 (|u|-r^*-1) \prod_{k \in u} h_k^2} ||\mathbf{h}||_2^{r^*+1}. \quad (8) $$

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

$$ \mathbb{E} \left[ \left( \widehat{\mathcal{D}^{|u|}f}(\mathbf{x}) - \mathcal{D}^{|u|}f(\mathbf{x}) \right)^2 \right] \le M_{|u|+1}^2 ||\mathbf{h}||_1^2 \zeta^2 \Gamma_{|u|+1}^2 + \frac{3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2}{N \zeta^2 (|u|-r^*-1) \prod_{k \in u} h_k^2} ||\mathbf{h}||_2^{2(r^*+1)}. \quad (9) $$

**Proof.** See Appendix D for the detailed proof. □

Theorem 2 provides the upper bounds of MSEs for the anisotropic case. Using a uniform bandwidth, that is, $h_k = h$, reveals that such upper bounds clearly depend on the
dimensionality of the function of interest. Indeed, we can check that the upper bounds of the MSEs provided in Equations (8) and (9) become, respectively,

$$
\begin{aligned}
& \sigma^2 M_{|u|+1}^2 K_1^2 \Gamma_{|u|+1}^2 dh^2 + \frac{M_{r^*+1}^2 \Gamma_{r^*+1}^2 K_{2,r^*+1}}{N \sigma^2 (|u|-r^*-1) h^2 (|u|-r^*-1)} d^{\frac{r^*+1}{2}}, \\
& \zeta^2 M_{|u|+1}^2 \Gamma_{|u|+1}^2 d^2 h^2 + \frac{3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2 d^{r^*+1}}{N \zeta^2 (|u|-r^*-1) h^2 (|u|-r^*-1)}.
\end{aligned}
$$

By minimizing such upper bounds with respect to h, the optimal rates of convergence of the proposed estimators are derived in Corollary 3.

**Corollary 3.** Under the assumptions made in Theorem 2, if $r^* < |u| - 1$, then

$$
\mathbb{E} \left[ \left( \widehat{\mathcal{D}|u|f}(\mathbf{x}) - \mathcal{D}|u|f(\mathbf{x}) \right)^2 \right] = \mathcal{O} \left( N^{-\frac{1}{|u|-r^*}} d^{1+\frac{|u|-1}{|u|-r^*}} \right).
$$

**Proof.** See Appendix E for the detailed proof. □

The optimal rates of convergence obtained in Corollary 3 are far away from the parametric ones, and such rates decrease with $|u|$. Nevertheless, such optimal rates are function of $d^2$ for any $|u| \ge 2$ using $r^* = 1$. The maximum rate of convergence that can be reached is $N^{-1/2}d^{\frac{|u|+1}{2}}$ by taking $r^* = |u| - 2$.

To derive the optimal and parametric rates of convergence, let us now choose $r^* = |u| - 1$ and $h_k = h$ with $k = 1, \dots, d$. Thus, we can see that the second terms of the upper bounds of the MSEs (provided in Theorem 2) are function of $d^{|u|}$, but they are independent of $h$. This key observation leads to Corollary 4.

**Corollary 4.** Under the assumptions made in Theorem 2, if $r^* = |u| - 1$; $\zeta = 1/(d\Gamma_{|u|+1})$ and $h_k = h \propto N^{-\gamma/2}$ with $\gamma \in [1, 2]$, then we have

$$
\mathbb{E} \left[ \left( \widehat{\mathcal{D}|u|f}(\mathbf{x}) - \mathcal{D}|u|f(\mathbf{x}) \right)^2 \right] = \mathcal{O} \left( N^{-1} d^{|u|} \right).
$$

**Proof.** The proof is straightforward since $h^2 \propto N^{-\gamma}$ and $Nh \to \infty$ if $N \to \infty$. □

It is worth noting that the upper bound of the squared bias obtained in Corollary 4 does not depend on the dimensionality thanks to $\zeta$. Also, the optimal and parametric rates of convergence are reached by means of $LN = (|u| + 1)N$ model evaluations, and such model runs can still be used for computing $\mathcal{D}^{|v|}f(\mathbf{x})$ for every $v \subseteq \{1, \dots, d\}$ with $|v| \le |u|$. Based on the same assumptions, it appears that the results provided in Corollary 4 are much more convenient for $|u| = 1, 2$ in higher dimensions, while those obtained in Corollary 3 are well suited for higher dimensions and for higher values of $|u| \in \{1, \dots, d\}$.

For highly smooth functions and for large values of $|u|$ and $d$, we are able to derive intermediate rates of convergence of the estimator of $\mathcal{D}^{|u|}f$ (see Theorem 3). To that end, consider an integer $L'$, $L = r^* + L' + 1$, and denote with $[b]$ the largest integer that is less than $b$ for any real $b$.

**Theorem 3.** For an integer $r^* \le |u| - 2$, consider $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$, $r = 0, 1, \dots, r^*$, $|u|$, $|u| + 2$, $|u| + 4, \dots, |u| + 2(L' - 1)$. If $f \in \mathcal{H}_{|u|+2L'}$ and (A2) hold, then

$$
\mathbb{E} \left[ \left( \widehat{\mathcal{D}|u|f}(\mathbf{x}) - \mathcal{D}|u|f(\mathbf{x}) \right)^2 \right] \le \sigma^{4L'} M_{|u|+2L'}^2 K_{1,L'}^2 d^{L'} h^{4L'} + \frac{M_{r^*+1}^2 \Gamma_{r^*+1}^2 K_{2,r^*+1}}{N \sigma^2 (|u|-r^*-1) h^2 (|u|-r^*-1)} d^{\frac{r^*+1}{2}}.
$$
Moreover, if $V_k \sim \mathcal{U}(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, then

$$ \mathbb{E} \left[ \left( \widehat{\mathcal{D}^{|u|} f}(\mathbf{x}) - \mathcal{D}^{|u|} f(\mathbf{x}) \right)^2 \right] \le \xi^{4L'} M_{|u|+2L'}^2 \Gamma_{|u|+2L'}^2 d^{2L'} h^{4L'} + \frac{3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2 d^{r^*+1}}{N \xi^2 (|u|-r^*-1) h^2 (|u|-r^*-1)}. $$

For a given $0 \le \epsilon_{op} < 1$, taking $L' = \left[ \frac{(|u|-r^*-1)(1-\epsilon_{op})}{2\epsilon_{op}} \right]$ leads to

$$ \mathbb{E} \left[ \left( \widehat{\mathcal{D}^{|u|} f}(\mathbf{x}) - \mathcal{D}^{|u|} f(\mathbf{x}) \right)^2 \right] = \mathcal{O} \left( N^{-1+\epsilon_{op}} d^{|u|(1-\epsilon_{op})} \right). \quad (10) $$

**Proof.** Detailed proofs are provided in Appendix F. □

It turns out that the optimal rate of convergence derived in Theorem 3 is a trade-off between the sample size $N$ and the dimensionality $d$. For instance, when $\epsilon_{op} = 1/2$, the optimal rate becomes $\mathcal{O}(N^{-1/2}d^{|u|/2})$, which improves the rate obtained in Corollary 3, but under different assumptions.

**Remark 6.** Since the bias vanishes for the class of functions $\mathcal{B}_0$, taking $\sigma = 1/h$ yields $\mathbb{E} \left[ \left( \widehat{\mathcal{D}^{|u|} f}(\mathbf{x}) - \mathcal{D}^{|u|} f(\mathbf{x}) \right)^2 \right] = \mathcal{O}(N^{-1})$ (see Appendix G). Note that such an optimal rate of convergence is dimension-free.

### 3.4. Derivative-Based Emulators of Smooth Functions

Using Equation (1) and bearing in mind the estimators of the cross-derivatives provided in Section 3.3, this section aims at providing surrogates of smooth functions, also known as emulators. The general expression of the surrogate of $f$ is given below.

**Corollary 5.** For any $|u| \in \{1, \dots, d\}$, consider $\sum_{\ell=1}^{L=5} C_{\ell}^{(|u|)} \beta_{\ell}^r = \delta_{|u|,r}$, $r = 0, 1, \dots, r^* = \max(|u| - 1, 3)$, $|u|$. Assume that $f \in \mathcal{H}_{d+1}$ and (A1) and (A2) hold. Then, an approximation of $f$ at $\mathbf{x} \in \Omega$ is given by

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

The above plug-in estimator is consistent using the law of large numbers. For a given $\mathbf{x}$, it is worth noting that the choice of $G_j$ is arbitrary, provided that such distributions are supported on an open neighborhood of $\mathbf{x}$.

Often, the higher-order cross-partial derivatives or, equivalently, the higher-order interactions among the model inputs almost vanish, leading us to consider the truncated expressions. Given an integer $s$ with $0 < s \ll d$ and keeping in mind the ANOVA decomposition, consider the class of functions that admit at most the $s$-th-order interactions, that is, $A_s := \{ h : \mathbb{R}^d \to \mathbb{R} : h(\mathbf{x}) = \sum_{\substack{v \subseteq \{1, \dots, d\} \\ |v| \le s}} h_v(\mathbf{x}_v) \}$. Truncating the functional expansion is a standard practice within the ANOVA-community, that is, $s \ll d$ is assumed in higher dimensions [10,12]. For such a class of functions, requiring $f \in \mathcal{H}_{\alpha_s}$ with $\alpha_s \ge s$ is sufficient to derive our results. Thus, the truncated surrogate of $f$ is given by

$$ \hat{f}_s(\mathbf{x}) := \mathbb{E}[f(\mathbf{X}')] + \sum_{\substack{v \subseteq \{1, \dots, d\} \\ 0 < |v| \le s}} \mathbb{E}_{\mathbf{X}'} \left[ \widehat{\mathcal{D}^{|v|} f}(\mathbf{X}') \prod_{k \in v} \frac{G_k(\mathbf{X}'_k) - \mathbb{I}_{[\mathbf{X}'_k \ge x_k]}}{g_k(\mathbf{X}'_k)} \right]. $$
Under the assumptions made in Corollary 5, $\hat{f}_s(\mathbf{x})$ reaches the optimal and parametric rate of convergence for the class of functions $\mathcal{A}_3$. For instance, taking $s = 1$ leads to the first-order emulator of $f$, which relies only on the gradient information. Thus, $\hat{f}_{s=1}$ provides accurate estimates of additive models of the form $\sum_{j=1}^d h_j(x_j)$, where $h_j$s are given functions. Likewise, $\hat{f}_2$ allows for incorporating the second-order terms, but it requires the second-order cross-partial derivatives. Thus, it is relevant to find the class of functions $\mathcal{A}_s$ which contains the model of interest before building emulators. The following section deals with such issues.

## 4. Applications: Computing Sensitivity Indices

In high-dimensional settings, reducing the dimension of functions is often achieved by using screening measures, that is, measures that can be used for quickly identifying non-relevant input variables. Screening measures based on the upper bounds of the total sensitivity indices rely on derivatives [14,23–25,60,61]. This section aims at providing optimal computations of upper bounds of the total indices, followed by the computations of the main indices using derivatives.

By evaluating the function $f$ given by Equation (1) at a random vector $\mathbf{X}$ using $G_j = F_j$, one obtains a random vector of the model outputs. Generalized sensitivity indices, including Sobol's indices, rely on the variance-covariance of sensitivity functionals (SFs), which are also random vectors containing the information about the overall contributions of inputs [14,26,62]. The derivative-based expressions of SFs are given below (see [14] for more details). Given $u \subseteq \{1, \dots, d\}$, the interaction SF of the inputs $\mathbf{X}_u$ is given by

$$f_u(\mathbf{X}_u) := \mathbb{E}_{\mathbf{X}'_u} \left[ \mathcal{D}^{|u|} f(\mathbf{X}') \prod_{k \in u} \frac{F_k(\mathbf{X}'_k) - \mathbb{I}_{[\mathbf{X}'_k \ge \mathbf{X}_k]}}{\rho_k(\mathbf{X}'_k)} \right],$$

and the first-order SF of $\mathbf{X}_u$ is given by

$$f_u^{fo}(\mathbf{X}_u) := \sum_{\substack{v, v \subseteq u \\ |v| > 0}} f_v(\mathbf{X}_v).$$

Likewise, the total-interaction SF of $\mathbf{X}_u$ is given by [14]

$$f_u^{sup}(\mathbf{X}) := \sum_{\substack{v \subseteq \{1, \dots, d\} \\ u \subseteq v}} f_v(\mathbf{X}_v) = \mathbb{E}_{\mathbf{X}'_u} \left[ \mathcal{D}^{|u|} f(\mathbf{X}'_u, \mathbf{X}_u) \prod_{k \in u} \frac{F_k(\mathbf{X}'_k) - \mathbb{I}_{[\mathbf{X}'_k \ge \mathbf{X}_k]}}{\rho_k(\mathbf{X}'_k)} \right],$$

and the total SF of $\mathbf{X}_u$ is given as [14]

$$f_u^{tot}(\mathbf{X}) := \sum_{\substack{v \subseteq \{1, \dots, d\} \\ u \cap v \neq \emptyset}} f_v(\mathbf{X}_v) = \sum_{\substack{v, v \subseteq u \\ |v| > 0}} \mathbb{E}_{\mathbf{X}'_u} \left[ \mathcal{D}^{|v|} f(\mathbf{X}'_u, \mathbf{X}_u) \prod_{k \in v} \frac{F_k(\mathbf{X}'_k) - \mathbb{I}_{[\mathbf{X}'_k \ge \mathbf{X}_k]}}{\rho_k(\mathbf{X}'_k)} \right].$$

For a single input $X_j$, we have $f_j(X_j) = f_j^{fo}(X_j)$ and $f_j^{tot}(\mathbf{X}) = f_j^{sup}(\mathbf{X})$. Among similarity measures [28,29], taking the variance-covariances of SFs, that is, $\Sigma_u := \mathbb{V}[f_u(\mathbf{X}_u)]$, $\Sigma_u^{sup} := \mathbb{V}[f_u^{sup}(\mathbf{X})]$ and $\Sigma_u^{tot} := \mathbb{V}[f_u^{tot}(\mathbf{X})]$, leads to [14]

$$\Sigma_u = \mathbb{E} \left[ \mathcal{D}^{|u|} f(\mathbf{X}) \mathcal{D}^{|u|} f(\mathbf{X}') \prod_{k \in u} \frac{F_k(\min(\mathbf{X}_k', \mathbf{X}'_k)) - F_k(\mathbf{X}_k) F_k(\mathbf{X}'_k)}{\rho_k(\mathbf{X}_k) \rho_k(\mathbf{X}'_k)} \right]; \quad (12)$$

$$\Sigma_u \le \Sigma_u^{sup} \le \Sigma_u^{ub} := \frac{1}{2^{|u|}} \mathbb{E} \left[ \left( \mathcal{D}^{|u|} f(\mathbf{X}) \right)^2 \prod_{k \in u} \frac{F_k(\mathbf{X}_k) [1 - F_k(\mathbf{X}_k)]}{(\rho_k(\mathbf{X}_k))^2} \right]. \quad (13)$$
Thus, $\Sigma_u^{ub}$ is the upper bound of $\Sigma_u^{sup}$. Likewise, $\Sigma_j^{ub}$ is the upper bound of $\Sigma_j^{tot}$ (i.e., $\Sigma_j^{tot} \le \Sigma_j^{ub}$), and it can be used for screening the input variables.

To provide new expressions of the screening measures and the main induces in the following proposition, denote with $\mathbf{V}'$ an i.i.d. copy of $\mathbf{V}$, and assume that

**Assumption 3 (A3).** $f(\mathbf{X})$ has finite second-order moments.

**Proposition 1.** Under the assumptions made in Corollary 4, assume that (A1)–(A3) hold. Then,

$$ \Sigma_u = \sum_{\substack{\ell_1=1 \\ \ell_2=1}}^{L} C_{\ell_1}^{(|u|)} C_{\ell_2}^{(|u|)} \mathbb{E}[f(\mathbf{X} + \beta_{\ell_1} \mathbf{h} \mathbf{V}) f(\mathbf{X}' + \beta_{\ell_2} \mathbf{h} \mathbf{V}') \\ \times \prod_{k \in u} \frac{V_k V_k'}{(h_k^2 \sigma^4)} \left( \frac{F_k(\min(X_k, X_k')) - F_k(X_k)F_k(X_k')}{\rho_k(X_k)\rho_k(X_k')} \right)] + \mathcal{O}\left(\|\mathbf{h}\|_2^2\right); \quad (14) $$

$$ \Sigma_u^{ub} = \frac{1}{2^{|u|}} \sum_{\substack{\ell_1=1 \\ \ell_2=1}}^{L} C_{\ell_1}^{(|u|)} C_{\ell_2}^{(|u|)} \mathbb{E}[f(\mathbf{X} + \beta_{\ell_1} \mathbf{h} \mathbf{V}) f(\mathbf{X} + \beta_{\ell_2} \mathbf{h} \mathbf{V}') \\ \times \prod_{k \in u} \frac{V_k V_k'}{h_k^2 \sigma^4} \left( \frac{F_k(X_k)[1 - F_k(X_k)]}{(\rho_k(X_k))^2} \right)] + \mathcal{O}\left(\|\mathbf{h}\|_2^2\right). \quad (15) $$

**Proof.** Bearing in mind Equations (2), (12) and (13), the proof of Proposition 1 is straightforward. $\square$

The method of moments allows for deriving the estimators of $\Sigma_u$ and $\Sigma_u^{ub}$ for all $u \subseteq \{1, \dots, d\}$. For screening inputs of models, we provide the estimators of $\Sigma_j$ and $\Sigma_j^{ub}$ for any $j \in \{1, \dots, d\}$. To that end, we are given four independent samples, that is, $\{\mathbf{X}_i\}_{i=1}^N$ from $\mathbf{X}$, $\{\mathbf{X}'_i\}_{i=1}^N$ from $\mathbf{X}'$, $\{\mathbf{V}_i\}_{i=1}^N$ from $\mathbf{V}$, and $\{\mathbf{V}'_i\}_{i=1}^N$ from $\mathbf{V}'$. Consistent estimators of $\Sigma_j$ and $\Sigma_j^{ub}$ are, respectively, given by

$$ \hat{\Sigma}_j := \sum_{i=1}^{N} \sum_{\substack{\ell_1=1 \\ \ell_2=1}}^{L} \frac{C_{\ell_1}^{(1)} C_{\ell_2}^{(1)}}{N} f(\mathbf{X}_i + \beta_{\ell_1} \mathbf{h} \mathbf{V}_i) f(\mathbf{X}'_i + \beta_{\ell_2} \mathbf{h} \mathbf{V}'_i) \frac{V_{i,j} V'_{i,j}}{h_j^2 \sigma^4} \\ \left( \frac{F_j(\min(X_{i,j}, X'_{i,j})) - F_j(X_{i,j})F_j(X'_{i,j})}{\rho_j(X_{i,j})\rho_j(X'_{i,j})} \right); \quad (16) $$

$$ \widehat{\Sigma}_j^{ub} := \sum_{i=1}^{N} \sum_{\substack{\ell_1=1 \\ \ell_2=1}}^{L} \frac{C_{\ell_1}^{(1)} C_{\ell_2}^{(1)}}{2N} f(\mathbf{X}_i + \beta_{\ell_1} \mathbf{h} \mathbf{V}_i) f(\mathbf{X}_i + \beta_{\ell_2} \mathbf{h} \mathbf{V}'_i) \frac{V_{i,j} V'_{i,j}}{h_j^2 \sigma^4} \left( \frac{F_j(\mathbf{X}_{i,j})[1 - F_j(\mathbf{X}_{i,j})]}{(\rho_j(\mathbf{X}_{i,j}))^2} \right). \quad (17) $$

The above (direct) estimators require 3LN model runs for obtaining the estimates for any $j \in \{1, \dots, d\}$. Additionally to such estimators, we derive the plug-in estimators, which are relevant in the presence of given data about the estimates of the first-order derivatives. To provide such estimators, we denote with $\{\widehat{\mathcal{D}}^{\{j\}}|f(\mathbf{X}_i)\}_{i=1}^{N_1}$ a sample of $N_1 > 1$ known or estimates of the first-order derivatives (i.e., $|u| = 1$). Using Equation (2), such estimates are obtained by considering $L = 1$ or $L = 2$ or $L = 3$ and $N$. Keeping in mind
Equations (12) and (13) and the U-statistic theory for one sample, the plug-in estimator of the main index of $X_j$ is given by

$$ \hat{\Sigma}'_j := \frac{2}{N_1(N_1-1)} \sum_{1 \le i_1 < i_2 \le N_1} \widehat{D^{\{j\}}|f(X_{i_1})} \widehat{D^{\{j\}}|f(X_{i_2})} \frac{F_j(\min(X_{i_1,j}, X_{i_2,j})) - F_j(X_{i_1,j})F_j(X_{i_2,j})}{\rho_j(X_{i_1,j})\rho_j(X_{i_2,j})} . \quad (18) $$

Likewise, the plug-in estimator of the upper bound of the total index of $X_j$ is given by

$$ \widehat{\Sigma}'_{j}^{ub} := \frac{1}{2N_1} \sum_{i=1}^{N_1} \left[ \widehat{D^{\{j\}}|f(X_i)} \right]^2 \frac{F_j(X_{i,j}) [1 - F_j(X_{i,j})]}{[\rho_j(X_{i,j})]^2} . \quad (19) $$

Note that the plug-in estimators are consistent and require a total of $LN_1N_0$ model runs for computing such indices, where $N_0$ is the number of model runs used for computing the gradient of $f$ at $X_i$.

# 5. Illustrations: Screening and Emulators of Models

## 5.1. Test Functions

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

The Ishigami function includes three independent inputs following a uniform distribution on $[-\pi, \pi]$, and it is given by

$$ f(\mathbf{x}) = \sin(x_1) + 7 \sin^2(x_2) + 0.1 x_3^4 \sin(x_1). $$

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$.

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

The g-function [63] includes ten independent inputs following a uniform distribution on $[0, 1]$, and it is defined as follows:

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

Note that such a function is differentiable almost everywhere. According to the values of $\mathbf{a} = (a_j, j = 1, 2, \dots, d)$, this function has different properties [23]:

* If $\mathbf{a} = [0, 0, 6.52, 6.52, 6.52, 6.52, 6.52, 6.52, 6.52]^T$, the values of sensitivity indices are $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$. Thus, this function has a low effective dimension (function of type A), and it belongs to $\mathcal{A}_s$ with $s > 1$ (see Section 3.4).

* If $\mathbf{a} = [50, 50, 50, 50, 50, 50, 50, 50, 50]^T$, the first and total indices are given as follows: $S_j = S_{T_j} = 0.1, \forall j \in \{1, 2, \dots, d\}$. Thus, all inputs are important, but there is no interaction among these inputs. This function has a high effective dimension (function of type B). Note that it belongs to $\mathcal{A}_1$.

* If $\mathbf{a} = [0, 0, 0, 0, 0, 0, 0, 0, 0]^T$, the function belongs to the class of functions with important interactions among inputs. Indeed, we have $S_j = 0.02$ and $S_{T_j} = 0.27, \forall j \in \{1, 2, \dots, d\}$. All the inputs are relevant due to important interactions (function of type C). Then, this function belongs to $\mathcal{A}_s$ with $s \ge 2$.

## 5.2. Numerical Comparisons of Estimators

This section provides a comparison of the direct and plug-in estimators of the main indices and the upper bounds of the total indices using the test functions of Section 5.1. Different total budgets for the model evaluations are considered in this paper, that is, $N = 500, 1000, 1500, 2000, 3000, 5000, 10,000, 15,000, 20,000$. In the case of the plug-in estimators, we used $N_0 = 2d$. For generating different random values, Sobol's sequence
(scrambled = 3) from the R-package `randtoolbox` [64] is used. We replicated each estimation 30 times by randomly choosing the seed, and the reported results are the average of the 30 estimates.

Figures 1–4 show the mean squared errors related to the estimates of the main indices for the Ishigami function, the g-functions of type A, type B, and type C, respectively. Each figure depicts the results for $L = 2$ and $L = 3$. All the figures show the convergence and accuracy of our estimates using either $L = 2$ or $L = 3$.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Average of $d = 3$ mean squared errors using the Ishigami function (∘ the direct estimator (16) and + for the plug-in estimator (18)).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Average of $d = 10$ mean squared errors using the g-function of type A (∘ the direct estimator (16) and + for the plug-in estimator (18)).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Average of *d* = 10 mean squared errors using the g-function of type B (◦ the direct estimator (16) and + for the plug-in estimator (18)).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Average of *d* = 10 mean squared errors using the g-function of type C (◦ the direct estimator (16) and + for the plug-in estimator (18)).</figcaption></figure>

Likewise, Figures 5-8 show the mean squared gaps (differences) between the true total index and its estimated upper bound for the Ishigami function, the g-functions of type A, type B, and type C, respectively.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Average of *d* = 3 mean squared gaps using the Ishigami function (∘ the direct estimator (17) and + for the plug-in estimator (19)).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Average of *d* = 10 mean squared gaps using the g-function of type A (∘ the direct estimator (17) and + for the plug-in estimator (19)).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Average of *d* = 10 mean squared gaps using the g-function of type B (◦ the direct estimator (17) and + for the plug-in estimator (19)).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Average of *d* = 10 mean squared gaps using the g-function of type C (◦ the direct estimator (17) and + for the plug-in estimator (19)).</figcaption></figure>

It turns out that the plug-in estimators outperform the direct ones. Also, increasing the values of *L* gives the same results. Moreover, the direct estimators associated with *L* = 1 fail to provide accuracy estimates (we do not report such results here). On the contrary, the plug-in estimates using *L* = 1 are reported in Table 1 for the Ishigami function and in Table 2 for the three types of the g-function. Such results suggest considering *L* = 1 or *L* = 2 with $\beta_1 = -1, \beta_2 = 1$ for plug-in estimators when the total budget of model runs is small. For a larger budget of model runs, the direct estimators associated with *L* = 2 and $\beta_1 = -1, \beta_2 = 1$ can be considered as well in practice.
<table><caption>**Table 1.** Average of 30 estimates of the main indices and upper bounds of total indices for the Ishigami function using the plug-in estimators, $L = 1$, and 2000 model runs.</caption><thead><tr><th></th><th>X1</th><th>X2</th><th>X3</th></tr></thead><tbody><tr><td>S<sub>j</sub></td><td>0.249</td><td>0.318</td><td>-0.006</td></tr><tr><td>UB<sub>j</sub></td><td>1.420</td><td>4.872</td><td>0.711</td></tr></tbody></table>

<table><caption>**Table 2.** Average of 30 estimates of the main indices and upper bounds of total indices for the g-functions using the plug-in estimators, $L = 1$, and 2000 model runs.</caption><thead><tr><th></th><th>X1</th><th>X2</th><th>X3</th><th>X4</th><th>X5</th><th>X6</th><th>X7</th><th>X8</th><th>X9</th><th>X10</th></tr></thead><tbody><tr><td colspan="11" style="text-align: center;">Type A</td></tr><tr><td>S<sub>j</sub></td><td>0.330</td><td>0.324</td><td>0.006</td><td>0.005</td><td>0.006</td><td>0.006</td><td>0.005</td><td>0.005</td><td>0.006</td><td>0.005</td></tr><tr><td>UB<sub>j</sub></td><td>2.022</td><td>2.005</td><td>0.045</td><td>0.046</td><td>0.047</td><td>0.047</td><td>0.046</td><td>0.046</td><td>0.046</td><td>0.047</td></tr><tr><td colspan="11" style="text-align: center;">Type B</td></tr><tr><td>S<sub>j</sub></td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td><td>0.085</td></tr><tr><td>UB<sub>j</sub></td><td>0.362</td><td>0.363</td><td>0.363</td><td>0.362</td><td>0.363</td><td>0.362</td><td>0.363</td><td>0.362</td><td>0.363</td><td>0.363</td></tr><tr><td colspan="11" style="text-align: center;">Type C</td></tr><tr><td>S<sub>j</sub></td><td>0.028</td><td>0.028</td><td>0.032</td><td>0.032</td><td>0.035</td><td>0.041</td><td>0.031</td><td>0.030</td><td>0.036</td><td>0.034</td></tr><tr><td>UB<sub>j</sub></td><td>2.033</td><td>1.301</td><td>1.825</td><td>1.605</td><td>1.634</td><td>1.641</td><td>2.216</td><td>1.526</td><td>1.793</td><td>1.503</td></tr></tbody></table>

### 5.3. Emulations of the g-Function of Type B

Based on the results obtained in Section 5.2 (see Table 2), all the inputs are important in the case of the g-function of type B, meaning that the dimension reduction is not possible. Also, the estimated upper bounds suggest weak interactions among inputs. As expected, our estimated results confirm that the g-function of type B belongs to $\mathcal{A}_1$. Thus, an emulator of $f$ based only on the first-order derivatives is sufficient. Using this information, we have derived the emulator of that function (i.e., $\hat{f}_{s=1}(x)$) under the assumptions made in Corollary 5 ($r^* = 0, L \ge 2$) and using $G_j = F_j$ with $j = 1, \dots, d$. For a given $L$, we used 300 model runs to build the emulator, and Figure 9 depicts the approximations of that function (called predictions) at the sample points involved in the construction of the emulator. Note that the evaluations of $f$ at such sample points (called observations) are not directly used in the construction of such an emulator. It turns out from Figure 9 that $\hat{f}_{s=1}$ provides predictions that are in line with the observations, showing the accuracy of our emulator.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Predictions of g-function of type B using the emulator $\hat{f}_{s=1}$ versus observations using $L = 2$ and $L = 3$.</figcaption></figure>
## 6. Conclusions

In this paper, we firstly provided i) stochastic expressions of cross-partial derivatives of any order, followed by their biases, and ii) estimators of such expressions. Our estimators of the $|u|$-th cross-partial derivatives ($\forall u \subseteq \{1, ..., d\}$) reach the parametric rates of convergence (i.e., $O(N^{-1}d^{|u|})$) by means of a set of $L \ge |u| + 1$ constraints for the Hölder space of $\alpha$-smooth functions $H_\alpha$ with $\alpha > |u|$. Moreover, we showed that the upper bounds of the biases of such estimators do not suffer from the curse of dimensionality. Secondly, the proposed surrogates of cross-partial derivatives are used for deriving (i) new derivative-based emulators of simulators or surrogates of models, even when a large number of model inputs contribute to the model outputs, and (ii) new repressions of the main sensitivity indices and the upper bounds of the total sensitivity indices.

Numerical simulations confirmed the accuracy of our approaches for not only screening the input variables, but also for identifying the class of functions that contains our simulator of interest, such as the class of functions with important or no interaction among inputs. This relevant information allows for designing and building the appropriate emulators of functions. In the case of the g-function of type B, our emulator of this function (based only on the first-order derivatives) provided approximations or predictions that are in line with the observations.

For functions with important interactions or, equivalently, for higher-order cross-partial derivatives, further numerical schemes are necessary to increase the numerical accuracy of the computations of such derivatives and predictions. Such perspectives will be investigated in the near future as well as the computations of the total sensitivity indices using the proposed surrogates of derivatives. Moreover, there is a need for a theoretical investigation to expect derivation of the parametric rates of convergence of the above estimators that do not suffer from the course of dimensionality. Working in $\mathbb{C}$ rather than in $\mathbb{R}$ may be helpful.

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

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** I would like to thank the three reviewers for their comments that helped improve my manuscript.

**Conflicts of Interest:** The author declares no conflict of interest.

## Appendix A. Proof of Theorem 1

Firstly, as $\vec{i} := (i_1, ..., i_d) \in \mathbb{N}^d$, denote $||\vec{i}||_1 = i_1 + ... + i_d$ and $\vec{u} := (\mathbb{I}_u(1), ..., \mathbb{I}_u(d))$. 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_{p=0}^{\alpha} \sum_{||\vec{i}||_1=p} \frac{\mathcal{D}(\vec{i}) f(\mathbf{x})}{\vec{i}!} \beta_{\ell}^{p} (\mathbf{h}\mathbf{V})^{\vec{i}} + \mathcal{O}(||\beta_{\ell} \mathbf{h}\mathbf{V}||_1^{\alpha+1}).$$

Multiplying such an expansion by the constant $C_\ell^{(|u|)}$, and taking the sum over $\ell = 1, ..., L$, the expectation $E := \sum_{\ell=1}^L C_\ell^{(|u|)} \mathbb{E}[f(\mathbf{x} + \beta_\ell \mathbf{h}\mathbf{V}) \prod_{k \in u} \frac{V_k}{(h_k \sigma^2)}]$ becomes

$$E = \sum_{p \ge 0} \sum_{||\vec{i}||_1 = p} \frac{\mathcal{D}(\vec{i}) f(\mathbf{x})}{\vec{i}!} \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^{p} \right) \mathbb{E} \left[ \frac{(\mathbf{V})^{\vec{i}+\vec{u}} (\mathbf{h})^{\vec{i}-\vec{u}}}{\sigma^{2|u|}} \right].$$

We can see that $\mathbb{E}[(\mathbf{V})^{\vec{i}+\vec{u}}(\mathbf{h})^{\vec{i}-\vec{u}}] \ne 0$ iff $\vec{i} + \vec{u} = 2\vec{q}$, $\forall \vec{q} \in \mathbb{N}^d$. Equation $\vec{i} + \vec{u} = 2\vec{q}$ implies $i_k = 2q_k \ge 0$ if $k \notin u$ and $i_k = 2q_k - 1 \ge 0$ otherwise. Thus, using $i_k = 2q_k + 1$
when $k \in u$ is much more convenient, and it leads to $\vec{i} = 2\vec{q} + \vec{u}, \forall \vec{q} \in \mathbb{N}^d$, which also implies that $||\vec{i}||_1 \ge ||\vec{u}||_1$. We then obtain $\mathcal{D}^{|u|}f$ when $||\vec{q}||_1 = 0$ or $\vec{i} = \vec{u}$, and the fact that $\mathbb{E}[(\mathbf{V})^{2\vec{u}}] = \mathbb{E}[\Pi_{k \in u} V_k^2] = \sigma^{2|u|}$ by independence. We can then write

$$E = \mathcal{D}^{|u|} f(\mathbf{x}) \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^{|u|} \right) + \sum_{\substack{r \ge 1 \\ ||\vec{q}||_1 = r}} \frac{\mathcal{D}^{(2\vec{q}+\vec{u})} f(\mathbf{x})}{(2\vec{q}+\vec{u})!} \left( \sum_{\ell} C_{\ell}^{(|u|)} \beta_{\ell}^{r+|u|} \right) \mathbb{E} \left[ \frac{(\mathbf{V})^{2(\vec{q}+\vec{u})} (\mathbf{h})^{2\vec{q}}}{\sigma^{2|u|}} \right],$$

using the change of variable $r = p - ||\vec{u}||_1$. At this point, setting $L = 1, \beta_\ell = 1$, and $C_\ell^{(|u|)} = 1$ results in the approximation of $\mathcal{D}^{|u|}f(\mathbf{x})$ of order $\mathcal{O}(||\mathbf{h}||_2^2)$.

Secondly, for $L > 1$, the constraints $\sum_{\ell=1}^L C_\ell^{(|u|)} \beta_\ell^{r+|u|} = \delta_{0,r} r = 0, 2, \dots, 2(L-1)$ allow to eliminate some higher-order terms so as to reach the order $\mathcal{O}(||\mathbf{h}||_2^{2L})$. One can also use $\sum_{\ell=1}^L C_\ell^{(|u|)} \beta_\ell^{r+|u|} = \delta_{0,r} r = -|u|, \dots, -|u| + L - 2, 0$ to increase the accuracy of approximations, but keeping the order $\mathcal{O}(||\mathbf{h}||_2^2)$ when $-|u| + L - 2 < 0$.

## Appendix B. Proof of Corollary 1

Let $\vec{q} = (q_1, \dots, q_d) \in \mathbb{N}^d$, $\vec{u} := (\mathbb{I}_u(1), \dots, \mathbb{I}_u(d)) \in \mathbb{N}^d$, and consider the set $\alpha := \{2\vec{q} + \vec{u} : ||\vec{q}||_1 = L\}$. As $f \in \mathcal{H}_{|u|+2L}$, the expansion of $f(\mathbf{x} + \mathbf{h}\mathbf{V})$ gives

$$f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{||\vec{i}||_1=0}^{|u|+2L-1} \mathcal{D}^{(\vec{i})} f(\mathbf{x}) \beta_{\ell}^{||\vec{i}||_1} \frac{(\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!} + \sum_{\substack{\vec{i} \notin \alpha \\ ||\vec{i}||_1=|u|+2L}} \mathcal{D}^{(\vec{i})} f(\mathbf{x}) \beta_{\ell}^{|u|+2L} \frac{(\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!} + R_{|u|}(\mathbf{h}, \beta_{\ell}, \mathbf{V}),$$

with the remainder term $R_{|u|}(\mathbf{h}, \beta_\ell, \mathbf{V}) := \beta_\ell^{|u|+2L} \sum_{\substack{\vec{i} \in \alpha \\ ||\vec{i}||_1=|u|+2L}} \mathcal{D}^{(\vec{i})} f(\mathbf{x} + \beta_\ell \mathbf{h} \mathbf{V}) \frac{(\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!}$. Thus, $R_{|u|}(\mathbf{h}, \beta_\ell, \mathbf{V}) := \beta_\ell^{|u|+2L} \sum_{\substack{\vec{i} \in \alpha \\ ||\vec{i}||_1=L}} \frac{\mathcal{D}^{(2\vec{q}+\vec{u})} f(\mathbf{x}+\beta_\ell \mathbf{h} \mathbf{V})}{(2\vec{q}+\vec{u})!} (\mathbf{h} \mathbf{V})^{2\vec{q}+\vec{u}}$ and

$$\begin{aligned} R_{|u|}(\mathbf{h}, \beta_{\ell}, \mathbf{V}) &= \beta_{\ell}^{|u|+2L} \sum_{||\vec{q}||_1=L} \frac{\mathcal{D}^{(2\vec{q}+\vec{u})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V})}{(2\vec{q} + \vec{u})!} (\mathbf{h} \mathbf{V})^{2\vec{q} + \vec{u}} \\ &= \beta_{\ell}^{|u|+2L} (\mathbf{h} \mathbf{V})^{\vec{u}} \sum_{||\vec{q}||_1=L} \frac{\mathcal{D}^{(2\vec{q}+\vec{u})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V})}{(2\vec{q} + \vec{u})!} (\mathbf{h}^2 \mathbf{V}^2)^{\vec{q}}. \end{aligned}$$

Using $E := \sum_{\ell=1}^L C_\ell^{(|u|)} \mathbb{E}[f(\mathbf{x} + \beta_\ell \mathbf{h} \mathbf{V}) \prod_{k \in u} \frac{V_k}{h_k \sigma^2}]$, Theorem 1 implies that the absolute value of the bias $B := |E - \mathcal{D}^{|u|}f(\mathbf{x})| = |\sum_{\ell=1}^L C_\ell^{(|u|)} \mathbb{E}[R_{|u|}(\mathbf{h}, \beta_\ell, \mathbf{V}) \prod_{k \in u} \frac{V_k}{h_k \sigma^2}]|$ is given by

$$B \le \sum_{\ell=1}^{L} |C_{\ell}^{(|u|)} \beta_{\ell}^{|u|+2L}| M_{|u|+2L} \mathbb{E} \left[ ||\mathbf{h}^2 \mathbf{V}^2||_1^L \prod_{k \in u} \frac{V_k^2}{\sigma^2} \right],$$

$$\text{as } \left| \sum_{||\vec{q}||_1=L} \frac{\mathcal{D}^{(2\vec{q}+\vec{u})} f(\mathbf{x}+\beta_\ell \mathbf{h} \mathbf{V})}{(2\vec{q}+\vec{u})!} (\mathbf{h}^2 \mathbf{V}^2)^{\vec{q}} \right| \le M_{|u|+2L} ||\mathbf{h}^2 \mathbf{V}^2||_1^L.$$

$$\text{Using } R_k = V_k/\sigma, \text{ the results hold because } \mathbb{E}[||\mathbf{h}^2 \mathbf{V}^2||_1^L \prod_{k \in u} \frac{V_k^2}{\sigma^2}] = \sigma^{2L} \mathbb{E}[||\mathbf{h}^2 \mathbf{R}^2||_1^L \prod_{k \in u} R_k^2].$$

For $V_k \sim \mathcal{U}(-\zeta, \zeta)$ with $\zeta > 0$ and $k = 1, \dots, d$, we have $\mathbb{E}[||\mathbf{h}^2 \mathbf{V}^2||_1^L \prod_{k \in u} \frac{V_k^2}{\sigma^2}] \le \zeta^{2L} ||\mathbf{h}^2||_1^L$.

## Appendix C. Proof of Corollary 2

Let $\alpha := \{\vec{q} + \vec{u} : ||\vec{q}||_1 = 1\}$. As $f \in \mathcal{H}_{|u|+1}^M$, we can write
$$f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{\|\tilde{\mathbf{v}}\|_1=0}^{|u|} \mathcal{D}^{(\tilde{\mathbf{v}})} f(\mathbf{x}) \beta_{\ell}^{\|\tilde{\mathbf{v}}\|_1} \frac{(\mathbf{h} \mathbf{V})^{\tilde{\mathbf{v}}}}{\tilde{\mathbf{v}}!} + \sum_{\substack{\|\tilde{\mathbf{v}}\|_1=|u|+1 \\ \tilde{\mathbf{v}} \notin \alpha}} \mathcal{D}^{(\tilde{\mathbf{v}})} f(\mathbf{x}) \beta_{\ell}^{\|\tilde{\mathbf{v}}\|_1+1} \frac{(\mathbf{h} \mathbf{V})^{\tilde{\mathbf{v}}}}{\tilde{\mathbf{v}}!} + R_{|u|}(\mathbf{h}, \beta_{\ell}, \mathbf{V}),$$

with the remainder term $R_{|u|}(\mathbf{h}, \beta_{\ell}, \mathbf{V}) := \beta_{\ell}^{|u|+1} \sum_{\tilde{\mathbf{v}} \in \alpha} \mathcal{D}^{(\tilde{\mathbf{v}})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{(\mathbf{h} \mathbf{V})^{\tilde{\mathbf{v}}}}{\tilde{\mathbf{v}}!}$. Using $E := \sum_{\ell=1}^L C_{\ell}^{(|u|)} \mathbb{E}[f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \prod_{k \in u} \frac{V_k}{(h_k \sigma^2)}]$ and Theorem 1, the results hold by analogy to the proof of Corollary 1. Indeed, if $R_k := V_k/\sigma$, then

$$B \leq M_{|u|+1} \mathbb{E} \left[ \left\| \mathbf{h} \mathbf{V} \right\|_1 \prod_{k \in u} \frac{V_k^2}{\sigma^2} \right] \Gamma_{|u|+1} = \sigma M_{|u|+1} \mathbb{E} \left[ \left\| \mathbf{h} \mathbf{R} \right\|_1 \prod_{k \in u} R_k^2 \right] \Gamma_{|u|+1} \\ \leq \sigma M_{|u|+1} \left\| \mathbf{h} \right\|_2 \mathbb{E} \left[ \left\| \mathbf{R} \right\|_2 \prod_{k \in u} R_k^2 \right] \Gamma_{|u|+1}.$$

For $V_k \sim \mathcal{U}(-\xi, \xi)$ with $\xi > 0$ and $k = 1, \dots, d$, we have

$$B \leq M_{|u|+1} \mathbb{E} \left[ \left\| \mathbf{h} \mathbf{V} \right\|_1 \prod_{k \in u} \frac{V_k^2}{\sigma^2} \right] \Gamma_{|u|+1} \leq M_{|u|+1} \xi \left\| \mathbf{h} \right\|_1 \mathbb{E} \left[ \prod_{k \in u} \frac{V_k^2}{\sigma^2} \right] \Gamma_{|u|+1} = M_{|u|+1} \xi \Gamma_{|u|+1}.$$

## Appendix D. Proof of Theorem 2

As $f \in \mathcal{H}_{|u|+1}$ implies that $f \in \mathcal{H}_{r^*+1}$ with $r^* \le |u| - 1$, we have

$$\left| f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) - \sum_{\|\tilde{\mathbf{v}}\|_1=0}^{r^*} \mathcal{D}^{(\tilde{\mathbf{v}})} f(\mathbf{x}) \beta_{\ell}^{\|\tilde{\mathbf{v}}\|_1} \frac{(\mathbf{h} \mathbf{V})^{\tilde{\mathbf{v}}}}{\tilde{\mathbf{v}}!} \right| \leq M_{r^*+1} \|\beta_{\ell} \mathbf{h} \mathbf{V}\|_2^{r^*+1},$$

Using the fact that $\sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^r = 0$ for $r = 0, 1, \dots, r^*$, we can write

$$\sum_{\ell=1}^{L} C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{\ell=1}^{L} C_{\ell}^{(|u|)} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) - \sum_{\|\tilde{\mathbf{v}}\|_1=0}^{r^*} \mathcal{D}^{(\tilde{\mathbf{v}})} f(\mathbf{x}) \beta_{\ell}^{\|\tilde{\mathbf{v}}\|_1} \frac{(\mathbf{h} \mathbf{V})^{\tilde{\mathbf{v}}}}{\tilde{\mathbf{v}}!} \right],$$

which leads to $\left| \sum_{\ell=1}^L C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \right| \le \sum_{\ell=1}^L \left| C_{\ell}^{(|u|)} \beta_{\ell}^{r^*+1} \right| M_{r^*+1} \|\mathbf{h} \mathbf{V}\|_2^{r^*+1}$.

By taking the variance of the proposed estimator, we have

$$\begin{aligned} \mathbb{V}[\widehat{\mathcal{D}^{|u|} f(\mathbf{x})}] &:= \frac{1}{N} \mathbb{V} \left[ \left\{ \sum_{\ell=1}^{L} C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \right\} \prod_{k \in u} \frac{V_k}{(h_k \sigma^2)} \right] \\ &\leq \frac{1}{N} \mathbb{E} \left[ \left\{ \sum_{\ell=1}^{L} C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \right\}^2 \prod_{k \in u} \frac{V_k^2}{(h_k^2 \sigma^4)} \right] \\ &\leq \frac{M_{r^*+1}^2 \left( \sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{r^*+1} \right)^2}{N \prod_{k \in u} h_k^2} \mathbb{E} \left[ \frac{\|\mathbf{h} \mathbf{V}\|_2^{2(r^*+1)}}{\sigma^2 |u|} \prod_{k \in u} \frac{V_k^2}{\sigma^2} \right] \\ &\leq \frac{M_{r^*+1}^2 \left( \sum_{\ell=1}^L C_{\ell}^{(|u|)} \beta_{\ell}^{r^*+1} \right)^2}{N \sigma^2 (|u|-r^*-1) \prod_{k \in u} h_k^2} \mathbb{E} \left[ \left\| \mathbf{h}^2 \right\|_2^{r^*+1} \mathbb{E} \left[ \left\| \mathbf{R}^2 \right\|_2^{r^*+1} \prod_{k \in u} R_k^2 \right] \right], \end{aligned}$$

where $R_k = V_k/\sigma$.

$$\text{If } V_k \sim \mathcal{U}(-\xi, \xi), \quad \mathbb{V}[\widehat{\mathcal{D}^{|u|} f(\mathbf{x})}] \le \frac{3^{|u|} M_{r^*+1}^2 \left( \sum_{\ell=1}^L \left| C_{\ell}^{(|u|)} \beta_{\ell}^{r^*+1} \right| \right)^2}{N \xi^{2(|u|-r^*-1)} \prod_{k \in u} h_k^2} \left\| \mathbf{h} \right\|_2^{2(r^*+1)}.$$

The results hold using Corollary 2 and the fact that MSE = $B^2 + \mathbb{V}[\widehat{\mathcal{D}^{|u|} f(\mathbf{x})}]$.
## Appendix E. Proof of Corollary 3

Let $\eta := |u| - r^* - 1 > 0$, $F_0 := 3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2 d^{r^*+1}$ and $F_1 := M_{|u|+1}^2 \Gamma_{|u|+1}^2 d^2$. By minimizing $\xi^2 M_{|u|+1}^2 \Gamma_{|u|+1}^2 d^2 h^2 + \frac{3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2 d^{r^*+1}}{N \xi^2 (|u|-r^*-1) h^2 (|u|-r^*-1)}$, we obtain

$$h_{op} = \frac{1}{\xi} \left( \frac{\eta F_0}{F_1} \right)^{\frac{1}{2\eta+2}} N^{-\frac{1}{2\eta+2}}; \quad \frac{F_0}{N \xi^2 \eta h^2 \eta} = \frac{F_0}{N N^{-\frac{2\eta}{2\eta+2}}} \left( \frac{F_1}{\eta F_0} \right)^{\frac{2\eta}{2\eta+2}} = \left( \frac{F_0}{N} \right)^{\frac{1}{\eta+1}} \left( \frac{F_1}{\eta} \right)^{\frac{\eta}{\eta+1}};$$

and

$$d^{\frac{r^*+1}{|u|-r^*}} d^{\frac{2(|u|-r^*-1)}{|u|-r^*}} = d^{\frac{2|u|-r^*-1}{|u|-r^*}} = d^{1+\frac{|u|-1}{|u|-r^*}}.$$

## Appendix F. Proof of Theorem 3

The first two results hold by combining the biases obtained in Corollary 1 and the upper bounds of the variance provided in Theorem 2.

For the last result, let $\eta := |u| - r^* - 1 > 0$, $F_0 := 3^{|u|} M_{r^*+1}^2 \Gamma_{r^*+1}^2 d^{r^*+1}$ and $F_1 := M_{|u|+2L'}^2 \Gamma_{|u|+2L'}^2 d^{2L'}$. By minimizing the last upper bound, we obtain $h_{op} := \frac{1}{\xi} \left( \frac{\eta F_0}{2L' F_1} \right)^{\frac{1}{2\eta+4L'}} N^{-\frac{1}{2\eta+4L'}}$ and

$$\frac{F_0}{N \xi^2 \eta h^2 \eta} = \frac{F_0}{N N^{-\frac{2\eta}{2\eta+4L'}}} \left( \frac{2L' F_1}{\eta F_0} \right)^{\frac{2\eta}{2\eta+4L'}} = \left( \frac{F_0}{N} \right)^{\frac{2L'}{\eta+2L'}} \left( \frac{2L' F_1}{\eta} \right)^{\frac{\eta}{\eta+2L'}} \propto N^{-\frac{2L'}{|u|-r^*-1+2L'}} d^{\frac{2L'|u|}{|u|-r^*-1+2L'}}.$$

## Appendix G. On Remark 6

The variance is $\mathbb{V}[\widehat{D^{|u|}f(x)}] = \frac{1}{N} \mathbb{V}[\sum_{\ell=1}^L C_{\ell}^{(|u|)} f(x + \beta_{\ell} h V) \prod_{k \in u} \frac{V_k}{h_k \sigma^2}]$ and

$$\mathbb{V}[\widehat{D^{|u|}f(x)}] \le \frac{1}{N} \mathbb{E} \left[ \left( \sum_{\ell=1}^L C_{\ell}^{(|u|)} f(x + \beta_{\ell} h V) \right)^2 \prod_{k \in u} \frac{V_k^2}{h_k^2 \sigma^4} \right] \le \frac{\|f\|_{\infty}^2 \sum_{\ell_1=1}^L \sum_{\ell_2=1}^L C_{\ell_1}^{(|u|)} C_{\ell_2}^{(|u|)}}{N h^2 |u| \sigma^2 |u|}.$$

## References

1. Robbins, H.; Monro, S. A Stochastic Approximation Method. Ann. Math. Stat. 1951, 22, 400–407. [CrossRef]

2. Fabian, V. Stochastic approximation. In Optimizing Methods in Statistics; Elsevier: Amsterdam, The Netherlands, 1971; pp. 439–470.

3. Nemirovsky, A.; Yudin, D. Problem Complexity and Method Efficiency in Optimization; Wiley & Sons: New York, NY, USA, 1983.

4. Polyak, B.; Tsybakov, A. Optimal accuracy orders of stochastic approximation algorithms. Probl. Peredachi Inf. 1990, 2, 45–53.

5. Cristea, M. On global implicit function theorem. J. Math. Anal. Appl. 2017, 456, 1290–1302. [CrossRef]

6. Lamboni, M. Derivative formulas and gradient of functions with non-independent variables. Axioms 2023, 12, 845. [CrossRef]

7. Morris, M.D.; Mitchell, T.J.; Ylvisaker, D. Bayesian design and analysis of computer experiments: Use of derivatives in surface prediction. Technometrics 1993, 35, 243–255. [CrossRef]

8. Solak, E.; Murray-Smith, R.; Leithead, W.; Leith, D.; Rasmussen, C. Derivative observations in Gaussian process models of dynamic systems. In Advances in Neural Information Processing Systems 15; MIT Press: Cambridge, MA, USA, 2002.

9. Hoeffding, W. A class of statistics with asymptotically normal distribution. Ann. Math. Stat. 1948, 19, 293–325. [CrossRef]

10. Efron, B.; Stein, C. The jacknife estimate of variance. Ann. Stat. 1981, 9, 586–596. [CrossRef]

11. Sobol, I.M. Sensitivity analysis for non-linear mathematical models. Math. Model. Comput. Exp. 1993, 1, 407–414.

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

13. Saltelli, A.; Chan, K.; Scott, E. Variance-Based Methods, Probability and Statistics; John Wiley & Sons: Hoboken, NJ, USA, 2000.

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

15. 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]

16. 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]

17. Conti, S.; O'Hagan, A. Bayesian emulation of complex multi-output and dynamic computer models. J. Stat. Plan. Inference 2010, 140, 640–651. [CrossRef]
18. Haylock, R.G.; O'Hagan, A.; Bernardo, J.M. On inference for outputs of computationally expensive algorithms with uncertainty on the inputs. In *Bayesian Statistics 5: Proceedings of the Fifth Valencia International Meeting*; Oxford Academic: Oxford, UK, 1996; Volume 5, pp. 629–638.

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

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

21. Wahba, G. An introduction to (smoothing spline) anova models in rkhs with examples in geographical data, medicine, atmospheric science and machine learning. *arXiv* **2004**, arXiv:math/0410419.

22. 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]

23. 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]

24. 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]

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

26. Lamboni, M. Derivative-based generalized sensitivity indices and Sobol' indices. *Math. Comput. Simul.* **2020**, *170*, 236–256.
[CrossRef]

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

28. Lamboni, M. Measuring inputs-outputs association for time-dependent hazard models under safety objectives using kernels. *Int. J. Uncertain. Quantif.* **2024**, 1–17. [CrossRef]

29. Lamboni, M. Kernel-based measures of association between inputs and outputs using ANOVA. *Sankhya A* **2024**. [CrossRef]

30. Russi, T.M. *Uncertainty Quantification with Experimental Data and Complex System Models*; Spring: Berlin/Heidelberg, Germany, 2010.

31. Constantine, P.; Dow, E.; Wang, S. Active subspace methods in theory and practice: Applications to kriging surfaces. *SIAM J. Sci. Comput.* **2014**, *36*, 1500–1524. [CrossRef]

32. Zahm, O.; Constantine, P.G.; Prieur, C.; Marzouk, Y.M. Gradient-based dimension reduction of multivariate vector-valued functions. *SIAM J. Sci. Comput.* **2020**, *42*, A534–A558. [CrossRef]

33. Kucherenko, S.; Shah, N.; Zaccheus, O. *Application of Active Subspaces for Model Reduction and Identification of Design Space*; Springer: Berlin/Heidelberg, Germany, 2024; pp. 412–418.

34. Kubicek, M.; Minisci, E.; Cisternino, M. High dimensional sensitivity analysis using surrogate modeling and high dimensional model representation. *Int. J. Uncertain. Quantif.* **2015**, *5*, 393–414. [CrossRef]

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

36. Bates, D.; Watts, D. Relative curvature measures of nonlinearity. *J. Royal Stat. Soc. Ser. B* **1980**, *42*, 1–25. [CrossRef]

37. Guidotti, E. calculus: High-dimensional numerical and symbolic calculus in R. *J. Stat. Softw.* **2022**, *104*, 1–37. [CrossRef]

38. 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]

39. 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]

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

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

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

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

44. Wang, Z.; Navon, I.M.; Le Dimet, F.-X.; Zou, X. The second order adjoint analysis: Theory and applications. *Meteorol. Atmos. Phys.* **1992**, *50*, 3–20. [CrossRef]

45. Agarwal, A.; Dekel, O.; Xiao, L. Optimal algorithms for online convex optimization with multi-point bandit feedback. In Proceedings of the 23rd Conference on Learning Theory, Haifa, Israel, 27–29 June 2010; pp. 28–40.

46. 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; Feldman, V., Rakhlin, A., Shamir, O., Eds.; Volume 49, pp. 257–283.

47. Akhavan, A.; Pontil, M.; Tsybakov, A.B. *Exploiting Higher Order Smoothness in Derivative-Free Optimization and Continuous Bandits*, NIPS'20; Curran Associates Inc.: Red Hook, NY, USA, 2020.

48. Lamboni, M. Optimal and efficient approximations of gradients of functions with nonindependent variables. *Axioms* **2024**, *13*, 426.
[CrossRef]

49. Patelli, E.; Pradlwarter, H. Monte Carlo gradient estimation in high dimensions. *Int. J. Numer. Methods Eng.* **2010**, *81*, 172–188.
[CrossRef]
50. Prashanth, L.; Bhatnagar, S.; Fu, M.; Marcus, S. Adaptive system optimization using random directions stochastic approximation. *IEEE Trans. Autom. Control.* **2016**, *62*, 2223-2238.

51. Agarwal, N.; Bullins, B.; Hazan, E. Second-order stochastic optimization for machine learning in linear time. *J. Mach. Learn. Res.* **2017**, *18*, 4148-4187.

52. Zhu, J.; Wang, L.; Spall, J.C. Efficient implementation of second-order stochastic approximation algorithms in high-dimensional problems. *IEEE Trans. Neural Netw. Learn. Syst.* **2020**, *31*, 3087-3099. [CrossRef] [PubMed]

53. Zhu, J. Hessian estimation via stein's identity in black-box problems. In *Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference, Online, 15-17 August 2022; Bruna, J., Hesthaven, J., Zdeborova, L., Eds.; Volume 145 of Proceedings of Machine Learning Research, PMLR*, pp. 1161-1178.

54. Erdogdu, M.A. Newton-stein method: A second order method for glms via stein's lemma. In *Advances in Neural Information Processing Systems*; Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., Garnett, R., Eds.; *Curran Associates, Inc.*: Red Hook, NY, USA, 2015; Volume 28.

55. Stein, C.; Diaconis, P.; Holmes, S.; Reinert, G. Use of exchangeable pairs in the analysis of simulations. *Lect.-Notes-Monogr. Ser.* **2004**, *46*, 1-26.

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

57. Strichartz, R. *A Guide to Distribution Theory and Fourier Transforms, Studies in Advanced Mathematics*; CRC Press: Boca, FL, USA, 1994.

58. Rawashdeh, E. A simple method for finding the inverse matrix of Vandermonde matrix. *Math. Vesn.* **2019**, *71*, 207-213.

59. Arafat, A.; El-Mikkawy, M. A fast novel recursive algorithm for computing the inverse of a generalized Vandermonde matrix. *Axioms* **2023**, *12*, 27. [CrossRef]

60. Morris, M. Factorial sampling plans for preliminary computational experiments. *Technometrics* **1991**, *33*, 161-174. [CrossRef]

61. Roustant, O.; Barthe, F.; Iooss, B. Poincaré inequalities on intervals-application to sensitivity analysis. *Electron. J. Stat.* **2017**, *11*, 3081-3119. [CrossRef]

62. Lamboni, M. Multivariate sensitivity analysis: Minimum variance unbiased estimators of the first-order and total-effect covariance matrices. *Reliab. Eng. Syst. Saf.* **2019**, *187*, 67-92. [CrossRef]

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

64. Dutang, C.; Savicky, P. *R Package*, version 1.13. Randtoolbox: Generating and Testing Random Numbers; The R Foundation: Vienna, Austria, 2013.

**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.