METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

arXiv:1701.04671v1 [math.ST] 17 Jan 2017

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

1. Introduction
We consider a Gaussian regression model
(1)

Y = m(X) + σε,

where X is a d random vector with a known distribution PX = P1 × . . . × Pd on X a compact
subset of Rd , and ε is independent of X, and distributed as a N (0, 1) variable. The variance σ 2 is
unknown and the number of variables d may be large. The function m is a complex and unknown
function from Rd to R. It may present strong non-linearities and high order interaction effects
between its coordinates. On the basis of a n-sample (Yi , Xi ), i = 1, . . . , n, we aim to construct
metamodels and perform sensitivity analysis in order to determine the infuence of each variable
or group of variables on the output.
Our approach combines methods for sensitivity analysis of a complex model and statistical
tools for sparse non-parametric estimation in multivariate Gaussian regression model. It rests
on the construction of a metamodel for aproximating the Hoeffding-Sobol decomposition of the
function m. This metamodel belongs to a reproducing kernel Hilbert spaces constructed as a
direct sum of Hilbert spaces leading to a functional ANOVA decomposition involving variables
and interactions between them. The estimation of the metamodel is carried out via a penalized
least-squares minimization allowing to select the subsets of variables X that contribute to predict
the output Y . Finally, the estimated metamodel allows to estimate the sensitivity indices of m.
A lot of work has been done around meta-modelling and sensitivity indices estimation.
For a complete account on global sensitivity analysis, see for example the book by Saltelli et
al. [29]. Let us briefly present the context of usual global sensitivity analysis. Suppose that we
are able to calculate the ouputs z of a model m for n realizations of the input vector X, such
that zi = m(Xi ) for i = 1, . . . , n. Starting from the values (zi , Xi ), i = 1, . . . n, the objectives
of meta-modelling and global sensitivity analysis are to quantify the influence of some subsets
of the variables X on the output z or are to approximate the function m by what is called a
metamodel. This metamodel helps to understand the behavior of the model, or allows to speed
up future calculation using it in place of the original model m.
In particular when the inputs variables X are independent, if m is square integrable, one may
consider the classical Hoeffding-Sobol decomposition [30, 33] that leads to write m according to
1

2

SYLVIE HUET AND MARIE-LUCE TAUPIN

its ANOVA functionnal expansion:
(2)

m(x) = m0 +

X

mv (xv ),

v∈P

where P denotes the set of parts of {1, . . . , d} with dimension 1 to d and where for all x ∈ Rd ,
xv denotes the vector with components xj for j ∈ v. The functions mv are centered and
orthogonal
in L2 (PX ) leading to the following decomposition of the variance of m: Var (m(x)) =
P
v Var (mv (xv )). The Sobol sensitivity indices, introduced by Sobol [30], are defined for any
group of variables xv , v ∈ P by
Var (mv (xv ))
Sv =
.
Var (m(x))
They quantify the contribution of a subset of variables x to the output m(x). Several approaches
are available for estimating these sensitivity indices, see for example Iooss and Lemaı̂tre [15] for
a recent review. Among all of them, let us consider the one based on metamodel construction that allows to directly obtain the sensitivity indices. Generally one consider metamodels
that correspond to an ANOVA-type decomposition, and that are candidate to approximate the
Hoeffding decomposition of m. The ANOVA-type decomposition leads to consider functions
defined as follows:
X
fv (xv ), EPX fv (xv ) = EPX fv (xv )fv′ (xv′ ) = 0 ∀v, v ′ ∈ P
(3)
f : X → R, f (x) = f0 +
v∈P

for functions fv that are chosen to belong to some functionnal spaces. The truncated (such that
to keep polynomials with degree less than some integer) polynomial Chaos construction was
proposed by Blatman and Sudret [4]. This truncated expansion is used as an approximation of
the Hoeffding decomposition of m. The authors propose a method for truncating the polynomial
Chaos expansion and then an algorithm based on least angle regression for selecting the terms
in the expansion. For the same purpose, Gu and Wu [13] propose an algorithm based on the
hierarchy principle (lower order effects are more likely to be important than higher order effects)
and on the heredity principle (interaction can be active only if one or all of its parent effects are
also active). This approach joins the one proposed by Bach [1] for variable selection based on
hierarchical kernel learning.
Durrande et al. [11] propose to approximate m by functions belonging to a reproducing kernel
Hilbert space (RKHS). The RKHS is constructed as a direct sum of Hilbert spaces leading to
a functional ANOVA decomposition (see Equation (3)), such that the projection of m onto the
RKHS is an approximation of the Hoeffding decomposition of m.
Following Lin and Zhang [20], Touzani and Busby [32] propose an algorithm to calculate
the penalized least-square estimator of m on the RKHS space, where the least-square criteria
is penalized by the sum of the norms of m on each Hilbert subspace. This group-lasso type
procedure allows both to select and calculate the non-zero terms in the functional ANOVA
decomposition.
Our objective is to propose an estimator of a metamodel which will approximate the Hoeffding
decomposition of m considering a Gaussian regression model defined at Equation (1) and to
deduce from this estimated metamodel, estimators for the sensitivity indices of m. Contrary to
the usual setting of sensitivity analysis where m(Xi ) is available, only the observations Y are
available, which leads us to consider the nonparametric multivariate regression setting.
Let us briefly describe the methods related to this regression setting and review their theoretical properties, starting with papers assuming an univariate additive decomposition for the
function m in the context of high-dimensional sparse additive models.
Pd Precisely, denote by
1-add
F
, the set of functions f defined on X such that f (x) = f0 + a=1 fa (xa ) where f0 is a
constant, and where for a = 1, . . . , d, the functions fa are centered and square-integrable with

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

3

respect to Pa . For each function f , the set Sf of indices a ∈ {1, . . . d} such that fa is not
identically zero is called the actice set of f .
Ravikumar et al. [27] propose a group-lasso procedure where each function fa is approximated
by its truncated decomposition on a basis of functions. They provide an algorithm and, assuming
that the function m belongs to the set F 1-add and that Sm is sparse, prove the consistency of
the active set and of the risk of the estimator of m.
Meier et al. [22] propose to combine a sparsity penalty (group-lasso) and a smoothness penalty
(ridge) for estimating m(x). Considering the fixed design framework, they establish some oracle
properties of the empirical risk for estimating the projection of m onto the set of univariate
additive functions F 1-add .
Raskutti et al. [26] consider the case where each univariate function fa belongs to a RKHS
and as Meier et al. combine a sparsity and a smoothness penalty. Assuming that the d variables
X are independent, they derive upper bounds for the integrated and the empirical risks, as
well as a lower bound for the integrated risk over spaces of sparse additive models whose each
component is bounded with respect to the RKHS norm.
Additive sparse modelling is too restrictive in practical settings because it does not take
into account interactions between variables that may affect the relationship between Y and
X. The generalization of additive smoothing splines to interaction smoothing splines leading
to an ANOVA-type decomposition (see Equation (3)) was proposed by several authors (see for
example Wahba [34], Friedman [12], Wahba et al. [35]).
To control smoothness and to enforce sparsity in the ANOVA-type decomposition, Gunn
and Kandola [14] propose to consider the ANOVA decomposition as a weighted linear sum of
kernels and to use a lasso penalty on the weights to select the terms in the decomposition as
well as a ridge penalty to ensure smoothness of the kernel expansion. The COSSO proposed by
Lin and Zhang [20] is based on smoothness penalty defined as the sum of the RKHS-norms of
the functions fv . The authors study existence and rate of convergence of the estimator. In a
more general framework, where the function m is written as a linear span of a large number of
kernels, Koltchinskii and Yuan [17] established oracle inequalities on the excess risk assuming
that the function m has a sparse representation (the set of v ∈ P such that fv∗ is non zero is
sparse). The authors generalized their results to a penalty function that combines sparsity and
smoothness [18], as proposed by Meier et al. [22] and Raskutti et al. [26].
Recently Kandasamy and Wu [16] proposed an estimator called SALSA, based on a ridge
penalty, where the ANOVA-type decomposition is restricted to set v ∈ P such that |v| ≤ Dmax .
The authors propose to choose Dmax via a cross-validation procedure.
Our contributions. Using the functionnal ANOVA-type decomposition as proposed by Durrande et al. [11], we propose an estimator of a metamodel which approximates the Hoeffding
decomposition of m. Following the most recent works in the framework of nonparametric estimation of sparse additive models, we propose a penalized least-square estimator where the
penalty function enforces both the sparsity and the smoothness of the terms in the decomposition. We show that our estimator satisfies an oracle inequality with respect to the empirical
and integrated risks.
Our procedure allows both to select and estimate the terms in the ANOVA decomposition,
and therefore, to select the sensitivity indices that are non-zero and estimate them. In particular
it makes possible to estimate Sobol indices of high order, a point known to be difficult in practice.
Finally, using convex optimization tools, we develop an algorithm (available on request) in
R [25], for calculating the estimator. A simulation study shows the good performances of our
estimator in practice.
The paper is organised as follows: the procedure for estimating a metamodel is presented in
Section 2.1. The theoretical properties of the metamodel estimator are stated in Theorem 2.1 and
Corollary 2.1 whose proofs are posponed in Section 5. Then we briefly explain in Section 2.3, the

4

SYLVIE HUET AND MARIE-LUCE TAUPIN

RKHS construction based on ANOVA kernels, and present the estimators of the Sobol indices
in Section 2.5. Section 3 is devoted to the calculation of the estimator and Section 4 to the
simulation study.
2. Meta-modelling and Sobol indices estimation
2.1. Procedure for meta-modelling. We propose to consider a functionnal space based on
the tensorial product of Reproducing Kernel Hilbert spaces (RKHS), and to approximate the
unknown function m by its projection denoted f ∗ on such a space. Precisely if P denotes the
set of parts of {1, . . . , d} with dimensions 1 to d, then we consider the space H of functions f
from X ⊂ Rd to R such that
X
f (x) = f0 +
fv (xv ),
v∈P

where for all x ∈ Rd , xv denotes the vector with components xj for j ∈ v, f0 a constant,
and where each function fv : R|v| → R belongs to a RKHS Hv associated to a kernel denoted
kv and equipped with a scalar product h·, ·iHv . The RKHS’s Hv , v ∈ P are constructed in
order to ensure that the functions in Hv are centered, and to ensure orthogonality between the
spaces Hv , v ∈ P with respect to L2 (PX
P) (the construction of such RKHS spaces is detailed in
Section 2.3). It follows that f ∗ = f0∗ + v∈P fv∗ which minimizes
km − f k2L2 (PX ) = EX (m(X) − f (X))2

over functions f ∈ H, is an approximation of the Hoeffding decomposition of m. Therefore if
the Hoeffding decomposition of m is written as in Equation (2), each function fv∗ approximates
the function mv . Since P is the set of parts of {1, . . . , d}, the number of functionsfv∗ is related
to the cardinality of P = 2d − 1 that may be huge. The construction is thus associtated to a
selection strategy.
The selection of fv∗ in f ∗ is based on a ridge-group-sparse type procedure which minimizes
the penalized least-squares criteria over functions f ∈ H. The least-squares criteria is penalized
in order to both select few terms in the additive decomposition of f over sets v ∈ P, and to
favour smoothness of the estimated fv . The ridge regularization is ensured by controling the
norm of fv in the Hilbert space Hv for all v, and the group-sparse regularization is strengthened
by controling the empirical norm of fv , defined as
v
u n
u1 X
fv2 (Xv,i ).
kf kn = t
n
i=1

P

For any f ∈ H such that f = f0 +
L(f ) be defined as
(4)

v∈P

fv , and for some tuning parameters (µv , γv , v ∈ P), let

!2
n
X
X
X
1X
µv kfv kHv +
γv kfv kn .
Yi − f 0 −
fv (Xv,i ) +
L(f ) =
n
i=1

v∈P

v∈P

v∈P

For some positive R, let us define the set of functions F
(
)
X
(5)
F = f such that f = f0 +
fv , fv ∈ Hv , kfv kHv < R .
v∈P

Then the estimator fb is defined as
(6)

fb = argmin {L(f ), f ∈ F} .

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

5

2.2. Theoretical result: oracle inequality for metamodel. In this section we state the
oracle inequality for estimating a metamodel which approximates the Hoeffding decomposition
of the unknown function m. Let us
Pbegin with notation and assumptions.
For a function f ∈ H, f = f0 + v∈P fv , we denote by Sf its support and |Sf | its cardinality.
More precisely
(7)

Sf = {v ∈ P, fv 6= 0} .

We consider RKHS spaces Hv , v ∈ P satisfying the following assumptions:
• for all fv ∈ Hv , EX fv2 (X) < ∞,
• for all fv ∈ Hv , fv′ ∈ Hv′ , EX fv (X) = 0 and EX fv (X)fv′ (X) = 0,
• there exists R′ > 0 such that

(8)

∀fv ∈ Hv kfv k∞ = sup {|fv (X)|, X ∈ X } ≤ R′ .

For each v ∈ P, let ωv,k , for k ≥ 1 be the eigenvalues of the operator associated to the self
reproducing kernel kv , arranged in the decreasing order. Let us define the function Qn,v (t), for
positive t, as follows:
s
5X
min(t2 , ωv,k ),
(9)
Qn,v (t) =
n
k≥1

and for some ∆ > 0 let νn,v be defined as follows

νn,v = inf t such that Qn,v (t) ≤ ∆t2 .
(10)

For each v ∈ P, νn,v refers to the so-called critical univariate rate, the minimax-optimal rate for
L2 (PX )-estimation of a single univariate function in the hilbert space Hv (e.g. Mendelson [23]).
Our choices of regularization parameters and rates are specified in terms of the quantities:
n
o
p
λn,v = max νn,v , d/n
(11)


X
X
d2n (f ) = max 
γv2 ,
µv  .
(12)
v∈Sf

v∈Sf

Theorem 2.1. Let us consider the regression model defined at Equation (1). Let (Yi , Xi ),
i = 1, . . . , n be a n-sample with the same law as (Y, X). Let fb be defined by (6).
If there exist constants Cl , l = 1, 2, 3, C1 ≥ 1, and 0 < η < 1 such that if the following
conditions are satisfied:

r 
µv
γv
,
, nλ2n,v ≥ −C2 log λn,v ,
(13)
for all v ∈ P, λn,v ≤ min
C1
C1
and
for all f ∈ F, d2n (f ) ≤ 1,
then, with probability greater than 1 − η,


 

X
d|Sf |  2 
2
νn,v
+
km − fbk2n ≤ C3 inf km − f k2n + 
.
σ

f ∈F 
n
v∈Sf

6

SYLVIE HUET AND MARIE-LUCE TAUPIN

Comments.
• The result is non asymptotic in the sense that it is shown for any (n, d). Nevertheless, the
upper bound is relevant when the infimum is reached for functions f whose decomposition
in H is sparse, and when d is small face to n. In fact, the coefficient d occuring in the rate
d|Sf |/n comes from the logarithm of the cardinality of P equal to log(2d − 1). When d
is large, it may be judicious to limit the decomposition of functions in H, to interactions
of limited order, so that the number of terms in the decomposition stays of the order
log(d). We will comment the term linked to the convergence rate on the space H in
Section 2.4.
• The term km−f k2n refers to the usual bias term quantifying the approximation properties
of the Hilbert space H as a distance between the true m and f , its approximation into
H.
• Koltchinskii et Yuan [18] considered what is called the multiple kernel learning problem,
where the functions in H have an additive representation over kernel spaces. They
do not assume that the variable X are independent, nor that the kernel spaces satisfy
an orthogonality condition. In return, they assume that some decomposability type
properties are satisfied, and they introduce some characteristics related to the degree of
“dependence” of the kernel spaces.
• In the particular case when the decomposition is limited to the main effects of the
variables, then the problem comes back to the classical nonparametric additive model.
The theoretical properties of the estimator based on a ridge-group-sparse type procedure
have already been established (see for example Meier et al. [22], and Raskutti et al. [26]).
• Weights in the penalty terms may be of interest for applications. The theoretical result
highlights that the tuning parameters (µv , γv ) should depend on the decreasing of the
eigenvalues of the kernel defining Hv . Besides, we may be interested by introducing
weigths that favor small order interaction terms.
• Because we aim to approximate the Hoeffding decomposition of m, we need to have
orthogonality between the spaces Hv , v ∈ P. This condition, required by our objective,
is also a key point in the proof of the Theorem, when the problem is to compare the
euclidean norm of functions in H with the norm in L2 (PX ). At this step we need to
assume that d2n (f ) ≤ 1 to conclude.
• We do not assume that the functions in F are uniformly globally bounded, that is that
the sup{|f (x)|, x ∈ X } is bounded by a constant that does not depend on f . Instead
we assume that each function within the unit ball of the Hilbert space Hv is uniformly
bounded by a constant
P multiple of its Hilbert norm. In fact, the functions f in the space
F, written as f0 + v fv satisfy that for each group v, fv is uniformly bounded. This
assumption is easily satisfied
p as soon as the kernel kv is bounded on the compact set X .
Indeed, kfv k∞ ≤ supX∈X kv (Xv , Xv )kfv kHv . We refer to [26] for a discussion on that
subject and a comparison with the work of Koltchinskii et Yuan [18].
• The constants R, R′ and ∆ in Equation (5), respectively (8) and (10), could depend on
v. Everything in the proof would be similar.
• Assuming that nλ2n,v ≥ −C2 log λn,v allows to control the probability of the union of |P|
√
√
events. This is a mild condition, satisfied for λn,v = Kn / n for Kn of the order log n.

The following corollary gives an upper bound of the risk with respect to the L2 (PX ) norm. Its
proof is given Section 5.2 page 20.

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

7

Corollary 2.1. Under the assumptions of Theorem 2.1, we have that, with probability greater
than 1 − η, for some constant C4 ,


 


X
d|S
|
f  2
2
νn,v
+
km − fbk2L2 (PX ) ≤ C4 inf km − f k2n + km − f k2L2 (PX ) + 
.
σ

f ∈F 
n
v∈Sf

Sketch of proof of Theorem 2.1. We give here a sketch of the proof and we postpone to Section 5
for complete statements. In particular, we denote by C constants that vary from an equation
to the other, and we assume that R = 1, for all v such that
∀f ∈ F, ∀v ∈ P, kfv kHv ≤ 1.

(14)

The proof of Theorem 2.1 starts in the same way as the proof of Theorem 1 in Raskutti et
al. [26]. Nevertheless it differs in several points, in particular because the terms occuring in
the decomposition of functions in H depend on several variables and thus are not independent.
Indeed, fv (Xv ) and fv′ (Xv′ ) are not independent as soon as the groups v and v ′ share some of
the variables Xa , a = 1, . . . , d. Moreover, we do not assume that the function m is in F.
Starting from the definition of fb, some simple calculation (see Equation (27)) give that for
all f ∈ F
n

X
1X
2
2
b
γv kfbv − fv kn + µv kfbv − fv kHv .
εi (fb(Xi ) − f (Xi )) +
Ckm − f kn ≤ km − f kn +
n
i=1

v∈Sf

P

If we set g = fb − f , then g ∈ H, g = g0 + v gv , with gv = fbv − fv , and for each v, kgv kHv ≤ 2
(see Equation (14)).
The main problem is now to control the empirical process. For each v, letting λn,v as in (11),
we state (see Lemma 5.1, page 17) that, with high probability,
if kgv kn ≤ λn,v kgv kHv

then

if kgv kn > λn,v kgv kHv

then

n
X

i=1
n
X
i=1

εi gv (Xv,i ) ≤ Cnλ2n,v kgv kHv
εi gv (Xv,i ) ≤ Cnλn,v kgv kn .

Therefore, if for all v, µv and γv satisfy Equation (13), we deduce that with high probability
(setting g = fb − f )

X
X
γv kfbv kn + µv kfbv kHv .
(γv kgv kn + µv kgv kHv ) +
Ckm − fbk2n ≤ km − f k2n
v∈Sf

v∈S
/ f

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

X
X
(γv kgv kn + µv kgv kHv ) .
γv kfbv kn + µv kfbv kHv ≤ C
v∈S
/ f

v∈Sf

Putting the things together, and noting again that kgv kHv ≤ 2, we obtain he following upper
bound
X
( µv + γv kgv kn ) .
Ckm − fbk2n ≤ km − f k2n +
v∈Sf

P
P
The last important step consists in comparing v∈Sf kgv kn to k v∈Sf gv kn . More precisely, it
can be shown (see lemma 5.3, page 18) that for all v ∈ P, with high probability, we have
kgv kn ≤ 2kgv kL2 (PX ) + γv .

8

SYLVIE HUET AND MARIE-LUCE TAUPIN

P
P
Using the orthogonality assumption between the spaces Hv , we have v∈Sf kgv k2L2 (PX ) = k v gv k2L2 (PX ) ,
and thus we get
X
X
X
(fbv − fv )2 k2L2 (PX ) .
γv2 + k
µv +
Ckm − fbk2n ≤ km − f k2n +
v∈Sf

v∈Sf

v∈Sf

Finally it remains to consider different cases according to the rankings of kfb− f k2L2 (PX ) , kfb− f k2n
P
and v∈Sf µv + γv2 to get the result of Theorem 2.1.

2.3. Practical construction of the RKHS. One of the key point is to construct the space H
such that the terms of the decomposition of a function f in H correspond to its Hoeffding-Sobol
decomposition (see Sobol [31] and Van der Vaart [33], p. 157).
Let us recall that the Hoeffding decomposition of the function m consists in writting m as in
Equation (2) where for all v, v ′ in P,
EX (mv (Xv )) = EX (mv (Xv )mv′ (Xv′ )) = 0.

Therefore any function f ∈ H is a candidate to approximate the Hoeffding decomposition of m.
The construction of spaces H, based on ANOVA kernels, was given by Durrande et al. [11].
Let X = X1 × . . . × Xd be a compact subset of Rd . For each coordinate a ∈ {1, · · · , d}, we
choose a RKHS Ha and its associated kernel ka defined on the set Xa ⊂ R such that the two
following properties are satisfied
(1) ka : p
Xa × Xa → R is Pa × Pa measurable
(2) EPa ka (Xa , Xa ) < ∞
⊥

The RKHS Ha may be decomposed as Ha = H0a ⊕ H1a , where

H0a = {fa ∈ Ha , EPa (fa (Xa )) = 0}
H1a = {fa ∈ Ha , fa (Xa ) = C} ,

the kernel k0a associated to the RKHS H0a being defined as follows (see Berlinet et ThomasAgnan [3]):
EU ∼Pa (ka (xa , U ))EU ∼Pa (ka (x′a , U ))
k0a (xa , x′a ) = ka (xa , x′a ) −
.
E(U,V )∼Pa ×Pa ka (U, V )
The ANOVA kernel is defined as
d
Y
X
Y

′
k(x, x ) =
1 + k0a (xa , x′a ) = 1 +
kv (xv , x′v ), where kv (xv , x′v ) =
k0a (xa , x′a ),
a=1

a∈v

v∈P

and the corresponding RKHS

H=

⊗da=1



⊥

1 ⊕ H0a



=1+

X

v∈P

Hv ,

where the RKHS Hv is associated with kernel kv .
According to this construction, any function f ∈ H satisfies
X
f (x) = hf, k(x, ·)iH = f0 +
fv (x),
v∈P

where fv (x) = hf, kv (x, ·)iH depends on xv only. For all v ∈ P, fv (xv ) is centered and for all
v ′ 6= v, fv (xv ) and fv′ (xv′ ) are uncorrelated. We get thus the Hoeffding decomposition of f ,
which is an approximation of the Hoeffding decomposition of m.
Let us remark that the decomposition based on the truncated polynomial Chaos expansion,
used for sensitivity analysis (see for example Blatman and Sudret [4]), is based on the distribution
of X, and only the choice of the truncation handles the smoothness of the approximation. The

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

9

construction of the RKHS spaces described above, allows to consider functionnal spaces that
suit well to the smoothness of the function m, irrespectively of the distribution of X. Indeed,
the kernels k0,a depend on the distribution of X only for calculating the projection onto the
space of constant functions.
2.4. Rate of convergence. Let us come back
oracle inequality given by Theorem 2.1
P to the
2 . For the sake of simplicity let us consider
and discuss the rate of convergence given by v∈Sf νn,v
the case where the variables X1 , . . . , Xd have the same distribution P1 onQX1 ⊂ R, and where
the unidimensionnal kernels k0a are all identical, such that kv (xv , x′v ) = a∈v k0 (xa , x′a ). The
kernel k0 admits an eigen expansion given by
X
k0 (x, x′ ) =
ω0,ℓ ζℓ (x)ζℓ (x′ )
ℓ≥1

where the eigenvalues ω0,ℓ are non negative and ranged in the decreasing order, and where the
ζℓ are the associated eigen functions, orthonormal with respect to L2 (P1 ). Therefore the kernel
kv admits the following expansion
X

′

kv (xv , xv ) =

|v|
Y

ω0,ℓa

ℓ=(ℓ1 ...ℓ|v| ) a=1

|v|
Y

ζℓa (xa )

| {z } a=1
| {z
ωv,ℓ

ζv,ℓ (xv )

|v|
Y

′

ζℓa (xa ) .

} |a=1 {z

ζv,ℓ (x′v )

}

Consider now the case where the eigenvalues ω0,ℓ are decreasing at a rate ℓ−2α for some α > 1/2.
Q|v|
Then ωv,ℓ is of order ( a=1 ℓa )−2α and
5X
min(t2 , ωv,ℓ )
Q2n,v (t) =
n
ℓ

≤

5 2X
C(|v|, α)
I(ℓ1 × . . . × ℓ|v| ≤ t−1/α ) +
t
n
n
ℓ

5
C(|v|, α)
t
+
n
n
where for any set E, I(E) denotes the characteristic function on E, and C(|v|, α) is a constant
depending on |v| and α. Therefore, the rate νn,v defined at Equation (10) is of order n−α/(2α+|v|)
which is the usual nonparametric rate for |v|-dimensionnal functions. It follows that the RKHS
space H should be chosen such that the unknown function m is well approximated by sparse
functions in H with low order of interactions.
≤

|v|
2− α

2.5. Sensitivity analysis. The orthogonality between two terms in the Hoeffding decomposition given at Equation (2) leads to the additive decomposition of the variance of m(x):
X
Var (m(x)) =
Var (mv (xv )) .
v∈P

Each of these variance terms are related to Sobol indices. For example, the Sobol indice linked
with the interaction between variables xv is defined as
Sv =

Var (mv (xv ))
,
Var (m(x))

or the global Sobol indices for the variable xa is
P
v⊃{a} Var (mv (xv ))
.
Ga =
Var (m(x))

10

SYLVIE HUET AND MARIE-LUCE TAUPIN

Thanks to the orthogonality property of functions in H, the variance of m(x) will be estimated
by


X
d (mv (xv )) , where Var
d (mv (xv )) = EX fbv2 (Xv ) = kfbv k2 2
d (m(x)) =
Var
Var
L (PX ) .
v∈P

d (mv (xv )) with the variance of mv (xv ). Thanks to the following inequality
Let us compare Var
kfbv kL2 (PX ) − kmv kL2 (PX ) ≤ kfbv − mv kL2 (PX ) ≤ kfb − mkL2 (PX ) ,

and to Corollary 2.1, we get the following result

d (mv (xv )) ≤ kfb − mk2 2
if mv ≡ 0 then Var
L (PX )

if kmv (xv )kL2 (PX ) ≥ c > 0 then

d (mv (xv ))
Var
− 1 ≤ kfb − mk2L2 (PX ) /c.
Var (mv (xv ))

In practice, in order to avoid calculating the variance of fbv2 (Xv ), one may use an estimator
based on the empirical variances of functions fbv . Precisely, if fbv,· is the mean of the fbv (Xv,i ),
for i = 1, . . . , n, then
(15)

d emp (mv (xv )) =
Var

2
1 Xb
fv (Xv,i ) − fbv,· .
n−1
n

i=1

3. Calculation of the estimator
The functional minimization problem described at Equation (6) is equivalent to a parametric minimization problem. Indeed, we know that if H is a RKHS associated with a kernel
k : X ×P
X → R, then for all (x1 , . . . , xn ) ∈ XPn , and for all (α1 , . . . , αn ) ∈ Rn , the function
n
n
′
2
it can be
f (·) =
i=1 αi k(xi , ·) is in H and kf kH =
i,i′ =1 αi αi′ k(xi , xi′ ). In particular,
P
shown that the solution to our minimization problem is written as f = f0 + v∈P fv where
P
fv (·) = ni=1 θvi kv (Xvi , ·) for some parameter θ ∈ Rn|P| with components (θv,i , i = 1, . . . , n, v =
1, . . . , |P|).
Let k · k denotes the usual euclidean norm in Rn . For each v ∈ P, let Kv be the n × n matrix
b be the
with components (Kv )i,i′ = kv (Xvi , Xvi′ ) that satisfies t(K 1/2 )K 1/2 = K. Let fb0 and θ
minimizer of the following penalized least-squares criteria:
X
X
1 X
1
Kv θ v k2 + √
γv kKv θ v k +
µv kKv1/2 θ v k.
(16)
C(f0 , θ) = kY − f0 1n −
n
n
v∈P

v∈P

Then the estimator fb defined at Equation (6) satisfies
fb(x) = fb0 +

X

v∈P

fbv (xv ) with fbv (xv ) =

v∈P

n
X
i=1

θbv,i kv (Xv,i , xv ).

b using a block
Because C(f0 , θ) is a convex and separable criteria, we propose to calculate θ
coordinate descent algorithm described in the following section.
1/2
Note that the estimator fb defined at Equation (6) should satisfy kfv kHv = kKv θ v k ≤ R.
Usually we have no idea of the value of this upper-bound in practice, and we propose to remove
this contraint in the optimization procedure. Nevertheless, if one wants to consider such an
additional constraint, the problem can be solved at the price of some additional complication,
considering a Lagrangian method for example.

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

11

3.1. Algorithm. We will assume that for all v ∈ P the matrices Kv are strictly definite positive.
If it not the case, one modify Kv by Kv + ξIn where ξ is a small positive value, in order to ensure
positive definiteness.
Using a coordinate descent procedure, we minimize the criteria C(f0 , θ) along each group
v at a time. At each step of the algorithm, the criteria is minimized as a function of the
current block’s parameters, while the parameters values for the other blocks are fixed to their
current values. The procedure is repeated until convergence, considering for example that the
convergence is obtained if the norm of the difference between two consecutive solutions is small
enough. See for exemple Boyd et al. [7] for optimization in such context.
For the sake of simplicity, we consider the minimization of the following criteria:
X
X
X
µ′v kKv1/2 θ v k.
(17)
C ′ (f0 , θ) = kY − f0 −
Kv θ v k2 +
γv′ kKv θ v k +
v∈P

v∈P

v∈P

√
= nγv and µ′v = nµv , this is exactly the same criteria as the one defined at
Taking
Equation (16).
Let us begin with the constant term f0 . Because the penalty function does not depend on f0 ,
minimizing C ′ (f0 , θ) with respect to f0 for fixed values of θ leads to
γv′

f 0 = Y· −

(18)

n
XX
v

(Kv θ v )i /n,

i=1

where Y· denotes the mean of Y and (Kv θ v )i denotes the i-th component of Kv θ v . In what
follows, we consider a group v, and fix the values of the parameters for all the other groups. We
describe the algorithm and postpone the proofs in Section 5.6.
Let us first consider the case where both µ′v and γv′ are non zero. If ∂Cv′ denotes the subdifferential of C ′ (f0 , θ) with respect to θ v , we need to solve 0 ∈ ∂Cv′ , which is equivalent to
− 2Kv Rv + 2Kv2 θ v + γv′ sv + µ′v tv = 0,

(19)
where

Rv = Y − f 0 −
and where sv and tv satisfy:
if θ v = 0
if θ v 6= 0

X

Kw θ w

w6=v

kKv−1 sv k ≤ 1, and kKv−1/2 tv k ≤ 1
Kv θ v
Kv2 θ v
.
, and tv =
sv =
1/2
kKv θ v k
kKv θv k

The first task is to obtain necessary and sufficient conditions for which the solution θ v = 0 is
the optimal one. Let
n
o
(20) J(t) = k2Rv − µ′v Kv−1 tk2 , and J ∗ = argmin J(t), for t ∈ Rn such that kKv−1/2 tk ≤ 1 .

Then the solution to Equation (19) is zero if and only if J ∗ ≤ γv′2 . Calculating J ∗ is a ridge
regression problem that can be easily solved (see Propositions 5.2 and 5.3 in Section 5.6).
If the solution to Equation (19) is not θ v 6= 0, the problem is to solve the subgradient equation:
!−1
γv′
µ′v
Rv .
In + Kv +
Kv
(21)
θv =
1/2
2kKv θ v k
2kKv θ v k

Because θ v appears in both sides of the equation, a numerical procedure is needed (see Proposition 5.4).

12

SYLVIE HUET AND MARIE-LUCE TAUPIN

Other cases.
• If all the µ′v are equal to 0, the parameters θ v are not identifiable.
• If all the γv′ are equal to 0, then we have to solve a classical group-lasso problem with
1/2
respect to the parameters θ ′v defined as θ ′v = Kv θv for all v ∈ P.
′
′
• Let v such that µv = 0, γv 6= 0 and assume that at least one of the µ′w is non zero for
w ∈ P, w 6= v. Then it is shown in Proposition 5.1 that θv = 0 if and only if 2kRv k ≤ γv′ .
1/2
• In the same way, if γv′ = 0, and µ′v 6= 0, then θ v = 0 if and only if 2kKv Rv k ≤ µ′v .

. Finally the algorithm is the following:

(1) Start with an initial value θ = θ 0 .
(2) Calculate f0 using Equation (18). For the group v, calculate Rv and determine if the
group v should be excluded or not either by solving the problem defined at Equation (20)
if µ′v 6= 0 and γv′ 6= 0, or directly if one of them equals 0. If it is the case, set θ v = 0. If
not, solve Equation (21) to obtain θ v .
(3) Iterate step 2. over all the groups v.
(4) Iterate step 2. and 3. until convergence.
3.2. Choice of the tuning parameters. For each value of the tuning parameters (µ′v , γv′ ), v ∈
P, the algorithm provides an estimate of m and of the Sobol indices. The problem for choosing
these parameters values is crucial. We propose to restrict this choice by considering tuning
parameters proportionnal to known weights: for all v ∈ P, µ′v = µωv and γv′ = γζv , where the
weights ωv and ζv are fixed. For example, one can take weights that increase with the cardinal
of v in order to favour effects with small interaction order between variables. Or, according to
2 and ζ = ν
the theoretical result given at Theorem 2.1, we can choose ωv = νbn,v
bn,v , where νbn,v
v
is an estimate of νn,v based on the eigenvalues of the matrix Kv . Any other choice, depending
on the problem of interest, may be relevant.
Once the weights are chosen, we estimate m, on the basis of a learning data set (Yi , Xi ), i =
1, . . . n, for a grid of values of (µ, γ). We first set γ = 0, and calculate µmax the smallest value
of µ such that the solution to the minimization of
X
X
ωv kKv1/2 θ v k,
Kv θ v k2 + µ
kY − f0 −
v∈P

v∈P

is θ v = 0 for all v ∈ P. Then we can consider µℓ = µmax 2−ℓ for ℓ ∈ {1, . . . , ℓmax }, as a grid of
values for µ. The grid of values for γ may be chosen after few attempts.
For choosing the final estimator, say fb, we suppose that we have at our disposal a testing
data set (YiT , XTi ), i = 1, . . . nT , and we propose two procedures.
Proc. GS: The first one, called the Group-Sparse procedure uses the testing data set for
estimating the prediction error. Precisely, for each value of (µ, γ) in the grid, let fb(µ,γ) (·)
be the estimation of m obtained with the learning data set. Then
nT
2
1 X T
Yi − fb(µ,γ) (XTi )
PE(µ, γ) = T
n
i=1

estimates the prediction error, and we propose to choose the pair (µ, γ) that minimizes
PE(µ, γ), say (b
µ, γ
b). Finally the estimator, denoted fbGS is defined as fbGS = fb(bµ,bγ) .
Proc. rdg: The second procedure, called the ridge procedure, considers the collection of
when (µ, γ) belongs to the grid. For each
supports composed of the different Sfb
(µ,γ)
support S in this collection, we estimate f using a ridge procedure assuming that the

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

rdg
fλ,S

13

rdg
support of f is S: for a given λ, fλ,S
is defined as follows:


!2
n

1 X
X
X
X
kfv k2Hv , f = f0 +
f v , f v ∈ Hv .
Yi − f 0 −
fv (Xv,i ) + λ
= argmin

n
i=1

v∈S

v∈S

v∈S

We choose a grid of values for λ, and for each S in the collection and λ in the grid, we
estimate the prediction error PE(λ, S) defined as follows:
nT
2
1 X T
rdg
Yi − fλ,S
(XTi ) .
PE(λ, S) = T
n
i=1

b
For each S in the collection, let λ(S)
be the minimizer of PE(λ, S) when λ varies in
b
b
the grid, and let S be the minimizer of PE(λ(S),
S), then the estimator denoted fbrdg is
rdg
rdg
b
defined as f
= fb b b.
λ(S),S

If a testing data set is not available, we can use the classical V-fold cross validation either to
estimate PE(µ, γ) or to estimate PE(λ, S).
4. Simulation study
In order to evaluate the performances of our method for estimating a meta-model and sensitivity indices of a function m we carried out a simulation study. We consider the g-function of
Sobol defined on [0, 1]d as
d
Y
|4xa − 2| + ca
m(x) =
, ca > 0,
1
+
c
a
a=1

whose Sobol indices can be expressed analytically (see Saltelli et al. [28]). Following Durrande
et al. [11], we take d = 5 and (c1 , c2 , c3 , c4 , c5 ) = (0.2, 0.6, 0.8, 100, 100). The lower the value
of ca , the more significant the variable xa . The variables Xa , a = 1, . . . , d are independent
and uniformly distributed on [0, 1]. We consider the regression model Yi = m(Xi ) + σεi , for
i = 1, . . . , n, with N (0, 1) independent error terms εi .
For each simulation, we generate three independent data sets as follows : a Latin Hypercube
Sample of the inputs is simulated to give the matrix X with n rows and d = 5 columns, and a
n-sample of independent centered Gaussian variable with variance 1 is simulated. This operation
is repeated three times in order to obtain the learning and testing data sets and a third data
set for estimating the estimators performances. As explained in Section 3.2, we choose optimal
values of the tuning parameters, (µ, γ) by minimizing a prediction error PE, and get an estimator
fb of m, as well as estimates of the Sobol indices:
d (mv (xv ))
Var
.
Sbv = P
d (mw (xw ))
Var
w∈P

Let Ωv be the matrix whose components satisfy
Y

EU ∼Pa k0a (U, Xa,i )k0a (U, Xa,i′ ) for i, i′ = 1, . . . , n.
(Ωv )i,i′ =
a∈v

d (mv (xv )) is calculated as follows:
The estimator Var

d (mv (xv )) = α
Var
bTv Ωv α
bv .

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

14

SYLVIE HUET AND MARIE-LUCE TAUPIN

σ=0
n = 50 n = 100 n = 200 n = 50
Proc. GS 0.814
0.920
0.959
0.737
0.976
0.989
0.763
Proc. rdg 0.874
Table 1. Estimated coefficient of determination R2
and σ, with the Matérn kernel.

σ = 0.2
n = 100 n = 200
0.835
0.889
0.854
0.892
for different values of n

σ=0
σ = 0.2
n = 50 n = 100 n = 200 n = 50 n = 100 n = 200
Proc. GS 0.033
0.0137
0.0139
0.051
0.028
0.020
0.0009
0.0007
0.042
0.022
0.013
Proc. rdg 0.011
Table 2. Estimated empirical risk ER for different values of n and σ with the
Matérn kernel.

Performance indicators. To evaluate the performances of our method for estimating a metamodel, we use the classical coefficient of determination R2 estimated using the third data set
(YiP , XPi ), i = 1, . . . , n:
2
Pn  P
b(XP )
Y
−
f
i
i
i=1
R2 = 1 − Pn
2 .
P
P
i=1 Yi − Y·
Moreover we calculate the empirical risk ER = km − fbk2n . For each simulation s, we get Rs2 and
ERs and we report the means of these quantities over all simulations.
Similarly, for each v, and each simulation s, we get Sbv,s and we report its mean, Sbv,· , its
estimated standard-error, and to sum up the behaviour of our procedure for estimating the
sensitivity indices, we estimate the global error, denoted GE, defined as follows
X
GE =
(Sbv,s − Sv )2 .
v

In order to assess the performances of our procedure for selecting the greatest Sobol Indices,
precisely those that are greater than some small quantity as ρ = 10−4 , we calculate for each
group v ∈ P, the percentage of simulations for which v is in the support of the estimator fb.
Then we average these quantities, on one hand over groups v such that Sv > ρ, and on the other
hand, over groups v such that Sv ≤ ρ. Let us denote these quantities pSelSv >ρ and pSelSv ≤ρ
respectively.

Simulation design. We present the results for n ∈ {50, 100, 200}, σ ∈ {0, 0.2}. For all a =
1, . . . , d, the kernels ka are the same: we considered the Brownian kernel, kb (x, x′ ) = 1 +
min {x, x′ }, the Matérn kernel, km (x, x′ ) = (1 + 2|x − x′ |) exp(−2|x − x′ |), and the Gaussian
kernel, kg (x, x′ ) = exp(x − x′ )2 . For each (µ, γ) in a grid of values, the estimator fbµ,γ is defined
as the minimizer of the criteria given at Equation (17) taking ωv = ζv = 1 for all v ∈ P. In
order to save computation time, we restrict the optimisation to sets v such that |v| ≤ 3. Some
preliminary simulations showed that the terms coresponding to |v| ≥ 4 are nearly always equal
to 0.
Choosing the tuning parameters. Let us begin with the comparison of the two methods proposed
for choosing the final estimator, see Section 3.2. The results are given in Tables 1 and 2. It
appears that the procedure based on the ridge estimator after selection of the groups outperforms
the method based on the group-sparse estimator. As expected both methods perform better
when n increases, and when σ = 0.

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

15

σ=0
σ = 0.2
n = 50 n = 100 n = 200 n = 50 n = 100 n = 200
Proc. GS
1.91
0.79
0.45
2.41
1.16
0.54
0.10
0.03
1.50
0.47
0.15
Proc. rdg 0.80
Table 3. Estimated global error GE × 100 for different values of n and σ with
the Matérn kernel.

v = {1}
v = {2}
v = {3} v = {1, 2} v = {1, 3} v = {2, 3} v = {1, 2, 3} sum
S.I.
43.3
24.3
19.2
5.63
4.45
2.50
0.579
99.98
Proc. GS 50.1 (6.2) 26.5 (5.4) 20.9 (4.9) 0.69 (0.8) 0.63 (1.1) 0.51 (0.9) 0.02 (0.07) 99.29
Proc. rdg 45.4 (4.3) 25.3 (3.4) 20.3 (3.5) 3.08 (2.5) 2.18 (2.1) 1.44 (1.8)
0.09 (0.5) 98.66
Table 4. The first line of the table gives the true values of the Sobol indices
×100 greater than 10−2 , as well as their sum in the last columns. The following
lines give the mean of the estimators as well as their standard-error (in parenthesis), calculated over 100 simulations, for n = 50, and σ = 0 with the Matérn
kernel.

v = {1}
v = {2}
v = {3} v = {1, 2} v = {1, 3} v = {2, 3} v = {1, 2, 3} sum
S.I.
43.3
24.3
19.2
5.63
4.45
2.50
0.579
99.98
Proc. GS 47.5 (5.4) 26.2 (4.9) 19.8 (3.7) 2.35 (1.5) 1.45 (1.2) 0.84 (0.8)
0.03 (0.1) 98.95
0.05 (0.3) 99.02
Proc. rdg 43.5 (3.9) 25.0 (3.6) 19.7 (2.8) 4.85 (1.7) 3.39 (1.5) 2.02 (1.3)
Table 5. The first line of the table gives the true values of the Sobol indices
×100 greater than 10−2 , as well as their sum in the last columns. The following
lines give the mean of the estimators as well as their standard-error (in parenthesis), calculated over 100 simulations, for n = 100, and σ = 0.2 with the Matérn
kernel.

SI < ρ SI ≥ ρ v = {1} v = {2} v = {3} v = {1, 2} v = {1, 3} v = {2, 3} v = {1, 2, 3}
Proc. GS
17.9
68
100
100
100
72
66
67
9
6.3
51
100
100
100
72
62
47
4
Proc. rdg
Table 6. The first two columns give respectively pSelSv ≤ρ and pSelSv >ρ . The
last columns give the values of pSelv for each group v such that Sv > ρ. Results
for n = 50 and σ = 0 with the Matérn kernel.

Similarly, the Sobol indices are better estimated, in the sense of the global error, with the
procedure based on the ridge estimate of the metamodel, see Table 3. The means of the estimators for the Sobol indices greater than ρ = 10−4 are given in Tables 4 and 5. It appears
than S{1} is over-estimated using the procedure based on the group-sparse estimator, leading to
under-estimate the Sobol indices associated with interactions of order 2. This tendancy is much
less pronounced with the procedure based on the ridge estimator.
Let us now consider the performances of the procedure for selecting the non zero Sobol indices.
In Tables 6 and 7 we report the percentages of simulations for which the Sobol indices smaller
(respectively greater) than ρ are selected, and for which each of Sobol indice greater than ρ is
selected. From these results, we conclude that the procedure based on the ridge estimator is
more strict for selecting non-zero Sobol indices.

16

SYLVIE HUET AND MARIE-LUCE TAUPIN

SI < ρ SI ≥ ρ v = {1} v = {2} v = {3} v = {1, 2} v = {1, 3} v = {2, 3} v = {1, 2, 3}
Proc. GS
38
85
100
100
100
100
99
98
18
6
58
100
100
100
98
92
77
3
Proc. rdg
Table 7. The first two columns give respectively pSelSv ≤ρ and pSelSv >ρ . The
last columns give the values of pSelv for each group v such that Sv > ρ. Results
for n = 100 and σ = 0.2 with the Matérn kernel.
n = 50, σ = 0
Matérn Gaussian Brownian
Proc. GS 0.033
0.054
0.027
0.024
0.025
Proc. rdg 0.011
Table 8. Estimated empirical risk
cording to the kernel choice.

n = 100, σ = 0.2
mixed Matérn Gaussian Brownian mixed
0.028
0.033
0.054
0.027
0.028
0.011
0.011
0.024
0.025
0.023
ER: Performances of the procedures ac-

n = 50, σ = 0
n = 100, σ = 0.2
Matérn Gaussian Brownian mixed Matérn Gaussian Brownian mixed
Proc. GS
1.91
2.49
2.19
1.88
1.16
1.90
1.69
1.16
Proc. rdg
0.80
1.39
1.40
0.85
0.48
0.83
0.73
0.51
Table 9. Estimated global error GE × 100: Performances of the procedures
according to the kernel choice.

Comparing different kernels. Finally we compare the performances of the procedures for different
kernels, see Tables 8 and 9. The means of the estimated empirical risk, ER, and of the global
error for estimating the sensitivity indices, GE, are calculated for each kernel. It appears that
the Matérn kernel gives the best results, except for the case n = 50 and σ = 0 where the
emptirical risk of fbGS is smaller for the Brownian kernel.
In practice, one may want to choose the kernel according to the smallest prediction error.
For that purpose, we propose to calculate the estimators fbGS and/or fbrdg for each kernel, as
described in Section 3.2, as well as their associated prediction errors. Then we choose the kernel
for which the prediction error is minimized. The results are reported under the column “mixed”
in Tables 8 and 9. It appears that the estimated empirical risks for this “mixed” procedure are
nearly equal to the minimum estimated empirical risks over the different kernels.
5. Proofs
Recall that we cconsider the regression model defined at Equation (1), where X has distribution PX = P1 × . . . × Pd defined on X a compact subset of Rd , and ε is distributed as N (0, σ 2 ).
We denote by PX,ε the distribution of (X, ε). We observe a n sample (Yi , Xi ), i = 1, . . . , n with
law PX,ε .
The notation and the procedure are given in Sections 2.1 and 2.2.
Let us add on few notations that will be used along the proofs.
For v ∈ P we denote |v| the cardinal of v. For a function φ : R|v| 7→ R, we denote Vn,ε the
empirical process defined by
n
1X
(22)
εi φ(Xv,i ).
Vn,ε (φ) =
n
i=1

For the sake of simplicity we assume σ = 1. Moreover, we set R = R′ = 1, see (5) and (8).
Consequently, for any function f ∈ H, kfv kHv ≤ 1, and kfv k∞ ≤ kf kHv . The proofs can be

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

17

done exactly in the same way by considering the general case. In the proofs, the k · kL2 (PX ) norm
will be denoted by k · k2 .
5.1. Proof of Theorem 2.1. The proof is based on four main lemmas proved in Section 5.4.
In Section 5.3 other lemmas used all along the proof are stated. Their proof are postponed to
Section 5.5.
Let us first establish inequalities that will be used in the following. Let f ∈ H and v ∈ Sf
(see (7)).
Using that for any v ∈ Sf , and any norm k · k in Hv , kfv k − kfbv k ≤ kfv − fbv k and that for
any v ∈
/ Sf , kfv k = 0, we get that
X
X
X
X
µv kfbv kHv ,
(23)
µv kfv − fbv kHv −
µv kfv kHv −
µv kfbv kHv ≤
v∈P

(24)

v∈P

X

v∈P

γv kfv kn −

X

v∈P

v∈Sf

γv kfbv kn ≤

X

v∈Sf

v∈Sfc

γv kfv − fbv kn −

X

v∈Sfc

γv kfbv kn .

Combining (23), and (24), to the fact that for any function f ∈ H, L(fb) ≤ L(f ), we obtain
that
km − fbk2n ≤ km − f k2n + B

with

i
i
Xh
Xh

µv kfbv kHv + γv kfbv kn .
µv kfbv − fv kHv + γv kfbv − fv kn −2
(25) B = 2Vn,ε fb−f +2
v∈Sfc

v∈Sf

If km − f k2n ≥ B, we immediately get the result since in that case
X
X
γv2 .
µv +
km − fbk2n ≤ 2km − f k2n ≤ 2km − f k2n +
v∈Sf

v∈Sf

If km − f k2n < B, we get that
(26)
(27)

kfb − mk2n ≤ 2B

i
Xh

µv kfbv − fv kHv + γv kfbv − fv kn .
≤ 4|Vn,ε fb − f | + 4
v∈Sf


The control of the empirical process |Vn,ε fb − f | is given by the following lemma (proved in
Section 5.4.1, page 22).
Lemma 5.1. Let Vn,ε be defined in (22). For any f in F, we consider the event T defined as
n
o

(28) T = ∀f ∈ F, ∀v ∈ P, |Vn,ε fbv − fv | ≤ κλ2n,v kfbv − fv kHv + κλn,v kfbv − fv kn ,

where the quantities λn,v are defined by Equation (11) and where κ = 10 + 4∆. Then, for some
positive constants c1 , c2 ,
X
exp(−nc2 λ2n,v ).
PX,ε (T ) ≥ 1 − c1
v∈P

Conditionning on T , Inequality (27) becomes
i
i
Xh
Xh
µv kfbv − fv kHv + γv kfbv − fv kn ,
kfb − mk2n ≤ 4κ
λ2n,v kfbv − fv kHv + λn,v kfbv − fv kn + 4
v∈P

v∈Sf

18

SYLVIE HUET AND MARIE-LUCE TAUPIN

which may be decomposed as follows
X
X

[κλn,v + γv ] kfbv − fv kn +
kfb − mk2n ≤ 4
κλ2n,v + µv kfbv − fv kHv + 4
v∈Sf

4

X

v∈Sfc

κλ2n,v kfbv − fv kHv + 4

X

v∈Sfc

v∈Sf

κλn,v kfbv − fv kn .

If we choose C1 ≥ 2κ in Theorem 2.1, then κλ2n,v ≤ µv /2 and κλn,v ≤ γv /2 and the previous
inequality becomes
i
i
Xh
Xh
µv kfbv kHv + γv kfbv kn .
µv kfbv − fv kHv + γv kfbv − fv kn + 6
(29)
kfb − mk2n ≤ 6
v∈Sfc

v∈Sf

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

Lemma 5.2. For any f ∈ F, under the assumptions of Theorem 2.1 with C1 ≥ 2κ, conditionnally on T , see (28), we have
X
X
X
X
γv kfbv kn ≤ 3
µv kfbv kHv +
µv kfbv − fv kHv + 3
γv kfbv − fv kn .
v∈Sfc

v∈Sfc

v∈Sf

v∈Sf

Hence, by combining (29) and Lemma (5.2) we obtain
X

kfb − mk2n ≤ 24
µv kfbv − fv kHv + γv kfbv − fv kn .
v∈Sf

For each v, kfbv − fv kHv ≤ 2 (because the functions fbv et fv belong to the class F with R = 1,
see (5)), and consequently, for some constant C,



X
X
γv kfbv − fv kn .
µv +
(30)
kfb − mk2n ≤ C


v∈Sf

v∈Sf

P
P
To finish the proof it remains to compare the two quantities v∈Sf kfbv −fv k2n and k v∈Sf fbv −
P
P
fv k2n . For that purpose we show that k v∈Sf fbv − fv kn is less than k v∈Sf fbv − fv k22 plus an
additive term coming from concentration results (see the Lemma given below). Next, thanks to
P
P
the orthogonality of the spaces Hv with respect to L2 (PX ), k v∈Sf fbv −fv k22 = v∈Sf kfbv −fv k22 .
P
To conclude, it remains to consider several cases, according to the rankings of k v∈Sf fbv − fv k22 ,
P
k v∈Sf fbv − fv k2n , and d2n (f ). This is the subject of the following lemma whose proof is given
in Section 5.4.3, page 24.
Lemma 5.3. For f ∈ H, let A be the event
o
n
(31)
A = ∀f ∈ F, ∀v ∈ P, kfbv − fv kn ≤ 2kfbv − fv k2 + γv .

Then, for some positive constant c2 ,

PX,ε (A) ≥ 1 −

X
v

exp(−nc2 γv2 ).

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

19

On the set A, Inequality (30) provides that, for all K > 0

X
1 b
(32)
µv + 2γv kfbv − fv k2 + γv2
kf − mk2n ≤
C
v∈Sf

X
1 b
2
2
≤
µv + (1 + K)γv + kfv − fv k2 ,
(33)
K
v∈Sf

≤

(34)

X

v∈Sf


1 X b
µv + (1 + K)γv2 + k
fv − fv k22 .
K
v∈Sf

Inequality (33) uses the inequality 2ab ≤ K1 a2 + Kb2 for all positive K, and Inequality (34) uses
the orthogonality with respect to L2 (PX ).
P
In the following we have to consider several cases, according to the rankings of k v∈P fbv −
P
fv k2 , k v∈P fbv − fv kn and dn (f ) defined as in (12). More precisely, we consider three cases
P
P
Case 1: k v∈P fbv − fv k2 ≤ k v∈P fbv − fv kn .
P
P
Case 2: k v∈P fbv − fv kn ≤ k v∈P fbv − fv k2 ≤ dn (f )
P
P
P
fbv − fv kn ≤ k
fbv − fv k2 and dn (f ) ≤ k
fbv − fv k2 .
Case 3: k
v∈P

v∈P

v∈P

Case 1: From (34), for any f ∈ H, we get
X

1
1 b
kf − mk2n ≤
µv + (1 + K)γv2 + kfb − f k2n .
C
K
v∈Sf

Hence, using that for all K ′ > 0, kfb− f k2n ≤ (1 + K ′ )kfb− mk2n + (1 + 1/K ′ )kf − mk2n , we obtain
for a suitable choice of K ′ , say 1 + K ′ < K/C, that, for some positive constant C ′ ,


X
X
kfb − mk2n ≤ C ′ kf − mk2n +
µv +
γv2  .
v∈Sf

v∈Sf

This shows the result in Case 1.
Case 2: Inequality (34) becomes

X

1 b
1
kf − mk2n ≤
µv + (1 + K)γv2 + d2n (f ),
C
K
v∈Sf
o
nP
P
2
which gives the expected result since d2n (f ) = max
v∈Sf γv .
v∈Sf µv ,

P
P
P
Case 3: Recall that in this case, k v∈P fbv − fv kn ≤ k v∈S fbv − fv k2 and dn (f ) ≤ k v∈S fbv −
fv k2 . This case is solved by applying the following Lemma
(proved in Section 5.4.4, page 25),
√
b
b
which states that with high probability, kf − f k2 ≤ 2kf − f kn .
P
Lemma 5.4. Let f = v fv ∈ F with
P support Sf , dn (f ) be defined by (12), and let G(f ) be
the class of functions written as g = v∈P gv , such that kgv kHv ≤ 2 satisfying for all f ∈ F
X
X
X
X
4γv kgv kn
4µv kgv kHv +
C1
µv kgv kHv +
γv kgv kn ≤
v∈P

C2

X

v∈Sf

C3

v∈P

γv kgv kn ≤ 2

kgkn ≤ kgk2 .

X

v∈Sf

v∈Sf

v∈Sf

γv kgv k2 +

X

v∈Sf

γv2

20

SYLVIE HUET AND MARIE-LUCE TAUPIN

Then the event



kgk22
2
kgkn ≥
, kgk2 ≥ dn (f )
2
P
has probability greater than 1 − exp(−nκ v∈Sf λ2n,v ), where κ = 10 + 4∆ (see Lemma 5.1).

√
Note that Assumption nλ2n,v ≥ −C2 log(λn,v ) implies that λn,v = Kn,v / n with Kn,v → ∞.
P
Then, if f is such that |Sf | ≥ 1, exp(−nc3 v∈Sf λ2n,v ) ≤ exp(−c3 minv∈P Kn,v ). If f is such
that |Sf | = 0, then Condition C1 is not satisfied except if gv = 0 for all v ∈ P. Because we will
apply Lemma 5.4 to gv = fbv − fv , this event has probability 0. Therefore the event
)
(
2
X
kgk
2
, kgk22 ≥ dn (f )
(35) C = ∀f ∈ F, such that g =
(fbv − fv ) ∈ G(f ), and kgk2n ≥
2
v

has probability greater than 1 − η/3 for some 0 < η < 1.
P
Conditionning on the events T and A (defined by (28) and (31)), v (fbv − fv ) belongs to the
set G(f ). According to (34), we conclude in the same way as in Case 1.
Finally, it remains to quantify PX,ε (T ∩ A ∩ C). Following
Lemma 5.1, and Lemma 5.3,
P
T , respectively A, has probability greater than 1 − c1 v∈P exp(−nc2 λ2n,v ), respectively 1 −
P
2
2
v exp(−nγv ). Each of these probabilities is greater 1 − η/3 thanks to the assumption nλn,v ≥
−C2 log λn,v .
5.2. Proof of Corrolary 2.1. We start from Theorem 2.1 which states that with high probability,



X 
X
γv2 ,
µv +
kfb − mk2n ≤ C kf − mk2n +


v∈Sf

v∈Sf

and use that for all θ > 0,

(36)

1
kfb − mk22 ≤ (1 + θ)kfb − f k22 + (1 + )km − f k22 .
θ

For dn defined by (12), we consider once again the three cases defined page 19.
Case 1: According to (36), we get the result since
1
kfb − mk22 ≤ (1 + θ)d2n (f ) + (1 + )km − f k22 .
θ

Case 2: We directly obtain that

1
kfb − mk22 ≤ (1 + θ)d2n (f ) + (1 + )km − f k22 .
θ

P
P
P
Case 3: Recall that in this case, k v∈P fbv − fv kn ≤ k v∈P fbv − fv k2 and dn (f ) ≤ k v∈S fbv −
fv k2 . Apply Lemma 5.4 (page 19) and conclude that conditionning on the events T and A,
P
P
defined by (28) and (31), then v∈P fbv −fv belongs to G(f ) defined in Lemma 5.4 and v∈P fbv −
fv satisfyies Conditions C1-C3. Now, conditionning on the event C we get the result since
X
X
fbv − fv k2n .
fbv − fv k22 ≤ 2k
k
v∈P

v∈P



METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

21

5.3. Intermediate Lemmas. For v ∈ P, let Hv be the RKHS associated to the self reproducing
kernel kv . Let Qn,v and νn,v and be defined by Equations (9) and (10). For any function gv ∈ Hv ,
let Vn,ε be defined at Equation (22) and consider the following processes
(37)
(38)

Wn,2,v (t) = sup {|Vn,ε (gv )| , kgv kHv ≤ 2, kgv k2 ≤ t}
Wn,n,v (t) = sup {|Vn,ε (gv )| , kgv kHv ≤ 2, kgv kn ≤ t} .

Lemma 5.5. If EX,ε denotes the expectation with respect to the distribution of (X, ε), we have
for all t > 0,
EX,ε Wn,2,v (t) ≤ Qn,v (t).
Its proof is given in Section 5.5.1 page 28.
Lemma 5.6. Let b > 0 and let G(t) be the following class of functions:

(39)

G(t) = {gv ∈ Hv , kgv kHv ≤ 2, kgv k2 ≤ t, kgv k∞ ≤ b} .

Let Ωv,t be the event defined as


bt
(40)
.
Ωv,t = sup {|kgv k2 − kgv kn |, gv ∈ G(t)} ≤
2

Then for any t ≥ νn,v , the event Ωv,t has probability greater than 1 − exp(−c2 nt2 ), for some
positive constant c2 .
Its proof is given in Section 5.5.2, page 29.
Lemma 5.7. For any function gv ∈ Hv satisfying kgv kHv ≤ 2, kgv k∞ ≤ b and kgv k2 ≥ t, for
all t ≥ νn,v and b ≥ 1, the event




b
b
kgv k2 ≤ kgv kn ≤ 1 +
kgv k2
1−
2
2

has probabilty greater than 1 − exp(−c2 nt2 ) for some positive constant c2 .
Its proof is given in Section 5.5.3, page 30.

Lemma 5.8. If Eε denotes the expectation with respect to the distribution of ε, we have




nδ2
PX,ε Wn,n,v (t) − Eε Wn,n,v (t) ≥ δt ≤ 4 exp −
(41)
.
2
Its proof is given in Section 5.5.4, page 30.

Lemma 5.9. Conditionnaly on the space Ωv,t defined by (40), we have the two following
inequalities:




nδ2
(42)
PX,ε Wn,2,v (t) − Eε Wn,2,v (t) ≥ δt ≤ 4 exp −
8




nx2
(43)
.
PX Eε Wn,2,v (t) − EX,ε Wn,2,v (t) ≥ x ≤ exp −
Qn,v (t)
Its proof is given in Section 5.5.5, page 31.

Lemma 5.10. Let λn,v be defined at Equation (11), ∆ at Equation (10) and κ = 10 + 4∆.
Conditionnaly on the space Ωv,λn,v defined at Equation (40), for some positive constants c1 , c2 ,
with probability greater than 1 − c1 exp(−c2 nλ2n,v ), we have
(44)

Wn,n,v (λn,v ) ≤ κλ2n,v and Eε Wn,n,v (λn,v ) ≤ κλ2n,v .

Its proof is given in Section 5.5.6, page 32.

22

SYLVIE HUET AND MARIE-LUCE TAUPIN

5.4. Proofs of Lemma 5.1 to 5.4:
5.4.1. Proof of Lemma 5.1 (page 17). For f ∈ F and v ∈ P, let gv = fbv − fv . Note that
kgv kHv ≤ 2. Let us show that


(45)
|Vn,ε (gv )| ≤ κ λ2n,v kgv kHv + λn,v kgv kn .
We start by writing that

(46)

|Vn,ε (gv )| = kgv kHv Vn,ε



gv
kgv kHv



≤ kgv kHv Wn,n,v

 kg k 
v n
.
kgv kHv

Consider the two following cases:
Case A: kgv kn ≤ λn,v kgv kHv
Case B: kgv kn > λn,v kgv kHv
Case A: Since kgv kn ≤ λn,v kgv kHv , we have
Wn,n,v

 kg k 
v n
≤ Wn,n,v (λn,v ).
kgv kHv

We then apply Lemma 5.10, page 21, and conclude that (45) holds in Case A for each v ∈ P
since, with high probability
(47)

|Vn,ε (gv )| ≤ κλ2n,v kgv kHv ≤ κλ2n,v kgv kHv + κλn,v kgv kn .

Case B: Consider now the case kgv kn > λn,v kgv kHv and let us show that for any v ∈ P,


kgv kn
≤ κλn,v kgv kn .
Wn,n,v
kgv kHv
Let rv be a deterministic number such that rv > λn,v . Our first step relies on the study of the
process Wn,n,v (rv ), for rv > λn,v . In that case we state two results:
R1 For any deterministic rv ≥ λn,v , with probability greater than 1 − c1 exp(−c2 nλ2n,v ),
(48)

Wn,n,v (rv ) ≤ (1 + κ)rv λn,v .

R2 Inequality (48) continues to hold for random rv of the form
rv =

kgv kn
.
kgv kHv

Combining these two points implies that, with probability greater than 1 − c1 exp(−c2 nλ2n,v ),
kgv kHv Wn,n,v

 kg k 
v n
≤ κkgv kn λn,v .
kgv kHv

Consequently, in Case B, according to (26), for each v, Inequality (45) holds because
|Vn,ε (gv )| ≤ κkgv kn λn,v ≤ κλ2n,v kgv kHv + κλn,v kgv kn .
This ends up the proof of Lemma 5.1.

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

23

Proof of R1. Taking t = rv and δ = λn,v in (41), with probability greater than 1−4 exp(−nλ2n,v ),
we have
Wn,n,v (rv ) ≤ Eε [Wn,n,v (rv )] + rv λn,v .

Next we prove that for some positive rv , with probability greater than 1 − exp(−ncλ2n,v ), we
have
Eε Wn,n,v (rv ) ≤ κrv λn,v .

(49)

Let νbn,v defined as the smallest solution of Eε [Wn,n,v (t)] ≤ κt2 . For Wn,n,v , defined by (37), we
write
rv
Eε Wn,n,v (rv ) =
Eε sup {|Vn,ε (gv )|, kgv kHv ≤ νbn,v /rv , kgv kn ≤ νbn,v }
νbn,v
rv
rv
Eε Wn,n,v (b
νn,v ) ≤
κb
ν 2 = κrv νbn,v .
≤
νbn,v
νbn,v n,v

Besides, Lemma 5.10 stated that on the event Ωv,λn,v , Eε Wn,n,v (λn,v ) ≤ κλ2n,v . It follows from
the definition of νbn,v
P, and 2Lemma 5.6, that νbn,v ≤ λn,v for all v ∈ P with probability greater
than 1 − exp(−nc2 v∈P λn,v ). Consequently, for any deterministic rv such that rv ≥ λn,v , (48)
is satisfied with high probability.
Proof of R2. Let us prove R2 by using a peeling-type argument. Our aim is to prove that (48)
holds for any rv of the form
kgv kn
.
rv =
kgv kHv

Since kgv k∞ /kgv kHv ≤ 1, we have kgv kn /kgv kHv ≤ 1. We thus restrict ourselves to rv satisfying
rv = kgv kn /kgv kHv with kgv kn /kgv kHv ∈ (λn,v , 1].
We start by splitting the interval (λn,v , 1] into M disjoint intervals such that (λn,v , 1] =
k−1 λ
k
c
∪M
n,v , 2 λn,v ], for some M that will be chosen later. Consider the event D defined as
k=1 (2
follows:
kg v kn
∈ (λn,v , 1]}.
D c = {∃v ∈ P and ∃g v , such that |Vn,ε (g v )| ≥ κλn,v kg v kn , with
kg v kH
We prove that, for some positive constants c1 , c2 ,

P [D c ] ≤ c1 exp(−c2 nλ2n,v ).
For g v ∈ D c , let k be the integer in {1, · · · , M }, such that
2k−1 λn,v ≤

kg v kn
≤ 2k λn,v .
kg v kH⊑

This k satisfies


 kg k 
v n
≥ kg v kHv Wn,n,v
≥ |Vn,ε (g v )| ≥ κλn,v kg v kn .
kg v kHv Wn,n,v 2k λn,v
kgv kHv

Therefore, we get

Wn,n,v (2k λn,v ) ≥ κλn,v

kg v kn
λn,v k
≥ κλ2n,v 2k−1 ≥ κ
2 λn,v .
kg v kHv
2

By taking rv = 2k λn,v in (48), we have


λn,v k
k
2 λn,v ≤ c1 exp(−c2 nλ2n,v ).
P Wn,n,v (2 λn,v ) ≥ κ
2

24

SYLVIE HUET AND MARIE-LUCE TAUPIN

Now let us write D c as follows:

[ 
kg v kn
kg v kn
c
k−1
k
D =
with
∈ (2 λn,v , 2 λn,v ) .
∃v, ∃ g v such that|Vn,ε (g v )| ≥ κλn,v
kg v kHv
kg v kH
k=1,··· ,M

The set D c has probability smaller than c1 M exp(−c2 nλ2n,v ). If we choose M such that log M ≤
(c2 /2)nλ2n,v , then the probability of the set T is greater than
X
c2
1−
c1 exp − nλ2n,v .
2
v∈P

It follows that R2 is proved which ends up the proof of Lemma 5.1.



5.4.2. Proof of Lemma 5.2 (page 18). Starting from (26) with B defined by Equation (25), we
write
1 b
kf − mk2n ≤ 2|Vn,ε (fb − f v)| +
2
X
X
[µv kfbv − fv kHv + γv kfbv − fv kn ] − 2
[µv kfbv kHv + γv kfbv kn ].
2
v∈S c

v∈Sf

On the event T defined in (28) we have
X
X
1 b
kf − mk2n ≤ 2κ
λ2n,v kfbv − fv kHv + 2κ
λn,v kfbv − fv kn +
2
v∈P
j∈P
X
X
[µv kfbv − fv kHv + γv kfbv − fv kn ] − 2
[µv kfbv kHv + γv kfbv kn ].
2
v∈S c

v∈Sf

Rearranging the terms we obtain that
X
X
1 b
(2κλn,v + 2γv )kfbv − fv kn +
(2κλ2n,v + 2µv )kfbv − fv kHv +
kf − mk2n ≤
2
v∈Sf
v∈Sf
X
X
(2κλn,v − 2γv )kfbv kn .
(2κλ2n,v − 2µv )kfbv kHv +
v∈Sfc

v∈Sfc

Now, thanks to Assumption (13) with C1 ≥ 2κ we have 2κλ2n,v ≤ µv and 2κλn,v ≤ γv and
Lemma 5.2 is shown since
X
X
1
0 ≤ kfb − mk2n ≤ 3
µv kfbv − fv kHv + 3
kfbv − fv kn −
2
v∈Sf
v∈Sf
X
X
γv kfbv kn .
µv kfbv kHv −
v∈Sfc

v∈Sfc



5.4.3. Proof of Lemma 5.3 (page 18): Let us consider the following two cases:
• kfbv − fv k2 ≤ γv . We apply Lemma 5.6 (page 21) to the function gv = fbv − fv . It satisfies
gv ∈ G(γv ) with b = 2 (recall that k·k∞ ≤ k·kHv ). Moreover, γv ≥ C1 λn,v ≥ C1 νnv ≥ νn,v
as soon as C1 ≥ 1.
It follows that, for some positive c2 , with probability greater than 1 − exp(−nc2 γv2 )
kfbv − fv kn ≤ kfbv − fv k2 + γv

• kfbv − fv k2 ≥ γv . We apply Lemma 5.7 (page 21) to the function gv = fbv − fv with b = 2.
It follows that, for some positive c2 , with probability greater than 1 − exp(−nc2 γv2 ),
kfbv − fv kn ≤ 2kfbv − fv k2 .

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

25


5.4.4. Proof of Lemma 5.4 (page 19): Let dn (f ) be defined by (12), and let G(f ) and G ′ (f ) be
the following sets
(
)
X
G(f ) =
g=
gv , satisfying kgv kHv ≤ 2, and Conditions C1, C2, C3 ,
v∈P

′

G (f ) = {g ∈ G(f ), such that kgk2 = dn (f )} .

For the sake of simplicity we remove the dependency in f in the notation. Let us consider the
two events B and B ′ defined as follows:
B ′ = {∀h ∈ G ′ , khk2n ≥ d2n /2} and B = {∀h ∈ G, khk2n ≥ khk22 /2, and khk2 ≥ dn }.

Let us first remark that B ′ is included into B: if h ∈ B ′ , then h ∈ G, khk2 = dn and khk2n ≥ d2n /2.
It follows that khk2n ≥ khk22 /2 and khk2 ≥ dn . Therefore Lemma 5.4 is proved if B ′ holds with
high probability. Consider
Zn (G ′ ) = sup {d2n − kgk2n }.
g∈G ′

has probability greater than 1 − c1 exp(−nc2 d2n ).
We show that the event Zn
≤
Let us briefly recall the notion of covering numbers for a totally bounded metric space (G, ρ),
consisting of a set G and a metric ρ defined from G × G into R+ . A δ-covering set of G is a
collection of functions f 1 , · · · , f N such that for all f ∈ G there exists k ∈ {1, 2, · · · , N } such that
ρ(f, f k ) ≤ δ. The δ-covering number N (δ, G, ρ) is the cardinality of the smallest δ-covering set.
A propercovering restricts the covering to use only elements in the set G. The proper covering
number denoted Npr (δ, G, ρ) satisfies
(G ′ )

(50)

d2n /2

N (δ, G, ρ) ≤ Npr (δ, G, ρ) ≤ N (δ/2, G, ρ).

Let us now consider a dn /8-covering of (G ′ , k · kn ), so that, for all g in G ′ there exists g k such
that kg − g k kn ≤ dn /8. The associated proper covering number is
Npr = Npr (dn /8, G ′ , k · kn ).

(51)

Now, for all g ∈ G ′ , T1 = kgk k2n − kgk2n , and T2 = d2n (f ) − kg k k2n , we write
d2n − kgk2n = T1 + T2 .

The proof is splitted into four steps:
Step 1 The first step consists in showing that

d2n
.
4
Step 2 The second step consists in proving that, for Npr given at Equation (51) and for some
constant C,




2
k 2
2
PX
max [dn − kg kn ] ≥ dn /4 ≤ exp log Npr − Cnd2n .
T1 = kg k k2n − kgk2n ≤

(52)

k∈{1,··· ,Npr }

Step 3 The third step concerns the control of Npr : we show the following result
!2
64
Eε sup |Vn,ε (g)| .
log Npr ≤ n
dn g∈G ′
Step 4 The last step consists in bounding from above the Gaussian complexity:
Eε sup |Vn,ε (g)| ≤
g∈G ′

20κ 2
d .
C1 n

26

SYLVIE HUET AND MARIE-LUCE TAUPIN

Let us conclude the proof of the lemma before proving these four steps.
Putting together Steps 3 and 4, if C1 is large enough, say C1 > 64 × 20κ/C, then Step 2 states
that






d2n
2
k 2
2
PX T2 ≥
max [dn − kg kn ] ≥ dn /4 ≤ exp − κnd2n .
≤ PX
4
k∈{1,··· ,Npr }
Now, we have




d2n
≤ exp − κnd2n .
4
g 1 ,··· ,g N
P
We conclude the proof of the lemma, by noting that d2n ≥ C12 v λ2n,v (see (11), (12) and (13)).


PX Zn (G ′ ) ≤ d2n /2 = PX



max {d2n − kgk k2n } ≥

Proof of Step 1: We start by writing that

i
1 Xh k
=
(g (Xi ))2 − (g(Xi ))2
n
i=1
v
u n
u1 X
k
≤ kg − gkn t
[gk (Xi ) + g(Xi )]2 .
n
n

kgk k2n

−

kgk2n

i=1

Using that (a + b)2 ≤ 2a2 + 2b2 , g ∈ G ′ , and g satisfies Condition C3, we get
n

1 X k (i)
[g (X ) + g(X (i) )]2 ≤ 2kgk k2n + 2kgk2n ≤ 4d2n .
n
i=1

Besides, the covering set is constructed such that kg k − gkn ≤ dn /8. It follows that Step 1 is
proved.
Proof of Step 2: We prove that for some constant C,






d2n
nd2n (f )
d2n
2
k 2
PX T2 ≥
≤ PX max {dn − kg kn } ≥
≤ exp log Npr − C
.
1≤k≤Npr
4
4
1 + dn (f ) + d2n (f )
As gk ∈ G ′ , dn = kgk k2 . Then
o
n
d2n − kgk k2n =
max
1≤k≤Npr

max [kgk k22 − kg k k2n ].

1≤k≤Npr

Applying Theorem 3.5 in Chung and Lu [9] we have that for all positive λ
#
" n


X
λ2
k
2
k
2
(g (Xi )) ≤ nE(g (Xi )) − λ ≤ exp −
PX
.
2nE(gk (X))4
i=1

Taking λ = nd2n /4 and using that kg k k22 = d2n we get




nd4n
d2n
2
k 2
.
≤ exp −
PX {dn − kg kn } ≥
4
32E(gk (X))4

It follows that

 X
Npr




d2n
nd4n
nd4n
2
k 2
PX max
(53)dn − kg kn ≥
≤
exp
log
N
−
.
≤
exp −
pr
1≤k≤Npr
4
32E(gk (X))4
32 maxk E(gk (X))4
k=1

It remains to calculate EX g4 (X) for g ∈ G ′ . Precisely we show the following result:

Eg 4 (X) ≤ 2d2n 2 + 11d2n + 4dn .

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

27

P
This result comes from the property of the RKHS H: indeed g ∈ H is written g = v∈P gv
where the functions gv are centered and orthogonal in L2 (PX ). Therefore Eg 4 (X) is the sum of
the following terms:
 X
X
4
4
A1 =
EX gv (Xv ),
EX gv2 (Xv )gv2′ (Xv′ ),
A2 =
2
v∈P
v6=v′
  X
  X
4
4
2
A3 =
EX gv1 (Xv1 )gv2 (Xv2 )gv3 (Xv3 ), A4 =
EX gv31 (Xv1 )gv2 (Xv2 )
3
3
v1 6=v2 6=v3
v1 6=v2
 
X
4
EX gv1 (Xv1 )gv2 (Xv2 )gv3 (Xv3 )gv4 (Xv4 ).
A5 =
1
v1 6=v2 6=v3 6=v4

Using the Cauchy-Schwartz inequality and the fact that kgv k∞ ≤ kgv kHv ≤ 2, and kgk2 = dn
(because g ∈ G ′ ), we get that A1 is proportionnal to d2n , A2 , A3 , A5 to d4n , and A4 to d3n . For
example
X
A1 =
EX gv4 (Xv ) ≤ kgv k2∞ kgv k22 ≤ 4d2n .
v∈P

After calculation of the terms Ai , since d2n (f ) is assumed to be smaller than one, we get that

(54)

max EX (gk (X))4 ≤ 4d2n (1 + dn + (11/2)d2n ) ≤ 34d2n .
k

Step 2 is proved by combining (53) and (54).
Proof of Step 3: Let Npr be defined at Equation (51). We prove that
r
64
log Npr
≤
Eε sup |Vn,ε (g)|.
n
dn
g∈G ′
We start from (50) and write that
log Npr (dn /8, G ′ , k · kn ) ≤ log N (dn /16, G ′ , k · kn ).
Using the Sudakov minoration (see Pisier [24]) we have that for all positive ω
"
#
√
p
n
4
log N (ω, G ′ , k · kn ) ≤
Eε sup |Vn,ε (g)| .
ω
g∈G ′
Hence by taking ω = dn /16, Step 3 is proved.

Proof of Step 4: The last step consists in bounding from above the Gaussian complexity Eε supg∈G ′ |Vn,ε (g)|.
This control is performed by using Lemma 5.5 (page 21). According to Inequality (45),
#
"
X
X
λ2n,v kgv kHv +
λn,v kgv kn ,
|Vn,ε (g)| ≤ κ
v∈P

v∈P

with λn,v defined by Equation (11) satisfying C1 λn,v ≤ γv and C1 λ2n,v ≤ µv for all v ∈ P. It
follows
#
"
X
X
X
2
λn,v kgv kHv +
λn,v kgv kn
sup
|Vn,ε (gv )| ≤ κ sup
g∈G ′ v∈P

g∈G ′

≤

v∈P

"

v∈P

X
X
κ
sup
µv kgv kHv +
γv kgv kn
C1 g∈G ′
v∈P

v∈P

#

28

SYLVIE HUET AND MARIE-LUCE TAUPIN

and according to Condition C1,
sup

X

g∈G ′ v∈P

|Vn,ε (gv )| ≤

4κ
C1

≤

4κ
C1

sup

X

g∈G ′ v∈S

2

X

µv kgv kHv + sup

X

g∈G ′ v∈S

µv + sup

X

g∈G ′ v∈S

v∈S

γv kgv kn

!

γv kgv kn

!

,

,

because kgv kHv ≤ 2. Now, according to Condition C2, and using that 2ab ≤ a2 + b2 , we get
"
#
X
X
X
X
4κ
2
2
2
µv + 2
γv + sup
kgv k2
sup
|Vn,ε (gv )| ≤
′
C1
g∈G
g∈G ′ v∈P
v∈S
v∈S
v∈S
#
"
X
4κ
≤
2
(µv + γv2 ) + d2n ,
C1
v∈S
P
the last inequality coming from the fact for all g ∈ G ′ , kgk22 = d2n ≥ v∈S kgv k22 .
Finally, thanks to (12), we get
X
20κ 2
d .
sup
|Vn,ε (gv )| ≤
C1 n
g∈G ′
v∈P

5.5. Proofs of intermediate Lemmas.
5.5.1. Proof of Lemma 5.5 (page 21): Let us write that the kernel kv is written as :
X
kv (xv , yv ) =
ωv,k φv,k (xv )φv,k (yv )
k

is an orthonormal basis of L2 (Pv ), where Pv =
where
Let us consider the class of functions K(t) defined as
(φv,k )∞
k=1

X
i

ai φv,i , with kgv k2Hv =

In the following, we set µk,v (t) = min
(55)

a∈v

Pa .

K(t) = {gv ∈ Hv , kgv kHv ≤ 2, kgv k2 ≤ t} .

It comes that
gv =

Q

X
k



X
X a2
i
≤ 4, and kgv k22 =
a2i ≤ t2
ωv,i
i

i

t2 , ω

k,v

. Hence

a2k
1 X 2 X a2k
1
≤ 2
= 2 kgv k22 + kgv k2Hv ≤ 5,
ak +
µk,v (t)
t
ωk,v
t
k

k

as soon as gv ∈ K(t).
Now, let us prove the lemma:

n

EX,ε Wn,2,v (t) = EX,ε sup

g∈K(t)

1X X
εi
aℓ φv,ℓ (xv,i )
n
i=1

ℓ

n

= EX,ε sup

g∈K(t)

X
a
1X
p ℓ
εi
n
µv,ℓ (t)

v
u
X
√ u
5tEX,ε
≤
ℓ

i=1

ℓ

n

1X
εi
n
i=1

q

q

µv,ℓ (t)φv,ℓ (xv,i )
!2

µv,ℓ (t)φv,l (xv,i )

.



METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

29

The last inequality follows from the Cauchy-Schwartz inequality and Inequality (55). Now,
simple calculation leads to
s
√
1X
µv,ℓ (t),
EX,ε Wn,2,v (t) ≤ 5
n
ℓ


p
√
√
5.5.2. Proof of Lemma 5.6 (page 21): Using that
a − b ≤ |a − b|, we get
q
|kgv k2 − kgv kn | ≤
kgv k22 − kgv k2n .

Hence

 


b2 t 2
bt
2
2
⊂ kgv k2 − kgv kn ≥
.
kgv k∞ ≤ b, |kgv k2 − kgv kn | ≥
2
4
The centered process
n

kgv k22

−

kgv k2n

1X 2
gv (Xv,i ) − E(gv2 (Xv )) ,
n

=

i=1

satisfies a concentration inequality given, for example, by Theorem 2.1 in Bartlett et al. [2] : if
C is a class of functions f such that kf k∞ ≤ B and Ef (X) = 0, and if there exists γ > 0 such
that for every f ∈ C, Var f (X) ≤ γ 2 . Then for every x > 0, with probability at least 1 − e−x ,


 r


 
n
n

X
X
1
1 x
1
2x
1
f (Xj ) ≤ inf 2(1 + α)E sup
f (Xj )  +
γ+B
+
.
(56) sup
α>0 
n
3 α n
f ∈C n
f ∈C n
j=1

j=1

For any t > 0, for G(t) defined by (39), let us consider the class of functions C(t) defined as
follows

C(t) = f such that f = gv2 − E(gv2 ), with gv ∈ G(t) .

Note that if f ∈ C(t), EX f (Xv ) = 0 and kf k∞ ≤ b2 . We have to study
γ 2 (t) = sup EX
gv ∈G(t)

It is easy to see that

2
gv2 (X) − kgv k22 ) and Γ(t) = EX

sup
gv ∈G(t)

kgv k2n − kgv k22

!

.

γ 2 (t) ≤ b2 sup EX (gv (X) + kgv k2 )2 ≤ 4b2 t2 .
gv ∈G(t)

Let ζi be independent and identically random variables Rademacher distributed and let EX,ζ
denotes the expectation with respect to the law of (X, ζ). By a symmetrization argument,
n

Γ(t) ≤ 2EX,ζ sup

gv ∈G(t)

1X 2
ζi gv (Xi ) .
n
i=1

Since kgv k∞ ≤ b, applying the contraction principal (see Ledoux-Talagrand [19]) we get that,
for Qn,v (t) defined by (9),
n

EX,ζ sup

gv ∈G(t)

1X 2
ζi gv (Xi )
n
i=1

n

≤ 4bEX,ζ sup

gv ∈G(t)

≤ 4bQn,v (t).

1X
ζi gv (Xi )
n
i=1

30

SYLVIE HUET AND MARIE-LUCE TAUPIN

The last inequality was proved by Mendelson [23], Theorem 41 (see the proof of Lemma 5.5).
Now, thanks to (56) we get that for all x > 0, with probability greater than 1 − e−x
(
r

 )
1 x
2x
2 1
2
2
.
2bt + b
+
sup kgv kn − kgv k2 ≤ inf 16(1 + α)bQn,v (t) +
α>0
n
3 α n
gv ∈G(t)
2

Taking x = c2 nt2 , t ≥ νn , we have that with probability greater than 1 − e−c2 nt
 


√
1
2
2
2
2 1
+
sup kgv kn − kgv k2 ≤ inf t 16(1 + α)b∆ + 2c2 4b + b
c2 .
α>0
3 α
gv ∈G(t)
p
The infimum of the right hand side is reached in α = c2 b/16∆, and equals
p
√
b2 c2
+ 8 ∆c2 b3/2 + 4(4∆ + 2c2 )b.
3

The constants ∆ and c2 should satisfy that this infimum is strictly smaller than b2 /4. For
example, if 16∆ < b/8, it remains to choose c2 small enough such that
√ 

√
b
c2
2c2
b
+
+ 4 2c2 < .
3
2
8

5.5.3. Proof of Lemma 5.7 (page 21): Let t > νn,v and h be defined as h = tgv /kgv k2 . If
gv satisfies the assumptions of the lemma, then h satisfies khk2 = t, khkH ≤ 2 and khk∞ ≤ b.
Applying Lemma 5.6 (page 21) to the function h, we obtain that for all t ≥ νn,v , with probability
greater than 1 − exp(−c2 nt2 ), we have |t − khkn | ≤ bt/2 for all h ∈ G(t). This concludes the
proof of the lemma.
5.5.4. Proof of Lemma 5.8 (page 21): The proof of Lemma 5.8 is based on an isoperimetric
inequality for Gaussian processes (Borell [5] or Cirel’son et al. [10]) as it is stated in Theorem
(3.8), page 61 in Massart [21]. Let us recall this inequality:
Lemma 5.11. Let P be the Gaussian probability measure on Rn and let φ be a function from
Rn to R, and kφkL its Lipschitz semi-norm:
|f (x) − f (y)|
kφkL = sup √
.
nkx − ykn
x6=y
Let Φ be the cumulative distribution of the standard Gaussian distribution. Then for any u,


u
(57)
.
P (|f − EP f | ≥ u ) ≤ 4Φ
kφkL
We apply Lemma 5.11 to φ(ε1 , . . . , εn ) = Wn,n,v (t). By Cauchy-Schwarz Inequality, kφkL =
√
t/ n. It follows that Lemma 5.8 is proved since







(δt)2 
nδ2
PX,ε (|Wn,n,v (t) − Eε Wn,n,v (t)| ≥ δt) ≤ 4 exp −  2 ≤ 4 exp −
.


2
 2 √t

n



METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

31

5.5.5. Proof of Lemma 5.9 (page 21): We start with the proof of (42) in Lemma 5.9 by applying
once again Lemma 5.11 given above, to the function φ(ε) = φ(ε1 , . . . , εn ) = Wn,2,v (t). On the
event Ωv,t defined by (40), we have kgv kn ≤ bt/2+kgv k2 . Besides if kgv kHv ≤ 2, then kgv k∞ ≤ 2.
Therefore applying Lemma 5.6 with b = 2, we get that if kgv k2 ≤ t,
φ(ε) − φ(ε′ ) ≤ kgv kn kε − ε′ kn ≤ 2tkε − ε′ kn ,
√
leading to kφkL = 2t/ n. It follows that (42) in Lemma 5.9 is proved since









nδ2
(δt)2 
c
.
PX,ε {|Wn,2,v (t) − Eε Wn,2,v (t)| ≥ δt} ∩ Ωv,t ≤ 4 exp −  2 ≤ 4 exp −

8

 2 √2t 
n

We now come to the proof of (43) in Lemma 5.9 using a Poissonian inequality for self-bounded
processes (see Boucheron et al. [6]) and Theorem 5.6, p 158 in Massart [21]). Let us recall it in
the particular case we are interested in:
Theorem 5.1. Let X1 , · · · , Xn be n iid variables. For i ∈ {1, · · · , n} let X(−i) = (X1 , . . . , Xi−1 , Xi+1 , . . . , Xn ).
Let Z be a nonnegative and bounded measurable function of X = (X1 , · · · , Xn ). Assume that
for all i ∈ {1, · · · , n}, there exists a measurable function Zi of X(−i) such that 0 < Z − Zi ≤ 1,

P
and ni=1 (Z − Zi ) ≤ Z. Then, for all x > 0, we have P {Z ≥ E(Z) + x} ≤ exp −x2 /2E(Z) .
We apply this result to Z defined as

Z = Z(X1 , · · · , Xn ) = nEε Wn,2,v (t) = nEε sup {|Vn,ε (gv )|, kgv k2 ≤ t, kgv kHv ≤ 2} .

The variable Z is positive, and because the distribution of (ε1 , . . . , εn ) is symmetric, we have
that
Z = Eε sup {nVn,ε (gv ), kgv k2 ≤ t, kgv kHv ≤ 2} .
Let τ be the function in Hv such that Z = Eε nVn,ε(τ ) (note that τ depends on (X1 , . . . , Xn )
and on (ε1 , . . . , εn )), and let
X
εj gv (Xj ).
Zi = Eε sup
gv

j6=i

We show that Z and Zi satisfy the assumptions of Theorem 5.1:


X
X
Z − Zi = Eε εi τ (Xi ) +
εj gv (Xj )
εj τ (Xj ) − sup
gv

j6=i

j6=i

r
2
1
,
≤ Eε εi τ (Xi ) ≤ √ Eε sup |τ (X)| ≤
π
2π
x∈X
where the last inequality comes from the fact that supx∈X |τ (X)| ≤ kτ kHv ≤ 2. Moreover
Z − Zi ≥ 0 since


n
n
X
X
εj gv (Xj ) = Eε Eεi sup
εj gv (Xj )
Z = Eε sup
gv



Finally we have:

gv

j=1

≥ Eε sup Eεi
gv

n
X
j=1



j=1

εj gv (Xj ) = Zi .



n
n
n
n
X
X
X
X
X


Eε εi τ (Xi ) +
(Z − Zi ) =
εj τ (Xj ) − sup
εj gv (Xj ) ≤
Eε εi τ (Xi ) = Z.
i

i=1

j6=i

g−v

j6=i

i=1

32

SYLVIE HUET AND MARIE-LUCE TAUPIN

Therefore, following Theorem 5.1, we get that for all postive u


h
u2
ui
≤ exp −
.
PX,ε Eε Wn,2,v (t) − EX,ε Wn,2,v (t) ≤
n
EX,ε Wn,2,v (t)

As EX,ε Wn,2,v (t) ≤ Qn,v (t), see Lemma 5.5 page 21, we get the expected result since for all
positive x


nx2
.
PX [Eε Wn,2,v (t) ≥ EX,ε Wn,2,v (t) + x] ≤ exp −
Qn,v (t)



5.5.6. Proof of Lemma 5.10 (page 21): From Lemma 5.8, page 21 with t = λn,v = δ, we get that
with probability greater than 1 − 4 exp(−nλ2n,v /2),
Eε Wn,n,v (λn,v ) ≤ λ2n,v + Wn,n,v (λn,v ).

The next step consists in comparing Wn,n,v (λn,v ) and Wn,2,v (2λn,v ). Recall that λn,v ≥ νn,v ,
see (11). Let gv such that kgv kn ≤ λn,v .
• When kgv k2 ≤ λn,v , according to Lemma 5.6 (page 21), taking b = 2 , since since
kgv kn ≤ λn,v , we get that with probability greater than 1 − exp(−c2 nλ2n,v ),
kgv kn − λn,v ≤ kgv k2 ≤ kgv kn + λn,v ≤ 2λn,v .

• When kgv k2 ≥ t, we apply Lemma 5.7 (page 21) with b = 2. For any function gv such
that kgv k∞ ≤ 2, and kgv k2 ≥ λn,v , we have kgv k2 ≤ 2kgv kn ≤ 2λn,v .
This implies that, with probability greater than 1 − exp(−c2 nλ2n,v ) we have
Wn,n,v (λn,v ) ≤ Wn,2,v (2λn,v ).

We now study the process Wn,2,v (λn,v ). By applying (42) in Lemma 5.9, page 21, with
δ = t = λn,v we get that with probability greater than 1 − 4 exp(−nλ2n,v /8)
It follows that

Wn,2,v (λn,v ) ≤ λ2n,v + Eε (Wn,2,v (λn,v )).

Eε Wn,n,v (λn,v ) ≤ λ2n,v + Wn,n,v (λn,v ))

≤ λ2n,v + Wn,2,v (2λn,v ))

≤ 5λ2n,v + Eε (Wn,2,v (2λn,v )).

Next, we apply (43) in Lemma 5.9, with t = 2λn,v and x = 4λ2n,v . We get that
Eε Wn,2,v (2λn,v ) ≤ 4λ2n,v + EX,ε (Wn,2,v (2λn,v )),

with probability greater than

nλ4n,v
1 − 2 exp −16
Qn,v (2λn,v )

!

4nλ2n,v
≥ 1 − 2 exp −
∆

!

.

The last inequality comes from the definition of νn,v , see (10), and from the fact that λn,v ≥ νn,v ,
see (11).
Putting everything together, we get that with probability greater than 1 − c1 exp(−c2 nλ2n,v )
for some positive constants c1 , c2 ,
Eε Wn,n,v (λn,v ) ≤ 9λ2n,v + EX,ε (Wn,2,v (2λn,v ))

≤ 9λ2n,v + Qn,v (2λn,v ), thanks to Lemma 5.5, page 21,
≤ 9λ2n,v + 4∆λ2n,v .

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

33

Applying once again Lemma 5.8, page 21, we get that



Wn,n,v (λn,v ) ≤ Eε Wn,n,v (λn,v ) + λ2n,v ≤ 10 + 4∆ λ2n,v .

This ends the proof of the lemma by taking κ = 10 + 4∆.


5.6. Proof of Propositions 5.1-5.4. We consider the minimization of C ′ (f0 , θ) given at Equation (17). Because C ′ (f0 , θ) is convex and separable, we use a block coordinate descent algorithm,
group v by group v. We refer to Bubeck [8] for a review on convex optimization.
In what follows, the group v is fixed, and for given values of f0 and θ w , w 6= v, we look for
the minimizer of C ′ with respect to θ v : Setting
X
Rv = Y − f 0 −
Kw θ w ,
w6=v

we aim at minimizing with respect to θ v

′

′

Q(θ v ) = kRv − Kv θ v k2 + γv kKv θ v k + µv kKv1/2 θ v k.

If γv′ = µ′v = 0, then θ v = Kv−1 Rv is the solution. In what follows we consider the case where
at least one of both is non zero.
If ∂Qv denotes the subdifferential of Q(θ v ) with respect to θ v , we need to solve 0 ∈ ∂Qv . Let
us recall that for all v ∈ P the matrices Kv are symmetric and strictly definite positive.
Let us begin with the calculation of the subdifferential of kKv θv k with respect to θv : if θv 6= 0,
we have
Kv2 θv
∂kKv θv k
=
,
∂θv
kKv θv k
and if θv = 0 the subdifferential is the set of x ∈ Rn such that kKv−1 xk ≤ 1.
Therefore if θv 6= 0,
∂Qv = −2Kv Rv + 2Kv2 θv + γv′

Kv θv
Kv2 θv
,
+ µ′v
1/2
kKv θv k
kKv θv k

while if θv = 0,
n
o
∂Qv = −2Kv Rv + γv′ t + µ′v s where t, s ∈ Rn such that kKv−1 tk ≤ 1, kKv−1/2 sk ≤ 1 .
Let us begin with the case µ′v = 0.

Proposition 5.1. Let (ρ)+ denotes the positive part of ρ ∈ R. If µ′v = 0,


γv′
θv = 1 −
K −1 Rv .
k2Rv k + v
Proof of Proposition 5.1. The problem comes to minimize
′

U (β v ) = kRv − β v k2 + γv kβ v k

and to take θ v = Kv−1 β v . The subdifferential of U is given by
β
∂Uv (β v ) = −2Rv + 2β v + γv′ v if β v =
6 0
kβ v k

∂Uv (0) = −2Rv + γv′ t where ktk ≤ 1 .

We get

βv =



γv′
1−
k2Rv k



Rv .
+



34

SYLVIE HUET AND MARIE-LUCE TAUPIN

Let µ′v > 0, the following proposition gives a necessary and sufficient condition for θ v = 0 to
be the minimizer of Q. For the sake of clarity let us recall the definitions of J and J ∗ given in
Section 3.1: For all x ∈ Rn
n
o
J(x) = k2Rv − µ′v Kv−1 xk2 , and J ∗ = min J(x), for x such that kKv−1/2 xk ≤ 1 .
Proposition 5.2. The minimizer of Q(θ v ) is θ v = 0 if and only if J ∗ ≤ γv′2 .

Proof of Proposition 5.2.
(1) Let us assume that θ v = 0 is the minimizer of Q(θ v ). Then 0 ∈ ∂Q implies that there
−1/2 ∗
s k ≤ 1 and such that γv′ t∗ + µ′v s∗ =
exists t∗ and s∗ such that kKv−1 t∗ k ≤ 1 and kKv
2Kv Rv .
If γv′ > 0, then

1p
1
Kv−1 t∗ = ′ 2Rv − µ′v Kv−1 s∗ , and kKv−1 t∗ k = ′ J(s∗ ).
γv
γv

Because J ∗ ≤ J(s∗ ) and kKv−1 t∗ k ≤ 1, we get that J ∗ ≤ γv′2 .
If γv′ = 0, then µ′v s∗ = 2Kv Rv and J(s∗ ) = J ∗ = 0.
(2) Let us now assume that J ∗ ≤ γv′2 . Note that minimizing the convex function J(s) over
−1/2
the convex set kKv
sk ≤ 1 has always a solution. Let us denote this solution by s∗ .
′
′
∗
If γv > 0, let t = (2Kv Rv − µ′v s∗ )/γv . Then −2Kv Rv + γv′ t∗ + µ′v s∗ = 0, and
′
−1/2 ∗
s k ≤ 1, and kKv−1 t∗ k = J(s∗ )/γv ≤ 1.
−2Kv Rv + γv′ t∗ + µ′v s∗ ∈ ∂Q(0) since kKv
Therefore θ v = 0 is the minimizer of Q.
−1/2 ∗
s k ≤ 1. Therefore
If γv′ = 0, and J ∗ = 0, then −2Kv Rv + µ′v s∗ = 0, and kKv
θ v = 0 is the minimizer of Q.


Proposition 5.3. Let µ′v > 0 and θ v be the minimizer of Q.
1/2

(1) If k2Kv Rv k ≤ µ′v , then θ v = 0
(2) If not, for ρ > 0, let
−1
θ(ρ) = 2µ′v µ′2
v Kv + ρIn

−1

Kv−1/2 Rv ,

and let ρ∗ defined as kθ(ρ∗ )k = 1. Then J(θ(ρ∗ )) ≤ γv′2 if and only if θ v = 0,
−1/2

Proof of Proposition 5.3. Minimizing J(x) under the constraint kKv
xk ≤ 1 is equivalent to
−1/2
1/2
′
2
2
minimizing K(β) = k2Rv − µv Kv
βk under the constraint kβk ≤ 1. Let β ∗ = 2Kv Rv /µ′v .
∗
∗
Then K(β ) = 0, which is smaller than γv′2 , and if kβ k ≤ 1, following Proposition 5.2, we get
θ v = 0.
If kβ ∗ k > 1, we have to solve a ridge regression problem by minimizing K(β) + ρkβk2 for
some positive ρ. The solution is given by θ(ρ). Let us note that kθ(ρ)k decreases to 0 when ρ
tends to infinity and that its maximum is kθ(0)k = kβ ∗ k. Therefore if kβ ∗ k > 1, there exists
ρ∗ such that kθ(ρ∗ )k = 1. Following Proposition 5.2, Proposition 5.3 is proved. A numerical
procedure can be used for calculating this ρ∗ .

Let us now consider the case where θ v is non zero. It should satisfy the subgradient condition
∂Qv = 0 which leads to
!−1
2
Kv θv
′
2
′ Kv θv
+ µv
2Kv Rv
(58)
θv = 2Kv + γv
1/2
kKv θv k
kKv θv k
which is equivalent to Equation (21).

METAMODEL CONSTRUCTION FOR SENSITIVITY ANALYSIS

35

Proposition 5.4. For all ρ1 , ρ2 > 0 let
θ(ρ1 , ρ2 ) = (Kv + ρ1 Kv + ρ2 In )−1 Rv .
If µ′v > 0, there exists a non zero solution to Equation (21) if and only if there exists ρ1 , ρ2 > 0
such that

γv′ = 2ρ1 kKv θ(ρ1 , ρ2 )k
(59)
1/2
µ′v = 2ρ2 kKv θ(ρ1 , ρ2 )k
Then θ v = θ(ρ1 , ρ2 ).

Proof of Proposition 5.4. If there exists a non zero solution to Equation (21) then kKv θ v k
1/2
and kKv θ v k are non zero because Kv is definite positive. Let ρ1 = γv′ /2kKv θ v k and ρ2 =
1/2
µ′v /2kKv θ v k, then
!−1
γv′
µ′v
θ(ρ1 , ρ2 ) = Kv +
Kv +
In
Rv = θ v ,
1/2
2kKv θ v k
2kKv θ v k

and, for such ρ1 , ρ2 , Equation (59) is satisfied.
Conversely, if there exist ρ1 , ρ2 such that Equation (59) is satisfied, then necessarily kKv θ(ρ1 , ρ2) )k
1/2

1/2

and kKv θ(ρ1 , ρ2) )k are non zero and ρ1 = γv′ /2kKv θ(ρ1 , ρ2 )k and ρ2 = µ′v /2kKv θ(ρ1 , ρ2 )k.
Then
!−1
µ′v
γv′
In
Rv ,
Kv +
θ(ρ1 , ρ2 ) = Kv +
1/2
2kKv θ(ρ1 , ρ2 )k
2kKv θ(ρ1 , ρ2 )k

which is exactly Equation (21) calculated in θ v = θ(ρ1 , ρ2 ).



References
[1] Francis Bach. High-Dimensional Non-Linear Variable Selection through Hierarchical Kernel Learning. working paper or preprint, September 2009.
[2] Peter L. Bartlett, Olivier Bousquet, and Shahar Mendelson. Local Rademacher complexities. Ann. Statist.,
33(4):1497–1537, 2005.
[3] A. Berlinet and C. Thomas-Agnan. Reproducing Kernel Hilbert Spaces in Probability and Statistics. Springer
US, 2003.
[4] G. Blatman and B. Sudret. Adaptive sparse polynomial chaos expansion based on least angle regression.
Journal of computational Physics, 230:2345–2367, 2011.
[5] Christer Borell. The Brunn-Minkowski inequality in Gauss space. Invent. Math., 30:207–216, 1975.
[6] Stéphane Boucheron, Gabor Lugosi, and Pascal Massart. A sharp concentration inequality with applications.
Random Structures & Algorithms, 16(3):277–292, 2000.
[7] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and
statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn., 3(1):1–
122, 2011.
[8] S. Bubeck. Convex optimization: Algorithms and complexity. volume 8 of Foundation and trends in Machine
Learning, pages 231–357. now, 2015.
[9] Fan Chung and Linyuan Lu. Concentration inequalities and martingale inequalities: a survey. Internet Math.,
3(1):79–127, 2006.
[10] B.S. Cirel’son, I.A. Ibragimov, and V.N. Sudakov. Norms of Gaussian sample functions. Proc. 3rd JapanUSSR Symp. Probab. Theory, Taschkent 1975, Lect. Notes Math. 550, 20-41 (1976)., 1976.
[11] N. Durrande, D. Ginsbourger, O. Roustant, and L. Carraro. {ANOVA} kernels and {RKHS} of zero mean
functions for model-based sensitivity analysis. Journal of Multivariate Analysis, 115(0):57 – 67, 2013.
[12] Jerome H. Friedman. Multivariate adaptive regression splines. Ann. Statist., 19(1):1–141, 1991. With discussion and a rejoinder by the author.
[13] L. Gu and F.J. Wu. A unified framework for uncertainty and sensitivity analysis of computational models
with many input. In SIMUL 2014 : The Sixth International Conference on Advances in System Simulation.
IARIA XPS Press, 2014.
[14] SR Gunn and JS Kandola. Structural modelling with sparse kernels. MACHINE LEARNING, 48(1-3):137–
163, 2002.

36

SYLVIE HUET AND MARIE-LUCE TAUPIN

[15] Bertrand Iooss and Paul Lemaı̂tre. A review on global sensitivity analysis methods. In C. Meloni and
G. Dellino, editors, Uncertainty management in Simulation-Optimization of Complex Systems: Algorithms
and Applications. Springer, 2015.
[16] K. Kandasamy and Y. Yu. Additive approximations in high dimensional nonparametric regression via the
salsa. In Proceedings of the 33rd International Conference on Machine Learning, New-York, volume 48.
JMLR: W & CP, 2016.
[17] V. Koltchinskii and M. Yuan. Sparse recovery in large ensembles of kernel machines. In Proceeding of the
21st Annual Conference on Learning Theory, Helsinki, pages 229–238. 2008.
[18] Vladimir Koltchinskii and Ming Yuan. Sparsity in multiple kernel learning. The Annals of Statistics,
38(6):3660–3695, 12 2010.
[19] M. Ledoux and M. Talagrand. Probability in Banach Spaces: Isoperimetry and Processes. Ergebnisse der
Mathematik und ihrer Grenzgebiete. U.S. Government Printing Office, 1991.
[20] Yi Lin and Hao Helen Zhang. Component selection and smoothing in multivariate nonparametric regression.
Ann. Statist., 34(5):2272–2297, 2006.
[21] P. Massart. Concentration inequalities and model selection, volume 1896 of Lecture Notes in Mathematics.
Springer, Berlin, 2007. Lectures from the 33rd Summer School on Probability Theory held in Saint-Flour,
July 6–23, 2003, With a foreword by Jean Picard.
[22] Lukas Meier, Sara van de Geer, and Peter Bühlmann. High-dimensional additive modeling. The Annals of
Statistics, 37(6B):3779–3821, 12 2009.
[23] Shahar Mendelson. Geometric parameters of kernel machines. In Computational learning theory (Sydney,
2002), volume 2375 of Lecture Notes in Comput. Sci., pages 29–43. Springer, Berlin, 2002.
[24] Gilles Pisier. The volume of convex bodies and Banach space geometry, volume 94 of Cambridge Tracts in
Mathematics. Cambridge University Press, Cambridge, 1989.
[25] R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for
Statistical Computing, Vienna, Austria, 2008. ISBN 3-900051-07-0.
[26] Garvesh Raskutti, Martin J. Wainwright, and Bin Yu. Minimax-optimal rates for sparse additive models over
kernel classes via convex programming. J. Mach. Learn. Res., 13:389–427, 2012.
[27] Pradeep Ravikumar, John Lafferty, Han Liu, and Larry Wasserman. Sparse additive models. J. R. Stat. Soc.
Ser. B Stat. Methodol., 71(5):1009–1030, 2009.
[28] A. Saltelli, K. Chan, and E.M. (eds) Scott. Sensitivity analysis. J. Wiley & Sons, 2000.
[29] Andrea Saltelli, Marco Ratto, Terry Andres, Francesca Campolongo, Jessica Cariboni, Debora Gatelli,
Michaela Saisana, and Stefano Tarantola. Global sensitivity analysis: The primer. Chichester: John Wiley &amp; Sons. , 2008.
[30] I. M. Sobol. Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp., 1(4):407–
414, 1993.
[31] I.M Sobol. Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates.
Mathematics and Computers in Simulation, 55(13):271 – 280, 2001. The Second {IMACS} Seminar on Monte
Carlo Methods.
[32] Samir Touzani and Daniel Busby. Smoothing spline analysis of variance approach for global sensitivity analysis
of computer codes. RELIABILITY ENGINEERING & SYSTEM SAFETY, 112:67–81, APR 2013.
[33] A. W. van der Vaart. Asymptotic statistics. Cambridge series in statistical and probabilistic mathematics.
Cambridge University Press, Cambridge (UK), New York (N.Y.), 1998. Autre tirage : 2000 (dition broche),
2005, 2006, 2007.
[34] Grace Wahba. Spline models for observational data, volume 59 of CBMS-NSF Regional Conference Series in
Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1990.
[35] Grace Wahba, Yuedong Wang, Chong Gu, Ronald Klein, and Barbara Klein. Smoothing spline ANOVA for
exponential families, with application to the Wisconsin Epidemiological Study of Diabetic Retinopathy. Ann.
Statist., 23(6):1865–1895, 1995.
Unit MaIAGE, INRA Jouy-en-Josas, France;
E-mail address: sylvie.huet@inra.fr
Laboratoire LaMME, UMR CNRS 8071- USC INRA, Université d’Evry Val d’Essonne, France
E-mail address: marie-luce.taupin@univ-evry.fr

