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

HAL
open science

# Mollified Active Subspace

Romain Verdière, Clémentine Prieur, Olivier Zahm

► To cite this version:

Romain Verdière, Clémentine Prieur, Olivier Zahm. Mollified Active Subspace. 2025. hal-05194926v2

HAL Id: hal-05194926

https://inria.hal.science/hal-05194926v2

Preprint submitted on 1 Aug 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
Mollified Active Subspace*

Romain Verdière†, Clémentine Prieur †, and Olivier Zahm †

**Abstract.** High-dimensional computational models are often expensive to evaluate, making surrogate modeling a critical tool for efficient analysis. Furthermore, constructing accurate surrogates in the small-sample regime remains a major challenge due to the curse of dimensionality. A viable strategy is to reduce the input dimension by identifying informative low-dimensional feature. Gradient-based methods such as Active Subspaces (AS) achieve this by analyzing the spectrum of the second moment of the gradient of the model. However, the performance of AS deteriorates in the presence of high-frequency, low-amplitude oscillations, which makes gradient evaluation non-informative for dimension reduction and leads to suboptimal feature selection. In this work, we introduce the Mollified Active Subspace (MAS) method, which mitigates these issues by applying an anisotropic mollifying operator to the model. We derive a tractable upper bound on the approximation error involving the mollified model and show that this bound can improve upon the original AS bound, enabling more reliable feature selection. Our analysis reveals a fundamental trade-off between smoothing strength and approximation fidelity, and we provide practical guidelines for selecting mollification parameters. Numerical experiments on a range of models with oscillatory behavior demonstrate that MAS consistently outperforms AS in feature selection leading to better surrogate accuracy.

**Key words.** Gradient based dimension reduction, Active Subspace, Mollifiers, High dimensional function approximation, Surrogate modeling.

**MSC codes.** 41A30, 41A63, 65D15, 81Q35

# 1. Introduction.
Modern computational models in science and engineering often involve a large number of input parameters and are expensive to evaluate, both in terms of time and computational resources. To alleviate this burden, these models are frequently replaced with surrogate models. However, building accurate surrogates for high-dimensional models remains a challenging task, especially in the small-sample regime where only a limited number of model evaluations is available. Traditional approximation tools such as polynomials, wavelets, or neural networks often struggle in this setting due to the curse of dimensionality. A common strategy is to reduce the input dimension before to construct the surrogate. This involves identifying a *feature map* which maps the high-dimensional input parameters to a low-dimensional variable which explains, in some sense, the model output. In global sensitivity analysis [11], for instance, the feature map is taken as a coordinate selection operator. Various sensitivity measures have been proposed to identify the relevant coordinates, the most commonly used are Sobol' indices and the derivative-based global sensitivity measure (DGSM) [27, 20], to cite just a few. DGSMs are typically cheaper to compute than Sobol' indices and, using Poincaré inequalities, they provide bounds for Sobol' indices [26]. The *Active Subspace* method [30, 9, 8] permits to identify *linear feature maps* from gradients eval-

*Submitted to the editors 07/31/2025.

**Funding:** This work was funded by the French National Research Agency (ANR) under the MODENA project (ANR-21-CE46-0006-01) and the BALTEEC Chair within the MIAI Cluster project (ANR-23-IACL-0006).

†Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France (romain.verdiere@inria.fr, clementine.prieur@univ-grenoble-alpes.fr, olivier.zahm@inria.fr).
uations of the model, assuming they are available. AS has been successfully applied in various domains to build efficient surrogate models [19, 15, 10]. Nonlinear feature have also been investigated, including approaches based on polynomial feature maps [5] and diffeomorphism-based transformations [28, 32], which allow for more flexible representations of the input-output relationship. All these methods require gradients of the model to estimate the Poincaré bounds, which can be computed via automatic differentiation [13, 2] or the adjoint-state method [24], typically at a cost comparable to evaluating the model itself.

In this paper, we focus on gradient-based methods for linear feature map detection. The standard AS method defines the feature maps using the dominant eigenvectors of the AS matrix

$$ (1.1) \qquad H(u) = \mathbb{E} \left[ \nabla u(\mathbf{X}) \nabla u(\mathbf{X})^{\top} \right] \in \mathbb{R}^{d \times d}, $$

where $u: \mathbb{R}^d \to \mathbb{R}$ denotes the model of interest and $\mathbf{X}$ the input random vector taking values in $\mathbb{R}^d$, $d \gg 1$. For simplicity, we assume that $\mathbf{X}$ follows a standard multivariate Gaussian distribution $N(0, I_d)$. The Gaussian setting enables meaningful theoretical analysis, particularly for single- and multiple-index models [21, 3]. It is a fairly standard assumption in the theoretical machine learning literature [25]. Working with the standard normal distribution is not limiting, since any normal distribution can be standardised, as explained in [7]. It is shown in [30, 31, 23] that the AS method consists in minimizing an error bound obtained via Poincaré inequalities. As mentioned in [30, Section 5.1.3], high-frequency, low-amplitude oscillations in the model may be overemphasized in the gradient-based bound, leading to irrelevant feature maps. We propose to utilize a mollifying operator $P_{t,M}(\cdot)$ in order to smooth out the high-frequencies of the model before applying the AS method. This mollifying operator is defined as

$$ P_{t,M}(u)(x) = \mathbb{E} \left[ u \left( e^{-tM}x + \sqrt{I_d - e^{-2tM}} \mathbf{Z} \right) \right], $$

where $\mathbf{Z} \sim N(0, I_d)$ is independent on $\mathbf{X}$ and where the scalar $t \ge 0$ and the positive matrix $M \in \mathbb{R}^{d \times d}$ are some parameters to be tuned. We demonstrate both theoretically and numerically that using the matrix $H(P_{t,M}(u))$ in place of $H(u)$ can improve the Active Subspace bound, and this improvement leads to better and more robust feature map constructions. This new bound decomposes as the sum of two terms. The first is the mean squared error between $u$ and $P_{t,M}(u)$, hence it vanished when $t \to 0$. The second corresponds to the AS bound applied to the mollified model $P_{t,M}(u)$, hence it is small when $t \to \infty$ (indeed, $H(P_{t,M}(u)) \to 0$ as $P_{t,M}(u)(x)$ tends to a constant function when $t \to \infty$). This highlights a fundamental trade-off between mollification and approximation accuracy: stronger smoothing $t \gg 1$ reduces the gradient variance but may also degrade the approximation. We explore various choices for the matrix $M$ and provide a practical method for selecting the mollification parameter $t$. Based on this analysis, we propose the *Mollified Active Subspace method*, which improves upon classical Active Subspace by reducing selection errors and yielding sharper upper bounds for models with oscillatory behavior.

The paper is organized as follows. In Section 2, we review the Active Subspace (AS) method and illustrate its limitations through a simple example. Subsection 3.2 analyzes key properties of the mollification operator $P_{t,M}(u)$, and Subsection 3.3 presents our main
theoretical result. In Subsection 4.1, we discuss practical strategies for selecting the mollification parameters $t$ and $M$. Subsection 4.2 describes how to estimate the MAS upper bound efficiently by reusing gradient snapshots from the AS computation. The complete Mollified Active Subspace procedure is detailed in Subsection 4.3. Finally, Subsection 4.4 illustrates the differences between MAS and AS on an analytical example, and Section 5 presents numerical comparisons on two benchmark problems.

**2. Active Subspace method and its limitations.** Linear dimension reduction consists in approximating a given high dimensional function $u: \mathbb{R}^d \to \mathbb{R}$ by a composed function $f \circ g$, where $f: \mathbb{R}^m \to \mathbb{R}$ is the profile function and where $g: \mathbb{R}^d \to \mathbb{R}^m$ is a linear feature map of the form $g(x) = U_m^\top x$ for some matrix $U_m \in \mathbb{R}^{d \times m}$ with orthogonal columns. The Active Subspace (AS) method consists in constructing the matrix $U_m$ whose columns are the first eigenvectors of the AS matrix $H(u)$ as in (1.1). In [30], it is shown that the AS method consists in minimizing an upper bound of the following reconstruction error

$$ \min_{f: \mathbb{R}^m \to \mathbb{R} \text{ measurable}} \mathbb{E}[(u(\mathbf{X}) - f(U_m^\top \mathbf{X}))^2] = \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}])^2], $$

for any $U_m$. Here the function $f$ which realizes the minimum is the conditional expectation $f: z \mapsto \mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}] = z$. Specifically when $\mathbf{X} \sim \mathcal{N}(0, I_d)$ is a standard Gaussian random vector, we have the following result.

**Proposition 2.1.** (Active Subspace [30]) Let $u: \mathbb{R}^d \to \mathbb{R}$ be a differentiable function such that $\mathbb{E}[u(\mathbf{X})^2] < +\infty$ and $\mathbb{E}[\|\nabla u(\mathbf{X})\|^2] < +\infty$, where $\mathbf{X} \sim \mathcal{N}(0, I_d)$. Then, for any matrix $U_m \in \mathbb{R}^{d \times m}$ with orthonormal columns, we have

$$ (2.1) \qquad \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}])^2] \le \mathbb{E}[\|(I_d - U_m U_m^\top) \nabla u(\mathbf{X})\|_2^2]. $$

Moreover, for $\{(\lambda_i, v_i), 1 \le i \le n\}$ denoting the set of eigenpairs of $H(u) = \mathbb{E}[\nabla u(\mathbf{X}) \nabla u(\mathbf{X})^\top]$ with $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_d$, we have

$$ (2.2) \qquad \min_{\substack{U_m \in \mathbb{R}^{d \times m} \\ U_m^\top U_m = I_m}} \mathbb{E}[\|(I_d - U_m U_m^\top) \nabla u(\mathbf{X})\|_2^2] = \sum_{i=m+1}^d \lambda_i, $$

and a solution to the above minimization problem is $U_m^* = [v_1, \dots, v_m]$.

The Active Subspace method consists in identifying $U_m$ by minimizing the right-hand side of (2.1). In particular, Equation (2.2) reveals that models with a rapidly decaying spectrum of $H(u)$ possess a low intrinsic dimension. This minimize-then-optimize strategy is effective when the bound in (2.1) is tight. While this holds true for smoothly varying models (for instance, (2.1) becomes an equality for affine functions $u$), Inequality (2.1) can be loose for models with strong oscillations, where gradients do not reliably capture the model global variations. Applying the AS method to these models can lead to selection errors in the feature map, as illustrated on the following toy model.

Let

$$ \begin{array}{l} u: \mathbb{R}^2 \quad \to \mathbb{R} \\ (x_1, x_2) \quad \mapsto \sin(x_1) + \frac{1}{6}\sin(10x_2). \end{array} $$
To simplify our analysis we suppose that $m = 1$ and that $U_m$ is either equal to $U_{\{1\}} = (\frac{1}{0})$ or $U_{\{2\}} = (\frac{0}{1})$. This corresponds to a coordinate selection strategy where $U_{\{1\}}$ selects $x_1$ and $U_{\{2\}}$ selects $x_2$. In this case the reconstruction errors and the upper-bounds read:

$$ \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_{\{1\}}^{\top}\mathbf{X}])^2] = \frac{1}{72}(1 - e^{-200}) \approx 0.14 $$

$$ \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_{\{2\}}^{\top}\mathbf{X}])^2] = \frac{1}{2}(1 - e^{-2}) \approx 0.43 $$

$$ \mathbb{E}[\|(I_d - U_{\{1\}}U_{\{1\}}^{\top})\nabla u(\mathbf{X})\|_2^2] = \frac{100}{72}(1 + e^{-200}) \approx 1.39 $$

$$ \mathbb{E}[\|(I_d - U_{\{1\}}U_{\{2\}}^{\top})\nabla u(\mathbf{X})\|_2^2] = \frac{1}{2}(1 + e^{-2}) \approx 0.57. $$

See Appendix G for calculus details. Minimizing the AS bound results in selecting $x_2$ as the most active coordinate (since $0.57 < 1.39$), and gives a rather good estimation of the resulting error (since $0.43 < 0.57$). However, selecting $x_2$ is clearly suboptimal compared to selecting $x_1$ (since $0.14 < 0.43$).

To avoid such selection errors, we propose to identify $U_m$ using the following Active Sub-space matrix:

$$ H(P_{t,M}(u)) = \mathbb{E}[\nabla P_t(u)(\mathbf{X})\nabla P_t(u)(\mathbf{X})^{\top}], $$

where given $t > 0$ the operator $P_t$ is defined by:

$$ P_t(u)(x) = \mathbb{E}[u(e^{-t}x + \sqrt{1 - e^{-2t}\mathbf{Z}})], $$

for $\mathbf{Z} \sim \mathcal{N}(0, I_d)$ independent of $\mathbf{X}$. The operator $P_t(\cdot)$ regularizes $u$ based on the parameter $t$, helping to mitigate model oscillations and prevent selection errors like those described above. This regularizing effect is illustrated on Figure 2, where we plot $u$ and $P_t(u)$ for $x_1 = x_2$ across different values of $t$. Figure 1 plots the Poincaré upperbound of $P_t(u)$ for $U_m = U_{\{1\}}$, $U_m = U_{\{2\}}$ as a function of $t$. Using $P_t(u)$ instead of $u$ allows to correct the selection error as soon as $t \ge 5 \times 10^{-3}$, however there is no guarantee that $\mathbb{E}[\|(I_d - U_mU_m^{\top})\nabla P_t(u)(\mathbf{X})\|_2^2]$ bounds the reconstruction error $\mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_m^{\top}\mathbf{X}])^2]$ (the blue curve goes below the green dashed line in Figure 1).

In the following sections, we derive in Proposition 3.6 an upper bound on the reconstruction error $\mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_m^{\top}\mathbf{X}])^2]$ that involves $P_t(u)$ instead of $u$. Minimizing this bound helps to prevent selection errors while also reducing the original Active Subspace upper bound for certain oscillatory models. Additionally, we provide practical techniques for selecting the parameter $t$.

**Remark 2.2.** In the statistical learning literature, the multi-index model [14, 16, 4] poses the relation $\mathbf{Y} = f(U_m^{*\top}\mathbf{X}) + \zeta$ for some unknown index space $U_m^* \in \mathbb{R}^{d \times m}$ and link function $f: \mathbb{R}^m \to \mathbb{R}$ with $\zeta$ a noise term such that $\mathbb{E}[\zeta|\mathbf{X}] = 0$. [14, 16, 4] study the statistical convergence of $U_m$ to $U_m^*$ with respect to the number of samples, or to the number of iterations of the stochastic gradient descent. Recent work in this community [29] uses Smoothed Gradient Outer Products (SGOP) to estimate the index space $U_m^*$ with a favorable convergence rate. Gradient Outer Products correspond to the Active Subspace matrix $H(u)$ while SGOP corresponds to the Active Subspace matrix of a mollified version of the model ($H(P_{t,M}(u))$ in our work). The framework we consider in the present work is different in the sense that we do not assume that $\mathbb{E}[\mathbf{Y}|\mathbf{X}]$ is of the form $f(U_m^{\top}\mathbf{X})$. We rather seek for a certified approximation of $u(\mathbf{X})$ in the form $f(U_m^{\top}\mathbf{X})$.
<figure><img src="image_3.png" /><figcaption>Figure 1: $P_t(u)$ AS bound as a function of $t$. $P_t(\cdot)$ corrects the selection error for $t \ge 5 \times 10^{-3}$, however $E[\|(I_d - U_m U_m^\top) \nabla P_t(u)(\mathbf{X})\|_2^2]$ no longer bounds the reconstruction error for all $t$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Plots of $u$ and $P_t(u)$ for $x_1 = x_2$ and for different values of $t$.</figcaption></figure>

**3. Mollified Active Subspace (MAS).** This section presents the main result of the paper, stated in Proposition 3.6. We begin by establishing a set of preliminary results on the mollifying operator. These results are then used to derive a new bound on the reconstruction error $E[(u(\mathbf{X}) - E[u(\mathbf{X})|U_m^\top \mathbf{X}])^2]$.

**3.1. Notations and definitions.** Let us start with the definition of the anisotropic mollification operator acting on the Hilbert space of function $L_\mu^2(\mathbb{R}^d) = \{u : \mathbb{R}^d \to \mathbb{R} : E[u(\mathbf{X})^2] < \infty\}$, where $\mathbf{X} \sim \mu = \mathcal{N}(0, I_d)$, endowed with the scalar product $\langle u, v \rangle = E[u(\mathbf{X})v(\mathbf{X})]$.

**Definition 3.1 (Anisotropic mollification operator).** For any $t \ge 0$ and any symmetric semi-definite positive matrix $M \in \mathbb{R}^{d \times d}$, we define the operator $P_{t,M} : L_\mu^2(\mathbb{R}^d) \to L_\mu^2(\mathbb{R}^d)$ by

$$ (3.1) \qquad P_{t,M}(u)(x) = E \left[ u \left( e^{-tM} x + \sqrt{I_d - e^{-2tM}} Z \right) \right], $$
for any $u \in L_{\mu}^2(\mathbb{R}^d)$, where $Z \sim \mathcal{N}(0, I_d)$. Here $e^A$ is the matrix exponential such that for $A \in \mathbb{R}^{d \times d}$, $e^A = \sum_{k=0}^{\infty} \frac{t^k}{k!} A^k$ and, for $B \in \mathbb{R}^{d \times d}$ a positive semi-definite matrix, $\sqrt{B}$ is the unique symmetric matrix which verifies $\sqrt{B}^2 = B$.

When $M = I_d$, the isotropic mollification operators $\{P_{t,I_d}\}_{t \ge 0}$ form the Ornstein-Uhlenbeck semigroup, see [1]. When $M \ne \lambda I_d$, the operator $P_{t,M}$ performs an anisotropic mollification. In particular with $M = \text{diag}(\lambda_1, \dots, \lambda_d)$, each canonical direction of $\mathbb{R}^d$ is molidified according to $\lambda_i t$. In any case, $t = 0$ leads to $P_{t,M}(u) = u$. We now introduce notation for the reconstruction error and its bound.

**Definition 3.2.** For any function $u \in L_{\mu}^2(\mathbb{R}^d)$ with $\nabla u \in (L_{\mu}^2(\mathbb{R}^d))^d$ and any matrix $U_m \in \mathbb{R}^{d \times m}$ with orthonormal columns, we let

$$ (3.2) \qquad \mathcal{E}(U_m, u) := \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}])^2], $$

$$ (3.3) \qquad \mathcal{J}(U_m, u) := \mathbb{E}[\| (I_d - U_m U_m^\top) \nabla u(\mathbf{X}) \|_2^2]. $$

As stated in Proposition 2.1, the reconstruction error $\mathcal{E}(U_m, u)$ of $u(\mathbf{X})$ knowing $U_m^\top \mathbf{X}$ is bounded by the Poincaré bound $\mathcal{J}(U_m, u)$. In this section, we aim to derive a novel bound for $\mathcal{E}(U_m, u)$ in terms of $\mathcal{J}(U_m, P_{t,M}(u))$.

**3.2. Expansions on Hermite Basis.** The Hermite polynomials provide a convenient and powerful tool for our analysis. We briefly recall below their definition and key properties, see [6] for a comprehensive treatment.

The univariate probabilist's Hermite polynomials $h_n(\cdot)$ of order $n$, defined by

$$ (3.4) \qquad h_n(\xi) = (-1)^n e^{\frac{\xi^2}{2}} \frac{d^n}{d\xi^n} e^{-\frac{\xi^2}{2}}, $$

satisfy the following properties for $Z \sim \mathcal{N}(0, 1)$:

(Pa) $\mathbb{E}[h_n(Z)h_m(Z)] = n!\delta_{n,m}$,

(Pb) $h'_n(\xi) = nh_{n-1}(\xi)$,

(Pc) $P_t h_n(\xi) = e^{-nt}h_n(\xi)$, where $P_t: f \mapsto P_t f(\cdot) = \mathbb{E}[f(e^{-t}\cdot) + \sqrt{1-e^{-2t}}Z]$ is the Ornstein-Uhlenbeck semigroup acting on $L_{\mu}^2(\mathbb{R})$,

(Pd) $(h_n)_{n \in \mathbb{N}}$ is an orthogonal basis of $L_{\mu}^2(\mathbb{R}^d)$.

We consider now the multivariate Hermite polynomials defined as the tensorization of normalized univariate Hermite polynomials as follows. Let $\mu = \mathcal{N}(O, I_d)$, for any $\alpha \in \mathbb{N}^d$, we define $e_{\alpha} \in L_{\mu}^2(\mathbb{R}^d)$ as

$$ (3.5) \qquad e_{\alpha}(x_1, \dots, x_d) = \prod_{i=1}^{d} \frac{h_{\alpha_i}(x_i)}{\sqrt{\alpha_i!}}, $$

for any $x \in \mathbb{R}^d$. Then $(e_{\alpha})_{\alpha \in \mathbb{N}^d}$ is an orthonormal basis of $L_{\mu}^2(\mathbb{R}^d)$ and we can write:

$$ (3.6) \qquad u(x) = \sum_{\alpha \in \mathbb{N}^d} u_{\alpha} e_{\alpha}(x), \quad \text{with } u_{\alpha} := \mathbb{E}[u(\mathbf{X})e_{\alpha}(\mathbf{X})] \text{ and } \mathbf{X} = \mathcal{N}(O, I_d). $$

Similarly, we can expand the conditional expectation $\mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}]$ on the Hermite polynomial basis. The proof of the next lemma is given in Appendix A.
<figure><img src="image_5.png" /><figcaption>Figure 3: L²-normalized Hermite polynomials $\psi_n = \frac{h_n}{\sqrt{n!}}$ of degrees 1, 5 and 20. Polynomials of higher degree have stronger oscillatory behaviors.</figcaption></figure>

**Lemma 3.3.** Let $u \in L_\mu^2(\mathbb{R}^d)$. For any matrix $U_m \in \mathbb{R}^{d \times m}$ with orthonormal columns, let $U_\perp \in \mathbb{R}^{d \times (d-m)}$ be an orthonormal complement to $U_m$ so that $U = [U_m, U_\perp] \in \mathbb{R}^{d \times d}$ is a unitary matrix. Then we have

$$ (3.7) \qquad \mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X}] = \sum_{\substack{\alpha \in \mathbb{N}^d \text{ s.t.} \\ \alpha_{m+1}, \dots, \alpha_d = 0}} u_\alpha^U e_\alpha(U^\top \mathbf{X}), $$

where $u_\alpha^U := \mathbb{E}[u(\mathbf{X})e_\alpha(U^\top \mathbf{X})]$.

The next lemma, which proof is given in Appendix B, expresses $\mathcal{E}(U_m, u)$ and $\mathcal{J}(U_m, u)$ in terms of the coefficients of $u$ on the Hermite Basis. The proof mainly relies on (Pa) and (Pb).

**Lemma 3.4.** Let $u \in L_\mu^2(\mathbb{R}^d)$ and $U_m \in \mathbb{R}^{d \times m}$ be a matrix with orthonormal columns. Let $U_\perp \in \mathbb{R}^{d \times (d-m)}$ be any orthonormal complement to $U_m$ so that $U = [U_m, U_\perp] \in \mathbb{R}^{d \times d}$ is a unitary matrix. Then we have

$$ (3.8) \qquad \mathcal{E}(U_m, u) = \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \mathbf{1}_{|\alpha_\perp|>0}, $$

$$ (3.9) \qquad \mathcal{J}(U_m, u) = \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 |\alpha_\perp|, $$

where $u_\alpha^U = \mathbb{E}[u(\mathbf{X})e_\alpha(U^\top \mathbf{X})]$, $|\alpha_\perp| = \alpha_{m+1} + \dots + \alpha_d$ and where $\mathbf{1}(\cdot)$ is the indicator function such that $\mathbf{1}_{|\alpha_\perp|>0} = 0$ if $|\alpha_\perp| = 0$ and $\mathbf{1}_{|\alpha_\perp|>0} = 1$ otherwise.

Let us note that Lemma 3.4 yields a simple proof for Proposition 2.1: since $|\alpha_\perp| \in \mathbb{N}$ we have $\mathbf{1}_{|\alpha_\perp|>0} \le |\alpha_\perp|$ so that $\mathcal{E}(U_m, u) \le \mathcal{J}(U_m, u)$. Also, Equations (3.8) and (3.9) allow to better appreciate the gap between the reconstruction error $\mathcal{E}(\Pi_m, u)$ and its bound $\mathcal{J}(\Pi_m, u)$:
replacing $\mathbf{1}_{|\alpha_{\perp}|>0}$ with $|\alpha_{\perp}|$ yields a loose bound when $(u_{\alpha}^{U})^2 = \mathcal{O}(1)$ for large $|\alpha_{\perp}|$. Oscillatory functions tend to have larger high order coefficients (Hermite polynomials of higher degree have stronger oscillatory behaviors, see Figure 3) and can fool the AS method as illustrated in the previous section.

The next lemma, whose proof is given in Appendix C exploits (Pc) to explicit the decomposition of $P_{t,M}(u)$ on Hermite polynomials basis.

**Lemma 3.5.** Let $u \in L_{\mu}^{2}(\mathbb{R}^{d})$, $t \ge 0$ and $M \in \mathbb{R}^{d \times d}$ be a symmetric semi-definite positive matrix. Consider the eigenvalue decomposition $M = Q\text{diag}(\lambda)Q^{\top}$, where $Q \in \mathbb{R}^{d \times d}$ is a unitary matrix and $\lambda = (\lambda_{1}, \dots, \lambda_{d})$ the vector containing its eigenvalues $\lambda_{i} \ge 0$. Then, for any $u \in L_{\mu}^{2}(\mathbb{R}^{d})$, we have

$$ (3.10) \qquad P_{t,M}(u)(x) = \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^{Q} e^{-\alpha^{\top} \lambda t} e_{\alpha}(Q^{\top} x), $$

where $u_{\alpha}^{Q} = \mathbb{E}[u(\mathbf{X})e_{\alpha}(Q^{\top}\mathbf{X})]$.

Comparing (3.6) with (3.10), we see the regularizing effect of the operator $P_{t,M}$ which turns oscillatory models into more regular ones, as illustrated on Figure 2. Indeed, the term $e^{-\alpha^{\top}\lambda t}$ in (3.10) decreases exponentially fast with $\alpha^{\top}\lambda$, which diminishes drastically the effect of high-order terms.

**3.3. Main result.** Here we leverage the previous lemmas to show that the reconstruction error $\mathcal{E}(\Pi_{m}, u)$ can be upper-bounded by a sum of two terms: one which is the Poincaré upper-bound of the mollified model $\mathcal{J}(U_{m}, P_{t,M}(u))$ and one that accounts for the approximation error made by replacing $u$ with $P_{t,M}(u)$.

**Proposition 3.6.** Let $u \in L_{\mu}^{2}(\mathbb{R}^{d})$ such that $\nabla u \in (L_{\mu}^{2}(\mathbb{R}^{d}))^{d}$ and let $U_{m} \in \mathbb{R}^{d \times m}$ be a matrix with orthonormal columns. Furthermore, let $t \ge 0$ and $M \in \mathbb{R}^{d \times d}$ be a nonzero symmetric semi-definite positive matrix. Assume that the columns of $U_{m}$ are eigenvectors of the matrix $M$. Then

$$ (3.11) \qquad \mathcal{E}(U_m, u) \le \mathcal{J}_{t,M}(U_m, u) := \mathcal{J}(U_m, P_{t,M}(u)) + \frac{1 - e^{-2rt}}{r} \mathbb{E}[\|\nabla u(\mathbf{X})\|_{M}^{2}], $$

where $r > 0$ is the smallest nonzero eigenvalue of $M$, $\mathbf{X} \sim \mu = \mathcal{N}(0, I_{d})$ and $\|\nabla u(\mathbf{X})\|_{M}^{2} = \nabla u(\mathbf{X})^{\top} M \nabla u(\mathbf{X})$.

*Proof.* Because the columns of $U_{m}$ are eigenvectors of the matrix $M$, there exists a unitary matrix of the form $U = [U_{m}, U_{\perp}] \in \mathbb{R}^{d \times d}$ such that $M = U\text{diag}(\lambda)U^{\top}$, where $\lambda = (\lambda_{1}, \dots, \lambda_{d})$ with $\lambda_{i} \ge 0$. Then, Lemma 3.5 permits to write $P_{t,M}(u)(x) = \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^{U} e^{-\alpha^{\top} \lambda t} e_{\alpha}(U^{\top} x)$. Applying Lemma 3.4 to $P_{t,M}(u)$ thus yields $\mathcal{E}(U_{m}, P_{t,m}u) = \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^{U})^2 \mathbf{1}_{|\alpha_{\perp}|>0} e^{-2\alpha^{\top}\lambda t}$ and
then

$$
\begin{align*}
\mathcal{E}(U_m, u) &\stackrel{(3.8)}{=} \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \mathbf{1}_{|\alpha_\perp|>0} \\
&= \mathcal{E}(U_m, P_{t,M}(u)) + \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \mathbf{1}_{|\alpha_\perp|>0} (1 - e^{-2\alpha^\top \lambda t}) \\
&\le \mathcal{J}(U_m, P_{t,M}(u)) + \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \mathbf{1}_{|\alpha_\perp|>0} (\alpha^\top \lambda) \frac{1 - e^{-2\alpha^\top \lambda t}}{\alpha^\top \lambda},
\end{align*}
$$

where we used the Active Subspace bound as given by Proposition 2.1. Next, because $x \mapsto \frac{1-e^{-2xt}}{x}$ is decreasing for all $t \ge 0$, we deduce that $\frac{1-e^{-2\alpha^\top \lambda t}}{\alpha^\top \lambda} \le \frac{1-e^{-2rt}}{r}$ for all $\alpha \in \mathbb{N}^d \setminus \{0\}$. Then we obtain

$$
\begin{align*}
\mathcal{E}(U_m, u) &\le \mathcal{J}(U_m, P_{t,M}(u)) + \frac{1-e^{-2rt}}{r} \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \mathbf{1}_{|\alpha_\perp|>0} (\alpha^\top \lambda) \\
&\le \mathcal{J}(U_m, P_{t,M}(u)) + \frac{1-e^{-2rt}}{r} \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 (\alpha^\top \lambda).
\end{align*}
$$

It remains to show $\sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 (\alpha^\top \lambda) = \mathbb{E}[\|\nabla u(\mathbf{X})\|_M^2]$. As in the proof of Lemma 3.4, let $v : y \mapsto u(Uy)$ and $\mathbf{Y} = U\mathbf{X}$. As $M = U\mathrm{diag}(\lambda)U^\top$, we have

$$
\begin{align}
\mathbb{E}[\|\nabla u(\mathbf{X})\|_M^2] &= \mathbb{E}[\|\nabla v(\mathbf{Y})\|_{\mathrm{diag}(\lambda)}^2] \nonumber \\
&= \sum_{i=1}^d \lambda_i \mathbb{E}[(\partial_i \nabla v(\mathbf{Y}))^2] \nonumber \\
&\stackrel{\text{(B.1)}}{=} \sum_{i=1}^d \lambda_i \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \alpha_i \nonumber \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 (\alpha^\top \lambda). \tag{3.12}
\end{align}
$$

This concludes the proof.

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

Before exploiting the bound stated in Proposition 3.6 in the next section, let us note that the new upper bound $\mathcal{J}_{t,M}(U_m, u)$ for $\mathcal{E}(U_m, u)$ converges to the standard Active Subspace upper-bound $\mathcal{J}(U_m, u)$ when $t \to 0$. An important question is whether this new bound provides an improvement over the standard AS bound, and under what conditions. While the answer remains unclear in general, the following proposition, whose proof is given in Appendix D, provides sufficient (but not necessary) conditions under which $\mathcal{J}_{t,M}(U_m, u) \le \mathcal{J}(U_m, u)$ is satisfied.

**Proposition 3.7.** Under the assumptions of Proposition 3.6 and for $t$ small enough, we have

$$
(3.13) \qquad \mathcal{J}_{t,M}(U_m, u) = \mathcal{J}(U_m, u) - 2t \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \alpha^\top \lambda (|\alpha_\perp| - 1) + \mathcal{O}(t^2),
$$
where $\lambda$ is the vector containing the eigenvalues of $M$ and where $U$ is the unitary matrix such that $M = U \text{diag}(\lambda)U^\top$. If we further assume $M = I_d - U_m U_m^\top$, with $U_m$ a matrix formed with $m$ first columns of $U$ and $m < d$ then

$$ (3.14) \qquad \mathcal{J}_{t, I_d - U_m U_m^\top} (U_m, u) \le \mathcal{J}(U_m, u), $$

holds for any $t \ge 0$.

**Remark 3.8.** Choosing $M$ depending on $U_m$ is difficult to use in practice, the main issue is that minimizing the objective $U_m \mapsto \mathcal{J}_{t,M}(U_m)(U_m, u)$ can no longer be done in closed form. This contrasts with the case where $M$ is independent of $U_m$, in which case minimizing $U_m \mapsto \mathcal{J}_{t,M}(U_m, u) = \text{trace}((I_d - U_m U_m^\top)H(P_{t,M}u))$ reduces to computing the leading eigenvectors of $H(P_{t,M}u)$.

Let us note that Proposition 3.6 requires a rather strong condition between $U_m$ and $M$, specifically that $U_m$ are eigenvectors of the matrix $M$. This is equivalent to say that $U_m^\top M U_m$ is a diagonal matrix. This condition is trivially satisfied when $M = I_d$, which corresponds to isotropic mollification. Proposition E.1, given in Appendix E, is a variant of Proposition 3.6 that does not require this assumption between $M$ and $U_m$, at the cost of yielding a slightly looser bound.

**4. MAS practical implementation and procedure.** This section describes how to apply the Mollified Active Subspace method. We first explain how to select the mollifying parameters and estimate the mollified active subspace matrix. We then provide a practical implementation procedure for the method.

**4.1. Choosing $M$ and $t$.** We propose the following choices for $M$.

* $M = I_d$ is the simplest way to choose $M$, all the input variables are mollified equally. Also, the condition of Proposition 3.6 that the columns of $U_m$ are eigenvector of $M$ is satisfied for all $U_m$.

* $M = U_{m_0}^{\text{AS}}(U_{m_0}^{\text{AS}})^\top$, where $U_{m_0}^{\text{AS}} \in \mathbb{R}^{d \times m_0}$ contains the first $m_0$ eigenvectors of $H(u)$. The number of components $m_0$ can be chosen as the smallest number such that $\mathcal{J}(U_{m_0}^{\text{AS}}, u) \le \varepsilon$ for some prescribed $\varepsilon > 0$. The intuition behind this choice of $M$ is to mollify only along the directions pre-selected by the AS method. This would help identify whether unnecessary components were selected, that is, whether the AS method was misled by some oscillatory behavior in $u$. For Proposition 3.6 to hold, the matrix $U_m$ needs to contain eigenvectors of $U_{m_0}^{\text{AS}}(U_{m_0}^{\text{AS}})^\top$.

* $M = H(u)$ yields an anisotropic mollification operator that regularizes more the “active directions” as given by the Active Subspace matrix $H(u)$, the strength of the mollification being proportional to the corresponding eigenvalue. This strategy is a variant of the previous choice which avoids the prior selection of active components. For Proposition 3.6 to hold, the matrix $U_m$ needs to contain eigenvectors of $H(u)$.

Next, we discuss the choice of $t$. For a given matrix $M \in \mathbb{R}^{d \times d}$ and for a fixed tolerance $\varepsilon > 0$, we aim at choosing $t$ such that $\mathcal{J}_{t,M}(U_m, u) \le \varepsilon$. However $\mathcal{J}_{t,M}(U_m, u)$ cannot be computed without knowing $t$, therefore we propose to set $t$ such that $\frac{1-e^{-2rt}}{r} \mathbb{E}[\|\nabla u(\mathbf{X})\|_M^2] = \alpha\varepsilon$ where $\alpha \in (0, 1)$. This way, if there exists $U_m$ with $m \le m_0$ such that $\mathcal{J}(U_m, P_{t,M}(u)) \le (1 - \alpha)\varepsilon$,
then $\mathcal{J}_{t,M}(U_m, u) \le \epsilon$. Using this rule we obtain the following expression:

$$ (4.1) \qquad t_{\alpha} = -\frac{1}{2r} \ln \left( 1 - r \frac{\alpha \epsilon}{\mathbb{E}[\|\nabla u(\mathbf{X})\|_{M}^{2}]} \right). $$

In practice, we use $\alpha \in \{\frac{1}{8}, \frac{1}{4}, \frac{3}{8}, \frac{1}{2}\}$, see Subsection 4.4 and section 5.

**4.2. Estimation of $\mathcal{J}_{t,M}(U_m, u)$**. Let us rewrite $\mathcal{J}_{t,M}(U_m, u)$ defined in (3.11) as follows

$$ \mathcal{J}_{t,M}(U_m, u) = \text{trace} \left( (I_d - U_m U_m^\top) H(P_{t,M}(u)) \right) + \frac{1 - e^{-2rt}}{r} \text{trace}(MH(u)), $$

where $H(u) = \mathbb{E}[\nabla u(\mathbf{X})\nabla u(\mathbf{X})^\top]$ is the AS matrix of $u$ and $H(P_{t,M}(u)) = \mathbb{E}[\nabla P_{t,M}(u)(\mathbf{X})\nabla P_{t,M}(u)(\mathbf{X})^\top]$ the AS matrix of $P_{t,M}(u)$. The above expression is not easy to work with in practice, as we do not have direct access to the gradients of $P_{t,M}(u)$. The next proposition provides an alternative expression for $H(P_{t,M}(u))$ that facilitates its numerical approximation.

**Proposition 4.1.** Let $u \in L_{\mu}^2(\mathbb{R}^d)$ such that $\nabla u \in (L_{\mu}^2(\mathbb{R}^d))^d$, $t \ge 0$ and $M \in \mathbb{R}^{d \times d}$ a symmetric semi-definite positive matrix. We have

$$ (4.2) \qquad H(P_{t,M}(u)) = e^{-tM} \mathbb{E} \left[ \nabla u(\mathbf{Y}) \nabla u(\mathbf{Y}')^\top \right] e^{-tM}, $$

where $(\mathbf{Y}, \mathbf{Y}') \sim \mathcal{N}(0, \Gamma)$, where $\Gamma = \begin{pmatrix} I_d & e^{-2tM} \\ e^{-2tM} & I_d \end{pmatrix}$.

Based on (4.2), we can use the following Monte-Carlo estimator for $H(P_{t,M}(u))$:

$$ \hat{H}^N(P_{t,M}(u)) = \frac{1}{N} \sum_{i=1}^{N} e^{-tM} \left( \nabla u(\mathbf{Y}^{(i)}) \nabla u(\mathbf{Y}^{(i)})^\top \right)_{\text{sym}} e^{-tM}, $$

where $(\mathbf{Y}^{(1)}, \mathbf{Y}^{(1)}), \dots, (\mathbf{Y}^{(N)}, \mathbf{Y}^{(N)})$ are $N$ independent copies of $(\mathbf{Y}, \mathbf{Y}') \sim \mathcal{N}(0, \Gamma)$, and where $(A)_{\text{sym}} = (A + A^\top)/2$. However, if one disposes of a Monte-Carlo estimation

$$ (4.3) \qquad \hat{H}^N(u) = \frac{1}{N} \sum_{i=1}^{N} \nabla u(\mathbf{X}^{(i)}) \nabla u(\mathbf{X}^{(i)})^\top $$

of $H(u)$ with $\mathbf{X}^{(i)} \stackrel{iid}{\sim} \mathcal{N}(0, I_d)$, one would prefer the following estimator which permits to recycle the already-computed gradient snapshots $\nabla u(\mathbf{X}^{(i)})$:

$$ (4.4) \qquad \hat{H}^{N,N'}(P_{t,M}(u)) = e^{-tM} \frac{1}{N} \sum_{i=1}^{N} \left( \nabla u(\mathbf{X}^{(i)}) \frac{1}{N'} \sum_{j=1}^{N'} \nabla u(\mathbf{Y}^{(i,j)})^\top \right)_{\text{sym}} e^{-tM}, $$

where, for all $1 \le i \le N$, $\mathbf{Y}^{(i,j)}$ are $j$ i.i.d realization of $(\mathbf{Y}'|(\mathbf{Y} = \mathbf{X}^{(i)}) \sim \mathcal{N}(e^{-2tM}\mathbf{X}^{(i)}, I_d - e^{-4tM})$. Note that the above estimator requires $(N' + 1)N$ gradient evaluations.
**Algorithm 4.1 MAS procedure**

**Require:** Training samples $(\mathbf{X}^{(i)}, \nabla u(\mathbf{X}^{(i)}))_{1\le i\le N}$ where $\mathbf{X}^{(i)} \stackrel{\text{iid}}{\sim} \mathcal{N}(0, I_d)$

1: Compute the Active Subspace matrix estimate $\hat{H}^N(u) = \frac{1}{N} \sum_{i=1}^N \nabla u(\mathbf{X}^{(i)}) \nabla u(\mathbf{X}^{(i)})^\top$.

2: Compute the singular value decomposition of $\hat{H}^N(u) = \hat{U}\text{diag}(\hat{\lambda})\hat{U}^\top$ where $\hat{U} \in \mathbb{R}^{d \times d}$ is a unitary matrix and $\hat{\lambda} = (\hat{\lambda}_1, \dots, \hat{\lambda}_d)$ with $\hat{\lambda}_1 \ge \hat{\lambda}_2 \ge \dots \ge 0$.

3: For a given threshold $\epsilon > 0$, find the smallest $m_0 \in \mathbb{N}$ for which $\hat{\mathcal{J}}^N(\hat{U}_{m_0}, u) := \sum_{i=m_0+1}^d \hat{\lambda}_i \le \epsilon$.

4: Choose matrix $M \in \{I_d; \hat{H}^N(u); \hat{U}_{m_0}\hat{U}_{m_0}^\top\}$ and parameter $t$ as in (4.1) for $\alpha \in \{\frac{1}{8}, \frac{1}{4}, \frac{3}{8}, \frac{1}{2}\}$.

5: For $1 \le i \le N$ and $1 \le j \le N'$, sample $\mathbf{Y}^{(i,j)}$ where $\mathbf{Y}^{(i,j)} \stackrel{\text{iid}}{\sim} \mathcal{N}(e^{-2tM}\mathbf{X}^{(i)}, I_d - e^{-4tM})$

6: Compute the Monte Carlo estimate $\hat{H}^{N,N'}(P_{t,M}(u))$ using (4.4).

7: **if** $M = I_d$ **then**

8: Compute the singular value decomposition of $\hat{H}^{N,N'}(P_{t,M}(u)) = \hat{U}^{t,M}\text{diag}(\hat{\lambda}^{t,M})(\hat{U}^{t,M})^\top$ where $\hat{\lambda}_1^{t,M} \ge \hat{\lambda}_2^{t,M} \ge \dots \ge 0$ are the eigenvalues of $\hat{H}^{N,N'}(P_{t,M}(u))$.

9: For $m \le m_0$ compute $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_m^{t,M}, u) = \sum_{i=m+1}^d \hat{\lambda}_i^{t,M} + (1 - e^{-2t}) \sum_{i=1}^d \hat{\lambda}_i$ which is a minimum of the MAS bound.

10: **if** there exists $m \le m_0$ such that $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_m^{t,M}, u) \le \hat{\mathcal{J}}^N(\hat{U}_{m_0}, u)$ **then**

11: Set $g(x) = (\hat{U}_m^{t,M})^\top x$.

12: **else**

13: Keep $g(x) = (\hat{U}_{m_0})^\top x$.

14: **end if**

15: **else**

16: Compute the singular value decomposition of $M$ such that $M = V\text{diag}(\mu)V^\top$ where $V := [v_1, \dots, v_d] \in \mathbb{R}^{d \times d}$ is the orthogonal matrix which columns $v_i$ are eigenvectors of $M$ and where $\mu = (\mu_1, \dots, \mu_d)$ with $\mu_1 \ge \mu_2 \ge \dots \ge 0$.

17: For $m \le m_0$ and $U_m = \{U_m := [u_1, \dots, u_m] \in \mathbb{R}^{d \times m} \mid \{u_1, \dots, u_m\} \subseteq \{v_1, \dots, v_d\}\}$ compute $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_m^{t,M}, u) = \min_{U_m \in U_m} \text{trace}((I_d - U_m U_m^\top)\hat{H}^{N,N'}(P_{t,M}(u))) + \frac{1-e^{-2rt}}{r} \sum_{i=1}^d \mu_i \hat{\lambda}_i$, where $r = \min_{\mu_i > 0} \mu_i$.

18: **if** there exists $m \le m_0$ such that $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_m^{t,M}, u) \le \hat{\mathcal{J}}^N(\hat{U}_{m_0}, u)$ **then**

19: Set $g(x) = (\hat{U}_m^{t,M})^\top x$.

20: **else**

21: Keep $g(x) = (\hat{U}_{m_0})^\top x$.

22: **end if**

23: **end if**
**4.3. Procedure.** Algorithm 4.1 details the practical procedure of the mollified Active Subspace method. Steps 1 to 3 correspond to the usual Active Subspace dimension reduction strategy. The following steps allow to correct some of the potential selection errors.

**4.4. Illustration on an analytical example.** In this section, we study a simple analytical toy model in which selection errors may occur. This model allows us to highlight the advantages of our strategy compared to Active Subspace, as well as to compare different choices for $t$ and $M$ proposed in Subsection 4.1. For $a, b \in \mathbb{R}$, we set:

$$ (4.5) \qquad u(x_1, x_2) = ax_1 + \frac{b}{\sqrt{120}}(x_2^5 - 10x_2^3 + 15x_2) $$

This choice is convenient as it yields the closed-form expressions $e_{1,0}(x_1, x_2) = x_1$ and $e_{0,5}(x_1, x_2) = \frac{x_2^5 - 10x_2^3 + 15x_2}{\sqrt{120}}$, enabling direct computation of all relevant quantities. Using the orthogonality of Hermite polynomials, we obtain the following Active Subspace matrix:

$$ H(u) = \begin{pmatrix} a^2 & 0 \\ 0 & 25b^2 \end{pmatrix}. $$

Since $H(u)$ is already diagonal, its eigenvector basis is the canonical basis, and its eigenvalues are $a^2$ and $25b^2$. Let $U_{\{1\}}$ and $U_{\{2\}}$ be such that $[U_{\{1\}}, U_{\{2\}}] = I_2$, so that $U_{\{1\}}^{\top}x = x_1$ and $U_{\{2\}}^{\top}x = x_2$. The values of the reconstruction error and the Poincaré upper-bound are then given by:

<table><thead><tr><th></th><th>$U_m = U_{\{1\}}$</th><th>$U_m = U_{\{2\}}$</th></tr></thead><tbody><tr><td>$\mathcal{E}(U_m, u)$</td><td>$b^2$</td><td>$a^2$</td></tr><tr><td>$\mathcal{J}(U_m, u)$</td><td>$25b^2$</td><td>$a^2$</td></tr></tbody></table>

For $a = 1$ and $b = 1/6$, minimizing the upper bound correctly selects $x_1$, leading to the optimal reconstruction error. However, for $a = 1$ and $b = 1/4$, minimizing the upper bound leads to selecting $x_2$, which is the worst possible since $\mathcal{E}(U_{\{2\}}, u) = 1$ whereas $\mathcal{E}(U_{\{1\}}, u) = 1/16$. We now demonstrate how the mollified Active Subspace strategy mitigates this selection error.

Using Lemma 3.5, for $M = \text{diag}(m_1, m_2) \in \mathbb{R}^{2 \times 2}$ a diagonal matrix, we obtain

$$ \begin{aligned} P_{t,M}(u)(x) &= aP_{t,M}(e_{1,0})(x) + bP_{t,M}(e_{0,5})(x) \\ &= ae^{-m_1 t}e_{1,0}(x) + be^{-5m_2 t}e_{0,5}(x). \end{aligned} $$

Using the orthogonality of Hermite polynomials, we deduce

$$ H(P_{t,M}(u)) = \begin{pmatrix} a^2 e^{-2m_1 t} & 0 \\ 0 & 25b^2 e^{-10m_2 t} \end{pmatrix} $$

Since $H(P_{t,M}(u))$ is diagonal, its eigenvector basis is the canonical basis and its eigenvalues are $\lambda_1^t = a^2 e^{-2m_1 t}$ and $\lambda_2^t = 25b^2 e^{-10m_2 t}$. Setting $r = \min(m_1, m_2)$ we obtain the bounds:

$$ \begin{aligned} \mathcal{J}_{t,M}(U_{\{1\}}, u) &= 25b^2 e^{-10m_2 t} + \frac{1 - e^{-2rt}}{r}(a^2 m_1 + 25b^2 m_2) \\ \mathcal{J}_{t,M}(U_{\{2\}}, u) &= a^2 e^{-2m_1 t} + \frac{1 - e^{-2rt}}{r}(a^2 m_1 + 25b^2 m_2). \end{aligned} $$
<table><caption>Table 1: Values of $t_α$ (as defined in Subsection 4.1) for $m_0 = 1$ and for different choices of the matrix $M$ and the parameters $a, b$. Here we set $\epsilon = \mathcal{J}(U_{m_0}, u)$. Bold values correspond to values of $t_α$ which allows to select the right coordinate.</caption><thead><tr><th rowspan="2">t<sub>α</sub></th><th colspan="3">a = 1, b = 1/4</th><th colspan="3">a = 1, b = 1/6</th></tr><tr><th>M = I<sub>2</sub></th><th>M = H(u)</th><th>M = U<sub>m<sub>0</sub></sub>U<sup>T</sup><sub>m<sub>0</sub></sub></th><th>M = I<sub>2</sub></th><th>M = H(u)</th><th>M = U<sub>m<sub>0</sub></sub>U<sup>T</sup><sub>m<sub>0</sub></sub></th></tr></thead><tbody><tr><td>t<sub>1/2</sub></td><td>0.109</td><td><strong>0.078</strong></td><td><strong>0.347</strong></td><td><strong>0.115</strong></td><td><strong>0.128</strong></td><td><strong>0.213</strong></td></tr><tr><td>t<sub>3/8</sub></td><td>0.079</td><td>0.058</td><td><strong>0.235</strong></td><td><strong>0.083</strong></td><td><strong>0.094</strong></td><td><strong>0.151</strong></td></tr><tr><td>t<sub>1/4</sub></td><td>0.051</td><td>0.038</td><td><strong>0.144</strong></td><td><strong>0.054</strong></td><td><strong>0.061</strong></td><td><strong>0.095</strong></td></tr><tr><td>t<sub>1/8</sub></td><td>0.025</td><td>0.018</td><td><strong>0.067</strong></td><td><strong>0.026</strong></td><td><strong>0.030</strong></td><td><strong>0.045</strong></td></tr></tbody></table>

For this example, $M = I_2$, $M = H(u)$ or $M = U_{m_0}U_{m_0}^\top$ (where $U_{m_0}$ represents the original Active Subspace projector) respectively correspond to $(m_1, m_2) = (1, 1)$; $(m_1, m_2) = (a^2, 25b^2)$ or $(m_1, m_2) = (\mathbf{1}_{a^2>25b^2}, \mathbf{1}_{a^2\le 25b^2})$; where $\mathbf{1}(\cdot)$ is the indicator function.

Figure 4 shows the upper bounds computed using the Active Subspace (AS) and Mollified Active Subspace (MAS) methods as functions of $t$, for various choices of $M$ and different values of $a$ and $b$. In these plots, the blue lines represent the upper bound obtained by selecting the first variable, and the orange lines represent the one obtained by selection the second variable. The MAS method provides a better upper bound when one of the MAS curves goes below both AS bounds and the MAS approach is capable of correcting selection errors when $t$ lies within the gray-shaded region. In Figures 4a to 4c, the choice of $a$ and $b$ leads to a selection error. In this case, the MAS method successfully corrects the error when $M = U_{m_0}U_{m_0}^\top$ or $M = H(u)$. Among these, $M = U_{m_0}U_{m_0}^\top$ appears to be more effective, as it offers a wider range of $t$ values that yield correct variable selection, along with a tighter bound. In contrast, in Figures 4d to 4f, the choice of $a$ and $b$ does not result in a selection error; hence, both AS and MAS select the same variable (in this case, $x_1$).

Finally, Table 1 reports the values of $t_α$, as defined in Subsection 4.1 for $\epsilon = \mathcal{J}(U_{m_0}, u)$ and for various choices of $M$. Bold values indicate values of $t$ that allow to select the right coordinate. For this example, the choice $\alpha = \frac{1}{2}$ is optimal, as it allows to correct the selection error across all considered cases.

**5. Numerical examples.** Let $(\mathbf{X}^{(1)}, \dots, \mathbf{X}^{(N)})$ denote $N$ i.i.d samples of $\mathbf{X} \sim \mathcal{N}(0, I_d)$ and for $1 \le i \le N'$ let $(\mathbf{Y}^{(i,1)}, \dots, \mathbf{Y}^{(i,N')})$ denote $N'$ i.i.d samples of $(\mathbf{Y}'|(\mathbf{X} = \mathbf{X}^{(i)}) \sim \mathcal{N}(e^{-2tM}\mathbf{X}^{(i)}, I_d - e^{-4tM})$. For a given rank $m < d$, using the Monte Carlo estimator (4.3) and (4.4), we can estimate the AS and the MAS upper bounds as follows:

$$ \begin{aligned} \hat{\mathcal{J}}^N(U_m, u) &= \mathrm{tr}((I_d - U_m U_m^\top) \hat{\mathcal{H}}^N(u)) \\ \hat{\mathcal{J}}_{t,M}^{N,N'}(U_m, u) &= \mathrm{tr}((I_d - U_m U_m^\top) \hat{\mathcal{H}}^{N,N'}(P_{t,M}(u))) + \frac{1 - e^{-2rt}}{r} \mathrm{tr}(M \hat{\mathcal{H}}^N(u)). \end{aligned} $$

Minimizing these estimators over $U_m$, as outlined in Algorithm 4.1, allows to compute the feature map $g(x) = U_m^\top x$. The corresponding profile function $f$ is modeled by a fully connected neural network with three hidden layers, each containing 20 neurons and ReLU activation
<figure><img src="image_7.png" /><figcaption>Figure 4: $\mathcal{J}(U_m, u)$ (AS) and $\mathcal{J}_{t,M}(U_m, u)$ (MAS) according to $t$ and for different values of $M$. The first line corresponds to the case $a = 1$ and $b = 1/4$ and the second line corresponds to the case $a = 1$ and $b = 1/6$. Here $U_{\{1\}}$ is the projector onto $x_1$ and $U_{\{2\}}$ is the projector onto $x_2$.</figcaption></figure>

functions. For a given feature map $g$, the network is trained by minimizing the empirical mean squared error:

$$ \frac{1}{n_x} \sum_{i=1}^{n_x} (u(x_i) - f \circ g(x_i))^2 $$

using stochastic gradient descent with the Adam optimizer.

**5.1. Sum of sines.** In this section, we consider the sum-of-sines model, which is known to be prone to selection error for certain frequency configurations $\omega$ [30]. Given vectors $a = (a_1, \dots, a_d) \in \mathbb{R}^d$ and $\omega = (\omega_1, \dots, \omega_d) \in \mathbb{R}^d$, the model is defined as:

$$ u(x_1, \dots, x_d) = \sum_{i=1}^{d} a_i \sin(\omega_i x_i) $$

**Example from Section 2.** We compare the Active Subspace (AS) and Mollified Active Subspace (MAS) methods on the toy example introduced in Section 2:

$$ u_2: \begin{array}{l} \mathbb{R}^2 \\ (x_1, x_2) \end{array} \rightarrow \begin{array}{l} \mathbb{R} \\ \sin(x_1) + \frac{1}{6} \sin(10x_2), \end{array} $$
<table><caption>Table 2: Mean ± standard deviation of the approximation error and upper bound for AS and MAS on the model $u_2$. MAS results are reported for different choices of the matrix $M$ and mollification parameter $t$. The upper bound column corresponds to $\hat{\mathcal{J}}^N(\hat{U}_{m_0}, u)$ for AS and to $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_{m_0}^{t,M}, u)$ for MAS, as defined in Algorithm 4.1. The approximation error column reports the empirical reconstruction error between $u$ and $f \circ g$, evaluated on a test set of 1000 samples.</caption><thead><tr><th>Method</th><th><i>M</i></th><th><i>t</i></th><th>Approximation Error</th><th>Upper Bound</th></tr></thead><tbody><tr><td>AS</td><td>–</td><td>–</td><td>0.446 ± 0.004</td><td>0.562 ± 0.021</td></tr><tr><td rowspan="3">MAS</td><td rowspan="3"><i>I</i><sub>2</sub></td><td><i>t</i><sub>1/2</sub></td><td>0.018 ± 0.004</td><td>0.338 ± 0.097</td></tr><tr><td><i>t</i><sub>3/8</sub></td><td>0.017 ± 0.003</td><td>0.164 ± 0.074</td></tr><tr><td><i>t</i><sub>1/4</sub></td><td>0.035 ± 0.031</td><td>0.136 ± 0.315</td></tr><tr><td rowspan="3">MAS</td><td rowspan="3"><i>Î</i><sup><i>N</i></sup>(<i>u</i>)</td><td><i>t</i><sub>1/2</sub></td><td>0.016 ± 0.001</td><td>0.267 ± 0.077</td></tr><tr><td><i>t</i><sub>3/8</sub></td><td>0.016 ± 0.001</td><td>0.211 ± 0.123</td></tr><tr><td><i>t</i><sub>1/4</sub></td><td>0.016 ± 0.002</td><td>0.062 ± 0.202</td></tr><tr><td rowspan="3">MAS</td><td rowspan="3"><i>Î</i><sub><i>m</i><sub>0</sub></sub><i>Î</i><sub><i>m</i><sub>0</sub></sub><sup>†</sup></td><td><i>t</i><sub>1/2</sub></td><td>0.014 ± 0.001</td><td>0.222 ± 0.103</td></tr><tr><td><i>t</i><sub>3/8</sub></td><td>0.015 ± 0.001</td><td>0.243 ± 0.114</td></tr><tr><td><i>t</i><sub>1/4</sub></td><td>0.015 ± 0.001</td><td>0.123 ± 0.187</td></tr></tbody></table>

which corresponds to a particular instance of the sum-of-sines model with $d = 2$, $a = (1, 1/6)$, and $\omega = (1, 10)$.

Since $d = 2$, we fix the target dimension to $m_0 = 1$ and consider only rank-1 projectors. In Table 2, we report the mean and standard deviation errors computed for 10 independent runs for both AS and MAS. Each training dataset consists of $N = 500$ and $N' = 10$ samples. The upper bound column corresponds to $\hat{\mathcal{J}}^N(\hat{U}_{m_0}, u)$ for AS and to $\hat{\mathcal{J}}_{t,M}^{N,N'}(\hat{U}_{m_0}^{t,M}, u)$ for MAS, as defined in Algorithm 4.1. The approximation error column reports the empirical reconstruction error between $u$ and $f \circ g$, evaluated on a test set of 1000 samples. For the MAS method, we consider three choices for the matrix $M$: the identity matrix $M = I_2$, the estimated AS matrix $M = \hat{\mathcal{H}}^N(u)$, and the rank-1 projection $M = \hat{U}_{m_0}\hat{U}_{m_0}^\top$, where $\hat{U}_{m_0}$ is the first active direction estimated using $\hat{\mathcal{H}}^N(u)$. The parameter $t_\alpha$ is selected according to the procedure described in Subsection 4.1, using $\epsilon = \hat{\mathcal{J}}^N(\hat{U}_{m_0}, u)$.

The results shown in Table 2 demonstrate that the MAS method consistently corrects the selection error observed in AS, regardless the choice of $M$ and $t$. In all cases, MAS achieves a near-optimal reconstruction error (see Section 2 for theoretical justification).

**Example in dimension d=8.** We compare the Active Subspace (AS) and Mollified Active Subspace (MAS) methods on the 8-dimensional sum-of-sines model $u_8$ with parameters:

$$ 
\begin{aligned}
a &= \left( 1, \frac{1}{2}, \frac{1}{3}, \frac{1}{4}, \frac{1}{6}, \frac{1}{8}, \frac{1}{10}, \frac{1}{10} \right), \\
\omega &= (1, 4, 5, 6, 7, 9, 2, 7).
\end{aligned}
$$
Figure 5 shows the AS and MAS upper bounds as a function of the lattent dimension m, for different choices of M and sample size $n_y$. Figure 6 and Figure 7 display the corresponding approximation error, defined as the empirical mean squared error between $u_8$ and $f \circ g$, evaluated on a test set of 10,000 samples. In Figure 5 the solid line represents the mean value over 10 runs of the algorithm, while the shaded area denotes the regions defined as [mean ± standard deviation]. In Figure 6 and Figure 7, the solid lines correspond to independent runs of the algorithm, each training dataset for the MAS method consists of $N = 500$ samples and either $N' = 10$ or $N' = 3$, depending on the experiment. The matrix M is set to $I_8$, $\hat{H}^N(u)$, or $\hat{U}_{m_0} \hat{U}_{m_0}^\top$, where $\hat{U}_{m_0}$ is the matrix which columns are the first $m_0 = 6$ active directions.

The results in Figures 5 and 6 indicate that MAS consistently outperforms AS for every choice of M, achieving significantly lower upper bounds and approximation errors up to $m = 5$. The best approximation errors are observed for $m = 3$ or $m = 4$, as increasing m makes learning the profile function f more challenging, even when the feature map g is well estimated. We note that even with $N' = 3$ samples, MAS achieves substantial improvements over AS. According to Figure 6, the best choice for M seems to be $M = I_8$ for this example. In Figure 7 we compare AS and MAS with a similar computational budget: we compute g using 4000 training points for the AS method and $N = 1000$, $N' = 3$ for the MAS method (which also results in 4000 training points). The profile function f is trained using 4000 training points for AS and 1000 training points for MAS. Figure 7 shows that MAS for $M = I_8$ and $M = \hat{H}^N(u)$ provides a lower approximation error than AS for $m \in \{3, 4, 5\}$ while preserving similar performances for $m \in \{1, 2\}$.

## 5.2. Wave propagation model.
In this section we consider a one-dimensional wave equation as a simplified model for acoustic wave propagation in heterogeneous media. The scalar displacement field $v(s,t)$ is given by the following partial differential equation with homogeneous Dirichlet boundary conditions:

$$ (5.1) \qquad \begin{aligned} \frac{\partial^2 v}{\partial t^2}(s,t) &= \frac{\partial}{\partial s} \left( c^2(s;X) \frac{\partial v}{\partial s}(s,t) \right), \quad s \in (0,1), t \in (0,T] \\ v(0,t) &= v(1,t) = 0, \quad t \in (0,T] \end{aligned} $$

where $X \in [0, 1]^{2q}$ is a vector of parameters and where $c(s;X)$ denotes the spatially varying wave speed. We consider initial conditions corresponding to a localized initial pulse in $s = 0.5$:

$$ v(s, 0) = \exp(-100(s - 0.5)^2), \quad \frac{\partial v}{\partial t}(s, 0) = 0. $$

The wave speed is modeled as:

$$ c(s; X) = 1 + \sum_{k=0}^{q} \frac{c_k - 0.5}{2^{k+1}} \sin(10^k \omega_k s), $$

where $X = (c_1, \dots, c_q, \omega_1, \dots, \omega_q) \in [0, 1]^{2q}$. This formulation incorporates Fourier coefficients and frequencies as parameters of the Fourier expansion of the wave speed. This allows to select key Fourier parameters for approximating a given quantity of interest (QoI). Here, the QoI we consider is the wave field value at a fixed spatial location $s_0 \in (0, 1)$ and final time $T$, i.e.,
<figure><img src="image_8.png" /><figcaption>Figure 5: Active Subspace (AS) and Mollified Active Subspace (MAS) upper bound as function of the latent dimension $m$ for $u_8$. The solid line represents the mean value over 10 runs, while the shaded area denotes the region of mean ± standard deviation. Here $N = 500$ and the results are shown for different choices of mollifying parameters $t$, $M$ and samples size $N'$.</figcaption></figure>

$u(s_0, T)$. In particular, we set $q = 4$, $s_0 = 0.7$ and $T = 1$ and we consider the 8 dimensional computational model $u_{\text{wave}}(\mathbf{X}) = v(s_0, T, \mathbf{X})$ where $\mathbf{X} \sim \mathcal{U}([0, 1]^8)$. Because $\mathbf{X} \sim \mathcal{U}([0, 1]^8)$, we apply our algorithm to the normalized version of the model $\bar{u}_{\text{wave}}(\mathbf{Y}) = u_{\text{wave}} \circ F_{\mathcal{N}(0,1)}(\mathbf{Y})$ where $\mathbf{Y} \sim \mathcal{N}(0, I_d)$ and where $F_{\mathcal{N}(0,1)}(\mathbf{Y})$ denotes the cumulative distribution function of $\mathcal{N}(0, 1)$ applied elementwise to $\mathbf{Y}$. We solve (5.1) using the finite element method [12] and the gradient is approximated via finite differences [18] with a step size of $10^{-3}$.

Figure 5 shows the AS and MAS upper bounds and approximation errors as a function of the latent dimension $m$, for different sample size $N$. Here $M = I_8$ and $N' = 5$, the approximation error is defined as the empirical mean squared error between $u_{\text{wave}}$ and $f \circ g$, evaluated on a test set of 1000 samples. The parameter $t$ is calibrated using expression (4.1) with $m_0 = 4$, $\alpha = 1/2$ and $\epsilon = \mathcal{J}^N(\hat{U}_{m_0}, u)$. In both figures, the solid lines correspond to independent runs of the algorithm. The MAS algorithm consistently achieves a lower upper bound and a lower reconstruction error than AS up to $m = 5$. In Figures 8b and 8d we compare AS and MAS with a similar computational budget: we compute $g$ using 3000 training points for the AS method and $N = 500$, $N' = 5$ for the MAS method (which also result in 3000 training points). The profile function $f$ is trained using 3000 training points for AS and 500 for MAS.
<figure><img src="image_9.png" /><figcaption>Figure 6: Active Subspace (AS) and Mollified Active Subspace (MAS) approximation error as a function of the latent dimension *m* for *u*<sub>8</sub>. Each solid line represent one independent run. Results are shown for different choices of the mollification matrix *M* and sample size *N*′.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Active Subspace (AS) and Mollified Active Subspace (MAS) approximation error as a function of the latent dimension *m* for *u*<sub>8</sub> and for the same computational budget. Here we use 4000 training points for AS and *N* = 1000, *N*′ = 3 for MAS which also results in 4000 training points. Each solid line represent one independent run.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Error bound AS $N = 500$, MAS $N = 500$, $N' = 5$ (b) Error bound AS $N = 3000$, MAS $N = 500$, $N' = 5$</figcaption></figure>

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

<figure><img src="image_13.png" /><figcaption>(c) Approximation error $N = 500$, MAS $N = 500$ (d) Approximation error $N = 3000$, MAS $N = 500$, $N' = 5$</figcaption></figure>

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

Figure 8: Active Subspace (AS) and Mollified Active Subspace (MAS) upper bound and approximation error as a function of the latent dimension $m$ for the 8 dimensional ($d = 2q$) wave propagation model. Here $M = I_8$, $N = 500$ and $N' = 5$ for the MAS algorithm. Each solid line represents a run of the algorithm. Figures from the right column compare MAS to AS with a similar budget of 3000 points in total.

training points for MAS. Even in this case MAS provides better results than AS up to $m = 5$.

**6. Conclusion and perspectives.** In this work, we address a key limitation of the Active Subspace (AS) method: its sensitivity to high-frequency, low-amplitude oscillations in the model output. Such oscillations can make gradient evaluations non-informative, leading to suboptimal feature maps and degraded surrogate performance. To overcome this issue, we proposed the Mollified Active Subspace (MAS) method, which introduces an anisotropic smoothing operator to attenuate noisy oscillatory components while preserving the model's essential structure. Our theoretical analysis of MAS yields a computable upper bound on the reconstruction error based on the mollified model. Although this bound provably improves upon the AS bound under restrictive assumptions, our numerical results show that these assumptions are not necessary for the proposed MAS to outperform AS. Further research is
needed to characterize the class of models for which MAS consistently outperforms AS. While we suggest practical parameter choices for the mollification operator (specifically, the parameters $t$ and $M$), determining their optimal values remains an open problem. Nonetheless, numerical experiments demonstrate that MAS achieves more reliable feature selection than standard AS in models with complex, oscillatory behavior for $M = I_d$. This approach opens several avenues for future research, including extending the method to non-Gaussian settings, as well as adapting it to gradient-free scenarios by leveraging the mollification trick to estimate model gradients [22, 17].

**Acknowledgments.** R. Verdière and O. Zahm acknowledge support from the ANR JCJC project MODENA (ANR-21-CE46-0006-01). The integrality of the authors acknowledge support from the consortium in Applied Mathematics CIROQUO, gathering partners in technological research and academia in the development of advanced methods for Computer Experiments. C. Prieur and O. Zahm have benefited from a French government grant managed by the Agence Nationale de la Recherche under the France 2030 program, reference ANR-23-IACL-0006. C. Prieur would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Representing, calibrating & leveraging prediction uncertainty from statistics to machine learning, where part of the work on this paper was undertaken. This work was supported by EPSRC grant EP/Z000580/1.

## REFERENCES

[1] D. BAKRY, I. GENTIL, AND M. LEDOUX, *Analysis and geometry of Markov diffusion operators*, vol. 348, Springer Science & Business Media, 2013.

[2] A. G. BAYDIN, B. A. PEARLMUTTER, A. A. RADUL, AND J. M. SISKIND, *Automatic differentiation in machine learning: a survey*, Journal of Marchine Learning Research, 18 (2018), pp. 1–43.

[3] A. BIETTI, J. BRUNA, AND L. PILLAUD-VIVIEN, *On learning gaussian multi-index models with gradient flow*, 2023, https://arxiv.org/abs/2310.19793, https://arxiv.org/abs/2310.19793.

[4] A. BIETTI, J. BRUNA, AND L. PILLAUD-VIVIEN, *On learning gaussian multi-index models with gradient flow*, 2023, https://arxiv.org/abs/2310.19793, https://arxiv.org/abs/2310.19793.

[5] D. BIGONI, Y. MARZOUK, C. PRIEUR, AND O. ZAHM, *Nonlinear dimension reduction for surrogate modeling using gradient information*, Information and Inference: A Journal of the IMA, 11 (2022), pp. 1597–1639, https://doi.org/10.1093/imaiai/iaac006, https://doi.org/10.1093%2Fimaiai%2Fiaac006.

[6] V. I. BOGACHEV, *Gaussian measures*, no. 62, American Mathematical Soc., 1998.

[7] W. BRYC, *Normal distributions*, in The Normal Distribution: Characterizations with Applications, Springer, 1995, pp. 23–38.

[8] Q. CHEN, E. ARNAUD, R. BAPTISTA, AND O. ZAHM, *Coupled input-output dimension reduction: Application to goal-oriented bayesian experimental design and global sensitivity analysis*, 2025, https://arxiv.org/abs/2406.13425, https://arxiv.org/abs/2406.13425.

[9] P. G. CONSTANTINE, E. DOW, AND Q. WANG, *Active subspace methods in theory and practice: Applications to kriging surfaces*, SIAM Journal on Scientific Computing, 36 (2014), p. A1500–A1524, https://doi.org/10.1137/130916138, http://dx.doi.org/10.1137/130916138.

[10] C. CUI, K. ZHANG, T. DAULBAEV, J. GUSAK, I. OSELEDETS, AND Z. ZHANG, *Active subspace of neural networks: Structural analysis and universal attacks*, SIAM Journal on Mathematics of Data Science, 2 (2020), pp. 1096–1122.

[11] S. DA VEIGA, F. GAMBOA, B. IOOSS, AND C. PRIEUR, *Basics and trends in sensitivity analysis: Theory and practice in R*, SIAM, 2021.

[12] A. ERN AND J.-L. GUERMOND, *Theory and practice of finite elements*, vol. 159, Springer, 2004.

[13] A. GRIEWANK ET AL., *On automatic differentiation*, Mathematical Programming: recent developments
and applications, 6 (1989), pp. 83–107.

[14] M. HRISTACHE, A. JUDITSKY, J. POLZEHL, AND V. SPOKOINY, *Structure adaptive approach for dimension reduction*, Annals of Statistics, (2001), pp. 1537–1566.

[15] J. L. JEFFERSON, J. M. GILBERT, P. G. CONSTANTINE, AND R. M. MAXWELL, *Active subspaces for sensitivity analysis and dimension reduction of an integrated hydrologic model*, Computers & geosciences, 83 (2015), pp. 127–138.

[16] T. KLOCK, A. LANTERI, AND S. VIGOGNA, *Estimating multi-index models with response-conditional least squares*, Electronic Journal of Statistics, 15 (2021), pp. 589–629, https://doi.org/10.1214/20-EJS1785, https://doi.org/10.1214/20-EJS1785.

[17] D. KOZAK, S. BECKER, A. DOOSTAN, AND L. TENORIO, *A stochastic subspace approach to gradient-free optimization in high dimensions*, Computational Optimization and Applications, 79 (2021), p. 339–368, https://doi.org/10.1007/s10589-021-00271-w, http://dx.doi.org/10.1007/s10589-021-00271-w.

[18] R. J. LEVEQUE, *Finite difference methods for differential equations*, Draft version for use in AMath, 585 (1998), p. 112.

[19] T. W. LUKACZYK, P. CONSTANTINE, F. PALACIOS, AND J. J. ALONSO, *Active subspaces for shape optimization*, in 10th AIAA multidisciplinary design optimization conference, 2014, p. 1171.

[20] N. LÜTHEN, O. ROUSTANT, F. GAMBOA, B. IOOSS, S. MARELLI, AND B. SUDRET, *Global sensitivity analysis using derivative-based sparse poincaré chaos expansions*, International Journal for Uncertainty Quantification, 13 (2023).

[21] A. MOUSAVI-HOSSEINI, D. WU, AND M. A. ERDOGDU, *Learning multi-index models with neural networks via mean-field langevin dynamics*, 2025, https://arxiv.org/abs/2408.07254, https://arxiv.org/abs/2408.07254.

[22] S. OSHER, H. HEATON, AND S. WU FUNG, *A hamilton-jacobi-based proximal operator*, Proceedings of the National Academy of Sciences, 120 (2023), https://doi.org/10.1073/pnas.2220469120, http://dx.doi.org/10.1073/pnas.2220469120.

[23] M. T. PARENTE, J. WALLIN, AND B. WOHLMUTH, *Generalized bounds for active subspaces*, Electronic Journal of Statistics, 14 (2020), pp. 917–943.

[24] R.-E. PLESSIX, *A review of the adjoint-state method for computing the gradient of a functional with geophysical applications*, Geophysical Journal International, 167 (2006), pp. 495–503.

[25] Y. REN, E. NICHANI, D. WU, AND J. D. LEE, *Emergence and scaling laws in sgd learning of shallow neural networks*, 2025, https://arxiv.org/abs/2504.19983, https://arxiv.org/abs/2504.19983.

[26] O. ROUSTANT, F. BARTHE, AND B. IOOSS, *Poincaré inequalities on intervals–application to sensitivity analysis*, (2017).

[27] I. M. SOBOL AND S. KUCHERENKO, *Derivative based global sensitivity measures*, Procedia-Social and Behavioral Sciences, 2 (2010), pp. 7745–7746.

[28] R. VERDIÈRE, C. PRIEUR, AND O. ZAHM, *Diffeomorphism-based feature learning using Poincaré inequalities on augmented input space*, working paper or preprint, Dec. 2023, https://hal.science/hal-04364208.

[29] G. YUAN, M. XU, S. KPOTUFE, AND D. HSU, *Efficient estimation of the central mean subspace via smoothed gradient outer products*, 2024, https://arxiv.org/abs/2312.15469, https://arxiv.org/abs/2312.15469.

[30] O. ZAHM, P. CONSTANTINE, C. PRIEUR, AND Y. MARZOUK, *Gradient-based dimension reduction of multivariate vector-valued functions*, (2019), https://arxiv.org/abs/1801.07922, https://arxiv.org/abs/1801.07922.

[31] O. ZAHM, T. CUI, K. LAW, A. SPANTINI, AND Y. MARZOUK, *Certified dimension reduction in nonlinear bayesian inverse problems*, Mathematics of Computation, 91 (2022), pp. 1789–1835.

[32] G. ZHANG, J. ZHANG, AND J. HINKLE, *Learning nonlinear level sets for dimensionality reduction in function approximation*, 2019, https://doi.org/10.48550/ARXIV.1902.10652, https://arxiv.org/abs/1902.10652.
**Appendix A. Proof of Lemma 3.3.** Let $\mathbf{Y} = U^\top \mathbf{X}$ and $z = U_m^\top x$. We have

$$ 
\begin{aligned}
\mathbb{E}[u(\mathbf{X})|U_m^\top \mathbf{X} = U_m^\top x] &\stackrel{(3.6)}{=} \mathbb{E}\left[\sum_{\alpha \in \mathbb{N}^d} u_\alpha^U e_\alpha(U^\top \mathbf{X})|U_m^\top \mathbf{X} = U_m^\top x\right] \\
&= \sum_{\alpha \in \mathbb{N}^d} u_\alpha^U \mathbb{E}[e_\alpha(\mathbf{Y})|\mathbf{Y} \le m = z] \\
&\stackrel{(3.5)}{=} \sum_{\alpha \in \mathbb{N}^d} u_\alpha^U \left(\prod_{i=1}^m h_{\alpha_i}(z_i)\right) \left(\prod_{i=m+1}^d \mathbb{E}[h_{\alpha_i}(\mathbf{Y}_i)]\right) \\
&= \sum_{\substack{\alpha \in \mathbb{N}^d \text{ s.t.} \\ \alpha_{m+1}, \dots, \alpha_d = 0}} u_\alpha^U \left(\prod_{i=1}^m h_{\alpha_i}(z_i)\right) = \sum_{\substack{\alpha \in \mathbb{N}^d \text{ s.t.} \\ \alpha_{m+1}, \dots, \alpha_d = 0}} u_\alpha^U e_\alpha(U^\top x)
\end{aligned}
$$

where we used the fact that $\mathbb{E}[h_{\alpha_i}(\mathbf{Y}_i)] = 1$ when $\alpha_i = 0$ and $\mathbb{E}[h_{\alpha_i}(\mathbf{Y}_i)] = 0$ otherwise.

**Appendix B. Proof of Lemma 3.4.** Equation (3.8) is a direct consequence of (3.6) and (3.7). To show (3.9), let $v : y \mapsto u(Uy)$ and $\mathbf{Y} = U\mathbf{X}$. We have $\nabla v(y) = U^\top \nabla u(Uy)$ so that

$$ 
\begin{aligned}
\mathcal{J}(U_m, u) &= \mathbb{E}[\|(I_d - U_m U_m^\top) \nabla u(\mathbf{X})\|_2^2] \\
&= \mathbb{E}[\|(I_d - U_m U_m^\top) U \nabla v(\mathbf{Y})\|_2^2] \\
&= \sum_{i=m+1}^d \mathbb{E}[(\partial_i v(\mathbf{Y}))^2].
\end{aligned}
$$

By (3.6), we have $v(y) = u(Uy) = \sum_{\alpha \in \mathbb{N}^d} u_\alpha^U e_\alpha(y)$ so that

$$ 
\begin{aligned}
\mathbb{E}[(\partial_i v(\mathbf{Y}))^2] &= \mathbb{E}\left[\left(\sum_{\alpha \in \mathbb{N}^d} u_\alpha^U e_\alpha(\mathbf{Y})\right)^2\right] \\
\text{((Pb) & (3.5))} &= \mathbb{E}\left[\left(\sum_{\alpha \in \mathbb{N}^d} u_\alpha^U \sqrt{\alpha_i} e_{\alpha-i}(\mathbf{Y})\right)^2\right] \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U \sqrt{\alpha_i})^2 \mathbb{E}[e_{\alpha-i}(\mathbf{Y})^2] \\
\text{(B.1)} \qquad &= \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \alpha_i.
\end{aligned}
$$

We deduce that $\mathcal{J}(U_m, u) = \sum_{i=m+1}^d \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 \alpha_i = \sum_{\alpha \in \mathbb{N}^d} (u_\alpha^U)^2 |\alpha_\perp|$, which gives (3.9).

**Appendix C. Proof of Lemma 3.5.** Let $y = Q^\top x$. By (3.6) with $U = Q$, we have
$u(x) = \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q e_{\alpha}(Q^{\top}x)$ so that

$$
\begin{align*}
P_{t,M}(u)(x) &= \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q P_{t,M} e_{\alpha}(Q^{\top}x) \\
&\stackrel{(3.1)}{=} \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q \mathbb{E} \left[ e_{\alpha} \left( Q^{\top} (e^{-tM}x + \sqrt{I_d - e^{-2tM}Z}) \right) \right] \\
(M = Q\mathrm{diag}(\lambda)Q^{\top}) &= \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q \mathbb{E} \left[ e_{\alpha} \left( e^{-t\mathrm{diag}(\lambda)}y + \sqrt{(I_d - e^{-2t\mathrm{diag}(\lambda)}Z)} \right) \right] \\
&\stackrel{(3.5)}{=} \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q \prod_{i=1}^d \frac{1}{\sqrt{\alpha_i!}} \mathbb{E} \left[ h_{\alpha_i} \left( e^{-t\lambda_i}y_i + \sqrt{(1 - e^{-2t\lambda_i})Z_i} \right) \right] \\
&\stackrel{(Pc)}{=} \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q \prod_{i=1}^d \frac{1}{\sqrt{\alpha_i!}} e^{-\alpha_i \lambda_i t} h_{\alpha_i}(y_i) \\
&\stackrel{(3.5)}{=} \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q e^{-\alpha^{\top} \lambda t} e_{\alpha}(y) = \sum_{\alpha \in \mathbb{N}^d} u_{\alpha}^Q e^{-\alpha^{\top} \lambda t} e_{\alpha}(Q^{\top}x).
\end{align*}
$$

This concludes the proof.

**Appendix D. Proof of Proposition 3.7.** We can write

$$
\begin{align*}
\mathcal{J}_{t,M}(U_m, u) &= \mathcal{J}(U_m, P_{t,m}u) + \frac{1 - e^{-2rt}}{r} \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 (\alpha^{\top} \lambda) \\
(\text{Equations (3.9}&\&(3.10)) \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| e^{-2\alpha^{\top}\lambda t} + \frac{1 - e^{-2rt}}{r} \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 (\alpha^{\top} \lambda) \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| (1 - 2\alpha^{\top}\lambda t) + 2t \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 (\alpha^{\top} \lambda) + \mathcal{O}(t^2) \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| - 2t \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 \alpha^{\top} \lambda (|\alpha_{\perp}| - 1) + \mathcal{O}(t^2),
\end{align*}
$$

which gives (3.13). Assume now that $M = I - U_m U_m^{\top}$ so that $\lambda_i = 0$ when $i \le m$, $\lambda_i = 1$ otherwise. Then we have

$$
\begin{align*}
\mathcal{J}_{t,M}(U_m, u) &= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| e^{-2\alpha^{\top}\lambda t} + \frac{1 - e^{-2rt}}{r} \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 (\alpha^{\top}\lambda) \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| e^{-2|\alpha_{\perp}|t} + (1 - e^{-2t}) \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| \\
&= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}| + \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 \underbrace{|\alpha_{\perp}|(e^{-2|\alpha_{\perp}|t} - e^{-2t})}_{\le 0} \\
&\le \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^U)^2 |\alpha_{\perp}|
\end{align*}
$$
which is (3.14)

## Appendix E. Proposition E.1.

**Proposition E.1.** Let $u \in L_{\mu}^2(\mathbb{R}^d)$ such that $\nabla u \in (L_{\mu}^2(\mathbb{R}^d))^d$ and let $U_m \in \mathbb{R}^{d \times m}$ be a matrix with orthonormal columns. Furthermore, let $t \ge 0$ and $M \in \mathbb{R}^{d \times d}$ be a symmetric semi-definite positive matrix. Then we have

$$ (E.1) \qquad \mathcal{E}(U_m, u) \le 2\mathcal{J}(U_m, P_{t,M}u) + \min\left\{t; \frac{2}{r}\right\} \mathbb{E}[\|\nabla u(\mathbf{X})\|_M^2], $$

where $r > 0$ is the smallest nonzero eigenvalue of $M$, $\mathbf{X} \sim \mu = \mathcal{N}(0, I_d)$ and $\|\nabla u(\mathbf{X})\|_M^2 = \nabla u(\mathbf{X})^{\top} M \nabla u(\mathbf{X})$.

*Proof.* Consider the eigendecomposition of $M = Q\mathrm{diag}(\lambda)Q^{\top}$. We can write

$$ \begin{aligned} \mathcal{E}(U_m, u) &\le \mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[P_{t,M}(u)(\mathbf{X})|U_m^{\top}\mathbf{X}])^2] \\ &\le 2\mathbb{E}[(u(\mathbf{X}) - P_{t,M}(u)(\mathbf{X}))^2] + 2\mathbb{E}[(P_{t,M}(u)(\mathbf{X}) - \mathbb{E}[P_{t,M}(u)(\mathbf{X})|U_m^{\top}\mathbf{X}])^2] \\ &\le 2\mathbb{E}[(u(\mathbf{X}) - P_{t,M}(u)(\mathbf{X}))^2] + 2\mathcal{J}(U_m, P_{t,m}u) \end{aligned} $$

Next, using the expression (3.10) of $P_{t,M}(u)$, we have

$$ \begin{aligned} \mathbb{E}[(u(\mathbf{X}) - P_{t,M}(u)(\mathbf{X}))^2] &= \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^Q)^2 (1 - e^{-\alpha^{\top} \lambda t})^2 \\ &\le h_{\lambda}(t)^2 \sum_{\alpha \in \mathbb{N}^d} (u_{\alpha}^Q)^2 \alpha^{\top} \lambda \\ &\stackrel{(3.12)}{=} h_{\lambda}(t)^2 \mathbb{E}[\|\nabla u(\mathbf{X})\|_M^2] \end{aligned} $$

where $h_{\lambda}(\cdot)$ is defined by

$$ h_{\lambda}(t) = \sup_{\substack{\alpha \in \mathbb{N}^{d} \text{ s.t.} \\ \alpha^{\top} \lambda \neq 0}} \frac{1 - e^{-\alpha^{\top} \lambda t}}{\sqrt{\alpha^{\top} \lambda}} $$

It remains to show $h_{\lambda}(t)^2 \le \min\{t/2; 1/r\}$ in order to obtain the result (E.1). First, we have

$$ h_{\lambda}(t) = \sqrt{t} \sup_{\substack{\alpha \in \mathbb{N}^{d} \text{ s.t.} \\ \alpha^{\top} \lambda \neq 0}} \frac{1 - e^{-\alpha^{\top} \lambda t}}{\sqrt{\alpha^{\top} \lambda t}} \le \sqrt{t} \sup_{x>0} \frac{1 - e^{-x}}{\sqrt{x}} \le \sqrt{t/2}, $$

where we used the fact that $\frac{(1-e^{-x})^2}{x} \le \frac{1}{2}$ for all $x > 0$. Furthermore, by definition of $r$, we have $\alpha^{\top} \lambda \ge r$ for all $\alpha \in \mathbb{N}^d$ such that $\alpha^{\top} \lambda \neq 0$ so that

$$ h_{\lambda}(t) \le \sup_{\substack{\alpha \in \mathbb{N}^{d} \text{ s.t.} \\ \alpha^{\top} \lambda \neq 0}} \frac{1}{\sqrt{\alpha^{\top} \lambda}} = \frac{1}{\sqrt{r}}. $$

This yields $h_{\lambda}(t) \le \min\{\sqrt{t/2}; 1/\sqrt{r}\}$ and concludes the proof.

<figure><img src="image_15.png" /></figure>
Compared to (3.11), the bound in (E.1) is more general (because doesn't require any assumption on $M$), but it might be looser. For instance when $t \to 0$, the right-hand side of (E.1) goes to $2\mathcal{J}(U_m, u)$, which is twice the AS bound (2.1).

**Appendix F. Proof of Proposition 4.1.** Using the definition of $P_{t,M}(u)$ we obtain:

$$H(P_{t,M}(u)) = \mathbb{E}[\nabla P_{t,M}(u)(\mathbf{X})\nabla P_{t,M}(u)(\mathbf{X})^{\top}] \\ = e^{-tM}\mathbb{E}_{\mathbf{X},\mathbf{Z},\mathbf{Z}'}[\nabla u(e^{-tM}\mathbf{X} + \sqrt{I_d - e^{-2tM}\mathbf{Z}})\nabla u(e^{-tM}\mathbf{X} + \sqrt{I_d - e^{-2tM}\mathbf{Z}'})^{\top}]e^{-tM}$$

where $\mathbf{Z}, \mathbf{Z}' \sim \mathcal{N}(0, I_d)$ and $\mathbf{Z}, \mathbf{Z}'$ and $\mathbf{X}$ are mutually independent. We denote :

$$\mathbf{Y} = e^{-tM}\mathbf{X} + \sqrt{I_d - e^{-2tM}\mathbf{Z}}, \\ \mathbf{Y}' = e^{-tM}\mathbf{X} + \sqrt{I_d - e^{-2tM}\mathbf{Z}},$$

such that:

$$\begin{pmatrix} \mathbf{Y} \\ \mathbf{Y}' \end{pmatrix} = \underbrace{\begin{pmatrix} e^{-tM} & \sqrt{I_d - e^{-2tM}} & O_{\mathbb{R}^{d \times d}} \\ e^{-tM} & O_{\mathbb{R}^{d \times d}} & \sqrt{I_d - e^{-2tM}} \end{pmatrix}}_{:=B} \begin{pmatrix} \mathbf{X} \\ \mathbf{Z} \\ \mathbf{Z}' \end{pmatrix}.$$

Therefore, using that $(\mathbf{X}, \mathbf{Z}, \mathbf{Z}')^{\top} \sim \mathcal{N}(0, I_{3d})$ we obtain that $(\mathbf{Y}, \mathbf{Y}')$ is a Gaussian vector of mean 0 and covariance $\Gamma = BB^{\top} = \begin{pmatrix} I_d & e^{-2tM} \\ e^{-2tM} & I_d \end{pmatrix}$. The matrix $H(P_{t,\Lambda}(u))$ reads:

$$H(P_{t,M}(u)) = e^{-tM}\mathbb{E}_{\mathbf{Y},\mathbf{Y}'}[\nabla u(\mathbf{Y})\nabla u(\mathbf{Y}')^{\top}]e^{-tM}.$$

Finally, the expression of the conditional probability distribution $\mathbf{Y}'|(\mathbf{Y} = y)$ is a direct consequence of the conditional law formula for Gaussian vectors.

**Appendix G. Calculus details toy example.** We recall the toy model from Section 2 :

$$u: \begin{array}{ccc} \mathbb{R}^2 & \to & \mathbb{R} \\ (x_1, x_2) & \mapsto & \sin(x_1) + \frac{1}{6}\sin(10x_2). \end{array}$$

The gradient of $u$ reads:

$$\nabla u(x_1, x_2) = \begin{pmatrix} \cos(x_1) \\ \frac{10}{6}\cos(10x_2) \end{pmatrix}.$$

For $\mathbf{X} = (\mathbf{x}_1, \mathbf{x}_2) \sim \mathcal{N}(0, I_2)$, $U_{x_1} = (1,0)^{\top}$ and $U_{x_2} = (0,1)^{\top}$ we obtain the following approximation errors:

$$\mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_{x_1}^{\top}\mathbf{X}])^2] = \frac{1}{36}\mathbb{E}[\sin(10\mathbf{x}_2)^2] = \frac{1}{72}(1 - e^{-200})$$

$$\mathbb{E}[(u(\mathbf{X}) - \mathbb{E}[u(\mathbf{X})|U_{x_2}^{\top}\mathbf{X}])^2] = \mathbb{E}[\sin(\mathbf{x}_1)^2] = \frac{1}{2}(1 - e^{-2}),$$
and the Poincaré upper-bounds reads:

$$ \mathbb{E}[\|\nabla u(\mathbf{X})\|_2^2] = \mathbb{E}[\cos(\mathbf{x}_1)^2] + \frac{100}{36}\mathbb{E}[\cos(10\mathbf{x}_2)^2] = \frac{1}{2}(1+e^{-2}) + \frac{100}{72}(1+e^{-200}) $$

$$ \mathbb{E}[\|(I_2 - U_{x_1}U_{x_1}^\top)\nabla u(\mathbf{X})\|_2^2] = \frac{100}{36}\mathbb{E}[\cos(10\mathbf{x}_2)^2] = \frac{100}{72}(1+e^{-200}) $$

$$ \mathbb{E}[\|(I_2 - U_{x_2}U_{x_2}^\top)\nabla u(\mathbf{X})\|_2^2] = \mathbb{E}[\cos(\mathbf{x}_1)^2] = \frac{1}{2}(1+e^{-2}). $$

Finally we can express the mollifying operator $P_{t,I_d}(u)$ as follows:

$$ \begin{aligned} P_{t,I_d}(u)(x_1, x_2) &= \mathbb{E}[\sin(e^{-t}x_1 + \sqrt{1-e^{-2t}}z)] + \frac{1}{6}\mathbb{E}[\sin(10(e^{-t}x_2 + \sqrt{1-e^{-2t}}z))] \\ &= \sin(e^{-t}x_1)\mathbb{E}[\cos(\sqrt{1-e^{-2t}}z)] + \frac{1}{6}\sin(10e^{-t}x_2)\mathbb{E}[\cos(10\sqrt{1-e^{-2t}}z)] \\ &= \sin(e^{-t}x_1)e^{-\frac{1-e^{-2t}}{2}} + \frac{1}{6}\sin(10e^{-t}x_2)e^{-50(1-e^{-2t})}, \end{aligned} $$

for $z \sim \mathcal{N}(0, 1)$.

$$ \nabla P_{t,I_d}(u)(x_1, x_2) = \begin{pmatrix} \cos(e^{-t}x_1)e^{-\frac{1-e^{-2t}}{2}-t} \\ \frac{10}{6}\cos(10e^{-t}x_2)e^{-50(1-e^{-2t})-t} \end{pmatrix} $$

$$ \begin{aligned} \mathbb{E}[\|(I_d - U_{\{1\}}U_{\{1\}}^\top)\nabla P_{t,I_d}(u)(\mathbf{X})\|_2^2] &= \frac{100}{36}e^{-100(1-e^{-2t})-2t}\mathbb{E}[\cos(10e^{-t}\mathbf{x}_2)^2] \\ &= \frac{100}{72}e^{-100(1-e^{-2t})-2t}(1+e^{-200e^{-2t}}) \end{aligned} $$

$$ \begin{aligned} \mathbb{E}[\|(I_d - U_{\{2\}}U_{\{2\}}^\top)\nabla P_{t,I_d}(u)(\mathbf{X})\|_2^2] &= e^{-1+e^{-2t}-2t}\mathbb{E}[\cos(e^{-t}\mathbf{x}_1)^2] \\ &= e^{-1+e^{-2t}-2t}\frac{1+e^{-2e^{-2t}}}{2} \end{aligned} $$