Metamodel construction for sensitivity analysis

TL;DR

This paper introduces a new method for constructing a metamodel to perform sensitivity analysis on complex models using Gaussian regression, combining sparse estimation and functional ANOVA decomposition.

Contribution

It develops a novel approach that estimates both the metamodel and sensitivity indices within a Gaussian regression framework, incorporating penalized least-squares for variable selection.

Findings

The method accurately estimates sensitivity indices in simulations.

The approach provides theoretical guarantees via an oracle inequality.

Simulation results demonstrate effective variable subset selection.

Abstract

We propose to estimate a metamodel and the sensitivity indices of a complex model m in the Gaussian regression framework. Our approach combines methods for sensitivity analysis of complex models and statistical tools for sparse non-parametric estimation in multivariate Gaussian regression model. It rests on the construction of a metamodel for aproximating the Hoeffding-Sobol decomposition of m. This metamodel belongs to a reproducing kernel Hilbert space constructed as a direct sum of Hilbert spaces leading to a functional ANOVA decomposition. The estimation of the metamodel is carried out via a penalized least-squares minimization allowing to select the subsets of variables that contribute to predict the output. It allows to estimate the sensitivity indices of m. We establish an oracle-type inequality for the risk of the estimator, describe the procedure for estimating the metamodel and the sensitivity indices, and assess the performances of the procedure via a simulation study.

Tables

Table 1. Table 1. Estimated coefficient of determination R 2 superscript 𝑅 2 R^{2} for different values of n 𝑛 n and σ 𝜎 \sigma , with the Matérn kernel.

	$\sigma =0$			$\sigma =0.2$
	$n=50$	$n=100$	$n=200$	$n=50$	$n=100$	$n=200$
Proc. GS	0.814	0.920	0.959	0.737	0.835	0.889
Proc. rdg	0.874	0.976	0.989	0.763	0.854	0.892

Table 2. Table 2. Estimated empirical risk ER ER \mathrm{ER} for different values of n 𝑛 n and σ 𝜎 \sigma with the Matérn kernel.

	$\sigma =0$			$\sigma =0.2$
	$n=50$	$n=100$	$n=200$	$n=50$	$n=100$	$n=200$
Proc. GS	0.033	0.0137	0.0139	0.051	0.028	0.020
Proc. rdg	0.011	0.0009	0.0007	0.042	0.022	0.013

Table 3. Table 3. Estimated global error GE × 100 GE 100 \mathrm{GE}\times 100 for different values of n 𝑛 n and σ 𝜎 \sigma with the Matérn kernel.

	$\sigma =0$			$\sigma =0.2$
	$n=50$	$n=100$	$n=200$	$n=50$	$n=100$	$n=200$
Proc. GS	1.91	0.79	0.45	2.41	1.16	0.54
Proc. rdg	0.80	0.10	0.03	1.50	0.47	0.15

Table 4. Table 4. The first line of the table gives the true values of the Sobol indices × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} , as well as their sum in the last columns. The following lines give the mean of the estimators as well as their standard-error (in parenthesis), calculated over 100 simulations, for n = 50 𝑛 50 n=50 , and σ = 0 𝜎 0 \sigma=0 with the Matérn kernel.

	$v=\left\{1\right\}$	$v=\left\{2\right\}$	$v=\left\{3\right\}$	$v=\{1,2\}$	$v=\{1,3\}$	$v=\{2,3\}$	$v=\{1,2,3\}$	sum
S.I.	43.3	24.3	19.2	5.63	4.45	2.50	0.579	99.98
Proc. GS	50.1 (6.2)	26.5 (5.4)	20.9 (4.9)	0.69 (0.8)	0.63 (1.1)	0.51 (0.9)	0.02 (0.07)	99.29
Proc. rdg	45.4 (4.3)	25.3 (3.4)	20.3 (3.5)	3.08 (2.5)	2.18 (2.1)	1.44 (1.8)	0.09 (0.5)	98.66

Table 5. Table 5. The first line of the table gives the true values of the Sobol indices × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} , as well as their sum in the last columns. The following lines give the mean of the estimators as well as their standard-error (in parenthesis), calculated over 100 simulations, for n = 100 𝑛 100 n=100 , and σ = 0.2 𝜎 0.2 \sigma=0.2 with the Matérn kernel.

	$v=\left\{1\right\}$	$v=\left\{2\right\}$	$v=\left\{3\right\}$	$v=\{1,2\}$	$v=\{1,3\}$	$v=\{2,3\}$	$v=\{1,2,3\}$	sum
S.I.	43.3	24.3	19.2	5.63	4.45	2.50	0.579	99.98
Proc. GS	47.5 (5.4)	26.2 (4.9)	19.8 (3.7)	2.35 (1.5)	1.45 (1.2)	0.84 (0.8)	0.03 (0.1)	98.95
Proc. rdg	43.5 (3.9)	25.0 (3.6)	19.7 (2.8)	4.85 (1.7)	3.39 (1.5)	2.02 (1.3)	0.05 (0.3)	99.02

Table 6. Table 6. The first two columns give respectively p ​ S ​ e ​ l S v ≤ ρ p 𝑆 𝑒 subscript 𝑙 subscript 𝑆 𝑣 𝜌 {\mathrm{p}Sel}_{S_{v}\leq\rho} and p ​ S ​ e ​ l S v > ρ p 𝑆 𝑒 subscript 𝑙 subscript 𝑆 𝑣 𝜌 {\mathrm{p}Sel}_{S_{v}>\rho} . The last columns give the values of p ​ S ​ e ​ l v p 𝑆 𝑒 subscript 𝑙 𝑣 {\mathrm{p}Sel}_{v} for each group v 𝑣 v such that S v > ρ subscript 𝑆 𝑣 𝜌 S_{v}>\rho . Results for n = 50 𝑛 50 n=50 and σ = 0 𝜎 0 \sigma=0 with the Matérn kernel.

		$\mathrm{SI}\ge \rho$	$v=\left\{1\right\}$	$v=\left\{2\right\}$	$v=\left\{3\right\}$	$v=\{1,2\}$	$v=\{1,3\}$	$v=\{2,3\}$	$v=\{1,2,3\}$
Proc. GS	17.9	68	100	100	100	72	66	67	9
Proc. rdg	6.3	51	100	100	100	72	62	47	4

Table 7. Table 7. blablaThe first two columns give respectively p ​ S ​ e ​ l S v ≤ ρ p 𝑆 𝑒 subscript 𝑙 subscript 𝑆 𝑣 𝜌 {\mathrm{p}Sel}_{S_{v}\leq\rho} and p ​ S ​ e ​ l S v > ρ p 𝑆 𝑒 subscript 𝑙 subscript 𝑆 𝑣 𝜌 {\mathrm{p}Sel}_{S_{v}>\rho} . The last columns give the values of p ​ S ​ e ​ l v p 𝑆 𝑒 subscript 𝑙 𝑣 {\mathrm{p}Sel}_{v} for each group v 𝑣 v such that S v > ρ subscript 𝑆 𝑣 𝜌 S_{v}>\rho . Results for n = 100 𝑛 100 n=100 and σ = 0.2 𝜎 0.2 \sigma=0.2 with the Matérn kernel.

		$\mathrm{SI}\ge \rho$	$v=\left\{1\right\}$	$v=\left\{2\right\}$	$v=\left\{3\right\}$	$v=\{1,2\}$	$v=\{1,3\}$	$v=\{2,3\}$	$v=\{1,2,3\}$
Proc. GS	38	85	100	100	100	100	99	98	18
Proc. rdg	6	58	100	100	100	98	92	77	3

Table 8. Table 8. Estimated empirical risk ER ER \mathrm{ER} : Performances of the procedures according to the kernel choice.

	$n=50$, $\sigma =0$				$n=100$, $\sigma =0.2$
	Matérn	Gaussian	Brownian	mixed	Matérn	Gaussian	Brownian	mixed
Proc. GS	0.033	0.054	0.027	0.028	0.033	0.054	0.027	0.028
Proc. rdg	0.011	0.024	0.025	0.011	0.011	0.024	0.025	0.023

Table 9. Table 9. Estimated global error GE × 100 GE 100 \mathrm{GE}\times 100 : Performances of the procedures according to the kernel choice.

	$n=50$, $\sigma =0$				$n=100$, $\sigma =0.2$
	Matérn	Gaussian	Brownian	mixed	Matérn	Gaussian	Brownian	mixed
Proc. GS	1.91	2.49	2.19	1.88	1.16	1.90	1.69	1.16
Proc. rdg	0.80	1.39	1.40	0.85	0.48	0.83	0.73	0.51

Equations

$$Y=m(\mathbf{X})+\sigma \epsilon ,$$

$$m(\mathbf{x})=m_0+\sum _{v\in \mathcal{P}}{m}_v({\mathbf{x}}_v),$$

$$S_v=\frac{\mathrm{Var}\left(m_v({\mathbf{x}}_v)\right)}{\mathrm{Var}\left(m(\mathbf{x})\right)}.$$

$$f:\mathcal{X}\to \mathbb{R},f(\mathbf{x})=f_0+\sum _{v\in \mathcal{P}}{f}_v({\mathbf{x}}_v),{\mathbb{E}}_{P_{\mathbf{X}}}{f}_v({\mathbf{x}}_v)={\mathbb{E}}_{P_{\mathbf{X}}}{f}_v({\mathbf{x}}_v)f_{v^{\prime }}({\mathbf{x}}_{v^{\prime }})=0\forall v,v^{\prime }\in \mathcal{P}$$

$$m(\mathbf{x})=m_0+\sum _{v\in \mathcal{P}}{m}_v({\mathbf{x}}_v),$$

$$E_{\mathbf{X}}\left(m_v({\mathbf{X}}_v)\right)=E_{\mathbf{X}}\left(m_v({\mathbf{X}}_v)m_{v^{\prime }}({\mathbf{X}}_{v^{\prime }})\right)=0.$$

$${\mathcal{H}}_{0a}$$

$${\mathcal{H}}_{1a}$$

$$k_{0a}(x_a,x_a^{\prime })=k_a(x_a,x_a^{\prime })-\frac{{\mathbb{E}}_{U\sim P_a}(k_a(x_a,U)){\mathbb{E}}_{U\sim P_a}(k_a(x_a^{\prime },U))}{{\mathbb{E}}_{(U,V)\sim P_a\times P_a}{k}_a(U,V)}.$$

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\prod _{a=1}^d\left(1+k_{0a}({\mathbf{x}}_a,{\mathbf{x}}_a^{\prime })\right)=1+\sum _{v\in \mathcal{P}}{k}_v({\mathbf{x}}_v,{\mathbf{x}}_v^{\prime }),\text{mbox}where{k}_v({\mathbf{x}}_v,{\mathbf{x}}_v^{\prime })=\prod _{a\in v}{k}_{0a}(x_a,x_a^{\prime }),$$

$$\mathcal{H}={\otimes }_{a=1}^d\left(1\stackrel{\perp }{\oplus }{\mathcal{H}}_{0a}\right)=1+\sum _{v\in \mathcal{P}}{\mathcal{H}}_v,$$

$$f(\mathbf{x})=⟨f,k(\mathbf{x},\cdot ){⟩}_{\mathcal{H}}=f_0+\sum _{v\in \mathcal{P}}{f}_v(\mathbf{x}),$$

$$\Vert m-f{\Vert }_{{\mathbb{L}}^2(P_{\mathbf{X}})}^2=E_{\mathbf{X}}{\left(m(\mathbf{X})-f(\mathbf{X})\right)}^2$$

$$\Vert f{\Vert }_n=\sqrt{\frac{1}{n}\sum _{i=1}^n{f}_v^2({\mathbf{X}}_{v,i})}.$$

$$\mathcal{L}(f)=\frac{1}{n}\sum _{i=1}^n{\left(Y_i-f_0-\sum _{v\in \mathcal{P}}{f}_v({\mathbf{X}}_{v,i})\right)}^2+\sum _{v\in \mathcal{P}}{\mu }_v\Vert f_v{\Vert }_{{\mathcal{H}}_v}+\sum _{v\in \mathcal{P}}{\gamma }_v\Vert f_v{\Vert }_n.$$

$$\mathcal{F}=\left\{f\text{mbox}suchthatf=f_0+\sum _{v\in \mathcal{P}}{f}_v,f_v\in {\mathcal{H}}_v,\Vert f_v{\Vert }_{{\mathcal{H}}_v}\le 1\right\}.$$

$$\hat{f}=\mathrm{argmin}\left\{\mathcal{L}(f),f\in \mathcal{F}\right\}.$$

$$\mathrm{Var}\left(m(\mathbf{x})\right)=\sum _{v\in \mathcal{P}}\mathrm{Var}\left(m_v({\mathbf{x}}_v)\right).$$

$$S_v=\frac{\mathrm{Var}\left(m_v({\mathbf{x}}_v)\right)}{\mathrm{Var}\left(m(\mathbf{x})\right)},$$

$$G_a=\frac{{\sum }_{v\supseteq \left\{a\right\}}\mathrm{Var}\left(m_v({\mathbf{x}}_v)\right)}{\mathrm{Var}\left(m(\mathbf{x})\right)}.$$

$$\widehat{\mathrm{Var}}\left(m(\mathbf{x})\right)=\sum _{v\in \mathcal{P}}\widehat{\mathrm{Var}}\left(m_v({\mathbf{x}}_v)\right),\text{mbox}where\widehat{\mathrm{Var}}\left(m_v({\mathbf{x}}_v)\right)={\mathbb{E}}_{\mathbf{X}}\left({\hat{f}}_v^2({\mathbf{X}}_v)\right)=\Vert {\hat{f}}_v{\Vert }_{{\mathbb{L}}^2(P_{\mathbf{X}})}^2.$$

$${\widehat{\mathrm{Var}}}^{\mathrm{emp}}\left(m_v({\mathbf{x}}_v)\right)=\frac{1}{n-1}\sum _{i=1}^n{\left({\hat{f}}_v({\mathbf{X}}_{v,i})-{\hat{f}}_{v,\cdot }\right)}^2.$$

$$S_f=\left\{v\in \mathcal{P},f_v\ne 0\right\}.$$

$$\forall f_v\in {\mathcal{H}}_v\Vert f_v{\Vert }_{\infty }=\sup \left\{|f_v(\mathbf{X})|,\mathbf{X}\in \mathcal{X}\right\}\le R^{\prime }.$$

$$Q_{n,v}(t)=\sqrt{\frac{5}{n}\sum _{k\ge 1}\min (t^2,{\omega }_{v,k})},$$

$${\nu }_{n,v}=\inf \left\{t\text{mbox}suchthat{Q}_{n,v}(t)\le \Delta t^2\right\}.$$

$${\lambda }_{n,v}$$

$$\text{mbox}forallv\in \mathcal{P},{\lambda }_{n,v}\le \min \left\{\frac{{\gamma }_v}{C_1},\sqrt{\frac{{\mu }_v}{C_1}}\right\},n{\lambda }_{n,v}^2\ge -C_2\log {\lambda }_{n,v},$$

$$\text{mbox}forallf\in \mathcal{F},\max \left(\sum _{v\in S_f}{\gamma }_v^2,\sum _{v\in S_f}{\mu }_v\right)\le 1,$$

$$\Vert m-\hat{f}{\Vert }_n^2\le C_3\underset{f\in \mathcal{F}}{\inf }\left\{\Vert m-f{\Vert }_n^2+{\sigma }^2\sum _{v\in S_f}\left({\mu }_v+{\gamma }_v^2\right)\right\}.$$

$${\mathcal{F}}_r=\left\{f\text{mbox}suchthatf=f_0+\sum _{v\in \mathcal{P}}{f}_v,f_v\in {\mathcal{H}}_v,\Vert f_v{\Vert }_{{\mathcal{H}}_v}\le r_v\right\}.$$

Full text

{resume}

Nous considérons l’estimation d’un méta-modèle d’un modèle complexe $m$ à partir des observations d’un $n$-échantillon dans un modèle de régression Gaussien. Nous en déduisons une estimation des indices de sensibilité de $m$. Notre approche combine les méthodes d’analyse de sensibilité de modèles complexes et les outils statistiques de l’estimation non-paramétrique en régression multivariée. Elle repose sur la construction d’un méta-modèle qui approche la décomposition de Hoeffding-Sobol de $m$. Ce méta-modèle appartient à un espace de Hilbert à noyau reproduisant qui est lui-même la somme directe d’espaces de Hilbert, permettant ainsi une décomposition de type ANOVA. On en déduit des estimateurs des indices de sensibilité de $m$. Nous établissons une inégalité de type oracle pour le risque de l’estimateur, nous décrivons la procédure pour estimer le méta-modèle et les indices de sensibilité, et évaluons les performances de notre méthode à l’aide d’une étude de simulations.

\runningauthorsS. Huet and M.L. Taupin

1. Introduction

We consider a Gaussian regression model

[TABLE]

where $\mathbf{X}$ is a $d$ random vector with a known distribution $P_{\mathbf{X}}=P_{1}\times\ldots\times P_{d}$ on $\mathcal{X}$ a compact subset of $\mathbb{R}^{d}$, and $\varepsilon$ is independent of $\mathbf{X}$, and distributed as a $\mathcal{N}(0,1)$ variable. The variance $\sigma^{2}$ is unknown and the number of variables $d$ may be large. The function $m$ is a complex and unknown function from $\mathbb{R}^{d}$ to $\mathbb{R}$. It may present strong non-linearities and high order interaction effects between its coordinates. On the basis of a $n$-sample $(Y_{i},\mathbf{X}_{i}),i=1,\ldots,n$, we aim to construct metamodels and perform sensitivity analysis in order to determine the infuence of each variable or group of variables on the output.

Our approach combines methods for sensitivity analysis of a complex model and statistical tools for sparse non-parametric estimation in multivariate Gaussian regression model. It rests on the construction of a metamodel for aproximating the Hoeffding-Sobol decomposition of the function $m$. This metamodel belongs to a reproducing kernel Hilbert space constructed as a direct sum of Hilbert spaces leading to a functional ANOVA decomposition involving variables and interactions between them. The estimation of the metamodel is carried out via a penalized least-squares minimization allowing to select the subsets of variables $\mathbf{X}$ that contribute to predict the output $Y$. Finally, the estimated metamodel allows to estimate the sensitivity indices of $m$.

A lot of work has been done around meta-modelling and sensitivity indices estimation.

For a complete account on global sensitivity analysis, see for example the book by Saltelli et al. [35]. Let us briefly present the context of usual global sensitivity analysis. Suppose that we are able to calculate the ouputs $z$ of a model $m$ for $n$ realizations of the input vector $\mathbf{X}$, such that $z_{i}=m(\mathbf{X}_{i})$ for $i=1,\ldots,n$. Starting from the values $(z_{i},\mathbf{X}_{i}),i=1,\ldots n$, the objectives of meta-modelling and global sensitivity analysis are to approximate the function $m$ by what is called a metamode or are to quantify the influence of some subsets of the variables $\mathbf{X}$ on the output $z$. This metamodel helps to understand the behavior of the model, or allows to speed up future calculation using it in place of the original model $m$.

In particular when the inputs variables $\mathbf{X}$ are independent, if $m$ is square integrable, one may consider the classical Hoeffding-Sobol decomposition [36, 41] that leads to write $m$ according to its ANOVA functional expansion:

[TABLE]

where $\mathcal{P}$ denotes the set of parts of $\left\{1,\ldots,d\right\}$ with dimension 1 to $d$ and where for all $\mathbf{x}\in\mathbb{R}^{d}$, $\mathbf{x}_{v}$ denotes the vector with components $x_{j}$ for $j\in v$. The functions $m_{v}$ are centered and orthogonal in $\mathbb{L}^{2}(P_{\mathbf{X}})$ leading to the following decomposition of the variance of $m$: $\mathop{\rm Var}\nolimits\left(m(\mathbf{x})\right)=\sum_{v}\mathop{\rm Var}\nolimits\left(m_{v}(\mathbf{x}_{v})\right)$. The Sobol sensitivity indices, introduced by Sobol [36], are defined for any group of variables $\mathbf{x}_{v}$, $v\in\mathcal{P}$ by

[TABLE]

They quantify the contribution of a subset of variables $\mathbf{x}$ to the output $m(\mathbf{x})$. Several approaches are available for estimating these sensitivity indices, see for example Iooss and Lemaître [17] for a recent review. Among all of them, let us consider the one based on metamodel construction that allows to directly obtain the sensitivity indices. Generally one consider metamodels that correspond to an ANOVA-type decomposition, and that are candidate to approximate the Hoeffding decomposition of $m$. The ANOVA-type decomposition leads to consider functions defined as follows:

[TABLE]

for functions $f_{v}$ that are chosen to belong to some functionnal spaces. The polynomial Chaos construction, see for example Ghanem and Spanos [14], Soize and Ghanem [38], can be used to approximate the Hoeffding decomposition of $m$. This approach was considered by Blatman and Sudret [5] who propose a method for truncating (such that to keep polynomials with degree less than some integer) the polynomial Chaos expansion and then an algorithm based on least angle regression for selecting the terms in the expansion. For the same purpose, Gu and Wu [15] propose an algorithm based on the hierarchy principle (lower order effects are more likely to be important than higher order effects) and on the heredity principle (interaction can be active only if one or all of its parent effects are also active). This approach joins the one proposed by Bach [2] for variable selection based on hierarchical kernel learning.

Inspired by Touzani [39], Durrande et al. [12] propose to approximate $m$ by functions belonging to a reproducing kernel Hilbert space (RKHS). The RKHS is constructed as a direct sum of Hilbert spaces leading to a functional ANOVA decomposition (see Equation (3)), such that the projection of $m$ onto the RKHS is an approximation of the Hoeffding decomposition of $m$.

Following Lin and Zhang [25], Touzani and Busby [40] propose an algorithm to calculate the penalized least-square estimator of $m$ on the RKHS space, where the least-square criteria is penalized by the sum of the norms of $m$ on each Hilbert subspace. This group-lasso type procedure allows both to select and calculate the non-zero terms in the functional ANOVA decomposition.

Our objective is to propose an estimator of a metamodel which will approximate the Hoeffding decomposition of $m$ considering a Gaussian regression model defined at Equation (1) and to deduce from this estimated metamodel, estimators for the sensitivity indices of $m$. Contrary to the usual setting of sensitivity analysis where $m(\mathbf{X}_{i})$ is available, only the observations $Y$ are available, which leads us to consider the nonparametric multivariate regression setting.

Let us briefly describe the methods related to this regression setting and review their theoretical properties, starting with papers assuming an univariate additive decomposition for the function $m$ in the context of high-dimensional sparse additive models. Precisely, denote by $\mathcal{F}^{\textrm{1-add}}$, the set of functions $f$ defined on $\mathcal{X}$ such that $f(\mathbf{x})=f_{0}+\sum_{a=1}^{d}f_{a}(x_{a})$ where $f_{0}$ is a constant, and where for $a=1,\ldots,d$, the functions $f_{a}$ are centered and square-integrable with respect to $P_{a}$. For each function $f$, the set $S_{f}$ of indices $a\in\left\{1,\ldots d\right\}$ such that $f_{a}$ is not identically zero is called the actice set of $f$.

Ravikumar et al. [32] propose a group-lasso procedure where each function $f_{a}$ is approximated by its truncated decomposition on a basis of functions. They provide an algorithm and, assuming that the function $m$ belongs to the set $\mathcal{F}^{\textrm{1-add}}$ and that $S_{m}$ is sparse, prove the consistency of the active set and of the risk of the estimator of $m$.

Meier et al. [27] propose to combine a sparsity penalty (group-lasso) and a smoothness penalty (ridge) for estimating $m(\mathbf{x})$. Considering the fixed design framework, they establish some oracle properties of the empirical risk for estimating the projection of $m$ onto the set of univariate additive functions $\mathcal{F}^{\textrm{1-add}}$.

Raskutti et al. [31] consider the case where each univariate function $f_{a}$ belongs to a RKHS and as Meier et al. combine a sparsity and a smoothness penalty. Assuming that the $d$ variables $\mathbf{X}$ are independent, they derive upper bounds for the integrated and the empirical risks, as well as a lower bound for the integrated risk over spaces of sparse additive models whose each component is bounded with respect to the RKHS norm.

Additive sparse modelling is too restrictive in practical settings because it does not take into account interactions between variables that may affect the relationship between $Y$ and $\mathbf{X}$. The generalization of additive smoothing splines to interaction smoothing splines leading to an ANOVA-type decomposition (see Equation (3)) was proposed by several authors (see for example Wahba [42], Friedman [13], Wahba et al. [43]).

To control smoothness and to enforce sparsity in the ANOVA-type decomposition, Gunn and Kandola [16] propose to consider the ANOVA decomposition as a weighted linear sum of kernels and to use a lasso penalty on the weights to select the terms in the decomposition as well as a ridge penalty to ensure smoothness of the kernel expansion. The COSSO proposed by Lin and Zhang [25] is based on smoothness penalty defined as the sum of the RKHS-norms of the functions $f_{v}$. The authors study existence and rate of convergence of the estimator. In a more general framework, where the function $m$ is written as a linear span of a large number of kernels, Koltchinskii and Yuan [21] established oracle inequalities on the excess risk assuming that the function $m$ has a sparse representation (the set of $v\in\mathcal{P}$ such that $f^{*}_{v}$ is non zero is sparse). The authors generalized their results to a penalty function that combines sparsity and smoothness [22], as proposed by Meier et al. [27] and Raskutti et al. [31].

Recently Kandasamy and Wu [19] proposed an estimator called SALSA, based on a ridge penalty, where the ANOVA-type decomposition is restricted to set $v\in\mathcal{P}$ such that $|v|\leq D_{\max}$. The authors propose to choose $D_{\max}$ via a cross-validation procedure.

Our contributions

Using the functionnal ANOVA-type decomposition as proposed by Durrande et al. [12], we propose an estimator of a metamodel which approximates the Hoeffding decomposition of $m$. Following the most recent works in the framework of nonparametric estimation of sparse additive models, we propose a penalized least-square estimator where the penalty function enforces both the sparsity and the smoothness of the terms in the decomposition. We show that our estimator satisfies an oracle inequality with respect to the empirical and integrated risks.

Our procedure allows both to select and estimate the terms in the ANOVA decomposition, and therefore, to select the sensitivity indices that are non-zero and estimate them. In particular it makes possible to estimate Sobol indices of high order, a point known to be difficult in practice.

Finally, using convex optimization tools, we develop an algorithm (available on request) in $R$ [30], for calculating the estimator. A simulation study shows the good performances of our estimator in practice.

The paper is organised as follows: The RKHS construction based on ANOVA kernels and the procedure for estimating a metamodel are presented in Section 2. The estimators of the Sobol indices are given in Section 3. The theoretical properties of the metamodel estimator are stated in Theorem 4.1 and Corollaries 4.1 and 4.2 whose proofs are postponed in Sections 7 and 8. Section 5 is devoted to the calculation of the estimator and Section 6 to the simulation study.

2. Meta-modelling

We start from the Hoeffding decomposition (see Sobol [37] and Van der Vaart [41], p. 157) of the function $m$ that consists in writting $m$ as in Equation (2)

[TABLE]

where $\mathcal{P}$ denotes the set of parts of $\left\{1,\ldots,d\right\}$ with dimension 1 to $d$ and where for all $\mathbf{x}\in\mathbb{R}^{d}$, $\mathbf{x}_{v}$ denotes the vector with components $x_{j}$ for $j\in v$. For all $v,v^{\prime}$ in $\mathcal{P}$,

[TABLE]

We propose to consider a functionnal space based on the tensorial product of Reproducing Kernel Hilbert spaces (RKHS), and to approximate the unknown function $m$ by its projection denoted $f^{*}$ on such such RKHS space. One of the key point is to construct the space $\mathcal{H}$ such that the terms of the decomposition of a function $f$ in $\mathcal{H}$ correspond to its Hoeffding-Sobol decomposition.

2.1. RKHS construction

Let us describe the construction of spaces $\mathcal{H}$, based on ANOVA kernels, construction which was given by Durrande et al. [12].

Let $\mathcal{X}=\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{d}$ be a compact subset of $\mathbb{R}^{d}$. For each coordinate $a\in\{1,\cdots,d\}$, we choose a RKHS $\mathcal{H}_{a}$ and its associated kernel $k_{a}$ defined on the set $\mathcal{X}_{a}\subset\mathbb{R}$ such that the two following properties are satisfied

1. (1)

$k_{a}:\mathcal{X}_{a}\times\mathcal{X}_{a}\rightarrow\mathbb{R}$ is $P_{a}\times P_{a}$ measurable 2. (2)

$\mathbb{E}_{P_{a}}\sqrt{k_{a}(X_{a},X_{a})}<\infty$

The RKHS $\mathcal{H}_{a}$ may be decomposed as $\mathcal{H}_{a}=\mathcal{H}_{0a}\stackrel{{\scriptstyle\perp}}{{\oplus}}\mathcal{H}_{1a}$, where

[TABLE]

the kernel $k_{0a}$ associated to the RKHS $\mathcal{H}_{0a}$ being defined as follows (see Berlinet et Thomas-Agnan [4]):

[TABLE]

The ANOVA kernel is defined as

[TABLE]

and the corresponding RKHS

[TABLE]

where the RKHS $\mathcal{H}_{v}$ is associated with kernel $k_{v}$. According to this construction, any function $f\in\mathcal{H}$ satisfies

[TABLE]

where $f_{v}(\mathbf{x})=\langle f,k_{v}(\mathbf{x},\cdot)\rangle_{\mathcal{H}}$ depends on $\mathbf{x}_{v}$ only. For all $v\in\mathcal{P}$, $f_{v}(\mathbf{x}_{v})$ is centered and for all $v^{\prime}\neq v$, $f_{v}(\mathbf{x}_{v})$ and $f_{v^{\prime}}(\mathbf{x}_{v^{\prime}})$ are uncorrelated. We get thus the Hoeffding decomposition of $f$.

2.2. Approximating the Hoeffding decomposition of $m$

Let $f^{*}=f_{0}^{*}+\sum_{v\in\mathcal{P}}f_{v}^{*}$ which minimizes

[TABLE]

over functions $f\in\mathcal{H}$. This $f^{*}$ can be viewed as an approximation of $m$ and more specifically his Hoeffding decomposition is an approximation of the Hoeffding decomposition of $m$. Therefore if the Hoeffding decomposition of $m$ is written as in Equation (2), each function $f^{*}_{v}$ approximates the function $m_{v}$.

The idea is propose an estimator of $f^{*}$ as estimator of $m$.

2.3. Selection step

Since $\mathcal{P}$ is the set of parts of $\left\{1,\ldots,d\right\}$, the number of functions$f^{*}_{v}$ is related to the cardinality of $\mathcal{P}=2^{d}-1$ that may be huge. Our construction is thus associated to a selection strategy.

The selection of $f^{*}_{v}$ in $f^{*}$ is based on a ridge-group-sparse type procedure which minimizes the penalized least-squares criteria over functions $f\in\mathcal{H}$. The least-squares criteria is penalized in order to both select few terms in the additive decomposition of $f$ over sets $v\in\mathcal{P}$, and to favour smoothness of the estimated $f_{v}$. The ridge regularization is ensured by controling the norm of $f_{v}$ in the Hilbert space $\mathcal{H}_{v}$ for all $v$, and the group-sparse regularization is strengthened by controling the empirical norm of $f_{v}$, defined as

[TABLE]

For any $f\in\mathcal{H}$ such that $f=f_{0}+\sum_{v\in\mathcal{P}}f_{v}$, and for some tuning parameters $(\mu_{v},\gamma_{v},v\in\mathcal{P})$, let $\mathcal{L}(f)$ be defined as

[TABLE]

Let us define the set of functions $\mathcal{F}$

[TABLE]

Then the estimator $\widehat{f}$ is defined as

[TABLE]

Remark 2.1**.**

The construction of the RKHS spaces described above, allows to consider functionnal spaces that suit well to the smoothness of the function $m$, irrespectively of the distribution of $\mathbf{X}$. Indeed, the kernels $k_{0,a}$ depend on the distribution of $\mathbf{X}$ only for calculating the projection onto the space of constant functions. In comparison, the decomposition based on the truncated polynomial Chaos expansion, used for sensitivity analysis (see Blatman and Sudret [5]), is based on the distribution of $\mathbf{X}$, and only the choice of the truncation handles the smoothness of the approximation.

3. Sensitivity analysis

3.1. Sobol indices

Let us go back to the Hoeffding decomposition Equation (2). The orthogonality between two terms in this decomposition leads to the additive decomposition of the variance of $m(\mathbf{x})$:

[TABLE]

Each of these variance terms are related to Sobol indices [36]. For example, the Sobol indice linked with the interaction between variables $\mathbf{x}_{v}$ is defined as

[TABLE]

or the global Sobol indices for the variable $x_{a}$, $a\in\{1,\cdots,d\}$, is

[TABLE]

Those Sobol indices and global Sobol indices quantify the contribution of a subset of variables $\mathbf{x}$ to the output $m(\mathbf{x})$. As it is said in the introduction direct estimation of these Sobol indices may require lot of calculations. We consider here methods based on metamodels to directly obtain Sensitivity indices.

3.2. Estimation of Sobol indices

Thanks to the orthogonality property of functions in $\mathcal{H}$, the variance of $m(\mathbf{x})$ will be estimated by

[TABLE]

In practice, in order to avoid calculating the variance of $\widehat{f}_{v}(\mathbf{X}_{v})$, one may use an estimator based on the empirical variances of functions $\widehat{f}_{v}$. Precisely, if $\widehat{f}_{v,\cdot}$ is the mean of the $\widehat{f}_{v}(\mathbf{X}_{v,i})$, for $i=1,\ldots,n$, then

[TABLE]

One of the main contribution of this approach is to allow the estimation of Sobol indices of any order, whereas classical methods only deal with small order, generally less than two.

4. Theoretical result: oracle inequality for metamodel

In this section we state the oracle inequality for the estimated metamodel $\widehat{f}$ which approximates the Hoeffding decomposition of the unknown function $m$.

4.1. Notations and Assumptions

For a function $f\in\mathcal{H}$, $f=f_{0}+\sum_{v\in\mathcal{P}}f_{v}$, we denote by $S_{f}$ its support and $|S_{f}|$ its cardinality. More precisely

[TABLE]

We consider RKHS spaces $\mathcal{H}_{v},v\in\mathcal{P}$ satisfying the following assumptions:

• •

for all $f_{v}\in\mathcal{H}_{v}$, $E_{\mathbf{X}}f^{2}_{v}(\mathbf{X})<\infty$,

• •

for all $f_{v}\in\mathcal{H}_{v}$, $f_{v^{\prime}}\in\mathcal{H}_{v^{\prime}}$, $E_{\mathbf{X}}f_{v}(\mathbf{X})=0$ and $E_{\mathbf{X}}f_{v}(\mathbf{X})f_{v^{\prime}}(\mathbf{X})=0$,

• •

there exists $R^{\prime}>0$ such that

[TABLE]

For each $v\in\mathcal{P}$, let $\omega_{v,k}$, for $k\geq 1$ be the eigenvalues of the operator associated to the self reproducing kernel $k_{v}$, arranged in the decreasing order. Let us define the function $Q_{n,v}(t)$, for positive $t$, as follows:

[TABLE]

and for some $\Delta>0$ let $\nu_{n,v}$ be defined as follows

[TABLE]

For each $v\in\mathcal{P}$, $\nu_{n,v}$ refers to the so-called critical univariate rate, the minimax-optimal rate for $\mathbb{L}^{2}(P_{\mathbf{X}})$-estimation of a single univariate function in the hilbert space $\mathcal{H}_{v}$ (e.g. Mendelson [28]).

Our choices of regularization parameters and rates are specified in terms of the quantities:

[TABLE]

Theorem 4.1**.**

Let us consider the regression model defined at Equation (1). Let $(Y_{i},\mathbf{X}_{i})$, $i=1,\ldots,n$ be a $n$-sample with the same law as $(Y,\mathbf{X})$. Let $\widehat{f}$ be defined by (6).

If there exist constants $C_{l},l=1,2,3$, $C_{1}\geq 1$, and $0<\eta<1$ such that the following conditions are satisfied:

[TABLE]

and

[TABLE]

then, with probability greater than $1-\eta$,

[TABLE]

The result can be easily generalized to the minimization of $\mathcal{L}(f)$ over the space $\mathcal{G}$ defined as follows: For some positive constants $r_{v},v\in\mathcal{P}$,

[TABLE]

Indeed, we just have to consider the RKHS $\mathcal{H}^{{}^{\prime}}_{v}$ associated with the kernel $k^{{}^{\prime}}_{v}(\mathbf{x}_{v},\mathbf{x}^{{}^{\prime}}_{v})=r_{v}^{2}k_{v}(\mathbf{x}_{v},\mathbf{x}^{{}^{\prime}}_{v})$ in place of the RKHS $\mathcal{H}_{v}$. Then minimizing

[TABLE]

over the set

[TABLE]

is equivalent to minimizing $\mathcal{L}(f)$ over the set $\mathcal{F}_{r}$.

Let us now comment on the theorem.

• •

The term $\|m-f\|_{n}^{2}$ refers to the usual bias term quantifying the approximation properties of the Hilbert space $\mathcal{H}$ as a distance between the true $m$ and $f$, its approximation into $\mathcal{H}$.

• •

Koltchinskii et Yuan [22] considered what is called the multiple kernel learning problem, where the functions in $\mathcal{H}$ have an additive representation over kernel spaces. They do not assume that the variable $\mathbf{X}$ are independent, nor that the kernel spaces satisfy an orthogonality condition. In return, they assume that some decomposability type properties are satisfied, and they introduce some characteristics related to the degree of “dependence” of the kernel spaces.

• •

In the particular case when the decomposition is limited to the main effects of the variables, then the problem comes back to the classical nonparametric additive model. The theoretical properties of the estimator based on a ridge-group-sparse type procedure have already been established (see for example Meier et al. [27], and Raskutti et al. [31]).

• •

Weights in the penalty terms may be of interest for applications. The theoretical result highlights that the tuning parameters $(\mu_{v},\gamma_{v})$ should depend on the decreasing of the eigenvalues of the kernel defining $\mathcal{H}_{v}$. Besides, we may be interested by introducing weigths that favor small order interaction terms.

• •

Because we aim to approximate the Hoeffding decomposition of $m$, we need to have orthogonality between the spaces $\mathcal{H}_{v},v\in\mathcal{P}$. This condition, required by our objective, is also a key point in the proof of the Theorem, when the problem is to compare the euclidean norm of functions in $\mathcal{H}$ with the norm in $\mathbb{L}^{2}(P_{\mathbf{X}})$. At this step we need to assume that Assumption (15) holds to conclude.

• •

We do not assume that the functions in $\mathcal{F}$ are uniformly globally bounded, that is that the $\sup\{|f(\mathbf{x})|,\mathbf{x}\in\mathcal{X}\}$ is bounded by a constant that does not depend on $f$. Instead we assume that each function within the unit ball of the Hilbert space $\mathcal{H}_{v}$ is uniformly bounded by a constant multiple of its Hilbert norm. In fact, the functions $f$ in the space $\mathcal{F}$, written as $f_{0}+\sum_{v}f_{v}$ satisfy that for each group $v$, $f_{v}$ is uniformly bounded. This assumption is easily satisfied as soon as the kernel $k_{v}$ is bounded on the compact set $\mathcal{X}$. Indeed, $\|f_{v}\|_{\infty}\leq\sup_{\mathbf{X}\in\mathcal{X}}\sqrt{k_{v}(\mathbf{X}_{v},\mathbf{X}_{v})}\|f_{v}\|_{\mathcal{H}_{v}}$. We refer to [31] for a discussion on that subject and a comparison with the work of Koltchinskii et Yuan [22].

• •

Assuming that $n\lambda_{n,v}^{2}\geq-C_{2}\log\lambda_{n,v}$ allows to control the probability of the union of $|\mathcal{P}|$ events. This is a mild condition, satisfied for $\lambda_{n,v}=K_{n}/\sqrt{n}$ for $K_{n}$ of the order $\sqrt{\log n}$.

The following corollary gives an upper bound of the risk with respect to the $\mathbb{L}^{2}(P_{\mathbf{X}})$ norm. It is mainly a consequence of Theorem 4.1 (its proof is given Section 8.2 page 8.2).

Corollary 4.1**.**

Under the assumptions of Theorem 4.1, we have that, with probability greater than $1-\eta$, for some constant $C_{4}$,

[TABLE]

From this corollary we can compare the $\widehat{\mathop{\rm Var}\nolimits}\left(m_{v}(\mathbf{x}_{v})\right)$ (see Equation (7)) with the variance of $m_{v}(\mathbf{x}_{v})$. Thanks to the following inequality

[TABLE]

and to Corollary 4.1, we get the following result

[TABLE]

4.2. Rate of convergence

Corollary 4.2**.**

Under the same condition as Theorem 4.1, if $\gamma_{v}=c\lambda_{n,v}$ and $\mu_{v}=c\lambda_{n,v}^{2}$ for $c\geq C1$, then

[TABLE]

The result is non asymptotic in the sense that it is shown for any $(n,d)$. Nevertheless, the upper bound is relevant when the infimum is reached for functions $f$ whose decomposition in $\mathcal{H}$ is sparse, and when $d$ is small face to $n$. In fact, the coefficient $d$ occuring in the rate $d|S_{f}|/n$ comes from the logarithm of the cardinality of $\mathcal{P}$ equal to $\log(2^{d}-1)$. When $d$ is large, it may be judicious to limit the decomposition of functions in $\mathcal{H}$, to interactions of limited order, so that the number of terms in the decomposition stays of the order $\log(d)$.

Let us discuss the rate of convergence given by $\sum_{v\in S_{f}}\nu^{2}_{n,v}$. For the sake of simplicity let us consider the case where the variables $X_{1},\ldots,X_{d}$ have the same distribution $P_{1}$ on $\mathcal{X}_{1}\subset\mathbb{R}$, and where the unidimensionnal kernels $k_{0a}$ are all identical, such that $k_{v}(\mathbf{x}_{v},\mathbf{x}^{\prime}_{v})=\prod_{a\in v}k_{0}(x_{a},x^{\prime}_{a})$. The kernel $k_{0}$ admits an eigen expansion given by

[TABLE]

where the eigenvalues $\omega_{0,\ell}$ are non negative and ranged in the decreasing order, and where the $\zeta_{\ell}$ are the associated eigen functions, orthonormal with respect to $\mathbb{L}^{2}(P_{1})$. Therefore the kernel $k_{v}$ admits the following expansion

[TABLE]

Consider now the case where the eigenvalues $\omega_{0,\ell}$ are decreasing at a rate $\ell^{-2\alpha}$ for some $\alpha>1/2$. It can be shown, see Section 8.3, that the rate $\nu_{n,v}$ defined at Equation (12) is bounded above by a term of order $n^{-\alpha/(2\alpha+1)}(\log n)^{\gamma}$, where $\gamma\geq(|v|-1)\alpha/(2\alpha-1)$. Note that in this particular case, the rate of convergence depends on $|v|$ through the logarithmic term, and that up to this logarihmic term the rate of convergence has the same order than the usual nonparametric rate for unidimensionnal functions. It follows that the RKHS space $\mathcal{H}$ should be chosen such that the unknown function $m$ is well approximated by sparse functions in $\mathcal{H}$ with low order of interactions.

5. Calculation of the estimator

The functional minimization problem described at Equation (6) is equivalent to a parametric minimization problem. Indeed, we know that if $\mathcal{H}$ is a RKHS associated with a kernel $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, then for all $(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})\in\mathcal{X}^{n}$, and for all $(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{R}^{n}$, the function $f(\cdot)=\sum_{i=1}^{n}\alpha_{i}k(\mathbf{x}_{i},\cdot)$ is in $\mathcal{H}$ and $\|f\|^{2}_{\mathcal{H}}=\sum_{i,i^{\prime}=1}^{n}\alpha_{i}\alpha_{i^{\prime}}k(\mathbf{x}_{i},\mathbf{x}^{\prime}_{i^{\prime}})$. In particular, it can be shown that the solution to our minimization problem is written as $f=f_{0}+\sum_{v\in\mathcal{P}}f_{v}$ where, according to the representer Theorem (see Kimeldorf and Whahba [20]), $f_{v}(\cdot)=\sum_{i=1}^{n}\theta_{vi}k_{v}(\mathbf{X}_{vi},\cdot)$ for some parameter $\boldsymbol{\theta}\in\mathbb{R}^{n|\mathcal{P}|}$ with components $(\theta_{v,i},i=1,\ldots,n,v=1,\ldots,|\mathcal{P}|)$.

Let $\|\cdot\|$ denotes the usual euclidean norm in $\mathbb{R}^{n}$. For each $v\in\mathcal{P}$, let $K_{v}$ be the $n\times n$ matrix with components $(K_{v})_{i,i^{\prime}}=k_{v}(\mathbf{X}_{vi},\mathbf{X}_{vi^{\prime}})$ that satisfies $t(K^{1/2})K^{1/2}=K$. Let $\widehat{f}_{0}$ and $\widehat{\boldsymbol{\theta}}$ be the minimizer of the following penalized least-squares criteria:

[TABLE]

Then the estimator $\widehat{f}$ defined at Equation (6) satisfies

[TABLE]

Because $C(f_{0},\boldsymbol{\theta})$ is a convex and separable criteria, we propose to calculate $\widehat{\boldsymbol{\theta}}$ using a block coordinate descent algorithm described in the following section.

Note that the estimator $\widehat{f}$ defined at Equation (6) should satisfy $\|f_{v}\|_{\mathcal{H}_{v}}=\|K^{1/2}_{v}\boldsymbol{\theta}_{v}\|\leq 1$, or generally $\|f_{v}\|_{\mathcal{H}_{v}}\leq r_{v}$ for some positive $r_{v}$, see (16). Usually we have no idea of the value of this upper-bound in practice, and we propose to remove this contraint in the optimization procedure. Nevertheless, if one wants to consider such an additional constraint, the problem can be solved at the price of some additional complication, considering a Lagrangian method for example, see Section 8.7 for more details.

5.1. Algorithm

We will assume that for all $v\in\mathcal{P}$ the matrices $K_{v}$ are strictly definite positive. If it is not the case, one modifies $K_{v}$ by $K_{v}+\xi I_{n}$ where $\xi$ is a small positive value, in order to ensure positive definiteness.

Using a coordinate descent procedure, we minimize the criteria $C(f_{0},\boldsymbol{\theta})$ along each group $v$ at a time. At each step of the algorithm, the criteria is minimized as a function of the current block’s parameters, while the parameters values for the other blocks are fixed to their current values. The procedure is repeated until convergence, considering for example that the convergence is obtained if the norm of the difference between two consecutive solutions is small enough. See for exemple Boyd et al. [8] for optimization in such context.

For the sake of simplicity, we consider the minimization of the following criteria:

[TABLE]

Taking $\gamma^{\prime}_{v}=\sqrt{n}\gamma_{v}$ and $\mu^{\prime}_{v}=n\mu_{v}$, this is exactly the same criteria as the one defined at Equation (17).

Let us begin with the constant term $f_{0}$. Because the penalty function does not depend on $f_{0}$, minimizing $C^{\prime}(f_{0},\boldsymbol{\theta})$ with respect to $f_{0}$ for fixed values of $\boldsymbol{\theta}$ leads to

[TABLE]

where $Y_{\cdot}$ denotes the mean of $\mathbf{Y}$ and $\left(K_{v}\boldsymbol{\theta}_{v}\right)_{i}$ denotes the $i$-th component of $K_{v}\boldsymbol{\theta}_{v}$. In what follows, we consider a group $v$, and fix the values of the parameters for all the other groups. We describe the algorithm and postpone the proofs in Section 8.7.

Let us first consider the case

where both $\mu^{\prime}_{v}$ and $\gamma^{\prime}_{v}$ are non zero

If $\partial C^{\prime}_{v}$ denotes the subdifferential of $C^{\prime}(f_{0},\boldsymbol{\theta})$ with respect to $\boldsymbol{\theta}_{v}$, we need to solve $0\in\partial C^{\prime}_{v}$, which is equivalent to

[TABLE]

where

[TABLE]

and where $s_{v}$ and $t_{v}$ satisfy:

[TABLE]

The first task is to obtain necessary and sufficient conditions for which the solution $\boldsymbol{\theta}_{v}=0$ is the optimal one. Let

[TABLE]

Then the solution to Equation (20) is zero if and only if $J^{*}\leq\gamma^{\prime 2}_{v}$. Calculating $J^{*}$ is a ridge regression problem that can be easily solved (see Propositions 8.2 and 8.3 in Section 8.7).

If the solution to Equation (20) is not $\boldsymbol{\theta}_{v}=0$, the problem is to solve the subgradient equation:

[TABLE]

Because $\boldsymbol{\theta}_{v}$ appears in both sides of the equation, a numerical procedure is needed (see Proposition 8.4).

Other cases

• •

If all the $\mu^{\prime}_{v}$ are equal to 0, the parameters $\boldsymbol{\theta}_{v}$ are not identifiable.

• •

If all the $\gamma^{\prime}_{v}$ are equal to 0, then we have to solve a classical group-lasso problem with respect to the parameters $\boldsymbol{\theta}^{\prime}_{v}$ defined as $\boldsymbol{\theta}^{\prime}_{v}=K_{v}^{1/2}\boldsymbol{\theta}_{v}$ for all $v\in\mathcal{P}$.

• •

Let $v$ such that $\mu^{\prime}_{v}=0$, $\gamma^{\prime}_{v}\neq 0$ and assume that at least one of the $\mu^{\prime}_{w}$ is non zero for $w\in\mathcal{P}$, $w\neq v$. Then it is shown in Proposition 8.1 that $\boldsymbol{\theta}_{v}=0$ if and only if $2\|R_{v}\|\leq\gamma^{\prime}_{v}$.

• •

In the same way, if $\gamma^{\prime}_{v}=0$, and $\mu^{\prime}_{v}\neq 0$, then $\boldsymbol{\theta}_{v}=0$ if and only if $2\|K_{v}^{1/2}R_{v}\|\leq\mu^{\prime}_{v}$.

Finally the algorithm is the following:

1. (1)

Start with an initial value $\boldsymbol{\theta}=\boldsymbol{\theta}_{0}$. 2. (2)

Calculate $f_{0}$ using Equation (19). For the group $v$, calculate $\mathbf{R}_{v}$ and determine if the group $v$ should be excluded or not either by solving the problem defined at Equation (21) if $\mu^{\prime}_{v}\neq 0$ and $\gamma^{\prime}_{v}\neq 0$, or directly if one of them equals 0. If it is the case, set $\boldsymbol{\theta}_{v}=0$. If not, solve Equation (22) to obtain $\boldsymbol{\theta}_{v}$. 3. (3)

Iterate step 2. over all the groups $v$. 4. (4)

Iterate step 2. and 3. until convergence.

5.2. Choice of the tuning parameters

For each value of the tuning parameters $(\mu^{\prime}_{v},\gamma^{\prime}_{v}),v\in\mathcal{P}$, the algorithm provides an estimate of $m$ and of the Sobol indices. The problem for choosing these parameters values is crucial. We propose to restrict this choice by considering tuning parameters proportionnal to known weights: for all $v\in\mathcal{P}$, $\mu^{\prime}_{v}=\mu\omega_{v}$ and $\gamma^{\prime}_{v}=\gamma\zeta_{v}$, where the weights $\omega_{v}$ and $\zeta_{v}$ are fixed. For example, one can take weights that increase with the cardinal of $v$ in order to favour effects with small interaction order between variables. Or, according to the theoretical result given at Theorem 4.1, we can choose $\omega_{v}=\widehat{\nu}^{2}_{n,v}$ and $\zeta_{v}=\widehat{\nu}_{n,v}$, where $\widehat{\nu}_{n,v}$ is an estimate of $\nu_{n,v}$ based on the eigenvalues of the matrix $K_{v}$. Any other choice, depending on the problem of interest, may be relevant.

Once the weights are chosen, we estimate $m$, on the basis of a learning data set $(Y_{i},\mathbf{X}_{i}),i=1,\ldots n$, for a grid of values of $(\mu,\gamma)$. We first set $\gamma=0$, and calculate $\mu_{\max}$ the smallest value of $\mu$ such that the solution to the minimization of

[TABLE]

is $\boldsymbol{\theta}_{v}=0$ for all $v\in\mathcal{P}$. Then we can consider $\mu_{\ell}=\mu_{\max}2^{-\ell}$ for $\ell\in\left\{1,\ldots,\ell_{\max}\right\}$, as a grid of values for $\mu$. The grid of values for $\gamma$ may be chosen after few attempts.

For choosing the final estimator, say $\widehat{f}$, we suppose that we have at our disposal a testing data set $(Y^{T}_{i},\mathbf{X}^{T}_{i}),i=1,\ldots n^{T}$, and we propose two procedures.

**Proc. GS: **

The first one uses the testing data set for estimating the prediction error. Precisely, for each value of $(\mu,\gamma)$ in the grid, let $\widehat{f}_{(\mu,\gamma)}(\cdot)$ be the estimation of $m$ obtained with the learning data set. Then

[TABLE]

estimates the prediction error, and we propose to choose the pair $(\mu,\gamma)$ that minimizes $\mathrm{PE}(\mu,\gamma)$, say $(\widehat{\mu},\widehat{\gamma})$. Finally the estimator, denoted $\widehat{f}^{\mathrm{GS}}$ is defined as $\widehat{f}^{\mathrm{GS}}=\widehat{f}_{(\widehat{\mu},\widehat{\gamma})}$. In the following, we will refer to this procedure as the Group-Sparse procedure.

**Proc. rdg: **

Doing the parallel with the inconsistency of the lasso for estimating the support of the parameters in the classical regression problem, we propose to choose the tuning parameter that minimizes the risk of the ridge estimator over the support estimated by the ridge-group-sparse procedure. Indeed, if the tuning parameter is chosen to minimize the prediction error, the lasso is not consistent for support estimation (see [24] for example). One idea to overcome this problem, is to choose the tuning parameter that minimizes the risk of the Gauss-lasso estimator which is calculated in two steps: For a given value of the tuning parameter, the estimation of the support of the parameter is estimated using a lasso procedure, then the least square estimator over this support is calculated. When the objective is support estimation, some numerical simulations [33] and theoretical results  [18] suggest that it may be more advisable not to apply the selection schemes based on prediction risk to the lasso estimators, but rather to the Gauss-lasso estimators. Our procedure, called the ridge procedure, applies the same idea in the framework of sparse nonparametric estimation. Precisely it considers the collection of supports composed of the different $S_{\widehat{f}_{(\mu,\gamma)}}$ when $(\mu,\gamma)$ belongs to the grid. For each support $S$ in this collection, we estimate $f$ using a ridge procedure assuming that the support of $f$ is $S$: for a given $\lambda$, $f^{\mathrm{rdg}}_{\lambda,S}$ is defined as follows:

[TABLE]

We choose a grid of values for $\lambda$, and for each $S$ in the collection and $\lambda$ in the grid, we estimate the prediction error $\mathrm{PE}(\lambda,S)$ defined as follows:

[TABLE]

For each $S$ in the collection, let $\widehat{\lambda}(S)$ be the minimizer of $\mathrm{PE}(\lambda,S)$ when $\lambda$ varies in the grid, and let $\widehat{S}$ be the minimizer of $\mathrm{PE}(\widehat{\lambda}(S),S)$, then the estimator denoted $\widehat{f}^{\mathrm{rdg}}$ is defined as $\widehat{f}^{\mathrm{rdg}}=f^{\mathrm{rdg}}_{\widehat{\lambda}(\widehat{S}),\widehat{S}}$.

If a testing data set is not available, we can use the classical V-fold cross validation (see [1] for example) either to estimate $\mathrm{PE}(\mu,\gamma)$ or to estimate $\mathrm{PE}(\lambda,S)$.

6. Simulation study

In order to evaluate the performances of our method for estimating a meta-model and sensitivity indices of a function $m$ we carried out a simulation study. We consider the $g$-function of Sobol defined on $[0,1]^{d}$ as

[TABLE]

whose Sobol indices can be expressed analytically (see Saltelli et al. [34]). Following the simulation experiment proposed by Durrande et al. [12], we take $d=5$ and $(c_{1},c_{2},c_{3},c_{4},c_{5})=(0.2,0.6,0.8,100,100)$. The lower the value of $c_{a}$, the more significant the variable $x_{a}$. The variables $X_{a},a=1,\ldots,d$ are independent and uniformly distributed on $[0,1]$. We consider the regression model $Y_{i}=m(\mathbf{X}_{i})+\sigma\varepsilon_{i}$, for $i=1,\ldots,n$, with $\mathcal{N}(0,1)$ independent error terms $\varepsilon_{i}$.

Simulation design

We present the results for $n\in\left\{50,100,200\right\}$, $\sigma\in\left\{0,0.2\right\}$. For all $a=1,\ldots,d$, the kernels $k_{a}$ are the same: we considered the Brownian kernel, $k^{b}(x,x^{\prime})=1+\min\left\{x,x^{\prime}\right\}$, the Matérn kernel, $k^{m}(x,x^{\prime})=(1+2|x-x^{\prime}|)\exp(-2|x-x^{\prime}|)$, and the Gaussian kernel, $k^{g}(x,x^{\prime})=\exp(x-x^{\prime})^{2}$.

For each simulation, we generate three independent data sets as follows : a Latin Hypercube Sample of the inputs is simulated to give the matrix $\mathbf{X}$ with $n$ rows and $d=5$ columns, and a $n$-sample of independent centered Gaussian variable with variance 1 is simulated. This operation is repeated three times in order to obtain the learning and testing data sets and a third data set for estimating the estimators performances. As explained in Section 5.2, we choose optimal values of the tuning parameters, $(\mu,\gamma)$ by minimizing a prediction error $\mathrm{PE}$, and get an estimator $\widehat{f}$ of $m$, as well as estimates of the Sobol indices:

[TABLE]

Let $\Omega_{v}$ be the matrix whose components satisfy

[TABLE]

The estimator $\widehat{\mathop{\rm Var}\nolimits}\left(m_{v}(\mathbf{x}_{v})\right)$ is calculated as follows:

[TABLE]

We also propose to estimate these quantities by their empirical variances as in Equation (8).

Performance indicators

To evaluate the performances of our method for estimating a meta-model, we use the classical coefficient of determination $R^{2}$ estimated using the third data set $(Y_{i}^{P},\mathbf{X}^{P}_{i}),i=1,\ldots,n$:

[TABLE]

Moreover we calculate the empirical risk $\mathrm{ER}=\|m-\widehat{f}\|_{n}^{2}$. For each simulation $s$, we get $R^{2}_{s}$ and $\mathrm{ER}_{s}$ and we report the means of these quantities over all simulations.

Similarly, for each $v$, and each simulation $s$, we get $\widehat{S}_{v,s}$ and we report its mean, $\widehat{S}_{v,\cdot}$, its estimated standard-error, and to sum up the behaviour of our procedure for estimating the sensitivity indices, we estimate the global error, denoted $\mathrm{GE}$, defined as follows

[TABLE]

In order to assess the performances of our procedure for selecting the greatest Sobol Indices, precisely those that are greater than some small quantity as $\rho=10^{-4}$, we calculate for each group $v\in\mathcal{P}$, the percentage of simulations for which $v$ is in the support of the estimator $\widehat{f}$. Then we average these quantities, on one hand over groups $v$ such that $S_{v}>\rho$, and on the other hand, over groups $v$ such that $S_{v}\leq\rho$. Let us denote these quantities ${\mathrm{pSel}}_{S_{v}>\rho}$ and ${\mathrm{pSel}}_{S_{v}\leq\rho}$ respectively.

For each $(\mu,\gamma)$ in a grid of values, the estimator $\widehat{f}_{\mu,\gamma}$ is defined as the minimizer of the criteria given at Equation (18) taking $\omega_{v}=\zeta_{v}=1$ for all $v\in\mathcal{P}$. In order to save computation time, we restrict the optimisation to sets $v$ such that $|v|\leq 3$. Some preliminary simulations showed that the terms coresponding to $|v|\geq 4$ are nearly always equal to 0.

Choosing the tuning parameters

Let us begin with the comparison of the two methods proposed for choosing the final estimator, see Section 5.2. The results are given in Tables 1 and 2. It appears that the procedure based on the ridge estimator after selection of the groups outperforms the method based on the group-sparse estimator. As expected both methods perform better when $n$ increases, and when $\sigma=0$.

Similarly, the Sobol indices are better estimated, in the sense of the global error, with the procedure based on the ridge estimate of the metamodel, see Table 3. The means of the estimators for the Sobol indices greater than $\rho=10^{-4}$ are given in Tables 4 and 5. It appears than $S_{\{1\}}$ is over-estimated using the procedure based on the group-sparse estimator, leading to under-estimate the Sobol indices associated with interactions of order 2. This tendancy is much less pronounced with the procedure based on the ridge estimator.

Let us now consider the performances of the procedure for selecting the non zero Sobol indices. In Tables 6 and 7 we report the percentages of simulations for which the Sobol indices smaller (respectively greater) than $\rho$ are selected, and for which each of Sobol index greater than $\rho$ is selected. From these results, we conclude that the procedure based on the ridge estimator is more strict for selecting non-zero Sobol indices.

Comparing different kernels

Finally we compare the performances of the procedures for different kernels, see Tables 8 and 9. The means of the estimated empirical risk, $\mathrm{ER}$, and of the global error for estimating the sensitivity indices, $\mathrm{GE}$, are calculated for each kernel. It appears that the Matérn kernel gives the best results, except for the case $n=50$ and $\sigma=0$ where the empirical risk of $\widehat{f}^{\mathrm{GS}}$ is smaller for the Brownian kernel.

In practice, one may want to choose the kernel according to the smallest prediction error. For that purpose, we propose to calculate the estimators $\widehat{f}^{\mathrm{GS}}$ and/or $\widehat{f}^{\mathrm{rdg}}$ for each kernel, as described in Section 5.2, as well as their associated prediction errors. Then we choose the kernel for which the prediction error is minimized. The results are reported under the column “mixed” in Tables 8 and 9. It appears that the estimated empirical risks for this “mixed” procedure are nearly equal to the minimum estimated empirical risks over the different kernels.

7. Sketch of proof of Theorem 4.1

We give here a sketch of the proof and we postpone to Section 8 for complete statements. In particular, we denote by $C$ constants that vary from an equation to the other, and we assume that $\sigma=1$.

The proof of Theorem 4.1 starts in the same way as the proof of Theorem 1 in Raskutti et al. [31]. Nevertheless it differs in several points, in particular because the terms occuring in the decomposition of functions in $\mathcal{H}$ depend on several variables and thus are not independent. Indeed, $f_{v}(\mathbf{X}_{v})$ and $f_{v}^{\prime}(\mathbf{X}_{v^{\prime}})$ are not independent as soon as the groups $v$ and $v^{\prime}$ share some of the variables $X_{a},a=1,\ldots,d$. Moreover, we do not assume that the function $m$ is in $\mathcal{F}$.

Starting from the definition of $\widehat{f}$, some simple calculation (see Equation (28)) give that for all $f\in\mathcal{F}$

[TABLE]

If we set $g=\widehat{f}-f$, then $g\in\mathcal{H}$, $g=g_{0}+\sum_{v}g_{v}$, with $g_{v}=\widehat{f}_{v}-f_{v}$, and for each $v$, $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$.

The main problem is now to control the empirical process. For each $v$, letting $\lambda_{n,v}$ as in (13), we state (see Lemma 8.1, page 8.1) that, with high probability,

[TABLE]

Therefore, if for all $v$, $\mu_{v}$ and $\gamma_{v}$ satisfy Equation (14), we deduce that with high probability (setting $g=\widehat{f}-f$)

[TABLE]

Besides we can express the decomposability property of the penalty as follows (see lemma 8.2, page 8.2): with high probability (in the set where the empirical process is controled as stated above),

[TABLE]

Putting the things together, and noting again that $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$, we obtain the following upper bound

[TABLE]

The last important step consists in comparing $\sum_{v\in S_{f}}\|g_{v}\|_{n}$ to $\|\sum_{v\in S_{f}}g_{v}\|_{n}$. More precisely, it can be shown (see lemma 8.3, page 8.3) that for all $v\in\mathcal{P}$, with high probability, we have

[TABLE]

Using the orthogonality assumption between the spaces $\mathcal{H}_{v}$, we have $\sum_{v\in S_{f}}\|g_{v}\|^{2}_{\mathbb{L}^{2}(P_{\mathbf{X}})}=\|\sum_{v}g_{v}\|^{2}_{\mathbb{L}^{2}(P_{\mathbf{X}})}$, and thus we get

[TABLE]

Finally it remains to consider different cases according to the rankings of $\|\widehat{f}-f\|^{2}_{\mathbb{L}^{2}(P_{\mathbf{X}})}$, $\|\widehat{f}-f\|^{2}_{n}$ and $\sum_{v\in S_{f}}\mu_{v}+\gamma^{2}_{v}$ to get the result of Theorem 4.1.

8. Proofs

Recall that we cconsider the regression model defined at Equation (1), where $\mathbf{X}$ has distribution $P_{\mathbf{X}}=P_{1}\times\ldots\times P_{d}$ defined on $\mathcal{X}$ a compact subset of $\mathbb{R}^{d}$, and $\varepsilon$ is distributed as $\mathcal{N}(0,\sigma^{2})$. We denote by $P_{\mathbf{X},\boldsymbol{\varepsilon}}$ the distribution of $(\mathbf{X},\boldsymbol{\varepsilon})$. We observe a $n$ sample $(Y_{i},\mathbf{X}_{i}),i=1,\ldots,n$ with law $P_{\mathbf{X},\boldsymbol{\varepsilon}}$.

The notation and the procedure are given in Sections 2 and 4.

Let us add on few notations that will be used along the proofs.

For $v\in\mathcal{P}$ we denote $|v|$ the cardinal of $v$. For a function $\phi:\mathbb{R}^{|v|}\mapsto\mathbb{R}$, we denote $V_{n,\varepsilon}$ the empirical process defined by

[TABLE]

For the sake of simplicity we assume $\sigma=1$. Moreover, we set $R^{\prime}=1$, see (10). Consequently, for any function $f\in\mathcal{H}$, $\|f_{v}\|_{\mathcal{H}_{v}}\leq 1$, and $\|f_{v}\|_{\infty}\leq\|f\|_{\mathcal{H}_{v}}$. The proofs can be done exactly in the same way by considering the general case. In the proofs, the $\|\cdot\|_{\mathbb{L}^{2}(P_{\mathbf{X}})}$ norm will be denoted by $\|\cdot\|_{2}$.

8.1. Proof of Theorem 4.1

The proof is based on four main lemmas proved in Section 8.5. In Section 8.4 other lemmas used all along the proof are stated. Their proof are postponed to Section 8.6.

Let us first establish inequalities that will be used in the following. Let $f\in\mathcal{H}$ and $v\in S_{f}$ (see (9)).

Using that for any $v\in S_{f}$, and any norm $\|\cdot\|$ in $\mathcal{H}_{v}$, $\|f_{v}\|-\|\widehat{f}_{v}\|\leq\|f_{v}-\widehat{f}_{v}\|$ and that for any $v\notin S_{f}$, $\|f_{v}\|=0$, we get that

[TABLE]

Combining (24), and (25), to the fact that for any function $f\in\mathcal{H}$, $\mathcal{L}(\widehat{f})\leq\mathcal{L}(f)$, we obtain that

[TABLE]

with

[TABLE]

If $\|m-f\|_{n}^{2}\geq B$, we immediately get the result since in that case

[TABLE]

If $\|m-f\|_{n}^{2}<B$, we get that

[TABLE]

The control of the empirical process |V_{n,\varepsilon}\big{(}\widehat{f}-f\big{)}| is given by the following lemma (proved in Section 8.5.1, page 8.5.1).

** Lemma 8.1****.**

Let $V_{n,\varepsilon}$ be defined in (23). For any $f$ in $\mathcal{F}$, we consider the event $\mathcal{T}$ defined as

[TABLE]

where the quantities $\lambda_{n,v}$ are defined by Equation (13) and where $\kappa=10+4\Delta$. Then, for some positive constants $c_{1},c_{2}$,

[TABLE]

Conditionning on $\mathcal{T}$, Inequality (28) becomes

[TABLE]

which may be decomposed as follows

[TABLE]

If we choose $C_{1}\geq\kappa$ in Theorem 4.1, then $\kappa\lambda_{n,v}^{2}\leq\mu_{v}/2$ and $\kappa\lambda_{n,v}\leq\gamma_{v}/2$ and the previous inequality becomes

[TABLE]

Next we use the decomposability property of the penalty expressed in the following lemma (proved in Section 8.5.2 page 8.5.2).

** Lemma 8.2****.**

For any $f\in\mathcal{F}$, under the assumptions of Theorem 4.1 with $C_{1}\geq\kappa$, conditionnally on $\mathcal{T}$, see (29), we have

[TABLE]

Hence, by combining (30) and Lemma (8.2) we obtain

[TABLE]

For each $v$, $\|\widehat{f}_{v}-f_{v}\|_{\mathcal{H}_{v}}\leq 2$ (because the functions $\widehat{f}_{v}$ et $f_{v}$ belong to the class $\mathcal{F}$, see (5)), and consequently, for some constant $C$,

[TABLE]

To finish the proof it remains to compare the two quantities $\sum_{v\in S_{f}}\|\widehat{f}_{v}-f_{v}\|_{n}^{2}$ and $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{n}^{2}$. For that purpose we show that $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{n}$ is less than $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{2}^{2}$ plus an additive term coming from concentration results (see the Lemma given below). Next, thanks to the orthogonality of the spaces $\mathcal{H}_{v}$ with respect to $\mathbb{L}^{2}(P_{\mathbf{X}})$, $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{2}^{2}=\sum_{v\in S_{f}}\|\widehat{f}_{v}-f_{v}\|_{2}^{2}$. To conclude, it remains to consider several cases, according to the rankings of $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{2}^{2}$, $\|\sum_{v\in S_{f}}\widehat{f}_{v}-f_{v}\|_{n}^{2}$, and $d^{2}(f)$. This is the subject of the following lemma whose proof is given in Section 8.5.3, page 8.5.3.

** Lemma 8.3****.**

For $f\in\mathcal{H}$, let $\mathcal{A}$ be the event

[TABLE]

Then, for some positive constant $c_{2}$,

[TABLE]

On the set $\mathcal{A}$, Inequality (31) provides that, for all $K>0$

[TABLE]

Inequality (8.1) uses the inequality $2ab\leq\frac{1}{K}a^{2}+Kb^{2}$ for all positive $K$, and Inequality (35) uses the orthogonality with respect to $\mathbb{L}^{2}(P_{\mathbf{X}})$.

In the following we have to consider several cases, according to the rankings of $\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$, $\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}$ and $d(f)$ defined as follows

[TABLE]

More precisely, we consider three cases

• Case 1:

$\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}.$

• Case 2:

$\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}\leq d(f)$

• Case 3:

$\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$ and $d(f)\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$.

Case 1: From (35), for any $f\in\mathcal{H}$, we get

[TABLE]

Hence, using that for all $K^{\prime}>0$,

[TABLE]

we obtain for a suitable choice of $K^{\prime}$, say $1+K^{\prime}<K/C$, that, for some positive constant $C^{\prime}$,

[TABLE]

This shows the result in Case 1.

Case 2: Inequality (35) becomes

[TABLE]

which gives the expected result since $d^{2}(f)=\max\left\{\sum_{v\in S_{f}}\mu_{v},\sum_{v\in S_{f}}\gamma_{v}^{2}\right\}$.

Case 3: Recall that in this case, $\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$ and $d(f)\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$. This case is solved by applying the following Lemma (proved in Section 8.5.4, page 8.5.4), which states that with high probability, $\|\widehat{f}-f\|_{2}\leq\sqrt{2}\|\widehat{f}-f\|_{n}$.

** Lemma 8.4****.**

Let $f=\sum_{v}f_{v}\in\mathcal{F}$ with support $S_{f}$, $d(f)$ be defined by (36), and let $\mathcal{G}(f)$ be the class of functions written as $g=\sum_{v\in\mathcal{P}}g_{v}$, such that $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$ satisfying for all $f\in\mathcal{F}$

[TABLE]

Then the event

[TABLE]

has probability greater than $1-\exp(-nc_{3}\sum_{v\in S_{f}}\lambda_{n,v}^{2})$ for some constant $c_{3}$.

Note that Assumption $n\lambda_{n,v}^{2}\geq-C_{2}\log(\lambda_{n,v})$ implies that $\lambda_{n,v}=K_{n,v}/\sqrt{n}$ with $K_{n,v}\rightarrow\infty$. Then, if $f$ is such that $|S_{f}|\geq 1$, $\exp(-nc_{3}\sum_{v\in S_{f}}\lambda_{n,v}^{2})\leq\exp(-c_{3}\min_{v\in\mathcal{P}}K_{n,v})$. If $f$ is such that $|S_{f}|=0$, then Condition C1 is not satisfied except if $g_{v}=0$ for all $v\in\mathcal{P}$. Because we will apply Lemma 8.4 to $g_{v}=\widehat{f}_{v}-f_{v}$, this event has probability 0. Therefore the event

[TABLE]

has probability greater than $1-\eta/3$ for some $0<\eta<1$.

Conditionning on the events $\mathcal{T}$ and $\mathcal{A}$ (defined by (29) and (32)), and according to Lemma 8.2, $\sum_{v\in\mathcal{P}}(\widehat{f}_{v}-{f}_{v})$ belongs to the set $\mathcal{G}(f)$. According to (35), we conclude in the same way as in Case 1.

Finally, it remains to quantify $P_{\mathbf{X},\boldsymbol{\varepsilon}}(\mathcal{T}\cap\mathcal{A}\cap\mathcal{C})$. Following Lemma 8.1, and Lemma 8.3, $\mathcal{T}$, respectively $\mathcal{A}$, has probability greater than $1-c_{1}\sum_{v\in\mathcal{P}}\exp(-nc_{2}\lambda_{n,v}^{2})$, respectively $1-\sum_{v}\exp(-n\gamma_{v}^{2})$. Each of these probabilities is greater $1-\eta/3$ thanks to the assumption $n\lambda_{n,v}^{2}\geq-C_{2}\log\lambda_{n,v}$.

8.2. Proof of Corrolary 4.1

We start from Theorem 4.1 which states that with high probability,

[TABLE]

and use that for all $\theta>0$,

[TABLE]

For $d$ defined by (36), we consider once again the three cases defined page 8.1.

Case 1: According to (37) and (39), we get the result since

[TABLE]

Case 2: We directly obtain that

[TABLE]

Case 3: Recall that in this case, $\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{n}\leq\|\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}\|_{2}$ and $d(f)\leq\|\sum_{v\in S}\widehat{f}_{v}-{f}_{v}\|_{2}$. Apply Lemma 8.4 (page 8.4) and conclude that conditionning on the events $\mathcal{T}$ and $\mathcal{A}$, defined by (29) and (32), then $\sum_{v\in\mathcal{P}}\widehat{f}_{v}-{f}_{v}$ belongs to $\mathcal{G}(f)$ defined in Lemma 8.4. Now, conditionning on the event $\mathcal{C}$ we get the result since

[TABLE]

$\Box$

8.3. Rate of convergence

Recall that we consider the case where the variables $X_{1},\ldots,X_{d}$ have the same distribution $P_{1}$ on $\mathcal{X}_{1}\subset\mathbb{R}$, and where the unidimensionnal kernels $k_{0a}$ are all identical.

In this context, our goal is to show that the rate $\nu_{n,v}$ defined at Equation (12) is bounded above by a term of order $n^{-\alpha/(2\alpha+1)}(\log n)^{\gamma}$, where $\gamma\geq(|v|-1)\alpha/(2\alpha-1)$.

We start from the fact that, $k_{v}(\mathbf{x}_{v},\mathbf{x}^{\prime}_{v})=\prod_{a\in v}k_{0}(x_{a},x^{\prime}_{a})$, with a kernel $k_{0}$ admitting an eigen expansion given by

[TABLE]

where the eigenvalues $\omega_{0,\ell}$ are non negative and ranged in the decreasing order at the rate $\ell^{-2\alpha}$ for some $\alpha>1/2$, and where the $\zeta_{\ell}$ are the associated eigen functions, orthonormal with respect to $\mathbb{L}^{2}(P_{1})$.

Therefore the kernel $k_{v}$ admits the following expansion

[TABLE]

According to this expansion the $\omega_{v,\boldsymbol{\ell}}$ are of order $(\prod_{a=1}^{|v|}\ell_{a})^{-2\alpha}$.

In order to control the rate $\nu_{n,v}$ defined as

[TABLE]

, we have to calculate an upper bound for $Q_{n,v}^{2}(t)$ We start with the following inequalities that hold up to some constant, for $t^{-1/\alpha}>1$

[TABLE]

Now let us mention that $\alpha>1/2$, $\sum_{j\geq 1}\frac{1}{j^{2\alpha}}$ is a constant that depends on $\alpha$. We thus focus on the first term in the right hand side of Equation (40).

Let $u=t^{-1/\alpha}\geq 1$ and let $B_{|v|}$ be defined as follows:

[TABLE]

Let us prove that

[TABLE]

Proof of Equation (41): First note that

[TABLE]

In the same way,

[TABLE]

More generally,

[TABLE]

Let $A_{|v|}$ be defined as follows:

[TABLE]

then we get $B_{|v|}\leq uA_{|v|-1}$. If we show that

[TABLE]

then Inequality (41) is proved.

Proof of Equation (42) :

[TABLE]

[TABLE]

In the same way we have

[TABLE]

And Bound (42) is proved.

Rate of convergence

Let us come back to the control of the rate $\nu_{n,v}=\inf\left\{t\mbox{ such that }Q_{n,v}(t)\leq\Delta t^{2}\right\}$. Thanks to (41) we obtain that, up to some constant that depends on $|v|$ and $\alpha$,

[TABLE]

It remains now to find $t$ such that, up to constant

[TABLE]

If $t=n^{-\beta}(\log(n))^{\gamma}$ with $\beta=\alpha/(1+2\alpha)$, $\gamma>0$, $\alpha>1/2$, then

[TABLE]

Therefore $\nu_{n,v}$ will be smaller than the infimum of $t$ such that

[TABLE]

which is satisfied if $\gamma\geq(|v|-1)\alpha/(2\alpha-1)$.

8.4. Intermediate Lemmas

For $v\in\mathcal{P}$, let $\mathcal{H}_{v}$ be the RKHS associated to the self reproducing kernel $k_{v}$. Let $Q_{n,v}$ and $\nu_{n,v}$ and be defined by Equations (11) and (12). For any function $g_{v}\in\mathcal{H}_{v}$, let $V_{n,\varepsilon}$ be defined at Equation (23) and consider the following processes

[TABLE]

** Lemma 8.5****.**

If $E_{\mathbf{X},\boldsymbol{\varepsilon}}$ denotes the expectation with respect to the distribution of $(\mathbf{X},\varepsilon)$, we have for all $t>0$,

[TABLE]

Its proof is given in Section 8.6.1 page 8.6.1.

** Lemma 8.6****.**

Let $b>0$ and let $\mathcal{G}(t)$ be the following class of functions:

[TABLE]

Let $\Omega_{v,t}$ be the event defined as

[TABLE]

Then for any $t\geq\nu_{n,v}$, the event $\Omega_{v,t}$ has probability greater than $1-\exp(-c_{2}nt^{2})$, for some positive constant $c_{2}$.

Its proof is given in Section 8.6.2, page 8.6.2.

** Lemma 8.7****.**

For any function $g_{v}\in\mathcal{H}_{v}$ satisfying $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$, $\|g_{v}\|_{\infty}\leq b$ and $\|g_{v}\|_{2}\geq t$, for all $t\geq\nu_{n,v}$ and $b\geq 1$, the event

[TABLE]

has probabilty greater than $1-\exp(-c_{2}nt^{2})$ for some positive constant $c_{2}$.

Its proof is given in Section 8.6.3, page 8.6.3.

** Lemma 8.8****.**

If $E_{\boldsymbol{\varepsilon}}$ denotes the expectation with respect to the distribution of $\varepsilon$, we have

[TABLE]

Its proof is given in Section 8.6.4, page 8.6.4.

** Lemma 8.9****.**

Conditionnaly on the space $\Omega_{v,t}$ defined by (46), we have the two following inequalities:

[TABLE]

Its proof is given in Section 8.6.5, page 8.6.5.

** Lemma 8.10****.**

Let $\lambda_{n,v}$ be defined at Equation (13), $\Delta$ at Equation (12) and $\kappa=10+4\Delta$. Conditionnaly on the space $\Omega_{v,\lambda_{n,v}}$ defined at Equation (46), for some positive constants $c_{1},c_{2}$, with probability greater than $1-c_{1}\exp(-c_{2}n\lambda_{n,v}^{2})$, we have

[TABLE]

Its proof is given in Section 8.6.6, page 8.6.6.

8.5. Proofs of Lemma 8.1 to

 8.4:

8.5.1. Proof of Lemma 8.1 (page 8.1)

For $f\in\mathcal{F}$ and $v\in\mathcal{P}$, let $g_{v}=\widehat{f}_{v}-f_{v}$. Note that $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$. Let us show that

[TABLE]

We start by writing that

[TABLE]

Consider the two following cases:

• Case A:

$\|g_{v}\|_{n}\leq\lambda_{n,v}\|g_{v}\|_{\mathcal{H}_{v}}$

• Case B:

$\|g_{v}\|_{n}>\lambda_{n,v}\|g_{v}\|_{\mathcal{H}_{v}}$

Case A: Since $\|g_{v}\|_{n}\leq\lambda_{n,v}\|g_{v}\|_{\mathcal{H}_{v}}$, we have

[TABLE]

We then apply Lemma 8.10, page 8.10, and conclude that (51) holds in Case A for each $v\in\mathcal{P}$ since, with high probability

[TABLE]

Case B: Consider now the case $\|g_{v}\|_{n}>\lambda_{n,v}\|g_{v}\|_{\mathcal{H}_{v}}$ and let us show that for any $v\in\mathcal{P}$,

[TABLE]

Let $r_{v}$ be a deterministic number such that $r_{v}>\lambda_{n,v}$. Our first step relies on the study of the process $W_{n,n,v}(r_{v}),$ for $r_{v}>\lambda_{n,v}$. In that case we state two results:

• R1

For any deterministic $r_{v}\geq\lambda_{n,v}$, with probability greater than $1-c_{1}\exp(-c_{2}n\lambda_{n,v}^{2})$,

[TABLE]

• R2

Inequality (54) continues to hold for random $r_{v}$ of the form

[TABLE]

Combining these two points implies that, with probability greater than $1-c_{1}\exp(-c_{2}n\lambda_{n,v}^{2})$,

[TABLE]

Consequently, in Case B, according to (52), for each $v$, Inequality (51) holds because

[TABLE]

This ends up the proof of Lemma 8.1.

Proof of R1

Taking $t=r_{v}$ and $\delta=\lambda_{n,v}$ in (47), with probability greater than $1-4\exp(-n\lambda_{n,v}^{2})$, we have

[TABLE]

Next we prove that for some positive $r_{v}$, with probability greater than $1-\exp(-nc\lambda_{n,v}^{2})$, we have

[TABLE]

Let $\widehat{\nu}_{n,v}$ defined as the smallest solution of $E_{\boldsymbol{\varepsilon}}[W_{n,n,v}(t)]\leq\kappa t^{2}$. For $W_{n,n,v}$, defined by (43), we write

[TABLE]

Besides, Lemma 8.10 stated that on the event $\Omega_{v,\lambda_{n,v}}$, $E_{\boldsymbol{\varepsilon}}W_{n,n,v}(\lambda_{n,v})\leq\kappa\lambda_{n,v}^{2}$. It follows from the definition of $\widehat{\nu}_{n,v}$, and Lemma 8.6, that $\widehat{\nu}_{n,v}\leq\lambda_{n,v}$ for all $v\in\mathcal{P}$ with probability greater than $1-\exp(-nc_{2}\sum_{v\in\mathcal{P}}\lambda_{n,v}^{2})$. Consequently, for any deterministic $r_{v}$ such that $r_{v}\geq\lambda_{n,v}$, (54) is satisfied with high probability.

Proof of R2

Let us prove R2 by using a peeling-type argument. Our aim is to prove that (54) holds for any $r_{v}$ of the form

[TABLE]

Since $\|g_{v}\|_{\infty}/\|g_{v}\|_{\mathcal{H}_{v}}\leq 1$, we have $\|g_{v}\|_{n}/\|g_{v}\|_{\mathcal{H}_{v}}\leq 1$. We thus restrict ourselves to $r_{v}$ satisfying $r_{v}=\|g_{v}\|_{n}/\|g_{v}\|_{\mathcal{H}_{v}}$ with $\|g_{v}\|_{n}/\|g_{v}\|_{\mathcal{H}_{v}}\in(\lambda_{n,v},1]$.

We start by splitting the interval $(\lambda_{n,v},1]$ into $M$ disjoint intervals such that $(\lambda_{n,v},1]=\cup_{k=1}^{M}(2^{k-1}\lambda_{n,v},2^{k}\lambda_{n,v}]$, for some $M$ that will be chosen later. Consider the event $\mathcal{D}^{c}$ defined as follows:

[TABLE]

We prove that, for some positive constants $c_{1},c_{2}$,

[TABLE]

For $\overline{g}_{v}\in\mathcal{D}^{c}$, let $\overline{k}$ be the integer in $\{1,\cdots,M\}$, such that

[TABLE]

This $\overline{k}$ satisfies

[TABLE]

Therefore, we get

[TABLE]

By taking $r_{v}=2^{\overline{k}}\lambda_{n,v}$ in (54), we have

[TABLE]

Now let us write $\mathcal{D}^{c}$ as follows:

[TABLE]

The set $\mathcal{D}^{c}$ has probability smaller than $c_{1}M\exp(-c_{2}n\lambda_{n,v}^{2})$. If we choose $M$ such that $\log M\leq(c_{2}/2)n\lambda_{n,v}^{2}$, then the probability of the set $\mathcal{T}$ is greater than

[TABLE]

It follows that R2 is proved which ends up the proof of Lemma 8.1. $\Box$

8.5.2. Proof of Lemma 8.2 (page 8.2).

Starting from (27) with $B$ defined by Equation (26), we write

[TABLE]

On the event $\mathcal{T}$ defined in (29) we have

[TABLE]

Rearranging the terms we obtain that

[TABLE]

Now, thanks to Assumption (14) with $C_{1}\geq\kappa$ we have $\kappa\lambda_{n,v}^{2}\leq\mu_{v}$ and $2\kappa\lambda_{n,v}\leq\gamma_{v}$ and Lemma 8.2 is shown since

[TABLE]

$\Box$

8.5.3. Proof of Lemma 8.3

(page 8.3):

Let us consider the following two cases:

• •

$\|\widehat{f}_{v}-{f}_{v}\|_{2}\leq\gamma_{v}$. We apply Lemma 8.6 (page 8.6) to the function $g_{v}=\widehat{f}_{v}-{f}_{v}$. It satisfies $g_{v}\in\mathcal{G}(\gamma_{v})$ with $b=2$ (recall that $\|\cdot\|_{\infty}\leq\|\cdot\|_{\mathcal{H}_{v}}$). Moreover, $\gamma_{v}\geq C_{1}\lambda_{n,v}\geq C_{1}\nu_{nv}\geq\nu_{n,v}$ as soon as $C1\geq 1$.

It follows that, for some positive $c_{2}$, with probability greater than $1-\exp(-nc_{2}\gamma_{v}^{2})$

[TABLE]

• •

$\|\widehat{f}_{v}-{f}_{v}\|_{2}\geq\gamma_{v}$. We apply Lemma 8.7 (page 8.7) to the function $g_{v}=\widehat{f}_{v}-{f}_{v}$ with $b=2$. It follows that, for some positive $c_{2}$, with probability greater than $1-\exp(-nc_{2}\gamma_{v}^{2})$,

[TABLE]

$\Box$

8.5.4. Proof of Lemma

8.4 (page 8.4):

Let $d(f)$ be defined by (36), and let $\mathcal{G}(f)$ and $\mathcal{G}^{\prime}(f)$ be the following sets

[TABLE]

Let us consider the two events $\mathcal{B}$ and $\mathcal{B}^{\prime}$ defined as follows:

[TABLE]

Let us first remark that $\mathcal{B}^{\prime}$ is included into $\mathcal{B}$: if $h\in\mathcal{B}^{\prime}$ , then $h\in\mathcal{G}$, $\|h\|_{2}=d(f)$ and $\|h\|_{n}^{2}\geq d^{(}f)2/2$. It follows that $\|h\|_{n}^{2}\geq\|h\|_{2}^{2}/2$ and $\|h\|_{2}\geq d(f)$. Therefore Lemma 8.4 is proved if $\mathcal{B^{\prime}}$ holds with high probability. Consider

[TABLE]

We show that the event $Z_{n}(\mathcal{G}^{\prime})\leq d(f)^{2}/2$ has probability greater than $1-c_{1}\exp(-nc_{2}d(f)^{2})$.

Let us briefly recall the notion of covering numbers for a totally bounded metric space $(\mathbb{G},\rho)$, consisting of a set $\mathbb{G}$ and a metric $\rho$ defined from $\mathbb{G}\times\mathbb{G}$ into $\mathbb{R}_{+}$. A $\delta$-covering set of $\mathbb{G}$ is a collection of functions $f^{1},\cdots,f^{N}$ such that for all $f\in\mathbb{G}$ there exists $k\in\{1,2,\cdots,N\}$ such that $\rho(f,f^{k})\leq\delta$.

The $\delta$-covering number $N(\delta,\mathbb{G},\rho)$ is the cardinality of the smallest $\delta$-covering set. A propercovering restricts the covering to use only elements in the set $\mathbb{G}$. The proper covering number denoted $N_{\rm{pr}}(\delta,\mathbb{G},\rho)$ satisfies

[TABLE]

Let us now consider a $d(f)/8$-covering of $(\mathcal{G}^{\prime},\|\cdot\|_{n})$, so that, for all $g$ in $\mathcal{G}^{\prime}$ there exists $g^{k}$ such that $\|g-g^{k}\|_{n}\leq d(f)/8.$ The associated proper covering number is

[TABLE]

Now, for all $g\in\mathcal{G}^{\prime}$, $T_{1}=\|g^{k}\|_{n}^{2}-\|g\|_{n}^{2},$ and $T_{2}=d^{2}(f)-\|g^{k}\|_{n}^{2}$, we write

[TABLE]

The proof is splitted into four steps:

1. Step 1

The first step consists in showing that

[TABLE] 2. Step 2

The second step consists in proving that, for $N_{\rm{pr}}$ given at Equation (57) and for some constant $C$,

[TABLE] 3. Step 3

The third step concerns the control of $N_{\rm{pr}}$: we show the following result

[TABLE] 4. Step 4

The last step consists in bounding from above the Gaussian complexity:

[TABLE]

Let us conclude the proof of the lemma before proving these four steps.

Putting together Steps 3 and 4, for $c_{3}<C$ and $C_{1}$ large enough, then Step 2 states that

[TABLE]

Now, we have

[TABLE]

We conclude the proof of the lemma, by noting that $d(f)^{2}\geq C_{1}^{2}\sum_{v}\lambda_{n,v}^{2}$ (see (13), (36) and (14)).

Proof of Step 1: We start by writing that

[TABLE]

Using that $(a+b)^{2}\leq 2a^{2}+2b^{2}$, $g\in\mathcal{G}^{\prime}$, and $g$ satisfies Condition C3, we get

[TABLE]

Besides, the covering set is constructed such that $\|g^{k}-g\|_{n}\leq d(f)/8$. It follows that Step 1 is proved.

Proof of Step 2: We prove that for some constant $C$,

[TABLE]

As $g^{k}\in\mathcal{G}^{\prime}$, $d=\|g^{k}\|_{2}$. Then

[TABLE]

Applying Theorem 3.5 in Chung and Lu [10] we have that for all positive $\lambda$

[TABLE]

Taking $\lambda=nd(f)^{2}/4$ and using that $\|g^{k}\|_{2}^{2}=d(f)^{2}$ we get

[TABLE]

It follows that

[TABLE]

It remains to calculate $E_{\mathbf{X}}g^{4}(\mathbf{X})$ for $g\in\mathcal{G}^{\prime}$. Precisely we show the following result:

[TABLE]

This result comes from the property of the RKHS $\mathcal{H}$: indeed $g\in\mathcal{H}$ is written $g=\sum_{v\in\mathcal{P}}g_{v}$ where the functions $g_{v}$ are centered and orthogonal in $\mathbb{L}^{2}(P_{\mathbf{X}})$. Therefore $\mathbb{E}g^{4}(\mathbf{X})$ is the sum of the following terms:

[TABLE]

Using the Cauchy-Schwartz inequality and the fact that $\|g_{v}\|_{\infty}\leq\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$, and $\|g\|_{2}=d(f)$ (because $g\in\mathcal{G}^{\prime}$), we get that $A_{1}$ is proportionnal to $d(f)^{2}$, $A_{2},A_{3},A_{5}$ to $d(f)^{4}$, and $A_{4}$ to $d(f)^{3}$. For example

[TABLE]

After calculation of the terms $A_{i}$, since $d(f)^{2}$ is assumed to be smaller than one, we get that

[TABLE]

Step 2 is proved by combining (59) and (60).

Proof of Step 3: Let $N_{\rm{pr}}$ be defined at Equation (57). We prove that

[TABLE]

We start from (56) and write that

[TABLE]

Using the Sudakov minoration (see Pisier [29]) we have that for all positive $\omega$

[TABLE]

Hence by taking $\omega=d(f)/16$, Step 3 is proved.

Proof of Step 4: The last step consists in bounding from above the Gaussian complexity $\mathbb{E}_{\varepsilon}\sup_{g\in\mathcal{G}^{\prime}}|V_{n,\varepsilon}(g)|$. This control is performed by using Lemma 8.5 (page 8.5). According to Inequality (51),

[TABLE]

with $\lambda_{n,v}$ defined by Equation (13) satisfying $C_{1}\lambda_{n,v}\leq\gamma_{v}$ and $C_{1}\lambda_{n,v}^{2}\leq\mu_{v}$ for all $v\in\mathcal{P}$. It follows

[TABLE]

and according to Condition C1,

[TABLE]

because $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$. Now, according to Condition C2, and using that $2ab\leq a^{2}+b^{2}$, we get

[TABLE]

the last inequality coming from the fact for all $g\in\mathcal{G}^{\prime}$, $\|g\|^{2}_{2}=d(f)^{2}\geq\sum_{v\in S}\|g_{v}\|_{2}^{2}$.

Finally, thanks to (36), we get

[TABLE]

$\Box$

8.6. Proofs of intermediate Lemmas.

8.6.1. Proof of Lemma 8.5 (page 8.5):

Let us write that the kernel $k_{v}$ is written as :

[TABLE]

where $(\phi_{v,k})_{k=1}^{\infty}$ is an orthonormal basis of $\mathbb{L}^{2}(P_{v})$, where $P_{v}=\prod_{a\in v}P_{a}$.

Let us consider the class of functions $\mathcal{K}(t)$ defined as

[TABLE]

It comes that

[TABLE]

In the following, we set $\mu_{k,v}(t)=\min\left\{t^{2},\omega_{k,v}\right\}$. Hence

[TABLE]

as soon as $g_{v}\in\mathcal{K}(t)$.

Now, let us prove the lemma:

[TABLE]

The last inequality follows from the Cauchy-Schwartz inequality and Inequality (61). Now, simple calculation leads to

[TABLE]

$\Box$

8.6.2. Proof of Lemma

8.6 (page 8.6):

Using that $\left|\sqrt{a}-\sqrt{b}\right|\leq\sqrt{|a-b|},$ we get

[TABLE]

Hence

[TABLE]

The centered process

[TABLE]

satisfies a concentration inequality given, for example, by Theorem 2.1 in Bartlett et al. [3] : if $\mathcal{C}$ is a class of functions $f$ such that $\|f\|_{\infty}\leq B$ and $Ef(\mathbf{X})=0$, and if there exists $\gamma>0$ such that for every $f\in\mathcal{C}$, $\mathop{\rm Var}\nolimits f(\mathbf{X})\leq\gamma^{2}$. Then for every $x>0$, with probability at least $1-e^{-x}$,

[TABLE]

For any $t>0$, for $\mathcal{G}(t)$ defined by (45), let us consider the class of functions $\mathcal{C}(t)$ defined as follows

[TABLE]

Note that if $f\in\mathcal{C}(t)$, $E_{\mathbf{X}}f(\mathbf{X}_{v})=0$ and $\|f\|_{\infty}\leq b^{2}$. We have to study

[TABLE]

It is easy to see that

[TABLE]

Let $\zeta_{i}$ be independent and identically random variables Rademacher distributed and let $E_{\mathbf{X},\zeta}$ denotes the expectation with respect to the law of $(\mathbf{X},\zeta)$. By a symmetrization argument,

[TABLE]

Since $\|g_{v}\|_{\infty}\leq b$, applying the contraction principal (see Ledoux-Talagrand [23]) we get that, for $Q_{n,v}(t)$ defined by (11),

[TABLE]

The last inequality was proved by Mendelson [28], Theorem 41 (see the proof of Lemma 8.5). Now, thanks to (62) we get that for all $x>0$, with probability greater than $1-e^{-x}$

[TABLE]

Taking $x=c_{2}nt^{2}$, $t\geq\nu_{n}$, we have that with probability greater than $1-e^{-c_{2}nt^{2}}$

[TABLE]

The infimum of the right hand side is reached in $\alpha=\sqrt{c_{2}b/16\Delta}$, and equals

[TABLE]

The constants $\Delta$ and $c_{2}$ should satisfy that this infimum is strictly smaller than $b^{2}/4$. For example, if $16\Delta<b/8$, it remains to choose $c_{2}$ small enough such that

[TABLE]

$\Box$

8.6.3. Proof of Lemma 8.7 (page 8.7):

Let $t>\nu_{n,v}$ and $h$ be defined as $h=tg_{v}/\|g_{v}\|_{2}$. If $g_{v}$ satisfies the assumptions of the lemma, then $h$ satisfies $\|h\|_{2}=t$, $\|h\|_{\mathcal{H}}\leq 2$ and $\|h\|_{\infty}\leq b$. Applying Lemma 8.6 (page 8.6) to the function $h$, we obtain that for all $t\geq\nu_{n,v}$, with probability greater than $1-\exp(-c_{2}nt^{2})$, we have $|t-\|h\|_{n}|\leq bt/2$ for all $h\in\mathcal{G}(t)$. This concludes the proof of the lemma.

8.6.4. Proof of Lemma

8.8 (page 8.8):

The proof of Lemma 8.8 is based on an isoperimetric inequality for Gaussian processes (Borell [6] or Cirel’son et al. [11]) as it is stated in Theorem (3.8), page 61 in Massart [26]. Let us recall this inequality:

** Lemma 8.11****.**

Let $P$ be the Gaussian probability measure on $\mathbb{R}^{n}$ and let $\phi$ be a function from $\mathbb{R}^{n}$ to $\mathbb{R}$, and $\|\phi\|_{L}$ its Lipschitz semi-norm:

[TABLE]

Let $\overline{\Phi}$ be the cumulative distribution of the standard Gaussian distribution. Then for any $u$,

[TABLE]

We apply Lemma 8.11 to $\phi(\varepsilon_{1},\ldots,\varepsilon_{n})=W_{n,n,v}(t)$. By Cauchy-Schwarz Inequality, $\|\phi\|_{L}=t/\sqrt{n}$. It follows that Lemma 8.8 is proved since

[TABLE]

$\Box$

8.6.5. Proof of Lemma 8.9 (page 8.9):

We start with the proof of (48) in Lemma 8.9 by applying once again Lemma 8.11 given above, to the function $\phi(\boldsymbol{\varepsilon})=\phi(\varepsilon_{1},\ldots,\varepsilon_{n})=W_{n,2,v}(t)$. On the event $\Omega_{v,t}$ defined by (46), we have $\|g_{v}\|_{n}\leq bt/2+\|g_{v}\|_{2}$. Besides if $\|g_{v}\|_{\mathcal{H}_{v}}\leq 2$, then $\|g_{v}\|_{\infty}\leq 2$. Therefore applying Lemma 8.6 with $b=2$, we get that if $\|g_{v}\|_{2}\leq t$,

[TABLE]

leading to $\|\phi\|_{L}=2t/\sqrt{n}$. It follows that (48) in Lemma 8.9 is proved since

[TABLE]

We now come to the proof of (49) in Lemma 8.9 using a Poissonian inequality for self-bounded processes (see Boucheron et al. [7]) and Theorem 5.6, p 158 in Massart [26]). Let us recall it in the particular case we are interested in:

Theorem 8.1**.**

Let $X_{1},\cdots,X_{n}$ be $n$ iid variables. For $i\in\{1,\cdots,n\}$ let $X_{(-i)}=(X_{1},\ldots,X_{i-1},X_{i+1},\ldots,X_{n})$. Let $Z$ be a nonnegative and bounded measurable function of $X=(X_{1},\cdots,X_{n})$. Assume that for all $i\in\{1,\cdots,n\}$, there exists a measurable function $Z_{i}$ of $X_{(-i)}$ such that $0<Z-Z_{i}\leq 1,$ and $\sum_{i=1}^{n}(Z-Z_{i})\leq Z$. Then, for all $x>0$, we have $P\{Z\geq E(Z)+x\}\leq\exp\left(-x^{2}/2E(Z)\right)$.

We apply this result to $Z$ defined as

[TABLE]

The variable $Z$ is positive, and because the distribution of $(\varepsilon_{1},\ldots,\varepsilon_{n})$ is symmetric, we have that

[TABLE]

Let $\tau$ be the function in $\mathcal{H}_{v}$ such that $Z=E_{\boldsymbol{\varepsilon}}nV_{n,\varepsilon}(\tau)$ (note that $\tau$ depends on $(\mathbf{X}_{1},\ldots,\mathbf{X}_{n})$ and on $(\varepsilon_{1},\ldots,\varepsilon_{n})$), and let

[TABLE]

We show that $Z$ and $Z_{i}$ satisfy the assumptions of Theorem 8.1:

[TABLE]

where the last inequality comes from the fact that $\sup_{x\in\mathcal{X}}|\tau(\mathbf{X})|\leq\|\tau\|_{\mathcal{H}_{v}}\leq 2$. Moreover $Z-Z_{i}\geq 0$ since

[TABLE]

Finally we have:

[TABLE]

Therefore, following Theorem 8.1, we get that for all postive $u$

[TABLE]

As $E_{\mathbf{X},\boldsymbol{\varepsilon}}W_{n,2,v}(t)\leq Q_{n,v}(t)$, see Lemma 8.5 page 8.5, we get the expected result since for all positive $x$

[TABLE]

$\Box$

8.6.6. Proof of Lemma 8.10 (page 8.10):

From Lemma 8.8, page 8.8 with $t=\lambda_{n,v}=\delta$, we get that with probability greater than $1-4\exp(-n\lambda_{n,v}^{2}/2)$,

[TABLE]

The next step consists in comparing $W_{n,n,v}(\lambda_{n,v})$ and $W_{n,2,v}(2\lambda_{n,v})$. Recall that $\lambda_{n,v}\geq\nu_{n,v}$, see (13). Let $g_{v}$ such that $\|g_{v}\|_{n}\leq\lambda_{n,v}$.

• •

When $\|g_{v}\|_{2}\leq\lambda_{n,v}$, according to Lemma 8.6 (page 8.6), taking $b=2$ , since since $\|g_{v}\|_{n}\leq\lambda_{n,v}$, we get that with probability greater than $1-\exp(-c_{2}n\lambda_{n,v}^{2})$,

[TABLE]

• •

When $\|g_{v}\|_{2}\geq t$, we apply Lemma 8.7 (page 8.7) with $b=2$. For any function $g_{v}$ such that $\|g_{v}\|_{\infty}\leq 2$, and $\|g_{v}\|_{2}\geq\lambda_{n,v}$, we have $\|g_{v}\|_{2}\leq 2\|g_{v}\|_{n}\leq 2\lambda_{n,v}$.

This implies that, with probability greater than $1-\exp(-c_{2}n\lambda_{n,v}^{2})$ we have

[TABLE]

We now study the process $W_{n,2,v}(\lambda_{n,v})$. By applying (48) in Lemma 8.9, page 8.9, with $\delta=t=\lambda_{n,v}$ we get that with probability greater than $1-4\exp(-n\lambda_{n,v}^{2}/8)$

[TABLE]

It follows that

[TABLE]

Next, we apply (49) in Lemma 8.9, with $t=2\lambda_{n,v}$ and $x=4\lambda^{2}_{n,v}$. We get that

[TABLE]

with probability greater than

[TABLE]

The last inequality comes from the definition of $\nu_{n,v}$, see (12), and from the fact that $\lambda_{n,v}\geq\nu_{n,v}$, see (13).

Putting everything together, we get that with probability greater than $1-c_{1}\exp(-c_{2}n\lambda_{n,v}^{2})$ for some positive constants $c_{1},c_{2}$,

[TABLE]

Applying once again Lemma 8.8, page 8.8, we get that

[TABLE]

This ends the proof of the lemma by taking $\kappa=10+4\Delta$.

$\Box$

8.7. Algorithm:

Propositions 8.1-8.4

We consider the minimization of $C^{\prime}(f_{0},\boldsymbol{\theta})$ given at Equation (18). Because $C^{\prime}(f_{0},\boldsymbol{\theta})$ is convex and separable, we use a block coordinate descent algorithm, group $v$ by group $v$. We refer to Bubeck [9] for a review on convex optimization.

In what follows, the group $v$ is fixed, and for given values of $f_{0}$ and $\boldsymbol{\theta}_{w},w\neq v$, we look for the minimizer of $C^{\prime}$ with respect to $\boldsymbol{\theta}_{v}$: Setting

[TABLE]

we aim at minimizing with respect to $\boldsymbol{\theta}_{v}$

[TABLE]

If $\gamma^{\prime}_{v}=\mu^{\prime}_{v}=0$, then $\boldsymbol{\theta}_{v}=K_{v}^{-1}\mathbf{R}_{v}$ is the solution. In what follows we consider the case where at least one of both is non zero.

If $\partial Q_{v}$ denotes the subdifferential of $Q(\boldsymbol{\theta}_{v})$ with respect to $\boldsymbol{\theta}_{v}$, we need to solve $0\in\partial Q_{v}$. Let us recall that for all $v\in\mathcal{P}$ the matrices $K_{v}$ are symmetric and strictly definite positive.

Let us begin with the calculation of the subdifferential of $\|K_{v}\theta_{v}\|$ with respect to $\theta_{v}$: if $\theta_{v}\neq 0$, we have

[TABLE]

and if $\theta_{v}=0$ the subdifferential is the set of $x\in\mathbb{R}^{n}$ such that $\|K_{v}^{-1}x\|\leq 1$.

Therefore if $\theta_{v}\neq 0$,

[TABLE]

while if $\theta_{v}=0$,

[TABLE]

Let us begin with the case $\mu^{\prime}_{v}=0$.

Proposition 8.1**.**

Let $(\rho)_{+}$ denotes the positive part of $\rho\in\mathbb{R}$. If $\mu^{\prime}_{v}=0$,

[TABLE]

Proof of Proposition 8.1

The problem comes to minimize

[TABLE]

and to take $\boldsymbol{\theta}_{v}=K_{v}^{-1}\boldsymbol{\beta}_{v}$. The subdifferential of $U$ is given by

[TABLE]

We get

[TABLE]

$\Box$

Let $\mu^{\prime}_{v}>0$, the following proposition gives a necessary and sufficient condition for $\boldsymbol{\theta}_{v}=0$ to be the minimizer of $Q$. For the sake of clarity let us recall the definitions of $J$ and $J^{*}$ given in Section 5.1: For all $x\in\mathbb{R}^{n}$

[TABLE]

Proposition 8.2**.**

The minimizer of $Q(\boldsymbol{\theta}_{v})$ is $\boldsymbol{\theta}_{v}=0$ if and only if $J^{*}\leq\gamma^{\prime 2}_{v}$.

Proof of Proposition 8.2

1. (1)

Let us assume that $\boldsymbol{\theta}_{v}=0$ is the minimizer of $Q(\boldsymbol{\theta}_{v})$. Then $0\in\partial Q$ implies that there exists $t^{*}$ and $s^{*}$ such that $\|K_{v}^{-1}t^{*}\|\leq 1$ and $\|K_{v}^{-1/2}s^{*}\|\leq 1$ and such that $\gamma^{\prime}_{v}t^{*}+\mu^{\prime}_{v}s^{*}=2K_{v}\mathbf{R}_{v}$.

If $\gamma^{\prime}_{v}>0$, then

[TABLE]

Because $J^{*}\leq J(s^{*})$ and $\|K_{v}^{-1}t^{*}\|\leq 1$, we get that $J^{*}\leq\gamma^{\prime 2}_{v}$.

If $\gamma^{\prime}_{v}=0$, then $\mu^{\prime}_{v}s^{*}=2K_{v}\mathbf{R}_{v}$ and $J(s^{*})=J^{*}=0$. 2. (2)

Let us now assume that $J^{*}\leq\gamma^{\prime 2}_{v}$. Note that minimizing the convex function $J(s)$ over the convex set $\|K_{v}^{-1/2}s\|\leq 1$ has always a solution. Let us denote this solution by $s^{*}$.

If $\gamma^{{}^{\prime}}_{v}>0$, let $t^{*}=(2K_{v}\mathbf{R}_{v}-\mu^{\prime}_{v}s^{*})/\gamma^{{}^{\prime}}_{v}$. Then $-2K_{v}\mathbf{R}_{v}+\gamma^{\prime}_{v}t^{*}+\mu^{\prime}_{v}s^{*}=0$, and $-2K_{v}\mathbf{R}_{v}+\gamma^{\prime}_{v}t^{*}+\mu^{\prime}_{v}s^{*}\in\partial Q(0)$ since $\|K_{v}^{-1/2}s^{*}\|\leq 1$, and $\|K_{v}^{-1}t^{*}\|=J(s^{*})/\gamma^{{}^{\prime}}_{v}\leq 1$. Therefore $\boldsymbol{\theta}_{v}=0$ is the minimizer of $Q$.

If $\gamma^{\prime}_{v}=0$, and $J^{*}=0$, then $-2K_{v}\mathbf{R}_{v}+\mu^{\prime}_{v}s^{*}=0$, and $\|K_{v}^{-1/2}s^{*}\|\leq 1$. Therefore $\boldsymbol{\theta}_{v}=0$ is the minimizer of $Q$.

$\Box$

Proposition 8.3**.**

Let $\mu^{\prime}_{v}>0$ and $\boldsymbol{\theta}_{v}$ be the minimizer of $Q$.

1. (1)

If $\|2K_{v}^{1/2}\mathbf{R}_{v}\|\leq\mu^{\prime}_{v}$, then $\boldsymbol{\theta}_{v}=0$ 2. (2)

If not, for $\rho>0$, let

[TABLE]

and let $\rho^{*}$ defined as $\|\boldsymbol{\theta}(\rho^{*})\|=1$. Then $J(\theta(\rho^{*}))\leq\gamma^{\prime 2}_{v}$ if and only if $\boldsymbol{\theta}_{v}=0$,

Proof of Proposition 8.3

Minimizing $J(x)$ under the constraint $\|K_{v}^{-1/2}x\|\leq 1$ is equivalent to minimizing $K(\boldsymbol{\beta})=\|2\mathbf{R}_{v}-\mu^{\prime}_{v}K_{v}^{-1/2}\boldsymbol{\beta}\|^{2}$ under the constraint $\|\boldsymbol{\beta}\|^{2}\leq 1$. Let $\boldsymbol{\beta}^{*}=2K_{v}^{1/2}\mathbf{R}_{v}/\mu^{\prime}_{v}$. Then $K(\boldsymbol{\beta}^{*})=0$, which is smaller than $\gamma^{\prime 2}_{v}$, and if $\|\boldsymbol{\beta}^{*}\|\leq 1$, following Proposition 8.2, we get $\boldsymbol{\theta}_{v}=0$.

If $\|\boldsymbol{\beta}^{*}\|>1$, we have to solve a ridge regression problem by minimizing $K(\boldsymbol{\beta})+\rho\|\boldsymbol{\beta}\|^{2}$ for some positive $\rho$. The solution is given by $\boldsymbol{\theta}(\rho)$. Let us note that $\|\boldsymbol{\theta}(\rho)\|$ decreases to 0 when $\rho$ tends to infinity and that its maximum is $\|\boldsymbol{\theta}(0)\|=\|\boldsymbol{\beta}^{*}\|$. Therefore if $\|\boldsymbol{\beta}^{*}\|>1$, there exists $\rho^{*}$ such that $\|\boldsymbol{\theta}(\rho^{*})\|=1$. Following Proposition 8.2, Proposition 8.3 is proved. A numerical procedure can be used for calculating this $\rho^{*}$. $\Box$

Let us now consider the case where $\boldsymbol{\theta}_{v}$ is non zero. It should satisfy the subgradient condition $\partial Q_{v}=0$ which leads to

[TABLE]

which is equivalent to Equation (22).

Proposition 8.4**.**

For all $\rho_{1},\rho_{2}>0$ let

[TABLE]

If $\mu^{\prime}_{v}>0$, there exists a non zero solution to Equation (22) if and only if there exists $\rho_{1},\rho_{2}>0$ such that

[TABLE]

Then $\boldsymbol{\theta}_{v}=\boldsymbol{\theta}(\rho_{1},\rho_{2})$.

Proof of Proposition 8.4

If there exists a non zero solution to Equation (22) then $\|K_{v}\boldsymbol{\theta}_{v}\|$ and $\|K_{v}^{1/2}\boldsymbol{\theta}_{v}\|$ are non zero because $K_{v}$ is definite positive. Let $\rho_{1}=\gamma^{\prime}_{v}/2\|K_{v}\boldsymbol{\theta}_{v}\|$ and $\rho_{2}=\mu^{\prime}_{v}/2\|K^{1/2}_{v}\boldsymbol{\theta}_{v}\|$, then

[TABLE]

and, for such $\rho_{1},\rho_{2}$, Equation (64) is satisfied.

Conversely, if there exist $\rho_{1},\rho_{2}$ such that Equation (64) is satisfied, then necessarily $\|K_{v}\boldsymbol{\theta}(\rho_{1},\rho_{2)})\|$ and $\|K_{v}^{1/2}\boldsymbol{\theta}(\rho_{1},\rho_{2)})\|$ are non zero and $\rho_{1}=\gamma^{\prime}_{v}/2\|K_{v}\boldsymbol{\theta}(\rho_{1},\rho_{2})\|$ and $\rho_{2}=\mu^{\prime}_{v}/2\|K_{v}^{1/2}\boldsymbol{\theta}(\rho_{1},\rho_{2})\|$. Then

[TABLE]

which is exactly Equation (22) calculated in $\boldsymbol{\theta}_{v}=\boldsymbol{\theta}(\rho_{1},\rho_{2})$. $\Box$

Taking into account that $\|K_{v^{1/2}}\boldsymbol{\theta}_{v}\|\leq r_{v}$

As already mentionned in Section 5, one may want to minimize $C(f_{0},\boldsymbol{\theta})$ under the additional constraint that $\|K_{v^{1/2}}\boldsymbol{\theta}_{v}\|\leq r_{v}$ for some positive constant $r_{v},v\in\mathcal{P}$.

For each group $v$, we have thus to minimize $Q(\boldsymbol{\theta}_{v})$ under the constraint $\|K_{v^{1/2}}\boldsymbol{\theta}_{v}\|\leq r_{v}$. We know that this problem is equivalent to minimize

[TABLE]

for some $\lambda$ that depends on $r_{v}$.

Let us first remark that, for a fixed $\gamma^{{}^{\prime}}_{v}$, and $\lambda\geq 0$, if $\widehat{\boldsymbol{\theta}}_{v}(\mu^{{}^{\prime}}_{v}+\lambda)$ minimizes $Q(\boldsymbol{\theta}_{v})+\lambda\|K_{v}^{1/2}\boldsymbol{\theta}_{v}\|$ with respect to $\boldsymbol{\theta}_{v}$, then $\|K^{1/2}_{v}\widehat{\boldsymbol{\theta}}_{v}(\mu^{{}^{\prime}}_{v}+\lambda)\|\leq\|K^{1/2}_{v}\widehat{\boldsymbol{\theta}}_{v}(\mu^{{}^{\prime}}_{v})\|$. It can be easily proved by writing

[TABLE]

Therefore, one can proceed as follows: calculate $\widehat{\boldsymbol{\theta}}_{v}$ for $\lambda=0$. If $\|K^{1/2}_{v}\widehat{\boldsymbol{\theta}}_{v}\|\leq r_{v}$, then one go to the next step of the algorithm. If $\|K^{1/2}_{v}\widehat{\boldsymbol{\theta}}_{v}\|>r_{v}$, one increases $\lambda$ untill $\|K^{1/2}_{v}\widehat{\boldsymbol{\theta}}_{v}\|=r_{v}$.