# Optimal and Efficient Approximations of Gradients of Functions with Nonindependent Variables

Matieyendou Lamboni

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

## To cite this version:

Matieyendou Lamboni. Optimal and Efficient Approximations of Gradients of Functions with Nonindependent Variables. Axioms, 2024, 13 (7), pp.426. 10.3390/axioms13070426. hal-05077907

HAL Id: hal-05077907

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

Submitted on 21 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 and Efficient Approximations of Gradients of Functions with Nonindependent Variables

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

**Abstract:** Gradients of smooth functions with nonindependent variables are relevant for exploring complex models and for the optimization of the functions subjected to constraints. In this paper, we investigate new and simple approximations and computations of such gradients by making use of independent, central, and symmetric variables. Such approximations are well suited for applications in which the computations of the gradients are too expansive or impossible. The derived upper bounds of the biases of our approximations do not suffer from the curse of dimensionality for any 2-smooth function, and they theoretically improve the known results. Also, our estimators of such gradients reach the optimal (mean squared error) rates of convergence (i.e., $O(N^{-1})$) for the same class of functions. Numerical comparisons based on a test case and a high-dimensional PDE model show the efficiency of our approach.

**Keywords:** dependent variables; gradients; high-dimensional models; optimal estimators; tensor metric of nonindependent variables

**MSC:** 26A24; 60H25; 62Gxx; 49Qxx

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

check for
updates

**Citation:** Lamboni, M. Optimal and Efficient Approximations of Gradients of Functions with Nonindependent Variables. *Axioms* **2024**, *13*, 426.
https://doi.org/10.3390/axioms13070426

Academic Editor: Delfim F. M.
Torres

Received: 12 May 2024

Revised: 16 June 2024

Accepted: 19 June 2024

Published: 25 June 2024

<figure><img src="image_4.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

Nonindependent variables arise when at least two variables do not vary independently, and such variables are often characterized by their covariance matrices, distribution functions, copulas, and weighted distributions (see, e.g., [1–7]). Recently, dependency models provide explicit functions that link these variables together by means of additional independent variables [8–12]. Models with nonindependent input variables, including functions subjected to constraints, are widely encountered in different scientific fields, such as data analysis, quantitative risk analysis, and uncertainty quantification (see, e.g., [13–15]).

Analyzing such functions requires being able to calculate or to compute their dependent gradients, that is, the gradients that account for the dependencies among the inputs. Recall that gradients are involved in (i) inverse problems and optimization (see, e.g., [16–20]), (ii) exploring complex mathematical models or simulators (see [21–28] for independent inputs and [9,15] for nonindependent variables); (iii) Poincaré inequalities and equalities [9,28–30], and recently in (iv) the derivative-based ANOVA (i.e., exact expansions) of functions [28]. While the first-order derivatives of functions with nonindependent variables have been derived in [9] for screening dependent inputs of high-dimensional models, the theoretical expressions of the gradients of such functions (dependent gradients) have been introduced in [15], thereby enhancing the difference between the gradients and the first-order partial derivatives when the input variables are dependent or correlated.

In high-dimensional settings and for time-demanding models, having an efficient approach for computing the dependent gradients provided in [15] using a few model evaluations is worth investigating. So far, the adjoint methods can provide the exact classical gradients for some classes of PDE/ODE-based models [31–36]. Additionally,

<footnote n="1">¹ Department DFR-ST, University of Guyane, 97346 Cayenne, France; matieyendou.lamboni@univ-guyane.fr or 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>
Richardson's extrapolation and its generalization considered in [37] provide accurate estimates of the classical gradients using a number of model runs that strongly depends on the dimensionality. In contrast, the Monte Carlo approach allows for computing the classical gradients using a number of model runs that can be very less than the dimensionality (i.e., $d \in \mathbb{N}$) [17,38,39]. The Monte Carlo approach is a consequence of the Stokes theorem, which claims that the expectation of a function evaluated at a random point about $\mathbf{x} \in \mathbb{R}^d$ is the gradient of a certain function. Such a property leads to randomized approximations of the classical gradients in derivative-free optimization or zero-order stochastic optimization (see [16,18–20] and the references therein). Such approximations are also relevant for applications in which the computations of the gradients are impossible [20].

Most of the randomized approximations of the classical gradients, including the Monte Carlo approach, rely on randomized kernels and/or random vectors that are uniformly distributed on the unit ball. The qualities of such approximations are often assessed by the upper bounds of the biases and the rates of convergence. The upper bounds provided in [19,20,40] depend on the dimensionality in general.

In this paper, we propose new surrogates of the gradients of smooth functions with nonindependent inputs and the associated estimators that comply with the following requirements:

* They are simple and applicable to a wide class of functions by making use of model evaluations at randomized points, which are only based on independent, central, and symmetric variables;

* They lead to a dimension-free upper bound of the bias, and they improve the best known upper bounds of the bias for the classical gradients;

* They lead to the optimal and parametric (mean squared error) rates of convergence;

* They are going to increase the computational efficiency and accuracy of the gradients estimates by means of a set of constraints.

The surrogates of the dependent gradients are derived in Section 3 by combining the properties of (i) the generalized Richardson extrapolation approach thanks to a set of constraints and (ii) the Monte Carlo approach based only on independent random variables that are symmetrically distributed about zero. Such expressions are followed by their order of approximations, biases, and a comparison with known results for the classical gradients. We also provide the estimators of such surrogates and their associated mean squared errors, including the rates of convergence for a wide class of functions (see Section 3.3). A number of numerical comparisons is considered so as to assess the efficiency of our approach. Section 4 presents comparisons of our approach to other methods and simulations based on a high-dimensional PDE (spatiotemporal) model with given autocollaborations among the initial conditions, which are considered in Section 5 to compare our approach to the adjoint-based methods. We conclude this work in Section 6.

## 2. Preliminaries

For an integer $d > 0$, let $\mathbf{X} := (X_1, \dots, X_d)$ be a random vector of continuous and nonindependent variables having $F$ as the joint cumulative distribution function (CDF) (i.e., $\mathbf{X} \sim F$). For any $j \in \{1, \dots, d\}$, we use $F_{X_j}$ or $F_j$ for the marginal CDF of $X_j$ and $F_j^{-1}$ for its inverse. Also, we use $(\sim j) := (1, \dots, j-1, j+1, \dots, d)$ and $\mathbf{X}_{\sim j} := (X_1, \dots, X_{j-1}, X_{j+1}, \dots, X_d)$. The equality (in distribution) $X \stackrel{d}{=} Z$ means that $X$ and $Z$ have the same CDF.

As the sample values of $\mathbf{X}$ are dependent, here we use $\frac{\partial f}{\partial x_k}$ for the formal partial derivative of $f$ with respect to $x_k$, that is, the partial derivative obtained by considering other inputs as constant or independent of $x_k$. Thus, $\nabla f := \left[ \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_d} \right]^T$ stands for the formal or classical gradient of $f$.
Given an open set $\Omega \subseteq \mathbb{R}^d$, consider a weak partial differentiable function $f: \Omega \to \mathbb{R}$ [41,42]. Given $\vec{i} := (i_1, \dots, i_d) \in \mathbb{N}^d$, denote $\mathcal{D}(\vec{i})f := (\prod_{k=1}^d \frac{\partial^{i_k}}{\partial x_k})f; (\mathbf{x})\vec{i} = \mathbf{x}\vec{i} := \prod_{k=1}^d x_k^{i_k}$, $\vec{i}! = 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}^n : \left| f(\mathbf{x}) - \sum_{0 \le i_1 + \dots + i_d \le \alpha - 1} \frac{\mathcal{D}(\vec{i}) f(\mathbf{y})}{\vec{i}!} (\mathbf{x} - \mathbf{y})^{\vec{i}} \right| \le M_{\alpha} ||\mathbf{x} - \mathbf{y}||_2^{\alpha} \right\}, $$

with $\alpha \ge 1$ and $M_\alpha > 0$. 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.

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

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

Also, denote $\mathbf{h}\mathbf{V} := (h_1V_1; \dots, h_dV_d)$; $\mathbf{h}^{-1}\mathbf{V} := (V_1/h_1; \dots, V_d/h_d)$ and $\beta_\ell \mathbf{h}\mathbf{V} := (\beta_\ell h_1V_1; \dots, \beta_\ell h_dV_d)$. The real $\beta_\ell$'s are used for controlling the order of approximations and the order of derivatives (i.e., $||\vec{i}||_1 = 1, 2$) that we are interested in. Finally, $h_j$'s are used to define a neighborhood of a sample point of $\mathbf{X}$ (i.e., $\mathbf{x}$). Thus, using $\beta_{\max} := \max(|\beta_1|, \dots, |\beta_L|)$ and keeping in mind the variance of $\beta_\ell h_j V_j$, we assume that $\forall j \in \{1, \dots, d\}$,

**Assumption 1.** $\beta_{\max}h_j\sigma \le 1/2$ or equivalently $0 < \beta_{\max}h_j|V_j| \le 1$ for bounded $V_j$'s.

### 3. Main Results

This section aims at providing new expressions of the gradient of a function with nonindependent variables and the associated order of approximations. We are also going to derive the estimators of such a gradient, including the optimal and parametric rates of convergence. Recall that the input variables are said to be nonindependent whenever there exist at least two variables $X_j, X_k$ such that the joint CDF $F_{j,k}(x_j, x_k) \ne F_j(x_j)F_k(x_k)$.

#### 3.1. Stochastic Expressions of the Gradients of Functions with Dependent Variables

Using the fact that $\mathbf{X} \sim F$, with $F(\mathbf{x}) \ne \prod_{j=1}^d F_j(x_j)$, we are able to model $\mathbf{X}$ as follows [8–12,14,43]:

$$ \begin{aligned} \mathbf{X}_{\sim j} &\stackrel{d}{=} r_j(\mathbf{X}_j, \mathbf{Z}_{\sim j}) \\ &= [r_{1,j}(\mathbf{X}_j, \mathbf{Z}_{\sim j}), \dots, r_{j-1,j}(\mathbf{X}_j, \mathbf{Z}_{\sim j}), r_{j+1,j}(\mathbf{X}_j, \mathbf{Z}_{\sim j}), \dots, r_{d,j}(\mathbf{X}_j, \mathbf{Z}_{\sim j})]^T, \end{aligned} \quad (1) $$

where $r_j : \mathbb{R}^d \to \mathbb{R}^{d-1}$; $X_j$; and $\mathbf{Z}_{\sim j} := (Z_1, \dots, Z_{j-1}, Z_{j+1}, \dots, Z_d)$ are independent. Moreover, we have $(X_j, \mathbf{X}_{\sim j}) \stackrel{d}{=} (X_j, r_j(X_j, \mathbf{Z}))$, and it is worth noting that the function $r_j$ is invertible with respect to $\mathbf{Z}_{\sim j}$ for continuous variables, that is,

$$ \mathbf{Z}_{\sim j} = r_j^{-1}(\mathbf{X}_{\sim j} | X_j). $$

Note that the formal Jacobian matrix of $g: \mathbb{R}^d \to \mathbb{R}^d$, $\mathbf{x} \mapsto \mathbf{x}$ is the identity matrix. As $\mathbf{x}$ is a sample value of $\mathbf{X}$, the dependent Jacobian of $g$ based on the above dependency function is clearly not the identity matrix due to the fact that such a matrix accounts for the
dependencies among the elements of $\mathbf{x}$. The dependent partial derivatives of $\mathbf{x}$ with respect to $x_j$ are then given by [9,15]

$$J^{(j)}(\mathbf{x}) := \frac{\partial \mathbf{x}}{\partial x_j} = \left[ \frac{\partial r_{1,j}}{\partial x_j} \cdots \underbrace{1}_{\text{jth position}} \cdots \frac{\partial r_{d,j}}{\partial x_j} \right]^T (x_j, r_j^{-1}(\mathbf{x}_{\sim j} | x_j)),$$

and the dependent Jacobian matrix becomes (see [15] for more details)

$$J^d(\mathbf{x}) := [J^{(1)}(\mathbf{x}), \dots, J^{(d)}(\mathbf{x})].$$

Moreover, the gradient of $f$ with nonindependent variables is given by [15]

$$\mathrm{grad}(f)(\mathbf{x}) := [J^d(\mathbf{x})^T J^d(\mathbf{x})]^{-1} \nabla f(\mathbf{x}) = G^{-1}(\mathbf{x}) \nabla f(\mathbf{x}), \quad (2)$$

with $G(\mathbf{x}) := J^d(\mathbf{x})^T J^d(\mathbf{x})$ being the tensor metric and $G^{-1}(\mathbf{x})$ being its generalized inverse. Based on the above framework, Theorem 1 provides the stochastic expression of $\mathrm{grad}(f)(\mathbf{x})$. In what follows, denote $\mathbb{I}_{\bullet} := [1, \dots, 1]^T \in \mathbb{R}^d$.

**Theorem 1.** Assume that $f \in \mathcal{H}_\alpha$, with $\alpha \ge 2L$, Assumption 1 holds and that the $\beta_\ell$s are distinct. Then, there exists $\alpha_1 \in \{1, \dots, L\}$ and reals coefficients $C_1, \dots, C_L$ such that

$$\mathrm{grad}(f)(\mathbf{x}) = G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L} C_\ell \mathcal{E} \left[ f(\mathbf{x} + \beta_\ell \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] + \mathcal{O} \left( \|\mathbf{h}\|_2^{2\alpha_1} \right) \mathbb{I}_{\bullet}. \quad (3)$$

**Proof.** See Appendix A for the detailed proof. $\square$

Using the Kronecker symbol $\delta_{1,r}$, the setting $L = 1, \beta_1 = 1, C_1 = 1$, or the constraints $\sum_{\ell=1}^{L=2} C_\ell \beta_\ell^r = \delta_{1,r}$; $r = 0, 1$ lead to the order of approximation $\mathcal{O}(\|\mathbf{h}\|_2^2)$, while the constraints $\sum_{\ell=1}^{L} C_\ell \beta_\ell^r = \delta_{1,r}$; $r = 1, 3, 5, \dots, 2L - 1$ allow for increasing that order up to $\mathcal{O}(\|\mathbf{h}\|_2^{2L})$. For distinct $\beta$s, the above constraints lead to the existence of the constants $C_1, \dots, C_L$. Indeed, some constraints rely on the Vandermonde matrix of the form

$$A_L := \begin{bmatrix} 1 & 1 & \cdots & 1 \\ \beta_1 & \beta_2 & \cdots & \beta_L \\ \beta_1^2 & \beta_2^2 & \cdots & \beta_L^2 \\ \vdots & \vdots & \ddots & \vdots \end{bmatrix},$$

which is invertible for distinct values of the $\beta_\ell$s (i.e., $\beta_{\ell_1} \neq \beta_{\ell_2}$), because the determinant $\det(A_L) = \prod_{1 \le \ell_1 < \ell_2 \le L} (\beta_{\ell_1} - \beta_{\ell_2})$.

**Remark 1.** For an even integer $L$, 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, there are other possibilities, provided that $\sum_{\ell=1}^{L} C_\ell \beta_\ell = 1$.

Beyond the strong assumption made on the functions in Theorem 1, and knowing that increasing $L$ will require more evaluations of $f$ at random points, we are going to derive the upper bounds of the biases of our appropriations under different structural assumptions on the deterministic functions $f$ and $V$, such as $f \in \mathcal{H}_\alpha$, with $\alpha > 1$. To that end, denote $\mathbf{R} := (R_1, \dots, R_d)$ as 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 let $\mathcal{R}_c$ be the set of such random vectors. Define

$$K_1 := \inf_{\mathbf{R} \in \mathcal{R}_c} \left\| \left\| G^{-1}(\mathbf{x}) \left| \mathbb{E}[\mathbf{R}^2 || \mathbf{R} ||_2] \right| \right\|_1 \right\};$$

$$K_2 := \inf_{\mathbf{R} \in \mathcal{R}_c} \left\| \left\| G^{-1}(\mathbf{x}) \left| \mathbb{E}[\mathbf{R}^2 || \mathbf{R} ||_2] \right| \right\|_2 \right\};$$

with $|G^{-1}|$ the matrix obtained by putting the entries of $G^{-1}$ in the absolute value. When $1 < \alpha \le 2$, only $L = 1$ or $L = 2$ can be considered for any function that belongs to $\mathcal{H}_\alpha$. To be able to derive the parametric rates of convergence, Corollary 1 starts providing the upper bounds of the bias when $L = 2$.

**Corollary 1.** Consider $\beta_1 = 1, \beta_2 = -1; C_1 = 1/2; \text{ and } C_2 = -1/2$. If $f \in \mathcal{H}_2$ and Assumption 1 hold, then there exists $M_2 > 0$ such that

$$\left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_1 \le \sigma M_2 K_1 ||\mathbf{h}||_2; \quad (4)$$

$$\left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \le \sigma M_2 K_2 ||\mathbf{h}||_2. \quad (5)$$

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

For a particular choice of $\mathbf{V}$, we obtain the results below.

**Corollary 2.** Consider $\beta_1 = 1, \beta_2 = -1; C_1 = 1/2; C_2 = -1/2$. If $V_k \sim \mathcal{U}(-\xi, \xi)$, where $\xi > 0, k = 1, \dots, d, f \in \mathcal{H}_2$, and Assumption 1 hold, then

$$\left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_1 \le M_2 \xi \left\| G^{-1}(\mathbf{x}) \|\mathbf{I}_{\bullet}\|_1 \|\mathbf{h}\|_1 \right\|_1; \quad (6)$$

$$\left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \le M_2 \xi \left\| G^{-1}(\mathbf{x}) \|\mathbf{I}_{\bullet}\|_2 \|\mathbf{h}\|_1 \right\|_2. \quad (7)$$

**Proof.** Since $|V_k| \le \xi$, we have $\|\mathbf{h}\mathbf{V}\|_1 \le \xi \|\mathbf{h}\|_1$, and the results hold using the following upper bounds: $M_2 \left\| G^{-1}(\mathbf{x}) \left| \mathbb{E}\left[ \frac{\mathbf{V}^2}{\sigma^2} ||\mathbf{h}\mathbf{V}||_1 \right] \right| \right\|_1$ and $M_2 \left\| G^{-1}(\mathbf{x}) \left| \mathbb{E}\left[ \frac{\mathbf{V}^2}{\sigma^2} ||\mathbf{h}\mathbf{V}||_1 \right] \right| \right\|_2$ obtained in Appendix B. □

It is worth noting that choosing $h_k = h$ and $\xi = 1/d^2$ leads to the dimension-free upper bound of the bias, that is,

$$\left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_1 \le \frac{M_2 h}{d} \left\| G^{-1}(\mathbf{x}) \|\mathbf{I}_{\bullet}\|_1 \right\|_1,$$

because $\|G^{-1}(\mathbf{x})\|\mathbf{I}_{\bullet}\|_2$ is a function of $d$ in general.

For the sequel of generality, Corollary 3 provides the bias of our approximations for highly smooth functions. To that end, define

$$K_{2,L} := \inf_{\mathbf{R} \in \mathcal{R}_c} \left\| \left\| G^{-1}(\mathbf{x}) \left| \mathbb{E}[\mathbf{R}^2 || \mathbf{R} ||_2^L] \right| \right\|_2 \right\};$$

$$K_3 := \sum_{\ell=1}^{L+1} \left| C_{\ell} \beta_{\ell}^{1+L} \right|.$$
**Corollary 3.** For an odd integer $L > 2$, consider $\sum_{\ell=1}^{L+1} C_{\ell} \beta_{\ell}^r = \delta_{1,r}$; $r = 0, 1, \dots, L$. If $f \in \mathcal{H}_{1+L}$ and Assumption 1 hold, then there exists $M_{1+L} > 0$ such that

$$ \left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L+1} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \le \sigma^L M_{1+L} K_{2,L} K_3 \|\mathbf{h}\|_2^L. $$

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

$$ \left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L+1} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \le \xi^L M_{1+L} \left\| G^{-1}(\mathbf{x}) \left| \mathbf{I}_{\bullet} \right| \right\|_2 \|\mathbf{h}\|_1^L K_3. $$

**Proof.** The proofs are similar to those of Corollary 1 (see Appendix B). □

In view of the results provided in Corollary 3, finding the $\beta$s and $C$s that minimize the quantity $K_3 = \sum_{\ell=1}^{L+1} |C_{\ell} \beta_{\ell}^{1+L}|$ might be helpful for improving the above upper bounds.

## 3.2. Links to Other Works for Independent Input Variables

Recall that for independent input variables, the matrix $|G^{-1}(\mathbf{x})|$ comes down to the identity matrix, and $\mathrm{grad}(f) = \nabla f$. Thus, Equation (7) becomes

$$ \left\| \nabla f(\mathbf{x}) - \sum_{\ell=1}^{L=2} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \le M_2 h, $$

when $\xi = \sqrt{d}/d^2$. Taking $\xi = \sqrt{d}/d$ leads to the upper bound $M_2hd$.

Other results about the upper bounds of the bias of the (formal) gradient approximations have been provided in [19,20] (and the references therein) under the same assumptions made on $f$ and the evaluations of $f$. Such results rely on a random vector $\mathbf{S}$ that is uniformly distributed on the unit ball and a kernel $K$. Under such a framework, the upper bound derived in [19,20] is

$$ \left\| \nabla f(\mathbf{x}) - \frac{d}{h} \mathbb{E}[f(\mathbf{x} + U h \mathbf{S}) \mathbf{S} K(U)] \right\|_2 \le 2\sqrt{2\alpha} d M_{\alpha} h^{\alpha-1}, $$

where $U \sim \mathcal{U}(-1, 1)$ is independent of $\mathbf{S}$. Therefore, our results improve the upper bound obtained in [19,20] when $\alpha = 2$, for instance.

## 3.3. Computation of the Gradients of Functions with Dependent Variables

Consider a sample of $\mathbf{V}$ given by $\{\mathbf{V}_i := (V_{i,1}, \dots, V_{i,d})\}_{i=1}^N$. Using Equation (3), the estimator of $\mathrm{grad}(f)(\mathbf{x})$ is derived as follows:

$$ \widehat{\mathrm{grad}}(f)(\mathbf{x}) := G^{-1}(\mathbf{x}) \frac{1}{N} \sum_{i=1}^{N} \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}_i) \frac{\mathbf{V}_i \mathbf{h}^{-1}}{\sigma^2}. $$

To assess the quality of such an estimator, it is common to use the mean squared error (MSE), including the rates of convergence. The MSEs are often used in statistics for determining the optimal value of $\mathbf{h}$ as well. Theorem 2 and Corollary 4 provide such quantities of interest. To that end, define

$$ K_4 := \inf_{\mathbf{R} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \mathbf{R} \mathbf{h}^{-1} \right\|_2^2 \left\| \mathbf{R}^2 \right\|_2 \right]. $$
**Theorem 2.** Consider $\beta_1 = 1, \beta_2 = -1$; $C_1 = 1/2$; and $C_2 = -1/2$. If $f \in \mathcal{H}_2$ and Assumption 1 hold, then

$$ \mathbb{E} \left[ \left\| \widehat{\mathrm{grad}}(f)(\mathbf{x}) - \mathrm{grad}(f)(\mathbf{x}) \right\|_2^2 \right] \le \sigma^2 M_2^2 K_2^2 \|\mathbf{h}\|_2^2 + \frac{M_1^2 K_4 \|\mathbf{h}^2\|_2}{N}. \quad (8) $$

Moreover, if $V_k \sim U(-\zeta, \zeta)$, with $\zeta > 0$, $k = 1, \dots, d$ and $\mathbf{R}_0 := \mathbf{V}/\sigma$, then

$$ \mathbb{E} \left[ \left\| \widehat{\mathrm{grad}}(f)(\mathbf{x}) - \mathrm{grad}(f)(\mathbf{x}) \right\|_2^2 \right] \le M_2^2 \zeta^2 \left\| G^{-1}(\mathbf{x}) \mathbb{I}_{\bullet} \right\|_2^2 \|\mathbf{h}\|_1^2 + \frac{3\sqrt{d} M_1^2 \|\mathbf{h}^2\|_2}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \mathbf{R}_0 \mathbf{h}^{-1} \right\|_2^2 \right]. \quad (9) $$

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

Using a uniform bandwidth, that is, $h_k = h$, with $k = 1, \dots, d$, the upper bounds of MSEs provided in Theorem 2 have simple expressions. Indeed, the upper bounds in Equations (8) and (9) become, respectively,

$$ \sigma^2 M_2^2 K_2^2 dh^2 + \frac{M_1^2 \sqrt{d}}{N} \inf_{\mathbf{R} \in \mathcal{R}_c} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \mathbf{R} \right\|_2^2 \left\| \mathbf{R}^2 \right\|_2 \right]; $$

$$ M_2^2 \zeta^2 \left\| G^{-1}(\mathbf{x}) \mathbb{I}_{\bullet} \right\|_2^2 d^2 h^2 + \frac{3d M_1^2}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \mathbf{R}_0 \right\|_2^2 \right]. $$

It comes out that the second terms of the above upper bounds do not depend on the bandwidth $h$. This key observation leads to the derivation of the optimal and parametric rates of convergence of the proposed estimator.

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

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

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

It is worth noting that the upper bound of the squared bias obtained in Corollary 4 does not depend on the dimensionality, thanks to the choice of $\zeta = d^{-3/2}$. But, the derived rate of convergence depends on $d^2$, thus meaning that our estimator suffers from the curse of dimensionality. In higher dimensions, an attempt to improve our results consists of controlling the upper bound of the second-order moment of the estimator through $\sum_{\ell=1}^L |C_\ell \beta_\ell|$.

**Remark 2.** For highly smooth functions (i.e., $f \in \mathcal{H}_{1+L}$, with $L > 3$) and under the assumptions made in Corollary 3, we can check that (see Appendix C)

$$ \mathbb{E} \left[ \left\| \widehat{\mathrm{grad}}(f)(\mathbf{x}) - \mathrm{grad}(f)(\mathbf{x}) \right\|_2^2 \right] \le \zeta^{2L} M_{1+L}^2 \left\| G^{-1}(\mathbf{x}) \mathbb{I}_{\bullet} \right\|_2^2 \|\mathbf{h}\|_1^{2L} K_3^2 + \frac{3\sqrt{d} M_1^2 \|\mathbf{h}^2\|_2}{N} \left( \sum_{\ell=1}^{L+1} |C_\ell \beta_\ell| \right)^2 \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \mathbf{R}_0 \mathbf{h}^{-1} \right\|_2^2 \right]. $$
### 4. Computations of the Formal Gradient of Rosenbrock's Function

For comparing our approach to (i) the finite differences method (FDM) using the R package numDeriv [44], with $h = 10^{-4}$, and (ii) the Monte Carlo (MC) approach provided in [17], with $h = 10^{-4}$, let us consider the Rosenbrock function given as follows: $\forall \mathbf{x} \in \mathbb{R}^d$,

$$r(\mathbf{x}) := \sum_{k=1}^{d-1} \left[ (1 - x_k)^2 + 100 (x_{k+1} - x_k^2)^2 \right].$$

The gradient of that function at $0$ is $\nabla r(0) = [-2, \dots, -2, 0]^T \in \mathbb{R}^{100}$ (see [17]). To assess the numerical accuracy of each approach, the following measure is considered:

$$Err := \frac{\| \nabla r(0) - \widehat{\nabla} r(0) \|_1}{\| \nabla r(0) \|_1},$$

where $\widehat{\nabla} r(0)$ is the estimated value of the gradient. Table 1 reports the values of Err for the three approaches. To obtain the results using our approach, we have used $h = 1/\sqrt{N}$, with $N$ as the sample size and $\zeta = 1/d^2 = 10^{-4}$, with $d = 100$. Also, the Sobol sequence has been used for generating the values of the $V_j$s, and the Gram–Schmidt algorithm is applied to obtain (perfect) orthogonal vectors for a given $N$.

<table><caption>**Table 1.** Values of Err for three different approximations of the formal gradients.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="6">Number of Total Model Evaluations (i.e., LN)</th></tr><tr><th>100</th><th>150</th><th>200</th><th>200</th><th>1000</th><th>1000</th></tr></thead><tbody><tr><td>FDM [44]</td><td>-</td><td>-</td><td>-</td><td>0.005</td><td>-</td><td>-</td></tr><tr><td>MC [17]</td><td>0.042</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="2">This paper</td><td>L = 1</td><td>L = 1</td><td>L = 1</td><td>L = 2</td><td>L = 1</td><td>L = 2</td></tr><tr><td>0.035</td><td>0.014</td><td>0.009</td><td>0.009</td><td>0.0020</td><td>0.00199</td></tr></tbody></table>

Based on Table 1, our approach provides efficient results compared to other methods. Since the FDM is not possible when $N < 2d = 200$, it comes out that our approach is quite flexible thanks to $L$ and the fact that the gradient can be computed for every value of $N$. Increasing $N$ improved our results, as expected.

### 5. Application to a Heat PDE Model with Stochastic Initial Conditions

#### 5.1. Heat Diffusion Model and Its Formal Gradient

Consider a time-dependent model $f(x, t)$ 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], \quad 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 \in \mathbb{R}_+$ represents the diffusion coefficient. It is common to consider $J(Z(x)) := \frac{1}{2} \int_0^T \int_0^{10} (f(x, t))^2 dxdt$ as the quantity of interest (QoI). The spatial discretisation consists in subdividing the spatial domain $[0, 1]$ in $d$ equally sized cells, which leads to $d$ initial conditions or inputs given by $Z(x_j)$ with $j = 1, \dots, d$. Given zero-mean random variables $(R_j, j = 1, \dots, d)$, assume that $X_j := Z(x_j) = \sin(2\pi x_j) + s_j R_j, j = 1, \dots, d$, where $s_j \in \mathbb{R}_+$ represents the inverse precision about our knowledge on the initial conditions. For the dynamic aspect, a time step of 0.025 is considered starting from 0 up to $T = 5$.
Given a direction $z(x)$ and the Gâteaux derivative $\check{f}(x, t) := \frac{\partial f}{\partial z(x)}$, the tangent linear model is derived as follows:

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

and we can check that the adjoint model (AM) (i.e., $f^a$) is given by

$$ \left\{ \begin{array}{ll} -\frac{\partial f^a}{\partial t} - D \frac{\partial^2 f^a}{\partial x^2} = 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 formal gradient of $J(Z(x))$ with respect to the inputs $Z(x)$ is $\nabla_Z J(Z(x)) = f^a(x, 0)$. Remark that the above gradient relies on $f^a(x, 0)$, and only one evaluation of such a function is needed.

## 5.2. Spatial Autocorrelations of Initial Conditions and the Tensor Metric

Recall that the above gradient is based on the assumption of independent input variables, thus suggesting that the initial conditions within different cells are uncorrelated. To account for the spatial autocorrelations between different cells, assume that the $d$ input variables follow the Gaussian process with the following autocorrelation function:

$$ \rho(X_{j_1}, X_{j_2}) = \left(\frac{1}{2}\right)^{|j_1 - j_2|} \mathbb{I}_{[0,3]}(|j_1 - j_2|); \quad \forall j_1, j_2 \in \{1, \dots, d\}, $$

where $\mathbb{I}_{[0,3]}(|j_1 - j_2|) = 1$ if $|j_1 - j_2| \in [0,3]$ and is zero otherwise. Such spatial autocorrelations lead to the correlation matrix of the form

$$ \mathcal{R} := \begin{bmatrix} 1 & 0.5 & 0.25 & 0.125 & 0 & 0 & 0 & \dots & 0 \\ 0.5 & 1 & 0.5 & 0.25 & 0.125 & 0 & 0 & \dots & 0 \\ 0.25 & 0.5 & 1 & 0.5 & 0.25 & 0.125 & 0 & \dots & 0 \\ 0.125 & 0.25 & 0.5 & 1 & 0.5 & 0.25 & 0.125 & 0 & \dots 0 \\ 0 & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \end{bmatrix}. $$

Using the same standard deviation $s_j = s$ leads to the following covariance matrix $\Sigma = s^2 \mathcal{R}$, and $\mathbf{X} = (X_1, \dots, X_d) \sim \mathcal{N}_d(\mu, \Sigma)$, with $\mu := (\sin(2\pi c_1), \dots, \sin(2\pi c_d))$, and $c_1, \dots, c_d$ being the centers of the cells. The associated dependency model is given below.

Consider the diagonal matrix $D_{\sim j} = \text{diag}(\Sigma_{1,1}, \dots, \Sigma_{j-1,j-1}, \Sigma_{j+1,j+1}, \dots, \Sigma_{d,d})$ and the Gaussian random vector $\mathbf{W} \sim \mathcal{N}_{d-1}(\mu_{\sim j}, D_{\sim j})$. Denote with $\Sigma^{(j)}$ the matrix obtained by moving the $j^{th}$ row and column of $\Sigma$ to the first row and column; $\mathcal{L}^{(j)}$ is the Cholesky factor of $\Sigma^{(j)}$, and $\mu^{(j)} := (\mu_j, \mu_1, \dots, \mu_{j-1}, \mu_{j+1}, \dots, \mu_d)$. We can see that $(X_j, \mathbf{X}_{\sim j}) \sim \mathcal{N}_d(\mu^{(j)}, \Sigma^{(j)})$, and the dependency model is given by [10]

$$ (X_j, \mathbf{X}_{\sim j}) = \mathcal{L}^{(j)} \left[ \frac{1}{\sqrt{\Sigma_{j,j}}} (X_j - \mathbb{E}[X_j]) \middle/ D_{\sim j}^{-1/2} (\mathbf{W} - \mu_{\sim j}) \right] + \mu^{(j)}; \quad j = 1, \dots, d. \quad (10) $$

Based on Equation (10), we have $\frac{\partial \mathbf{X}_{\sim j}}{\partial x_j} = \frac{\mathcal{L}^{(j)}}{\sqrt{\Sigma_{j,j}}} = \frac{\Sigma_{\sim 1,1}}{\Sigma_{j,j}} = \frac{\Sigma_{\sim j,j}}{\Sigma_{j,j}}$. Thus, we can deduce that $J^{(j)} = \frac{\Sigma_{\sim j,j}}{\Sigma_{j,j}}$, with $\Sigma_{\sim j,j}$ being the $j^{th}$ column of $\Sigma$, and the dependent Jacobian becomes
$J^d = [J^{(1)}, \dots, J^{(d)}] = \left[ \frac{\Sigma_{\bullet,1}}{\Sigma_{1,1}}, \dots, \frac{\Sigma_{\bullet,d}}{\Sigma_{d,d}} \right] = \frac{\Sigma}{s^2} = \mathcal{R}, \text{ since } \Sigma_{j,j} = s_j^2 = s^2, \text{ and } \Sigma = s^2 \mathcal{R}. \text{ The tensor metric is given by } G = \mathcal{R}^T \mathcal{R}.$

### 5.3. Comparisons between Exact Gradient and Estimated Gradients

For running the above PDE-based model using the R package deSolve [45], we are given $D = 0.0011$ and $s = 1.96$. The exact and formal gradient associated with the mean values of the initial conditions is obtained by running the corresponding adjoint model. For estimating the gradient using the proposed estimators, we consider $L = 2, 3$ and $N = 50, 100, 150, 200$. We also use $h = 1/\sqrt{N}$ and $V_j \sim N(0, 1)$, $j = 1, \dots, d = 50$. The Sobol sequence is used for generating the random values of the $V_j$s, and the Gram–Schmidt algorithm is applied to obtain perfect orthogonal vectors for a given $N$.

Figure 1 shows the comparisons between the estimated and the exact values of the formal gradient $\nabla f$ (i.e., $\rho(X_{j_1}, Z_{j_2}) = 0$) for $L = 1, 2$. Likewise, Figures 2 and 3 depict the dependent gradient $\text{grad}(f) = (\mathcal{R}^T \mathcal{R})^{-1} \nabla f$ and its estimation. The estimates of both gradients are in line with the exact values using only $NL = 50$ (respectively, $NL = 100$) model evaluations when $L = 1$ and $N = 50$ (respectively, $L = 1$ and $N = 100$ or $L = 2$ and $N = 50$). Increasing the values of $L$ and $N$ gives the same quasiperfect results for both the formal and dependent gradients (see Figure 3).

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Exact gradient versus estimated gradients using $L = 1$ (◦) and $L = 2$ (+) of the QoI by considering the inputs as independent (formal gradients).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Exact gradient versus estimated gradients using $L = 1$ (∘) and $L = 2$ (+) of the QoI by considering the autocorrelations among the inputs (dependent gradients).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Exact gradient versus estimated gradients using $L = 2$ (∘) and $L = 3$ (+) of the QoI by considering the autocorrelations among the inputs (dependent gradients).</figcaption></figure>

**6. Conclusions**

In this paper, we have proposed new, simple, and generic approximations of the gradients of functions with nonindependent input variables by means of independent, central, and symmetric variables and a set of constraints. It comes out that the biases of our approximations for a wide class of functions, such as 2-smooth functions, do not suffer from the curse of dimensionality by properly choosing the set of independent, central, and symmetric variables. For functions including only independent input variables, a
theoretical comparison has shown that the upper bounds of the bias of the formal gradient derived in this paper outperformed the best known results.

For computing the dependent gradient of the function of interest, we have provided estimators of such a gradient by making use of evaluations of that function at LN randomized points. Such estimators reach the optimal (mean squared error) rates of convergence (i.e., $O(N^{-1}d^2)$) for a wide class of functions. Numerical comparisons using a test case and simulations based on a PDE model with given autocollaborations among the initial conditions have shown the efficiency of our approach, even when $L = 1, 2$ constraints were used. Our approach is therefore flexible, thanks to $L$ and the fact that the gradient can be computed for every value of the sample size $N$ in general.

While the proposed estimators reach the parametric rate of convergence, note that the second-order moments of such estimators depend on $d^2$. An attempt to reach a dimension-free rate of convergence requires working in $\mathbb{C}$ rather than $\mathbb{R}$ when $L = 2$. In the future, it is worth investigating the derivation of the optimal rates of convergence that are dimension-free or (at least) are linear with respect to $d$ by considering $L > 3$ constraints. Also, combining such a promising approach with a transformation of the original space might be helpful for reducing the number of model evaluations in higher dimensions.

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

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

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

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Proof of Theorem 1

$$ \text{As } \vec{t} = (t_1, \dots, t_d), \text{ let } \vec{k} = \left( \begin{array}{c} 0, \dots, 0, \\ \multicolumn{1}{c}{\underbrace{\phantom{0, \dots, 0}}_{\text{k th position}}} \\ 0, \dots, 0 \end{array} \right) \in \mathbb{R}^d \text{ and } \vec{q} = (q_1, \dots, q_d) \in \mathbb{N}^d. \text{ Multiplying the Taylor expansion of } f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \text{ about } \mathbf{x}, \text{ that is,} $$

$$ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{p=0}^{m} \sum_{\|\vec{t}\|_1=p} \frac{\mathcal{D}^{(\vec{t})} f(\mathbf{x})}{\vec{t}!} \beta_{\ell}^{p} (\mathbf{h} \mathbf{V})^{\vec{t}} + \mathcal{O}(\|\beta_{\ell} \mathbf{h} \mathbf{V}\|_2^{m+1}), $$

by $\frac{\mathbf{V}h^{-1}}{\sigma^2} \in \mathbb{R}^d$ and the constant $C_\ell$, as well as taking the sum over $\ell = 1, \dots, L$, we can see that the expectation $E := \sum_{\ell=1}^L C_\ell \mathbb{E}[f(\mathbf{x} + \beta_\ell \mathbf{h} \mathbf{V}) \frac{\mathbf{V}h^{-1}}{\sigma^2}]$ becomes

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

Firstly, for a given $k \in \{1, \dots, d\}$ and by independence, we can see that

$$ \mathbb{E}[(\mathbf{V})^{\vec{t}}(\mathbf{h})^{\vec{t}} V_k h_k^{-1}] = \mathbb{E}[(\mathbf{V})^{\vec{t}+\vec{k}}(\mathbf{h})^{\vec{t}-\vec{k}}] \neq 0 $$

if $t_k = 2q_k + 1$; thus, $t_j = 2q_j$ for any $j \in \{1, \dots, d\} \setminus \{k\}$, with $q_k, q_j \in \mathbb{N}$, which implies that $\vec{t} = \vec{k} + 2\vec{q}$. Thus, one obtains $\frac{\partial f}{\partial x_k}$ when $\|\vec{t}\|_1 = \|\vec{k} + 2\vec{q}\|_1 = 1$, and the fact that $\mathbb{E}[V_k^2] = \sigma^2$; $\mathbb{E}[V_j] = 0$; and $\sum_{\ell} C_\ell \beta_\ell = 1$. At this point, by taking $k = 1, \dots, d$ and setting $L = 1, \beta_\ell = 1$ and $C_\ell = 1$, this results in the approximation of $\nabla f(\mathbf{x})$ of order $O(\|h\|_2^2)$, because when $\|\vec{t}\|_1 = 2$, $\mathbb{E}[(\mathbf{V})^{\vec{t}}(\mathbf{h})^{\vec{t}} V_k h_k^{-1}] = 0$.
Secondly, for $L > 1$ the constraints $\sum_{\ell=1}^{L} C_{\ell} \beta_{\ell}^{r+1} = \delta_{0,r}$ $r = 0, 2, \dots, 2(L-1)$, we eliminate some higher-order terms so as to reach the order $O(||h||_2^{2L})$. Using other constraints, we complete the proof, thus bearing in mind Equation (2).

## Appendix B. Proof of Corollary 1

For $\vec{q} = (q_1, \dots, q_d) \in \mathbb{N}^d; k \in \{1, \dots, d\}$ and $\vec{k} = \left( \begin{array}{c} 0, \dots, 0, \\ \underbrace{1}_{k \text{ th position}}, \dots, 0 \end{array} \right) \in \mathbb{R}^d$, consider $s_k := \{\vec{q} + \vec{k} : ||\vec{q}||_1 = 1\}$. As $f \in \mathcal{H}_2$, we can write

$$f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) = \sum_{||\vec{i}||_1=0}^{1} \mathcal{D}^{(\vec{i})} f(\mathbf{x}) \beta_{\ell}^{||\vec{i}||_1} \frac{(\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!} + \sum_{||\vec{i}||_1=2} \mathcal{D}^{(\vec{i})} f(\mathbf{x}) \frac{\beta_{\ell}^2 (\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!} + R_k(\mathbf{h}, \beta_{\ell}, \mathbf{V}),$$

with the remainder term

$$R_k(\mathbf{h}, \beta_{\ell}, \mathbf{V}) = \sum_{\substack{||\vec{i}||_1=2 \\ \vec{i} \in s_k}} \mathcal{D}^{(\vec{i})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\beta_{\ell}^2 (\mathbf{h} \mathbf{V})^{\vec{i}}}{\vec{i}!} = \sum_{\substack{||\vec{q}||_1=1 \\ \vec{q} \in s_k}} \mathcal{D}^{(\vec{k}+\vec{q})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\beta_{\ell}^2 (\mathbf{h} \mathbf{V})^{\vec{q}+\vec{k}}}{(\vec{k}+\vec{q})!} \\ = h_k V_k \beta_{\ell}^2 \sum_{||\vec{q}||_1=1} \mathcal{D}^{(\vec{k}+\vec{q})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{(\mathbf{h} \mathbf{V})^{\vec{q}}}{(\vec{k}+\vec{q})!}.$$

Denote $R_k^0 := \sum_{||\vec{q}||_1=1} \mathcal{D}^{(\vec{k}+\vec{q})} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{(\mathbf{h} \mathbf{V})^{\vec{q}}}{(\vec{k}+\vec{q})!}$, and remark that $|R_k^0| \le M_2 ||\mathbf{h} \mathbf{V}||_1$. Using Theorem 1, we can see that the absolute value of the bias, that is,

$$B := \left\| \mathrm{grad}(f)(\mathbf{x}) - G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L} C_{\ell} \mathbb{E} \left[ f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_1 \text{ is given by}$$

$$B = \left\| G^{-1}(\mathbf{x}) \mathbb{E} \left[ \nabla f(\mathbf{x}) - \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_1 \\
= \left\| G^{-1}(\mathbf{x}) \sum_{\ell=1}^{L} C_{\ell} \frac{\beta_{\ell}^2}{\sigma^2} \mathbb{E} \left[ V_1^2 R_1^0, \dots, V_d^2 R_d^0 \right]^T \right\|_1 \\
\le \sum_{\ell=1}^{L} |C_{\ell}| \frac{\beta_{\ell}^2}{\sigma^2} \left\| G^{-1}(\mathbf{x}) \mathbb{E} \left[ V_1^2 R_1^0, \dots, V_d^2 R_d^0 \right]^T \right\|_1 \\
\le \sum_{\ell=1}^{L} |C_{\ell}| \beta_{\ell}^2 M_2 \left\| G^{-1}(\mathbf{x}) \mathbb{E} \left[ \frac{\mathbf{V}^2}{\sigma^2} ||\mathbf{h} \mathbf{V}||_1 \right] \right\|_1,$$

using the expansion of the product between matrices.

Using the same reasoning and taking the Euclidean norm, we obtain

$$B_2 = \left\| G^{-1}(\mathbf{x}) \mathbb{E} \left[ \nabla f(\mathbf{x}) - \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} \mathbf{h} \mathbf{V}) \frac{\mathbf{V} \mathbf{h}^{-1}}{\sigma^2} \right] \right\|_2 \\
\le \sum_{\ell=1}^{L} |C_{\ell}| \beta_{\ell}^2 M_2 \left\| G^{-1}(\mathbf{x}) \mathbb{E} \left[ \frac{\mathbf{V}^2}{\sigma^2} ||\mathbf{h} \mathbf{V}||_1 \right] \right\|_2.$$

The results hold using $\mathbf{R} := \mathbf{V}/\sigma$.
## Appendix C. Proof of Theorem 2

Firstly, remark that $MSE := \mathbb{E}[||\widehat{grad}(f)(\mathbf{x}) - grad(f)(\mathbf{x})||_2^2]$ is given by

$$MSE = \mathbb{E}[||\widehat{grad}(f)(\mathbf{x}) - \mathbb{E}[\widehat{grad}(f)(\mathbf{x})]||_2^2 + ||\mathbb{E}[\widehat{grad}(f)(\mathbf{x})] - grad(f)(\mathbf{x})||_2^2].$$

Since the bias $\mathbb{E}[||\mathbb{E}[\widehat{grad}(f)(\mathbf{x})] - grad(f)(\mathbf{x})||_2^2]$ has been derived in previous Corollaries, we are going to treat the second-order moment.

Secondly, since $f \in \mathcal{H}_2$ implies that $f \in \mathcal{H}_1$, we have $|f(\mathbf{x} + \beta_{\ell}h\mathbf{V}) - f(\mathbf{x})| \le M_1||\beta_{\ell}h\mathbf{V}||_2$. Also, since $\sum_{\ell=1}^L C_{\ell} = 0$, we then have

$$\sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) = \sum_{\ell=1}^{L} C_{\ell} [f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) - f(\mathbf{x})],$$

which leads to $|\sum_{\ell=1}^L C_{\ell}^{(|u|)} f(\mathbf{x} + \beta_{\ell} h \mathbf{V})| \le \sum_{\ell=1}^L |C_{\ell} \beta_{\ell}| M_1 ||h\mathbf{V}||_2$ and

$$Q(\mathbf{x}) := G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) = G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \sum_{\ell=1}^{L} C_{\ell} [f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) - f(\mathbf{x})]. \quad (A1)$$

Thirdly, using (3), we can see that $\mathbb{E}[Q(\mathbf{x})] = \mathbb{E}[\widehat{grad}(f)(\mathbf{x})]$. Bearing in mind the definition of the Euclidean norm and the variance, the centered second-order moment, that is, $V_{grad} := \mathbb{E}[||\widehat{grad}(f)(\mathbf{x}) - \mathbb{E}[\widehat{grad}(f)(\mathbf{x})]||_2^2]$ is given by

$$ \begin{aligned} V_{grad} &\le \frac{1}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) - \mathbb{E}[\widehat{grad}(f)(\mathbf{x})] \right\|_2^2 \right] \\ &\le \frac{1}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \sum_{\ell=1}^{L} C_{\ell} f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) \right\|_2^2 \right] \\ &\stackrel{(A1)}{=} \frac{1}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \sum_{\ell=1}^{L} C_{\ell} \{f(\mathbf{x} + \beta_{\ell} h \mathbf{V}) - f(\mathbf{x})\} \right\|_2^2 \right] \\ &\le \frac{1}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \right\|_2^2 \|h\mathbf{V}\|_2^2 \right] M_1^2 \left( \sum_{\ell=1}^{L} |C_{\ell} \beta_{\ell}| \right)^2 \\ &\le \frac{1}{N} \mathbb{E} \left[ \left\| G^{-1}(\mathbf{x}) \frac{\mathbf{V}\mathbf{h}^{-1}}{\sigma^2} \right\|_2^2 \|\mathbf{V}^2\|_2 \right] \|h^2\|_2 M_1^2 \left( \sum_{\ell=1}^{L} |C_{\ell} \beta_{\ell}| \right)^2 \end{aligned} $$

when bearing in mind the Hölder inequality. The results hold using $R := V/\sigma$ and the fact that when $V_k \sim U(-\xi, \xi)$, $||V^2||_2 \le \sqrt{d\xi^2}$, and $\sigma^2 = \xi^2/3$.

## References

1. Rosenblatt, M. Remarks on a Multivariate Transformation. Ann. Math. Statist. 1952, 23, 470-472. [CrossRef]

2. Nataf, A. Détermination des distributions dont les marges sont données. Comptes Rendus L'Académie Des Sci. 1962, 225, 42-43.

3. Joe, H. Dependence Modeling with Copulas; Chapman & Hall/CRC: London, UK, 2014.

4. McNeil, A.J.; Frey, R.; Embrechts, P. Quantitative Risk Management; Princeton University Press: Princeton, NJ, USA; Oxford, UK, 2015.

5. Navarro, J.; Ruiz, J.M.; Aguila, Y.D. Multivariate weighted distributions: A review and some extensions. Statistics 2006, 40, 51-64. [CrossRef]

6. Sklar, A. Fonctions de Rpartition à n Dimensions et Leurs Marges. Publ. l'Institut Stat. L'Université Paris 1959, 8, 229-231.

7. Durante, F.; Ignazzi, C.; Jaworski, P. On the class of truncation invariant bivariate copulas under constraints. J. Math. Anal. Appl. 2022, 509, 125898. [CrossRef]
8. Skorohod, A.V. On a representation of random variables. *Theory Probab. Appl.* **1976**, *21*, 645-648.

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

10. Lamboni, M. Efficient dependency models: Simulating dependent random variables. *Math. Comput. Simul.* **2022**, *200*, 199-217. [CrossRef]

11. Lamboni, M. On exact distribution for multivariate weighted distributions and classification. *Methodol. Comput. Appl. Probab.* **2023**, *25*, 41. [CrossRef]

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

13. Kucherenko, S.; Klymenko, O.; Shah, N. Sobol' indices for problems defined in non-rectangular domains. *Reliab. Eng. Syst. Saf.* **2017**, *167*, 218-231. [CrossRef]

14. Lamboni, M. On dependency models and dependent generalized sensitivity indices. *arXiv* **2021**, arXiv2104.12938.

15. Lamboni, M. Derivative Formulas and Gradient of Functions with Non-Independent Variables. *Axioms* **2023**, *12*, 845. [CrossRef]

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

17. Patelli, E.; Pradlwarter, H. Monte Carlo gradient estimation in high dimensions. *Int. J. Numer. Methods Eng.* **2010**, *81*, 172-188. [CrossRef]

18. Agarwal, A.; Dekel, O.; Xiao, L. Optimal Algorithms for Online Convex Optimization with Multi-Point Bandit Feedback. In *Proceedings of the The 23rd Conference on Learning Theory, COLT 2010*, Haifa, Israel, 27-29 June 2010; pp. 28-40.

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

20. Akhavan, A.; Pontil, M.; Tsybakov, A.B. *Exploiting Higher Order Smoothness in Derivative-Free Optimization and Continuous Bandits*; Red Hook, NY, USA, 2020; NIPS '20. Available online: https://arxiv.org/abs/2006.07862 (accessed on 18 June 2024).

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

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

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

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

25. Fruth, J.; Roustant, O.; Kuhnt, S. Total interaction index: A variance-based sensitivity index for second-order interaction screening. *J. Stat. Plan. Inference* **2014**, *147*, 212-223. [CrossRef]

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

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

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

29. Bobkov, S. Isoperimetric and Analytic Inequalities for Log-Concave Probability Measures. *Ann. Probab.* **1999**, *27*, 1903-1921. [CrossRef]

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

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

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

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

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

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

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

37. Guidotti, E. Calculus: High-Dimensional Numerical and Symbolic Calculus in R. *J. Stat. Softw.* **2022**, *104*, 1-37. [CrossRef]

38. Ancell, B.; Hakim, G.J. Comparing Adjoint- and Ensemble-Sensitivity Analysis with Applications to Observation Targeting. *Mon. Weather Rev.* **2007**, *135*, 4117-4134. [CrossRef]

39. Pradlwarter, H. Relative importance of uncertain structural parameters. Part I: algorithm. *Comput. Mech.* **2007**, *40*, 627-635. [CrossRef]

40. Polyak, B.; Tsybakov, A. Optimal accuracy orders of stochastic approximation algorithms. *Probl. Peredachi Inf.* **1990**, *26*, 45-53.
41. 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.

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

43. Lamboni, M. *Kernel-based Measures of Association Between Inputs and Outputs Using ANOVA*. *Sankhya A* **2024**. [CrossRef]

44. Gilbert, P.; Varadhan, R. *R-Package numDeriv: Accurate Numerical Derivatives*; CRAN Repository. 2019. Available online: http://optimizer.r-forge.r-project.org/ (accessed on 18 June 2024).

45. Soetaert, K.; Petzoldt, T.; Setzer, R.W. *Solving Differential Equations in R: Package deSolve*. *J. Stat. Softw.* **2010**, *33*, 1–25. [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.