RKHSMetaMod: An R package to estimate the Hoeffding decomposition of a complex model by solving RKHS ridge group sparse optimization problem

TL;DR

The paper introduces RKHSMetaMod, an R package that estimates the Hoeffding decomposition of complex models using RKHS ridge group sparse optimization, enabling sensitivity analysis and Sobol index estimation.

Contribution

It presents a novel R package implementing a penalized least-squares method for estimating Hoeffding decompositions in complex models, combining kernel methods with sparse optimization.

Findings

Efficient estimation of Hoeffding decomposition terms.

Ability to select non-zero Sobol indices.

Enhanced sensitivity analysis capabilities.

Abstract

In this paper, we propose an R package, called RKHSMetaMod, that implements a procedure for estimating a meta-model of a complex model. The meta-model approximates the Hoeffding decomposition of the complex model and allows us to perform sensitivity analysis on it. It belongs to a reproducing kernel Hilbert space that is constructed as a direct sum of Hilbert spaces. The estimator of the meta-model is the solution of a penalized empirical least-squares minimization with the sum of the Hilbert norm and the empirical L^2-norm. This procedure, called RKHS ridge group sparse, allows both to select and estimate the terms in the Hoeffding decomposition, and therefore, to select and estimate the Sobol indices that are non-zero. The RKHSMetaMod package provides an interface from R statistical computing environment to the C++ libraries Eigen and GSL. In order to speed up the execution time and optimize the storage memory, except for a function that is written in R, all of the functions of this package are written using the efficient C++ libraries through RcppEigen and RcppGSL packages. These functions are then interfaced in the R environment in order to propose a user-friendly package.

Figures

Tables

Table 1. Table 1: List of the reproducing kernels used to construct the RKHS ℋ ℋ \mathcal{H} . The range parameters r 𝑟 r in the Gaussian and Matérn 3 / 2 3 2 3/2 kernels are assumed to be fixed and set as 1 / 2 1 2 1/2 and 3 / 2 3 2 \sqrt{3}/2 , respectively. The value 1 1 1 is added to the Brownian kernel to relax the constraint of nullity at the origin.

Kernel type	Mathematics formula for $u\in {ℝ}^n,v\in ℝ$	RKHSMetaMod name
Gaussian	$k_a(u,v)=\exp \left(-{\Vert u-v\Vert }^2/2r^2\right)$	"gaussian"
Matérn 3/2	$k_a(u,v)=\left(1+\sqrt{3}\left|u-v\right|/r\right)\exp \left(-\sqrt{3}\left|u-v\right|/r\right)$	"matern"
Brownian	$k_a(u,v)=\min (u,v)+1$	"brownian"
Quadratic	$k_a(u,v)={\left(u^{T}v+1\right)}^2$	"quad"
Linear	$k_a(u,v)=u^{T}v+1$	"linear"

Table 2. Table 2: Example 1 : The columns of the table correspond to the different datasets with n ∈ { 50 , 100 , 200 } 𝑛 50 100 200 n\in\{50,100,200\} and d = 5 𝑑 5 d=5 . Each line of the table, from up to down, gives the value of GPE obtained for each dataset associated with the "matern", "brownian" and "gaussian" kernels, respectively.

$n$	$50$	$100$	$200$
$GP{E}_m$	$0.13$	$0.07$	$0.03$
$GP{E}_b$	$0.14$	$0.10$	$0.05$
$GP{E}_g$	$0.15$	$0.10$	$0.07$

Table 3. Table 3: Example 1 : The columns of the table correspond to the different datasets with n ∈ { 50 , 100 , 200 } 𝑛 50 100 200 n\in\{50,100,200\} and d = 5 𝑑 5 d=5 . Each line of the table, from up to down, gives the value of MSE obtained for each dataset associated with the "matern", "brownian" and "gaussian" kernels, respectively.

$n$	$50$	$100$	$200$
$MS{E}_m$	$75.12$	$46.72$	$28.22$
$MS{E}_b$	$110.71$	$84.99$	$41.06$
$MS{E}_g$	$78.22$	$94.67$	$67.02$

Table 4. Table 4: Example 1 : The first line of the table gives the true values of the Sobol indices × 100 absent 100 \times 100 . The second line gives the mean of the estimated empirical Sobol indices ( S ^ v , . = ∑ r = 1 N r S ^ v , r / N r \widehat{S}_{v,.}=\sum_{r=1}^{N_{r}}\widehat{S}_{v,r}/N_{r} ) × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} calculated over fifty simulations for n = 200 𝑛 200 n=200 and "matern" kernel. The sum of the Sobol indices is displayed in the last column.

$v$	$\left\{1\right\}$	$\left\{2\right\}$	$\left\{3\right\}$	$\{1,2\}$	$\{1,3\}$	$\{2,3\}$	$\{1,2,3\}$	sum
$S_v$	43.30	24.30	19.20	5.63	4.45	2.50	0.57	99.98
${\hat{S}}_{v,.}$	46.10	26.33	20.62	2.99	2.22	1.13	0.0	99.39

Table 5. Table 5: Example 2 : Obtained prediction errors in step 1.

${\mu }_g$	$1.304$	$0.652$	$0.326$	$0.163$	$0.081$	$0.041$	$0.020$	$0.010$	$0.005$
$\gamma =0$	$0.197$	$0.156$	$0.145$	$0.097$	$0.063$	$0.055$	$0.056$	$0.063$	$0.073$

Table 6. Table 6: Example 2 : Obtained prediction errors in step 2.

$\mu$	$0.081/\sqrt{n}$	$0.041/\sqrt{n}$	$0.020/\sqrt{n}$
$\gamma =0.2$	$0.153$	$0.131$	$0.119$
$\gamma =0.1$	$0.098$	$0.079$	$0.072$
$\gamma =0.01$	$0.065$	$0.054$	$0.053$
$\gamma =0.005$	$0.064$	$0.054$	$0.054$

Table 7. Table 7: Example 2 : The estimated empirical Sobol indices × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} . The last two columns show ∑ v S ^ v subscript 𝑣 subscript ^ 𝑆 𝑣 \sum_{v}\widehat{S}_{v} and RE, respectively.

$v$	$\left\{1\right\}$	$\left\{2\right\}$	$\left\{3\right\}$	$\{1,2\}$	$\{1,3\}$	$\{2,3\}$	$\{1,2,3\}$	sum	RE
${\hat{S}}_v$	$42.91$	$25.50$	$20.81$	$4.40$	$3.84$	$2.13$	$0.00$	$99.60$	$1.64$

Table 8. Table 8: Example 3 : The kernel used is "matern". The execution time for the functions RKHSgrplasso and pen _ _ \_ MetMod is displayed in each row for two pairs of values of tuning parameters ( μ 1 = μ m ​ a ​ x / ( n × 2 7 ) , γ = 0.01 ) formulae-sequence subscript 𝜇 1 subscript 𝜇 𝑚 𝑎 𝑥 𝑛 superscript 2 7 𝛾 0.01 (\mu_{1}=\mu_{max}/(\sqrt{n}\times 2^{7}),\gamma=0.01) on up, and ( μ 2 = μ m ​ a ​ x / ( n × 2 8 ) , γ = 0.01 ) formulae-sequence subscript 𝜇 2 subscript 𝜇 𝑚 𝑎 𝑥 𝑛 superscript 2 8 𝛾 0.01 (\mu_{2}=\mu_{max}/(\sqrt{n}\times 2^{8}),\gamma=0.01) on below. In the column | S f ^ | subscript 𝑆 ^ 𝑓 |S_{\widehat{f}}| , the number of the active groups associated with each estimated RKHS meta-model is displayed.

$(n,d)$	calc$\_$Kv	mu$\_$max	RKHSgrplasso	pen$\_$MetMod	$\left|S_{\hat{f}}\right|$	sum
(100,5)	0.09s	0.01s	1s	2s	18	$\sim$ 3s
			2s	3s	19	$\sim$ 5s
(500,10)	33s	9s	247s	333s	39	$\sim$ 10min
			599s	816s	64	$\sim$ 24min
(1000,10)	197s	53s	959s	1336s	24	$\sim$ 42min
			2757s	4345s	69	$\sim$ 2h
(2000,10)	1498s	420s	3984s	4664s	12	$\sim$ 2h:56min
			12951s	22385s	30	$\sim$ 10h:20min
(5000,10)	34282s	6684s	38957s	49987s	11	$\sim$ 36h:05min
			99221s	111376s	15	$\sim$ 69h:52min

Table 9. Table 9: Example 3 : Obtained prediction errors.

$n$	$({\mu }_{max}/\left(\sqrt{n}\times 2^7\right),\gamma )$	$({\mu }_{max}/(\sqrt{n}\times 2^8,\gamma )$
$2000$	$0.052$	$0.049$
$5000$	$0.049$	$0.047$

Table 10. Table 10: Example 3 : The estimated empirical Sobol indices × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} associated with each estimated RKHS meta-model is printed. The last two columns show ∑ v S ^ v subscript 𝑣 subscript ^ 𝑆 𝑣 \sum_{v}\widehat{S}_{v} and RE, respectively. We have μ 1 = μ m ​ a ​ x / ( n × 2 7 ) subscript 𝜇 1 subscript 𝜇 𝑚 𝑎 𝑥 𝑛 superscript 2 7 \mu_{1}=\mu_{max}/(\sqrt{n}\times 2^{7}) , μ 2 = μ m ​ a ​ x / ( n × 2 8 ) subscript 𝜇 2 subscript 𝜇 𝑚 𝑎 𝑥 𝑛 superscript 2 8 \mu_{2}=\mu_{max}/(\sqrt{n}\times 2^{8}) and γ = 0.01 𝛾 0.01 \gamma=0.01 .

$n$	$v$	$\left\{1\right\}$	$\left\{2\right\}$	$\left\{3\right\}$	$\{1,2\}$	$\{1,3\}$	$\{2,3\}$	$\{1,2,3\}$	sum	RE
$2000$	${\hat{S}}_{v;({\mu }_1,\gamma )}$	$45.54$	$24.78$	$21.01$	$3.96$	$3.03$	$1.65$	$0.00$	$99.97$	$2.12$
	${\hat{S}}_{v;({\mu }_2,\gamma )}$	$45.38$	$25.07$	$19.69$	$4.36$	$3.66$	$1.79$	$0.00$	$99.95$	$1.79$
$5000$	${\hat{S}}_{v;({\mu }_1,\gamma )}$	$44.77$	$25.39$	$20.05$	$4.49$	$3.38$	$1.90$	$0.00$	$99.98$	$1.81$
	${\hat{S}}_{v;({\mu }_2,\gamma )}$	$43.78$	$24.99$	$19.56$	$5.43$	$3.90$	$2.32$	$0.00$	$99.98$	$1.29$

Table 11. Table 11: Example 5 : The true values of the Sobol indices × 100 absent 100 \times 100 greater than 10 − 2 superscript 10 2 10^{-2} are given in the first raw. The estimated Sobol indices associated with the RKHS meta-model ( S ^ v subscript ^ 𝑆 𝑣 \hat{S}_{v} ) and the Kriging based meta-model ( S ^ k ​ m v subscript ^ 𝑆 𝑘 subscript 𝑚 𝑣 \widehat{S}_{km_{v}} ) are given in second and third rows, respectively.

$v$	$\left\{1\right\}$	$\left\{2\right\}$	$\left\{3\right\}$	$\left\{4\right\}$	$\{1,2\}$	$\{1,3\}$	$\{1,4\}$	$\{2,3\}$	$\{2,4\}$	$\{1,2,3\}$	$\{1,2,4\}$	sum
$S_v$	71.62	17.90	2.37	0.72	5.97	0.79	0.24	0.20	0.06	0.07	0.02	99.96
${\hat{S}}_v$	75.78	17.42	1.71	0.47	4.00	0.05	0.07	0.28	0.09	0.00	0.00	99.87
${\hat{S}}_{k{m}_v}$	71.18	15.16	1.42	0.44	0.00	0.00	0.00	0.00	0.00	0.00	0.00	88.20

Table 12. Table 12: List of the input arguments of the RKHSMetMod function.

Input parameter	Description
Y	Vector of the response observations of size $n$.
X	Matrix of the input observations with $n$ rows and $d$ columns. Rows correspond to the observations and columns correspond to the variables.
kernel	Character, indicates the type of the kernel (see Table 1) chosen to construct the RKHS $\mathscr{H}$.
Dmax	Integer, between $1$ and $d$, indicates the maximum order of interactions considered in the RKHS meta-model: Dmax$=1$ is used to consider only the main effects, Dmax$=2$ to include the main effects and the second-order interactions, and so on.
gamma	Vector of non-negative scalars, values of the tuning parameter $\gamma$ in decreasing order. If $\gamma =0$ the function solves the RKHS group lasso optimization problem and for $\gamma >0$ it solves the RKHS ridge group sparse optimization problem.
frc	Vector of positive scalars. Each element of the vector sets a value to the tuning parameter $\mu$: $\mu ={\mu }_{\text{max}}/\left(\sqrt{n}\times \text{frc}\right)$. The value ${\mu }_{\text{max}}$ (see Equation (16)) is calculated inside the program.
verbose	Logical. Set as TRUE to print: the group $v$ for which the correction of the Gram matrix $K_v$ is done (see Section 3.1), and for each pair of the tuning parameters $(\mu ,\gamma )$: the number of current iteration, active groups and convergence criteria. It is set as FALSE by default.

Table 13. Table 13: List of the arguments of the output "Meta-Model" of RKHSMetMod function.

Output parameter	Description
intercept	Scalar, estimated value of intercept.
teta	Matrix with vMax rows and $n$ columns. Each row of the matrix is the estimated vector ${\theta }_v$ for $v=1,\dots ,$vMax.
fit.v	Matrix with $n$ rows and vMax columns. Each row of the matrix is the estimated value of $f_v=K_v{\theta }_v$.
fitted	Vector of size $n$, indicates the estimator of $m$.
Norm.n	Vector of size vMax, estimated values for the empirical $L^2$-norm.
Norm.H	Vector of size vMax, estimated values for the Hilbert norm.
supp	Vector of active groups.
Nsupp	Vector of the names of the active groups.
SCR	Scalar equals to ${\Vert Y-f_0-{\sum }_v{K}_v{\theta }_v\Vert }^2$.
crit	Scalar indicates the value of the penalized criterion.
gamma.v	Vector of size vMax, coefficients of the empirical $L^2$-norm.
mu.v	Vector of size vMax, coefficients of the Hilbert norm.
iter	List of two components: maxIter, and the number of iterations until the convergence is achieved.
convergence	TRUE or FALSE. Indicates whether the algorithm has converged or not.
RelDiffCrit	Scalar, value of the first convergence criterion at the last iteration, i.e. ${\Vert {\theta }_{lastIter}-{\theta }_{lastIter-1}/{\theta }_{lastIter-1}\Vert }^2$.
RelDiffPar	Scalar, value of the second convergence criterion at the last iteration, i.e. $cri{t}_{lastIter}-cri{t}_{lastIter-1}/cri{t}_{lastIter-1}$.

Equations

$$Y=m(X)+\epsilon ,$$

$$m(X)=m_0+\sum _{a=1}^d{m}_a(X_a)+\sum _{a<a^{\prime }}{m}_{a,a^{\prime }}(X_a,X_{a^{\prime }})+...+m_{1,...,d}(X).$$

$$m_0$$

$$m_{a,a^{\prime }}(X_a,X_{a^{\prime }})$$

$$m(X)=m_0+\sum _{v\in \mathcal{P}}{m}_v(X_v).$$

$$\text{var}(m(X))=\sum _{v\in \mathcal{P}}\text{var}(m_v(X_v)).$$

$$S_v=\frac{\text{var}(m_v(X_v))}{\text{var}(m(X))}.$$

$$m(X)=\sum _{j=0}^{\infty }{h}_j{\varphi }_j(X),$$

$$f(X)=<f,k(X,.){>}_{\mathcal{H}}=f_0+\sum _{v\in \mathcal{P}}{f}_v(X_v),$$

$$\gamma \Vert f{\Vert }_{n,1}+\mu \Vert f{\Vert }_{\mathcal{H},1},$$

$$\gamma \Vert f{\Vert }_n+\mu \Vert f{\Vert }_{\mathcal{H}},\text{mbox}with\Vert f{\Vert }_n=\sum _{v\in \mathcal{P}}\Vert f_v{\Vert }_n,\Vert f{\Vert }_{\mathcal{H}}=\sum _{v\in \mathcal{P}}\Vert f_v{\Vert }_{{\mathcal{H}}_v}.$$

$$\displaystyle\mathcal{L}(f)=\frac{1}{n}\sum_{i=1}^{n}\Big{(}Y_{i}-f_{0}-\sum_{v\in\mathcal{P}}f_{v}(X_{vi})\Big{)}^{2}+\sum_{v\in\mathcal{P}}\gamma_{v}\|f_{v}\|_{n}+\sum_{v\in\mathcal{P}}\mu_{v}\|f_{v}\|_{\mathcal{H}_{v}},$$

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

$$f_v(.)=\sum _{i=1}^n{\theta }_{vi}{k}_v(X_{vi},.),\text{mbox}and\Vert f_v{\Vert }_{{\mathcal{H}}_v}^2=\sum _{i,i^{\prime }=1}^n{\theta }_{vi}{\theta }_{v{i}^{\prime }}{k}_v(X_{vi},X_{v{i}^{\prime }}).$$

$$C(f_0,\theta )=\Vert Y-f_0{I}_n-\sum _{v\in \mathcal{P}}{K}_v{\theta }_v{\Vert }^2+\sqrt{n}\sum _{v\in \mathcal{P}}{\gamma }_v\Vert K_v{\theta }_v\Vert +n\sum _{v\in \mathcal{P}}{\mu }_v\Vert K_v^{1/2}{\theta }_v\Vert .$$

$$\hat{f}(X)={\hat{f}}_0+\sum _{v\in \mathcal{P}}{\hat{f}}_v(X_v)\text{mbox}with{\hat{f}}_v(X_v)=\sum _{i=1}^n{\hat{\theta }}_{vi}{k}_v(X_{vi},X_v).$$

$$C(f_0,\theta )=\Vert Y-f_0{I}_n-\sum _{v\in \mathcal{P}}{K}_v{\theta }_v{\Vert }^2+\sqrt{n}\gamma \sum _{v\in \mathcal{P}}\Vert K_v{\theta }_v\Vert +n\mu \sum _{v\in \mathcal{P}}\Vert K_v^{1/2}{\theta }_v\Vert .$$

$$C_g(f_0,\theta )=\Vert Y-f_0{I}_n-\sum _{v\in \mathcal{P}}{K}_v{\theta }_v{\Vert }^2+n\mu \sum _{v\in \mathcal{P}}\Vert K_v^{1/2}{\theta }_v\Vert ,$$

$${\mu }_g=\sqrt{n}\mu .$$

$$\displaystyle\mu_{\text{max}}=\max_{v}\Big{(}2\|K_{v}^{1/2}(Y-\bar{Y})\|\Big{)}/\sqrt{n}.$$

$$\text{ErrPred}(\mu ,\gamma )=\sum _{i=1}^{n^{\text{test}}}(Y_i^{\text{test}}-{\hat{f}}_{(\mu ,\gamma )}(X_i^{\text{test}}){)}^2/n^{\text{test}},$$

$${\hat{f}}_{(\mu ,\gamma )}(X^{\text{test}})={\hat{f}}_0+\sum _{v\in S_{\hat{f}}}\sum _{i=1}^n{\hat{\theta }}_{vi}{k}_v(X_{vi},X_v^{\text{test}}).$$

$$\text{var}(\hat{f}(X))=\sum _{v\in \mathcal{P}}\text{var}({\hat{f}}_v(X_v)),\text{mbox}where\forall v\in \mathcal{P},\text{var}({\hat{f}}_v(X_v))=E_X({{\hat{f}}_v}^2(X_v))=\Vert {\hat{f}}_v{\Vert }_2^2.$$

$$\widehat{\text{var}}({\hat{f}}_v(X_v))=\frac{1}{n-1}\sum _{i=1}^n({\hat{f}}_v(X_{vi})-{\hat{f}}_{v.}{)}^2.$$

$${\hat{S}}_v=\widehat{\text{var}}({\hat{f}}_v(X_v))/\sum _{v\in \mathcal{P}}\widehat{\text{var}}({\hat{f}}_v(X_v)).$$

$$k_{0a}=\min (X_a,X_a^{\prime })+1-(3/4)(1+X_a-X_a^2/2)(1+X_a^{\prime }-X_a^{\prime 2}/2).$$

$$k_{0a}$$

$$m(X)=\prod _{a=1}^d\frac{|4X_a-2|+c_a}{1+c_a},c_a>0.$$

$$\forall v\in \mathcal{P},S_v=\frac{1}{D}\prod _{a\in v}{D}_a,D_a=\frac{1}{3(1+c_a{)}^2},D=\prod _{a=1}^d(D_a+1)-1.$$

$${\mu }_{(1:5)}={\mu }_{max}/(\sqrt{n}\times 2^{(2:6)}),{\gamma }_{(1:5)}=(0.2,0.1,0.01,0.005,0).$$

$$GPE=\frac{1}{N_r}\sum _{r=1}^{N_r}\frac{1}{N}\sum _{i=1}^N({\hat{f}}_r(X_i^{third})-m(X_i^{third}){)}^2.$$

Full text

Keywords: Meta model, Hoeffding decomposition, Ridge Group Sparse penalty, Reproducing Kernel Hilbert Spaces.

1 Introduction

Consider a phenomenon described by a model $m$ depending on $d$ input variables $X=(X_{1},...,X_{d})$. This model $m$ from $\mathbb{R}^{d}$ to $\mathbb{R}$ may be a known model that is calculable in all points of $X$, i.e. $Y=m(X)$, or it may be an unknown regression model defined as follows:

[TABLE]

where the error $\varepsilon$ is assumed to be centered with a finite variance, i.e. $E(\varepsilon)=0$ and $\mbox{var}(\varepsilon)<\infty$. The components of $X$ are independent with a known law $P_{X}=\prod_{a=1}^{d}P_{X_{a}}$ on $\mathcal{X}$, a subset of $\mathbb{R}^{d}$. The number $d$ of components of $X$ may be large. The model $m$ may present high complexity as strong non-linearities and high order interaction effects, and it is assumed to be square-integrable, i.e. $m\in L^{2}(\mathcal{X},P_{X})$. Based on the data points $\{(X_{i},Y_{i})\}_{i=1}^{n}$, we estimate a meta-model that approximates the Hoeffding decomposition of $m$. This meta-model belongs to a reproducing kernel Hilbert space (RKHS), which is constructed as a direct sum of the Hilbert spaces leading to an additive decomposition including variables and interactions between them [\citeauthoryearDurrande, Ginsbourger, Roustant, and CarraroDurrande et al.2013]. The estimator of the meta-model is calculated by minimizing an empirical least-squares criterion penalized by the sum of two penalty terms: the Hilbert norm and the empirical norm [\citeauthoryearHuet and TaupinHuet and Taupin2017]. This procedure allows us to select the subsets of variables $X_{1},...,X_{d}$ that contribute to predict $Y$. The estimated meta-model is used to perform sensitivity analysis, and therefore, to determine the influence of each variable and groups of them on the output variable $Y$.

In the classical framework of sensitivity analysis, the function $m$ is calculable in all points of $X$, and one may use the method of [\citeauthoryearSobolSobol1993] to perform sensitivity analysis on $m$. Let us briefly introduce this method: The independency between the components of $X$ leads to write the function $m$ according to its Hoeffding decomposition ([\citeauthoryearSobolSobol1993], [\citeauthoryearVan der VaartVan der Vaart1998]):

[TABLE]

The terms in this decomposition are defined using the conditional expected values:

[TABLE]

These terms are known as the constant term, main effects, interactions of order two and so on. Let $\mathcal{P}$ be the set of all subsets of $\{1,...,d\}$ with dimension $1$ to $d$. For all $v\in\mathcal{P}$ and $X\in\mathcal{X}$, let $X_{v}$ be the vector with components $X_{a}$, $a\in v$. For a set $A$ let $|A|$ be its cardinality, and for all $v\in\mathcal{P}$, let $m_{v}:\mathbb{R}^{|v|}\rightarrow\mathbb{R}$ be the function associated with $X_{v}$ in Equation (2). Then Equation (2) can be expressed as follows:

[TABLE]

This decomposition is unique, all the terms $m_{v}$, $v\in\mathcal{P}$ are centered, and they are orthogonal with respect to $L^{2}(\mathcal{X},P_{X})$. The functions $m$ and $m_{v}$, $v\in\mathcal{P}$ in Equation (3) are square-integrable. As any two terms of decomposition (3) are orthogonal, by squaring (3) and integrating it with respect to the distribution of $X$, a decomposition of the variance of $m(X)$ is obtained as follows:

[TABLE]

The Sobol indices associated with the group of variables $X_{v}$, $v\in\mathcal{P}$ are defined by:

[TABLE]

For each $v$, $S_{v}$ expresses the fraction of variance of $m(X)$ explained by $X_{v}$. For all $v\in\mathcal{P}$, when $|v|=1$, the $S_{v}$s are referred to as the first order indices; when $|v|=2$, i.e. $v=\{a,a^{\prime}\}$ and $a\neq a^{\prime}$, they are referred to as the second order indices or the interaction indices of order two (between $X_{a}$ and $X_{a^{\prime}}$); and the same holds for $|v|>2$.

The total number of the Sobol indices to be calculated is equal to $|\mathcal{P}|=2^{d}-1$, which raises exponentially with the number $d$ of the input variables. When $d$ is large, the evaluation of all the indices can be computationally demanding and even not reachable. In practice, only the indices of order not higher than two are calculated. However, only the first and second order indices may not provide a good information on the model sensitivities. In order to provide better information on the model sensitivities, [\citeauthoryearHomma and SaltelliHomma and Saltelli1996] proposed to calculate the first order and the total indices defined as follows:

Let $\mathcal{P}_{a}\subset\mathcal{P}$ be the set of all the subsets of $\{1,...,d\}$ including $a$, then $S_{T_{a}}=\sum_{v\in\mathcal{P}_{a}}S_{v}$. For all $a\in\{1,...,d\}$, $S_{T_{a}}$ denotes the total effect of $X_{a}$. It expresses the fraction of variance of $m(X)$ explained by $X_{a}$ alone and all the interactions of it with the other variables. The total indices allow us to rank the input variables with respect to the amount of their effect on the output variable. However, they do not provide complete information on the model sensitivities as do all the Sobol indices.

The classical computation of the Sobol indices is based on the Monte Carlo methods (see for example: [\citeauthoryearSobolSobol1993] for the main effect and interaction indices, and [\citeauthoryearSaltelliSaltelli2002] for the main effect and total indices). For models that are expensive to evaluate, the Monte Carlo methods lead to a high computational burden. Moreover, in the case where $d$ is large, $m$ is complex and the calculation of the variances (see Equation (4)) is numerically complicated or not possible (as in the case where the model $m$ is unknown) the methods described above are not applicable. Another approach consists in approximating $m$ by a simplified model, called a meta-model, which is much faster to evaluate and to perform sensitivity analysis on it. Beside the approximations of the Sobol indices of $m$ at a lower computational cost, a meta-model provides a deeper view of the input variables effects on the model output. Among the meta-modelling methods proposed in the literature, the expansion based on polynomial Chaos ([\citeauthoryearWienerWiener1938], [\citeauthoryearSchoutensSchoutens2000]) can be used to approximate the Hoeffding decomposition of $m$ [\citeauthoryearSudretSudret2008]. The principle of the polynomial Chaos is to project $m$ onto a basis of orthonormal polynomials. The polynomial Chaos expansion of $m$ is written as [\citeauthoryearSoize and GhanemSoize and Ghanem2004]:

[TABLE]

where $\{h_{j}\}_{j=0}^{\infty}$ are the coefficients, and $\{\phi_{j}\}_{j=0}^{\infty}$ are the multivariate orthonormal polynomials associated with $X$ which are determined according to the distribution of the components of $X$. In practice, expansion (6) shall be truncated for computational purposes, and the model $m$ may be approximated by $\sum_{j=0}^{v_{max}}h_{j}\phi_{j}(X)$, where $v_{max}$ is determined using a truncation scheme. The Sobol indices are obtained then by summing up the squares of the suitable coefficients. [\citeauthoryearBlatman and SudretBlatman and Sudret2011] proposed a method for truncating the polynomial Chaos expansion and an algorithm based on the least angle regression for selecting the terms in the expansion. In this approach, according to the distribution of the components of $X$, a unique family of orthonormal polynomials $\{\phi_{j}\}_{j=0}^{\infty}$ is determined. However, this family may not be necessarily the best functional basis to approximate $m$ well.

Gaussian Process (GP) can also be used to construct meta-models as highlighted in [\citeauthoryearWelch, Buck, Sacks, Wynn, Mitchell, and MorrisWelch et al.1992], [\citeauthoryearOakley and O’HaganOakley and O’Hagan2004], [\citeauthoryearKleijnenKleijnen2007, \citeauthoryearKleijnenKleijnen2009], [\citeauthoryearMarrel, Iooss, Laurent, and RoustantMarrel et al.2009], [\citeauthoryearDurrande, Ginsbourger, and RoustantDurrande et al.2012], [\citeauthoryearLe Gratiet, Cannamela, and IoossLe Gratiet et al.2014]. The principle is to consider that the prior knowledge about the function $m(X)$, can be modelled by a GP $\mathcal{Z}(X)$ with a mean $m_{\mathcal{Z}}(X)$ and a covariance kernel $k_{\mathcal{Z}}(X,X^{\prime})$. To perform sensitivity analysis from a GP model, one may replace the true model $m(X)$ with the mean of the conditional GP and deduce the Sobol indices from it. A review on the meta-modelling based on polynomial Chaos and GP is presented in [\citeauthoryearLe Gratiet, Marelli, and SudretLe Gratiet et al.2017].

[\citeauthoryearDurrande, Ginsbourger, Roustant, and CarraroDurrande et al.2013] considered a class of the functional approximation methods similar to the GP and obtained a meta-model that satisfies the properties of the Hoeffding decomposition. They proposed to approximate $m$ by functions belonging to a RKHS $\mathcal{H}$ which is a direct sum of the Hilbert spaces. Their RKHS $\mathcal{H}$ is constructed in a way that the projection of $m$ onto $\mathcal{H}$, denoted $f^{*}$, is an approximation of the Hoeffding decomposition of $m$. The function $f^{*}$ is defined as the minimizer over the functions $f\in\mathcal{H}$ of the criterion $E_{X}(m(X)-f(X))^{2}.$

Let $<.,.>_{\mathcal{H}}$ be the scalar product in $\mathcal{H}$, let also $k$ and $k_{v}$ be the reproducing kernels associated with the RKHS $\mathcal{H}$ and the RKHS $\mathcal{H}_{v}$, respectively. The properties of the RKHS $\mathcal{H}$ insures that any function $f\in\mathcal{H}$, $f:\mathcal{X}\subset\mathbb{R}^{d}\rightarrow\mathbb{R}$ is written as the following decomposition:

[TABLE]

where $f_{0}$ is constant, and $f_{v}:\mathbb{R}^{|v|}\rightarrow\mathbb{R}$ is defined by $f_{v}(X)=<f,k_{v}(X,.)>_{\mathcal{H}}$. For more details on the RKHS construction and the definition of the Hilbert norm see Section 6.1.1 in 6.1.

For all $v\in\mathcal{P}$, the functions $f_{v}(X_{v})$ are centered and for all $v\neq v^{\prime}$, the functions $f_{v}(X_{v})$ and $f_{v^{\prime}}(X_{v^{\prime}})$ are orthogonal with respect to $L^{2}(\mathcal{X},P_{X})$. Therefore, the decomposition of the function $f$ presented in Equation (7) is its Hoeffding decomposition. As the function $f^{*}$ belongs to the RKHS $\mathcal{H}$, it is decomposed as its Hoeffding decomposition, $f^{*}=f^{*}_{0}+\sum_{v\in\mathcal{P}}f^{*}_{v}$, and each function $f^{*}_{v}$ approximates the function $m_{v}$ in Equation (3). The number of terms $f^{*}_{v}$ that should be estimated in the Hoeffding decomposition of $f^{*}$ is equal to $|\mathcal{P}|=2^{d}-1$, which may be huge since it rises very quickly by increasing $d$. In order to deal with this problem, in the regression framework, one may estimate $f^{*}$ by a sparse meta-model $\widehat{f}\in\mathcal{H}$. To this end, the estimation of $f^{*}$ is done on the basis of $n$ observations by minimizing a least-squares criterion suitably penalized in order to deal with both the non-parametric nature of the problem and the possibly large number of functions that have to be estimated. In the classical framework of the sensitivity analysis one may calculate a sparse approximation of $f^{*}$ using least-squares penalized criterion as it is done in the non-parametric regression framework. In order to obtain a sparse solution of a minimization problem, the penalty function should enforce the sparsity. There exists various ways of enforcing sparsity for a minimization (maximization) problem, see for example [\citeauthoryearHastie, Tibshirani, and WainwrightHastie et al.2015] for a review. Some methods, such as the Sparse Additive Models (SpAM) procedure ([\citeauthoryearRavikumar, Lafferty, Liu, and WassermanRavikumar et al.2009], [\citeauthoryearLiu, Wasserman, and LaffertyLiu et al.2009]) are based on a combination of the $l_{1}$-norm with the empirical $L^{2}$-norm: $\|f\|_{n,1}=\sum_{a=1}^{d}\|f_{a}\|_{n},$ where $\|f_{a}\|_{n}^{2}=\sum_{i=1}^{n}f_{a}^{2}(X_{ai})/n,$ is the squared empirical $L^{2}$-norm of the univariate function $f_{a}$. The Component Selection and Smoothing Operator (COSSO) method developed by [\citeauthoryearLin and ZhangLin and Zhang2006] enforces sparsity using a combination of the $l_{1}$-norm with the Hilbert norm: $\|f\|_{\mathcal{H},1}=\sum_{a=1}^{d}\|f_{a}\|_{\mathcal{H}_{a}}$. Instead of focusing on only one penalty term, one may consider a more general family of estimators, called doubly penalized estimator, which is obtained by minimizing a criterion penalized by the sum of two penalty terms. [\citeauthoryearRaskutti, Wainwright, and YuRaskutti et al.2009, \citeauthoryearRaskutti, Wainwright, and YuRaskutti et al.2012] proposed a doubly penalized estimator, which is the solution of the minimization of a least-squares criterion penalized by the sum of a sparsity penalty term and a combination of the $l_{1}$-norm with the Hilbert norm:

[TABLE]

where $\gamma,\mu\in\mathbb{R}$ are the tuning parameters that should be suitably chosen.

[\citeauthoryearMeier, van de Geer, and BühlmannMeier et al.2009] proposed a related family of estimators, based on the penalization with the empirical $L^{2}$-norm. Their penalty function is the sum of the sparsity penalty term, $\|f\|_{n,1}$, and a smoothness penalty term. [\citeauthoryearHuet and TaupinHuet and Taupin2017] considered the same approximation functional spaces as [\citeauthoryearDurrande, Ginsbourger, Roustant, and CarraroDurrande et al.2013], and obtained a doubly penalized estimator of a meta-model which approximates the Hoeffding decomposition of $m$. Their estimator is the solution of the least-squares minimization penalized by the penalty function defined in Equation (8) adapted to the multivariate setting,

[TABLE]

This procedure, called RKHS ridge group sparse, estimates the groups $v$ that are suitable for predicting $f^{*}$, and the relationship between $f^{*}_{v}$ and $X_{v}$ for each group. The obtained estimator, called RKHS meta-model, is used then to estimate the Sobol indices of $m$. This approach renders it possible to estimate the Sobol indices for all groups in the support of the RKHS meta-model, including the interactions of possibly high order, a point known to be difficult in practice.

In this paper, we introduce an R package, called RKHSMetaMod, that implements the RKHS ridge group sparse procedure. The functions of this package allows us to:

• (1)

calculate the reproducing kernels and their associated Gram matrices (see Section 3.1),

• (2)

implement the RKHS ridge group sparse procedure and a special case of it, called the RKHS group lasso procedure (when $\gamma=0$ in the penalty function (9)), in order to estimate the terms $f^{*}_{v}$ in the Hoeffding decomposition of the meta-model $f^{*}$ leading to an estimation of the function $m$ (see Section 3.2),

• (3)

choose the tuning parameters $\mu$ and $\gamma$ (see Equation (9)), using a procedure that leads to obtain the best RKHS meta-model in terms of the prediction quality,

• (4)

estimate the Sobol indices of the function $m$ (see Section 2.3).

The current version of the package supports uniformly distributed input variables on $\mathcal{X}=[0,1]^{d}$. However, it could easily be adapted to a different distribution of input variables by making a small modification to one of its functions (see Remark 3 of Section 3.1).

Let us give a brief overview of the related existing statistical packages to RKHSMetaMod. The R package sensitivity is designed to implement sensitivity analysis methods and provides the approaches for numerical calculation of the Sobol indices. In particular, Kriging method can be used to reduce the number of observations in global sensitivity analysis. The function sobolGP of the package sensitivity builds a Kriging based meta-model using the function km of the package DiceKriging [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012], and estimates its Sobol indices. This procedure can also be done using the function km and the function fast99 of the package sensitivity (see Section 4.5. of [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012]). In this case, the idea is once again to build a Kriging based meta-model using the function km and then estimate its Sobol indices using the function fast99. In both cases the true function is substituted by a Kriging based meta-model and then its Sobol indices are estimated. In the sobolGP function, the Sobol indices are estimated by Monte Carlo integration, while the fast99 function estimates them using the extended-FAST method [\citeauthoryearSaltelli, Tarantola, and ChanSaltelli et al.1999]. To reduce the computational burden when dealing with large datasets and complex models, in RKHSMetaMod package, we propose to use the empirical variances to estimate the Sobol indices (see Section 2.3). Besides, the estimation of the Sobol indices in the RKHSMetaMod package is done based on the RKHS meta-model which is a sparse estimator. It is beneficial since instead of calculating the Sobol indices of all groups $v\in\mathcal{P}$, only the Sobol indices associated with the groups in the support of the RKHS meta-model are computed (see Section 2.3). Moreover, the functions sobolGP and fast99 provide the estimation of the first order and the total Sobol indices only, while the procedure in the RKHSMetaMod package makes it possible to estimate the high order Sobol indices. The R packages DiceKriging and DiceOptim (Deep Inside Computer Experiments Kriging/Optim) [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012] implement the Kriging based meta-models to estimate complex models in the high dimensional context. They provide different GP (Kriging) models corresponding to the Gaussian, Matérn, Exponential and Power-Exponential correlation functions. The estimation of the parameters of the correlation functions in these packages relies on the global optimizer with gradient genoud algorithm of the package rgenoud [\citeauthoryearMebane and SekhonMebane and Sekhon2011]. These packages do not implement any method of the sensitivity analysis themselves. However, some authors (see Section 4.5. of [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012] for example) perform sensitivity analysis on their estimated meta-models by employing the functions of the package sensitivity. The R package RobustGaSP (Robust Gaussian Stochastic Process) [\citeauthoryearGu, Palomo, and BergerGu et al.2019] approximates a complex model by a GP meta-model. This package implements marginal posterior mode estimation of the GP model parameters. The estimation method in this package insures the robustness of the parameter estimation in the GP model, and allows one also to identify input variables that have no effect on the variability of the function under study. The R package mlegp (maximum likelihood estimates of Gaussian processes) [\citeauthoryearDancik and DormanDancik and Dorman2008] provides functions to implement both meta-modelling approaches and sensitivity analysis methods. It obtains maximum likelihood estimates of GP model for the output of costly computer experiments. The GP models are built either on the basis of the Gaussian correlation function or on the basis of the first degree polynomial trend. The sensitivity analysis methods implemented in this package include Functional Analysis of Variance (FANOVA) decomposition, plot functions to obtain diagnostic plots, main effects, and second order interactions. The prediction quality of the meta-model depends on the quality of the estimation of its parameters and more precisely the estimation of parameters in the correlation functions [\citeauthoryearKennedy and O’HaganKennedy and O’Hagan2000]. The maximum likelihood estimation of these parameters often produce unstable results, and as a consequence, the obtained meta-model may have an inferior prediction quality ([\citeauthoryearGu, Wang, and BergerGu et al.2018], [\citeauthoryearGuGu2019]). The RKHSMetaMod package is devoted to the meta-model estimation on the RKHS $\mathcal{H}$. It implements the convex optimization algorithms to calculate meta-models; provides the functions to compute the prediction error of the obtained meta-models; performs the sensitivity analysis on the obtained meta-models and more precisely calculate their Sobol indices. The convex optimization algorithms used in this package are all written using C++ libraries, and are adapted to take into account the problem of high dimensionality in this context. This package is available from the Comprehensive R Archive Network (CRAN) [\citeauthoryearKamariKamari2019].

The organization of the paper is as follows: In the next Section, we describe the estimation method. In Section 3, we present in details the algorithms used in the RKHSMetaMod package to obtain the RKHS meta-model. In Section 6.2, an overview of the package functions as well as a brief documentation of them are given. Section 4 includes two parts: In the first part, Section 4.1, the performance of the RKHSMetaMod package functions is validated through a simulation study. In the second part, Section 4.2, the comparison in terms of the predictive accuracy between the RKHS meta-model and the Kriging based meta-models from RobustGaSP [\citeauthoryearGu, Palomo, and BergerGu et al.2019] and DiceKriging [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012] packages is given through two examples.

2 Estimation method

In this Section, we present: the RKHS ridge group sparse and the RKHS group lasso procedures (see 2.1), the strategy of choosing the tuning parameters in the RKHS ridge group sparse algorithm (see 2.2), and the calculation of the empirical Sobol indices of the RKHS meta-model (see 2.3).

2.1 RKHS ridge group sparse and RKHS group lasso procedures

Let us denote by $n$ the number of observations. The dataset consists of a vector of $n$ observations $Y=(Y_{1},...,Y_{n})$, and a $n\times d$ matrix of features $X$ with components $(X_{ai},i=1,...,n,a=1,...,d)\in\mathbb{R}^{n\times d}.$ For some tuning parameters $\gamma_{v}$, $\mu_{v}$, $v\in\mathcal{P}$, the RKHS ridge group sparse criterion is defined by,

[TABLE]

where $X_{v}$ represents the matrix of variables corresponding to the $v$-th group, i.e. $X_{v}=(X_{vi},i=1,...,n,v\in\mathcal{P})\in\mathbb{R}^{n\times|\mathcal{P}|},$ and where $\|f_{v}\|_{n}$ is the empirical $L^{2}$-norm of $f_{v}$ defined by the sample $\{X_{vi}\}_{i=1}^{n}$ as $\|f_{v}\|_{n}=\sqrt{\sum_{i=1}^{n}f_{v}^{2}(X_{vi})/n}.$

The penalty function in the criterion (10) is the sum of the Hilbert norm and the empirical norm, which allows us to select few terms in the additive decomposition of $f$ over sets $v\in\mathcal{P}$. Moreover, the Hilbert norm favours the smoothness of the estimated $f_{v}$, $v\in\mathcal{P}$.

Let $\mathcal{F}=\{f:\>f=f_{0}+\sum_{v\in\mathcal{P}}f_{v},\>f_{v}\in\mathcal{H}_{v},\>\|f_{v}\|_{\mathcal{H}_{v}}\leq r_{v},\>r_{v}\in\mathbb{R}^{+}\}$ be the set of functions. Then the RKHS meta-model is defined by,

[TABLE]

According to the Representer Theorem [\citeauthoryearKimeldorf and WahbaKimeldorf and Wahba1970], the non-parametric functional minimization problem described above is equivalent to a parametric minimization problem. Indeed, the solution of the minimization problem (11) belonging to the RKHS $\mathcal{H}$ is written as $f=f_{0}+\sum_{v\in\mathcal{P}}f_{v}$, where for some matrix $\theta=(\theta_{vi},i=1,...,n,v\in\mathcal{P})\in\mathbb{R}^{n\times|\mathcal{P}|}$ we have for all $v\in\mathcal{P}$,

[TABLE]

Let $\|.\|$ be the Euclidean norm (called also $L^{2}$-norm) in $\mathbb{R}^{n}$, and for each $v\in\mathcal{P}$, let $K_{v}$ be the $n\times n$ Gram matrix associated with the kernel $k_{v}(.,.)$, i.e. $(K_{v})_{i,i^{\prime}}=k_{v}(X_{vi},X_{vi^{\prime}})$. Let also $K^{1/2}$ be the matrix that satisfies $t(K^{1/2})K^{1/2}=K$, and let $\widehat{f}_{0}$ and $\widehat{\theta}$ be the minimizers of the following penalized least-squares criterion:

[TABLE]

Then the estimator $\widehat{f}$ defined in Equation (11) satisfies,

[TABLE]

Remark 1**.**

The constraint $\|f_{v}\|_{\mathcal{H}_{v}}\leq r_{v}$ is crucial for theoretical properties, but the value of $r_{v}$ is generally unknown and has no practical usefulness. In this package, it is not taken into account in the parametric minimization problem.

For each $v\in\mathcal{P}$, let $\gamma^{\prime}_{v}$ and $\mu^{\prime}_{v}$ be the weights that are chosen suitably. We define $\gamma_{v}=\gamma\times\gamma^{\prime}_{v}$ and $\mu_{v}=\mu\times\mu^{\prime}_{v}$ with $\gamma,\mu\in\mathbb{R}^{+}$.

Remark 2**.**

This formulation simplifies the choice of the tuning parameters since instead of tuning $2\times|\mathcal{P}|$ parameters $\gamma_{v}$ and $\mu_{v}$, $v\in\mathcal{P}$, only two parameters $\gamma$ and $\mu$ are tuned. Moreover, the weights $\gamma^{\prime}_{v}$ and $\mu^{\prime}_{v}$, $v\in\mathcal{P}$ may be of interest in practice. For example, one can take weights that increase with the cardinal of $v$ in order to favour the effects with small interaction order between variables.

For the sake of simplicity, in the rest of this paper for all $v\in\mathcal{P}$ the weights $\gamma^{\prime}_{v}$ and $\mu^{\prime}_{v}$ are assumed to be set as one, and the RKHS ridge group sparse criterion is then expressed as follows:

[TABLE]

If we consider only the second part of the penalty function in the criterion (13) ( i.e. set $\gamma=0$), we obtain the RKHS group lasso criterion as follows:

[TABLE]

which is a group lasso criterion [\citeauthoryearYuan and LinYuan and Lin2006] up to a scale transformation.

In the RKHSMetaMod package, the RKHS ridge group sparse algorithm is initialized using the solutions obtained by solving the RKHS group lasso algorithm. Indeed, the penalty function in the RKHS group lasso criterion (14) insures the sparsity in the solution. Therefore, for a given value of $\mu$, by implementing the RKHS group lasso algorithm (see Section 3.2.1), a RKHS meta-model with few terms in its additive decomposition is obtained. The support and the coefficients of a RKHS meta-model which is obtained by implementing RKHS group lasso algorithm will be denoted by $\widehat{S}_{\widehat{f}_{\text{Group Lasso}}}$ and $\widehat{\theta}_{\text{Group Lasso}}$, respectively. From now on the tuning parameter in the RKHS group lasso criterion will be denoted by:

[TABLE]

2.2 Choice of the tuning parameters

While dealing with an optimization problem of a criterion of the form (13), one of the essential steps is to choose appropriately the tuning parameters. Cross-validation is generally used for that purpose. Nevertheless in the context of high-dimensional complex models, the computational time for a cross-validation procedure may be prohibitively high. Therefore we propose a procedure based on a single testing data set:

• •

we first choose, a grid of values of the tuning parameters $\mu$ and $\gamma$;

Let $\mu_{\text{max}}$ be the smallest value of $\mu_{g}$ (see Equation (15)), such that the solution to the minimization of the RKHS group lasso problem for all $v\in\mathcal{P}$ is $\theta_{v}=0$. We have,

[TABLE]

In order to set up the grid of values of $\mu$, one may find $\mu_{\text{max}}$ and then a grid of values of $\mu$ could be defined by $\mu_{l}=\mu_{\text{max}}/(\sqrt{n}\times 2^{l})$ for $l\in\{1,...,l_{\text{max}}\}.$ The grid of values of $\gamma$ is chosen by the user.

• •

next, for the grid of values of $\mu$ and $\gamma$, we calculate a sequence of estimators. Each estimator associated with the pair $(\mu,\gamma)$ in the grid of values of $\mu$ and $\gamma$, denoted by $\widehat{f}_{(\mu,\gamma)}$, is the solution of the RKHS ridge group sparse optimization problem or the RKHS group lasso optimization problem if $\gamma=0$.

• •

finally, the obtained estimators $\widehat{f}_{(\mu,\gamma)}$ are evaluated using a testing dataset, $\{(Y^{\text{test}}_{i},X^{\text{test}}_{i})\}_{i=1}^{n^{\text{test}}}.$ The prediction error associated with the estimator $\widehat{f}_{(\mu,\gamma)}$ is calculated by,

[TABLE]

where for $S_{\widehat{f}}$ being the support of the estimator $\widehat{f}_{(\mu,\gamma)}$ we have,

[TABLE]

The pair $(\widehat{\mu},\widehat{\gamma})$ with the smallest value of the prediction error is chosen, and the estimator $\widehat{f}_{(\widehat{\mu},\widehat{\gamma})}$ is considered as the best estimator of the function $m$, in terms of the prediction error.

In the RKHSMetaMod package, the algorithm to calculate a sequence of the RKHS meta-models, the value of $\mu_{\text{max}}$, and the prediction error are implemented as RKHSMetMod, mu$\_$max, and PredErr functions, respectively. These functions are described in Section 6.2, and illustrated in Example 1, Example 2, and Examples 1, 2, 3, respectively.

2.3 Estimation of the Sobol indices

The variance of the function $m$ is estimated by the variance of the estimator $\widehat{f}$. As the estimator $\widehat{f}$ belongs to the RKHS $\mathcal{H}$, it admits the Hoeffding decomposition and,

[TABLE]

In order to reduce the computational cost in practice, one may estimate the variances of $\widehat{f}_{v}(X_{v})$, $v\in\mathcal{P}$ by their empirical variances. Let $\widehat{f}_{v.}$ be the empirical mean of $\widehat{f}_{v}(X_{vi})$, $i=1,...,n$, then:

[TABLE]

For the groups $v$ that do not belong to the support of $\widehat{f}$, we have $\widehat{S}_{v}=0$ and for the groups $v$ that belong to the support of $\widehat{f}$, the estimators of the Sobol indices of $m$ are defined by,

[TABLE]

In the RKHSMetaMod package, the algorithm to calculate the empirical Sobol indices $\widehat{S}_{v}$, $v\in\mathcal{P}$ is implemented as SI$\_$emp function. This function is described in Section 6.2.2 and illustrated in Examples 1, 2, 3.

3 Algorithms

The RKHSMetaMod package implements two optimization algorithms: the RKHS ridge group sparse (see Algorithm 2) and the RKHS group lasso (see Algorithm 1). These algorithms rely on the Gram matrices $K_{v}$, $v\in\mathcal{P}$ that have to be positive definite. Therefore, the first and essential step in this package is to calculate these matrices and insure their positive definiteness. This step is detailed in an algorithm that is described in the next Section. The second step is to estimate the RKHS meta-model. In the RKHSMetaMod package, two different objectives based on different procedures are considered to calculate this estimator:

The RKHS meta-model with the best prediction quality.

The procedure to calculate the RKHS meta-model with the best prediction quality has been described in Section 2.2: a sequence of values of the tuning parameters $(\mu,\gamma)$ is considered, and the RKHS meta-models associated with each pair of the values of $(\mu,\gamma)$ are calculated. For $\gamma=0$, the RKHS meta-model is obtained by solving the RKHS group lasso optimization problem, while for $\gamma\neq 0$ the RKHS ridge group sparse optimization problem is solved to calculate the RKHS meta-model. The obtained RKHS meta-models are evaluated by considering a new dataset and the RKHS meta-model with the minimum value of the prediction error is chosen as the best estimator.

• 2.

The RKHS meta-model with at most $qmax$ groups in its support, i.e. $|S_{\widehat{f}}|\leq qmax$.

First, the tuning parameter $\gamma$ is set as zero. Then, a value of $\mu$ for which the number of groups $v\in\mathcal{P}$ in the solution of the RKHS group lasso optimization problem is equal to $qmax$, is computed. This value of $\mu$ will be denoted by $\mu_{qmax}$. Finally, the RKHS meta-models containing at most $qmax$ groups in their support are obtained by implementing the RKHS ridge group sparse algorithm for a grid of values of $\gamma\neq 0$ and $\mu_{qmax}$. This procedure is described in more details in Section 3.2.3.

3.1 Calculation of the Gram matrices

The available kernels in the RKHSMetaMod package are: Gaussian kernel, Matérn $3/2$ kernel, Brownian kernel, quadratic kernel and linear kernel. The usual presentation of these kernels is given in Table 1.

The choice of the kernel that is done by the user determines the functional approximation space. For a chosen kernel, the algorithm to calculate the Gram matrices $K_{v}$, $v\in\mathcal{P}$ in the RKHSMetaMod package, is implemented as calc$\_$Kv function. This algorithm is based on three essential points:

• (1)

Modify the chosen kernel:

In order to satisfy the conditions of constructing the RKHS $\mathcal{H}$ described in Section 6.1.1, these kernels are modified according to Equation (19). Let us take the example of the Brownian kernel:

The RKHS associated with the Brownian kernel $k_{a}(X_{a},X_{a}^{\prime})=\min(X_{a},X_{a}^{\prime})+1$ is well known to be $\mathcal{H}_{a}=\{f:[0,1]\rightarrow\mathbb{R}\mbox{ is absolutely continuous, and }f(0)=0,\>\int_{0}^{1}{f^{\prime}(X_{a})}^{2}dX_{a}<\infty\},$ with the inner product $\langle f,h\rangle_{\mathcal{H}_{a}}=\int_{0}^{1}f^{\prime}(X_{a})h^{\prime}(X_{a})dX_{a}.$ Easy calculations lead to obtain the Brownian kernel as follows,

[TABLE]

The RKHS associated with the kernel $k_{0a}$ is the set $\mathcal{H}_{0a}=\{f\in\mathcal{H}_{a}:\>\int_{0}^{1}f(X_{a})dX_{a}=0\}$, and we have $\mathcal{H}=\mathbbm{1}+\sum_{v\in\mathcal{P}}\mathcal{H}_{v}=\{f:[0,1]^{d}\rightarrow\mathbb{R}:\>f=f_{0}+\sum_{v\in\mathcal{P}}f_{v}(X_{v}),\mbox{ with }f_{v}\in\mathcal{H}_{v}\}.$

Remark 3**.**

In the current version of the package, we consider the input variables $X=(X_{1},...,X_{d})$ that are uniformly distributed on $[0,1]^{d}$. In order to consider the input variables that are not distributed uniformly, it suffices to modify a part of the function calc$\_$Kv related to the calculation of the kernels $k_{0a}$, $a=1,...,d$. For example, for $X=(X_{1},...,X_{d})$ being distributed with law $P_{X}=\prod_{a=1}^{d}P_{a}$ on $\mathcal{X}=\bigotimes_{a=1}^{d}\mathcal{X}_{a}\subset\mathbb{R}^{d}$, the kernel $k_{0a}$ associated with the Brownian kernel is calculated as follows,

[TABLE]

The other parts of the function calc$\_$Kv remain unchanged.

• (2)

Calculate the Gram matrices $K_{v}$ for all $v$:

First, for all $a=1,...,d$, the Gram matrices $K_{a}$ associated with kernels $k_{0a}$ are calculated using Equation (19), $(K_{a})_{i,i^{\prime}}=k_{0a}(X_{ai},X_{ai^{\prime}}).$ Then, for all $v\in\mathcal{P}$, the Gram matrices $K_{v}$ associated with kernel $k_{v}=\prod_{a\in v}k_{0a}$ are computed by $K_{v}=\bigodot_{a\in v}K_{a},$

• (3)

Insure the positive definiteness of the matrices $K_{v}$:

The output of the function calc$\_$Kv is one of the input arguments of the functions associated with the RKHS group lasso and the RKHS ridge group sparse algorithms. Throughout these algorithms we need to calculate the inverse and the square root of the matrices $K_{v}$. In order to avoid the numerical problems and insure the invertibility of the matrices $K_{v}$, it is mandatory to have these matrices positive definite. One way to render the matrices $K_{v}$ positive definite is to add a nugget effect to them. That is, to modify matrices $K_{v}$ by adding a diagonal with a constant term, i.e. $K_{v}+\mbox{epsilon}\times I_{n}$. The value of epsilon is computed based on the data and through a part of the algorithm of the function calc$\_$kv. Let us briefly explain this part of the algorithm:

For each group $v\in\mathcal{P}$, let $\lambda_{v,i},i=1,...,n$ be the eigenvalues associated with the matrix $K_{v}$. Set $\lambda_{v,\text{max}}={\text{max}}_{i}\lambda_{v,i}$ and $\lambda_{v,\text{min}}={\text{min}}_{i}\lambda_{v,i}$. For some fixed value of tolerance "tol", and for each matrix $K_{v}$, if $\lambda_{v,\text{min}}<\lambda_{v,\text{max}}\times\text{tol}$: then, the eigenvalues of $K_{v}$ are replaced by $\lambda_{v,i}+\text{epsilon}$, with "epsilon" being equal to $\lambda_{v,\text{max}}\times$"tol". The value of "tol" is set as $1e^{-8}$ by default, but one may consider a smaller or a greater value for it depending on the kernel chosen and the value of $n$.

The function calc$\_$Kv is described in Section 6.2.2 and illustrated in Example 2.

3.2 Optimization algorithms

The RKHS meta-model is the solution of one of the optimization problems: the minimization of the RKHS group lasso criterion presented in Equation (14) (if $\gamma=0$), or the minimization of the RKHS ridge group sparse criterion presented in Equation (13) (if $\gamma\neq 0$). In the following, the algorithms to solve these optimization problems are presented.

3.2.1 RKHS group lasso

A popular technique for doing group wise variable selection is group lasso. With this procedure, depending on the value of the tuning parameter $\mu$, an entire group of predictors may drop out of the model. An efficient algorithm for solving group lasso problem is the classical block coordinate descent algorithm ([\citeauthoryearBoyd, Parikh, Chu, Peleato, and EcksteinBoyd et al.2011], [\citeauthoryearBubeckBubeck2015]). Following the idea of [\citeauthoryearFuFu1998], [\citeauthoryearYuan and LinYuan and Lin2006] implemented a block wise descent algorithm for the group lasso penalized least-squares under the condition that the model matrices in each group are orthonormal. A block coordinate (gradient) descent algorithm for solving the group lasso penalized logistic regression is then developed by [\citeauthoryearMeier, van de Geer, and BühlmannMeier et al.2008]. This algorithm is implemented in the R package grplasso available from CRAN [\citeauthoryearMeierMeier2020]. [\citeauthoryearYang and ZouYang and Zou2015] proposed a unified algorithm named group wise majorization descent for solving the general group lasso learning problems by assuming that the loss function satisfies a quadratic majorization condition. The implementation of their work is done in the gglasso R package available from CRAN [\citeauthoryearYang, Zou, and BhatnagarYang et al.2020].

In order to solve the RKHS group lasso optimization problem, we use the classical block coordinate descent algorithm. The minimization of criterion $C_{g}(f_{0},\theta)$ (see Equation (14)) is done along each group $v$ at a time. At each step of the algorithm, the criterion $C_{g}(f_{0},\theta)$ 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. This procedure leads to Algorithm 1 (see 6.1 for more details on this procedure).

In the RKHSMetaMod package, the Algorithm 1 is implemented as RKHSgrplasso function. This function is described in Section 6.2.2.

3.2.2 RKHS ridge group sparse

In order to solve the RKHS ridge group sparse optimization problem, we propose an adapted block coordinate descent algorithm. This algorithm is provided in two steps:

• Step $1$

Initialize the input parameters by the solutions of the RKHS group lasso algorithm for each value of the tuning parameter $\mu$, and implement the RKHS ridge group sparse algorithm through the active support of the RKHS group lasso solutions until it achieves convergence. This step is provided in order to decrease the execution time. In fact, instead of implementing the RKHS ridge group sparse algorithm over the set of all groups $\mathcal{P}$, it is implemented only over the groups in the support of the solution of the RKHS group lasso algorithm, $\widehat{S}_{\widehat{f}_{\text{Group Lasso}}}$.

• Step $2$

Re-initialize the input parameters with the obtained solutions of Step $1$ and implement the RKHS ridge group sparse algorithm through all groups in $\mathcal{P}$ until it achieves convergence. This second step makes it possible to verify that no group is missing in the output of Step $1$.

This procedure leads to Algorithm 2 (see 6.1 for more details on this procedure).

In the RKHSMetaMod package the Algorithm 2 is implemented as pen$\_$MetMod function. This function is described in Section 6.2.2 and illustrated in Example 2.

3.2.3 RKHS meta-model with at most $qmax$ groups in its support

By considering some prior information about the data, one may be interested in a RKHS meta-model $\widehat{f}$ with the number of groups in its support not greater than some "$qmax$". In order to obtain such an estimator, we provide the following procedure in the RKHSMetaMod package:

• •

First, the tuning parameter $\gamma$ is set as zero and a value of $\mu$ for which the solution of the RKHS group lasso algorithm, Algorithm 1, contains exactly $qmax$ groups in its support is computed. This value is denoted by $\mu_{qmax}$.

• •

Then, the RKHS ridge group sparse algorithm, Algorithm 2, is implemented by setting the tuning parameter $\mu$ equal to $\mu_{qmax}$ and a grid of values of the tuning parameter $\gamma>0$.

This procedure leads to Algorithm 3.

This algorithm is implemented in the RKHSMetaMod package, as function RKHSMetMod$\_$qmax (see Section "Main RKHSMetaMod functions" (supplementary materials) for more details on this function).

Remark 4**.**

As both terms in the penalty function of criterion (13) enforce sparsity to the solution, the estimator obtained by solving the RKHS ridge group sparse associated with the pair of the tuning parameters $(\mu_{qmax},\gamma>0)$ may contain a smaller number of groups than the solution of the RKHS group lasso optimization problem (i.e. the RKHS ridge group sparse with $(\mu_{qmax},\gamma=0)$). And therefore, the estimated RKHS meta-model contains at most "$qmax$" groups in its support.

4 RKHSMetaMod through examples

4.1 Simulation study

Let us consider the g-function of Sobol [\citeauthoryearSaltelli, Chan, and ScottSaltelli et al.2009] in the Gaussian regression framework, i.e. $Y=m(X)+\varepsilon$. The error term $\varepsilon$ is a centered Gaussian random variable with variance $\sigma^{2}$, and the function $m$ is the g-function of Sobol defined over $[0,1]^{d}$ by,

[TABLE]

The Sobol indices of the g-function can be expressed analytically:

[TABLE]

Set $c_{1}=0.2$, $c_{2}=0.6$, $c_{3}=0.8$ and $(c_{a})_{a>3}=100$. With these values of coefficients $c_{a}$, the variables $X_{1},X_{2}$ and $X_{3}$ explain $99.98\%$ of the variance of the function $m(X)$ (see Table 4).

In this Section, three examples are presented. In all examples, the value of Dmax is set as three. Example 1 illustrates the use of the RKHSMetMod function by considering three different kernels, "matern", "brownian", and "gaussian" (see Table 1), and three datasets of $n\in\{50,100,200\}$ observations and $d=5$ input variables. The larger datasets with $n\in\{1000,2000,5000\}$ observations and $d=10$ input variables are studied in Examples 2 and 3.

In each example, two independent datasets are generated: $(X,Y)$ to estimate the meta-models, and $(XT,YT)$ to estimate the prediction errors. The design matrices $X$ and $XT$ are the Latin Hypercube Samples of the input variables that are generated using maximinLHS function of the package lhs available at CRAN [\citeauthoryearCarnellCarnell2021]:

R> library(lhs); X <- maximinLHS(n, d); XT <- maximinLHS(n, d)

The response variables $Y$ and $YT$ are calculated as $Y=m(X)+\varepsilon$ and $YT=m(XT)+\varepsilon_{T}$, where $\varepsilon$, and $\varepsilon_{T}$ are centered Gaussian random variables with $\sigma^{2}=(0.2)^{2}$:

R> a <- c(0.2, 0.6, 0.8, 100, 100, 100, 100, 100, 100, 100)[1:d]

R> g=1; for (i in 1:d) g = g*(abs(4*X[,i]-2)+a[i])/(1+a[i])

R> sigma <- 0.2

R> epsilon <- rnorm(n, 0, sigma*sigma); Y <- g + epsilon

R> gT=1; for (i in 1:d) gT = gT*(abs(4*XT[,i]-2)+a[i])/(1+a[i])

R> epsilonT <- rnorm(n, 0, sigma*sigma); YT <- gT + epsilonT

Example 1**.**

RKHS meta-model estimation using RKHSMetMod function:

In this example, three datasets of $n$ points maximinLHS over $[0,1]^{d}$ are generated with $n\in\{50,100,200\}$ and $d=5$, and a grid of five values for each of the tuning parameters $\mu$ and $\gamma$ is considered as follows:

[TABLE]

For each dataset, the experiment is repeated $N_{r}=50$ times. At each repetition, the RKHS meta-models associated with the pair of the tuning parameters $(\mu,\gamma)$ are estimated using the RKHSMetMod function:

R> Dmax <- 3; kernel <- "matern" # kernel <- "brownian" # kernel <- "gaussian"

R> gamma <- c(0.2, 0.1, 0.01, 0.005, 0); frc <- c(4,8,16,32,64)

R> res <- RKHSMetMod(Y, X, kernel, Dmax, gamma, frc, FALSE)

These meta-models are evaluated using a testing dataset. The prediction errors are computed for them using the PredErr function. The RKHS meta-model with minimum prediction error is chosen to be the best estimator for the model. Finally, the Sobol indices are computed for the best RKHS meta-model using the function SI$\_$emp:

R> Err <- PredErr(X, XT, YT, mu, gamma, res, kernel, Dmax)

R> SI <- SI_emp(res, Err)

The vector "mu" is the values of the tuning parameter $\mu$ that are calculated throughout the function RKHSMetMod (see argument "frc" in Table 12). It could be recovered from the output of the RKHSMetMod function as follows:

R> mu <- vector()

R> l <- length(gamma); for(i in 1:length(frc))mu[i] <- res[[(i-1)*l+1]]$mu

The performances of this method for estimating a meta-model are evaluated by considering a third dataset $(m(X^{third}_{i}),X_{i}^{third})$, $i=1,...,N$, with $N=1000$. The global prediction error is calculated as follows:

Let $\widehat{f}_{r}(.)$ be the best RKHS meta-model obtained in the repetition $r$, $r=1,...,N_{r}$, then

[TABLE]

The values of $GPE$ obtained for different kernels and values of $n$ are given in Table 2.

As expected the value of $GPE$ decreases as $n$ increases. The lowest values of $GPE$ are obtained when using the "matern" kernel.

In order to sum up the behaviour of the procedure for estimating the Sobol indices, we consider the mean square error (MSE) criterion obtained by $\sum_{v}(\sum_{r=1}^{N_{r}}(\widehat{S}_{v,r}-S_{v})^{2}/{N_{r}})$, where for each group $v$, $S_{v}$ denotes the true values of the Sobol indices, and $\widehat{S}_{v,r}$ is the empirical Sobol indices of the best RKHS meta-model in repetition $r$.

The obtained values of MSE for different kernels and values of $n$ are given in Table 3.

As expected, the values of MSE are smaller for larger values of $n$. The smallest values are obtained when using "matern" kernel.

The means of the empirical Sobol indices of the best RKHS meta-models through all repetitions for $n=200$ and "matern" kernel are displayed in Table 4.

It appears that the estimated Sobol indices are close to the true ones, nevertheless, they are overestimated for the main effects, i.e. groups $v\in\{\{1\},\{2\},\{3\}\}$, and underestimated for the interactions of order two and three, i.e. groups $v\in\{\{1,2\},\{1,3\},\{2,3\},\{1,2,3\}\}$.

Note that the strategy of choosing the tuning parameters is based on the minimization of the prediction error of the estimated meta-model, which may not minimize the error of estimating the Sobol indices.

Taking into account the results obtained for this Example 1, the calculations in the rest of the examples is done using only the "matern" kernel.

Example 2**.**

A time-efficient strategy to obtain the "optimal" tuning parameters when dealing with large data sets:

A dataset of $n$ points maximinLHS over $[0,1]^{d}$ with $n=1000$ and $d=10$ is generated. First, we use functions calc$\_$Kv and mu$\_$max to compute the eigenvalues and eigenvectors of the positive definite matrices $K_{v}$, and the value of $\mu_{max}$, respectively:

R> kernel <- "matern"; Dmax <- 3

R> Kv <- calc_Kv(X, kernel, Dmax, TRUE, TRUE)

R> mumax <- mu_max(Y, Kv$kv)

Then, we consider the two following steps:

Set $\gamma=0$ and, $\mu_{(1:9)}=\mu_{max}/(\sqrt{n}\times 2^{(2:10)}).$ Calculate the RKHS meta-models associated with the values of $\mu_{g}=\mu\times\sqrt{n}$ by using the function RKHSgrplasso. Gather the obtained RKHS meta-models in a list, "res$\_$g" (while this job could be done with the function RKHSMetMod by setting $\gamma=0$, in this example, we use the function RKHSgrplasso in order to avoid the re-calculation of $K_{v}$’s at the next step). Thereafter, for each estimator in res$\_$g, the prediction error is calculated by considering a new dataset and using the function PredErr. The value of $\mu$ with the smallest error of prediction in this step is denoted by $\mu_{i}$. Let us implement this step:

For a grid of values of $\mu_{g}$, a sequence of the RKHS meta-models are calculated and gathered in the "res$\_$g" list:

R> frc <- c(4,8,16,32,64,128,256,512,1024)

R> mu$\_$g <- mumax/frc

R> res$\_$g <- list();resg <- list()

R> for(i in 1:length(mu$\_$g))

R> resg[[i]] <- RKHSgrplasso(Y,Kv, mu$\_$g[i],1000,FALSE)

R> res$\_$g[[i]] <- list("mu$\_$g"=mu$\_$g,"gamma"=0,"MetaModel"=resg[[i]])

R>

Output res$\_$g contains nine RKHS meta-models and they are evaluated using a testing dataset:

R> gamma <- c(0); Err$\_$g <- PredErr(X, XT, YT, mu$\_$g, gamma, res$\_$g, kernel, Dmax)

The prediction errors of the RKHS meta-models obtained in this step are displayed in Table 5.

It appears that the minimum prediction error corresponds to the solution of the RKHS group lasso algorithm with $\mu_{g}=0.041$, so $\mu_{i}=0.041/\sqrt{n}$.

• 2.

Choose a smaller grid of values of $\mu$, $(\mu_{(i-1)},\mu_{i},\mu_{(i+1)})$, and set a grid of values of $\gamma>0$. Calculate the RKHS meta-models associated with each pair of the tuning parameters $(\mu,\gamma)$ by the function pen$\_$MetMod. Calculate the prediction errors for the new sequence of the RKHS meta-models using the function PredErr. Compute the empirical Sobol indices for the best estimator. Let us go back to the implementation of the example and apply this step $2$:

The grid of the values of $\mu$ in this step is $(0.081,0.041,0.020)/\sqrt{n}.$ The RKHS meta-models associated with this grid of the values of $\mu$ are gathered in a new list "resgnew". We set $\gamma_{(1:4)}=(0.2,0.1,0.01,0.005)$, and we calculate the RKHS meta-models for this new grid of the values of $(\mu,\gamma)$ using pen$\_$MetMod function:

R> gamma <- c(0.2, 0.1, 0.01, 0.005); mu <- c(mu$\_$g[5], mu$\_$g[6], mu$\_$g[7])/sqrt(n)

R> resgnew <- list()

R> resgnew[[1]] <- resg[[5]]; resgnew[[2]] <- resg[[6]]; resgnew[[3]] <- resg[[7]]

R> res <- pen$\_$MetMod(Y, Kv, gamma, mu, resgnew, 0, 0)

The output "res" is a list of twelve RKHS meta-models. These meta-models are evaluated using a new dataset, and their prediction errors are displayed in Table 6.

The minimum prediction error is associated with the pair $(0.020/\sqrt{n},0.01)$, and the best RKHS meta-model is then $\widehat{f}_{(0.020/\sqrt{n},0.01)}$.

The performance of this procedure for estimating the Sobol indices is evaluated using the relative error (RE) defined as follows:

For each $v$, let $S_{v}$ be the true value of the Sobol indices displayed in Table 4 and $\widehat{S}_{v}$ be the estimated empirical Sobol indices. Then

[TABLE]

In Table 7 the estimated empirical Sobol indices, their sum, and the value of RE are displayed.

The obtained RE for each group $v$ is smaller than $1.64\%$, therefore, the estimated Sobol indices in this example are very close to the true values of the Sobol indices displayed in the first row of Table 4.

Example 3**.**

Dealing with larger datasets:

Two datasets of $n$ points maximinLHS over $[0,1]^{d}$ with $n\in\{2000,5000\}$ and $d=10$ are generated. In order to obtain one RKHS meta-model associated with one pair of the tuning parameters $(\mu,\gamma)$, the number of coefficients to be estimated is equal to $n\times$vMax$=n\times 175$. Table 8 gives the execution time for different functions used throughout the Examples 1-3. In all examples we used a cluster of computers with: 2 Intel Xeon E5-2690 processors (2.90GHz) and 96Gb Ram (6x16Gb of memory 1600MHz).

As we can see, the execution time increases fast as $n$ increases. In Figure 1 the plot of the logarithm of the time (in seconds) versus the logarithm of $n$ is displayed for the functions calc$\_$Kv, mu$\_$max, RKHSgrplasso and pen$\_$MetMod.

It appears that, the algorithms of these functions are of polynomial time $O(n^{\alpha})$ with $\alpha\backsimeq 3$ for the functions calc$\_$Kv and mu$\_$max, and $\alpha\backsimeq 2$ for the functions RKHSgrplasso and pen$\_$MetMod.

Taking into account the results obtained for the prediction error and the values of $(\widehat{\mu},\widehat{\gamma})$ in Example 2, in this example, only two values of the tuning parameter $\mu_{(1:12)}=\mu_{max}/(\sqrt{n}\times 2^{(7:8)})$, and one value of the tuning parameter $\gamma=0.01$ are considered. The RKHS meta-models associated with the pair of values $(\mu_{i},\gamma)$, $i=1,2$ are estimated using the RKHSMetMod function:

R> kernel <- "matern"; Dmax <- 3

R> gamma <- c(0.01); frc <- c(128,256)

R> res <- RKHSMetMod(Y, X, kernel, Dmax, gamma, frc, FALSE)

The prediction error and the empirical Sobol indices are then calculated for the obtained meta-models using the functions PredErr and SI$\_$emp:

R> mu <- vector(); mu[1] <- res[[1]]$mu; mu[2] <- res[[2]]$mu

R> Err <- PredErr(X, XT, YT, mu, gamma, res, kernel, Dmax)

R> SI <- SI$\_$emp(res, NULL)

The result of the prediction errors associated with the obtained estimators for two different values of $n$ are displayed in Table 9.

For $n$ equal to $5000$, we obtained smaller values of the prediction error, so as expected, the prediction quality improves by increasing the number of the observations $n$. Table 10 gives the estimated empirical Sobol indices as well as their sum and the values of RE (see Equation (18)).

Comparing the values of RE, we can see that the empirical Sobol indices are better estimated for n equal to $5000$, so as expected, the estimation of the Sobol indices is better for larger values of $n$.

In Figure 2 the result of the estimation quality and the Sobol indices for the dataset with $n$ equal to $5000$, $d$ equal to $10$, and $(\mu_{2},\gamma)$ are displayed.

The line $y=x$ in red crosses the cloud of points as long as the values of the g-function are smaller than three. When the values of the g-function are greater than three, the estimator $\widehat{f}$ tends to underestimate the g-function.

4.2 Comparison examples

This section includes two examples. In the first example we reproduce an example from paper [\citeauthoryearGu, Palomo, and BergerGu et al.2019] and compare the prediction quality of the RKHS meta-model with the GP (Kriging) based meta-models from the RobusGaSP [\citeauthoryearGu, Palomo, and BergerGu et al.2019] and DiceKriging [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012] packages.

The objective is to evaluate the quality of the RKHS estimated meta-model and to compare it with methods recently proposed for approximating complex models. In the first example we consider one-dimensional model and focuse on the comparison between the true model and the estimated meta-model. In the second example we reproduce an example from paper [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012] which allows us to compare the prediction quality of the RKHS meta-model with the Kriging based meta-model from DiceKriging package, as well as the estimation quality of the Sobol indices in our package with the well-known package sensitivity.

For the sake of comparison between the three methods, the meta-models are calculated using the same experimental design and outputs and the same kernel function available in three packages is used. However, in packages RobustGaSP and DiceKriging the range parameter $r$ (see Table 1) in the kernel function is estimated by marginal posterior modes with the robust parameterization and by MLE with upper and lower bounds, respectively, while it is assumed to be fixed and set as $\sqrt{3}/2$ in the RKHSMetaMod package.

Example 4**.**

"The modified sine wave function" [\citeauthoryearGu, Palomo, and BergerGu et al.2019].

We consider the $1$-dimensional modified sine wave function defined over $[0,1]$ by:

[TABLE]

The same experimental design as described in [\citeauthoryearGu, Palomo, and BergerGu et al.2019] is considered: the design matrix $X$ is a sequence of $12$ equally spaced points on $[0,1]$, and the response variable $Y$ is calculated as $Y=m(X)$:

R> X <- as.matrix(seq(0,1,1/11)); Y <- sinewave(X)

where sinewave function is defined in [\citeauthoryearGu, Palomo, and BergerGu et al.2019]. We build the GP based meta-models by the RobustGaSP and the DiceKriging packages using the constant mean function and kernel Matérn $3/2$:

R> library(RobustGaSP)

R> res.rgasp <- rgasp(design=X, response=Y, kernel$\_$type="matern$\_$3$\_$2")

R> library(DiceKriging)

R> res.km <- km(design=X, response=Y, covtype="matern3$\_$2")

As $d=1$, we have $Dmax=1$. We consider the grid of values of $\mu_{(1:9)}=\mu_{max}/(\sqrt{n}\times 2^{(2:10)})$ and $\gamma_{(1:5)}=(0.2,0.1,0.01,0.005,0)$. The RKHS meta-models associated with the pair of values $(\mu_{i},\gamma_{j})$, $i=1,\cdots,9$, $j=1,\cdots,5$ are estimated using the RKHSMetMod function:

R> kernel <- "matern"; Dmax <- 1

R> gamma <- c(0.2, 0.1, 0.01, 0.005,0)

R> frc <- c(4,8,16,32,64,128,256,512,1024)

R> res <- RKHSMetMod(Y, X, kernel, Dmax, gamma, frc, FALSE)

Given a testing dataset $(XT,YT)$, the prediction errors associated with the obtained RKHS meta-models are calculated using the PredErr function, and the best RKHS meta-model is chosen to be the estimator of the model $m(X)$:

R> XT <- as.matrix(seq(0,1,1/11)); YT <- sinewave(XT)

R> Err <- PredErr(X, XT, YT, mu, gamma, res, kernel, Dmax)

To compare these three estimators in terms of the prediction quality, we perform prediction on $100$ test points, equally spaced in $[0,1]$:

R> predict$\_$X <- as.matrix(seq(0,1,1/99))

R> #prediction with the GP based meta-models:

R> rgasp.predict <- predict(res.rgasp, predict$\_$X)

R> km.predict <- predict(res.km, predict$\_$X, type=’UK’)

R> #prediction with the best RKHS meta-model:

R> res.predict <- prediction(X, predict$\_$X, kernel, Dmax, res, Err)

The prediction results are plotted in Figure 3. The black circles that correspond to the prediction from the RKHSMetMod package are closer to the real output than the green and the blue circles corresponding to the predictive means from the RobustGaSP and DiceKriging packages.

The meta-model results are plotted in Figure 4. The prediction from the RKHSMetaMod package plotted as the black curve is much more accurate as an estimate of the true function (plotted in red) than the predictive mean from the RobustGaSP and DiceKriging packages plotted as the blue and green curves, respectively. As already noted by [\citeauthoryearGu, Palomo, and BergerGu et al.2019], for that sine wave example, the meta-model from the DiceKriging package "degenerates to the fitted mean with spikes at the design points".

Example 5**.**

"A standard SA 8-dimensional example" [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012].

We consider the $8$-dimensional g-function of Sobol implemented in the package sensitivity: the function $m(X)$ as defined in Equation (17) with coefficients $c_{1}=0,\>c_{2}=1,\>c_{3}=4.5,\>c_{4}=9,\>(c_{a})_{a=5,6,7,8}=99$. With these values of coefficients $c_{a}$, the variables $X_{1}$, $X_{2}$, $X_{3}$ and $X_{4}$ explain $99.96\%$ of the variance of the function $m(X)$ (see Table 11).

We consider the same experimental design as described in [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012]: the design matrices $X$ and $XT$ are $80$-point optimal Latin Hypercube Samples of the input variables generated by the optimumLHS function of package lhs, and the response variables $Y$ and $YT$ are calculated as $Y=m(X)$, and $YT=m(XT)$ using sobol.fun function of the package sensitivity:

R> n <- 80; d <- 8

R> library(lhs); X <- optimumLHS(n, d); XT <- optimumLHS(n, d)

R> library(sensitivity); Y <- sobol.fun(X); YT <- sobol.fun(XT)

Let us first consider the RKHS meta-model method. We set Dmax$=3$, and we consider the grid of values of $\mu_{(1:9)}=\mu_{max}/(\sqrt{n}\times 2^{(2:10)})$, and $\gamma_{(1:5)}=(0.2,0.1,0.01,0.005,0)$. The RKHS meta-models associated with the pair of values $(\mu_{i},\gamma_{j})$, $i=1,\cdots,9$, $j=1,\cdots,5$ are estimated using the RKHSMetMod function:

R> kernel <- "matern"; Dmax <- 3

R> gamma <- c(0.2, 0.1, 0.01, 0.005,0)

R> frc <- c(4,8,16,32,64,128,256,512,1024)

R> res <- RKHSMetMod(Y, X, kernel, Dmax, gamma, frc, FALSE)

Given the testing dataset $(XT,YT)$, the prediction errors associated with the obtained RKHS meta-models are calculated using PredErr function, and the best RKHS meta-model is chosen to be the estimator of the model $m(X)$. Finally, the Sobol indices are computed for the best RKHS meta-model using the function SI$\_$emp:

R> Err <- PredErr(X, XT, YT, mu, gamma, res, kernel, Dmax)

R> SI <- SI$\_$emp(res, Err)

Secondly, let us build the GP based meta-model. We use the km function of the package DiceKriging with the constant mean function and kernel Matérn $3/2$:

R> library(DiceKriging)

R> res.km <- km(design = X, response = Y, covtype = "matern3$\_$2")

The Sobol indices associated with the estimated GP based meta-model are calculated using fast99 function of the package sensitivity:

R> SI.km <- fast99(model = kriging.mean, factors = d, n = 1000,

• q = "qunif", q.arg = list(min = 0, max = 1), m = res.km)

where kriging.mean function is defined in [\citeauthoryearRoustant, Ginsbourger, and DevilleRoustant et al.2012].

The result of the estimation with the best RKHS meta-model and the Kriging based meta-model is drawn in Figure 5.

The black circles that correspond to the best RKHS meta-model are closer to the real output than the blue circles corresponding to the GP based meta-model from the DiceKriging package. Another way to evaluate the prediction quality of the estimated meta-models is to consider the mean square error of the fitted meta-model computed by $\sum_{i=1}^{80}(m(X_{i})-\hat{f}(X_{i}))^{2}/80$. We obtained $3.96\%$ and $0.07\%$ for the Kriging based meta-model and the RKHS meta-model, respectively, which confirms the good behavior of the RKHS meta-model.

The estimated Sobol indices associated with the RKHS meta-model and the Kriging based meta-model are given in Table 11.

As shown, with RKHS meta-model, we obtained non-zero values for the interactions of order two. Concerning the main effects, excepting the first one, the estimated Sobol indices with the RKHS meta-model are closer to the true ones. However, the interactions of order three are ignored by both meta-models. For a general comparison of the estimation quality of the Sobol indices, one may consider the criterion RE defined in Equation (18), which is equal to $7.95$ for the Kriging based meta-model, and $5.59$ for the RKHS meta-model. Comparing the values of RE, we can point out that the Sobol indices are better estimated with the RKHS meta-model in that model.

5 Summary and discussion

In this paper, we proposed an R package, called RKHSMetaMod, that estimates a meta-model of a complex model $m$. This meta-model belongs to a reproducing kernel Hilbert space constructed as a direct sum of Hilbert spaces [\citeauthoryearDurrande, Ginsbourger, Roustant, and CarraroDurrande et al.2013]. The estimation of the meta-model is carried out via a penalized least-squares minimization allowing both to select and estimate the terms in the Hoeffding decomposition, and therefore, to select the Sobol indices that are non-zero and estimate them [\citeauthoryearHuet and TaupinHuet and Taupin2017]. This procedure makes it possible to estimate the Sobol indices of high order, a point known to be difficult in practice. Using the convex optimization tools, RKHSMetaMod package implements two optimization algorithms: the minimization of the RKHS ridge group sparse criterion (13) and the RKHS group lasso criterion (14). Both of these algorithms rely on the Gram matrices $K_{v}$, $v\in\mathcal{P}$ and their positive definiteness. Currently, the package considers only uniformly distributed input variables. If one is interested by another distribution of the input variables, it suffices to modify the calculation of the kernels $k_{0a}$, $a=1,...,d$ (see Equation (19)) in the function calc$\_$Kv of this package (see Remark 3). The available kernels in the RKHSMetaMod package are: Gaussian kernel (with the fixed range parameter $r=1/2$), Matérn kernel (with the fixed range parameter $r=\sqrt{3}/2$), Brownian kernel, quadratic kernel and linear kernel (see Table 1). With regard to the problem being under study, one may consider other kernels or kernels with different values of the range parameter $r$ and add them easily to the list of the kernels in the calc$\_$Kv function. For the large values of $n$ and $d$ the calculation and storage of the eigenvalues and the eigenvectors of all the Gram matrices $K_{v}$, $v\in\mathcal{P}$ require a lot of time and a very large amount of memory. In order to optimize the execution time and also the storage memory, except for a function that is written in R, all of the functions of RKHSMetaMod package are written using the efficient C++ libraries through RcppEigen and RcppGSL packages. These functions are then interfaced with the R environment in order to propose a user friendly package.

6 Supplementary materials

6.1 Appendix

Preliminary 1**.**

*For $F(x)=\|Ax\|$, where $A$ is a symmetric matrix that does not depend on x, we have: $\partial F(x)=\{{A^{2}x/\|Ax\|}\}$ if $x\neq 0$, and $\partial F(x)=\{w\in\mathbb{R}^{n},\quad\|A^{-1}w\|\leq 1\}$ if $x=0$. *

Preliminary 2**.**

Let $F:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a convex function. We have the following first order optimality condition:

[TABLE]

This results from the fact that $F(y)\geq F(\widehat{x})+<0,y-\widehat{x}>$ for all $y\in\mathbb{R}^{n}$ in both cases [\citeauthoryearGiraudGiraud2014].

6.1.1 RKHS construction

We begin this Section with a brief introduction to the RKHS. Let $\mathcal{H}$ be a Hilbert space of real valued functions on a set $\mathcal{X}$. The space $\mathcal{H}$ is a RKHS if for all $X\in\mathcal{X}$ the evaluation functionals $L_{X}:\mathcal{H}(f)\rightarrow\mathbb{R}(f(X))$ are continuous. The Riesz representation Theorem ensures the existence of an unique element $k_{X}(.)$ in $\mathcal{H}$ verifying the property that $\forall X\in\mathcal{X},\>\forall f\in\mathcal{H},\>f(X)=L_{X}(f)=\langle f,k_{X}\rangle_{\mathcal{H}},$ where $\langle.,.\rangle_{\mathcal{H}}$ denotes the inner product in $\mathcal{H}$. It follows that for all $X$ ,$X^{\prime}$ in $\mathcal{X}$, and $k_{X}(.)$, $k_{X^{\prime}}(.)$ in $\mathcal{H}$, we have $k_{X}(X^{\prime})=L_{X^{\prime}}(k_{X})=\langle k_{X},k_{X^{\prime}}\rangle_{\mathcal{H}}.$ This allows to define the reproducing kernel of $\mathcal{H}$ as $k:\mathcal{X}\times\mathcal{X}((X,X^{\prime}))\rightarrow\mathbb{R}(k_{X}(X^{\prime})).$ The reproducing kernel $k(X,X^{\prime})$ is positive definite since it is symmetric, and for any $n\in\mathbb{N}$, $\{X_{i}\}_{i=1}^{n}\in\mathcal{X}$ and $\{c_{i}\}_{i=1}^{n}\in\mathbb{R}$, we have:

[TABLE]

For more background on RKHS, we refer to various standard references such as [\citeauthoryearAronszajnAronszajn1950], [\citeauthoryearSaitohSaitoh1988], and [\citeauthoryearBerlinet and Thomas-AgnanBerlinet and Thomas-Agnan2003].

In this work, the idea is to construct a RKHS $\mathcal{H}$ such that any function $f$ in $\mathcal{H}$ is decomposed as its Hoeffding decomposition, and therefore, any function $f$ in $\mathcal{H}$ is a candidate to approximate the Hoeffding decomposition of $m$. To do so, the method of [\citeauthoryearDurrande, Ginsbourger, Roustant, and CarraroDurrande et al.2013] as described below is used.

Let $\mathcal{X}=\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{d}$ be a subset of $\mathbb{R}^{d}$. For each $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:

• (i)

$k_{a}:\mathcal{X}_{a}\times\mathcal{X}_{a}\rightarrow\mathbb{R}$ is $P_{a}\otimes P_{a}$ measurable,

• (ii)

$E_{X_{a}}\sqrt{k_{a}(X_{a},X_{a})}<\infty$.

The property (ii) depends on the kernel $k_{a}$, $a=1,...,d$ and the distribution of $X_{a}$, $a=1,...,d$. It is not very restrictive since it is satisfied, for example, for any bounded kernel.

The RKHS $\mathcal{H}_{a}$ can be decomposed as a sum of two orthogonal sub-RKHS, $\mathcal{H}_{a}=\mathcal{H}_{0a}\stackrel{{\scriptstyle\perp}}{{\oplus}}\mathcal{H}_{1a},$ where $\mathcal{H}_{0a}$ is the RKHS of zero mean functions, $\mathcal{H}_{0a}=\{f_{a}\in\mathcal{H}_{a}:\>E_{X_{a}}(f_{a}(X_{a}))=0\},$ and $\mathcal{H}_{1a}$ is the RKHS of constant functions, $\mathcal{H}_{1a}=\{f_{a}\in\mathcal{H}_{a}:\>f_{a}(X_{a})=C\}.$ The kernel $k_{0a}$ associated with the RKHS $\mathcal{H}_{0a}$ is defined by:

[TABLE]

Let $k_{v}(X_{v},X^{\prime}_{v})=\prod_{a\in v}k_{0a}(X_{a},X^{\prime}_{a}),$ then the ANOVA kernel $k(.,.)$ is defined as follows:

[TABLE]

For $\mathcal{H}_{v}$ being the RKHS associated with the kernel $k_{v}$, the RKHS associated with the ANOVA kernel is then defined by $\mathcal{H}=\prod_{a=1}^{d}(\mathbbm{1}\stackrel{{\scriptstyle\perp}}{{\oplus}}\mathcal{H}_{0a})=\mathbbm{1}+\sum_{v\in\mathcal{P}}\mathcal{H}_{v},$ where $\perp$ denotes the $L^{2}$ inner product. According to this construction, any function $f\in\mathcal{H}$ satisfies decomposition (7), $f(X)=\langle f,k(X,.)\rangle_{\mathcal{H}}=f_{0}+\sum_{v\in\mathcal{P}}f_{v}(X_{v}),$ which is the Hoeffding decomposition of $f$.

The regularity properties of the RKHS $\mathcal{H}$ constructed as described above, depend on the set of the kernels ($k_{a}$, $a=1,...,d$). This method allows us to choose different approximation spaces independently of the distribution of the input variables $X_{1},...,X_{d}$, by choosing different sets of kernels. While, as mentioned earlier, in the meta-modelling approach based on the polynomial Chaos expansion, according to the distribution of the input variables $X_{1},...,X_{d}$, a unique family of orthonormal polynomials $\{\phi_{j}\}_{j=0}^{\infty}$ is determined. Here, the distribution of the components of $X$ occurs only for the orthogonalization of the spaces $\mathcal{H}_{v}$, $v\in\mathcal{P}$, and not in the choice of the RKHS, under the condition that properties (i) and (ii) are satisfied. This is one of the main advantages of this method compared to the method based on the truncated polynomial Chaos expansion where the smoothness of the approximation is handled only by the choice of the truncation [\citeauthoryearBlatman and SudretBlatman and Sudret2011].

6.1.2 RKHS group lasso algorithm

We consider the minimization of the RKHS group lasso criterion given by,

[TABLE]

We begin with the constant term $f_{0}$. The ordinary first derivative of the function $C_{g}(f_{0},\theta)$ at $f_{0}$ is equal to:

[TABLE]

and therefore,

[TABLE]

where $(K_{v}\theta_{v})_{i}$ denotes the i-th component of $K_{v}\theta_{v}$.

The next step is to calculate $\widehat{\theta}=\text{argmin}_{\theta\in\mathbb{R}^{n\times|\mathcal{P}|}}C_{g}(f_{0},\theta).$ Since $C_{g}(f_{0},\theta)$ is convex and separable, we use a block coordinate descent algorithm, group $v$ by group $v$. In the following, we fix a group $v$, and we find the minimizer of $C_{g}(f_{0},\theta)$ with respect to $\theta_{v}$ for the given values of $f_{0}$ and $\theta_{w}$, $w\neq v$. Set

[TABLE]

where

[TABLE]

We aim to minimize $C_{g,v}(f_{0},\theta_{v})$ with respect to $\theta_{v}$. Let $\partial C_{g,v}$ be the sub-differential of $C_{g,v}(f_{0},\theta_{v})$ with respect to $\theta_{v}$:

[TABLE]

The first order optimality condition (see Preliminary (2)) ensures the existence of $\widehat{t}_{v}\in\partial\|K_{v}^{1/2}\theta_{v}\|$ fulfilling,

[TABLE]

Using the sub-differential definition (see Preliminary 1), we obtain:

[TABLE]

and,

[TABLE]

Let $\widehat{\theta}_{v}$ be the minimizer of $C_{g,v}$. The sub-differential equations above give the two following cases:

Case 1. If $\widehat{\theta}_{v}=0$, then there exists $\widehat{t}_{v}\in\mathbb{R}^{n}$ such that $\|K_{v}^{-1/2}\widehat{t}_{v}\|\leq 1$ and it fulfils Equation (22), $2K_{v}R_{v}=\sqrt{n}\mu_{g}\widehat{t}_{v}$. Therefore, the necessary and sufficient condition for which the solution $\widehat{\theta}_{v}=0$ is the optimal one is $2\|K_{v}^{1/2}R_{v}\|/\sqrt{n}\leq\mu_{g}.$

Case 2. If $\widehat{\theta}_{v}\neq 0$, then $\widehat{t}_{v}=K_{v}\widehat{\theta}_{v}/\|K_{v}^{1/2}\widehat{\theta}_{v}\|$ and it fulfils Equation (22),

[TABLE]

We obtain then,

[TABLE]

Since $\widehat{\theta}_{v}$ appears in both sides of the Equation (23), a numerical procedure is needed:

Proposition 1**.**

For $\rho>0$ let $\theta(\rho)=(K_{v}+\rho I_{n})^{-1}R_{v}$. There exists a non-zero solution to Equation (23) if and only if there exists $\rho>0$ such that

[TABLE]

Then $\widehat{\theta}_{v}=\theta(\rho)$.

Proof

If there exists a non-zero solution to Equation (23), then $\|K_{v}^{1/2}\widehat{\theta}_{v}\|\neq 0$ since $K_{v}$ is positive definite. Take $\rho={\sqrt{n}\mu_{g}/2\|K_{v}^{1/2}\widehat{\theta}_{v}\|},$ then $\theta(\rho)=(K_{v}+\frac{\sqrt{n}\mu_{g}}{2\|K_{v}^{1/2}\widehat{\theta}_{v}\|}I_{n})^{-1}R_{v}=\widehat{\theta}_{v},$ and, for such $\rho$ Equation (24) is satisfied. Conversely, if there exists $\rho>0$ such that Equation (24) is satisfied, then $\|K_{v}^{1/2}\theta(\rho)\|\neq 0$ and $\rho=\frac{\sqrt{n}\mu_{g}}{2\|K_{v}^{1/2}\theta(\rho)\|}.$ Therefore, $\theta(\rho)=(K_{v}+\frac{\sqrt{n}\mu_{g}}{2\|K_{v}^{1/2}\theta(\rho)\|}I_{n})^{-1}R_{v},$ which is Equation (23) calculated in $\widehat{\theta}_{v}=\theta(\rho)$. $\Box$

Remark 5**.**

Define $y(\rho)=2\rho\|K_{v}^{1/2}\theta(\rho)\|-\sqrt{n}\mu_{g}$ with $\theta(\rho)=(K_{v}+\rho I_{n})^{-1}R_{v}$, then $y(\rho)=0$ has a unique solution, denoted $\widehat{\rho}$, which leads to calculate $\widehat{\theta}(\widehat{\rho})$.

Proof

For $\rho=0$ we have $y(0)=-\sqrt{n}\mu_{g}<0$, since $\mu_{g}>0$; and for $\rho\rightarrow+\infty$ we have $y(\rho)>0$, since $\|K_{v}^{1/2}(\frac{K_{v}}{\rho}+I_{n})^{-1}R_{v}\|\rightarrow\|K_{v}^{1/2}R_{v}\|$ and $\|2K_{v}^{1/2}R_{v}\|>\sqrt{n}\mu_{g}$.

Moreover, we have $y(\rho)=2\|(\frac{I_{n}}{\rho}+k_{v}^{-1})^{-1}k_{v}^{-1/2}R_{v}\|-\sqrt{n}\mu_{g}$, which is equal to $2(X^{T}A^{-2}X)^{1/2}-\sqrt{n}\mu_{g}$ with $A=(I_{n}/\rho+k_{v}^{-1})$ and $X=k_{v}^{-1/2}R_{v}$. The first derivative of $y(\rho)$ in $\rho$ is obtained by $\frac{\partial y(\rho)}{\partial\rho}=(X^{T}A^{-2}X)^{-1/2}\frac{\partial(X^{T}A^{-2}X)}{\partial\rho}$. Finally, by simple calculations we get,

[TABLE]

Therefore $y(\rho)$ is an increasing function of $\rho$, and the proof is complete. $\Box$

In order to calculate $\rho$ and so $\widehat{\theta}_{v}=\theta(\rho)$, we use Algorithm 4 which is a part of the RKHS group lasso Algorithm 1 when $\widehat{\theta}_{v}\neq 0$.

6.1.3 Computational cost

The complexity for the matrices $K_{v}$, $v\in\mathcal{P}$ is equal to $n^{3}$, which is given by the singular value decomposition to get eigenvalues and eigenvectors of each $K_{v}$. Supposing that the matrices $K_{v}$, $v\in\mathcal{P}$ was first created and are already stored, the complexity for the constant term $f_{0}$ is given by the second term in equation (20), which is equal to $n^{2}$. Given $\widehat{\rho}$, the complexity for $\widehat{\theta}_{v}$, $v\in\mathcal{P}$, is given by the backsolving of $(K_{v}+\widehat{\rho}I_{n})\widehat{\theta}_{v}=R_{v}$ to get $\widehat{\theta}_{v}$, which is equal to $n^{2}$. The computation of $\widehat{\rho}$ is done using Brent-Dekker method implemented as function gsl_root_fsolver_brent from [\citeauthoryearGalassi, Davies, Theiler, Gough, Jungman, Alken, Booth, Rossi, and UlerichGalassi et al.2018]. This method combines an interpolation strategy with the bisection algorithm and takes $O(m)$ iterations to converge, where $m$ is the number of steps that the bisection algorithm would take.

6.1.4 RKHS ridge group sparse algorithm

We consider the minimization of the RKHS ridge group sparse criterion:

[TABLE]

The constant term $f_{0}$ is estimated as in the RKHS group lasso algorithm. In order to calculate $\widehat{\theta}=\text{argmin}_{\theta\in\mathbb{R}^{n\times|\mathcal{P}|}}C(f_{0},\theta)$, we use once again the block coordinate descent algorithm group $v$ by group $v$. In the following, we fix a group $v$, and we find the minimizer of $C(f_{0},\theta)$ with respect to $\theta_{v}$ for given values of $f_{0}$ and $\theta_{w}$, $w\neq v$. We aim at minimizing with respect to $\theta_{v}$,

[TABLE]

where $R_{v}$ is defined by (21).

Let $\partial C_{v}$ be the sub-differential of $C_{v}(f_{0},\theta_{v})$ with respect to $\theta_{v}$,

[TABLE]

According to the first order optimality condition (see Preliminary 2), we know that there exists $\widehat{s}_{v}\in\partial\|K_{v}\theta_{v}\|$ and $\widehat{t}_{v}\in\partial\|K_{v}^{1/2}\theta_{v}\|$ such that,

[TABLE]

The sub-differential definition (see Preliminary 1) gives:

[TABLE]

and,

[TABLE]

Let $\widehat{\theta}_{v}$ be the minimizer of the $C_{v}(f_{0},\theta_{v})$. Using the sub-differential equations above, the estimator $\widehat{\theta}_{v}$, $v\in\mathcal{P}$ is obtained following the two cases below:

Case 1. If $\widehat{\theta}_{v}=0$, then there exists $\widehat{s}_{v}\in\mathbb{R}^{n}$ such that $\|K_{v}^{-1}\widehat{s}_{v}\|\leq 1$ and it fulfils Equation (25), $2K_{v}R_{v}-n\mu\widehat{t}_{v}=\sqrt{n}\gamma\widehat{s}_{v},$ with $\widehat{t}_{v}\in\mathbb{R}^{n}$, $\|K_{v}^{-1/2}\widehat{t}_{v}\|\leq 1$. Set $J(\widehat{t}_{v})=\|2R_{v}-n\mu K_{v}^{-1}\widehat{t}_{v}\|,$ and, $J^{*}=\text{argmin}_{\widehat{t}_{v}\in\mathbb{R}^{n}}\{J(\widehat{t}_{v}),\text{ such that }\|K_{v}^{-1/2}\widehat{t}_{v}\|\leq 1\}.$ Then the solution to Equation (25) is zero if and only if $J^{*}\leq\gamma$.

Case 2. If $\widehat{\theta}_{v}\neq 0$, then we have $\widehat{s}_{v}=K_{v}^{2}\widehat{\theta}_{v}/\|K_{v}\widehat{\theta}_{v}\|$, and $\widehat{t}_{v}=K_{v}\widehat{\theta}_{v}/\|K_{v}^{1/2}\widehat{\theta}_{v}\|$ fulfilling Equation (25), $2K_{v}(R_{v}-K_{v}\widehat{\theta}_{v})=\sqrt{n}\gamma\frac{K_{v}^{2}\widehat{\theta}_{v}}{\|K_{v}\widehat{\theta}_{v}\|_{2}}+n\mu\frac{K_{v}\widehat{\theta}_{v}}{\|K_{v}^{1/2}\widehat{\theta}_{v}\|},$ that is,

[TABLE]

In this case the calculation of $\widehat{\theta}_{v}$ needs a numerical algorithm.

Proposition 2**.**

(Proposition 8.4 in [\citeauthoryearHuet and TaupinHuet and Taupin2017]) For $\rho_{1},\rho_{2}>0$, let $\theta(\rho_{1},\rho_{2})=(K_{v}+\rho_{1}K_{v}+\rho_{2}I_{n})^{-1}R_{v}$. If $\mu>0$, there exists a non zero solution to Equation (26) if and only if there exists $\rho_{1},\rho_{2}>0$ such that $\gamma=\frac{2\rho_{1}}{\sqrt{n}}\|K_{v}\theta(\rho_{1},\rho_{2})\|,\>\mbox{and}\>\mu=\frac{2\rho_{2}}{n}\|K^{1/2}_{v}\theta(\rho_{1},\rho_{2})\|.$ Then $\widehat{\theta}_{v}=\theta(\rho_{1},\rho_{2})$.

Proof

The proof is given in [\citeauthoryearHuet and TaupinHuet and Taupin2017].

6.1.5 Computational cost

The complexity for the matrices $K_{v}$, $v\in\mathcal{P}$ and the constant term $f_{0}$ is the same as for RKHS group lasso algorithm. Given $\widehat{\rho}_{1},\widehat{\rho}_{2}$, the complexity for $\widehat{\theta}_{v}$, $v\in\mathcal{P}$, is given by the backsolving of $(K_{v}+\rho_{1}K_{v}+\rho_{2}I_{n})\widehat{\theta}_{v}=R_{v}$ to get $\widehat{\theta}_{v}$, which is equal to $n^{2}$. The computation of $\widehat{\rho}_{1},\widehat{\rho}_{2}$, is insured using a combination of three methods: first we implement a modified version of Newton’s method, if it does not achieve the convergence second we implement a version of the Hybrid algorithm and if it does not achieve the convergence, third we implement a version of the discrete Newton algorithm called Broyden algorithm. These methods are implemented as functions gsl_multiroot_fdfsolver_gnewton, gsl_multiroot_fsolver_hybrids, and gsl_multiroot_fsolver_broyden from [\citeauthoryearGalassi, Davies, Theiler, Gough, Jungman, Alken, Booth, Rossi, and UlerichGalassi et al.2018], respectively. These methods converge fast in general, and for $n$ big enough their complexity is dominated by $n^{2}$.

6.2 Overview of the RKHSMetaMod functions

In the R environment, one can install and load the RKHSMetaMod package by using the following commands:

Example 6**.**

install.packages("RKHSMetaMod") library("RKHSMetaMod")

The optimization problems in this package are solved using block coordinate descent algorithm which requires various computational algorithms including generalized Newton, Broyden and Hybrid methods. In order to gain the efficiency in terms of the calculation time and be able to deal with high dimensional problems, the computationally efficient tools of C++ packages Eigen [\citeauthoryearGuennebaud, Jacob, et al.Guennebaud et al.2010] and GSL [\citeauthoryearGalassi, Davies, Theiler, Gough, Jungman, Alken, Booth, Rossi, and UlerichGalassi et al.2018] via RcppEigen [\citeauthoryearBates and EddelbuettelBates and Eddelbuettel2013] and RcppGSL [\citeauthoryearEddelbuettel and FrancoisEddelbuettel and Francois2019] packages are used in the RKHSMetaMod package. For different examples of usage of RcppEigen and RcppGSl functions see the work by [\citeauthoryearEddelbuettelEddelbuettel2013].

The complete documentation of RKHSMetaMod package is available from CRAN [\citeauthoryearKamariKamari2019]. Here, a brief documentation of some of its main and companion functions is presented in the next two Sections.

6.2.1 Main RKHSMetaMod functions

Let us begin by introducing some notations. For a given Dmax$\in\mathbb{N}$, let $\mathcal{P}_{\text{Dmax}}$ be the set of parts of $\{1,...,d\}$ with dimension $1$ to Dmax. The cardinal of $\mathcal{P}_{\text{Dmax}}$ is denoted by $\text{vMax}=\sum_{j=1}^{\text{Dmax}}\binom{d}{j}.$

RKHSMetMod function:

For a given value of Dmax and a chosen kernel (see Table 1), this function calculates the Gram matrices $K_{v}$, $v\in\mathcal{P}_{\text{Dmax}}$, and produces a sequence of estimators $\widehat{f}$ associated with a given grid of values of tuning parameters $\mu,\gamma$, i.e. the solutions to the RKHS ridge group sparse (if $\gamma\neq 0$) or the RKHS group lasso problem (if $\gamma=0$). Table 12 gives a summary of all the input arguments of the RKHSMetMod function as well as the default values for non-mandatory arguments.

The RKHSMetMod function returns a list of $l$ components, with $l$ being equal to the number of pairs of the tuning parameters $(\mu,\gamma)$, i.e. $l=|\text{gamma}|\times|\text{frc}|$. Each component of the list is a list of three components "mu", "gamma" and "Meta-Model":

• •

mu: value of the tuning parameter $\mu$ if $\gamma>0$, or $\mu_{g}=\sqrt{n}\times\mu$ if $\gamma=0$.

• •

gamma: value of the tuning parameter $\gamma$.

• •

Meta-Model: a RKHS ridge group sparse or RKHS group lasso object associated with the tuning parameters mu and gamma. Table 13 gives a summary of all arguments of the output "Meta-Model" of RKHSMetMod function.

RKHSMetMod$\_$qmax function:

For a given value of Dmax and a chosen kernel (see Table 1), this function calculates the Gram matrices $K_{v}$, $v\in\mathcal{P}_{\text{Dmax}}$; determines $\mu$, reffered to as $\mu_{qmax}$, for which the number of groups in the support of the RKHS group lasso solution is equal to $qmax$; and produces a sequence of estimators $\widehat{f}$ associated with the tuning parameter $\mu_{qmax}$ and a grid of values of the tuning parameter $\gamma$. All the estimators $\widehat{f}$ produced by this function have at most $qmax$ groups in their support. This function has the following input arguments:

• $-$

$Y$, $X$, kernel, Dmax, gamma, verbose (see Table 12).

• $-$

qmax: integer, the maximum number of groups in the support of the obtained solution.

• $-$

rat: positive scalar, to restrict the minimum value of $\mu$ considered in Algorithm 3,

[TABLE]

where the value of $\mu_{\text{max}}$ is given by Equation (16) and is calculated inside the program.

• $-$

Num: integer, to restrict the number of different values of the tuning parameter $\mu$ to be evaluated in the RKHS group lasso algorithm until it achieves $\mu_{qmax}$. For example, if Num equals $1$, the program is implemented for three different values of $\mu\in[\mu_{\text{min}},\mu_{\text{max}})$:

[TABLE]

where $|\widehat{S}_{\widehat{f}(\mu_{1})_{\text{Group Lasso}}}|$ is the number of groups in the support of the solution of the RKHS group lasso problem, Algorithm 1, associated with $\mu_{1}$.

If Num$>1$, the path to cover the interval $[\mu_{\text{min}},\mu_{\text{max}})$ is detailed in Algorithm 3.

The RKHSMetMod$\_$qmax function returns a list of three components "mus", "qs", and "MetaModel":

• •

mus: vector of all values of $\mu_{i}$ in Algorithm 3.

• •

qs: vector with the same length as mus. Each element of the vector shows the number of groups in the support of the RKHS meta-model obtained by solving RKHS group lasso problem for an element in mus.

• •

MetaModel: list with the same length as the vector gamma. Each component of the list is a list of three components "mu", "gamma" and "Meta-Model":

• –

mu: value of $\mu_{qmax}$.

• –

gamma: element of the input vector gamma associated with the estimated "Meta-Model".

• –

Meta-Model: a RKHS ridge group sparse or RKHS group lasso object associated with the tuning parameters mu and gamma (see Table 13).

6.2.2 Companion functions

calc$\_$Kv function:

For a given value of Dmax and a chosen kernel (see Table 1), this function calculates the Gram matrices $K_{v}$, $v\in\mathcal{P}_{\text{Dmax}}$, and returns their associated eigenvalues and eigenvectors. This function has,

• •

four mandatory input arguments:

• –

$Y$, $X$, kernel, Dmax (see Table 12).

• •

three facultative input arguments:

• –

correction: logical, set as TRUE to make correction to the matrices $K_{v}$ (see Section 3.1). It is set as TRUE by default.

• –

verbose: logical, set as TRUE to print the group for which the correction is done. It is set as TRUE by default.

• –

tol: scalar to be chosen small, set as $1e^{-8}$ by default.

The calc$\_$Kv function returns a list of two components "kv" and "names.Grp":

• •

kv: list of vMax components, each component is a list of,

• –

Evalues: vector of eigenvalues.

• –

Q: matrix of eigenvectors.

• •

names.Grp: vector of group names of size vMax.

RKHSgrplasso function:

For a given value of the tuning parameter $\mu_{g}$, this function fits the solution to the RKHS group lasso optimization problem by implementing Algorithm 1. This function has,

• •

three mandatory input arguments:

• –

$Y$ (see Table 12).

• –

Kv: list of the eigenvalues and the eigenvectors of the positive definite Gram matrices $K_{v}$ for $v=1,...,$vMax and their associated group names (output of the function calc$\_$Kv).

• –

mu: positive scalar indicates the value of the tuning parameter $\mu_{g}$ defined in Equation (15).

• •

two facultative input arguments:

• –

maxIter: integer, to set the maximum number of loops through all groups. It is set as $1000$ by default.

• –

verbose: logical, set as TRUE to print the number of current iteration, active groups and convergence criterion. It is set as FALSE by default.

This function returns a RKHS group lasso object associated with the tuning parameter $\mu_{g}$. Its output is a list of $13$ components:

• •

intercept, teta, fit.v, fitted, Norm.H, supp, Nsupp, SCR, crit, MaxIter, convergence, RelDiffCrit, and RelDiffPar (see Table 13).

mu$\_$max function:

This function calculates the value $\mu_{\text{max}}$ defined in Equation (16). It has two mandatory input arguments: the response vector $Y$, and the list matZ of the eigenvalues and eigenvectors of the positive definite Gram matrices $K_{v}$ for $v=1,...,$vMax. This function returns the $\mu_{max}$ value.

pen$\_$MetMod function:

This function produces a sequence of the RKHS meta-models associated with a given grid of values of the tuning parameters $\mu,\gamma$. Each RKHS meta-model in the sequence is the solution to the RKHS ridge group sparse optimization problem (obtained by implementing Algorithm 2) associated with a pair of values of $(\mu,\gamma)$ in the grid of values of $\mu,\gamma$. This function has,

• •

seven mandatory input arguments:

• –

$Y$ (see Table 12).

• –

gamma: vector of positive scalars. Values of the penalty parameter $\gamma$ in decreasing order.

• –

Kv: list of the eigenvalues and the eigenvectors of the positive definite Gram matrices $K_{v}$ for $v=1,...,$vMax and their associated group names (output of the function calc$\_$Kv).

• –

mu: vector of positive scalars. Values of the tuning parameter $\mu$ in decreasing order.

• –

resg: list of the RKHS group lasso objects associated with the components of "mu", used as initial parameters at Step $1$.

• –

gama_v and mu_v: vector of vMax positive scalars. These two inputs are optional. They are provided to associate the weights to the two penalty terms in the RKHS ridge group sparse criterion (13). In order to consider no weights, i.e. all the weights are equal to one, we set these two inputs to scalar zero.

• •

three facultative input arguments:

• –

maxIter: integer, to set the maximum number of loops through initial active groups at Step $1$ and maximum number of loops through all groups at Step $2$. It is set as $1000$ by default.

• –

verbose: logical, set as TRUE to print for each pair of the tuning parameters $(\mu,\gamma)$, the number of current iteration, active groups and convergence criterion. It is set as FALSE by default.

• –

calcStwo: logical, set as TRUE to execute Step $2$. It is set as FALSE by default.

The function pen$\_$MetMod returns a list of $l$ components, with $l$ being equal to the number of pairs of the tuning parameters $(\mu,\gamma)$. Each component of the list is a list of three components "mu", "gamma" and "Meta-Model":

• •

mu: positive scalar, an element of the input vector "mu" associated with the estimated "Meta-Model".

• •

gamma: positive scalar, an element of the input vector "gamma" associated with the estimated "Meta-Model".

• •

Meta-Model: a RKHS ridge group sparse object associated with the tuning parameters mu and gamma (see Table 13).

PredErr function:

By considering a testing dataset, this function calculates the prediction errors for the obtained RKHS meta-models. This function has eight mandatory input arguments:

• $-$

$X$, gamma, kernel, Dmax (see Table 12).

• $-$

$XT$: matrix of observations of the testing dataset with $n^{\text{test}}$ rows and $d$ columns.

• $-$

$YT$: vector of response observations of the testing dataset of size $n^{\text{test}}$.

• $-$

mu: vector of positive scalars. Values of the tuning parameter $\mu$ in decreasing order.

• $-$

res: list of the estimated RKHS meta-models for the learning dataset associated with the tuning parameters $(\mu,\gamma)$ (it could be the output of one of the functions RKHSMetMod or pen$\_$MetMod).

Note that, the same kernel and Dmax have to be chosen as the ones used for the learning dataset.

The function PredErr returns a matrix of the prediction errors. Each element of the matrix corresponds to the prediction error of one RKHS meta-model in "res".

prediction function:

This function calculates the predicted values for a new dataset based on the best RKHS meta-model estimator. It has six input arguments:

• $-$

X, kernel, Dmax (see Table 12).

• $-$

Xnew: matrix of new observations with $n^{\text{new}}$ rows and $d$ columns.

• $-$

res: list of the estimated RKHS meta-models for a learning dataset associated with the tuning parameters $(\mu,\gamma)$ (it could be the output of one of the functions RKHSMetMod, RKHSMetMod$\_$qmax or pen$\_$MetMod).

• $-$

Err: matrix of the prediction errors associated with the RKHS meta-models in "res" (output of the function PredErr).

The function prediction returns a vector of the predicted values based on the best RKHS meta-model estimator in "res". This function is available at GitHub.

SI$\_$emp function:

For each RKHS meta-model $\widehat{f}$, this function calculates the empirical Sobol indices for all the groups that are active in the support of $\widehat{f}$. This function has two input arguments:

• $-$

res: list of the estimated meta-models using RKHS ridge group sparse or RKHS group lasso algorithms (it could be the output of one of the functions RKHSMetMod, RKHSMetMod$\_$qmax or pen$\_$MetMod).

• $-$

ErrPred: matrix or NULL. If matrix, each element of the matrix corresponds to the prediction error of a RKHS meta-model in "res" (output of the function PredErr). Set as NULL by default.

The empirical Sobol indices are then calculated for each RKHS meta-model in "res", and a list of vectors of the Sobol indices is returned. If the argument "ErrPred" is the matrix of the prediction errors, the vector of empirical Sobol indices is returned for the best RKHS meta-model in the "res".