Gabriel Sarazin
email: gabriel.sarazin@cea.fr

Amandine Marrel
email: amandine.marrel@cea.fr

Sebastien Da Veiga
email: sebastien.da-veiga@ensai.fr

Vincent Chabridon
email: vincent.chabridon@edf.fr

New insights into the feature maps of Sobolev kernels: application in global sensitivity analysis

Keywords: 2020 Mathematics Subject Classification. 34A30, 46E35, 47G10, 62G05, 62G10 Sensitivity analysis, Hilbert-Schmidt independence criterion, cross-covariance operators, Mercer kernels, kernel feature maps, kernel integral operators, ordinary differential equations

come

Submitted November 29, 2023.

Introduction

Kernel methods

Kernels, to be understood here as symmetric and positive definite functions of two arguments, are essential tools in probability and statistics. In probability, they are rather called covariance functions and they are used in stochastic calculus to characterize the properties of random fields and processes [START_REF] Hristopulos | Random Fields for Spatial Data Modeling[END_REF]. In statistics, they are the backbone of numerous non-linear methods designed to carry out machine learning tasks, whether supervised (e.g. Gaussian process regression [START_REF] Williams | Gaussian Processes for Machine Learning, volume 2 of Adaptive Computation and Machine Learning series[END_REF], support vector classification [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF], denoising with smoothing splines [START_REF] Wahba | Spline Models for Observational Data[END_REF][START_REF] Wang | Smoothing Splines: Methods and Applications[END_REF]) or unsupervised (e.g. spectral clustering [START_REF] Bach | Learning spectral clustering[END_REF], kernel k-means [START_REF] Dhillon | Kernel k-means: spectral clustering and normalized cuts[END_REF], kernel PCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF]). Kernels have also enabled the development of advanced methods in areas such as novelty detection [START_REF] Campbell | Kernel methods: a survey of current techniques[END_REF][START_REF] Campbell | A linear programming approach to novelty detection[END_REF], Bayesian quadrature [START_REF] Bach | On the equivalence between kernel quadrature rules and random feature expansions[END_REF][START_REF] Briol | Frank-Wolfe Bayesian quadrature: probabilistic integration with theoretical guarantees[END_REF], independent component analysis [START_REF] Bach | Kernel independent component analysis[END_REF], hypothesis testing (especially independence testing [START_REF] Gretton | A kernel statistical test of independence[END_REF] or two-sample testing [START_REF] Gretton | A kernel two-sample test[END_REF]), optimal experimental design [START_REF] Mutny | Experimental design for linear functionals in reproducing kernel Hilbert spaces[END_REF], optimal quantization [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF] and many others. The ubiquity of kernels in statistical applications is mainly explained by their ability to introduce non-linearity where standard methods are only able to handle linearity. From a mathematical viewpoint, a kernel can be described in two different but equivalent ways.

First, in virtue of the Moore-Aronszajn theorem [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], every kernel is related to a unique reproducing kernel Hilbert space (RKHS) which is generated from the only knowledge of the kernel and inherits most of its properties (separability, measurability, boundedness, continuity, differentiability) [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF][START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF]. In particular, the link between a kernel and its RKHS is summarized by a reproducing property, which often leads to speak of reproducing kernels (instead of kernels). The RKHS offers a sound theoretical framework where some challenging problems find easy solutions. A classic example is provided by kernel ridge regression [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF] where a penalized least-square score has to be minimized over a high-dimensional class of candidate functions [START_REF] Gu | Smoothing Spline ANOVA Models[END_REF][START_REF] Wahba | Spline Models for Observational Data[END_REF]. If an RKHS is chosen, the representer theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF] indicates that any minimizer can be expressed as a weighted linear combination of functions picked from the RKHS.

Then, beyond the abstract existence of the RKHS, there is always a transformation (called feature map) that allows to rewrite the kernel as the inner product in a specific Hilbert space (called feature space) between the representatives (called feature functions) of the two initial arguments [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF]. The role played by a kernel within an algorithm is sometimes easier to understand by adopting the feature viewpoint. A typical example is the case of binary classification with support vector machines [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. In presence of non-separable data (i.e. when linear boundaries fail to separate the two classes), the idea is to use a kernel in order to transport the data into a feature space where they become linearly separable. The knowledge of the kernel feature map then allows to clearly identify the transformations applied to the data to make them separable.

Global sensitivity analysis [START_REF] Veiga | Basics and Trends in Sensitivity Analysis: Theory and Practice in R[END_REF][START_REF] Saltelli | Global Sensitivity Analysis: the Primer[END_REF] is another area where kernel methods overcome the limitations of the standard methodology and where the feature maps allow to understand how kernels perform the expected task.

Global sensitivity analysis

In many industrial fields, reliability and risk assessment is based on the joint use of computer simulation and uncertainty management. For highly sophisticated systems, the physical phenomena involved in accidental scenarios are modeled by numerical simulators which compute one (or several) output(s) of interest from a large number of uncertain input parameters. The generic Monte Carlo approach for uncertainty propagation relies on input-output evaluations of the simulator. As a complement, global sensitivity analysis (GSA) ambitions to quantify the influence of input uncertainties on the output uncertainty [START_REF] Borgonovo | Sensitivity Analysis: an Introduction for the Management Scientist[END_REF][START_REF] Saltelli | Sensitivity Analysis in Practice: a Guide to Assessing Scientific Models[END_REF]. When computer-based experiments are time-consuming, the maximum feasible number of runs may be very limited. The main challenge of GSA is then to deliver accurate sensitivity estimates from the few available simulation data.

The most prevailing approach to perform GSA is to apportion the output variance between all subsets of input variables [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]. Variance allocation is notably guided by the ANOVA1 decomposition of the multivariate function representing the simulator [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF]. The resulting Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] are very popular because their definition (as percentages of the output variance) is very intuitive. However, the computational burden required to accurately estimate them is often prohibitive. To bypass this problem, it is often preferred to adopt a pairwise approach based on the estimation and comparison of dependence measures [START_REF] Borgonovo | A new uncertainty importance measure[END_REF][START_REF] Borgonovo | A common rationale for global sensitivity measures and their estimation[END_REF][START_REF] Borgonovo | Sensitivity analysis: a review of recent advances[END_REF]. This can be achieved in many different ways, but kernel methods are certainly the most attractive option, as they notably allow to define the Hilbert Schmidt Independence Criterion (HSIC) [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF]. Computing the HSIC between two random variables (each previously assigned an appropriate kernel) amounts to using a dissimilarity measure to compare their joint distribution and their hypothetical distribution under independence (within an RKHS built from the two selected kernels). Interestingly, and unlike the numerator of high-order Sobol' indices, the HSIC can be estimated accurately, even from a limited number of samples. Hence, to meet the needs of GSA when simulation data are expensive, the HSIC can be applied to all input-output pairs [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF], and the resulting sensitivity measures are called HSIC indices. Among their many advantages, HSIC indices are often praised for their ability to characterize independence because, assuming that all kernels are characteristic, an HSIC index is equal to zero if and only if the output variable is independent from the input variable.

Despite all this, HSIC indices suffer from a lack of interpretability, and this hinders their dissemination to a wider audience. There are actually two problems. Firstly, as HSIC indices do not arise from an ANOVA-like decomposition, their sum is not equal to one, and they cannot be universally upper bounded (i.e. any bound can be questioned with exotic choices of kernels and/or distributions). Secondly, the nature of the information captured by HSIC indices is not directly accessible, unless adopting the cross-covariance viewpoint [START_REF] Baker | Joint measures and cross-covariance operators[END_REF][START_REF] Vakhania | Probability Distributions on Banach Spaces[END_REF]. A prerequisite to understand finely the way HSIC indices measure sensitivity is the prior knowledge of the feature maps nested in the kernels chosen for the input and output variables. Many long-established results may be found in the dedicated literature [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF][START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF][START_REF] Kanagawa | Gaussian processes and kernel methods: a review on connections and equivalences[END_REF][START_REF] Novak | Reproducing kernels of Sobolev spaces on R d and applications to embedding constants and tractability[END_REF], especially for Gaussian and Matérn kernels. However, for a given kernel, there is no guarantee that a fully analytical feature map can be easily found. This must be discussed on a case-by-case basis, which makes the interpretation of HSIC indices very dependent on the selected kernels.

New generation of HSIC indices based on Sobolev kernels

In response to the first interpretation problem, an ANOVA-like decomposition has been recently set up for HSIC indices [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF]. This HSIC-ANOVA decomposition allows to define the sensitivity indices of the same name. A contribution (defined in the HSIC style) can be assigned to each subset of inputs, and the sum of all contributions is then equal to one. This is very similar to Sobol' indices, except that all HSIC-ANOVA contributions can be estimated from a small amount of input-output data.

Unfortunately, the HSIC-ANOVA decomposition is obtained at the cost of strong restrictions on the input variables (assumed independent) and the input kernels (subject to additional orthogonality constraints). The input and output kernels must also be characteristic in order to preserve the ability to detect independence. Finding kernels that satisfy all these constraints is a major obstacle to the implementation of the HSIC-ANOVA methodology. The solution promoted in [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] is to apply a preliminary transformation to the input variables so that they all follow the standard uniform distribution2 , and then to assign them some specific Sobolev kernels (namely the unanchored Sobolev kernels).

Sobolev kernels have already proved to be efficient for spline regression [START_REF] Craven | Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF][START_REF] Gu | Smoothing Spline ANOVA Models[END_REF][START_REF] Wahba | Spline Models for Observational Data[END_REF]. They are also of great theoretical use to demonstrate the (strong) tractability of some quasi-Monte Carlo techniques in very high dimension [START_REF] Dick | Liberating the weights[END_REF][START_REF] Novak | Intractability results for integration and discrepancy[END_REF][START_REF] Sloan | When are quasi-Monte Carlo algorithms efficient for high dimensional integrals[END_REF]. Still from a theoretical perspective, they may be used to construct Hilbert spaces of multivariate functions in which the ANOVA components can be computed with simplified formulas [START_REF] Kuo | On decompositions of multivariate functions[END_REF][START_REF] Wasilkowski | Finite-order weights imply tractability of linear multivariate problems[END_REF]. However, in the context of GSA, the use of Sobolev kernels raises a dilemma. On the one hand, Sobolev kernels are necessary to compute HSIC-ANOVA indices (which may be seen as a more interpretable collection of HSIC indices). On the other hand, this gain in interpretability is counterbalanced by a loss of transparency regarding the way sensitivity is measured. Therefore, the main objective of this paper is to investigate the mathematical properties of Sobolev kernels, especially those which have immediate consequences in GSA. For this, several points deserve to be clarified.

(Q1) We need to identify at least one orthonormal basis (ONB) for each Sobolev RKHS. The knowledge of the basis functions is indeed the key to understanding the preliminary transformations applied to the input variables when HSIC-ANOVA indices are computed with Sobolev kernels.

(Q2) We need to extract at least one easily interpretable feature map for each Sobolev kernel. Of all possible feature maps, those arriving in ℓ 2 -spaces are particularly attractive because they often lead to an ONB of the RKHS.

(Q3) We need to know if Sobolev kernels are characteristic. This will determine the ability of HSIC-ANOVA indices to characterize independence.

Organization of the document

This paper is divided into seven sections which can be split into two parts. The first part, consisting of Sections 2 to 4, introduces all useful concepts before explaining why the acquisition of new knowledge about Sobolev kernels gives even more credit to the HSIC-ANOVA approach.

• Section 2 provides the essential ideas of the theory of reproducing kernels. The notion of feature map is introduced and its importance in the understanding of kernel action is emphasized. In particular, Section 2.4 presents two different strategies that can be used to identify an ONB of the RKHS associated to a given kernel. The first one (based on Mercer features) is well-known whereas the second one (based on ℓ 2 -linearly independent features) is a first contribution of this work.

• Section 3 says a few words about Sobolev kernels. Their links with Sobolev spaces are detailed and their most relevant properties are presented. The question (Q3) is addressed at this step because it is ultimately quite simple after introducing all the necessary concepts. In addition, the unanchored Sobolev kernels are introduced as a specific family of Sobolev kernels built upon Bernoulli polynomials and parameterized by an integer r ≥ 1 (accounting for the order of smoothness in the associated Sobolev space). Since the unanchored Sobolev spaces are among the few kernels fulfilling all the requirements of the HSIC-ANOVA approach, the study will be restricted to them.

• Section 4 is then focused on kernel-based GSA. The advantages of HSIC indices (over Sobol' indices) are highlighted. The added value of the newly-developed HSIC-ANOVA framework is justified and the problems caused by the use of Sobolev kernels are clearly stated. In this section, the information captured by HSIC-based sensitivity measures is thoroughly examined from the perspective of kernel feature maps. It is notably explained how basic manipulations on kernels sometimes allow to rewrite HSIC indices as series of explicit covariance terms.

In the second part of this work, Sections 5 to 8 strive to answer the questions (Q1) and (Q2). These four sections are therefore the core contribution of this work.

• In Section 5, a simulation-based method is applied to Sobolev kernels in order to get a first idea of the feature map arising from Mercer's theorem.

• In Section 6, the eigenvalue problem (solved numerically in Section 5) is transformed into a boundary value problem. In particular, it is shown that this equivalent version of the original problem can be solved analytically only for r = 1.

• In Section 7, a sub-kernel decomposition of Sobolev kernels is used to disclose a feature map of a different type which is a little less informative than Mercer's but holds whatever is r ≥ 2.

• Section 8 deals with the limit Sobolev kernel, obtained in the hypothetical situation where r = ∞.

Section 9 provides some concluding remarks. The key findings are summarized and they are reformulated in terms of the original questions. Throughout this work, the proofs which are sufficiently short are inserted in the body of the text. On the contrary, the long and more technical ones are reported to Appendix F.

Basic reminders on kernels

In this first section, the objective is to introduce briefly the few elements from the theory of reproducing kernels that will be used in this work. Section 2.1 is dedicated to recall some basic definitions. Then, Section 2.2 deals with the particular links between kernels and probability distributions, in particular through the notions of orthogonal kernels and kernel mean embeddings. In Section 2.3, the fundamental concept of feature map is put forward and some examples are provided. The specific case of Mercer kernels is carefully examined in Section 2.4 as it will be of great importance in the forthcoming developments.

Kernels and their reproducing kernel Hilbert spaces

In this work, for the sake of simplicity, X will always denote an interval in R but most definitions and results can be extended to locally compact spaces in R d (with d ∈ N * ) or to even more general spaces. In the following, R X denotes the space of all real-valued functions defined on X . Definition 2.1. A function K : X × X → R is said to be a kernel if it is symmetric and positive definite.

• Symmetry: ∀ x, x ′ ∈ X , K(x, x ′ ) = K(x ′ , x) ,
• Positive definiteness:

∀ n ≥ 1, ∀ x 1 , . . . , x n ∈ X , ∀ c 1 , . . . , c n ∈ R, n i=1 n j=1 c i c j K(x i , x j ) ≥ 0 .
After adopting a matrix viewpoint, the positive definiteness property may be rewritten as follows:

c T K n c ≥ 0 with K n := [K(x i , x j )] 1≤i,j≤n and c := [c i ] 1≤i≤n .
The n-by-n matrix K n is called the Gram matrix (or more simply the kernel matrix).

Definition 2.2. Let H be a Hilbert space 3 (of functions h : X → R) with inner product ⟨•, •⟩ H . For any x ∈ X , the evaluation functional at x is the linear form defined by

L x : h ∈ H → h(x) ∈ R.
Then, H is said to be an RKHS if all evaluation functionals are continuous:

∀ x ∈ X , ∃ C x > 0 such that ∀ h ∈ H, |h(x)| ≤ C x ∥h∥ H .
Definition 2.3. Let K : X × X → R be a kernel and let H be an RKHS (of functions h : X → R). K is said to be a reproducing kernel of H if the two following conditions are verified:

• Embedding property:

∀ x ∈ X , K(•, x) ∈ H , • Reproducing property: ∀ x ∈ X , ∀ h ∈ H, h(x) = ⟨h, K(•, x)⟩ H .
Remark 2.4. In Definition 2.3, the two points are crucial. Together, they notably lead to the fundamental property verified by K: In fact, any RKHS admits a reproducing kernel and any kernel is reproducing for an RKHS. More precisely, there is a one-to-one mapping between kernels and RKHSs. For a given RKHS, the Riesz-Fréchet representation theorem [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] (see Theorem 5.5, p. 135) applied to all evaluation functionals allows to construct a reproducing kernel and it is the only one for this RKHS. Conversely, for a given kernel K, the linear span of the functions K(•, x) is a pre-Hilbert which can be completed into a Hilbert space. By uniqueness of the completion, this Hilbert space is the only RKHS having K as reproducing kernel [START_REF] Aronszajn | Theory of reproducing kernels[END_REF].

∀ x, x ′ ∈ X , K(x, x ′ ) = ⟨K(•, x) , K(•, x ′ )⟩ H . (2.1)
Theorem 2.5 (Moore-Aronszajn). For any kernel K : X × X → R, there is a unique Hilbert space for which K is a reproducing kernel.

A detailed proof is provided in [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] (see Theorem 3,p. 19).

Remark 2.6. If K 1 and K 2 are two kernels on X , it is straightforward to see that K 1 + K 2 is also a kernel on X . The RKHS induced by K 1 + K 2 is not so simple to characterize [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] (see Theorem 5, p. 24). In the particular case where H 1 ∩ H 2 = {0}, it is exactly the direct sum H 1 ⊕ H 2 . If K 1 and K 2 are defined on two (different) intervals X 1 and X 2 , it can be proved 4 that the tensor product K 1 ⊗ K 2 is a kernel on X 1 × X 2 . The related RKHS (denoted by H 1 ⊗ H 2 ) is obtained after completing the linear span of all the functions h 1 ⊗ h 2 where h 1 ∈ H 1 and h 2 ∈ H 2 [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] (see Theorem 13, p. 31).

Links between kernels and distributions

For any interval X ⊆ R, the space of all Borel 5 probability measures on X is denoted by M + 1 (X ).

Orthogonal kernels

Definition 2.7. A kernel k : X × X → R is said to be orthogonal 6 with respect to ν ∈ M + 1 (X ) if one has:

∀ x ∈ X , X k(ξ, x) dν(ξ) = 0 . (2.2)
For fixed ν ∈ M + 1 (X ), let k be an orthogonal kernel and let F be the associated RKHS. By definition, all the functions k(•, x) with x ∈ X have zero mean. Due to the denseness of the pre-Hilbert space F pre := Span({k(•, x) with x ∈ X }) in F, the zero-mean property can actually be extended to any function f ∈ F. As a result, F does not contain any non-zero constant function, which may be written as R ∩ F = {0}. In this context, R does not represent the real line but the space of all constant functions on X . This space is an RKHS, and more precisely the RKHS induced by the constant kernel (x, x ′ ) → 1.

Definition 2.8. A kernel

K : X × X → R is said to be ANOVA (with respect to ν) if K = 1 + k where k : X × X → R is an orthogonal kernel (with respect to ν). Let K = 1 + k be an ANOVA kernel. As it is known that R ∩ F = {0}, one can write H = R ⊕ F.
For most widespread parametric families of distributions, no ANOVA kernel can be found in the literature. The only exception is the standard uniform distribution. In this particular case, some Sobolev kernels naturally satisfy the orthogonality constraint. This point will be further discussed in Section 3. 4 This directly follows from the Schur product theorem [START_REF] Styan | Hadamard products and multivariate statistical analysis[END_REF] (see Theorem 3.1, p. 221) which guarantees that the Hadamard product of two Gram matrices is indeed a symmetric positive semi-definite matrix. 5 Let B(X ) be the Borel σ-algebra of X , i.e. the σ-algebra generated by the open sets of X . Any probability measure ν defined on B(X ) is called a Borel probability measure on X . 6 In order to establish a clear distinction with other kernels, the letters k and F will be used (instead of K and H) for denoting orthogonal kernels and their RKHSs.

Embedding in L p -spaces

Let K : X × X → R be a kernel and let ν ∈ M + 1 (X ). If the random variable K(X, X) with X ∼ ν verifies some finite-moment assumptions, it can easily be shown that H ⊂ R X is contained in an L p -space related to ν:

∀ p ≥ 1, E ν K(X, X) p /2 < ∞ =⇒ H ⊂ L p (X , ν) . (2.3) Thus, it is enough to assume that E ν [K(X, X)] < ∞ for H to be part of L 2 (X , ν). In particular, if K is a bounded kernel, H is contained in any space L p (X , ν) with p ≥ 1 and ν ∈ M + 1 (X ).

Kernel embedding of probability distributions

Under the only assumption that E ν K(X, X) < ∞, the kernel mean embedding µ ν ∈ H of the probability distribution ν can be defined as:

∀ x ∈ X , µ ν (x) = E ν [K(X, x)] = X K(ξ, x) dν(ξ) .
(2.4)

The embedding mechanism allows to define a dissimilarity measure on M + 1 (X ) which is called the maximum mean discrepancy (MMD) since [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]:

∀ ν 1 , ν 2 ∈ M + 1 (X ), MMD(ν 1 , ν 2 ) := ∥µ ν1 -µ ν2 ∥ H .
For two given probability measures ν 1 and ν 2 , the MMD must be understood as the distance between their images µ ν1 and µ ν2 in a function space where the two images are well-defined and the distance can be easily computed (or at least approximated from sample data). After reverting to the initial definition of the norm ∥•∥ H in terms of the inner product ⟨•, •⟩ H , an alternative formula of the MMD can be derived:

∀ ν 1 , ν 2 ∈ M + 1 (X ), MMD 2 (ν 1 , ν 2 ) = E ν1⊗ν1 [K(X, X ′ )] + E ν2⊗ν2 [K(X, X ′ )] -2 E ν1⊗ν2 [K(X, X ′ )] . (2.5)
In the first (resp. second) term, X and X ′ are two independent copies of ν 1 (resp. ν 2 ). By way of comparison, in the third term, X and X ′ are still independent but they are, this time, distributed according to ν 1 and ν 2 .

Definition 2.9. Let K : X × X → R be a kernel with RKHS denoted by H. The kernel K is said to be characteristic to

M + 1 (X ) if the map ν ∈ M + 1 (X ) → µ ν ∈ H is injective. Most of
the time, it is quite difficult to determine whether a given kernel is characteristic or not. A necessary condition for a kernel to be characteristic is to generate an infinite-dimensional RKHS. For some specific types of kernels, there also exist sufficient conditions (no longer involving the notion of kernel mean embedding) [START_REF] Sriperumbudur | On the relation between universality, characteristic kernels and RKHS embedding of measures[END_REF][START_REF] Sriperumbudur | Kernel choice and classifiability for RKHS embeddings of probability distributions[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]. A recent overview of existing results is also provided in [START_REF] Simon-Gabriel | Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions[END_REF].

Remark 2.10. Let K 1 : X 1 ×X 1 → R and K 2 : X 2 ×X 2 →
R be two kernels which are respectively characteristic to M + 1 (X 1 ) and M + 1 (X 2 ). Then, the tensor product kernel

K 1 ⊗ K 2 is characteristic to M + 1 (X 1 × X 2 )
. Against all odds, the generalization of this result to more than two kernels is not so simple. Indeed, it is true for translation-invariant kernels [START_REF] Sriperumbudur | On the relation between universality, characteristic kernels and RKHS embedding of measures[END_REF] (see Proposition 8, p. 777) but false in general [START_REF] Szabó | Characteristic and universal tensor product kernels[END_REF] (see Example 2, p. 9).

Kernels and their feature maps

Definition 2.11. For a given kernel K : X × X → R, let us assume that there exists a Hilbert space G and a map ψ : X → G such that:

∀ x, x ′ ∈ X , K(x, x ′ ) = ⟨ψ(x), ψ(x ′ )⟩ G . (2.6)
ψ is called a feature map, G a feature space and any object ψ(x) a feature function (or simply a feature).

Having access to a feature map ψ : X → G provides a different overview of the kernel K. In particular, if ψ is sufficiently explicit, the kernel behavior (in other words, the mathematical operations performed on the provided data when K is evaluated) may be understood at a deeper level. Given two points x and x ′ in X , computing the value K(x, x ′ ) amounts to:

• creating two feature functions ψ(x) and ψ(x ′ ) living in a possibly very sophisticated feature space G,

• comparing these two feature functions with the metric existing in G. When G is infinite-dimensional (for instance when G is a sequence space or a function space), the metric in G is expected to capture a much richer and more subtle information than the initial metric in X . On the one side, from a mathematical viewpoint, computing K(x, x ′ ) amounts to applying preliminary transformations to x and x ′ before comparing their respective images ψ(x) and ψ(x ′ ) in a higher-dimensional feature space. On the other side, from a numerical viewpoint, everything remains tractable because the value K(x, x ′ ) is computed with a simple analytical formula. Thus, the feature space (hidden within the kernel structure) can be leveraged without even knowing the feature map ψ and the inner product ⟨•, •⟩ G . Kernel evaluations are indeed sufficient to virtually create pairs of feature functions and measure their similarity in the feature space. This property, known as the kernel trick, explains in large part the great popularity of kernel methods [START_REF] Hofmann | Kernel methods in machine learning[END_REF][START_REF] Vert | A primer on kernel methods[END_REF].

Mathematically speaking, there always exists a feature map. Indeed, for any given kernel K, the fundamental property stated in Eq. (2.1) always holds and therefore always offers a naive way of satisfying Eq. (2.6): 7 because the feature space is exactly the RKHS induced by K. Unfortunately, this feature map is not very informative in most cases. In fact, it is often difficult to understand intuitively the true mathematical nature of the operations that allow to transform a point x ∈ X into the function K(•, x) = [θ(x)](•) ∈ H ⊂ R X , especially when H is unknown or partly implicit.

∀ x, x ′ ∈ X , K(x, x ′ ) = ⟨θ(x), θ(x ′ )⟩ H with θ(x) := K(•, x) ∈ H . θ : X → H is called the canonical feature map
Thankfully, the canonical feature map is not the only existing feature map. Depending on the kernel, there may exist many other feature maps and the associated feature spaces may have very different natures. Sometimes, the use of basic tools in mathematical analysis is sufficient to extract a feature map from the analytical expression of the kernel under study, as illustrated in the example below.

Example 2.12. The Gaussian kernel K G : R 2 → R (with bandwidth parameter γ > 0) is defined by:

∀ x, x ′ ∈ R, K G (x, x ′ ) = exp - 1 2 x -x ′ γ 2 .
(2.7)

Remember that the Taylor series expansion of exp(•) holds everywhere on R because the associated radius of convergence is infinite. Hence, the initial expression of K G can thus be rewritten as follows:

K G (x, x ′ ) = exp - 1 2 x γ 2 exp - 1 2 x ′ γ 2 exp x x ′ γ 2 with exp x x ′ γ 2 = ∞ k=0 1 k! x γ k x γ k .
This points the way towards a feature map ψ G from R to the Hilbert space ℓ 2 (N) of all square-summable sequences indexed by N:

K G (x, x ′ ) = ⟨ψ G (x), ψ G (x ′ )⟩ ℓ 2 (N) where ψ G (x) = [ ψ k (x) ] k≥0 ∈ ℓ 2 (N) with ψ k (x) := 1 √ k! exp - 1 2 x γ 2 x γ k . (2.8)
As the feature space is now ℓ 2 (N), ψ G can be interpreted as a catalogue of reference transformations applied to the data. Here, the catalogue consists of an infinite number of damped polynomial features.

Beyond mathematical tricks, for certain types of kernels, there also exists a well-adapted framework to address feature extraction. The case of Mercer kernels is covered in the next section.

Focus on Mercer kernels

Mercer expansions and their L 2 -orthogonal feature maps

From now on, the space of all continuous functions on X is denoted by C(X ). Definition 2.13. A kernel K : X × X → R is said to be a Mercer kernel if X is compact and K is continuous.

Generally speaking, the RKHS associated to any bounded and continuous kernel K is composed of bounded and continuous functions [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] (see Lemma 4.28,p. 128), meaning that H ⊂ C(X ).

Definition 2.14. A Mercer kernel

K : X × X → R is said to be universal if H is uniformly dense in C(X ): ∀ f ∈ C(X ), ∀ ϵ > 0, ∃ h ϵ ∈ H such that ∥f -h ϵ ∥ ∞ := sup x∈X |f (x) -h ϵ (x)| < ϵ .
Remark 2.15. Several variants of Definition 2.14 may be found in the literature [START_REF] Simon-Gabriel | Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions[END_REF][START_REF] Sriperumbudur | On the relation between universality, characteristic kernels and RKHS embedding of measures[END_REF]. In particular, a continuous kernel K defined on R is said to be c 0 -universal if H is uniformly dense in the space C 0 (R) of all continuous functions vanishing at infinity. Remark 2.16. For a Mercer kernel, being universal is a sufficient condition to be characteristic to M + 1 (X ) [START_REF] Sriperumbudur | On the relation between universality, characteristic kernels and RKHS embedding of measures[END_REF] (see Theorem 13, p. 778).

Let K : X × X → R be a Mercer kernel and let ν be a probability measure (with support X ). K and ν are the only two elements needed to define the following integral transform:

T K : L 2 (X , ν) -→ L 2 (X , ν) f -→ T K f with [T K f ](x) := E ν [K(x, X) f (X)] = X K(x, ξ) f (ξ) dν(ξ) .
(2.9)

The linear operator T K is called the kernel integral operator. Changing either K or ν gives rise to another operator with possibly very different properties. The integral operators built from Mercer kernels and Borel probability measures have many remarkable properties. Only those which are directly useful for the rest of this work are recalled hereafter.

T K is a Hilbert-Schmidt operator (i.e. the Hilbert-Schmidt norm of T K is finite). In particular, one has:

∥T K ∥ 2 HS := i≥1 j≥1 |⟨T K e i , e j ⟩ L 2 | 2 = i≥1 j≥1 |⟨K, e i ⊗ e j ⟩ L 2 | 2 = ∥K∥ 2 L 2 < ∞ , ( 2.10)
where (e i ) i≥1 denotes any possible ONB of L 2 (X , ν). The first equality is obtained by combining Eq. (2.9) and Fubini's theorem in L 2 (X 2 , ν ⊗2 ). The second one follows from Parseval's identity after noting that the tensorized system (ϕ i ⊗ ϕ j ) i,j≥1 is an ONB of L 2 (X 2 , ν ⊗2 ). In addition, ∥K∥ L 2 is always a finite norm because K is a Mercer kernel (and is therefore bounded on X ).

Remark 2.17. In light of Eq. (2.9) and (2.10), it is enough to assume that K ∈ L 2 (X 2 , ν ⊗2 ) for T K to be both well-defined and Hilbert-Schmidt. However, only the case of Mercer kernels is considered here because it provides a nice framework where all relevant results can be stated without restriction.

In addition to being Hilbert-Schmidt, T K is positive, compact and self-adjoint. These properties come from the fact that T K is built from a Mercer kernel. At this point, the spectral theory of linear operators allows to go a step further. In fact, the spectral theorem for compact self-adjoint operators [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] (see Theorem A.5.13, p. 505) provides an eigendecomposition of T K which is the cornerstone to demonstrate Mercer's theorem [START_REF] König | Eigenvalue Distribution of Compact Operators[END_REF] (see Theorem 3.a.1, p. 145).

Theorem 2.18 (Mercer). Let K : X × X → R be a Mercer kernel and let ν be a probability measure with support X . Let T K : L 2 (X , ν) → L 2 (X , ν) denote the resulting kernel integral operator. Then, there exists an ONB of L 2 (X , ν) denoted by (ϕ i ) i≥1 which is only composed of eigenfunctions of T K . The associated eigenvalues (λ i ) i≥1 are all non-negative and the eigenfunctions corresponding to positive eigenvalues are continuous on X . In addition, K can be decomposed as follows:

∀ x, x ′ ∈ X , K(x, x ′ ) = ∞ i=1 λ i ϕ i (x) ϕ i (x ′ ) where λ 1 ≥ λ 2 ≥ . . . ≥ 0 , (2.11)
and the convergence of the series is absolute and uniform.

Eq. (2.11) is called the Mercer decomposition (or the Mercer expansion, or the Mercer representation) of K with respect to ν. The rank of K is defined as the number of positive eigenvalues (counted with multiplicity) in the eigenspectrum λ(T K ) := (λ i ) i≥1 . If K only has a finite number of positive eigenvalues, it is said to be a degenerate kernel (or a finite-rank kernel).

Remark 2.19. Using the ONB of eigenfunctions (ϕ i ) i≥1 in Eq. (2.10) leads to ∥T K ∥ 2 HS = i≥1 λ 2 i < ∞. For a Mercer kernel, a little more is known about the decay speed of the eigenvalues. In fact, T K is actually a trace-class operator [START_REF] Brislawn | Kernels of trace class operators[END_REF], which implies that Tr(T

K ) = i≥1 λ i < ∞. Hence, one has (λ i ) i≥1 ∈ ℓ 1 (N * ).
Remark 2.20. For a fixed probability measure ν, the Mercer decomposition established in Eq. (2.11) is not uniquely defined. To be more precise, the eigenspectrum and the eigenspaces of T K are invariant, but there may be many different ways to pick the eigenfunctions from the eigenspaces, especially if some eigenspaces have dimension greater than or equal to 2. Theorem 2.18 may be useful for feature extraction because it provides a feature map φ from X to ℓ 2 (N * ):

K(x, x ′ ) = ⟨φ(x), φ(x ′ )⟩ ℓ 2 (N * ) where φ(x) := λ i ϕ i (x) i≥1 .
(2.12)

φ is called the Mercer feature map8 of K (with respect to ν). Let us insist on the fact that φ strongly depends on ν (either for the decay speed of the eigenvalues or for the shapes of the eigenfunctions). If ν is replaced by another probability measure, a completely different collection of Mercer features could emerge from T K . Theorem 2.18 must not be seen as a miracle solution because the features revealed in Eq. (2.12) are only defined implicitly. In fact, the ability to extract a Mercer feature map φ relies on the ability to solve the infinite-dimensional eigenvalue problem defined by:

T K ϕ = λ ϕ with ϕ ∈ L 2 (X , ν) and λ > 0 .
(2.13)

For common characteristic kernels and probability distributions, this eigenvalue problem seldom has a closedform solution. A list of the few rares examples where this happens is given in Appendix B. For a given Mercer kernel K, the knowledge of an explicit Mercer decomposition allows to characterize the related RKHS thanks to the eigenvalues and eigenfunctions of T K . Theorem 2.21. In regard of the Mercer expansion stated in Eq. (2.11), the RKHS induced by K is:

H = h ∈ R X h(•) = ∞ i=1 a i λ i ϕ i (•) with (a i ) i≥1 ∈ ℓ 2 (N * ) , (2.14)
with inner product:

⟨•, •⟩ H : H × H -→ R h 1 (•) = ∞ i=1 a i √ λ i ϕ i (•) , h 2 (•) = ∞ i=1 b i √ λ i ϕ i (•) -→ ∞ i=1 a i b i .
(2.15)

The system (

√ λ i ϕ i ) i≥1 is therefore an ONB of H.
The reader is referred to [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] (see Theorem 4.51, pp. 150-151) for the detailed proof.

Remark 2.22. With Eq. (2.14), it can be seen that the size of the RKHS is directly linked to the decay rate of the eigenvalues (λ i ) i≥1 . The slower the decay rate, the larger the RKHS. Thus, for two kernels having the same basis of eigenfunctions (ϕ i ) i≥1 , the RKHS induced by the kernel having the fastest decay rate is contained in the RKHS induced by the kernel having the slowest decay rate.

Kernel expansions leading to non-orthogonal feature maps

Theorem 2.21 is actually a specific case of a more general result which applies to any kernel (whether Mercer or not) for which a series expansion is known. Indeed, if a kernel can be decomposed as a series of symmetric and separable functions, a feature map can be directly identified and a feature-based characterization of the associated RKHS follows from this decomposition. This result is stated in the theorem below.

Theorem 2.23. Let X ⊆ R be an interval and let K : X × X → R be a kernel. Let us assume that there exist a countable set I and a system (g i ) i∈I of ℓ 2 -linearly independent functions such that:

∀ x, x ′ ∈ X , K(x, x ′ ) = i∈I g i (x) g i (x ′ ) .
(2.16)

Then, the RKHS induced by K is:

H = h ∈ R X h(•) = i∈I a i g i (•) with (a i ) i∈I ∈ ℓ 2 (I) , (2.17)
with inner product:

⟨•, •⟩ H : H × H -→ R h 1 (•) = i∈I a i g i (•) , h 2 (•) = i∈I b i g i (•) -→ i∈I a i b i .
(2.18)

The system (g i ) i∈I is therefore an ONB of H.

The reader is referred to Appendix F.1 for the detailed proof. It is much inspired from what is done in [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] (see the proof of Theorem 4.51, pp. 150-151) but it is more general since no assumption is made on K. Remark 2.24. An important point of Theorem 2.23 lies in the fact that the system (g i ) i∈I has to be ℓ 2 -linearly independent. As pointed out in [START_REF] Saliani | ℓ 2 -linear independence for the system of integer translates of a square integrable function[END_REF], there exist several non-equivalent ways of defining linear independence in an infinite-dimensional Hilbert space. Among them, ℓ 2 -linear independence means that the implication:

i∈I a i g i (•) = 0 =⇒ ∀ i ∈ I, a i = 0 (2.19)
holds for any square-summable sequence (a i ) i∈I ∈ ℓ 2 (I) such that the series i∈I a i g i (•) converges everywhere on X . By way of comparison, the system (g i ) i∈I is said to be ω-independent if Eq. (2. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]) is verified for any sequence leading to a convergent series. For this reason, the ℓ 2 -linear independence of the features (g i ) i∈I is the weakest assumption one can make for Eq. (2.18) to define an inner product in H. In particular, for any function h(•) = i∈I a i g i (•) belonging to H, this assumption ensures that the sequence (a i ) i∈I is unique in ℓ 2 (I). It also guarantees that the series i∈I a i b i in Eq. (2.18) has always finite sum.

For a given kernel, Theorem 2.23 reveals that finding a decomposition in the same form as Eq. (2.16) is enough to extract a feature map, characterize the RKHS and derive an ONB. Mercer's theorem is nothing but a particular case of this general result, where the prior choice of a probability measure allows to define an integral operator (making L 2 -orthogonal features emerge). Hence, the search for a feature map can be done in two different ways:

• One may try to write K as a sum (or series) of symmetric and separable functions using various expansion methods.

• One may try to find a probability measure ν ∈ M + 1 (X ) for which the eigenvalue problem in Eq. (2.13) can be solved analytically. The first strategy seems simpler but requires the right intuition about how to transform the kernel formula. In comparison, the framework laid by Mercer's theorem, although not offering a miracle solution, clearly indicates a class of eigenvalue problems whose solutions disclose features.

The next section is devoted to a particular family of Mercer kernels which is derived from Sobolev spaces.

Sobolev kernels

In this section, the main goal is to provide a quick overview of Sobolev kernels. First, Section 3.1 describes the Hilbertian setting from which Sobolev kernels originate. In particular, it will be shown that they emerge from the Sobolev spaces H r ([0, 1]) after using well-adapted inner products. For r ≥ 2, we will see in Section 3.2 that it not so simple to obtain a kernel with closed-form expression unless considering the unanchored Sobolev spaces. Some general properties of Sobolev kernels will eventually be put forward in Section 3.3.

Throughout this section, when nothing is mentioned, the reference probability measure is the standard uniform distribution. This remark mainly concerns orthogonality properties and Mercer expansions.

Sobolev spaces and their reproducing kernels

Before we start, let us take time to clarify some notations (although they are, admittedly, very standard in mathematical analysis).

• C([0 , 1]) is the space of all continuous functions on [0 , 1].

• C k ([0 , 1]
) is the space of all k-times continuously differentiable functions on [0 , 1].

• C ∞ ([0 , 1]
) is the space of all infinitely differentiable functions on [0 , 1].

• C ∞ 0 ([0 , 1]
) is the space of all infinitely differentiable functions on [0 , 1] such that ϕ(0) = ϕ(1) = 0. • For any k ≥ 1 and for any h ∈ C k ([0 , 1]), h [k] is the k-th (classic) derivative of f on [0 , 1]. The notations h ′ := h [1] and h ′′ := h [2] will also be used in some situations that lend themselves well to this. Now, let us consider h ∈ L 1 ([0 , 1]). Remember that there exists a unique solution g h ∈ L 1 ([0 , 1]) to the following integral equation:

1 0 g h (x) ϕ(x) dx = - 1 0 h(x) ϕ ′ (x) dx with ϕ ∈ C ∞ 0 ([0 , 1]). (3.1)
This solution is called the weak derivative of h and is often denoted by D 1 h. Thus, the use of an appropriate set of smooth test functions allows to define a generalized notion of derivative in L 1 ([0 , 1]). Further technical details on the mathematical foundations of weak derivatives may be found in [START_REF] Adams | Sobolev spaces[END_REF] (see Sections 1.55-1.62, pp. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Briol | Frank-Wolfe Bayesian quadrature: probabilistic integration with theoretical guarantees[END_REF][START_REF] Brislawn | Kernels of trace class operators[END_REF][START_REF] Campbell | Kernel methods: a survey of current techniques[END_REF].

In particular, if the bounded interval [0 , 1] is replaced by a possibly unbounded domain Ω ∈ R n (for n ≥ 1), everything said above remains valid, provided you replace L 1 ([0 , 1]) by the space L 1 loc (Ω) of all locally integrable functions on Ω.

Since

D 1 h ∈ L 1 ([0 , 1]
) by construction, the process can be iterated to define D 2 h, then D 3 h, and finally any weak derivative D k h of order k ≥ 1. Knowing that L 2 ([0 , 1]) ⊂ L 1 ([0 , 1]), the weak differentiation mechanism can therefore be applied infinitely many times in L 2 ([0 , 1]) and all the resulting weak derivatives belong to L 1 ([0 , 1]). However, there is no guarantee that these functions also belong to L 2 ([0 , 1]). For this reason, the subspace of L 2 ([0 , 1]) where the weak derivatives up to order r ≥ 1 remain in L 2 ([0 , 1]) is a very specific subspace of L 2 ([0 , 1]) called the Sobolev space of order r (on [0 , 1] and for the L 2 -norm):

H r ([0 , 1]) := h ∈ R [0 ,1] ∀ 0 ≤ k ≤ r, D k h ∈ L 2 ([0 , 1]) . (3.2)
The integer parameter r directly controls the level of smoothness within H r ([0 , 1]). Indeed, even if not necessarily intuitive, the integrability conditions imposed on the weak derivatives D k h determine how many times h can be differentiated (in the classic meaning of the word). This point is precisely the subject of the Sobolev embedding theorem [START_REF] Adams | Sobolev spaces[END_REF] (see Theorem 4.12, p. 85) which ensures that H r ([0 , 1]) is continuously embedded in C r-1 ([0 , 1]) for any r ≥ 1. This notably shows that the functions contained in the Sobolev space H r ([0 , 1]) are at least continuous. Hence, there is no need to consider the quotient space of H r ([0 , 1]) for the almost-everywhere equality relation. Now, let us see which inner products are suitable to get the most out of Sobolev spaces. Of course, a naive attempt could be to equip H r ([0 , 1]) with the L 2 -inner product but this would unfortunately not produce an RKHS. A more specific inner product must be used instead:

∀ h 1 , h 2 ∈ H r ([0 , 1]), ⟨h 1 , h 2 ⟩ H r := r k=0 1 0 D k h 1 (x) D k h 2 (x) dx . (3.3)
The inner product ⟨•, •⟩ H r truly accounts for the entire specificity of the functions in H r ([0 , 1]). It is therefore considered as the standard inner product. The resulting induced norm ∥•∥ H r involves the L 2 -norms of all the weak derivatives (from h = D 0 h to D r h). If endowed with this norm, H r ([0 , 1]) becomes an RKHS [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] (see Theorem 121, p. 276) with reproducing kernel denoted by K r . This RKHS (resp. its kernel) is better known as the standard Sobolev RKHS (resp. standard Sobolev kernel) of order r. For r = 1, it was demonstrated in [START_REF] Duc-Jacquet | Approximation des fonctionnelles linéaires sur les espaces Hilbertiens autoreproduisants[END_REF] (see Chapter IV) that K 1 has a closed-form expression:

∀ x, x ′ ∈ [0 , 1], K 1 (x, x ′ ) = 2e e 2 -1 cosh [ min(x, x ′ )] cosh [1 -max(x, x ′ )] , (3.4)
where cosh(•) denotes the hyperbolic cosine function. Moreover, it was later shown in [START_REF] Thomas-Agnan | Computing a family of reproducing kernels for statistical applications[END_REF] (see Corollary 2, p. 27) that K 1 admits a Mercer decomposition based on sinusoidal eigenfunctions (more precisely cosine functions with increasing frequencies) and eigenvalues decaying at a polynomial rate of 1/k 2 (see Appendix B.3.3). When r ≥ 2, unlike what happens for r = 1, there is no explicit formula for the reproducing kernel K r of the standard Sobolev RKHS. In the hope of obtaining a kernel with a closed-form expression even for r ≥ 2, a well-known strategy is to replace ⟨

•, •⟩ H r by another inner product defined on H r ([0 , 1]). Let ⟨•, •⟩ H be one of them. The associated induced norm ∥•∥ H is said to be equivalent to ∥•∥ H r if: ∃ 0 < c 1 < c 2 < ∞ such that ∀ h ∈ H r ([0 , 1]), c 1 ∥h∥ H r ≤ ∥h∥ H ≤ c 2 ∥h∥ H r . (3.5)
By Definition 2.2, if an inner product ⟨•, •⟩ H indeed leads to an equivalent norm, H r ([0 , 1]) equipped with ⟨•, •⟩ H remains an RKHS. This new RKHS is denoted by H in order to avoid any confusion with the standard Sobolev RKHS. According to Theorem 2.5, the reproducing kernel K H (associated to H) is unique and different from the standard kernel K r . Since K r and K H come out from the same Sobolev space, they can be considered as equivalent Sobolev kernels. This terminology is strictly delineated in [START_REF] Oates | Minimum discrepancy methods in uncertainty quantification[END_REF].

Definition 3.1. Let K : [0 , 1] × [0 , 1] → R be a kernel with associated RKHS denoted by H. Then, K is said to be a Sobolev kernel (of order r) if ∥•∥ H and ∥•∥ H r are two equivalent norms on H r ([0 , 1]).
The next section focuses on the particular case of the unanchored Sobolev RKHSs. Just as the standard Sobolev RKHSs, they are built from the Sobolev spaces H r ([0 , 1]). The only difference lies in the choice of the inner product used to build the Hilbertian structure.

Unanchored Sobolev spaces

A simple way to define an alternative inner product on H r ([0 , 1]) is to take:

⟨h 1 , h 2 ⟩ H r Sob := r-1 k=0 1 0 D k h 1 (x) dx 1 0 D k h 2 (x) dx + 1 0 D r h 1 (x) D r h 2 (x) dx . (3.6)
The induced norm ∥•∥ H r Sob is equivalent to the standard Sobolev norm ∥•∥ H r . The two inequalities in Eq. (3.5) are indeed satisfied by ∥•∥ H r Sob . The right-hand side is a direct consequence of the Cauchy-Schwarz inequality (when it is applied in ∥h∥ 2

H r

Sob

to majorize the mean values of the derivatives). As regards the left-hand side, it results from a recursive use of the Poincaré-Wirtinger inequality [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] (see Comment 3.A on Chapter 9, p. 312). In particular, this inequality allows to derive successive upper bounds of ∥h∥ 2 H r which are gradually freed from the L 2 -norms of the low-order derivatives. Definition 3.2. The RKHS obtained when H r ([0 , 1]) is endowed with the inner product defined in Eq. (3.6) is called the unanchored Sobolev space (of order r) and denoted by H r Sob . Contrary to what was deplored for the reproducing kernel K r of the standard Sobolev space, the reproducing kernel K r Sob of the unanchored Sobolev space has an easily computable expression at all orders. Theorem 3.3. For any r ≥ 1, the reproducing kernel of the unanchored Sobolev space H r Sob is given by:

∀ x, x ′ ∈ [0 , 1], K r Sob (x, x ′ ) := 1 + k r Sob (x, x ′ ) = r k=0 B k (x) B k (x ′ ) (k!) 2 + (-1) r+1 (2r)! B 2r (|x -x ′ |) , (3.7)
where (B k ) k≥0 denotes the sequence of Bernoulli polynomials.

For the detailed proof, the interested reader is invited to consult [START_REF] Gu | Smoothing Spline ANOVA Models[END_REF] (see Section 2.3.3, p. [START_REF] Dick | Proof techniques in quasi-Monte Carlo theory[END_REF][START_REF] Dick | Approximation of analytic functions in Korobov spaces[END_REF][START_REF] Dick | Lattice rules for nonperiodic smooth integrands[END_REF][START_REF] Dick | Liberating the weights[END_REF]. Unlike the scattered information found in older works [START_REF] Craven | Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF][START_REF] Wahba | Spline Models for Observational Data[END_REF], this proof is self-contained and easy to read.

Remark 3.4. The kernels (K r

Sob ) r≥1 are all ANOVA. This can be checked very easily with the properties of Bernoulli polynomials (recalled in Appendix A.1), especially the zero-mean property (see Appendix A.1.6) and the symmetry properties (see Appendix A.1.4).

Properties of Sobolev kernels

Several general results regarding Sobolev kernels (in the meaning of Definition 3.1) are now highlighted because they are essential to understand the rest of this work. Proposition 3.5. All Sobolev kernels are Mercer kernels.

Proof. Let K be a Sobolev kernel (of order r ≥ 1). The Sobolev embedding theorem yields H r ([0 , 1]) ⊂ C([0 , 1]). Then, Lemma 4.28 in [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] allows to justify that K is continuous. □ Proposition 3.6. All Sobolev kernels are characteristic to M + 1 ([0 , 1]).

Proof. Let K be a Sobolev kernel (of order r ≥ 1). According to Proposition 3.5 and Remark 2.16, it is sufficient to prove that K is universal. In other words, the only point to justify is the uniform denseness of H r ([

0 , 1]) in C([0 , 1]). With Eq. (3.2), it is obvious that H r ([0 , 1]
) contains all the polynomial functions. According to the Stone-Weierstrass theorem, the algebra of polynomials is uniformly dense in C([0 , 1]). Hence, the same is true for H r ([0 , 1]). Therefore, K is universal, and thus characteristic to

M + 1 ([0 , 1]). □ Proposition 3.7. A Sobolev kernel is not always ANOVA. Proof. A counterexample is provided by the kernel K 1 anch (x, x ′ ) = 1+k 1 anch (x, x ′ ) where k 1 anch (x, x ′ ) := min(x, x ′ ). K 1
anch is the reproducing kernel of the RKHS (denoted by [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] (see Example 23,p. 322). This RKHS is sometimes called the Sobolev space (of order 1) anchored at 0 [START_REF] Griebel | Sparse grids for higher dimensional problems[END_REF][START_REF] Hickernell | On tractability of weighted integration over bounded and unbounded regions in R s[END_REF][START_REF] Kuo | On decompositions of multivariate functions[END_REF]. This explains why the subscript "anch" is used here to denote the RKHS and its reproducing kernel. Coming back to the proof, the kernel k 1 anch cannot be orthogonal since it is positive almost everywhere on [0 , 1]

H 1 anch ) obtained when H 1 ([0 , 1]) is equipped with the inner product ⟨h 1 , h 2 ⟩ H 1 anch := h 1 (0) h 2 (0) + ⟨ h ′ 1 , h ′ 2 ⟩ L 2
2 . Therefore, K 1 anch is not an ANOVA kernel. □ Theorem 3.8. If K is a Sobolev kernel (of order r ≥ 1), then the eigenvalues of the kernel integral operator

T K : L 2 ([0 , 1]) → L 2 ([0 , 1]) verify λ k = O(1/k 2r ).
This result is mentioned in many recent papers [START_REF] Bach | On the equivalence between kernel quadrature rules and random feature expansions[END_REF][START_REF] Wei | Early stopping for kernel boosting algorithms: a general analysis with localized complexities[END_REF][START_REF] Zhang | Divide and conquer kernel ridge regression[END_REF] but the original proof dates back to much older works [START_REF] Birman | Piecewise-polynomial approximations of functions of the classes W α p[END_REF][START_REF] Birman | Estimates of singular numbers of integral operators[END_REF] where more general conclusions on Sobolev spaces are enunciated.

Several illustrations of Theorem 3.8 can be found in Appendix B.3. Polynomial eigendecays are one major singularity of Sobolev kernels. For comparison, Gaussian kernels are characterized by exponential eigendecays and thus correspond to much smaller RKHSs. Now that all the necessary mathematical concepts have been clearly introduced, it is time to explain what kernels are used for in sensitivity analysis.

Kernel-based global sensitivity analysis

The main ambition of this section is to provide a brief overview of kernel methods in global sensitivity analysis (GSA) and to explain why the feature maps of Sobolev kernels provide valuable insights into the newly developed HSIC-ANOVA framework. First, Section 4.1 deals with HSIC indices and highlights their ability to characterize independence. Then, Section 4.2 explains why the (Mercer) feature maps of the input and output kernels help identify the (most important) dependence patterns captured by HSIC indices. Finally, Section 4.3 focuses on the HSIC-ANOVA decomposition. The reasons for the use of Sobolev kernels are detailed and the questions raised by this choice are put forward.

Sensitivity measures based on the Hilbert-Schmidt independence criterion

Let us take the usual notations in the field of uncertainty quantification. In this context, an output of interest Y is computed by a numerical simulator g : X 1 × . . . × X d → Y which is given a set of possibly dependent random variables X 1 , . . . , X d gathered in a random vector X. In terms of modeling, the simulator is handled as a deterministic black-box function. Each random input X i takes its values in X i ⊆ R and follows a probability distribution P Xi ∈ M + 1 (X i ). In addition, a continuous kernel K i : X i ×X i → R (with RKHS denoted by H i ) is assigned to X i . Similar mathematical objects and notations are adopted for the output variable Y .

Definition and estimation of HSIC indices

As recalled in [START_REF] Borgonovo | Sensitivity analysis: a review of recent advances[END_REF][START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF], many existing sensitivity measures seek to quantify the discrepancy between the joint input-output distribution P XiY and the product of marginal distributions P Xi ⊗ P Y (representing hypothetical independence between X i and Y ). The key step is thus the choice of a dissimilarity measure on M + 1 (X i × Y). For example, taking the total-variation distance leads to Borgonovo's indices [START_REF] Borgonovo | A new uncertainty importance measure[END_REF]. Despite real advantages, this method involves a delicate preliminary step in which all input-output densities need to be estimated from the available data [START_REF] Derennes | A nonparametric importance sampling estimator for moment independent importance measures[END_REF]. To avoid this, a kernel strategy consists in using the MMD (see Section 2.2.3) related to the tensor product kernel K i ⊗ K Y . The discrepancy between P XiY and P Xi ⊗ P Y is thus measured through the distance of their respective images in the tensor product RKHS H i ⊗ H Y . This kernel-based approach amounts to applying the Hilbert-Schmidt independence criterion (HSIC) to all input-output pairs:

∀ 1 ≤ i ≤ d, HSIC(X i , Y ) := MMD 2 (P XiY , P Xi ⊗ P Y ) = ∥µ P X i Y -µ P X i ⊗ P Y ∥ 2 Hi⊗H Y . (4.1)
The resulting sensitivity measures are called HSIC indices. Although they were initially intended for variable selection in machine learning [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF][START_REF] Gretton | Kernel methods for measuring independence[END_REF], they have been increasingly used over the past few years in GSA [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF][START_REF] Marrel | The ICSCREAM methodology: identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF].

The alternative formula of the MMD given in Eq. (2.5) allows to rewrite each index HSIC(X i , Y ) as a sum of three expectations:

HSIC(X i , Y ) = E [K i (X i , X ′ i ) K Y (Y, Y ′ )] + E [K i (X i , X ′ i )] E [K Y (Y, Y ′ )] -2 E [K i (X i , X ′ i ) K Y (Y, Y ′′ )] , (4.2)
where (X i , Y ), (X ′ i , Y ′ ) and (X ′′ i , Y ′′ ) are three random pairs following the joint input-output distribution P XiY while being independent of each other. Based on Eq. (4.2), the quantity HSIC(X i , Y ) can be expressed as a single expectation (see Appendix C.2.1), and more precisely as the expectation of a symmetric function involving four independent copies of the input-output vector Z := (X i , Y ). This paves the way to an estimator of HSIC(X i , Y ) in the form of a U-statistic or a V-statistic [START_REF] Gretton | A kernel statistical test of independence[END_REF][START_REF] Song | Supervised feature selection via dependence estimation[END_REF]. Further technical details related to the construction of H U i (U-statistic estimator) and H V i (V-statistic estimator) are postponed to Appendix C.2.2. In the following, when a result is valid for both H U i and H V i , no distinction will be made between them and the HSIC estimator will be merely denoted by H i .

From a theoretical perspective, an important point regarding the efficiency of inference is the existence of a central limit theorem (CLT) that governs the asymptotic behavior of H i [START_REF] Gretton | A kernel statistical test of independence[END_REF] (see Theorem 1, p. 4). This notably shows that H i converges to the exact value HSIC(X i , Y ) at the rate of 1/n. A similar convergence rate can also be achieved for most Sobol' index estimators, as highlighted by the CLTs established in [START_REF] Janon | Asymptotic normality and efficiency of two Sobol index estimators[END_REF] (see Sections 3.1 and 4.2, pp. 345-348) and more recently in [START_REF] Gamboa | Global sensitivity analysis: a novel generation of mighty estimators based on rank statistics[END_REF] (see Section 4.2, p. 2351). Although the convergence speed of an HSIC estimator (taken alone) seems comparable to that of any Sobol' index estimator, it is actually much easier to estimate all HSIC indices than all Sobol' indices. Indeed, the entire collection of HSIC estimates can be directly computed from a single Monte Carlo design of size n, whereas the computation of all Sobol' estimates asks for a Pick-Freeze design composed of at least n (d + 1) samples [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]. Hence, HSIC indices are particularly attractive to tackle large-scale problems because the number of input-output evaluations required to reach sufficient accuracy (on the entire collection of sensitivity indices) does not increase with the dimensionality.

Characterization of independence with HSIC indices

From Eq. (4.1), it is clear that one has HSIC(X i , Y ) = 0 in case of independence between X i and Y . The converse is false in general but becomes true if both K i and K Y are sufficiently sophisticated.

Proposition 4.1. Let K i : X i × X i → R and K Y : Y × Y → R be two kernels such that K i is characteristic to M + 1 (X i ) and K Y is characteristic to M + 1 (Y).
Then, one has:

X i ⊥ Y ⇐⇒ HSIC(X i , Y ) = 0 .
Proof. It is straightforward in view of the definitions given so far. As already mentioned in Remark 2.10, the tensor product kernel

K i ⊗ K Y is characteristic to M + 1 (X i × Y). The map ν ∈ M + 1 (X i × Y) → µ ν ∈ H 1 ⊗ H Y is
therefore injective and this ends the proof. □ Proposition 4.1 is fundamental in the prospect of building a test of independence. Indeed, a plethora of test procedures have been developed on the basis of this property [START_REF] Albert | Adaptive test of independence based on HSIC measures[END_REF][START_REF] Chwialkowski | A kernel independence test for random processes[END_REF][START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF][START_REF] El Amri | More powerful HSIC-based independence tests, extension to space-filling designs and functional data[END_REF][START_REF] El Amri | Optimized HSIC-based tests for sensitivity analysis: Application to thermalhydraulic simulation of accidental scenario on nuclear reactor[END_REF]. Generally speaking, a test of independence is a statistical procedure intended to make a choice between:

(H 0 ) : X i ⊥ Y
vs. (H 1 ) : X i and Y are dependent variables.

If K i and K Y are two characteristic kernels, taking the HSIC allows to transform the initial problem of hypothesis test into a simpler one:

(H 0 ) : HSIC(X i , Y ) = 0 vs. (H 1 ) : HSIC(X i , Y ) > 0.
From here, any estimator of HSIC(X i , Y ) with known distribution under (H 0 ) is a suitable test statistic. In particular, the estimators H U i and

H V i (see Appendix C.2.
2) are two typical examples. Their asymptotic distributions under (H 0 ) were first derived in [START_REF] Gretton | A kernel statistical test of independence[END_REF] (see Theorem 2, p. 4) and several refinements have been proposed since then [START_REF] Zhang | Kernel-based conditional independence test and application in causal discovery[END_REF][START_REF] Zhang | Large-scale kernel methods for independence testing[END_REF]. Naturally, the acceptance (resp. rejection) region of (H 0 ) is composed of the smallest (resp. largest) values taken by the test statistic H i . When the sample size is small, asymptotic test procedures are no longer valid and non-parametric test procedures based on random permutations must be used instead [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF][START_REF] Gretton | A kernel statistical test of independence[END_REF].

Connections between feature maps and HSIC indices

Beyond their original mathematical construction based on kernel mean embeddings, HSIC indices can also be rewritten as generalized covariance operators [START_REF] Baker | Joint measures and cross-covariance operators[END_REF][START_REF] Vakhania | Probability Distributions on Banach Spaces[END_REF]. In fact, it is proved in [START_REF] Gretton | Kernel methods for measuring independence[END_REF] that HSIC(X i , Y ) = ∥C XiY ∥ 2 HS where the cross-covariance operator C XiY : H Y → H i is defined by:

∀ h i ∈ H i , ∀ h Y ∈ H Y , ⟨C XiY h Y , h i ⟩ Hi = Cov(h i (X i ), h Y (Y )) .
After replacing ∥•∥ HS by its definition, one has:

∀ 1 ≤ i ≤ d, HSIC(X i , Y ) = k l |Cov (u ik (X i ) , v l (Y )) | 2 with (u ik ) k an ONB of H i , (v l ) l an ONB of H Y . ( 4.3)
With the above reformulation, HSIC(X i , Y ) appears to be an aggregation of covariance terms obtained after scanning all basis directions in the two RKHSs. Quantifying with Cov(•, •) the linear dependence between u ik (X i ) and v l (Y ) amounts to quantifying a part of the non-linear dependence between X i and Y , more precisely the non-linear dependence pattern characterized by the pair of preliminary transformations (u ik , v l ).

Remark 4.2. In Eq. (4.3), the indexation sets associated to k and l are deliberately not specified. This is done to emphasize the fact that indexing will depend on the selected kernels. Indeed, if a kernel as finite (resp. infinite) rank, the induced RKHS has finite (resp. infinite) dimension, and the associated ONB is composed of a finite (resp. infinite) number of basis functions.

Eq. ( 4.3) reveals that the identification of the dependence patterns captured by the HSIC (when built with K i and K Y ) simply requires the knowledge of ONBs of H i and H Y . For a given RKHS, it has already been explained in Section 2.4 that an ONB can be directly obtained from the kernel expression by finding a Mercer expansion (see Theorem 2.21) or any other series expansion based on symmetric and separable functions (see Theorem 2.23). Depending on the kernel, this may be more or less difficult.

Example 4.3. Let us imagine that HSIC(X i , Y ) is computed with two Gaussian kernels (with scale parameters respectively denoted by γ i and γ Y ). This is the most common situation in practice. A series expansion of the Gaussian kernel K G was provided in Example 2.12 to reveal the hidden existence (within K G ) of infinitely many damped polynomial features. Theorem 2.23 can then be used to justify that the features exactly form an ONB of the RKHS (induced by K G ). With this in mind, the generalized covariance defined in Eq. (4.3) becomes:

HSIC(X i , Y ) = ∞ k=0 ∞ l=0 |Cov (u ik (X i ) , v l (Y )) | 2 with              ∀ k ≥ 0, u ik (x i ) := 1 √ k! exp - 1 2 x i γ i 2 x i γ i k , ∀ l ≥ 0, v l (y) := 1 √ l! exp - 1 2 y γ Y 2 y γ Y l .
Now, let (u ik ) k and (v l ) l be the functions obtained after an L 2 -normalization (with respect to P Xi and P Y ) of the basis functions extracted from H i and H Y . With these new notations, Eq. (4.3) becomes:

HSIC(X i , Y ) = k l ∥u ik ∥ 2 L 2 ∥v l ∥ 2 L 2 |Cov (u ik (X i ) , v l (Y )) | 2 , (4.4)
and it is now possible to find a uniform upper bound for all covariance terms:

0 ≤ |Cov (u ik (X i ) , v l (Y )) | 2 ≤ V(u ik (X i )) V(v l (Y )) = 1 -E P X i [u ik (X i )] 2 1 -E P Y [v l (Y )] 2 ≤ 1 .
Hence, each covariance term in Eq. (4.4) lies in [0 , 1] and it is weighted by a coefficient equal to

∥u ik ∥ 2 L 2 ∥v l ∥ 2 L 2 .
This means that the HSIC captures many different dependence patterns but they are not equally weighted and the weighting system depends on the encountered distributions P Xi and P Y .

• For each input variable X i , the L 2 -norms of the basis functions (u ik ) k can be calculated by hand or approximated by numerical integration because P Xi is provided.

• For the output variable Y , since P Y is unknown, the L 2 -norms of the basis functions (v l ) l can only be estimated from the available output samples. It can be easily proved that the two sequences (∥u ik ∥ L 2 ) k and (∥v l ∥ L 2 ) l are square summable and thus vanish at infinity. As a result, only a small number of dependence patterns actually count towards the final value of HSIC(X i , Y ). Furthermore, an enhanced (if not optimal) version of Eq. (4.4) can be achieved if K i and K Y have explicit Mercer decompositions (for P Xi and P Y respectively):

∀ x i , x ′ i ∈ X i , K i (x i , x ′ i ) = k λ k ϕ ik (x i ) ϕ ik (x ′ i ) with λ i1 ≥ λ i2 ≥ . . . > 0 , ∀ y, y ′ ∈ Y, K Y (y, y ′ ) = l µ l ψ l (y) ψ l (y ′ ) with µ 1 ≥ µ 2 ≥ . . . > 0 .
Theorem 2.21 allows to rewrite Eq. (4.3) with the eigenvalues and eigenfunctions of T Ki and T K Y :

HSIC(X i , Y ) = k l λ ik µ l |Cov (ϕ ik (X i ) , ψ l (Y )) | 2 . (4.5)
Since the eigenfunctions are already L 2 -normalized, Eq. (4.5) is already of the same form as Eq. (4.4) and no renormalization effort is necessary. Therefore, each coefficient λ ik µ l must be interpreted as the weight assigned to the dependence pattern characterized by the pair (ϕ ik , ψ l ). Just as the L 2 -norms in Eq. (4.4), the two sequences of eigenvalues (λ ik ) k and (µ l ) l are square summable and vanish at infinity. They even decrease faster since the eigenfunctions (ϕ ik ) k and (ψ l ) l are orthogonal in the L 2 -sense, contrary to the features (u ik ) k and (v l ) l in Eq. (4.4). In short, the Mercer feature map of K i (resp. K Y ) with respect to P Xi (resp. P Y ) provides the most relevant description of how X i (resp. Y ) is transformed when HSIC(X i , Y ) is computed. In particular, Mercer features have two major advantages over any other collection of features:

• They are naturally ranked by increasing order of influence.

• Their L 2 -norms have a faster decay rate, which results in a sparser representation of HSIC(X i , Y ).

An ANOVA framework for HSIC indices

Despite their many advantages (low estimation cost, characterization of independence, almost no limiting assumption on the input and output distributions), HSIC indices also suffer from major shortcomings. Above all, as the sum of all HSIC indices is not equal to 1, they cannot be interpreted as percentages of influence. Moreover, in the absence of a universal bound for HSIC indices, it is difficult to know what is meant by small and large values. To remedy this difficulty, the ANOVA framework (until now reserved for Sobol' indices) has been freshly extended to kernel-based sensitivity measures [START_REF] Barr | A generalized kernel method for global sensitivity analysis[END_REF][START_REF] Barr | Kernel-based global sensitivity analysis obtained from a single data set[END_REF][START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF]. More specifically, an ANOVA-like decomposition has been set up for HSIC indices in [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] and it has immediately aroused much interest in the GSA community [START_REF] Fellmann | Kernel-based sensivity analysis for (excursion) sets[END_REF][START_REF] Novello | Making sense of dependence: efficient black-box explanations using dependence measure[END_REF][START_REF] Sarazin | Independence test based on total-order HSIC-ANOVA indices[END_REF]. However, this breakthrough was obtained at the cost of stronger assumptions (on both the input kernels and the input probability distributions). (A2) There is an ANOVA kernel

K i = 1 + k i (with RKHS H i = R ⊕ F i ) for each marginal distribution P Xi .
A multivariate ANOVA kernel can then be constructed for each subset of input variables:

∀ u := {u 1 , . . . , u s } ⊆ {1, . . . , d}, K u := K u1 ⊗ K u2 ⊗ . . . ⊗ K us .
(A3) The joint input distribution P X and the output distribution P Y are such that:

∀ u ⊆ {1, . . . , d}, E P Xu [K u (X u , X u )] < ∞ and E P Y [K Y (Y, Y )] < ∞ .
An HSIC index can then be computed for each subset of input variables:

∀ u ⊆ {1, . . . , d}, HSIC(X u , Y ) := MMD 2 (P XuY , P Xu ⊗ P Y ) .
The HSIC-ANOVA decomposition is then given by: The proof may be found in [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] (see Appendix A.3, pp. [START_REF] Derennes | A nonparametric importance sampling estimator for moment independent importance measures[END_REF][START_REF] Dhillon | Kernel k-means: spectral clustering and normalized cuts[END_REF][START_REF] Dick | Proof techniques in quasi-Monte Carlo theory[END_REF][START_REF] Dick | Approximation of analytic functions in Korobov spaces[END_REF][START_REF] Dick | Lattice rules for nonperiodic smooth integrands[END_REF]. It notably relies on earlier works dealing with ANOVA decompositions of multivariate functions in well-adapted RKHSs [START_REF] Kuo | On decompositions of multivariate functions[END_REF][START_REF] Wasilkowski | Finite-order weights imply tractability of linear multivariate problems[END_REF]. The assumption (A1) may seem too restrictive but it is in fact similar to what is demanded for Sobol' indices, at least in their most usual setting9 . Theorem 4.4 provides a rigorous decomposition of the quantity HSIC(X, Y ) into the sum of 2 d -1 HSIC terms, one per each subset X u of input variables. If taking u := {i} in Eq. (4.6), one has HSIC i = HSIC(X i , Y ). This suggests to renormalize HSIC indices in the following way:

HSIC(X, Y ) =
∀ 1 ≤ i ≤ d, S HSIC i := HSIC(X i , Y ) HSIC(X, Y ) . (4.7)
The resulting collection of sensitivity indices are called the first-order HSIC-ANOVA indices [START_REF] Veiga | Basics and Trends in Sensitivity Analysis: Theory and Practice in R[END_REF]. Indices of higher order can also be defined in the same spirit. In particular, the total-order HSIC-ANOVA indices are studied in [START_REF] Sarazin | Independence test based on total-order HSIC-ANOVA indices[END_REF]. Unlike the R 2 -HSIC indices (see Appendix C.1.2), introduced much earlier [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] and resulting from a different renormalization technique, HSIC-ANOVA indices can be regarded as percentage-like importance measures, especially because their sum is equal to 1.

Remark 4.5. Let us assume that the input and output kernels verify (A2) and (A3). In addition, if all kernels are characteristic, HSIC indices are able to characterize independence, and the same is true for the first-order HSIC-ANOVA indices S HSIC i which are proportional to them.

From a theoretical viewpoint, the HSIC-ANOVA methodology successfully combines the advantages of the HSIC and ANOVA frameworks. Now comes the question of its implementation in practice. To comply with the assumptions of Theorem 4.4 and Remark 4.5, the output kernel K Y must be characteristic and a characteristic ANOVA kernel K i = 1 + k i must be found for each input probability distribution P Xi ∈ M + 1 (X i ). As already evoked in Section 2.2.1, there is no characteristic ANOVA kernel for most parametric families, except for the standard uniform distribution. A well-known trick is to replace the initial variable X i by U i := F Xi (X i ) where F Xi denotes the cumulative distribution function of P Xi . Of course, the mathematical modeling of the numerical simulator is adapted accordingly and becomes g :

[0 , 1] d → Y.
As the problem can always be reformulated in this way, X 1 , . . . , X d are now assumed to follow standard uniform distributions. Then, among all possible kernel choices, the practice promoted in [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] is to use Sobolev kernels, and more precisely the unanchored Sobolev kernels (K r Sob ) r≥1 . They are several reasons for this choice: • They are ANOVA (see Remark 3.4) and characteristic (see Proposition 3.6).

• They have a simple analytical formula (unlike the orthogonalized kernels popularized in [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for modelbased sensitivity analysis[END_REF][START_REF] Ginsbourger | On ANOVA decompositions of kernels and Gaussian random field paths[END_REF]).

• No parameter tuning is necessary and the level of smoothness in H r

Sob can be adjusted with the integer parameter r ≥ 1. Up to now, HSIC indices have been mainly estimated with kernels drawn from Gaussian process regression (Gaussian kernels, Laplace kernels, Matérn kernels). To the best of our knowledge, the initiative of computing HSIC indices with Sobolev kernels is one major peculiarity of the HSIC-ANOVA framework. For a complete understanding of the dependence patterns captured by the first-order HSIC-ANOVA indices, one must return to the cross-covariance viewpoint. In particular, Eq. (4.3) becomes:

S HSIC i ∝ HSIC(X i , Y ) = k l |Cov (u ik (X i ) , v l (Y )) | 2 with (u ik ) k an ONB of H r Sob , (v l ) l
an ONB of H Y . Sob . This objective is pursued from Section 5 to Section 8. Of all possible feature maps, Mercer's is the most instructive because it provides the most synthetic description of the role played by K r Sob in the computation of S HSIC i . That is why Sections 5 and 6 only focus on this specific feature map. In Section 5, numerical experiments are first conducted and this surprisingly helps respond many theoretical answers.

Since the rest of this work is dedicated to the study of the unanchored Sobolev kernels, a simplified terminology is adopted. Sobolev kernels will only refer to the kernels (K r Sob ) r≥1 even if Definition 3.1 is more general. Likewise, Sobolev spaces will only refer to the unanchored Sobolev spaces (H r Sob ) r≥1 .

Numerical extraction of the Mercer feature maps

This section aims at using a numerical method called kernel feature analysis (KFA) to get a first idea of what the Mercer expansion of Sobolev kernels may be like. Since it is precised in Theorem 3.8 that the decay speed of the eigenvalues is always polynomial for Sobolev kernels, the only unknown here is the mathematical nature of the eigenfunctions. In Section 5.1, KFA is described in the general context of a Mercer kernel. In Section 5.2, it is applied to Sobolev kernels and valuable conclusions are drawn from this numerical study.

Key principles of KFA

Let K : X × X → R be a Mercer kernel, let ν be a probability measure with support X and let T K be the integral operator built from K and ν. By Theorem 2.18, K admits an eigendecomposition based on a sequence of L 2 -normalized eigenfunctions with positive eigenvalues. The general idea of KFA (originally called the Nyström method [START_REF] Spence | On the convergence of the Nyström method for the integral equation eigenvalue problem[END_REF][START_REF] Spence | Error bounds and estimates for eigenvalues of integral equations[END_REF]) is to estimate this eigendecomposition by solving a discretized version of the eigenvalue problem. Unlike the Galerkin method [START_REF] Krasnosel'skii | Approximate Solution of Operator Equations[END_REF][START_REF] Sloan | Iterated Galerkin method for eigenvalue problems[END_REF] that relies on a deterministic discretization of the eigenvalue equation T K ϕ = λ ϕ, KFA is a simulation-based approach. Only the key points of the method are outlined in what follows. For a step-by-step explanation, the reader is referred to [START_REF] Liang | Eigen-analysis of nonlinear PCA with polynomial kernels[END_REF] or to Appendix D.

Let f ∈ L 2 (X , ν) be an eigenfunction of T K associated to a positive eigenvalue λ. For now, no constraint is imposed on the L 2 -norm of f . It will be seen later how to recover a unit-norm function (as demanded by Theorem 2.18). The starting point to estimate both f and λ is to draw a n-sample x sim := (x i ) 1≤i≤n from ν. In each pointwise equality [T K f ](x i ) = λ f (x i ), the left-hand term is then replaced by its Monte Carlo approximation. The resulting system of n equations then includes n + 1 unknowns (namely λ and the n values of f ). The Gram matrix K n (built from x sim ) helps rewrite the system as a matrix equation:

K n v = (n λ) v with K n := [K(x i , x j )] 1≤i,j≤n ∈ R n×n and v := [f (x i )] 1≤i≤n ∈ R n .
(5.1)

The eigenvalue problem related to T K is thus transformed into a matrix eigenvalue problem which can be solved numerically (as long as n is not too large). Let (γ k ) 1≤k≤n denote the eigenvalues of K n (with γ 1 ≥ . . . ≥ γ n ≥ 0) and let (v k ) 1≤k≤n denote the corresponding eigenvectors. In light of Eq. (5.1), the k-th largest eigenvalue of T K can be estimated by λk = γ k /n. In addition, the eigenvector v k admits the latent structure [f k (x i )] 1≤i≤n where f k is an eigenfunction of T K associated to λ k (the k-th largest eigenvalue). For this reason, v k must be understood as a discrete estimate of f k . As standard eigenvalue algorithms are designed to return unit eigenvectors in R n , one has

∥f k ∥ L 2 ≈ 1/ √ n in practice. To bypass this problem, v k is replaced by w k := √ n v k which is now a discrete estimate of the L 2 -normalized eigenfunction ϕ k := f k /∥f k ∥ L 2 .
Then, for extrapolation purposes, it can be useful to construct an estimate of ϕ k that can be evaluated everywhere on the domain X . An interpolation method called the Nyström extension [START_REF] Fowlkes | Spectral grouping using the Nyström method[END_REF] can be used to build ϕ k : X → R from the only knowledge of γ k and v k (see Appendix D).

From a practical viewpoint, plotting the histogram of the eigenvalues ( λk ) 1≤k≤n allows to visualize how fast the eigenvalues decrease. In particular, if the eigendecay is very fast, the Mercer decomposition can be truncated after a few terms and the corresponding eigenfunctions stand for the most influential features. The eigenfunction estimates ( ϕ k ) 1≤k≤n are also worth plotting. They may help understand what the true eigenfunctions look like, or at least identify some of their properties (monotonicity, periodicity, upper and lower bounds). Sometimes the displayed curves have such distinctive shapes that one can easily guess the analytical expression of the eigenfunctions (see Figures 1, 3 and4). The purpose of the next section is to apply KFA to Sobolev kernels in order to extract relevant information about their eigenvalues and eigenfunctions, and from there to deduce information about their Mercer feature maps.

Application to Sobolev kernels

The procedure detailed in Section 5.1 is now applied to Sobolev kernels. This means that ν is the uniform distribution, X = [0 , 1] and K = k r Sob (for some r ≥ 1). The eigenvalue problem under study is then:

T k r Sob ϕ = λ ϕ with ϕ ∈ L 2 ([0 , 1]) and λ > 0 . (S r λ )
The notation (S r λ ) will be extensively used in all what follows, and sometimes declined to specific values of r. In particular, the letter S was chosen to indicate that the ambition is to extract a feature map by means of a spectral approach. For convenience, two additional notations are introduced:

• λ(T k r

Sob ) is the eigenspectrum of the integral operator T k r Sob . • For any fixed eigenvalue λ, E r (λ) is the eigenspace of λ, i.e. the subspace of L 2 ([0 , 1]) composed of all the eigenfunctions of T k r Sob associated to λ.

Application to the Sobolev kernel k 1

Sob

The Gram matrix K n is built from n = 500 sample points and the results are represented in Figure 1. The histogram plot in Figure 1 (a) allows to visually appreciate the decay speed of the estimated eigenvalues. In Figure 1 (b), the eigenvectors are plotted together with the Nyström approximations that they allow to construct. For each pair ( λk , w k ), the eigenvector coefficients are marked by colored dots, while the resulting Nyström extension ϕ k is represented as a solid line (of the same color) which interpolates these dots. To check if the eigendecay is truly polynomial, a simple trick is to apply a logarithmic transformation to the estimated eigenvalues. As can be observed in Figure 2 (a), the resulting estimates are then arranged along a decreasing straight line. The estimated slope β1 ≈ -2.0 means that λ k = O(1/k 2 ). This is fully consistent with Theorem 3.8, which gives significant credit to the KFA method.

Going back to Figure 1 (b), it can be observed that the eigenfunction estimates ( ϕ k ) 1≤k≤6 look like some sinusoidal functions. Beyond the general appearance of the curves, some other particularities can be observed. Firstly, the eigenfunctions correspond to an increasing number of half-periods. Secondly, there is no phase shift and the eigenfunctions take their largest amplitude at the bounds of the interval [0 , 1]. In light of these elements, all the functions t → √ 2 cos(kπt) with k ≥ 1 emerge as plausible eigenfunctions. This conjecture is indeed true, as confirmed by the proposition below.

Theorem 5.1. For any

k ≥ 1, the function c k : t ∈ [0 , 1] → √ 2 cos(kπt) is an L 2 -normalized eigenfunction of T k 1
Sob with associated eigenvalue λ k := 1/(kπ) 2 . The eigenspectrum λ(T k 1 Sob ) consists of the null eigenvalue 0 and the eigenvalues (λ k ) k≥1 . Accordingly, the Mercer expansion of k 1 Sob may be written as follows:

∀x, x ′ ∈ [0 , 1], k 1 Sob (x, x ′ ) = ∞ k=1 1 (kπ) 2 c k (x) c k (x ′ ) .
(5.2)

Proof. For any given x ∈ [0 , 1], [T k 1 Sob c k ](x) can be expressed as an integral over [0 , 1]. Then, basic calculations allow to show that [T k 1 Sob c k ](x) = c k (x)/(kπ) 2 . They are detailed in Appendix G.1 as supplementary material. This confirms that c k is an eigenfunction and this also indicates that λ k := 1/(kπ) 2 is the corresponding eigenvalue. In addition, since k 1 Sob is an orthogonal kernel, T k 1 Sob 1 = 0 with 1 denoting the constant function equal to 1. This proves that 0 ∈ λ(T k 1 Sob ) and that 1 ∈ E 1 (0). There is no other eigenvalue since the orthonormal system (ONS) defined by (c k ) k≥0 := {1 ; (c k ) k≥1 } is already an ONB of L 2 ([0 , 1]) as recalled in [START_REF] Hernández | A First Course on Wavelets[END_REF] (see Theorem 4.1, p. 21). By Theorem 2.18, the Mercer decomposition of k 1 Sob directly follows. □ Remark 5.2. It must be acknowledged that the Mercer expansion of k 1 Sob was already existing in the literature. A specific remark on the eigendecomposition of T k 1 Sob may be found in [START_REF] Dick | Proof techniques in quasi-Monte Carlo theory[END_REF] while a detailed proof is provided in [START_REF] Dick | Lattice rules for nonperiodic smooth integrands[END_REF] (see Lemma 1, p. 9). The proof technique is rather different from ours since the idea is to expand

k 1 Sob in an ONB of L 2 ([0 , 1]
2 ) obtained by tensorization of the ONB (c k ) k≥0 . In the continuation of this work, two additional proof techniques will be proposed. In Section 6, the Mercer decomposition is recovered by solving a boundary value problem. In Section 7, the trick consists in rewriting k 1 Sob only in terms of B 2 before using a Fourier series expansion. Sob is actually very similar to that of the standard Sobolev kernel K 1 (see Appendix B.3). The eigenfunctions {1 ; (c k ) k≥1 } are the same and the eigenvalues are asymptotically equivalent. At first sight, this may be surprising because the analytical expressions of K 1 and K 1 Sob are much different. However, from a theoretical viewpoint, this is pretty natural since the two kernels are related to the same function space H 1 ([0 , 1]) and they result from the use of two equivalent metrics. Furthermore, Theorem 3.8 guarantees that the eigenvalues of T K 1 and T K 1 Sob have the same decay speed, namely O(1/k 2 ).

Now that k 1

Sob has been carefully studied, the question is whether k 2 Sob has similar properties. 5.2.2. Application to the Sobolev kernel k 2

Sob

The estimation method is exactly the same and the results (for n = 500 sample points) are represented in Figure 3. The largest estimated eigenvalue λ1 is predominant to such an extent that it is almost equal to the total eigensum. In comparison with the histogram plot from Figure 1 (a), the eigendecay observed in Figure 3 (a) is much faster. This is not really surprising since Theorem 3.8 indicates that the new decay order is 1/k 4 . However, this situation could not have been fully anticipated only on the basis of Theorem 3.8 because the very start of the eigendecay could as well have been much slower than the final asymptotic trend.

Just as for r = 1, the logarithmic transformation is applied to the eigenvalue estimates. The results are shown in Figure 2 (b) and one can see that the log-transformed eigenvalue estimates perfectly match with a straight line. After being rounded, the estimated slope β1 = -4.44 yields λ k = O(1/k 4 ) and the decay rate stipulated in Theorem 3.8 is recovered.

As the largest eigenvalue λ 1 seems to be heavily dominating in λ(T k 2 Sob ), the Mercer expansion of k 2 Sob may be restricted to its first term. Numerically, one can check that λ1 ≈ 1/12. It can also be observed in Figure 3 (b) that ϕ 1 looks like the L 2 -normalized zero-mean linear function

P 1 : t ∈ [0 , 1] → 2 √ 3 (t -1/2). Hence, k 2
Sob can be approximated by:

k 2 Sob (x, x ′ ) ≈ λ 1 ϕ 1 (x) ϕ 1 (x ′ ) ≈ λ1 P 1 (x) P 1 (x ′ ) = x - 1 2 x ′ - 1 2 = B 1 (x) B 1 (x ′ ) =: k lin (x, x ′ ) . (5.3)
k lin is the dot-product kernel centered at the midpoint (1/2, 1/2) of [0 , 1] 2 . This kernel is not characteristic and it is therefore unsuitable for many tasks. The fact that k 2 Sob ≈ B 1 ⊗ B 1 also indicates that the numerical behavior of k 2 Sob is almost entirely driven by one single term of its definition. The strong similarity between k 2

Sob and k lin = B 1 ⊗ B 1 is an important conclusion whose implications will be further discussed.

The eigenfunction estimates ( ϕ k ) 1≤k≤6 shown in Figure 3 (b) deserve to be paid some attention. In fact, the shape of ϕ k suggests that the underlying eigenfunction ϕ k might be a polynomial function of degree k: ϕ 1 looks like a straight line, ϕ 2 like a parabola, ϕ 3 like a cubic curve and so on. Remember that the shifted Legendre polynomials (see Appendix A.2) are a family of orthogonal polynomials in L 2 ([0 , 1]). More precisely, it is the only family (P k ) k≥0 such that deg(P k ) = k and ⟨P k , P l ⟩ L 2 = δ kl . As a result, one can naturally wonder whether the eigenfunctions of T k 2 Sob are the shifted Legendre polynomials. Surprisingly, the conjecture is false. Therefore, the visual insights brought by KFA are misleading this time. This situation can be illustrated in the specific case of

P 1 = 2 √ 3 B 1 .
In fact, when T k 2 Sob is applied to B 1 ∝ P 1 , a polynomial of degree 5 comes out:

∀x ∈ [0 , 1], [T k 2 Sob B 1 ](x) = 1 12 708 720 x - 1 2 + 1 120 x 5 - 1 48 x 4 + 1 72 x 3 .
(5.4)

The calculation details are provided in Appendix G.2 as supplementary material. Eq. (5.4) shows that T k 2 Sob B 1 and B 1 cannot be proportional. Hence, P 1 cannot be an L 2 -normalized eigenfunction of T k 2 Sob . The same strategy could be used for P 2 but the integral calculations (which are already tedious for P 1 ) would become extremely tiresome. In Section 6, a much faster way to reach a more general conclusion will be presented.

The next objective is to study what happens when r ≥ 3.

Application to the Sobolev kernels k r

Sob with r ≥ 3 Performing KFA for different Sobolev kernels such that r ≥ 3 confirms that the eigenvalues of T k r Sob decrease with a polynomial rate of 1/k 2r . In addition, two important properties are always verified. Firstly, the largest eigenvalue is strongly predominant and it can be approximated by λ1 = 1/12. Secondly, the eigenfunction estimates ( ϕ k ) k≥1 look like the shifted Legendre polynomials (P k ) k≥1 . Based on these two remarks, all what was said in Section 5.2.2 can be repeated identically and the same conclusions can be drawn. As a result, the Sobolev kernels (k r Sob ) r≥2 and the dot-product kernel k lin are expected to show similar behaviors. This closeness is even accentuated because the eigendecay occurs faster (as r increases) and the truncation of the Mercer expansion in Eq. (5.3) is therefore all the more justified.

One major difference between k 2 Sob and k r Sob emerges as r becomes larger. This can be observed in Figure 4 where KFA is used to estimate the eigenfunctions associated to k 5 Sob . The eigenfunction estimates ( ϕ k ) 1≤k≤6 and the shifted Legendre polynomials (P k ) 1≤k≤6 seem to perfectly match. Although this is again a numerical illusion (as will be demonstrated in Section 6), this suggests that the shifted Legendre polynomials may be asymptotic eigenfunctions (in other words, the eigenfunctions of the limit integral operator obtained as r → ∞). Such a thorny question is not answered here but it will be further studied in Section 8.

As a conclusion, KFA was of great help in the study of the Mercer feature maps associated to Sobolev kernels. For r = 1, KFA provided the guidance to explicit the Mercer decomposition of k 1 Sob and thereby understand that the associated feature map only consists of sinusoidal features. For r = 2, KFA highlighted the fact that k r

Sob behaves like the dot-product kernel k lin centered at (1/2, 1/2) but the exact Mercer decomposition of k 2

Sob

could not be disclosed. More generally, for r ≥ 3, KFA showed that the eigenfunctions involved in the Mercer decomposition of k r Sob are very close to the shifted Legendre polynomials. To go beyond these first conclusions, further investigations on Sobolev kernels are conducted in the next section. In particular, the use of a differential approach provides new answers.

Extraction of Mercer feature maps with a differential approach

In this section, a differential approach is proposed to gain additional insights into the Mercer feature maps of Sobolev kernels (especially in the case where r ≥ 2). For this, the key idea is to differentiate the eigenvalue equation coming from Mercer's theorem and to transform the infinite-dimensional eigenvalue problem (S r λ ) into a boundary value problem, i.e. an ordinary differential equation (ODE) subject to a system of boundary conditions. Section 6.1 explains how to switch from one formulation of the problem to the other. Section 6.2 deals with the resolution of the resulting boundary value problem. Depending on the smoothness parameter r, a closed-form solution may exist or not.

Transformation of the eigenvalue problem into a boundary value problem

Before computing derivatives, the important issue of the differentiability of the eigenfunctions involved in (S r λ ) must be addressed. For fixed r ≥ 1, let ϕ ∈ L 2 ([0 , 1]) be an eigenfunction of T k r Sob . With Theorem 2.21, it is clear that ϕ belongs to the function space H r ([0 , 1]). Then, in virtue of the basic results recalled in Section 3.1, and especially the definition of H r ([0 , 1]) provided in Eq. (3.2), it can be said that:

(i) the weak derivative D k ϕ is well defined at any order k ≥ 1,

(ii) D k ϕ remains in L 2 ([0 , 1]) for all 0 ≤ k ≤ r, (iii) ϕ is r -1 times continuously differentiable on [0 , 1].
In short, ϕ is always at least continuous, with the worst case occurring for r = 1. It must be acknowledged that this simple conclusion could have been directly drawn from Theorem 2.18 (as k r Sob is continuous on [0 , 1] 2 whatever r ≥ 1), and independently of the Sobolev embedding theorem, which is (in a sense) too strong a result. The continuity of ϕ is an essential property. It is actually the key argument to demonstrate that ϕ is infinitely differentiable on [0 , 1]. As a consequence, ϕ is much smoother than first expected, and appears to be a very special function within H r ([0 , 1]). This point is clarified in the theorem below. Proposition 6.1. Let ϕ ∈ L 2 ([0 , 1]) be an eigenfunction of T k r Sob with eigenvalue λ > 0. Then, ϕ ∈ C ∞ ([0 , 1]) and it is a solution of the following ODE:

λ ϕ [2r] + (-1) r+1 ϕ = 0 . (E r λ )
The reader is referred to Appendix F.2 for the detailed proof. Knowing that ϕ ∈ C([0 , 1]), the Leibniz integral rule allows to justify that the eigenvalue equation λ ϕ = T k r Sob ϕ can be differentiated on both sides. This leads to an integral expression of λ ϕ [1] and the same technique can be repeated. The derivatives ϕ [1] , . . . , ϕ [2r] are computed recursively and the properties of Bernoulli polynomials make ϕ be a solution of the ODE (E r λ ). As ϕ [2r] ∝ ϕ, a simple proof by induction leads to ϕ ∈ C ∞ ([0 , 1]).

With the ODE (E r λ ) from Proposition 6.1, it can be proved much more easily that T k r Sob does not admit any polynomial eigenfunction, as stated in Corollary 6.2 below.

Corollary 6.2. Let ϕ ∈ L 2 ([0 , 1]) be an eigenfunction of T k r

Sob with eigenvalue λ > 0. Then, ϕ cannot be a polynomial function.

Proof. For the sake of contradiction, let us assume that T k r Sob admits one polynomial eigenfunction ϕ with eigenvalue λ > 0. Let P ∈ R[X] denote the polynomial expression associated to the polynomial function ϕ ∈ R [0,1] . According to Theorem 6.1, ϕ is one solution of the ODE (E r λ ). After switching to polynomial expressions, one has λ P [2r] = (-1) r P . Taking the degree on both sides of the previous equality leads to a contradiction. □

It is important to realize that Proposition 6.1 provides a necessary but not sufficient condition for being an eigenfunction. In particular, for any positive eigenvalue λ, the solution space associated to (E r λ ) is too big compared to what is expected for the solution space of (S r λ ). On the one side, the solution space associated to (S r λ ) is the eigenspace E r (λ). According to what was said in Section 5 (especially in Theorem 5.1), E r (λ) is a 1-dimensional linear subspace of L 2 ([0 , 1]). On the other side, the solution space associated to (E r λ ) is a 2r-dimensional linear subspace of L 2 ([0 , 1]). As a consequence, additional constraints must be added to (E r λ ) in order to restrain the solution space. This is the subject of the theorem below. Theorem 6.3. Let ϕ ∈ L 2 ([0 , 1]) and λ > 0. Then, the two following statements are equivalent:

(i) ϕ is an eigenfunction of the integral operator T k r Sob with eigenvalue λ. (ii) ϕ is a solution of the boundary value problem defined as:

λ ϕ [2r] + (-1) r+1 ϕ = 0 with          ϕ [r] (0) = ϕ [r] (1) = 0 ∀ 0 ≤ p ≤ r -2, (-1) r+p (ϕ [p] (1) -ϕ [p] (0)) = ϕ [2r-p-1] (0) ∀ 0 ≤ p ≤ r -2, ϕ [2r-p-1] (0) = ϕ [2r-p-1] (1) . (B r λ )
The reader is referred to Appendix F.3 for the detailed proof. As Theorem 6.3 is an extended version of Proposition 6.1, the proof does not repeat what is said in Appendix F.2 and rather focuses on how to obtain the boundary conditions and the converse statement. For the boundary conditions, the derivatives of the eigenvalue equation (already calculated for the purposes of Proposition 6.1) need to be evaluated at the bounds of [0 , 1] and the properties of Bernoulli polynomials eventually lead to the expected constraints. For the converse statement, one must write T k r Sob ϕ = (-1) r λ T k r Sob ϕ [2r] and repeat 2r times a trick that consists in using an integration by parts before selecting an appropriate boundary condition to enable simplifications. As many analytical techniques are available to solve linear ODEs, the reformulation of (S r λ ) as (B r λ ) opens the way for a possibly faster and simpler analytical resolution. This question is discussed in the next section.

Analytical resolution of the boundary value problem

For r = 1, the boundary value problem (B r λ ) becomes:

λ ϕ ′′ + ϕ = 0 with ϕ ′ (0) = ϕ ′ (1) = 0 and λ > 0 . (B 1 λ )
The resolution is straightforward (see Appendix E.1) and fortunately allows to recover the Mercer expansion already established in Theorem 5.1.

For r ≥ 2, solving the homogeneous linear ODE (E r λ ) simply boils down to finding the complex roots of the characteristic polynomial χ r (z) := λ z 2r + (-1) r+1 . To achieve this, one can write:

χ r (z) = 0 ⇐⇒ λ z 2r = (-1) r ⇐⇒ 2r √ λ z 2r = i 2r ⇐⇒ ζ 2r = 1 with ζ := -i 2r √ λ z .
After using the 2r-th roots of unity in C, a solution set containing 2r elements can be derived:

χ r (z) = 0 ⇐⇒ ζ ∈ e 2iπ( j 2r ) with 0 ≤ j ≤ 2r -1 ⇐⇒ z ∈ ξ e 2iπ( j 2r + 1 4 ) with ξ := 1 2r √ λ and r 2 ≤ j ≤ 3r 2 -1 .
where ⌊•⌋ and ⌈•⌉ are respectively the floor and ceiling functions. For the sake of convenience, the roots of χ r are indifferently denoted by (z i ) 1≤i≤2r in what follows. Any solution of the ODE (E r λ ) can then be expressed as: Since all coefficients in χ r are real numbers, there also exists a real-valued variant of Eq. ( 6.1) but it would be tedious and useless to formalize it. Thus, the construction of the solution space associated to (E r λ ) is not a problem, especially if a basis of complex-valued solutions is used to do so. On the contrary, it is much harder to restrain this solution space according to the boundary conditions specified in (B r λ ). In fact, when the boundary conditions are rewritten in terms of the general solution given in Eq. (6.1), a non-linear system depending on both ξ and w 1 , . . . , w 2r comes out and it cannot be solved analytically. Remark 6.5. For the Laplace kernel and the uniform distribution on a symmetric interval [-a , a], the same pitfall is encountered (see Appendix B.2). The eigenvalue problem is transformed into a boundary value problem which cannot be solved completely [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF] (see Section 2.3.3, pp. 27-45). In particular, the Mercer expansion remains partly implicit because the eigenvalues rely on the unknown solutions of a transcendental equation.

∀t ∈ [0 , 1], ϕ(t) =
Let us take a close look at what happens for r = 2. This time, the boundary value problem is given by: λ ϕ [4] -ϕ = 0 with ϕ ′′ (0) = ϕ ′′ (1) = 0 ϕ(1) -ϕ(0) = ϕ [3] (0) = ϕ [3] (1) and λ > 0 , (B 2 λ )

and any real-valued solution of the ODE may be written as:

ϕ(t) = αe ξt + βe -ξt + γ cos(ξt) + δ sin(ξt) with ξ := 1 4 √ λ and (α, β, γ, δ) ∈ R 4 . (6.2)
With this kind of general solution, the boundary conditions may be rearranged as a system of 2r = 4 equations:

M 2 (ξ) w = 0 where M 2 (ξ) :=     1 1 -1 0 e ξ e -ξ -cos(ξ) -sin(ξ) e ξ -1 1 -e -ξ sin(ξ) 1 -cos(ξ) e ξ -ξ 3 -1 e -ξ + ξ 3 -1 cos(ξ) -1 sin(ξ) + ξ 3     and w :=     α β γ δ     .
Contrary to what such a rearrangement might suggest, the equation M 2 (ξ) w = 0 is not a linear system because most coefficients in M 2 (ξ) depend on ξ (or equivalently on λ) which is unknown. Therefore, a system of 2r = 4 equations in 2r + 1 = 5 unknowns must be solved. However, as ξ is the only unknown involved in the parametrization of the square matrix M 2 (ξ), the specific structure adopted by the boundary conditions paves the way to a simple characterization of the eigenvalues:

• If M 2 (ξ) is an invertible matrix, the linear equation M 2 (ξ) w = 0 admits w = 0 as unique solution in R 4 . In this case, the solution space of (B 2 λ ) is restricted to ϕ = 0. Hence, λ = 1/ξ 4 is not an eigenvalue. • If M 2 (ξ) is a singular matrix, M 2 (ξ) w = 0 admits infinitely many non-zero solutions in R 4 . The functions resulting from the use of such coefficients in Eq. ( 6.2) are non-zero solutions of the boundary value problem (B 2 λ ). They thus live in the eigenspace E 2 (λ) associated to the eigenvalue λ = 1/ξ 4 . As a consequence, non-trivial eigenfunctions only arise from situations where the matrix M 2 (ξ) is singular. The eigenspectrum of T k 2 Sob can thus be found out after solving the equation:

η 2 (ξ) := det [M 2 (ξ)] = 0 .
Computing analytically the determinant of M 2 (ξ) yields:

∀ξ > 0, η 2 (ξ) = -2e -ξ 2 (e ξ -1) 2 sin(ξ) + (ξ 3 + e 2ξ (ξ 3 -2) + 2) cos(ξ) -2 (1 + ξ 3 e ξ -e 2ξ ) . (6.

3)

The roots of η 2 on (0 , +∞) do not have a closed-form expression. However, η 2 is asymptotically equivalent to the very simple function η ∞ 2 defined below.

Proposition 6.6. The asymptotic behavior of η 2 obeys the following equivalence relation:

η 2 (ξ) ∼ ξ→∞ η ∞ 2 (ξ) := -a 2 (ξ) b 2 (ξ) with    a 2 (ξ) := 2 ξ 3 e ξ b 2 (ξ) := cos(ξ)
.

Proof. After identifying that ξ 3 e 2ξ is the leading-order factor in Eq. ( 6.3), η ∞ 2 follows immediately:

η 2 (ξ) = -2ξ 3 e ξ        2(e ξ -1) 2 ξ 3 e 2ξ ---→ ξ→∞ 0 sin(ξ) + ξ 3 + e 2ξ (ξ 3 -2) + 2 ξ 3 e 2ξ ---→ ξ→∞ 1 cos(ξ) - 2(1 + ξ 3 e ξ -e 2ξ ) ξ 3 e 2ξ ---→ ξ→∞ 0        ∼ ξ→∞ -2 ξ 3 e ξ cos(ξ) .

□

A direct consequence of Proposition 6.6 is the following approximation formula for small eigenvalues:

η 2 (ξ) = 0 for ξ → ∞ ⇐⇒ η ∞ 2 (ξ) = 0 for ξ → ∞ ⇐⇒ cos(ξ) = 0 for ξ → ∞ ⇐⇒ ξ k = 2k+1 2 π for k → ∞ ⇐⇒ λ k = 1 2k+1 2 π 4 for k → ∞ . (6.4)
These asymptotic derivations are consistent with Theorem 3.8 which says that the eigenvalues decrease with a rate proportional to 1/k 4 in the case of K 2 Sob . In fact, Eq. (6.4) does even better than Theorem 3.8 because the eigendecay is henceforth known up to the multiplicative constant (and not only bounded).

For r = 3, the same method can be used to deduce the decay rate of the eigenvalues. The general solution of the ODE (E 3 λ ) is easy to obtain (see Appendix E.2 for further details) but it is still a bit more complicated to leverage the extra information brought by the boundary conditions. They can be reorganized as M 3 (ξ) w = 0 with w containing the 2r = 6 coefficients involved in the general solution of the ODE. Therefore, the eigenvalues of T k 3 Sob correspond to the roots of the function η 3 : ξ → det [M 3 (ξ)] on (0 , +∞). The intractable expression of M 3 (ξ) makes pointless the calculation by hand of the determinant involved in the definition of η 3 . Therefore, no counterpart of Eq. ( 6.3) is available when r = 3. However, an asymptotically equivalent function η ∞ 3 can be found, as shown below.

Proposition 6.7. The asymptotic behavior of η 3 obeys the following equivalence relation:

η 3 (ξ) ∼ ξ→∞ η ∞ 3 (ξ) := -a 3 (ξ) b 3 (ξ) with    a 3 (ξ) := 3 4 ξ 8 e √ 3 ξ b 3 (ξ) := sin(ξ)
.

The reader is referred to Appendix E.2 for the entire proof. Once again, an approximation formula for small eigenvalues can be derived:

η 3 (ξ) = 0 for ξ → ∞ ⇐⇒ η ∞ 3 (ξ) = 0 for ξ → ∞ ⇐⇒ sin(ξ) = 0 for ξ → ∞ ⇐⇒ ξ k = kπ for k → ∞ ⇐⇒ λ k = 1 (kπ) 6 for k → ∞ . ( 6.5)
Just as for r = 2, the approximation formula obtained after asymptotic derivations may be seen as a refined version of Theorem 3.8 which only allowed to affirm that λ k = O(1/k 6 ).

For r ≥ 4, the general solution of the ODE (E r λ ) can be derived in the same way as for r ∈ {2, 3}. However, the analytical resolution of the boundary value problem (B r λ ) cannot be pushed any further because the expression of M 2r (ξ) is too cumbersome. Therefore, no generalization of Eq. (6.4) and (6.5) might be obtained, unless using a symbolic math software.

As a conclusion, Section 6 focused on the transformation of the original eigenvalue problem (S r λ ) into the boundary value problem (B r λ ). A general solution can always be found for the associated ODE but it is then much harder to take into account the boundary conditions. A complete analytical resolution is possible for r = 1 and allows to recover the Mercer expansion of k 1 Sob . For r ≥ 2, the boundary conditions result in a non-linear system of equations without closed-form solution. This unfortunately means that the Mercer expansion of k r Sob remains partly implicit, and that we should not hope for better.

Beyond the expected conclusions, the differential approach also brings additional answers. On the one hand, it is now mathematically proved that the eigenfunctions cannot be of polynomial type. On the second hand, when the boundary conditions do not allow for an exact resolution, they can still be used to derive approximation formulas which confirm that the decay rates are polynomial with order of 1/k 2r (at least for r ∈ {1, 2, 3}). Despite those findings, the identification of a fully explicit feature map remains an open question. In the next section, an approach based on Fourier series provides complementary insights.

Feature maps based on a sub-kernel decomposition

Significant efforts were made in Sections 5 and 6 with the aim of discovering the Mercer decomposition of the Sobolev kernel k r Sob . The related feature map φ r Sob is of particular interest because its features take the form of mutually orthogonal functions in L 2 ([0 , 1]). Unfortunately, φ r Sob has a closed-form expression only for r = 1 and remains partly implicit in all other cases. However, the study of Sobolev kernels must not be stopped here.

The main objective of this section is therefore to identify another feature map ψ r Sob : [0 , 1] → ℓ 2 (N * ). To this end, Section 7.1 seeks to rewrite k r Sob as the sum of two kernels in order to split the initial problem into two simpler problems. Then, it is shown in Section 7.2 that an explicit feature map can be easily extracted for each of the two identified kernels. The feature map ψ r Sob obtained after merging the two collections of features offers a different view on Sobolev kernels and notably allows to derive an ONB of H r Sob (whatever is r ≥ 2). Finally, in Section 7.3, the links between the non-orthogonal feature map ψ r Sob and the Mercer feature map φ r Sob are discussed. In particular, the additional insights provided by ψ r Sob are leveraged to better understand the numerical results observed in Section 5.

Decomposition of Sobolev kernels into a sum of two sub-kernels

Let us consider the following decomposition of k r

Sob as the sum of two kernels:

k r Sob (x, x ′ ) = k r A (x, x ′ ) + k r B (x, x ′ ) with k r A := r i=1 B i (x) B i (x ′ ) (i!) 2 and k r B (x, x ′ ) := (-1) r+1 (2r)! B 2r (|x -x ′ |) .
It is trivial to see that k r A is a kernel as it is a positive linear combination of r symmetric separable functions. In the case of k r B , positive definiteness deserves further explanations. First, it should be noted that all Bernoulli polynomials (except B 1 ) have equal boundary values, i.e. they verify B n (0) = B n [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] for n ̸ = 1. Hence, each polynomial B n may be envisioned as the restriction on [0 , 1] of a 1-periodic continuous function defined on R. This suggests to compute the Fourier series expansion of B n . The Fourier coefficients of B n are actually easy to derive (see Appendix A.1.5) and one eventually has:

∀ n ≥ 2, ∀ x ∈ [0 , 1], B n (x) = (-2) n! ∞ k=1 cos(2kπx -nπ 2 ) (2kπ) n . (7.1)
After replacing B 2r by its Fourier series expansion, k r B turns into:

∀ x, x ′ ∈ [0 , 1], k r B (x, x ′ ) = 2 ∞ k=1 1 (2kπ) 2r [ cos(2kπx) cos(2kπx ′ ) + sin(2kπx) sin(2kπx ′ ) ] . (7.2)
Since k r B is expressed as a series of symmetric separable functions with positive coefficients, it is now obvious that it is a positive definite function on [0 , 1]

2 .

Remark 7.1. K r B = 1 + k r B is a well-known kernel in the literature [START_REF] Bach | On the equivalence between kernel quadrature rules and random feature expansions[END_REF][START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF][START_REF] Craven | Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF]. It is often introduced as the reproducing kernel of the following periodic Sobolev space:

H r per ([0 , 1]) := h ∈ R [0,1] h [k] ∈ L 2 ([0 , 1]) ∀ 0 ≤ k ≤ r h [k] (0) = h [k] (1) ∀ 0 ≤ k ≤ r -1 , (7.3)
when this one is equipped with the inner product:

⟨•, •⟩ H r per : H r per ([0 , 1]) × H r per ([0 , 1]) -→ R ( h 1 , h 2 ) -→ 1 0 h 1 (x) dx 1 0 h 2 (x) dx + 1 0 h [r] 1 (x) h [r] 2 (x) dx . (7.4) H r per ([0 , 1]) is obtained from H r ([0 , 1]
) by adding boundary conditions that force the first r -1 derivatives to have equal boundary values.

Identification of explicit feature maps

Explicit feature map for k r A

For the kernel k r A , a feature map ψ r A : [0 , 1] → R r immediately stands out:

k r A (x, x ′ ) = ⟨ψ r A (x), ψ r A (x ′ )⟩ R r with ψ r A (x) = B k (x) 1≤k≤r where B k (x) := B k (x) k! . (7.5)
The feature space is the Euclidean space R r and every feature function ψ r A (x) consists of r polynomial features. Remark 7.2. As Bernoulli polynomials are not orthogonal in L 2 ([0 , 1]), the definition of k r A is not a Mercer expansion of k r A with respect to the uniform distribution on [0 , 1]. To obtain one, there is no other option than solving analytically the eigenvalue problem related to T k r A . Whether it is straightforward to prove that:

T k r A ϕ = λ ϕ =⇒ ϕ ∈ Span({B 1 , . . . , B r })
, it is much more complicated to go further. In fact, tedious and error-prone hand calculations are required in order to determine which linear combinations of B 1 , . . . , B r correspond to the eigenfunctions of T k r A . Even for r = 3, this asks for considerable efforts. Consequently, there is no point in going further in this direction.

Explicit feature map for k r B

For the kernel k r B , a feature map ψ r B : [0 , 1] → ℓ 2 (N * ) stems from the series expansion established in Eq. (7.2):

k r B (x, x ′ ) = ⟨ψ r B (x), ψ r B (x ′ )⟩ ℓ 2 with ψ r B (x) := (c r 2k ) k≥1 ; (s r 2k ) k≥1 where c r 2k (x) := c 2k (x)/(2kπ) r s r 2k (x) := s 2k (x)/(2kπ) r and c 2k (x) := √ 2 cos(2kπx) s 2k (x) := √ 2 sin(2kπx) . (7.6)
Here, the feature space is the Hilbert space ℓ 2 (N * ) and every feature function ψ r B (x) is composed of an infinite number of sinusoidal features which are organized into pairs. From one pair to the next, the frequency f k := k increases whereas the amplitude A r k := √ 2/(2kπ) r decreases. The ONS {(c 2k ) k≥1 , (s 2k ) k≥1 } is nothing but the Fourier ONB of L 2 ([0 , 1]), after removing the constant function 1 from it. As a consequence, Eq. (7.2) is exactly the Mercer decomposition of k r B and the eigenvalues of T k r B are given by µ k := 1/(2kπ) 2r with k ≥ 1. The polynomial rate of the eigendecay is perfectly in line with Theorem 3.8 because k r B is the reproducing kernel of an infinite-dimensional sub-RKHS of H r ([0 , 1]). In the rest of this work, the notation ψ r B is replaced by φ r B to stress that this feature map is of Mercer's type. It is worth noting that the smoothness parameter r has a very targeted influence on the Mercer decomposition of k r B . One the one hand, the exponent of the power law describing the eigendecay is equal to 2r. On the other hand, r does not interfere in the definition of the eigenfunctions (c 2k ) k≥1 and (s 2k ) k≥1 . This means that all kernels k r B (with r ≥ 1) induce the same initial collection of sinusoidal features but they are not weighted in the same way from one kernel k r B to another. High-frequency features are always more penalized than low-frequency features, and the imbalance between them worsens as r increases.

Explicit feature map for k r

Sob

Replacing k r B by its Mercer expansion in the initial definition of k r Sob yields:

k r Sob (x, x ′ ) = r k=1 B k (x) B k (x ′ ) (k!) 2 + ∞ k=1 c 2k (x) c 2k (x ′ ) (2kπ) 2r + ∞ k=1 s 2k (x) s 2k (x ′ ) (2kπ) 2r (7.7) = r k=1 B k (x) B k (x ′ ) + ∞ k=1 c r 2k (x) c r 2k (x ′ ) + ∞ k=1 s r 2k (x) s r 2k (x ′ ) (7.8)
and a global feature map ψ r Sob : [0 , 1] → ℓ 2 (N * ) follows:

k r Sob (x, x ′ ) = ⟨ψ r Sob (x), ψ r Sob (x ′ )⟩ ℓ 2 with ψ r Sob (x) := [ ψ r A (x) ; ψ r B (x) ] = B k 1≤k≤r ; (c r 2k ) k≥1 ; (s r 2k ) k≥1 .
ψ r Sob is not a Mercer feature map because the identified features are not mutually orthogonal in L 2 ([0 , 1]) even though there exist orthogonality relations between some of them.

In return for its non-orthogonality, ψ r Sob has the great advantage of being totally explicit. In particular, ψ r Sob allows to clearly establish the coexistence of two different kinds of features within k r Sob . Moreover, the role played by the smoothness parameter r in the balance between the polynomial and sinusoidal features is now clear. When r is increased by one unit, the former polynomial features ( B k ) 1≤k≤r-1 remain unchanged, a new polynomial feature B r comes into play and the weights ( √ µ k ) k≥1 assigned to the sinusoidal features (c 2k ) k≥1 and (s 2k ) k≥1 are all multiplied by the factor 1/(2kπ). In short, the amount of polynomial features increases while the amplitude of sinusoidal features is gradually shrunk. This explains why KFA shows polynomial-like eigenfunctions for r = 2 and beyond (see Figures 3 to 4). The numerical results brought by KFA and the new insights brought by the analysis ψ r Sob are indeed consistent and complementary. From a theoretical viewpoint, the analytical expression of ψ r Sob indicates that there is a gradual evolution of the balance between polynomial and sinusoidal features as r becomes larger. From a numerical viewpoint, KFA reveals that this transition is rather abrupt and corresponds to the switch from r = 1 to r ≥ 2. Given ψ 1 Sob , inserting B 2 and replacing the initial weights 1/(2kπ) by smaller weights 1/(2kπ) 2 is thus enough to trigger a definitive reversal in the balance between the polynomial and sinusoidal features. This question is further investigated in Section 7.3.

Feature-based characterization of H r

Sob k r

Sob was rewritten in Eq. (7.8) as a series of symmetric and separable functions. In view of Theorem 2.23, one may monder whether the system defined by:

(g r k ) k := 1 ; ( B k ) 1≤k≤r ; (c r 2k ) k≥1 ; (s r 2k ) k≥1 = { 1 ; (f r k ) k } (7.9)
is an ONB of H r Sob or not. It can be easily proved that this system is not ω-independent (see Remark 2.24). Indeed, after taking n = 2 in Eq. (7.1), a series expansion of B 2 follows:

B 2 (•) = 1 2 B 2 (•) = 2 ∞ k=1 (2kπ) r-2 c r 2k (•) ,
which means that one can find a non-zero sequence of coefficients (γ k ) k making the series k γ k g r k (•) converge and be equal to zero everywhere on [0 , 1]. However, no square-summable sequence is able to do the same.

Proposition 7.3. The system (g r k ) k defined by Eq. (7.9) is ℓ 2 -linearly independent.

The reader is referred to Appendix F.4 for the detailed proof. A key step is to prove that H r ([0 , 1]) is the direct sum of the function spaces induced by the sub-kernels k r A and K r B . Now that Proposition 7.3 is stated, Theorem 2.23 can be applied rigorously. Therefore, it can be concluded that the system (g r k ) k is indeed an ONB of H r Sob .

Comparison with Mercer feature maps

The objective is now to compare the Mercer feature map φ r Sob (known for r = 1 and approximated for r ≥ 2) and the non-orthogonal feature map ψ r Sob (known whatever is r ≥ 1). φ r Sob involves features which are either purely sinusoidal (for r = 1) or pseudo-polynomial (for r ≥ 2). On the contrary, the closed-form expression of ψ r Sob shows that polynomial and sinusoidal features always coexist. To bridge the gap between the two feature maps, the L 2 -norms of the non-orthogonal features (g r k ) k must be computed. In fact, according to what was said in Section 4.2, and more precisely what was established in Eq. (4.4), the contribution of each input basis function g r k to the final value of the index S HSIC i defined in Eq. (4.7) can be read through its norm in L 2 ([0 , 1]). Moreover, it can also be shown that:

∥T K r Sob ∥ 2 HS = k≥1 (λ r k ) 2 ≥ k ∥g r k ∥ 4 L 2 , (7.10)
where (λ r k ) k≥1 are the eigenvalues of T K r Sob . For each feature g r k , the computation of its L 2 -norm allows to better understand its contribution to ∥T K r Sob ∥ 2 HS , and therefore its influence on the Mercer decomposition of K r Sob . For more details on how the above inequality is obtained, the reader is referred to Appendix F.5, where a proof is given in the general case of a Mercer kernel K : X × X → R.

Calculation of L 2 -norms

The properties of Bernoulli polynomials and trigonometric functions allow to derive the L 2 -norms of the features identified in Eq. (7.9). Of course, one immediately has ∥c r 2k ∥ L 2 = ∥s r 2k ∥ L 2 = 1/(2kπ) r . As regards the polynomial features ( B k ) 1≤k≤r , their L 2 -norms are expected to decrease extremely fast because of the presence of factorials in the denominators. In fact, the L 2 -norms of Bernoulli polynomials are closely related to Bernoulli numbers (see Appendix A.1.7) and asymptotic results (see Appendix A.1.2) are available for them:

∥B k ∥ 2 L 2 = |B 2k (0)| (k!) 2 (2k)! with |B 2k (0)| ∼ k→∞ 4 √ kπ k π e 2k .
Then, Stirling's approximation for factorials can be used to deduce an asymptotically-equivalent sequence:

k! ∼ k→∞ √ 2kπ k e k and thus ∥ B k ∥ L 2 = ∥B k ∥ L 2 k! = |B 2k (0)| (2k)! 1/2 ∼ k→∞ √ 2 (2π) k . (7.11)
The decay speed is therefore much lower than first thought. It can be checked numerically that the asymptotic approximation of ∥ B k ∥ L 2 is actually pretty accurate for small values of k (including k = 1). Hence, the less influential polynomial feature (i.e. B r ) and the most influential pair of sinusoidal features (i.e. c r 2 and s r 2 ) are weighted the same:

∥ B r ∥ L 2 ≈ √ 2/(2π) r and ∥c r 2 ∥ L 2 = ∥s r 2 ∥ L 2 = 1/(2π) r . This yields ∥ B r ∥ 2 L 2 ≈ ∥c r 2 ∥ 2 L 2 + ∥s r 2 ∥ 2 L 2
, which strengthens the idea of a competition between B r and (c r 2 , s r 2 ).

Application to the Sobolev kernel k 1

Sob

For r = 1, B 1 is the unique existing polynomial feature and it is facing the entire collection of sinusoidal features {(c 1 2k ) k≥1 ; (s 1 2k ) k≥1 }. In particular, the pair (c 1 2 , s 1 2 ) is directly challenging B 1 in terms of L 2 -norm. Because of this competition, the joint action of polynomial and sinusoidal features generates a specific collection of Mercer features, namely the collection of sinusoidal features identified in Theorem 5.1.

Remark 7.4. If the Mercer expansion of k 1

Sob had not be provided by Theorem 5.1, it could also have been obtained only with the tools from Section 7. The trick consists in rewriting k 1 Sob in the following way:

k 1 Sob (x, x ′ ) = B 1 (x) B 1 (x ′ ) + 1 2 B 2 (|x -x ′ |) = 1 2 x 2 + (x ′ ) 2 -max(x, x ′ ) + 1 3 = B 2 |x -x ′ | 2 + B 2 x + x ′ 2 .
Then, Eq. ( 7.1) applied to B 2 provides:

k 1 Sob (x, x ′ ) = 4 ∞ k=1 cos [kπ(x -x ′ )] (2kπ) 2 + 4 ∞ k=1 cos [kπ(x + x ′ )] (2kπ) 2 = ∞ k=1 c k (x) c k (x ′ ) (kπ) 2 ,
and this offers a third manner of accessing the Mercer decomposition of k 1 Sob .

Application to the Sobolev kernels k r

Sob with r ≥ 2 For r ≥ 2, the competition between B r and (c r 2 , s r 2 ) is still existing but it is overwhelmed by the lower-degree polynomial features B 1 , . . . , B r-1 which have much larger L 2 -norms. Since a small number of polynomial features contribute almost exclusively to ∥T K r Sob ∥ 2 HS , it is quite natural that the Mercer decomposition of k r Sob involves eigenfunctions which have much in common with the shifted Legendre polynomials (see Appendix A.2).

To go even further, it should be emphasized that the L 2 -norm of B 1 is the largest by far. This is consistent with the conclusion stated in Section 5.2 where it was pointed out that B 1 ⊗ B 1 is the leading term in k r Sob (as soon as r ≥ 2).

In this section, the objective was to identify a fully analytical and easily interpretable feature map. This was achieved by rewriting k r Sob as the sum of the two sub-kernels. The unified feature map ψ r Sob arising from this approach is composed of non-orthogonal features (in the L 2 -sense) but has the tremendous advantage of clearly establishing the coexistence of polynomial and sinusoidal features. The balance between both types of features is ruled by the smoothness parameter r. For r = 1, there is a true competition and this gives birth to a purely sinusoidal Mercer feature map. For r ≥ 2, the sinusoidal features are too heavily penalized and they become negligible in comparison with the polynomial features. Knowing this, one may think that the sinusoidal features will vanish asymptotically (i.e. as r → ∞). This intuition will be demonstrated in Section 8.

A glance at the asymptotic feature maps

The objective of this section is to understand how the behavior of Sobolev kernels (K r Sob ) r≥1 evolves when the smoothness parameter r ≥ 1 becomes extremely large. Findings from Section 7 indicate that the influence of sinusoidal features vanishes as r increases. This trend suggests that there will be no sinusoidal feature asymptotically. This intuition has now to be formalized theoretically.

In Section 8.1, it is shown that the sequence of Sobolev kernels converges uniformly to a limit kernel. The analytical expression of this kernel allows to identify an asymptotic feature map and to derive an ONB of the limit RKHS. In Section 8.2, further attention is paid to the asymptotic Mercer feature map. Once again, the shifted Legendre polynomials are proved not to be the unknown eigenfunctions.

Limit Sobolev kernel and associated RKHS

Since Bernoulli polynomials are continuous functions on [0 , 1], they are bounded and attain their bounds:

∀ n ≥ 1, -∞ < m n := min x∈[0 ,1] B n (x) ≤ max x∈[0 ,1] B n (x) =: M n < +∞ . ( 8.1)
Explicit formulas were provided by [START_REF] Lehmer | On the maxima and minima of Bernoulli polynomials[END_REF] in order to compute (or at least approximate) those bounds:

∀ n ≥ 1,        m n = 2 ζ(n) n! (2π) n if n ≡ 0 (mod 4) m n > - 2 n! (2π) n otherwise and            M n = 2 ζ(n) n! (2π) n if n ≡ 2 (mod 4) M n < 2 n! (2π) n otherwise , ( 8.2)
where ζ is the Riemann zeta function:

ζ : N \ {0, 1} -→ R n -→ ζ(n) := ∞ k=1 1 k n . As ζ is decreasing on N \ {0, 1}, one has 1 ≤ ζ(n) ≤ ζ(2) = π 2 /6
≤ 2 and this leads to:

∀ n ≥ 1, ∥B n ∥ ∞ = max x∈[0 ,1] |B n (x)| =: M + n ≤ 2 ζ(n) n! (2π) n ≤ 4 n! (2π) n . ( 8.3)
The above inequality may be seen as a slightly relaxed summary of Eq. (8.1) and (8.2) obtained after taking a loose bound of ζ(n). The term on the right-hand side is not the tightest possible upper bound for M + n but it is quite sufficient to prove the convergence results presented below.

∀ x, x ′ ∈ [0 , 1], k ∞ A (x, x ′ ) := lim r→∞ k r A (x, x ′ ) = ∞ i=1 B i (x) B i (x ′ ) (i!) 2 . (8.4) (b) The kernel sequence (k r B ) r≥1 converges uniformly to the zero kernel k ∞ B := lim r→∞ k r B = 0.
The reader is referred to Appendix F.6 for the detailed proof. This proposition allows to define the limit Sobolev kernel as

K ∞ Sob := 1 + k ∞ Sob = 1 + k ∞ A . The related RKHS is denoted by H ∞ Sob = R ⊕ F ∞ Sob since k ∞

Sob

remains an orthogonal kernel. As Bernoulli polynomials are not mutually orthogonal in L 2 ([0 , 1]), the closed-form expression in Eq. (8.4) is not the Mercer expansion of k ∞ Sob . However, a feature map ψ ∞ Sob : [0 , 1] → ℓ 2 (N * ) can still be identified:

k ∞ Sob (x, x ′ ) = ⟨ψ r Sob (x), ψ r Sob (x ′ )⟩ R r with ψ r Sob (x) = B k (x) k≥1 . ( 8.5)
It is composed of an infinite number of polynomial features having increasing degrees. ψ ∞ Sob must be regarded as the infinite-dimensional generalization of the feature map ψ r A extracted from the finite-rank kernel k r A . In the spirit of what was done in Section 7.2.4 with the polynomial and sinusoidal features extracted from k r Sob , the question here is whether the polynomial features ( B k ) k≥0 are the basis functions of the limit RKHS. Proposition 8.2. The system ( B k ) k≥0 is ℓ 2 -linearly independent.

The reader is referred to Appendix F.7 for the detailed proof. Theorem 2.23 can then be applied seamlessly to K ∞ Sob in order to show that the polynomial features ( B k ) k≥0 form an ONB of H ∞ Sob . Remark 8.3. If a sequence of kernels (K n ) n≥1 converges pointwise, the limit function K ∞ is also a kernel [START_REF] Christmann | Kernels and reproducing kernel Hilbert spaces[END_REF] (see Corollary 4.17, p. 119). However, the characterization of the limit RKHS is in general a delicate issue [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] (see Section 9, pp. 362-368). In the case of K ∞ Sob , the feature-based formalism proposed in Theorem 2.23 offers a very convenient way to describe H ∞ Sob . In particular, the function space can be expressed as:

H ∞ Sob = h ∈ R [0,1] h(•) = ∞ k=0 a k B k (•) k! with (a k ) k≥0 ∈ ℓ 2 (N) .
The feature-based viewpoint can also be used to write ⟨•, •⟩ H ∞ Sob . It would have been much more difficult to obtain the same results with the theory of Sobolev spaces of infinite order [START_REF] Dubinskii | Sobolev Spaces of Infinite Order and Differential Equations[END_REF][START_REF] Dubinskii | Sobolev spaces of infinite order[END_REF].

This first section has allowed to lay the foundations of the asymptotic framework. The limit Sobolev kernel K ∞ Sob has been identified and its properties have been carefully studied. In the next section, the connections between the Mercer feature map of K ∞ Sob and the shifted Legendre polynomials are investigated.

Conclusion

The declared objective of this work was to shed light on how the new generation of HSIC indices (called HSIC-ANOVA indices and based on the use of Sobolev kernels) measures sensitivity. More specifically, as the knowledge of the feature maps is the key to identify the non-linear transformations applied to the input and output variables, the most expected achievement was the extraction of at least one fully explicit feature map from each (unanchored) Sobolev kernel.

Three strategies have been explored: (i) KFA, (ii) the differential approach and (iii) the decomposition into sub-kernels. It was found that their conclusions overlap on many points, while remaining highly complementary in the information they provide. Regarding how the feature maps of Sobolev kernels change with r, the three approaches agree that only two cases need to be distinguished (namely r = 1 vs. r ≥ 2). The Sobolev kernel K 1 Sob is only composed of purely sinusoidal features. The eigenfunctions are cosine functions and the decay rate of the eigenvalues is 1/k 2 . Hence, if K 1 Sob is assigned to all input variables, the information captured by the first-order HSIC-ANOVA indices is now transparent:

S HSIC i ∝ HSIC(X i , Y ) = ∞ k=1 l 1 (kπ) 2 |Cov (c k (X i ) , v l (Y )) | 2 with    ∀ k ≥ 1, c k (x i ) = √ 2 cos(kπx i ) , (v l ) l an ONB of H Y .
In contrast, if r ≥ 2, K r Sob consists of polynomial-like features. The eigenfunctions estimated with the KFA method look like the shifted Legendre polynomials and the eigenvalues decrease much faster. From a numerical standpoint, K 2 Sob and the following Sobolev kernels behave almost like the dot-product kernel K lin because their Mercer feature maps are characterized by the predominance of one single linear feature. If K 2 Sob is used to implement the HSIC-ANOVA decomposition, the information captured by the first-order indices is now:

S HSIC i ∝ HSIC(X i , Y ) ≈ l |Cov (B 1 (X i ) , v l (Y )) | 2 with    B 1 (x i ) = x i -1 2 , (v l ) l an ONB of H Y .
An important lesson is therefore to avoid using a Sobolev kernel K r Sob of order r ≥ 2 within the HSIC. Even though all Sobolev kernels are characteristic, the computation of HSIC-ANOVA indices must be restricted to K 1 Sob . In a more general perspective, this prompts to question the importance of using characteristic kernels in GSA. Of course, characteristic kernels offer theoretical guarantees regarding the ability of the HSIC to characterize independence. However, in practice, this is not sufficient to ensure an efficient detection of statistical dependence. Using a characteristic kernel is a good practice but investigating the macroscopic kernel behavior is just as important, if not more so.

At the end of this work, one can say that most of the questions raised by the use of Sobolev kernels received a satisfactory answer. In light of the numerical and theoretical results, K 1 Sob appears to be a nice kernel that should be trusted to implement the HSIC-ANOVA decomposition. As the first-order HSIC-ANOVA indices are sufficient to characterize independence (see Remark 4.5), one may wonder about the added value of the totalorder indices. This question echoes the numerical study proposed in [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] where the first-order and total-order HSIC-ANOVA indices are almost equal in all presented test cases. The very nature of the additional information captured by higher-order indices will be investigated in future works.

Olivier Roustant 11 , Anthony Nouy 12 and Luc Pronzato 13 for sharing their expertise on the theory of reproducing kernel Hilbert spaces.

INSA Toulouse, Institut Mathématiques de Toulouse, Université Toulouse III -Paul Sabatier, France. École Centrale Nantes, Nantes Université, Laboratoire de Mathématiques Jean Leray, UMR CNRS 6629, France. Laboratoire I3S, Université Côte d'Azur -CNRS, Sophia Antipolis, France. This leads to distinguish three cases:

               B 1 (1) = -B 1 (0) = 1 2 ; ∀ n ≥ 0, B 2n (1) = B 2n (0) ̸ = 0 ; ∀ n ≥ 1, B 2n+1 (1) = B 2n+1 (0) = 0 .
For any n ≥ 2, the polynomial function B n : [0 , 1] → R can be seen as the restriction of a continuous 1-periodic function defined on R.

A.1.5. Fourier series expansions

Any Bernoulli polynomial (except B 1 ) is equal to its Fourier series on [0 , 1]. For B 1 , the equality only holds on (0 , 1) since B 1 cannot be extended as a continuous 1-periodic function. The complex and real versions of the Fourier series are respectively given by:

• ∀ 0 ≤ x ≤ 1, B 0 (x) = 1 ; • ∀ 0 < x < 1, B 1 (x) = - 1 2πi k∈Z * e 2πikx k = - ∞ k=1 sin(2kπx) kπ ; • ∀ n ≥ 2, ∀ 0 ≤ x ≤ 1, B n (x) = - n! (2πi) n k∈Z * e 2πikx k n = (-2) n! ∞ k=1 cos(2kπx -nπ 2 ) (2kπ) n .
A.1.6. Mean values

∀ n ≥ 1, 1 0 B n (x) dx = 0 A.1.7. Integrals • ∀ (i, j) ∈ N 2 , β ij := ⟨B i , B j ⟩ L 2 = 1 0 B i (x) B j (x) dx = (-1) i+1 i! j! (i + j)! B i+j (0) . • ∀ n ∈ N, ∥B n ∥ 2 L 2 = 1 0 B n (x) 2 dx = (n!) 2 (2n)! |B 2n (0)| . A.1.8. Upper bound ∀ n ≥ 1, ∥B n ∥ ∞ = max x∈[0 ,1] |B n (x)| =: M + n ≤ 4 n! (2π) n A.

Legendre polynomials

A few elements about Legendre polynomials are introduced in this section. More details on orthogonal families of polynomials are provided in [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] (see pp.

∀ k ≥ 0, ∥L k ∥ 2 L 2 = 1 -1 L k (t) 2 dt = 2 2k + 1 ̸ = 1 .
For any given interval [a , b], an ONB of L 2 ([a , b]) can be derived from Legendre polynomials. Such basis will be denoted by (P k ) k≥1 in this work. It is simply obtained by rescaling and renormalizing Legendre polynomials:

∀ t ∈ [a , b], σ(t) = 2t -(a + b) b -a and ∀ k ≥ 0, P k (t) := [L k • σ](t) ∥L k • σ∥ L 2 = 2k + 1 b -a L k 2t -(a + b) b -a .
∀t ∈ [0 , 1], ∀k ≥ 0, P k (t) := √ 2k + 1 L k (2t -1).

A.2.2. Explicit expressions for low degrees

The very first Legendre and shifted Legendre polynomials are given below:

• L 0 (x) = 1 ; • P 0 (x) = 1 ; • L 1 (x) = x ; • P 1 (x) = 2 x -1 ; • L 2 (x) = 1 2 3 x 2 -1 ; • P 2 (x) = 6 x 2 -6 x + 1 ; • L 3 (x) = 1 2 5 x 3 -3 x ; • P 3 (x) = 20 x 3 -30 x 2 + 12 x -1 ; • L 4 (x) = 1
∀ x, x ′ ∈ [0 , 1], K W (x, x ′ ) := min(x, x ′ ) and ν is U([0 , 1]) .
Eigenvalues and eigenfunctions:

∀ k ≥ 1, λ k = 1 2k-1 2 π 2 and ϕ k (x) = √ 2 sin 2k -1 2 πx .
Equivalent boundary value problem:

For any λ > 0,

T K W ϕ = λ ϕ ⇐⇒ λ ϕ ′′ + ϕ = 0 with ϕ(0) = 0 ϕ ′ (1) = 0 . B.3.2.
Covariance kernel of the Brownian bridge

∀ x, x ′ ∈ [0 , 1], K B (x, x ′ ) := min(x, x ′ ) -xx ′ and ν is U([0 , 1]) .
Eigenvalues and eigenfunctions:

∀ k ≥ 1, λ k = 1 (kπ) 2 and ϕ k (x) = √
2 sin(kπx) .

Equivalent boundary value problem:

For any λ > 0,

T K B ϕ = λ ϕ ⇐⇒ λ ϕ ′′ + ϕ = 0 with ϕ(0) = 0 ϕ(1) = 0 . B.3.3.
Reproducing kernel of the standard Sobolev space of order r = 1

∀ x, x ′ ∈ [0 , 1], K 1 (x, x ′ ) = 2e e 2 -1 cosh [ min(x, x ′ )] cosh [1 -max(x, x ′ )] and ν is U([0 , 1]) .
Eigenvalues and eigenfunctions:

• λ 0 = 1 and ϕ 0 (x) = 1 ;

• ∀ k ≥ 1, λ k = 1 1 + (kπ) 2 and ϕ k (x) = √ 2 cos(kπx) .

Equivalent boundary value problem:

For any λ > 0,

T K 1 ϕ = λ ϕ ⇐⇒ λ ϕ ′′ + (1 -λ) ϕ = 0 with ϕ ′ (0) = 0 ϕ ′ (1) = 0 . B.3.4.
Reproducing kernel of the unanchored Sobolev space of order r = 1

∀ x, x ′ ∈ [0 , 1], K 1 Sob (x, x ′ ) = 1 + k 1 Sob (x, x ′ ) = 1 + B 1 (x)B 1 (x ′ ) + 1 2 B 2 (|x -x ′ |) and ν is U([0 , 1]) .
Eigenvalues and eigenfunctions:

• λ 0 = 1 and ϕ 0 (x) = 1 ;

• ∀ k ≥ 1, λ k = 1 (kπ) 2 and ϕ k (x) = √
2 cos(kπx) .

Equivalent boundary value problem:

For any λ > 0,

T k 1 Sob ϕ = λ ϕ ⇐⇒ λ ϕ ′′ + ϕ = 0 with ϕ ′ (0) = 0 ϕ ′ (1) = 0 . B.3.5.
Reproducing kernel of the periodic Sobolev space of order r ≥ 1

∀ x, x ′ ∈ [0 , 1], K r per (x, x ′ ) = 1 + k r per (x, x ′ ) = 1 + (-1) r+1 (2r)! B 2r (|x -x ′ |) and ν is U([0 , 1]) .
Eigenvalues and eigenfunctions:

• λ 0 = 1 and ϕ 0 (x) = 1 ;

• ∀ k ≥ 1, λ 1k = 1 (2kπ)

Equivalent boundary value problem:

For any λ > 0, T k r per ϕ = λ ϕ ⇐⇒ λ ϕ [2r] + (-1) r+1 ϕ = 0 with

   ϕ [k] (0) = ϕ [k] (1) ∀ 0 ≤ k ≤ 2r -1 .

B.4. Related literature

If additional details are sought, the reader is kindly referred to the following papers:

• For K G : [START_REF] Williams | Gaussian Processes for Machine Learning, volume 2 of Adaptive Computation and Machine Learning series[END_REF] (see pp. 97-98).

• For K L : [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF] (see pp. 29-32).

• For K W : [START_REF] Hawkins | Some practical problems in implementing a certain sieve estimator of the Gaussian mean function[END_REF] (see pp. 487-492).

• For K B : [START_REF] Corlay | Functional quantization-based stratified sampling methods[END_REF] (see p. 7).

• For K 1 : [START_REF] Duc-Jacquet | Approximation des fonctionnelles linéaires sur les espaces Hilbertiens autoreproduisants[END_REF] (see p. IV.8) and [START_REF] Thomas-Agnan | Computing a family of reproducing kernels for statistical applications[END_REF] (see p. 27).

• For K 1 Sob : [START_REF] Dick | Approximation of analytic functions in Korobov spaces[END_REF] (see pp. 9-10). • For K r per : [START_REF] Wahba | Smoothing noisy data with spline functions[END_REF] (see p. 386).

Appendix C. Inference of kernel-based sensitivity measures

Let g : X 1 × . . . X d → Y be a numerical simulator. The input variables X 1 , . . . , X d and the output variable Y := g(X 1 , . . . , X d ) are assumed to be scalar 14 . They are respectively assigned the kernels K 1 , . . . , K d and K Y . In addition, for each pair Z i := (X i , Y ) composed of one given input and the common output, the two marginal distributions and the joint distribution are respectively denoted by P Xi , P Y and P XiY .

C.1. Kernel-based sensitivity measures

C.1.1. HSIC indices

The most convenient way to write HSIC indices is to use their formula based on expectations [START_REF] Gretton | Kernel methods for measuring independence[END_REF]. In fact, the HSIC index for the pair (X i , Y ) may be expressed as:

HSIC(X i , Y ) = E [K i (X i , X ′ i ) K Y (Y, Y ′ )] + E [K i (X i , X ′ i )] E [K Y (Y, Y ′ )] -2 E [K i (X i , X ′ i ) K Y (Y, Y ′′ )] , (C.1)
where

Z i = (X i , Y ), Z ′ i := (X ′ i , Y ′
) and Z ′′ i := (X ′′ i , Y ′′ ) are three independent random pairs following the same bivariate distribution P XiY .

C.1.2. R 2 -HSIC indices

The normalized HSIC index (also called the R 2 -HSIC index) between X i and Y is defined in the same fashion as Pearson's correlation coefficient (with the HSIC replacing the covariance operator):

R 2 (X i , Y ) := HSIC(X i , Y ) HSIC(X i , X i ) HSIC(Y, Y ) . (C.2)
Normalization has the desired effect since one has 0 ≤ R 2 (X i , Y ) ≤ 1. However, the proof of this result is much more complicated than a simple application of the Cauchy-Schwarz inequality (see [START_REF] Veiga | Basics and Trends in Sensitivity Analysis: Theory and Practice in R[END_REF] for more details).

C.1.3. HSIC-ANOVA indices

Provided that the input and output kernels satisfy all the assumptions of Theorem 4.4, an ANOVA-like decomposition exists for HSIC indices [START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF]. In particular, the HSIC between X := [X 1 , . . . , X d ] and Y can be decomposed as follows:

HSIC(X, Y ) = This allows to define the first-order and total-order HSIC-ANOVA indices in a similar way as Sobol' indices:

S HSIC i := HSIC(X i , Y ) HSIC(X, Y ) and T HSIC i := 1 - HSIC(X -i , Y ) HSIC(X, Y ) (C.3)
where X -i := [X j ] j̸ =i . Moreover, after realizing that HSIC(X, Y ) -HSIC(X -i , Y ) = u∋i HSIC(X u , Y ), it can be easily proved that 0 ≤ S HSIC i ≤ T HSIC i ≤ 1.

C.2. Estimation of the HSIC

Here, the objective is to show how to estimate the HSIC between an input variable X ∼ P X (equipped with the kernel K X ) and an output variable Y ∼ P Y (equipped with the kernel K Y ). Of course, everything can be generalized to a vector of input variables, thus allowing any of the kernel-based sensitivity measure defined in Section C.1 to be estimated.

C.2.1. Preliminary formalization

After grouping all its terms within a single expectation, the HSIC becomes:

HSIC(X, Y ) = E [K X (X, X) K Y (Y, Y ′ )] + E [K X (X, X ′ )] E [K Y (Y, Y ′ )] -2 E [K X (X, X ′ ) K Y (Y, Y ′′ )] = E [K X (X 1 , X 2 ) K Y (Y 1 , Y 2 ) + K X (X 1 , X 2 ) K Y (Y 3 , Y 4 ) -2 K X (X 1 , X 2 ) K Y (Y 1 , Y 3 )] = E [θ(Z 1 , Z 2 , Z 3 , Z 4 )]
where Z 1 := (X 1 , Y 1 ), Z 2 := (X 2 , Y 2 ), Z 3 := (X 3 , Y 3 ) and Z 4 := (X 4 , Y 4 ) are four independent random pairs following the same bivariate distribution P XY , and θ : Z 4 → R (with Z := X × Y) is the 4-argument function defined by:

θ : (z 1 , z 2 , z 3 , z 4 ) -→ K X (x 1 , x 2 ) K Y (y 1 , y 2 ) + K X (x 1 , x 2 ) K Y (y 3 , y 4 ) -2 K X (x 1 , x 2 ) K Y (y 1 , y 3 ) .
Since θ is not symmetric (i.e. it is not invariant under all possible permutations of its arguments), the theory of U-statistics and V-statistics [START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF] (see Chapters 5 and 6, pp. 171-242) cannot be applied directly. To circumvent this pitfall, θ needs to be replaced by the symmetric function:

θ(z 1 , z 2 , z 3 , z 4 ) := 1 4! σ∈S4 θ z σ(1) , z σ(2) , z σ(3) , z σ(4)
where S 4 is the set of all permutations of {1, 2, 3, 4}. Ultimately, one has HSIC(X, Y ) = E θ(Z 1 , Z 2 , Z 3 , Z 4 ) , and this expression is the key to build Monte Carlo estimators.

C.2.2. U-statistic and V-statistic estimators

The available data is a Monte Carlo design of experiments composed of n input-output samples:

Z obs := Z (i) 1≤i≤n = X (i) , Y (i) 1≤i≤n ∼ (P XY ) ⊗n .
The U-statistic estimator of HSIC(X, Y ) is given by:

H U := n 4 -1 1≤i1<...<i4≤n θ Z (i1) , Z (i2) , Z (i3) , Z (i4) = 1 (n) 2 i1̸ =i2 K X X (i1) , X (i2) K Y Y (i1) , Y (i2) . . . - 2 (n) 3 i1̸ =i2̸ =i3 K X X (i1) , X (i2) K Y Y (i1) , Y (i3) . . . + 1 (n) 4 i1̸ =...̸ =i4 K X X (i1) , X (i2) K Y Y (i3) , Y (i4) ,
with (n) p := n!/(n -p)! for any 0 ≤ p ≤ n. By way of comparison, the V-statistic estimator of HSIC(X, Y ) has a similar (although slightly simpler) expression:

H V := 1 n 4 1≤i1,...,i4≤n
θ Z (i1) , Z (i2) , Z (i3) , Z (i4) By default, the Euclidean norm of the eigenvector v k returned by the eigenvalue algorithm is equal to 1. Hence, the L 2 -norm of f k is decreasing with the sample size n:

∥f k ∥ 2 L 2 = X |f k (x)| 2 dν(x) ≈ 1 n n i=1 |f k (x i )| 2 ≈ 1 n n i=1 |v ik | 2 = 1 n ∥v k ∥ 2 R n = 1 n .
This means ∥f k ∥ L 2 ≈ 1/ √ n whereas the eigenfunction ϕ k involved in the Mercer decomposition must respect the unit-norm constraint in L 2 (X , ν). To bypass this problem, v k is replaced by w k := √ n v k and it is then possible to construct a discrete estimate of the L 2 -normalized eigenfunction

ϕ k := f k /∥f k ∥ L 2 : w ik ≈ ϕ k (x i ) .
Of course, one could decide to take -w k instead of w k . A rule has to be laid down to eliminate this last degree of freedom. For this, the simulated data in x sim are sorted in increasing order and the eigenvector w k is renumbered accordingly. There exists an integer m k such that the subsequence (w ik ) m k ≤i≤n is monotonous. The sign of w k is then chosen so that this subsequence is increasing. Thus, among the two possible functions, ϕ k is assumed to be the one which is increasing in the neighborhood of the upper bound. This decision rule applies in all the upcoming numerical experiments.

Remark D.2. The eigenvalue and eigenfunction estimates of T K are simply obtained by properly renormalizing the eigendecompositon of the Gram matrix K n :

K n = n k=1 γ k v k v T k = n k=1 γ k n √ n w k √ n w k T = n k=1 λk ϕ k (x sim ) ϕ k (x sim ) T .

D.3. Nyström extension

For extrapolation purposes, it would be interesting to construct a continuous estimate of the unknown eigenfunction ϕ k . To do this, an interpolation method called the Nyström extension [START_REF] Fowlkes | Spectral grouping using the Nyström method[END_REF] can be used. In this approach, for any given point x ∈ X , the value ϕ k (x) is estimated from the only knowledge of γ k and w k :

ϕ Nys k (x) := 1 γ k n i=1 w ik K(x, x i ) ≈ 1 λk X K(x, ξ) ϕ k (ξ) dν(ξ) ≈ ϕ k (x) , (D.2)
and this provides an overall estimator ϕ Nys k of the k-th eigenfunction ϕ k on X . Computing the value of ϕ Nys k (x) for any given prediction point x ∈ X only asks for n additional kernel evaluations between x and the points in x sim . As ϕ Nys k is expressed as a linear combination of n canonical features, ϕ Nys k inherits from the regularity of the kernel K. In particular, as K is a continuous kernel, ϕ Nys k is in turn a continuous function.

Remark D.3. Although it was not built to this end, note that the function ϕ Nys k interpolates the points x sim and w k . Indeed, one has:

∀ 1 ≤ j ≤ n, ϕ Nys k (x j ) = 1 γ k n i=1 w ik K(x j , x i ) = 1 γ k (K n w k ) j = 1 γ k (γ k w k ) j = w jk .
To go further, the boundary conditions of (B 3 λ ) need to be rewritten in terms of the general solution of the ODE. For the sake of convenience, Eq. (E.1) is rewritten in a more compact style:

ϕ(t) = w T Q(t) = 6 i=1 w i Q i (t | ξ) with w := [ w i ] 1≤i≤6 =         α β γ δ ϵ ζ         and Q(t) := [ Q i (t | ξ) ] 1≤i≤6 =              cos(ξt) sin(ξt) exp ( √ 3 2 ξt) cos ( 1 2 ξt) exp ( √ 3 2 ξt) sin ( 1 2 ξt) exp ( - √ 3 2 ξt) cos ( 1 2 ξt) exp ( - √ 3 2 ξt) sin ( 1 2 ξt)             
.

Assuming that λ = 1/ξ 6 is an eigenvalue of T k 3 Sob , the functions [Q i (• | ξ)]
1≤i≤6 form a basis of the 6-dimensional solution space of the ODE involved in (B 3 λ ). After rewriting the trigonometric functions as the real parts of complex-valued functions, the derivatives of the basis functions can be easily calculated:

∀k ≥ 0,                            Q [k] 1 (t | ξ) = ξ k cos (ξt + kπ 2 ) Q [k] 2 (t | ξ) = ξ k sin (ξt + kπ 2 ) Q [k] 3 (t | ξ) = ξ k exp ( √ 3 2 ξt) cos ( 1 2 ξt + kπ 6 ) Q [k] 4 (t | ξ) = ξ k exp ( √ 3 2 ξt) sin ( 1 2 ξt + kπ 6 ) Q [k] 5 (t | ξ) = ξ k exp ( - √ 3 2 ξt) cos ( 1 2 ξt + 5kπ 6 ) Q [k] 6 (t | ξ) = ξ k exp ( - √ 3 2 ξt) sin ( 1 2 ξt + 5kπ 6 )
.

These formulas allow to rearrange the boundary conditions as a system of equations which can be summarized by M 3 (ξ) w = 0. In order to provide the analytical expression of all coefficients, M 3 (ξ) is described with the help of its two half-matrices:

M 3 (ξ) = M A (ξ) M B (ξ) ∈ R 6×6 with M A (ξ) ∈ R 6×3 M B (ξ) ∈ R 6×3
.

They are respectively given by:

M A (ξ) =                    0 -1 0 sin(ξ) -cos(ξ) -e √ 3 2 ξ sin 1 2 ξ cos(ξ) -1 sin(ξ) e √ 3 2 ξ cos 1 2 ξ + 2π 3 + 1 2 -sin(ξ) cos(ξ) -1 e √ 3 2 ξ cos 1 2 ξ + 5π 6 + √ 3 2 cos(ξ) -1 sin(ξ) + ξ 5 e √ 3 2 ξ cos 1 2 ξ -1 - √ 3 2 ξ 5 -( sin(ξ) + ξ 3 ) cos(ξ) -1 e √ 3 2 ξ cos 1 2 ξ + π 6 - √ 3 2 + 1 2 ξ 3                   
and:

M B (ξ) =                    1 0 1 e √ 3 2 ξ cos 1 2 ξ -e - √ 3 2 ξ sin 1 2 ξ e - √ 3 2 ξ cos 1 2 ξ e √ 3 2 ξ sin 1 2 ξ + 2π 3 - √ 3 2 e - √ 3 2 ξ cos 1 2 ξ -2π 3 + 1 2 e - √ 3 2 ξ sin 1 2 ξ -2π 3 + √ 3 2 e √ 3 2 ξ sin 1 2 ξ + 5π 6 -1 2 e - √ 3 2 ξ cos 1 2 ξ + π 6 - √ 3 2 e - √ 3 2 ξ sin 1 2 ξ + π 6 -1 2 e √ 3 2 ξ sin 1 2 ξ + 1 2 ξ 5 e - √ 3 2 ξ cos 1 2 ξ -1 + √ 3 2 ξ 5 e - √ 3 2 ξ sin 1 2 ξ + 1 2 ξ 5 e √ 3 2 ξ sin 1 2 ξ + π 6 -1 2 - √ 3 2 ξ 3 e - √ 3 2 ξ cos 1 2 ξ + 5π 6 + √ 3 2 + 1 2 ξ 3 e - √ 3 2 ξ sin 1 2 ξ + 5π 6 -1 2 + √ 3 2 ξ 3                    . E.2.2. Proof of Proposition 6.7
It would be useless to seek to obtain a one-line expression of η 3 (ξ) = det [M 3 (ξ)] with a computation by hand of the determinant. To find the function η ∞ 3 which is asymptotically equivalent to η 3 , the trick is to factorize each row L i (with 1 ≤ i ≤ 6) and each column C j (with 1 ≤ j ≤ 6) in M 3 (ξ) by its leading-order term. In particular:

• L 5 must be factorized by ξ 5 .

• L 6 must be factorized by ξ 3 .

• C 3 and C 4 must be factorized by e The idea is then to find a simple function which is asymptotically equivalent to det [N 3 (ξ)]. A first step consists in taking asymptotically equivalent functions for all the coefficients of the matrix N 3 (ξ):

det [N 3 (ξ)] ∼ ξ→∞ 0 -1 0 e - √ 3 2 ξ 0 1 sin(ξ) -cos(ξ) -sin 1 2 ξ cos 1 2 ξ e - √ 3 2 ξ sin 1 2 ξ e - √ 3 2 ξ cos 1 2 ξ cos(ξ) -1 sin(ξ) cos 1 2 ξ + 2π 3 sin 1 2 ξ + 2π 3 1 2 √ 3 2 -sin(ξ) cos(ξ) -1 cos 1 2 ξ + 5π 6 sin 1 2 ξ + 5π 6 - √ 3 2 -1 2 cos(ξ) -1 ξ 5 1 cos 1 2 ξ ξ 5 sin 1 2 ξ ξ 5 √ 3 2 1 2 -1 cos(ξ) -1 ξ 3 cos 1 2 ξ + π 6 ξ 3 sin 1 2 ξ + π 6 ξ 3 1 2 √ 3 2
coefficients have simple expressions and most of them vanish asymptotically. To go further, the determinant calculation rules (and especially the cofactor expansion) must be applied cleverly. More specifically, at each step, the calculations must be carried out so that the coefficients which are asymptotically vanishing enable simplifications. It is a fastidious job but it must be completed patiently in order to obtain:

det [N ∞ 3 (ξ)] ∼ ξ→∞ 3 4 sin(ξ) .
The details of the step-by-step dimension reduction are not provided here for conciseness. This finally leads to the expected result:

η 3 (ξ) = det [M 3 (ξ)] = ξ 8 e √ 3ξ det [N 3 (ξ)] ∼ ξ→∞ ξ 8 e √ 3ξ det [N ∞ 3 (ξ)] ∼ ξ→∞ η ∞ 3 (ξ) = 3 4 sin(ξ) ξ 8 e √ 3ξ .
Theorem F.1 (Leibniz integral rule). Let I ⊆ R be an interval. Let a : I → R and b : I → R be two continuously differentiable functions. Let h : I × R -→ R such that:

• h is continuous on I × R, • h(•, ξ) is continuously differentiable on I for all ξ ∈ R, • h(x, •) is integrable on R for all x ∈ I.
Then, the function defined by:

H : I -→ R x -→ b(x) a(x) h(x, ξ) dξ ,
is continuously differentiable on I with:

∀ x ∈ I, H ′ (x) = b(x) a(x) ∂h ∂x (x, ξ) dξ + b ′ (x) h(x, b(x)) -a ′ (x) h(x, a(x)) . (F.1)
The main idea of our proof is to use the Leibniz integral rule to justify that the eigenvalue equation

[T k r Sob ϕ](•) = λ ϕ(•) is (at least) 2r times differentiable on [0 , 1]. The derivatives λ ϕ [k] = (T k r Sob ϕ) [k]
for all k ∈ {1, . . . , 2r} can then be computed with the help of Eq. (F.1). To begin, let us revert to the integral expression of T k r Sob ϕ:

λ ϕ(x) = T k r Sob ϕ (x) = 1 0 k r Sob (x, ξ) ϕ(ξ) dξ = 1 0 r k=1 B k (x)B k (ξ) (k!) 2 + (-1) r+1 (2r)! B 2r (|x -ξ|) ϕ(ξ) dξ = r k=1 1 0 B k (ξ) ϕ(ξ) dξ B k (x) (k!) 2 + (-1) r+1 (2r)! 1 0 B 2r (|x -ξ|) ϕ(ξ) dξ = r k=1 β k (k!) 2 B k (x) + (-1) r+1 (2r)! x 0 B 2r (x -ξ) ϕ(ξ) dξ + 1 x B 2r (ξ -x) ϕ(ξ) dξ = G A (x) + (-1) r+1 (2r)! G B (x) . (F.2) Of course, G A ∈ C ∞ ([0 , 1]
). Because of the recursive formula ruling the differentiation of Bernoulli polynomials (see Appendix A.1.3), the derivatives of G A can be expressed as follows:

     ∀ p ∈ {1, . . . , r}, G [p] A (x) = r k=p β k k!(k -p)! B k-p (x) ; ∀ p > r, G [p] A (x) = 0 . (F.3)
Regarding G B , it is a bit more complicated. It must be noted that Theorem F.1 applies to both integrals. In particular, ϕ is continuous on [0 , 1] because λ > 0 (according to Theorem 2.18). As a consequence, the functions ξ → B 2r (x -ξ) ϕ(ξ) and ξ → B 2r (ξ -x) ϕ(ξ) verify all the assumptions of Theorem F.1. Therefore, G B is continuously differentiable on [0 , 1] and one has:

G ′ B (x) = 2r x 0 B 2r-1 (x -ξ) ϕ(ξ) dξ + B 2r (0) ϕ(x) -2r x 0 B 2r-1 (ξ -x) ϕ(ξ) dξ -B 2r (0) ϕ(x) = 2r x 0 B 2r-1 (x -ξ) ϕ(ξ) dξ - 1 x B 2r-1 (x -ξ) ϕ(ξ) dξ .
This suggests that G ′ B is also continuously differentiable and G ′′ B can be obtained in the same way. In fact, a simple proof by induction (not given here for conciseness) allows to show that G B is (at least) 2r -1 times continuously differentiable on [0 , 1] with:

∀ 1 ≤ p ≤ 2r -1, G [p] B (x) = (2r)! (2r -p)! x 0 B 2r-p (x -ξ) ϕ(ξ) dξ + (-1) p 1 x B 2r-p (ξ -x) ϕ(ξ) dξ . (F.4)
The key arguments are the Leibniz integral rule, the derivatives of Bernoulli polynomials (see Appendix A.1.3) and the boundary values of Bernoulli polynomials (see Appendix A.1.4). For p = 2r -1, Eq. (F.4) becomes:

G [2r-1] B (x) = (2r)! x 0 B 1 (x -ξ) ϕ(ξ) dξ - 1 x B 1 (ξ -x) ϕ(ξ) dξ . Hence, G [2r-1] B is continuously differentiable on [0 , 1] and G [2r]
B can be computed with Eq. (F.1):

G [2r] B (x) = (2r)! x 0 ϕ(ξ) dξ + B 1 (0) ϕ(ξ) - 1 x ϕ(ξ) dξ -B 1 (0) ϕ(x) = (2r)! 1 0 ϕ(ξ) dξ -ϕ(x) = -(2r)! ϕ(ξ) since B 1 (0) = - 1 2
and ϕ has zero mean.

Finally, the eigenvalue equation [T k r

Sob ϕ](•) = λ ϕ(•) can be differentiated 2r times and one has:

λ ϕ [2r] (x) = (T k r Sob ϕ) [2r] (x) = G [2r] A (x) + (-1) r+1 (2r)! G [2r] B (x) = (-1) r ϕ(x) ,
which justifies that ϕ is indeed a solution of (E r λ ). Besides, the equality λ ϕ [2r] = (-1) r ϕ allows to prove that ϕ is infinitely smooth since:

ϕ ∈ C 2r ([0 , 1]) =⇒ ϕ [2r] ∈ C 2r ([0 , 1]) =⇒ ϕ ∈ C 4r ([0 , 1]) =⇒ ϕ [2r] ∈ C 4r ([0 , 1]) =⇒ . . . =⇒ ϕ ∈ C ∞ ([0 , 1]) .
Remark F.2. Eq. (F.3) and (F.4) will be very useful in the proof of Theorem 6.3 since they provide an integral expression of the first 2r -1 derivatives of ϕ. In short, one may write:

∀ p ∈ {0, . . . , 2r -1}, λ ϕ [p] (x) = G [p] A + (-1) r+1 (2r)! G [p] B = r k=p β k k!(k -p)! B k-p (x) . . . + (-1) r+1 (2r -p)! x 0 B 2r-p (x -ξ) ϕ(ξ) dξ . . . + (-1) p 1 x B 2r-p (ξ -x) ϕ(ξ) dξ . (F.5)
Note that the first term is equal to zero as soon as p > r.

F.3. Proof of Theorem 6.3

Let ϕ ∈ L 2 ([0 , 1]) and λ > 0. It must be proved that the two following statements are equivalent: (i) ϕ is an eigenfunction of the integral operator T k r Sob with eigenvalue λ. (ii) ϕ is a solution of the boundary value problem defined by:

λ ϕ [2r] + (-1) r+1 ϕ = 0 with          ϕ [r] (0) = ϕ [r] (1) = 0 ∀ 0 ≤ p ≤ r -2, (-1) r+p ϕ [p] (1) -ϕ [p] (0) = ϕ [2r-p-1] (0) ∀ 0 ≤ p ≤ r -2, ϕ [2r-p-1] (0) = ϕ [2r-p-1] (1) . (B r λ )
For convenience, specific labels are used to distinguish the three different types of boundary conditions: (C 1 ) ϕ [r] = ϕ [r] (1) = 0 ;

(C 2 ) ∀ 0 ≤ p ≤ r -2, (-1) r+p (ϕ [p] (1) -ϕ [p] (0)) = ϕ [2r-p-1] (0) ; (C 3 ) ∀ 0 ≤ p ≤ r -2, ϕ [2r-p-1] (0) = ϕ [2r-p-1] (1) .
First, let us prove that (i) =⇒ (ii) .

As ϕ is an eigenfunction of T k r Sob with λ > 0, Theorem 6.1 states that it is a solution of the ODE (E r λ ). Then, the integral formula given in Eq. (F.5) to express the first 2r -1 derivatives of T k r Sob ϕ can be used to demonstrate that ϕ verifies all the boundary conditions.

Boundary conditions (C 1 )

• λ ϕ [r]

(0) = β r r! + (-1) 2r+1 r! 1 0 B r (ξ) ϕ(ξ) dξ = β r r! - β r r! = 0 . • λ ϕ [r] (1) = β r r! + (-1) r+1 r! 1 0 B r (1 -ξ) ϕ(ξ) dξ = β r r! + (-1) r+1 r! 1 0 (-1) r B r (ξ) ϕ(ξ) dξ = β r r! - β r r! = 0 .
As a result, λ ϕ [r] (0) = λ ϕ [r] (1) = 0 and this leads to (C 1 ) since λ > 0.

Boundary conditions (C 3 )

For any q ∈ {r + 1, . . . , 2r -1}, one has:

• λ ϕ [q] (0) = (-1) r+1 (2r -q)! (-1) q 1 0 B 2r-q (ξ) ϕ(ξ) dξ = (-1) r+q+1 (2r -q)! β 2r-q .

• λ ϕ [q] (1) = (-1) r+1 (2r -q)! 1 0

B 2r-q (1-ξ) ϕ(ξ) dξ = (-1) r+1 (2r -q)! 1 0 (-1) 2r-q B 2r-q (ξ) ϕ(ξ) dξ = (-1) r+q+1 (2r -p)! β 2r-q .
When p ∈ {0, . . . , r -2}, note that q = 2r -p -1 ∈ {r + 1, . . . , 2r -1}. After operating the associated index shift, one has λ ϕ [2r-p-1] (0) = λ ϕ [2r-p-1] (1) for any p ∈ {0, . . . , r -2} and this leads to (C 3 ) since λ > 0.

Boundary conditions (C 2 )

For any p ∈ {0, . . . , r -2}, one has:

• λ ϕ [p] (0) = r k=p β k k!(k -p)! B k-p (0) + (-1) r+p+1 (2r -p)! β 2r-p . • λ ϕ [p] (1) = r k=p β k k!(k -p)! B k-p (1) + (-1) r+p+1 (2r -p)! β 2r-p .
Because of the boundary values of Bernoulli polynomials (see Appendix A.1.4), the difference of the two terms above can be expressed very simply:

λ (ϕ [p] (1) -ϕ [p] (0)) = r k=p β k k!(k -p)! (B k-p (1) -B k-p (0)) = β p+1 (p + 1)! .
To compute ϕ [2p-r-1] (0), the formula obtained at the previous step can be applied since 2r -p -1 ∈ {r + 1, . . . , 2r -1}:

ϕ [2p-r-1] (0) = (-1) 3r-p (p + 1)! β p+1 = (-1) r+p β p+1 (p + 1)! .
As a result, (-1) r+p λ (ϕ [p] (1) -ϕ [p] (0)) = λ ϕ [2r-p-1] (0) and this leads to (C 2 ) since λ > 0.

As a conclusion, ϕ verifies both the ODE and the boundary conditions, hence verifies (B r λ ). Now, let us prove that (ii) =⇒ (i) .

For this, ϕ is assumed to be a solution of (B r λ ). This means that ϕ verifies the equality λ ϕ [2r] = (-1) r ϕ and all the boundary conditions. Under these assumptions, it must be proved that ϕ is an eigenfunction of T k r Sob with eigenvalue λ. To begin, it can be proved that ϕ has zero mean since:

1 0 ϕ(ξ) dξ = (-1) r λ 1 0 ϕ [2r] (ξ) dξ = (-1) r λ ϕ [2r-1] (1) -ϕ [2r-1] (0) = 0 .
The final equality is provided by the boundary values of ϕ [2r-1] in (C 3 ). Coming back to the main objective, the key idea of our proof is to write 2r] ](•) and to transform T k r Sob ϕ [2r] by repeating integration by parts in order to demonstrate that ϕ verifies the eigenvalue equation T k r Sob ϕ = λ ϕ. Before doing this, let us introduce the following notation:

[T k r Sob ϕ](•) = (-1) r λ [T k r Sob ϕ [
∀ g : [0 , 1] → R, g(1) -g(0) = [g(t)] 1 t=0 .
This will be notably useful to summarize the result of an integration by parts. The integral expression of T k r Sob ϕ [2r] is now splitted into three terms:

[T k r Sob ϕ [2r] ](x) = 1 0 k r Sob (x, ξ) ϕ [2r] (ξ) dξ = 1 0   r j=1 B j (x)B j (x ′ ) (j!) 2 + (-1) r+1 (2r)! B 2r (|x -ξ|)   ϕ [2r] (ξ) dξ = r j=1 1 0 B j (ξ) ϕ [2r] (ξ) dξ B j (x) (j!) 2 . . . + (-1) r+1 (2r)! x 0 B 2r (x -ξ) ϕ [2r] (ξ) dξ + 1 x B 2r (ξ -x) ϕ [2r] (ξ) dξ = r j=1 I j (j!) 2 B j (x) + (-1) r+1 (2r)! [ J -(x) + J + (x) ] . (F.6)
Whether for I j , J -(x) or J + (x), the integrand is the product of ϕ [2r] by a Bernoulli polynomial. Therefore, after repeating integration by parts as many times as necessary, these integrals can be calculated, or at least transformed into expressions where integration only concerns ϕ.

Integration by parts to transform I j It can be shown that:

∀ 1 ≤ j ≤ r, ∀ 1 ≤ k ≤ j, I j = k-1 l=0 (-1) l j! (j -l)! B j-l (t) ϕ [2r-l-1] (t) 1 t=0 . . . + (-1) k j! (j -k)! 1 0 B j-k (ξ) ϕ [2r-k] (ξ) dξ . (F.7)
The proof (by induction) relies on the formula of the derivatives of Bernoulli polynomials (see Appendix A.1.3). Then, taking k = j in Eq. (F.7) yields:

∀ 1 ≤ j ≤ r, I j = j l=0 (-1) l j! (j -l)! B j-l (t) ϕ [2r-l-1] (t) 1 t=0 . (F.8)
Integration by parts to transform J -(x) It can be shown that:

∀ 1 ≤ k ≤ 2r, J -(x) = k-1 l=0 (2r)! (2r -l)! B 2r-l (x -t) ϕ [2r-l-1] (t) x t=0 . . . + (2r)! (2r -k)! x 0 B 2r-k (x -ξ) ϕ [2r-k] (ξ) dξ . (F.9)
Once again, the proof (by induction) relies on the formula of the derivatives of Bernoulli polynomials. Then, taking k = 2r in Eq. (F.9) yields:

J -(x) = 2r-1 l=0 (2r)! (2r -l)! B 2r-l (x -t) ϕ [2r-l-1] (t) x t=0 + (2r)! x 0 ϕ(ξ) dξ . (F.10)
Integration by parts to transform J + (x)

It can be shown that:

∀ 1 ≤ k ≤ 2r, J + (x) = k-1 l=0 (-1) l (2r)! (2r -l)! B 2r-l (t -x) ϕ [2r-l-1] (t) 1 t=x . . . + (-1) k (2r)! (2r -k)! 1 x B 2r-k (ξ -x) ϕ [2r-k] (ξ) dξ . (F.11)
The proof (by induction) is similar to what was done for J -(x). Then, taking k = 2r in Eq. (F.11) yields:

J + (x) = 2r-1 l=0 (-1) l (2r)! (2r -l)! B 2r-l (t -x) ϕ [2r-l-1] (t) 1 t=x + (2r)! 1 x ϕ(ξ) dξ . (F.12)
The summation of Eq. (F.10) and (F.12) leads to: To switch from Eq. (F.13) to Eq. (F.14), one needs to perform the index shift k = 2r -l and to remember that ϕ has zero mean. In addition:

J -(x) + J + (x) = 2r-1 l=0 (2r)! (2r -l)! B 2r-l (x -t) ϕ [2r-l-1] (t) x t=0 . . . +(-1) l B 2r-l (t -x) ϕ [2r-l-1] (t) 1 t=x + 1 0 ϕ(ξ) dξ (F.13) = 2r k=1 (2r)! k! B k (x -t) ϕ [k-1] (t) x t=0 + (-1) 2r-k B k (t -x) ϕ [k-1] (t) 1 t=x (F.14) = 2r k=1 (2r)! k! -B k (x) ϕ [k-1] (0) + (-1) k B k (1 -x) ϕ [k-1] (
• With (C 1 ) and the symmetry properties of Bernoulli polynomials (see Appendix A.1.4), the sum S 1 (x) can be simplified: (F.17)

S 1 (x) =
After bringing together Eq. (F.6), (F.8), (F.15), (F. [START_REF] Borgonovo | Sensitivity Analysis: an Introduction for the Management Scientist[END_REF]) and (F.17), one has: The unanchored Sobolev space (of order r) is the RKHS obtained when H r ([0 , 1]) is equipped with the inner product ⟨•, •⟩ H r Sob defined in Eq. (3.6). This RKHS is denoted by H r Sob and its reproducing kernel is:

T k r Sob ϕ [2r] (x) =
K r Sob (x, x ′ ) = 1 + k r Sob (x, x ′ ) = 1 + r k=1 B k (x)B k (x ′ ) (k!) 2 + (-1) r+1 (2r)! B 2r (|x -x ′ |) = 1 + k r A (x, x ′ ) + k r B (x, x ′ ) .
Preliminary results on the two kernels k r A and k r B need to be stated.

Lemma F.4. For any r ≥ 1, the following statements are true:

(a) k r A and k r B are orthogonal kernels. (b) The RKHS F r A induced by k r A is:

F r A = f ∈ R [0,1] f (•) = r k=1 a k B k (•) with (a k ) 1≤k≤r ∈ R r ,
with inner product:

⟨•, •⟩ F r A : F r A × F r A -→ R f 1 (•) = r k=1 a k B k (•) , f 2 (•) = r k=1 α k B k (•) -→ r k=1 a k α k .
(c) The RKHS F r B induced by k r B is:

F r B = f ∈ R [0,1] f (•) = ∞ k=1 p k c r 2k (•) + q k s r 2k (•) with (p k ) k≥1 ∈ ℓ 2 (N * ) (q k ) k≥1 ∈ ℓ 2 (N * ) ,
with inner product:

⟨•, •⟩ F r B : F r B × F r B -→ R f 1 (•) = ∞ k=1 p k c r 2k (•) + q k s r 2k (•) , f 2 (•) = ∞ k=1 π k c r 2k (•) + ρ k s r 2k (•) -→ ∞ k=1 p k π k + q k ρ k . (d) F r A ∩ F r B = {0}.
Proof. The statement (a) can be proved very easily by using the zero-mean property of Bernoulli polynomials (see Appendix A.1.6). As regards the statements (b) and (c), they are direct consequences of Theorems 2.23 and 2.21 respectively. Only the statement (d) deserves further clarifications. For any f ∈ F r A ∩ F r B , it must be proved that f = 0. Since f ∈ F r A , one can write: The question is now to examine which conditions (on a 2 , . . . , a r ) are required so that the sequence (γ k ) k≥1 is square summable. A simple factorization allows to see that: It must be proved that the system (g r k ) k defined by Eq. (7.9) is truly ℓ 2 -linearly independent. Let us consider a square-summable sequence (γ k ) k which is rearranged in the following way:

f (•) =
γ k =
(γ k ) k := (a k ) 0≤k≤r ; (p k ) k≥1 ; (q k ) k≥1 .
This rearrangement allows for a better matching of the coefficients in (γ k ) k and the functions in (g r k ) k . Assuming that (γ k ) k is square summable amounts to assuming that the subsequences (p k ) k≥1 and (q k ) k≥1 are both in ℓ 2 (N * ). Note that the series defined as:

S(•) = a 0 ∈R + r k=1 a k B k (•) ∈F r A + ∞ k=1 p k c r 2k (•) + q k s r 2k (•) ∈F r B = a 0 + f r A (•) + f r B (•)
belongs to H r ([0 , 1]) in virtue of Eq. (F.24). Now, let us assume that the square-summable sequence (γ k ) k is such that S = 0. As the zero function belongs simultaneously to R, F r A and F r B , one can write:

S(•) = 0 + 0 + 0 = a 0 + f r A (•) + f r B (•)
Because of the direct sum, the decomposition of S into a sum of three sub-functions belonging respectively to R, F r A and F r B is unique. This leads to: Eq. (F.25) implies a 1 = . . . = a r = 0 because the polynomials ( B k ) 1≤k≤r are linearly independent. Furthermore, after turning back to the definition of the sinusoidal features, Eq. (F.26) may be rewritten in the following way:

a 0 = 0 , f r A (•) =
∞ k=1 p k (2kπ) r c 2k (•) + q k (2kπ) r s 2k (•) = 0 .
(F.27) Remember that the system {1 ; (c 2k ) k≥1 ; (s 2k ) k≥1 } is the Fourier basis of L 2 ([0 , 1]) and is therefore ℓ 2 -linearly independent (because of Parseval's identity). Additionally, as the sequences with general terms p k /(2kπ) r and q k /(2kπ) r are in ℓ 2 (N * ), Eq. (F.27) cannot be verified unless p k = 0 and q k = 0 for all k ≥ 1. Finally, (γ k ) k is only composed of zero coefficients. This proves that the system (g r k ) k is ℓ 2 -linearly independent.

F.5. Additional details for Section 7.3

Let K : X × X → R be a Mercer kernel, let ν be a probability measure with support X and let T K be the integral operator built from K and ν. Just as in Theorem 2.18, (λ i ) i≥1 and (ϕ i ) i≥1 respectively denote the eigenvalues of T K and an ONB of L 2 (X , ν) only composed of eigenfunctions of T K . In addition, let us assume that there also exists a series expansion of K based on a system (g i ) i≥1 of non-orthogonal but ℓ 2 -linearly independent features. In light of these assumptions, K can be decomposed in two different ways:

∀ x, x ′ ∈ X , K(x, x ′ ) = i≥1 λ i ϕ i (x) ϕ i (x ′ ) = i≥1 g i (x) g i (x ′ ) .
(F.28) Theorem 2.21 and Theorem 2.23 ensure that ( √ λ i ϕ i ) i≥1 and (g i ) i≥1 are two possible ONBs of H. However, as the functions (g i ) i∈I are not mutually orthogonal (in the L 2 -sense), they cannot be eigenfunctions of T K . The difference between the two families of basis functions can be summarized by the following inequality:

∥T K ∥ 2 HS = i≥1 λ 2 i ≥ i≥1 ∥g i ∥ 4 L 2 ,
with equality holding if and only if the system (g i ) i≥1 is L 2 -orthogonal. Let us demonstrate this result.

For the first equality, just write the Hilbert-Schmidt norm of T K with the ONB of eigenfunctions (ϕ i ) i≥1 :

∥T K ∥ 2 HS = i≥1 j≥1 |⟨T K e i , e j ⟩ L 2 | 2 = i≥1 λ 2 i .
Then, remember that it was stated in Eq. (2.10) that ∥T K ∥ 2 HS = ∥K∥ 2 L 2 . If using the second decomposition in Eq. (F.28), ∥K∥ 2 L 2 can be computed differently:

∥K∥ 2 L 2 = i≥1 g i ⊗ g i 2 L 2 = i≥1 g i ⊗ g i , j≥1 g j ⊗ g j L 2 = i≥1 j≥1
⟨g i ⊗ g i , g j ⊗ g j ⟩ L 2 with the dominated convergence theorem,

= i≥1 j≥1 ⟨g i , g j ⟩ 2 L 2
by definition of the inner product in L 2 (X 2 , ν ⊗2 ),

≥ i≥1 ∥g i ∥ 4 L 2
by only taking the terms such that i = j.

Before using Lebesgue's dominated convergence theorem (to switch from the first to the second line), it must be proved that the absolute value of the integrand belongs to L 1 (X 2 , ν ⊗2 ): A ) r≥1 converges uniformly to k ∞ B = 0. For the two convergence results, the idea of the proof is the same. The difference between the reproducing kernel of order r and the proposed asymptotic kernel must be uniformly bounded on [0 , 1]

X ×X   i≥1 g i (x) g i (ξ)
2 by a constant C(r) converging to zero as r tends to ∞. This may be summarized by:

∀ x, x ′ ∈ [0 , 1], |k r ⋆ (x, x ′ ) -k ∞ ⋆ (x, x ′ )| ≤ C(r) ---→ r→∞ 0 ,
where ⋆ is either A or B.

Proof of (a)

First, it must be checked that the series k ∞ A is convergent for any point (x, x ′ ) ∈ [0 , 1]. This can be easily achieved with the help of Eq. (8.3):

|k ∞ A (x, x ′ )| = ∞ k=1 B k (x) B k (x ′ ) (k!) 2 ≤ ∞ k=1 M + k k! 2 ≤ 16 ∞ k=1 1 2π 2k = 16 4π 2 -1 < ∞ .

The function k ∞

A is thus well-defined everywhere on [0 , 1] 2 . In addition, the remainder of the series can be uniformly bounded:

|k ∞ A (x, x ′ ) -k r A (x, x ′ )| = ∞ k=r+1 B k (x) B k (x ′ ) (k!) 2 ≤ ∞ k=r+1 M + k k! 2 ≤ 16 4π 2 -1 1 2π 2r ---→ r→∞ 0 .
(F. [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] This justifies that (k r A ) r≥1 uniformly converges to k ∞ A on [0 , 1] 2 . It is trivial to see that k ∞ A is a kernel because it is expressed as a convergent series of symmetric and separable functions. More precisely, k r A is a continuous kernel because of the uniform limit theorem.

Proof of (b)

The same kind of reasoning must be applied but it is even simpler because there is only one single Bernoulli polynomial to bound: This justifies that (k r B ) r≥1 uniformly converges to the zero kernel on [0 , 1] 2 .

F.7. Proof of Proposition 8.2

First, it must be noted that the series k≥0 a k B k (•) is uniformly convergent for any sequence (a k ) k≥0 ∈ ℓ 2 (N). Indeed, with the Cauchy-Schwarz inequality and the same upper bounding technique as in Eq. (F.29), one has: In the above equation, the integral over [0 , 1] and the summation over N are freely interchanged because the series converges uniformly on [0 , 1]. With Eq. (F.30), one has a 0 = 0. Let us now explain why all remaining coefficients (a i ) i≥1 are also zero. As S is a uniformly convergent series of polynomials (defined on a bounded interval), it can be easily justified that S ∈ C ∞ ([0 , 1]). In addition, the derivatives of S can be computed through term-by-term differentiation:

∀ i ≥ 1, S [i] (•) = ∞ k=0 a k B [i] k (•) = ∞ k=i a k B k-i (•) = ∞ k=0 a k+i B k (•) = 0 .
With the analytical expressions of α, β, γ 0 , γ 2 (x), γ - 3 (x), γ + 3 (x) and γ 4 (x), Eq. (G.4), (G.5) and (G.6) yield:

T k 2 Sob B 1 (x) = α B 1 (x) +

Theorem 4 . 4 .

 44 It is assumed that:(A1) The input variables X 1 , . . . , X d are mutually independent.

  u⊆{1,...,d} HSIC u = u⊆{1,...,d} v⊆u (-1) |u|-|v| HSIC(X v , Y ) . (4.6)

(4. 8 )

 8 If K Y is the Gaussian kernel, an ONB of H Y has already been provided in Example 4.3. The identification of an ONB of H r Sob seems to be the only remaining technical deadlock. According to Theorems 2.21 and 2.23, a solution may be found through the extraction of a feature map from K r

Theorem 5 .

 5 1 allows to disclose the Mercer feature map of K 1 Sob = 1 + k 1 Sob . One can see that it is composed of one constant feature and an infinite number of purely sinusoidal features. The Mercer decomposition of K 1 Sob paves the way to a feature-based characterization of the unanchored Sobolev space H 1 Sob . Importantly, Theorem 2.21 ensures that the system {1 ; (c k ) k≥1 } with ck (•) := c k (•)/kπ is an ONB of H 1 Sob . Remark 5.3. The Mercer expansion of the unanchored Sobolev kernel K 1

Remark 6 . 4 .

 64 Theorem 6.3 is an important result since the equivalence between the eigenvalue problem (S r λ ) and the boundary value problem (B r λ ) is clearly established. Differentiating the eigenvalue equation T K ϕ = λ ϕ in order to obtain an equivalent boundary value problem is a well-known proof technique to disclose the Mercer expansion of a given kernel. Other examples of boundary value problems obtained in this way are provided in the Appendix B.

w

  i exp(z i t) with w := (w 1 , . . . , w 2r ) ∈ C 2r .(6.1)

Proposition 8 . 1 .

 81 For any r ≥ 1, let k r A and k r B be the kernels introduced in Section 7.1 to decompose k r Sob . (a) The kernel sequence (k r A ) r≥1 converges uniformly to the continuous kernel k ∞ A defined by:

  773-780). A.2.1. Definition Definition A.1. There exists a unique family (L k ) k≥0 of orthogonal polynomial functions in L 2 ([-1 , 1]) such that deg(L k ) = k and L k (1) = 1 for all k ≥ 0. They are called Legendre polynomials and they form a complete and orthogonal system in L 2 ([-1 , 1]).

Remark A. 2 .

 2 Legendre polynomials are not an ONB of L 2 ([-1 , 1]) because they do not have unit norm:

1 L 2 =

 12 The shift function σ is a linear mapping from [a , b] into [-1 , 1] that helps go back to the standard definition domain of Legendre polynomials. Using σ allows to preserve orthogonality among the polynomials (P k ) k≥1 while the multiplying factor ∥L k • σ∥ -(2k + 1)/(b -a) ensures that the polynomial functions (P k ) k≥0 are normalized in L 2 ([a , b]). In the specific case where a = 0 and b = 1, the polynomials (P k ) k≥1 are called the shifted Legendre polynomials:

8 35 x 4 -B. 3 . 1 .

 431 30 x 2 + 3 ; • P 4 (x) = 70 x 4 -140 x 3 + 90 x 2 -20 x + 1 .B.3. Sobolev kernels Covariance kernel of the Wiener process

  2r and ϕ 1k (x) = √ 2 cos(2kπx) ;• ∀ k ≥ 1, λ 2k = 1 (2kπ) 2r and ϕ 2k (x) = √ 2 sin(2kπx) .

  u⊆{1,...,d} HSIC u = u⊆{1,...,d} v⊆u (-1) |u|-|v| HSIC(X v , Y ) .

  has η 3 (ξ) = ξ 8 e √ 3ξ det [N 3 (ξ)] with N 3 (ξ) the matrix obtained after operating all factorizations.

1 )

 1 k (0) ϕ [k-1] (x) -(-1) k B k (0) ϕ [k-1] (x) = S 1 (x) + S 2 (x) .(F.15)

1 ] ( 1 )•

 11 -ϕ [k-1] (0) B k (x) . (F.16) With the boundary values of Bernoulli polynomials (see Appendix A.1.4), one can see that the sum S 2 (x) is actually composed of only one non-zero term: S 2 (x) = -(2r)! ϕ(x) .

F. 4 . 2 .

 42 ) 2 B k (x) + (-1) r+1 r k=1 1 k! ϕ [k-1] (1) -ϕ [k-1] (0) B k (x) + (-1) r ϕ(x)A few more details on the sub-kernel decomposition

  i! with (a i ) 1≤i≤r ∈ R r . Now, let us see why all coefficients are zero. First, let us assume that a 1 ̸ = 0 for the sake of contradiction.It is said in AppendixA.1.4 that B 1 (1) = -B 1 (0) = 1/2 and B k (1) = B k (0) for all k ≥ 2 and this leads to f (1) -f (0) = a 1 ̸ = 0. Therefore, f cannot belong to F rB because it only contains 1-periodic functions. Since there is a contradiction with f ∈ F r B , it must be that a 1 = 0. Then, let us write f in the formalism of F r B : -jπ/2) (2kπ) j with the Fourier series expansion of B j , j ϵ k (x) withϵ k := if j = 0 mod 4 -s 2k if j = 1 mod 4 c 2k if j = 2 mod 4 s 2k if j = 3 2kπ) r-j   ϵ k (x) (2kπ) r = ∞ k=1 γ k ϵ k (x) (2kπ) r .

a 2 F. 4 . 3 .

 243 (2kπ) r-2 if a 2 ̸ = 0 a 3 (2kπ) r-3 if a 2 = 0 but a 3 ̸ = 0 . . . . . . a r if a 2 = . . . = a r-1 = 0 but a r ̸ = 0 0 if a 2 = . . . = a r-1 = a r = 0 .The only way to make (γ k ) k≥1 be square summable is thus to take a 2 = . . . = a r = 0. This leads to f = 0 . □ Knowing (a), one has R ∩ F r A = {0} and R ∩ F r B = {0}. Since (d) also ensures that F r A ∩ F r B = {0}, the sum of the three sub-RKHSs is direct. With the help of Theorem F.3, it is then straightforward to obtain:H r ([0 , 1]) = R ⊕ F r A ⊕ F r B .(F.24) Proof of the main result

  (•) + q k s r 2k (•) = 0 . (F.26)

2 < 2 .F. 6 . 1

 2261 , x) K(ξ, ξ) dν(x) dν(ξ) = X K(x, x) dν(x) ∞ .The first line is obtained by applying the Cauchy-Schwarz inequality in ℓ 2 (N * ). The two equalities in the second line stem from Eq. (F.28) and Fubini's theorem. The final quantity is indeed finite because K is a Mercer kernel. Eventually, one has∥T K ∥ 2 HS = i≥1 λ 2 i = ∥K∥ 2 L 2 ≥ i≥1 ∥g i ∥ 4 L Proof of Proposition 8.It must be proved that: (a) The sequence (k r A ) r≥1 converges uniformly to the continuous kernel k ∞ A = k≥1 B k ⊗ B k . (b) The sequence (k r

1 ((

 1 |k r B (x, x ′ )| = 2r)! |B 2r (|x -x ′ |)| ≤ M + 2r

B

  us assume that there exist a sequence (a k ) k≥0 ∈ ℓ 2 (N) such that S(•) := k≥0 a k B k (•) = 0. Under this assumption, the mean value of the function S : [0 , 1] → R is equal to zero. The zero-mean property of Bernoulli polynomials (see Appendix A.1.6) leads to: k (x) dx = a 0 .(F.30)

=

  g A (x) + g B (x) .

( a )

 a Histogram of the estimated eigenvalues. (b) Curves of the estimated eigenfunctions.

Figure 1 .

 1 Figure 1. Estimation of the eigenvalues and eigenfunctions involved in the Mercer decomposition of the kernel k 1Sob . The KFA method is performed with n = 500 sample points.

( a )

 a Estimation of the decay rate for the kernel k 1 Sob after a logarithmic transformation. (b) Estimation of the decay rate for the kernel k 2 Sob after a logarithmic transformation.

Figure 2 .

 2 Figure 2. Regression-based eigendecay analysis for the kernel k r Sob (with r ∈ {1, 2}). The eigenvalues are first estimated by KFA (with n = 500 points).

( a )

 a Histogram of the estimated eigenvalues. (b) Curves of the estimated eigenfunctions.

Figure 3 .

 3 Figure 3. Estimation of the eigenvalues and eigenfunctions involved in the Mercer decomposition of the kernel k 2Sob . The KFA method is performed with n = 500 sample points.

( a )

 a Curves of the estimated eigenfunctions. (b) Approximation of the theoretical eigenfunctions.

Figure 4 .

 4 Figure 4. Brute-force estimation of the eigenfunctions involved in the Mercer decomposition of the kernel k 5Sob . The KFA method is performed with n = 3000 sample points.

ANOVA means ANalysis Of VAriance. This acronym is encountered in many different areas of statistics[START_REF] Archer | Sensitivity measures, ANOVA-like techniques and the use of bootstrap[END_REF][START_REF] Ginsbourger | On ANOVA decompositions of kernels and Gaussian random field paths[END_REF][START_REF] Gu | Smoothing Spline ANOVA Models[END_REF][START_REF] Scheffe | The Analysis of Variance[END_REF][START_REF] Stitson | Support vector regression with ANOVA decomposition kernels[END_REF].

In this work, the standard uniform distribution refers to the uniform distribution on [0 , 1]. For convenience, this distribution will sometimes be denoted by U ([0 , 1]) in equations.

In the rest of the paper, the Greek letter φ will be exclusively reserved for denoting the Mercer feature maps.

The original definition of Sobol' indices[START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] is based on the Sobol'-Hoeffding decomposition of the input-output numerical simulator g : X 1 × . . . × X d → Y. Since the uniqueness of this decomposition is only guaranteed in the presence of mutually independent inputs, Sobol' indices are well defined only under this assumption. Over the past fifteen years, there have been several attempts to extend Sobol' indices to correlated inputs[START_REF] Chastaing | Generalized Sobol' sensitivity indices for dependent variables: Numerical methods[END_REF][START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF][START_REF] Owen | On Shapley value for measuring importance of dependent inputs[END_REF] but none of them has proved sufficiently satisfactory.

This choice is made for the sake of simplicity. However, note that everything remains true when the random objects live in (possibly different) separable spaces[START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF].

https://www.wolframalpha.com

Acknowledgments

This research work is part of the SAMOURAI 10 project funded by the French National Research Agency (ANR-20-CE46-0013). This financial support is gratefully acknowledged. The authors would also like to thank

Asymptotic Mercer feature maps

The question is whether the Mercer decomposition of k ∞ Sob has a closed-form expression relying on basic math functions. Since k ∞ Sob = k ∞ A is the limit kernel of the sequence (k r A ) r≥1 , the Mercer expansion of these kernels could help guess the asymptotic Mercer decomposition. Unfortunately, it was explained in Section 7.2.1 that the eigenvalue problem related to T k r A is more or less difficult to solve depending on the value of r. The resolution is trivial for r ∈ {1, 2} because the kernel definitions

are almost Mercer representations (up to the normalization of B 1 and B 2 ). On the contrary, the eigenvalue equation is much more difficult to solve when r ≥ 3.

• For r = 3, the Mercer decomposition of k 3 A can be obtained after long and exhausting hand calculations. The three resulting eigenfunctions are polynomials (which is quite normal since the RKHS induced by k r A is the space of all zero-mean polynomials with degree at most 3) but they are not P 1 , P 2 and P 3 .

• For r ≥ 4, there is no point in trying to apply a brute-force approach based on hand calculations.

Therefore, the Mercer decomposition of k ∞ Sob cannot be derived from a general result on the Mercer decomposition of k r A . A simple method in order to obtain a definitive answer is to compute T k ∞ Sob P 1 and to prove that it cannot be proportional to P 1 . For the sake of convenience, the integral calculation is carried out with B 1 ∝ P 1 :

An integration by parts is used to switch from the first line to the second line. Then, the properties of Bernoulli numbers (see Appendix A.1.2) enable simplifications leading to the final expression. According to Proposition 8.2, the polynomials ( B k ) k≥1 are ℓ 2 -linearly independent. This allows to prove that T k ∞ Sob B 1 and B 1 cannot be proportional. As a consequence, B 1 is not an eigenfunction of T k ∞ Sob and neither is P 1 . With the same technique but longer calculations, it could be proved that P 2 is not an eigenfunction of T k ∞ Sob , and so on.

Once again, the shifted Legendre polynomials (P k ) k≥1 appear not to be the eigenfunctions of the integral operator under study. Thus, the asymptotic framework does not explain why those polynomials provide so accurate approximations of the eigenfunctions involved in the Mercer expansions of the kernels (k r A ) r≥1 and (k r Sob ) r≥2 . Identifying the Mercer expansion of k ∞ Sob thus remains an open question.

Appendix A. Remarkable properties of two families of polynomials

Bernoulli polynomials

This section brings together various results on Bernoulli polynomials which are notably used in the proofs of Appendix F. A similar table dedicated to Bernoulli polynomials can be found in [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] (see Section 23,. For a much more comprehensive overview of existing results on Bernoulli polynomials, it is strongly recommended to have a look at [START_REF] Kouba | Lecture notes: Bernoulli polynomials and applications[END_REF].

A.1.1. Explicit expressions for low degrees

Bernoulli polynomials (B n ) n≥0 are all defined on [0 , 1]. The analytical expressions of the first Bernoulli polynomials are detailed below:

A.1.2. Bernoulli numbers

For any n ≥ 0, the value B n (0) is called the n-th Bernoulli number. Here is a list of the first ones:

; B 9 (0) = 0 ; B 10 (0) = 5 66 ; B 11 (0) = 0 . . .

In particular, Bernoulli numbers verify the following properties:

A.1.4. Symmetries and boundary values

Eigenvalues and eigenfunctions:

(H k ) k≥1 are the Hermite polynomials [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] (see pages 773-780).

B.2. Laplace kernel

Equivalent boundary value problem:

For any λ > 0,

Eigenvalues and eigenfunctions:

.

The sequences (w 1k ) k≥1 and (w 2k ) k≥1 are the (unknown) solutions of the two following non-linear equations:

Appendix D. Kernel feature analysis

Let K : X × X → R be a Mercer kernel, let ν be a probability measure with support X and let T K be the integral operator built from K and ν. The objective here is to give additional details regarding the numerical procedure described in Section 5.1 to solve the eigenvalue problem:

Originally, the approach described below was known as the Nyström method [START_REF] Spence | On the convergence of the Nyström method for the integral equation eigenvalue problem[END_REF][START_REF] Spence | Error bounds and estimates for eigenvalues of integral equations[END_REF]. In more recent papers, it is also called kernel principal component analysis (kernel PCA) [START_REF] Liang | Eigen-analysis of nonlinear PCA with polynomial kernels[END_REF][START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF][START_REF] Schölkopf | Kernel principal component analysis[END_REF] or kernel feature analysis (KFA) [START_REF] Jiang | Accelerated kernel feature analysis[END_REF][START_REF] Smola | Sparse kernel feature analysis[END_REF]. In this work, it was decided to use the latter terminology.

D.1. Discretization of the initial eigenvalue problem

Let f ∈ L 2 (X , ν) be an eigenfunction of T K associated to λ > 0. For now, the idea is to solve the eigenvalue equation T K f = λ f without paying attention to the fact that the eigenfunction mentioned in Theorem 2.18 are expected to have unit L 2 -norm. This renormalization issue will be addressed later. The starting point is to draw a n-sample x sim := (x i ) 1≤i≤n from ν. For any x ∈ X , the pointwise equality [T K f ](x) = λ f (x) can be discretized with the Monte Carlo method:

∼ ν .

Taking x = x i for i ∈ {1, . . . , n} yields a system of n linear equations:

After denoting by v i := f (x i ) the unknown values taken by f at the simulated points, the above system of equations may be rewritten as a matrix equation based on the simulated Gram matrix K n :

The eigenvalue problem related to T K is thus transformed into a matrix eigenvalue problem which can be solved numerically (as long as n is not too large).

Remark D.1. Let ψ : X → F be one possible feature map of K. Knowing that K(x i , x j ) = ⟨ψ(x i ), ψ(x j )⟩ F , the Gram matrix K n can be envisioned as the covariance matrix of the n feature functions ψ(x i ) in the feature space F. Furthermore, since KFA is basically an eigenanalysis of K n , it is often called kernel PCA.

D.2. Estimation of the eigenvalues and eigenfunctions

Let (γ k ) 1≤k≤n denote the eigenvalues of K n and let (v k ) 1≤k≤n denote the corresponding eigenvectors. By convention, it is assumed that γ 1 ≥ . . . ≥ γ n ≥ 0. Note that the eigenvalues are all non-negative because K n is a symmetric positive semi-definite matrix. With Eq. (D.1), it is clear that the k-th largest eigenvalue of T K can be estimated by λk = γ k /n. Now, the question is how to estimate the eigenfunctions (ϕ k ) 1≤k≤n of T K from the eigenvectors (v k ) 1≤k≤n of the matrix K n . As already said when establishing Eq. (D.1), the k-th eigenvector v k := [v ik ] 1≤i≤n may be seen as a discrete estimate of an eigenfunction f k associated to λ k :

Appendix E. Analytical resolution of the boundary value problem E.1. Exact resolution for r = 1

In this case, the boundary value problem (B r λ ) is restricted to:

To begin with, any real-valued solution of this ODE may be written as:

The boundary conditions on ϕ ′ yield β = 0 and sin(ξ) = 0. Hence, there are countably infinitely many possible solutions for ξ and they are all given by ξ k = kπ (with k ≥ 1 since λ > 0). The sequence of eigenvalues (λ k ) k≥1 directly follows with λ k = 1/ξ 2 k = 1/(kπ) 2 . For each eigenvalue λ k , a pair (α k , β k ) ∈ R 2 is used to parametrize the corresponding eigenfunction ϕ k . One has β k = 0 but α k remains free of any constraint. This is natural since the solution space is expected to be an eigenspace (i.e. a linear subspace of dimension at least 1). More precisely, the eigenspace associated to λ k can be expressed as:

The solution space of (B 1 λ ) contains infinitely many solutions but the Mercer decomposition of k 1 Sob requires L 2 -normalized eigenfunctions. This leads to take α k = ± √ 2 and ϕ k : t ∈ [0 , 1] → √ 2 cos(kπt). Everything is therefore consistent with the Mercer expansion of k 1 Sob already found in Section 5.

E.2. Partial resolution for r = 3

For r = 3, the boundary value problem is given by: λϕ [6] + ϕ = 0 with

ϕ [3] (0) = ϕ [3] (1) = 0 ϕ(0) -ϕ(1) = ϕ [5] (0) = ϕ [5] (1) [4] (0) = ϕ [4] (1)

and λ > 0 .

Any real-valued solution of this ODE may be written as:

2 ξt) sin ( 1 2 ξt) . . .

(E.1)

The framed coefficients in N ∞ 3 (ξ) correspond to asymptotic approximations as opposed to all other coefficients which correspond to the original terms in N 3 (ξ). It is much more convenient to handle N ∞ 3 (ξ) because the Appendix F. Proofs

F.1. Proof of Theorem 2.23

The idea is to prove that K is indeed the reproducing kernel of the RKHS defined in Theorem 2.23. First, it must be checked that all functions in H are well-defined. For any h ∈ H, there exists (a i ) i∈I ∈ ℓ 2 (I) such that h(•) = i∈I a i g i (•) and the convergence of the series must be justified at any point x ∈ X . For this, simply apply the Cauchy-Schwarz inequality:

The vector space H from Eq. (2.17) is therefore composed of well-defined functions. The map ⟨•, •⟩ H : H×H → R defined in Eq. (2.18) is also well-defined:

• The sequences (a i ) i∈I and (b i ) i∈I used to write h 1 and h 2 are uniquely defined since the functions (g i ) i∈I are assumed to be ℓ 2 -linearly independent.

• The series i∈I a i b i is always convergent because the sequences (a i ) i∈I and (b i ) i∈I belong to ℓ 2 (I). From this, it can be easily proved that (H, ⟨•, •⟩ H ) is a Hilbert space. In particular, the normed space (H, ∥•∥ H ) inherits from the completeness of ℓ 2 (I). Now, it only remains to prove that K is a reproducing kernel for the constructed Hilbert space. For any x ∈ X , one has:

which allows to prove that K(•, x) belongs to H. In addition, the reproducing property is verified:

(H, ∥•∥ H ) is therefore a Hilbert space for which K is reproducing. By Theorem 2.5, we know that it is the only one, which ends the proof.

F.2. Proof of Proposition 6.1

) and that ϕ is a solution of the ODE recalled below:

First, let us prove that ϕ has zero mean:

The next step is to demonstrate that α 1 = . . . = α r = 0.

Nullity of α r

Given the definition of α r , it follows that:

Note that:

• Eq. (F.18) is obtained after replacing I r by the sum given in Eq. (F.8). Then, the term indexed by l = r is taken out of the sum and it is exactly the opposite of the last term.

• Eq. (F. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]) is obtained after taking the term indexed by l = r -1 out of the sum. In fact, this term is equal to zero because of (C 1 ).

• In Eq. (F.20), the sum indexed by even integers is equal to zero for two reasons. Firstly, all the associated Bernoulli polynomials have equal boundary values (see Appendix A.1.4). Secondly, all the functions ϕ [2r-l-1] have equal boundary values according to (C 3 ).

• In Eq. (F.20), the sum indexed by odd integers is equal to zero because all the associated Bernoulli polynomials have zero boundary values (see Appendix A.1.4).

Nullity of α k with k ∈ {1, . . . , r -1} Given the definition of α k , it follows that:

Note that:

• Eq. (F.21) is obtained after replacing I j by the sum given in Eq. (F.8). Then, the term indexed by l = r is taken out of the sum. This term is always equal to zero because of (C 1 ) when k = r -1 and because of (C 3 ) when 1 ≤ k ≤ r -2.

• Eq. (F.22) is obtained after taking the term indexed by l = r -1 out of the sum. The entire expression in the second line of Eq. (F. [START_REF] Campbell | Kernel methods: a survey of current techniques[END_REF]) is equal to zero because of (C 2 ) and (C 3 ).

• In Eq. (F.23), the sum indexed by even integers is equal to zero for two reasons. Firstly, all the associated Bernoulli polynomials have equal boundary values (see Appendix A. 1.4). Secondly, all the functions ϕ [2r-l-1] have equal boundary values according to (C 3 ).

• In Eq. (F.23), the sum indexed by odd integers is equal to zero because all the associated Bernoulli polynomials have zero boundary values (see Appendix A.1.4).

As a conclusion, having α 1 = . . . = α r = 0 finally leads to T k r Sob ϕ [2r]=(-1) r ϕ and then to T k r Sob ϕ = λ ϕ.

F.4. Proof of Proposition 7.3

Sum of kernels

In this proof, a fundamental theorem regarding the RKHS induced by the sum of two kernels is required.

Theorem F.3. Let X ⊆ R be an interval. Let K 1 : X × X → R and K 2 : X × X → R be two kernels with respective RKHSs denoted by H 1 and H 2 . The RKHS H 12 induced by the kernel K 1 + K 2 is:

If H 1 ∩ H 2 = {0}, one must write H 12 = H 1 ⊕ H 2 and the inner product is:

Appendix G. Supplementary material G.1. Proof of Theorem 5.1

It only remains to prove that T k 1 Sob c k (•) = c k (•)/(kπ) 2 for any k ≥ 1. This can be done by simply calculating the integral T k 1 Sob c k (x) for a given point x ∈ [0 , 1]. This is not particularly difficult but the calculations are tedious and they must be carried out patiently. The main calculation steps are provided below. First, one has: It can be proved that α = 1 (kπ) 2 (-1) k -1 .

Then, the integral β(x) may be divided into three terms: It can be proved that β 2 (x) = 2 (kπ) 2 (-1) k (1 -x) + x and β 0 = 0 .

To eliminate the absolute value in the integral expression of β 1 (x), one may write:

With the analytical expressions of α, β 0 , β - 1 (x), β + 1 (x) and β 2 (x), Eq. (G.1), (G.2) and (G.3) yield:

G.2. Additional details for Section 5.2.2

It only remains to prove that T k 2 Sob B 1 (•) = g A (•) + g B (•) as stated in Eq. (5.4). This can be done by simply calculating the integral [T k 2 Sob B 1 ](x) for a given point x ∈ [0 , 1]. This is not particularly difficult but the calculations are very tedious. To save time, most integrals were computed with the online integral calculator proposed by Wolfram|Alpha 15 . First, one has: To eliminate the absolute value in the integral expression of γ 3 (x), one may write: