Pierre Alquier
email: alquier@math.jussieu.fr

Lpma

Gérard Biau
email: gerard.biau@upmc.fr

Sparse Single-Index Model

Keywords: Single-index model, sparsity, regression estimation, PAC-Bayesian, oracle inequality, reversible jump Markov chain Monte Carlo method. 2010 Mathematics Subject Classification: 62G08, 62G05, 62G20

Let (X, Y ) be a random pair taking values in R p × R. In the socalled single-index model, one has Y = f ⋆ (θ ⋆T X) + W , where f ⋆ is an unknown univariate measurable function, θ ⋆ is an unknown vector in R d , and W denotes a random noise satisfying E[W |X] = 0. The single-index model is known to offer a flexible way to model a variety of high-dimensional real-world phenomena. However, despite its relative simplicity, this dimension reduction scheme is faced with severe complications as soon as the underlying dimension becomes

Introduction

Let D n = {(X 1 , Y 1 ), . . . , (X n , Y n )} be a collection of independent observations, distributed as a generic independent pair (X, Y ) taking values in R p ×R and satisfying EY 2 < ∞. Throughout, we let P be the distribution of (X, Y ), so that the sample D n is distributed according to P ⊗n . In the regression function estimation problem, the goal is to use the data D n in order to construct an estimate r n : R p → R of the regression function r

Y = f ⋆ (θ ⋆T X) + W, (1.1)
where f ⋆ is an unknown univariate measurable function. This is the celebrated single-index model, which is recognized as a particularly useful variation of the linear formulation and can easily be interpreted: The model changes only in the direction θ ⋆ , and the way it changes in this direction is described by the function f ⋆ . This model has applications to a variety of fields, such as discrete choice analysis in econometrics and dose response models in biometrics, where high-dimensional regression models are often employed. There are too many references to be included here, but the monographs of McCullagh and Nelder [START_REF] Mccullagh | Generalized Linear Models[END_REF] and Horowitz [START_REF] Horowitz | Semiparametric Methods in Econometrics[END_REF] together with the references [START_REF] Härdle | Optimal smoothing in single-index models[END_REF][START_REF] Ichimura | Semiparametric least squares (SLS) and weighted SLS estimation of single-index models[END_REF][START_REF] Delecroix | On semiparametric Mestimation in single-index regression[END_REF][START_REF] Dalalyan | A new algorithm for estimating the effective dimension-reduction subspace[END_REF][START_REF] Lopez | Single-index regression models with right-censored responses[END_REF] will provide the reader with good introductions to the general subject area.

One of the main advantages of the single-index model is its supposed ability to deal with the problem of high dimension (Bellman [6]). It is known that estimating the regression function is especially difficult whenever the dimension p of X becomes large. As a matter of fact, the optimal mean square convergence rate n -2k/(2k+p) for the estimation of a k-times differentiable regression function converges to zero dramatically slowly if the dimension p is large compared to k. This leads to an unsatisfactory accuracy of estimation for moderate sample sizes, and one possibility to circumvent this problem is to impose additional assumptions on the regression function. Thus, in particular, if r(x) = f ⋆ (θ ⋆T x) holds for every x ∈ R p , then the underlying structural dimension of the model is 1 (instead of p) and the estimation of r can hopefully be performed easier. In this regard, Gaïffas and Lecué show in [START_REF] Gaïffas | Optimal rates and adaptation in the singleindex model using aggregation[END_REF] that the optimal rate of convergence over the single-index model class is n -2k/(2k+1) (instead of n -2k/(2k+p) ), thereby answering a conjecture of Stone [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF].

Nevertheless, practical estimation of the link function f ⋆ and the index θ ⋆ still requires a degree of statistical smoothing. Perhaps the most common approach to reach this goal is to use a nonparametric smoother (for instance, a kernel or a local polynomial method) to construct an approximation fn of f ⋆ , then substitute fn into an empirical version R n (θ) of the mean square error R(θ) = E[Yf (θ T X)] 2 , and finally choose θn to minimize R n (θ) (see e.g. Härdle, Hall and Ichimura [START_REF] Härdle | Optimal smoothing in single-index models[END_REF] and Delecroix, Hristache and Patilea [START_REF] Delecroix | On semiparametric Mestimation in single-index regression[END_REF] where the procedure is discussed in detail). The rationale behind this type of two-stage approach, which is asymptotic in spirit, is that it produces a √ n-consistent estimate of θ, thereby devolving the difficulty to the simpler problem of computing a good estimate for the one-dimensional function f ⋆ . However, the relative simplicity of this strategy is accompanied by severe difficulties (overfitting) when the dimension p becomes larger than the number of observations n. Estimation in this setting (called "p larger than n" paradigm) is generally acknowledged as an important challenge in contemporary statistics, see e.g. the recent monograph of Bühlmann and van de Geer [START_REF] Bühlmann | Statistics for High-Dimensional Data[END_REF]. In fact, this drawback considerably reduces the ability of the single-index model to behave as an effective dimension reduction technique.

On the other hand, there is empirical evidence that many signals in highdimensional spaces admit a sparse representation. As an example, wavelet coefficients of images often exhibit exponential decay, and a relatively small subset of all wavelet coefficients allow for a good approximation of the original image. Such signals have few nonzero coefficients and can therefore be described as sparse in the signal domain (see for instance [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]). Similarly, recent advances in high-throughput technologies-such as array comparative genomic hybridization-indicate that, despite the huge dimensionality of problems, only a small number of genes may play a role in determining the outcome and be required to create good predictors ( [START_REF] Van't Veer | Gene expression profiling predicts clinical outcome of breast cancer[END_REF] for instance). Sparse estimation is playing an increasingly important role in the statistics and machine learning communities, and several methods have recently been developed in both fields, which rely upon the notion of sparsity (e.g. penalty methods like the Lasso and Dantzig selector, see [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Candès | The Dantzig selector: Statistical estimation when p is much larger than n[END_REF][START_REF] Bunea | Sparsity oracle inequalities for the Lasso[END_REF][START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF] and the references therein).

In the present document, we consider the single-index model (1.1) from a sparsity perspective, i.e., we assume that θ ⋆ has only a few coordinates different from 0. In the dimension reduction scenario we have in mind, the ambient dimension p can be very large, much larger than the sample size n, but we believe that the representation is sparse, i.e., that very few coordinates of θ ⋆ are nonzero. This assumption is helpful at least for two reasons: If p is large and the number of nonzero coordinates is small enough, then the model is easier to interpret and its efficient estimation becomes possible. Our setting is close in spirit of the approach of Cohen, Daubechies, DeVore, Kerkyacharian and Picard [START_REF] Cohen | Capturing ridge functions in high dimension from point queries[END_REF], who study approximation from queries of functions of the form f (θ T x), where θ is approximately sparse (in the sense that it belongs to a weak-ℓ p space). However, these authors dot not provide any statistical study of their model. Our modus operandi will rather rely on the so-called PAC-Bayesian approach, originally developed in the classification context by Shawe-Taylor and Williamson [START_REF] Shawe-Taylor | A PAC analysis of a Bayes estimator[END_REF], McAllester [START_REF] Mcallester | Some PAC-Bayesian theorems[END_REF] and Catoni [START_REF] Catoni | Statistical Learning Theory and Stochastic Optimization[END_REF][START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF]. This strategy was further investigated for regression by Audibert [START_REF] Audibert | Aggregated estimators and empirical complexity for least square regression[END_REF] and Alquier [START_REF] Alquier | PAC-Bayesian bounds for randomized empirical risk minimizers[END_REF] and, more recently, worked out in the sparsity framework by Dalalyan and Tsybakov [START_REF] Dalalyan | Aggregation by exponential weighting, sharp PAC-Bayesian bounds and sparsity[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF] and Alquier and Lounici [START_REF] Alquier | PAC-Bayesian bounds for sparse regression estimation with exponential weights[END_REF]. The main message of [START_REF] Dalalyan | Aggregation by exponential weighting, sharp PAC-Bayesian bounds and sparsity[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF][START_REF] Alquier | PAC-Bayesian bounds for sparse regression estimation with exponential weights[END_REF] is that aggregation with a properly chosen prior is able to deal nicely with the sparsity issue. Contrary to procedures such as the Lasso, the Dantzig selector and other penalized least square methods, which achieve fast rates under rather restrictive assumptions on the Gram matrix associated to the predictors, PAC-Bayesian aggregation requires only minimal assumptions on the model. Besides, it is computationally feasible even for a large p and exhibits good statistical performance.

The paper is organized as follows. In Section 2, we first set out some notation and introduce the single-index estimation procedure. Then we state our main result (Theorem 2.1), which offers a sparsity oracle inequality more powerful than the best known oracle inequalities for other common procedures of single-index recovery. Section 3 is devoted to the practical implementation of the estimate via a reversible jump Markov chain Monte Carlo (MCMC) algorithm, and to numerical experiments on both simulated and real-life data sets. In order to preserve clarity, proofs have been postponed to Section 4 and the description of the MCMC method in its full length is given in the Appendix Section 5.

Note finally that our techniques extend to the case of multiple-index models, of the form

Y = f ⋆ (θ ⋆T 1 X, . . . , θ ⋆T m X) + W,
where the underlying structural dimension m is supposed to be larger than 1 but substantially smaller than p. However, to keep things simple, we let m = 1 and leave the reader the opportunity to adapt the results to the more general situation m ≥ 1.

Sparse single-index estimation

We start this section with some notation and basic requirements.

Notation

Throughout the document, we suppose that the recorded data D n is generated according to the single-index model (1.1). More precisely, for each i = 1, . . . , n,

Y i = f ⋆ (θ ⋆T X i ) + W i ,
where f ⋆ is a univariate measurable function, θ ⋆ is a p-variate vector, and W 1 , . . . , W n are independent copies of W . We emphasize that it is implicitly assumed that the observations are drawn according to the true model under study. Therefore, this casts our problem in a nonparametric setting rather than in a classical PAC-Bayesian one.

Recall that, in model (1.1), E[W |X] = 0 and, consequently, that EW = 0. However, the distribution of W (in particular, the variance) may depend on X. We shall not precisely specify this dependence, and will rather require the following condition on the distribution of W .

Assumption N. There exist two positive constants σ and L such that, for all integers k ≥ 2,

E |W | k | X ≤ k! 2 σ 2 L k-2 .
Observe that Assumption N holds in particular if W = Φ(X)ε, where ε is a standard Gaussian random variable independent of X and Φ(X) is almost surely bounded.

Let θ 1 denote the ℓ 1 -norm of the vector θ = (θ 1 , . . . , θ p ) T , i.e., θ 1 = p j=1 |θ j |. Without loss of generality, it will be assumed throughout the document that the index θ ⋆ belongs to S p 1,+ , where S p 1,+ is the set of all θ ∈ R p with θ 1 = 1 and θ j(θ) > 0, where j(θ) is the smallest j ∈ {1, . . . , p} such that θ j = 0.

We will also require that the random variable X is almost surely bounded by a constant which, without loss of generality, can be taken equal to 1. Moreover, it will also be assumed that the link function f ⋆ is bounded by some known positive constant C. Thus, denoting by X ∞ the supremum norm of X and by f ⋆ ∞ the functional supremum norm of f ⋆ over [-1, 1], we set: Assumption B. The condition X ∞ ≤ 1 holds almost surely and there exists a positive constant C larger than 1 such that f ⋆ ∞ ≤ C. Remark 2.1 To keep a sufficient degree of clarity, no attempt was made to optimize the constants. In particular, the requirement C ≥ 1 is purely technical. It can easily be removed by replacing C by max[C, 1] throughout the document.

In order to approximate the link function f ⋆ , we shall use the vector space F spanned by a given countable dictionary of measurable functions {ϕ j } ∞ j=1 . Put differently, the approximation space F is the set of (finite) linear combinations of functions of the dictionary. Each ϕ j of the collection is assumed to be defined on [-1, 1] and to take values in [-1, 1]. To avoid getting into too much technicalities, we will also assume that each ϕ j is differentiable and such that, for some positive constant ℓ, ϕ ′ j ∞ ≤ ℓj. This assumption is satisfied by the (non-normalized) trigonometric system

ϕ 1 (t) = 1, ϕ 2j (t) = cos(πjt), ϕ 2j+1 (t) = sin(πjt), j = 1, 2, . . .
Finally, for any measurable f : R p → R and θ ∈ S p 1,+ , we let

R(θ, f ) = E Y -f (θ T X) 2
and denote by

R n (θ, f ) = 1 n n i=1 Y i -f (θ T X i ) 2
the empirical counterpart of R(θ, f ) based on the sample D n .

Estimation procedure

We are now in a position to describe our estimation procedure. The method which is presented here is inspired by the approach developed by Catoni in [START_REF] Catoni | Statistical Learning Theory and Stochastic Optimization[END_REF][START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF]. It strongly relies on the choice of a probability measure π on S p 1,+ ×F, called the prior, which in our framework should enforce the sparsity properties of the target regression function. With this objective in mind, we first let dπ(θ, f ) = dµ(θ)dν(f ),

i.e., we assume that the distribution over the indexes is independent of the distribution over the link functions. With respect to the parameter θ, we put

dµ(θ) = p i=1 10 -i I⊂{1,...,p},|I|=i p i -1 dµ I (θ) 1 -( 1 10 ) p , (2.1)
where |I| denotes the cardinality of I and dµ I (θ) is the uniform probability measure on the set

S p 1,+ (I) = {θ = (θ 1 , . . . , θ p ) ∈ S p 1,+ : θ j = 0 if and only if j / ∈ I}.
We see that S p 1,+ (I) may be interpreted as the set of "active" coordinates in the single-index regression of Y on X, and note that the prior on S p 1,+ is a convex combination of uniform probability measures on the subsets S p 1,+ (I). The weights of this combination depend only on the size of the active coordinate subset I. As such, the value |I| characterizes the sparsity of the model: The smaller |I|, the smaller the number of variables involved in the model. The factor 10 -i penalizes models of high dimension, in accordance with the sparsity idea.

The choice of the prior ν on F is more involved. To begin with, we define, for any positive integer M ≤ n and all Λ > 0,

B M (Λ) = (β 1 , . . . , β M ) ∈ R M : M j=1 j|β j | ≤ Λ and β M = 0 .
Next, we let F M (Λ) ⊂ F be the image of B M (Λ) by the map

Φ M : R M → F (β 1 , . . . , β M ) → M j=1 β j ϕ j .
It is worth pointing out that, roughly, Sobolev spaces are well approximated by F M (Λ) as M grows (more on this in Subsection 2.3). Finally, we define ν M (df ) on the set F M (C +1) as the image of the uniform measure on B M (C + 1) induced by the map Φ M , and take

dν(f ) = n M =1 10 -M dν M (f ) 1 -( 1 10 ) n . (2.2)
Some comments are in order here. First, we note that the prior π is defined on S p 1,+ × F n (C + 1) endowed with its canonical Borel σ-field. The choice of C + 1 instead of C in the definition of the prior support is essentially technical. This bound ensures that when the target f ⋆ belongs to F n (C), then a small ball around it is contained in F n (C + 1). It could be safely replaced by C + u n , where {u n } ∞ n=1 is any positive sequence vanishing sufficiently slowly as n → ∞. Next, the integer M should be interpreted as a measure of the "dimension" of the function f -the larger M , the more complex the function-and the prior ν adapts again to the sparsity idea by penalizing large-dimensional functions f . The coefficients 10 -i and 10 -M which appear in (2.1) and (2.2) show that more complex models have a geometrically decreasing influence. Note however that the value 10, which has been chosen because of its good practical results, is somehow arbitrary. It could be, in all generality, replaced by a more general coefficient α at the price of a more technical analysis. Finally, we observe that, for each f = M j=1 β j ϕ j ∈ F M (C + 1),

f ∞ ≤ M j=1 |β j | ≤ C + 1.
Now, let λ be a positive real number, called the inverse temperature parameter hereafter. The estimates θλ and fλ of θ ⋆ and f ⋆ , respectively, are simply obtained by randomly drawing

( θλ , fλ ) ∼ ρλ ,
where ρλ is the so-called Gibbs posterior distribution over S p 1,+ × F n (C + 1), defined by the probability density

dρ λ dπ (θ, f ) = exp [-λR n (θ, f )] exp [-λR n (θ, f )] dπ(θ, f )
.

[The notation dρ λ /dπ means the density of ρλ with respect to π.] The estimate ( θλ , fλ ) has a simple interpretation. Firstly, the level of significance of each pair (θ, f ) is assessed via its least square error performance on the data D n . Secondly, a Gibbs distribution with respect to the prior π enforcing those pairs (θ, f ) with the most empirical significance is assigned on the space S p 1,+ × F n (C + 1). Finally, the resulting estimate is just a random realization (conditional to the data) of this Gibbs posterior distribution.

Sparsity oracle inequality

For any I ⊂ {1, . . . , p} and any positive integer M ≤ n, we set

θ ⋆ I,M , f ⋆ I,M ∈ arg min (θ,f )∈S p 1,+ (I)×F M (C) R(θ, f ).
At this stage, it is very important to note that, for each M , the infimum f ⋆ I,M is defined on F M (C), whereas the prior charges a slightly bigger set, namely

F M (C + 1).
The main result of the paper is the following theorem. Here and everywhere, the wording "with probability 1δ" means the probability evaluated with respect to the distribution P ⊗n of the data D n and the conditional probability measure ρλ . Recall that ℓ is a positive constant such that ϕ ′ j ∞ ≤ ℓj.

Theorem 2.1 Assume that Assumption N and Assumption B hold. Set w = 8(2C + 1) max[L, 2C + 1] and take

λ = n w + 2 [(2C + 1) 2 + 4σ 2 ] . (2.3)
Then, for all δ ∈ ]0, 1[, with probability at least 1δ we have

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n R(θ ⋆ I,M , f ⋆ I,M ) -R(θ ⋆ , f ⋆ ) + M log(Cn) + |I| log(pn) + log 2 δ n ,
where Ξ is a positive constant, function of L, C, σ and ℓ only.

Remark 2.2 Interestingly enough, analysis of the estimate ( θλ , fλ ) is still possible when Assumption N is not satisfied. Indeed, even if Bernstein's inequality (see Lemma 4.1) is not valid, a recent paper by Seldin, Cesa-Bianchi, Laviolette, Auer, Shawe-Taylor and Peters [START_REF] Seldin | PAC-Bayesian analysis of the exploration-exploitation trade-off[END_REF] provides us with a nice alternative inequality assuming less restrictive assumptions. However, we would then suffer a loss in the upper bound of Theorem 2.1. It is also interesting to note that recent results by Audibert and Catoni [START_REF] Audibert | Robust linear least squares regression[END_REF] allow the study of PAC-Bayesian estimates without Assumption N. However, the results of these authors are valid for linear models only, and it is therefore not clear to what extent their technique can be transposed to our setting.

Theorem 2.1 can be given a simple interpretation. Indeed, we see that if there is a "small" I and a "small"

M such that R(θ ⋆ I,M , f ⋆ I,M ) is close to R(θ ⋆ , f ⋆ ), then R( θλ , fλ ) is also close to R(θ ⋆ , f ⋆ )
up to terms of order 1/n. However, if no such I or M exists, then one of the terms M log(Cn)/n and |I| log(pn)/n starts to dominate, thereby deteriorating the general quality of the bound. A good approximation with a "small" I is typically possible when θ ⋆ is sparse or, at least, when it can be approximated by a sparse parameter. On the other hand, a good approximation with a "small" M is possible if f ⋆ has a sufficient degree of regularity.

To illustrate the latter remark, assume for instance that {ϕ j } ∞ j=1 is the (nonnormalized) trigonometric system and suppose that the target f ⋆ belongs to the Sobolev ellipsoid, defined by

W k, 6C 2 π 2 = f ∈ L 2 ([-1, 1]) : f = ∞ j=1 β j ϕ j and ∞ j=1 j 2k β 2 j ≤ 6C 2 π 2
for some unknown regularity parameter k ≥ 2 (see, e.g., Tsybakov [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]).

Observe that, in this context, the approximation sets F M (C + 1) take the form

F M (C + 1) = f ∈ L 2 ([-1, 1]) : f = M j=1 β j ϕ j , M j=1 j|β j | ≤ C + 1 and β M = 0 .
It is important to note that the regularity parameter k is assumed to be unknown, and this casts our results in the so-called adaptive setting. The following additional assumption will be needed:

Assumption D. The random variable θ ⋆T X has a probability density on [-1, 1], bounded from above by a positive constant B.

Last, we let I ⋆ be the set I such that θ ⋆ ∈ S p 1,+ (I) and set

θ ⋆ 0 = |I ⋆ |.
Corollary 2.1 Assume that Assumption N, Assumption B and Assumption D hold. Suppose also that f ⋆ belongs to the Sobolev ellipsoid

W(k, 6C 2 /π 2 ),
where the real number k ≥ 2 is an (unknown) regularity parameter. Set w = 8(2C + 1) max[L, 2C + 1] and take λ as in (2.3). Then, for all δ ∈ ]0, 1[, with probability at least 1δ we have

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ ′ log(Cn) n 2k 2k+1 + θ ⋆ 0 log(pn) n + log 2 δ n , (2.4)
where Ξ ′ is a positive constant, function of L, C, σ, ℓ and B only.

As far as we are aware, all existing methods achieving rates of convergence similar to the ones provided by Corollary 2.1 are valid in an asymptotic setting only (p fixed and n → ∞). The strength of Corollary 2.1 is to provide a finite sample bound and to show that our estimate still behaves well in a nonasymptotic situation if the intrinsic dimension (i.e., the sparsity) is small with respect to n. To understand this remark, just assume that p is a function of n such that p → ∞ as n → ∞. Whereas a classical asymptotic approach cannot say anything clever about this situation, our bounds still provide some information, provided the model is sparse enough (i.e., θ ⋆ 0 is sufficiently small with respect to n).

We see that, asymptotically (p fixed and n → ∞), the leading term on the right-hand side of inequality (2.4) is (log

(n)/n) 2k 2k+1
. This is the minimax rate of convergence over a Sobolev class, up to a log(n) factor. However, when n is "small" and θ ⋆ is not sparse (i.e., θ ⋆ 0 is not "small"), the term θ ⋆ 0 log(pn)/n starts to emerge and cannot be neglected. Put differently, in large dimension, the estimation of θ ⋆ itself is a problem-this phenomenon is not taken into account by asymptotic studies.

It is worth mentioning that the approach developed in the present article does not offer any guarantee on the point of view of variable (feature) selection. To reach this objective, an interesting route to follow is the sufficient dimension reduction (SDR) method proposed by Chen, Zou and Cook [START_REF] Chen | Coordinate-independent sparse sufficient dimension reduction and variable selection[END_REF], which can be applied to the single-index model to estimate consistently the parameter θ ⋆ and perform variable selection in a sparsity framework. Note however that such results require strong assumptions on the distribution of the data.

Finally, it should be stressed that the choice of λ in Theorem 2.1 and Corollary 2.1 is not the best possible and may eventually be improved, at the price of a more technical analysis however.

Implementation and numerical results

A series of experiments was conducted, both on simulated and real-life data sets, in order to assess the practical capabilities of the proposed method and compare its performance with that of standard procedures. Prior to analysis, we first need to discuss its concrete implementation, which has been carried out via a Markov Chain Monte Carlo (MCMC) method.

Implementation via reversible jump MCMC

The use of MCMC methods has become a popular way to compute Bayesian estimates. For an introduction to the domain, one should refer to the comprehensive monograph of Marin and Robert [START_REF] Marin | Bayesian Core: A Practical Approach to Computational Bayesian Analysis[END_REF] and the references therein. Importantly, in this computational framework, an adaptation of the well-known Hastings-Metropolis algorithm to the case where the posterior distribution gives mass to several models of different dimensions was proposed by Green [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] under the name Reversible Jump MCMC (RJMCMC) method. In the PAC-Bayesian setting, MCMC procedures were first considered by Catoni [START_REF] Catoni | Statistical Learning Theory and Stochastic Optimization[END_REF], whereas Dalalyan and Tsybakov [START_REF] Dalalyan | Aggregation by exponential weighting, sharp PAC-Bayesian bounds and sparsity[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF] and Alquier and Lounici [START_REF] Alquier | PAC-Bayesian bounds for sparse regression estimation with exponential weights[END_REF] explore their practical implementation in the sparse context using Langevin Monte Carlo and RJMCMC, respectively. Regarding the single-index model, MCMC algorithms were used to compute Bayesian estimates by Antoniadis, Grégoire and McKeague [START_REF] Antoniadis | Bayesian estimation in single-index models[END_REF] and, more recently, by Wang [START_REF] Wang | Bayesian estimation and variable selection for single index models[END_REF], who develop a fully Bayesian method to analyse the single-index model. Our implementation technique is close in spirit to the one of Wang [START_REF] Wang | Bayesian estimation and variable selection for single index models[END_REF].

As a starting point for the approximate computation of our estimate, we used the RJMCMC method of Green [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF], which is in fact an adaptation of the Hastings-Metropolis algorithm to the case where the objective posterior probability distribution (here, ρλ ) assigns mass to several different models. The idea is to start from an initial given pair (θ (0) , f (0) ) ∈ S p 1,+ × F n (C + 1) and then, at each step, to iteratively compute (θ (t+1) , f (t+1) ) from (θ (t) , f (t) ) via the following chain of rules:

• Sample a random pair (τ (t) , h (t) ) according to some proposal conditional density k t ( . |(θ (t) , f (t) )) with respect to the prior π;

• Take (θ (t+1) , f (t+1) ) = (τ (t) , h (t) ) with probability α t (θ (t) , f (t) ) with probability 1 -α t ,
where

α t = min 1, dρ λ dπ (τ (t) , h (t) ) × k t (θ (t) , f (t) )|(τ (t) , h (t) ) dρ λ dπ (θ (t) , f (t) ) × k t ((τ (t) , h (t) )|(θ (t) , f (t) ))
.

This protocol ensures that the sequence {(θ (t) , f (t) )} ∞ t=0 is a Markov chain with invariant probability distribution ρλ (see e.g. Marin and Robert [START_REF] Marin | Bayesian Core: A Practical Approach to Computational Bayesian Analysis[END_REF]). A usual choice is to take k t ≡ k, so that the Markov chain is homogeneous. However, in our context, it is more convenient to let k t = k 1 if t is odd and k t = k 2 if t is even. Roughly, the effect of k 1 is to modify the index θ (t) while k 2 will essentially act on the link function f (t) . While the ideas underlying the proposal densities k 1 and k 2 are quite simple, a precise description in its full length turns out to be more technical. Thus, in order to preserve the readability of the paper, the explicit construction of k 1 and k 2 has been postponed to the Appendix Section 5.

From a theoretical point of view, it is clear that the implementation of our method requires knowledge of the constant C (the upper bound on f ⋆ ∞ ). A too small C will result in a smaller model, which is unable to perform a good approximation. On the other hand, a larger C induces a poor bound in Theorem 2.1. In practice, however, the influence of C turns out to be secondary compared to the impact of the parameter λ. Indeed, it was found empirically that a very large choice of C (e.g., C = 10 100 ) does not deteriorate the overall quality of the results, as soon as λ is appropriately chosen. This is the approach that was followed in the experimental testing process.

Besides, the time for the Markov chains to converge depends strongly on the ambient dimension p and the starting point of the simulations. When the dimension is small (typically, p ≤ 10), the chains converge fast and any value may be chosen as a starting point. In this case, we let the MCMC run 1000 steps and obtained satisfying results. On the other hand, when the dimension is larger (typically, p > 10), the convergence is very slow, in the sense that R n (θ (t) , f (t) ) takes a very long time to stabilize. However, using as a starting point for the chains the preliminary estimate θHHI (see below) significantly reduces the number of steps needed to reach convergence-we let the chains run 5000 steps in this context. Nevertheless, as a general rule, we encourage the users to inspect the convergence of the chains by checking if R n (θ (t) , f (t) ) is stabilized, and to run several chains starting from different points to avoid their attraction into local minima.

Simulation study

In this subsection, we illustrate the finite sample performance of the presented estimation method on three synthetic data sets and compare its predictive capabilities with those of three standard statistical procedures. In all our experiments, we took as dictionary the (non-normalized) trigonometric system {ϕ j } ∞ j=1 and denote accordingly the resulting regression function estimate defined in Section 2 by FFourier . In accordance with the order of magnitude indicated by the theoretical results, we set λ = 4n. This choice can undoubtedly be improved a bit but, as the numerical results show, it seems sufficient for our procedure to be fairly competitive.

The tested competing methods are the Lasso (Tibshirani [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]), the standard regression kernel estimate (Nadaraya [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Nadaraya | Remarks on nonparametric estimates for density functions and regression curves[END_REF] and Watson [START_REF] Watson | Smooth regression analysis[END_REF], see also Tsybakov [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]), and the estimation strategy discussed in Härdle, Hall and Ichimura [START_REF] Härdle | Optimal smoothing in single-index models[END_REF]. While the procedure of Härdle, Hall and Ichimura is specifically tailored for single-index models, the Lasso is designed to deal with the estimation of sparse linear models. On the other hand, the nonparametric kernel method is one of the best options when no obvious assumption (such as the single-index one) can be made on the shape of the targeted regression function.

We briefly recall that, for a linear model of the form Y = θ ⋆T X + W , the Lasso estimate takes the form FLasso (x) = θT Lasso x, where θLasso ∈ arg min

θ∈R p 1 n n i=1 Y i -θ T X i 2 + ξ p j=1 |θ j |
and ξ > 0 is a regularization parameter. Theoretical results (see e.g. Bunea, Tsybakov and Wegkamp [START_REF] Bunea | Sparsity oracle inequalities for the Lasso[END_REF]) indicate that ξ should be of the order ξ ⋆ = σ log(p)/n. Throughout, σ is assumed to be known, and we let ξ = ξ ⋆ /3, since this choice is known to give good practical results. The Nadaraya-Watson kernel estimate will be denoted by FNW . It is defined by

FNW (x) = n i=1 Y i K h (x -X i ) n i=1 K h (x -X i )
for some nonnegative kernel K on R p and K h (z) = K(z/h)/h. In the experiments, we let K be the Gaussian kernel K(z) = exp(-z T z) and chose the smoothing parameter h via a classical leave-one-out procedure on the grid G = {0.75 k , k = 0, . . . , ⌊log(n)⌋}, see, e.g., Györfi, Kohler, Krzyżak and Walk [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF] (notation ⌊.⌋ stands for the floor function). Finally, the estimation procedure advocated in Härdle, Hall and Ichimura [START_REF] Härdle | Optimal smoothing in single-index models[END_REF] takes the form

FHHI (x) = n i=1 Y i G ĥ θT HHI (x -X i ) n i=1 G ĥ θT HHI (x -X i )
for some kernel G on R, with G h (z) = G(z/h)/h and ĥ, θHHI ∈ arg min

h>0,θ∈R p n i=1 Y i - j =i Y j G h θ T (X j -X i ) j =i G h (θ T (X j -X i )) 2 .
All calculations were performed with the Gaussian kernel. We used the grid G for the optimization with respect to h, whereas the best search for θ was implemented via a pathwise coordinate optimization.

The various methods were tested for the general regression model

Y i = F (X i ) + W i , i = 1, . . . , n,
for three different choices of F (single-index or not) and two values of n, namely n = 50 and n = 100. In each of these models, the observations X i take values in R p , with p = 10 and p = 50, and have independent components uniformly distributed on [-1, 1]. The noise variables W 1 , . . . , W n are independently distributed according to a Gaussian N (0, σ 2 ), with σ = 0.2. It is worth pointing out that for n = 50 and p = 50, p and n are of the same order, which means that the setting is nonasymptotic. It is essentially in this case that the use of estimates tailored to sparsity, which reduce the variance, is expected to improve the performance over generalist methods.

On the other hand, the situation n = 100 and p = 10 is less difficult and mimics the asymptotic setting.

The three examined functions F (x), for x = (x 1 , . . . , x p ), were the following ones:

[Model 1] A linear model F Linear (x) = 2θ ⋆T x. [Model 2] A single-index function F SI (x) = 2(θ ⋆T x) 2 + θ ⋆T x. [Model 3] A purely nonparametric model F NP (x) = 2|x 2 | |x 1 | -x 3 3 ,
where, in the first and second model, θ ⋆ = (0.5, 0.5, 0, . . . , 0) For each experiment, a learning set of size n was generated to compute the estimates and their performance, in terms of mean square prevision error, was evaluated on a separate test set of the same size. The results are shown in Table 1 (p = 10) and Table 2 (p = 50). As each experiment was repeated 20 times, these tables report the median, the mean and the standard deviation (s.d.) of the prevision error of each procedure.

Some comments are in order. First, we note without surprise that:

1. The Lasso performs well in the linear setting [Model 1]. Interestingly, FFourier provides slightly better results than the single-indextailored estimate FHHI , especially for p = 50. This observation can be easily explained by the fact that FHHI does not integrate any sparsity information regarding the parameter θ ⋆ , whereas FFourier tries to focus on the dimension of the active coordinates, which is equal to 2 in this simulation. As a general finding, we retain that FFourier is the most robust of all the tested procedures.

Real data

The real-life data sets used in this second series of experiments are from two different sources. The first one, called AIR-QUALITY data (n = 111, p = 3), has been first used by Chambers, Cleveland, Kleiner and Tukey [START_REF] Chambers | Graphical Methods for Data Analysis[END_REF] and has been later considered as a benchmark in the study and comparison of single-index models (see, for example, Antoniadis, Grégoire and McKeague [START_REF] Antoniadis | Bayesian estimation in single-index models[END_REF] and Wang [START_REF] Wang | Bayesian estimation and variable selection for single index models[END_REF], among others). This data set originated from an environmental study relating n = 111 ozone concentration measures at p = 3 meteorological variables, namely wind speed, temperature and radiation. The data is available as a package in the software R [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF], which we employed in all the numerical experiments. The programs are available upon request from the authors.

The second category of data arises from the UC Irvine Machine Learning Repository http://archive.ics.uci.edu/ml, where the following packages have been downloaded from:

• AUTO-MPG (Quinlan [START_REF] Quinlan | Combining instance-based and model-based learning[END_REF], n = 392, p = 7).

• CONCRETE (Yeh [START_REF] Yeh | Modeling of strength of high-performance concrete using artificial neural networks[END_REF], n = 1030, p = 8).

• HOUSING (Harrison and Rubinfeld [START_REF] Harrison | Hedonic housing prices and the demand for clean air[END_REF], n = 508, p = 13).

• SLUMP-1, SLUMP-2 and SLUMP-3, which correspond to the concrete slump test data introduced by Yeh [START_REF] Yeh | Modeling slump flow of concrete using second-order regressions and artificial neural networks[END_REF] (n = 51, p = 7). Since there are 3 different output variables Y in the original data set, we created a single experiment for each of these variables (1 refers to the output "slump", 2 to the output "flow" and 3 to the output "28-day Compressive Strength").

• WINE-RED and WINE-WHITE (Cortez, Cerdeira, Almeida, Matos and Reis [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF], n = 1599, n = 4898, p = 11).

We refer to the above-mentioned references for a precise description of the meaning of the variables involved in these data sets. For homogeneity reasons, all data were normalized to force the input variables to lie in [-1, 1]-in accordance with the setting of our method-and to ensure that all output variables have standard deviation 0.5. In two data sets (AIR-QUALITY and AUTO-MPG) there were some missing values and the corresponding observations were simply removed.

For each method and each of the nine data sets, we randomly split the observations in a learning and a test set of equal sizes, computed the estimate on the learning set, evaluated the prediction error on the test set, and repeated this protocol 20 times. The results are summarized in Table 3.

We see that all the tested methods provide reasonable results on most data sets. The Lasso is very competitive, especially in the nonasymptotic framework. The estimation procedure FFourier offers outcomes which are similar to the ones of FHHI , with a slight advantage for the latter method however. Altogether, FFourier and FHHI provide the best performance in terms of prediction error in 6 out of 9 experiments. Besides, when it is not the best, the method FFourier is close to the best one, as for example in SLUMP-3 and WINE-RED. As an illustrative example, the plot of the resulting fit of our procedure to the data set AUTO-MPG is shown in Figure 1.

Clearly, all data sets under study have a dimension p which is small compared to n. To correct this situation, we ran the same series of experiments by adding some additional irrelevant dimensions to the data. Specifically, the observations were embedded into a space of dimension p × 4 by letting the new fake coordinates follow independent uniform [0, 1] random variables. The results are shown in Table 4. In this nonasymptotic framework, the method FHHI -which is not designed for sparsity-collapses, whereas FFourier takes a clear advantage over its competitors. In fact, it provides the best results in 3 out of 9 experiments (AUTO-MPG, CONCRETE and HOUSING).

Besides, when it is not the best, the method FFourier is very close to the best one, as for example in SLUMP-3 and WINE-RED.

Thus, as a general conclusion to this experimental section, we may say that our PAC-Bayesian oriented procedure has an excellent predictive ability, even in nonasymptotic/high-dimensional situations. It is fast, robust, and exhibits performance at the level of the gold standard Lasso.

Proofs

Preliminary results

Throughout this section, we let π be the prior probability measure on R p × F n (C + 1) equipped with its canonical Borel σ-field. Recall that F n (C + 1) ⊂ F and that, for each f ∈ F n (C + 1), we have

f ∞ ≤ C + 1. Besides, since E[Y |X] = f ⋆ (θ ⋆T X)
almost surely, we note once and for all that for all (θ, f

) ∈ S p 1,+ × F n (C + 1), R(θ, f ) -R(θ ⋆ , f ⋆ ) = E Y -f (θ T X) 2 -E Y -f ⋆ (θ ⋆T X) 2 = E f (θ T X) -f ⋆ (θ ⋆T X)
E (T i ) k + ≤ k! 2 vw k-2 .
Then, for any ζ ∈ ]0, 1/w[,

E exp ζ n i=1 [T i -ET i ] ≤ exp vζ 2 2(1 -wζ) .
Given a measurable space (E, E) and two probability measures µ 1 and µ 2 on (E, E), we denote by K(µ 1 , µ 2 ) the Kullback-Leibler divergence of µ 1 with respect to µ 2 , defined by

K(µ 1 , µ 2 ) =    log dµ 1 dµ 2 dµ 1 if µ 1 ≪ µ 2 , ∞ otherwise.
(Notation µ 1 ≪ µ 2 means "µ 1 is absolutely continuous with respect to µ 2 ".) Lemma 4.2 Let (E, E) be a measurable space. For any probability measure µ on (E, E) and any measurable function h :

E → R such that (exp •h)dµ < ∞, we have log (exp •h)dµ = sup m hdm -K(m, µ) , (4.1)
where the supremum is taken over all probability measures on (E, E) and, by convention, ∞ -∞ = -∞. Moreover, as soon as h is bounded from above on the support of µ, the supremum with respect to m on the right-hand side of (4.

λ ∈ 0, n w + [(2C + 1) 2 + 4σ 2 ]
.

Then, for all δ ∈ ]0, 1[ and any data-dependent probability measure ρ absolutely continuous with respect to π we have, with probability at least

1 -δ, R( θ, f ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ R n ( θ, f ) -R n (θ ⋆ , f ⋆ ) + log dρ dπ ( θ, f ) + log 1 δ λ ,
where the pair ( θ, f ) is distributed according to ρ.

Proof of Lemma 4.3. Fix θ ∈ S p 1,+ and f ∈ F n (C + 1). The proof starts with an application of Lemma 4.1 to the random variables

T i = -Y i -f (θ T X i ) 2 + Y i -f ⋆ (θ ⋆T X i ) 2 , i = 1, . . . , n.
Note that these random variables are independent, identically distributed, and that n i=1

ET 2 i = n i=1 E 2Y i -f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 = n i=1 E 2W i + f ⋆ (θ ⋆T X i ) -f (θ T X i ) 2 f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 ≤ n i=1 E 4W 2 i + (2C + 1) 2 f (θ T X i ) -f ⋆ (θ ⋆T X i )
Thus, by Assumption N,

n i=1 ET 2 i ≤ (2C + 1) 2 + 4σ 2 n i=1 E f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 ≤ v,
where we set

v = 2n[(2C + 1) 2 + 4σ 2 ] [R(θ, f ) -R(θ ⋆ , f ⋆ )] .
More generally, for all integers k ≥ 3,

n i=1 E (T i ) k + ≤ n i=1 E 2Y i -f (θ T X i ) -f ⋆ (θ ⋆T X i ) k f (θ T X i ) -f ⋆ (θ ⋆T X i ) k = n i=1 E 2W i + f ⋆ (θ ⋆T X i ) -f (θ T X i ) k f (θ T X i ) -f ⋆ (θ ⋆T X i ) k ≤ 2 k-1 n i=1 E 2 k |W i | k + (2C + 1) k (2C + 1) k-2 f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 .
In the last inequality, we used the fact that

|a+b| k ≤ 2 k-1 (|a| k +|b| k ) together with f (θ T X i ) -f ⋆ (θ ⋆T X i ) k = f (θ T X i ) -f ⋆ (θ ⋆T X i ) k-2 × f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 ≤ (2C + 1) k-2 f (θ T X i ) -f ⋆ (θ ⋆T X i ) 2 .
Therefore, by Assumption N,

n i=1 E (T i ) k + ≤ n i=1 2 2k-2 k!σ 2 L k-2 + 2 k-1 (2C + 1) k (2C + 1) k-2 [R(θ, f ) -R(θ ⋆ , f ⋆ )] = v × 2 2k-2 k!σ 2 L k-2 + 2 k-1 (2C + 1) k (2C + 1) k-2 [(2C + 1) 2 + 4σ 2 ] ≤ v × 8 k-2 k! max L k-2 , (2C + 1) k-2 (2C + 1) k-2 2 = k! 2 vw k-2 , with w = 8(2C + 1) max[L, 2C + 1].
Thus, for any inverse temperature parameter λ ∈ ]0, n/w[, taking ζ = λ/n, we may write by Lemma 4.1

E exp [λ (R(θ, f ) -R(θ ⋆ , f ⋆ ) -R n (θ, f ) + R n (θ ⋆ , f ⋆ ))] ≤ exp vλ 2 2n 2 (1 -wλ n )
.

Therefore, using the definition of v, we obtain

E exp λ - λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ, f ) -R(θ ⋆ , f ⋆ )) + λ (-R n (θ, f ) + R n (θ ⋆ , f ⋆ )) -log 1 δ ≤ δ.
Next, we use a standard PAC-Bayesian approach (Catoni [START_REF] Catoni | Statistical Learning Theory and Stochastic Optimization[END_REF][START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF], Audibert [START_REF] Audibert | Aggregated estimators and empirical complexity for least square regression[END_REF] and Alquier [START_REF] Alquier | PAC-Bayesian bounds for randomized empirical risk minimizers[END_REF]). Let us remind the reader that π is a prior probability measure on the set S p 1,+ × F n (C + 1). We have

E exp λ - λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ, f ) -R(θ ⋆ , f ⋆ )) + λ (-R n (θ, f ) + R n (θ ⋆ , f ⋆ )) -log 1 δ dπ(θ, f ) ≤ δ
and consequently, using Fubini's theorem,

E exp λ - λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ, f ) -R(θ ⋆ , f ⋆ )) + λ (-R n (θ, f ) + R n (θ ⋆ , f ⋆ )) -log 1 δ dπ(θ, f ) ≤ δ.
Therefore, for any data-dependent posterior probability measure ρ absolutely continuous with respect to π, adopting the convention ∞ × 0 = 0,

E exp λ - λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ, f ) -R(θ ⋆ , f ⋆ )) + λ (-R n (θ, f ) + R n (θ ⋆ , f ⋆ )) -log dρ dπ (θ, f ) -log 1 δ dρ(θ, f ) ≤ δ.
Recalling that P ⊗n stands for the distribution of the sample D n , the latter inequality can be more conveniently written as

E Dn∼P ⊗n E ( θ, f )∼ρ exp λ - λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) R( θ, f ) -R(θ ⋆ , f ⋆ ) + λ -R n ( θ, f ) + R n (θ ⋆ , f ⋆ ) -log dρ dπ ( θ, f ) -log 1 δ ≤ δ.
Thus, using the elementary inequality exp(λx) ≥ 1 R + (x) we obtain, with probability at most δ,

1 - λ [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) R( θ, f ) -R(θ ⋆ , f ⋆ ) ≥ R n ( θ, f ) -R n (θ ⋆ , f ⋆ ) + log dρ dπ ( θ, f ) + log 1 δ λ ,
where the probability is evaluated with respect to the distribution P ⊗n of the data D n and the conditional probability measure ρ. Put differently, letting

λ ∈ 0, n w + [(2C + 1) 2 + 4σ 2 ] ,
we have, with probability at least 1δ,

R( θ, f ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ R n ( θ, f ) -R n (θ ⋆ , f ⋆ ) + log dρ dπ ( θ, f ) + log 1 δ λ .
This concludes the proof of Lemma 4.3.

Lemma 4.4 Under the conditions of Lemma 4.3 we have, with probability at least 1δ,

R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) ≤ 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ(θ, f ) -R(θ ⋆ , f ⋆ ) + K (ρ, π) + log 1 δ λ .
Proof of Lemma 4.4. The beginning of the proof is similar to the one of Lemma 4.3. More precisely, we apply Lemma 4.1 with

T i = (Y i -f (θ T X i )) 2 - (Y i -f ⋆ (θ ⋆T X i ))
2 and obtain, for any inverse temperature parameter λ ∈ ]0, n/w[,

E exp [λ (R(θ ⋆ , f ⋆ ) -R(θ, f ) -R n (θ ⋆ , f ⋆ ) + R n (θ, f ))] ≤ exp vλ 2 2n 2 (1 -wλ n )
.

Thus, using the definition of v,

E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ ⋆ , f ⋆ ) -R(θ, f )) + λ (R n (θ, f ) -R n (θ ⋆ , f ⋆ )) -log 1 δ ≤ δ.
Integrating with respect to π leads to

E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ ⋆ , f ⋆ ) -R(θ, f )) + λ (R n (θ, f ) -R n (θ ⋆ , f ⋆ )) -log 1 δ dπ(θ, f ) ≤ δ whence, by Fubini's theorem, E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ ⋆ , f ⋆ ) -R(θ, f )) + λ (R n (θ, f ) -R n (θ ⋆ , f ⋆ )) -log 1 δ dπ(θ, f ) ≤ δ.
Thus, for any data-dependent posterior probability measure ρ absolutely continuous with respect to π,

E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ ⋆ , f ⋆ ) -R(θ, f )) + λ (R n (θ, f ) -R n (θ ⋆ , f ⋆ )) -log dρ dπ (θ, f ) -log 1 δ dρ(θ, f ) ≤ δ.
Therefore, by Jensen's inequality,

E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) (R(θ ⋆ , f ⋆ ) -R(θ, f )) + λ (R n (θ, f ) -R n (θ ⋆ , f ⋆ )) -log dρ dπ (θ, f ) -log 1 δ dρ(θ, f ) = E exp λ + λ 2 [(2C + 1) 2 + 4σ 2 ] n(1 -wλ n ) R(θ ⋆ , f ⋆ ) -R(θ, f )dρ(θ, f ) + λ R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) -K (ρ, π) -log 1 δ ≤ δ.
Consequently, by the elementary inequality exp(λx) ≥ 1 R + (x), we obtain, with probability at most δ,

R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) ≥ 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ(θ, f ) -R(θ ⋆ , f ⋆ ) + K (ρ, π) + log 1 δ λ .
Equivalently, with probability at least 1δ,

R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) ≤ 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ(θ, f ) -R(θ ⋆ , f ⋆ ) + K (ρ, π) + log 1 δ λ .

Proof of Theorem 2.1

The proof starts with an application of Lemma 4.3 with ρ = ρλ (the Gibbs distribution) as posterior distribution. More precisely, we know that, with probability larger than 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ R n ( θλ , fλ ) -R n (θ ⋆ , f ⋆ ) + log dρ λ dπ ( θλ , fλ ) + log 1 δ λ ,
where the probability is evaluated with respect to the distribution P ⊗n of the data D n and the conditional probability measure ρλ . Observe that

log dρ λ dπ ( θλ , fλ ) = log     exp -λR n ( θλ , fλ ) exp [-λR n (θ, f )] dπ(θ, f )     = -λR n ( θλ , fλ ) -log exp [-λR n (θ, f )] dπ(θ, f ).
Consequently, with probability at least 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 λ 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ -log exp [-λR n (θ, f )] dπ(θ, f ) -λR n (θ ⋆ , f ⋆ ) + log 1 δ .
Next, using Lemma 4.2 we deduce that, with probability at least 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf ρ R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) + K(ρ, π) + log 1 δ λ ,
where the infimum is taken over all probability measures on S p 1,+ ×F n (C +1). In particular, letting M(I, M ) be the set of all probability measures on S p 1,+ (I) × F M (C + 1), we have, with probability at least 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n inf ρ∈M(I,M ) R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) + K(ρ, π) + log 1 δ λ .
Next, observe that, for ρ ∈ M(I, M ),

K(ρ, π) = K(ρ, µ ⊗ ν) = K(ρ, µ I ⊗ ν M ) + log 1 -1 10 p 1 -1 10 n p |I| 10 -|I|-M ≤ K(ρ, µ I ⊗ ν M ) + log p |I| 10 -|I|-M . (4.2)
Therefore, with probability at least 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n inf ρ∈M(I,M ) R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) + K(ρ, µ I ⊗ ν M ) + log ( p |I| ) 10 -|I|-M + log 1 δ λ . (4.3)
By Lemma 4.4 and inequality (4.2), for any data-dependent distribution ρ ∈ M(I, M ), with probability at least 1δ,

R n (θ, f )dρ(θ, f ) -R n (θ ⋆ , f ⋆ ) ≤ 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ(θ, f ) -R(θ ⋆ , f ⋆ ) + K(ρ, µ I ⊗ ν M ) + log ( p |I| ) 10 -|I|-M + log 1 δ λ . (4.4)
Thus, combining inequalities (4.3) and (4.4), we may write, with probability at least 1 -2δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n inf ρ∈M(I,M ) 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ(θ, f ) -R(θ ⋆ , f ⋆ ) + 2 K(ρ, µ I ⊗ ν M ) + log ( p |I| ) 10 -|I|-M + log 1 δ λ . (4.5)
For any subset I of {1, . . . , p}, any positive integer M ≤ n and any η, γ ∈ ]0, 1/n], let the probability measure ρ I,M,η,γ be defined by

dρ I,M,η,γ (θ, f ) = dρ 1 I,M,η (θ)dρ 2 I,M,γ (f ), with dρ 1 I,M,η dµ I (θ) ∝ 1 [ θ-θ ⋆ I,M 1 ≤η] and dρ 2 I,M,γ dν M (f ) ∝ 1 [ f -f ⋆ I,M M ≤γ]
where, for f = M j=1 β j ϕ j ∈ F M (C + 1), we put

f M = M j=1 j|β j |.
With this notation, inequality (4.5) leads to Similarly, by technical Lemma 4.6,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n inf η,γ>0 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ, f )dρ I,M,η,γ (θ, f ) -R(θ ⋆ , f ⋆ ) + 2 K(ρ I,M,η,γ , µ I ⊗ ν M ) + log ( p |I| ) 10 -|I|-M + log 1 δ λ . ( 4
K(ρ 2 I,M,γ , ν M ) = M log C + 1 γ .
Putting all the pieces together, we are led to

K(ρ I,M,η,γ , µ I ⊗ ν M ) ≤ (|I| -1) log max |I|, 4 η + M log C + 1 γ . (4.8)
Finally, it remains to control the term R(θ, f )dρ I,M,η,γ (θ, f ).

To this aim, we write

R(θ, f )dρ I,M,η,γ (θ, f ) = E Y -f (θ T X) 2 dρ I,M,η,γ (θ, f ) = E Y -f ⋆ I,M (θ ⋆T I,M X) + f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) + f (θ ⋆T I,M X) -f (θ T X) 2 dρ I,M,η,γ (θ, f ) = R(θ ⋆ I,M , f ⋆ I,M ) + E f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) 2 + f (θ ⋆T I,M X) -f (θ T X) 2 + 2 Y -f ⋆ I,M (θ ⋆T I,M X) f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) + 2 Y -f ⋆ I,M (θ ⋆T I,M X) f (θ ⋆T I,M X) -f (θ T X) + 2 f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) f (θ ⋆T I,M X) -f (θ T X) dρ I,M,η,γ (θ, f ) := R(θ ⋆ I,M , f ⋆ I,M ) + A + B + C + D + E.
Computation of C By Fubini's theorem,

C = E 2 Y -f ⋆ I,M (θ ⋆T I,M X) f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) dρ I,M,η,γ (θ, f ) = E 2 Y -f ⋆ I,M (θ ⋆T I,M X) × f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) dρ 2 I,M,γ (f ) dρ 1 I,M,η (θ) .
By the triangle inequality, for f = M j=1 β j ϕ j and

f ⋆ I,M = M j=1 (β ⋆ I,M ) j ϕ j , it holds M j=1 j|β j | ≤ M j=1 j β j -(β ⋆ I,M ) j + M j=1 j (β ⋆ I,M ) j . Since f ⋆ I,M ∈ F M (C), we have M j=1 j|(β ⋆ I,M ) j | ≤ C, so that M j=1 j|β j | ≤ C + 1 as soon as f -f ⋆ I,M M ≤ 1. This shows that the set f = M j=1 β j ϕ j : f -f ⋆ I,M M ≤ γ
is contained in the support of ν M . In particular, this implies that ρ 2 I,M,γ is centered at f ⋆ I,M and, consequently,

f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) dρ 2 I,M,γ (f ) = 0.
This proves that C = 0.

Control of A Clearly,

A ≤ sup y∈R (f ⋆ I,M (y) -f (y) 2 dρ 2 I,M,γ (f ) ≤ γ 2 .
Control of B We have

B = E f (θ ⋆T I,M X) -f (θ T X) 2 dρ I,M,η,γ (θ, f ) ≤ E ℓ(C + 1)(θ ⋆T I,M -θ T )X 2 dρ 1 I,M,η (θ)
(by the mean value theorem)

≤ ℓ 2 (C + 1) 2 E X 2 ∞ θ ⋆ I,M -θ 2 1 dρ 1 I,M,η (θ) ≤ ℓ 2 (C + 1) 2 η 2 (by Assumption D). Control of E Write |E| ≤ 2 E f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) × f (θ ⋆T I,M X) -f (θ T X) dρ I,M,η,γ (θ, f ) ≤ 2 E f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) × ℓ(C + 1) (θ ⋆T I,M -θ T )X dρ I,M,η,γ (θ, f ) ≤ 2 E f ⋆ I,M (θ ⋆T I,M X) -f (θ ⋆T I,M X) 2 dρ I,M,η,γ (θ, f ) 1 2 E ℓ(C + 1)(θ ⋆T I,M -θ T )X 2 dρ I,M,η,γ (θ, f ) 1 2
(by the Cauchy-Schwarz inequality)

≤ 2 γ 2 1 2 ℓ 2 (C + 1) 2 η 2 1 2 = 2ℓ(C + 1)γη. Control of D Finally, D = 2 E Y -f ⋆ I,M (θ ⋆T I,M X) f (θ ⋆T I,M X) -f (θ T X) dρ I,M,η,γ (θ, f ) = 2 E Y -f ⋆ I,M (θ ⋆T I,M X) f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) (since f dρ 2 I,M,γ (f ) = f ⋆ I,M ) = 2E Y -f ⋆ I,M (θ ⋆T I,M X) f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) ≤ 2 E Y -f ⋆ I,M (θ ⋆T I,M X) 2 × E f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) 2 (by the Cauchy-Schwarz inequality) = 2 R(θ ⋆ I,M , f ⋆ I,M ) E f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) 2 .
The inequality

f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) ≤ ℓ(C + 1) (θ ⋆T I,M -θ T )X ≤ ℓ(C + 1) θ ⋆ I,M -θ 1 leads to f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) 2 ≤ ℓ 2 (C + 1) 2 θ ⋆ I,M -θ 1 dρ 1 I,M,η (θ) 2 .
Consequently,

f ⋆ I,M (θ ⋆T I,M X) -f ⋆ I,M (θ T X) dρ 1 I,M,η (θ) 2 ≤ ℓ 2 (C + 1) 2 η 2 ,
and therefore

D ≤ 2ℓ(C + 1)η R(0, 0)/2 ≤ √ 2ℓ(C + 1)η √ C 2 + σ 2 .
Thus, taking η = γ = 1/n and putting all the pieces together, we obtain

A + B + C + D + E ≤ Ξ 1 n ,
where Ξ 1 is a positive constant, function of C, σ and ℓ. Combining this inequality with (4.6)-(4.8) yields, with probability larger than 1 -2δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ 1 1 -λ[(2C+1) 2 +4σ 2 ] n-wλ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n 1 + λ [(2C + 1) 2 + 4σ 2 ] n -wλ R(θ ⋆ I,M , f ⋆ I,M ) -R(θ ⋆ , f ⋆ ) + Ξ 1 n + 2 M log(10(C + 1)n) + |I| log(40epn) + log 1 δ λ .

Choosing finally

λ = n w + 2 [(2C + 1) 2 + 4σ 2 ]
,

we obtain that there exists a positive constant Ξ 2 , function of L, C, σ and ℓ such that, with probability at least 1 -2δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ 2 inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n R(θ ⋆ I,M , f ⋆ I,M ) -R(θ ⋆ , f ⋆ ) + M log(10Cn) + |I| log(40epn) + log 1 δ n .
This concludes the proof of Theorem 2.1.

Proof of Corollary 2.1

We already know, by Theorem 2.1, that with probability at least 1δ,

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ inf I ⊂ {1, . . . , p} 1 ≤ M ≤ n R(θ ⋆ I,M , f ⋆ I,M ) -R(θ ⋆ , f ⋆ ) + M log(Cn) + |I| log(pn) + log 2 δ n .
By definition, for all (θ, f

) ∈ S p 1,+ (I) × F M (C), R(θ ⋆ I,M , f ⋆ I,M ) ≤ R(θ, f ). In particular, if I ⋆ is such that θ ⋆ ∈ S p 1,+ (I ⋆ ), then R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ inf 1 ≤ M ≤ n f ∈ F M (C) R(θ ⋆ , f ) -R(θ ⋆ , f ⋆ ) + M log(Cn) + |I ⋆ | log(pn) + log 2 δ n . (4.9)
Observe that, for any f ∈ F M (C),

R(θ ⋆ , f ) -R(θ ⋆ , f ⋆ ) = R p f θ ⋆T x -f ⋆ θ ⋆T x 2 dP(x, y) ≤ B 2 1 -1 [f (t) -f ⋆ (t)] 2 dt. Since f ⋆ ∈ L 2 ([-1, 1]), we may write f ⋆ = ∞ j=1 β ⋆ j ϕ j
and apply (4.9) with

f = M j=1 β ⋆ j ϕ j .
In order to do so, we just need to check that f ∈ F M (C), that is

M j=1 j|β ⋆ j | ≤ C.
But, by the Cauchy-Schwarz inequality,

M j=1 j|β ⋆ j | = M j=1 j k |β ⋆ j |j 1-k ≤ M j=1 j 2k (β ⋆ j ) 2 M j=1 j 2-2k . Thus, M j=1 j|β ⋆ j | ≤ π √ 6 M j=1 j 2k (β ⋆ j ) 2 (since, by assumption, k ≥ 2) ≤ C (since f ⋆ ∈ W(k, 6C 2 /π 2 )).
Next, with this choice of f ,

1 -1 [f (t) -f ⋆ (t)] 2 dt ≤ ΛM -2k
for some positive constant Λ depending only on k and C (see for instance Tsybakov [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]). Therefore, inequality (4.9) leads to

R( θλ , fλ ) -R(θ ⋆ , f ⋆ ) ≤ Ξ inf 1≤M ≤n ΛM -2k + M log(Cn) + |I ⋆ | log(pn) + log 2 δ n . ( 4

Some technical lemmas

Lemma 4.5 For any subset I of {1, . . . , p}, any positive integer M ≤ n and any η ∈ ]0, 1/n], let the probability measure ρ 1 I,M,η be defined by

dρ 1 I,M,η dµ I (θ) ∝ 1 [ θ-θ ⋆ I,M 1 ≤η] . Then K(ρ 1 I,M,η , µ I ) ≤ (|I| -1) log max |I|, 4 η .
Proof. For simplicity, we assume that I = {1, . . . , |I|}. Up to a permutation of the coordinates, the proof remains valid for any subset I of {1, . . . , p}.

Still for simplicity, we let θ denote θ ⋆ I,M . By a symmetry argument, it can be assumed that θ has nonnegative coordinates-this just means that θ is arbitrarily fixed in one of the 2 |I|-1 faces of S p 1,+ (I). We denote by FA this face and note that

FA =    θ ∈ (R + ) |I| × {0} p-|I| : |I| j=1 θ j = 1    .
Finally, without loss of generality, we suppose that the largest coordinate in θ is θ1 , and let ν be the uniform probability measure on FA, defined by

dν dµ I (θ) = 2 |I|-1 1 [θ∈F A] .
Set u = min(1/|I|, η/2), and let T 2 = ( θ1u, θ2 + u, θ3 , . . . , θ|I| , 0, . . . , 0), T 3 = ( θ1u, θ2 , θ3 + u, . . . , θ|I| , 0, . . . , 0), . . . . . .

T |I| = ( θ1u, θ2 , θ3 , . . . , θ|I| + u, 0, . . . , 0).

Note that u ≤ 1/|I| ≤ θ1 . Therefore, for each j, all the coordinates of T j are nonnegative. Obviously T j 1 = 1, so that, for all j, T j ∈ FA. Denoting by K the convex hull of the set { θ, T 2 , . . . , T |I| }, we also have K ⊂ FA. Next, observe that T j -θ 1 = 2u ≤ η, which implies K ⊂ {θ ∈ R p : θ -θ 1 ≤ η}.

Clearly,

K(ρ 1 I,M,η , µ I ) = log     1
1 [ θ-θ ⋆ I,M 1 ≤η] dµ I (θ)     ≤ log     1
1 [θ∈F A] 1 [ θ-θ ⋆ I,M 1 ≤η] dµ I (θ)     .
Thus,

K(ρ 1 I,M,η , µ I ) ≤ log     2 |I|-1 1 [ θ-θ ⋆ I,M 1 ≤η] dν(θ)     ≤ log     2 |I|-1 1 [θ∈K] dν(θ)     .
Observe that K is homothetic to FA, by a factor of u. This means that

Annex: Description of the MCMC algorithm

This annex is intended to make thoroughly clear the specification of the proposal conditional densities k 1 and k 2 introduced in Section 3.

Notation

To provide explicit formulas for the conditional densities k 1 ((τ, h)|(θ, f )) and k 2 ((τ, h)|(θ, f )), we first set

f = m f j=1 β f,j ϕ j and h = m h j=1 β h,j ϕ j ,
where it is recalled that {ϕ j } ∞ j=1 denotes the (non-normalized) trigonometric system. We let I (respectively, J) be the set of nonzero coordinates of the vector θ (respectively, τ ), and denote finally by θ I (respectively, τ J ) the vector of dimension |I| (respectively, |J|) which contains the nonzero coordinates of θ (respectively, |τ |). Recall that all densities are defined with respect to the prior π, which is made explicit in Subsection 2.2.

For a generic h ∈ F m h (C + 1), given τ ∈ S p 1,+ and s > 0, we let the density dens s (h|τ, m h ) with respect to π be defined by

Y i - m h j=1 b j ϕ j (τ T X i ) 2 .
In the experiments, we fixed s = 0.1. Note that simulating with respect to dens s (h|τ, m h ) is an easy task, since one just needs to compute a least square estimate and then draw from a truncated Gaussian distribution.

Description of k 1

We take It is noteworthy that when we change the parameter from θ to τ , then we also change the function from f to h. Thus, with this procedure, the link function h is more "adapted" to τ and the subsequent move is more likely to be accepted in the Hastings-Metropolis algorithm.

In the case where we are to remove one component, k 1,-is given by

k 1,-((τ, h)|(θ, f )) = j∈I c j 1 [τ =θ -j ] dens s (h|τ, m f ),
where θ -j is just obtained from θ by setting the j-th component to 0 and by renormalizing the parameter in order to have θ -j 1 = 1. We set .

The idea is that smaller components are more likely to be removed than larger ones. Finally, the density k 1,+ takes the form

  (x) = E[Y |X = x]. In the classical parametric linear model, one assumesY = θ ⋆T X + W, where θ ⋆ = (θ ⋆ 1 , . . . , θ ⋆ p ) T ∈ R p and E[W |X] = 0. Here r(x) = θ ⋆T x = p j=1 θ ⋆ j x jis a linear function of the components of x = (x 1 , . . . , x p ) T . More generally, we may define

Figure 1 :

 1 Figure 1: AUTO-MPG example: Estimated link function by the method FFourier .

  I,M,η,γ , µ I ⊗ ν M ) = K(ρ 1 I,M,η ⊗ ρ 2 I,M,γ , µ I ⊗ ν M ) = K(ρ 1 I,M,η , µ I ) + K(ρ 2 I,M,γ , ν M ). By technical Lemma 4.5, we know thatK(ρ 1 I,M,η , µ I ) ≤ (|I| -1) log max |I|,4 η .

. 10 )

 10 Letting ⌈.⌉ be the ceiling function and choosing M = ⌈(n/ log(Cn)) 1 2β+1 ⌉ in (4.10) concludes the proof.

1|I|- 1 ≤Lemma 4 . 6 2 I 1 [

 14621 [θ∈K] dν(θ) = u |I|-1 .Consequently, we obtainK(ρ 1 I,M,η , µ I ) ≤ log 2 u (|I| -1) log max |I|, 4 η .For any subset I of {1, . . . , p}, any positive integer M ≤ n and any γ ∈ ]0, 1/n], let the probability measure ρ 2 I,M,γ be defined bydρ ,M,γ dν M (f ) ∝ 1 [ f -f ⋆ I,M M ≤γ]where, forf = M j=1 β j ϕ j ∈ F M (C + 1), we put f -f ⋆ I,M M ≤γ] (f ) ζ , where ζ = 1 [ f -f ⋆ I,M M ≤γ] (f )dν M (f ).It easily follows, using the fact that the support ofρ 2 I,M,γ is included in the set {f ∈ F M (C +1) : f -f ⋆ I,M ≤ γ}, that K(ρ 2 I,M,γ , ν M ) = log(1/ζ). Note that ζ = 1 [ f -f ⋆ I,M M ≤γ] (f )dν M (f ) = 1 [ M j=1 j|β j -(β ⋆ I,M ) j |≤γ] (β)1 [ M j=1 j|β j |≤C+1] (β)dβ 1 [ M j=1 j|β j |≤C+1] (β)dβ ,where the second equality is true since ν M is (the image of) the uniform probability measure on {β ∈ R M : M j=1 j|β j | ≤ C + 1}. This impliesK(ρ 2 I,M,γ , ν M ) = log j=1 j|β j |≤C+1] (β)dβ 1 [ M j=1 j|β j -(β ⋆ I,M ) j |≤γ] (β)1 [ M j=1 j|β j |≤C+1] (β)dβ j -(β ⋆ I,M ) j + M j=1 j (β ⋆ I,M ) j . Since f ⋆ I,M ∈ F M (C), we have M j=1 j|(β ⋆ I,M ) j | ≤ C, so that 1 [ M j=1 j|β j |≤C+1] ≥ 1 [ M j=1 j|β j -(β ⋆ I,M ) j |≤γ]as soon as γ ≤ 1. We conclude that K(ρ 2 I,M,γ , ν M ) = log

  dens s (h|τ, m h ) j -βj (τ, m h ) j | ≤ C + 1 ,where the βj (τ, m h ) are the empirical least square coefficients given by βj (τ, m h ) j=1,...,m h∈ arg minb∈R m h n i=1

k 1 (

 1 •|(θ, f )) = 2k 1,= (•|(θ, f )) + k 1,+ (•|(θ, f )) 3 1 [|I|=1] + k 1,-(•|(θ, f )) + 2k 1,= (•|(θ, f )) + k 1,+ (•|(θ, f )) 4 1 [1<|I|<p] + k 1,-(•|(θ, f )) + 2k 1,= (•|(θ, f )) 3 1 [|I|=p] .Roughly, the idea is that k 1,-tries to remove one component in θ, k 1,= keeps the same number of components, whereas k 1,+ adds one component. The density k 1,= takes the formk 1,= ((τ, h)|(θ, f )) = k 1,= (τ |θ)dens s (h|τ, m f ).The density k 1,= (.|θ) is the density of τ when J = I andτ I = θ I + E θ I + E 1 sgn (θ I + E) j(θ I +E),where E = (E 1 , . . . , E |I| ) and the E i are independent random variables uniformly distributed in [-δ, δ]. Throughout, the value of δ was fixed at 0.5.

  c j = exp (-|θ j |) 1 [|θ j |<δ] ℓ∈I exp (-|θ ℓ |) 1 [|θ ℓ |<δ]

k 1 ,

 1 + ((τ, h)|(θ, f )) = j / ∈I c ′ j 1 [τ -j =θ] 1 [|τ j |<δ] 2δ dens s (h|τ, m f ).

Table 1 :

 1 T . Thus, in [Model 1] and [Model 2], even if the ambient dimension is large, the intrinsic dimension of the model is in fact equal to 2. Numerical results for the simulated data, with n = 50 and n = 100, p = 10. The characters in bold indicate the best performance.

	n = 50	p = 10	FFourier	FHHI	FLasso	FNW
	F Linear	median 0.061	0.063	0.046	0.293
		mean	0.061	0.063	0.047	0.290
		s.d.	0.016	0.014	0.011	0.063
	F SI	median 0.050	0.067	0.307	0.198
		mean	0.069	0.080	0.338	0.208
		s.d.	0.081	0.057	0.082	0.072
	F NP	median 0.375	0.405	0.830	0.354
		mean	0.402	0.407	0.890	0.336
		s.d.	0.166	0.110	0.176	0.006
	n = 100 p = 10	FFourier	FHHI	FLasso	FNW
	F Linear	median 0.053	0.051	0.042	0.227
		mean	0.056	0.050	0.043	0.237
		s.d.	0.011	0.006	0.004	0.044
	F SI	median 0.047	0.052	0.332	0.209
		mean	0.049	0.053	0.337	0.218
		s.d.	0.009	0.012	0.063	0.045
	F NP	median 0.305	0.343	0.793	0.333
		mean	0.321	0.338	0.833	0.324
		s.d.	0.092	0.042	0.145	0.041

Table 2 :

 2 Numerical

	2. The single-index methods FFourier and FHHI are the best ones when
	the targeted regression function really involves a single-index model
	[Model 2].
	3. The kernel method gives good results in the purely nonparametric set-
	ting [Model 3].

results for the simulated data, with n = 50 and n = 100, p = 50. The characters in bold indicate the best performance.

Table 3 :

 3 Numerical results for the real-life data sets. The characters in bold indicate the best performance.

	Data set		FFourier	FHHI	FLasso	FNW
	AIR QUALITY median 0.117	0.099	0.107	0.129
	n = 111	mean	0.128	0.096	0.113	0.130
	p = 3	s.d.	0.044	0.029	0.029	0.035
	AUTO-MPG	median 0.044	0.049	0.070	0.068
	n = 392	mean	0.051	0.050	0.072	0.069
	p = 7	s.d.	0.017	0.006	0.011	0.009
	CONCRETE	median 0.089	0.087	0.106	0.094
	n = 1030	mean	0.091	0.087	0.107	0.094
	p = 8	s.d.	0.008	0.003	0.005	0.004
	HOUSING	median 0.074	0.059	0.086	0.086
	n = 508	mean	0.076	0.061	0.085	0.088
	p = 11	s.d.	0.015	0.013	0.012	0.016
	SLUMP-1	median 0.289	0.171	0.201	0.208
	n = 51	mean	0.244	0.187	0.213	0.226
	p = 7	s.d.	0.062	0.050	0.049	0.047
	SLUMP-2	median 0.219	0.196	0.172	0.215
	n = 51	mean	0.216	0.194	0.171	0.213
	p = 7	s.d.	0.053	0.025	0.019	0.022
	SLUMP-3	median 0.065	0.070	0.053	0.116
	n = 51	mean	0.073	0.079	0.052	0.126
	p = 7	s.d.	0.033	0.027	0.010	0.026
	WINE-RED	median 0.173	0.171	0.183	0.171
	n = 1599	mean	0.174	0.170	0.174	0.183
	p = 11	s.d.	0.009	0.008	0.007	0.010
	WINE-WHITE median 0.191	0.187	0.185	0.184
	n = 4898	mean	0.202	0.188	0.186	0.185
	p = 11	s.d.	0.045	0.003	0.004	0.004

Table 4 :

 4 Numerical results for the real-life data sets augmented with noise variables. The characters in bold indicate the best performance.

	Augmented data set		FFourier	FHHI	FLasso	FNW
	AIR QUALITY	median 0.172	0.272	0.164	0.281
	n = 111	mean	0.244	0.291	0.163	0.291
	p = 12	s.d.	0.163	0.116	0.038	0.046
	AUTO-MPG	median 0.043	0.062	0.085	0.202
	n = 392	mean	0.044	0.072	0.086	0.203
	p = 28	s.d.	0.009	0.018	0.008	0.014
	CONCRETE	median 0.087	0.093	0.113	0.245
	n = 1030	mean	0.087	0.094	0.112	0.094
	p = 32	s.d.	0.007	0.008	0.005	0.009
	HOUSING	median 0.071	0.199	0.092	0.226
	n = 508	mean	0.075	0.181	0.095	0.227
	p = 44	s.d.	0.023	0.084	0.013	0.018
	SLUMP-1	median 0.270	0.426	0.276	0.271
	n = 51	mean	0.290	0.409	0.274	0.262
	p = 44	s.d.	0.101	0.079	0.055	0.042
	SLUMP-2	median 0.276	0.332	0.195	0.253
	n = 51	mean	0.285	0.349	0.198	0.254
	p = 44	s.d.	0.075	0.063	0.043	0.034
	SLUMP-3	median 0.079	0.371	0.061	0.372
	n = 51	mean	0.082	0.361	0.058	0.279
	p = 28	s.d.	0.025	0.079	0.013	0.031
	WINE-RED	median 0.178	0.222	0.172	0.245
	n = 1599	mean	0.176	0.226	0.174	0.246
	p = 44	s.d.	0.085	0.033	0.006	0.029
	WINE-WHITE	median 0.199	0.239	0.187	0.252
	n = 4898	mean	0.204	0.256	0.188	0.260
	p = 11	s.d.	0.091	0.041	0.005	0.019

  Pythagora's theorem). We start with four technical lemmas. Lemma 4.1 is a version of Bernstein's inequality, whose proof can be found in Massart [31, Chapter 2, inequality (2.21)]. Lemma 4.2 is a classical result, whose proof can be found, for example, in Catoni [13, page 4]. For a random variable Z, the notation (Z) + means the positive part of Z.

2

(Lemma 4.1 Let T 1 , . . . , T n be independent real-valued random variables. Assume that there exist two positive constants v and w such that, for all integers k ≥ 2, n i=1

  [START_REF] Alquier | PAC-Bayesian bounds for randomized empirical risk minimizers[END_REF] is reached for the Gibbs distribution g given by

	dg dµ	(e) =	exp [h(e)] (exp •h)dµ	, e ∈ E.
	Lemma 4.3 Assume that Assumption N holds. Set w = 8(2C + 1) max[L,
	2C + 1] and take

(since E[W i |X i ] = 0).

Acknowledgments. The authors thank three referees for valuable comments and insightful suggestions, which lead to a substantial improvement of the paper. They also thank John O'Quigley for his careful reading of the article.

 [START_REF] Alquier | PAC-Bayesian bounds for randomized empirical risk minimizers[END_REF]Corresponding author 2 Research partially supported by the French "Agence Nationale pour la Recherche" under grant ANR-09-BLAN-0128 "PARCIMONIE". 3 Research partially supported by the French "Agence Nationale pour la Recherche" under grant ANR-09-BLAN-0051-02 "CLARA". 4 Research carried out within the INRIA project "CLASSIC" hosted by Ecole Normale Supérieure and CNRS.

1

We set

where (X i ) j denotes the j-th component of X i . In words, the idea is that a new nonzero coordinate in θ is more likely to be interesting in the model if the corresponding feature is correlated with the current residual.

Description of k 2

In the same spirit, we let the conditional density k 2 be defined by

We choose

With this choice, m h = m f +1, which means that the proposal density tries to add one coefficient in the expansion of h, while leaving θ unchanged. Finally


Olivier Bokanowski

Hasnaa And
email: hasnaa.zidani@ensta.fr.

Zidani

MINIMAL TIME PROBLEMS WITH MOVING TARGETS AND OBSTACLES

Keywords: Minimal time problem, moving targets, time-dependent state constraints, motion planning, obstacle avoidance, Hamilton-Jacobi-Bellman equations, level set method, reachable set (attainable set)

We consider minimal time problems governed by nonlinear systems under general time dependent state constraints and in the twoplayer games setting. In general, it is known that the characterization of the minimal time function, as well as the study of its regularity properties, is a difficult task in particular when no controllability assumption is made. In addition to these difficulties, we are interested here to the case when the target, the state constraints and the dynamics are allowed to be time-dependent.

We introduce a particular "reachability" control problem, which has a supremum cost function but is free of state constraints. This auxiliary control problem allows to characterize easily the backward reachable sets, and then, the minimal time function, without assuming any controllability assumption. These techniques are linked to the well known level-set approaches. Partial results of the study have been published recently by the authors in SICON. Here, we generalize the method to more complex problems of moving target and obstacle problems.

Our results can be used to deal with motion planning problems with obstacle avoidance.

Introduction

The regularity and characterization of the minimal time function is a widely studied topic (see [START_REF] Cannarsa | Convexity properties of the minimum time function[END_REF], [START_REF] Colombo | Some new regularity properties for the minimal time function[END_REF], [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], [START_REF] Cardaliaguet | Optimal times for constrained nonlinear control problems without local controllability[END_REF][START_REF] Cardaliaguet | Pursuit differential games with state constraints[END_REF] [5] and references therein).

When there is no state constraint, and under some local metric properties around the target, the minimal time function T is the unique continuous viscosity solution of an Hamilton-Jacobi (HJ) equation. See [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] in the oneplayer game setting and for corresponding Hamilton-Jacobi-Bellman (HJB) equation satisfied by T .

In the presence of state constraints K (where K is a non-empty closed subset of R d ), the minimum time function (and more generally the value function of an optimal control problem) is not necessarily continuous, unless some controllability assumption is made on the boundary of K.

In the work of [START_REF] Cardaliaguet | Pursuit differential games with state constraints[END_REF], no controllability assumption is assumed and a characterization is obtained involving non-smooth analysis.

In the present paper, we follow some ideas developed in [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF] in order to treat minimum time function and capture basin for two-player games. We are mainly interested by the case of a moving target and time-dependent state constraints, and the dynamics will be time-dependent as well.

We refer to [START_REF] Varaiya | On the existence of solutions to a differential game[END_REF][START_REF] Elliott | The existence of value in differential games[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Soravia | Pursuit-evasion problems and viscosity solutions of Isaacs equations[END_REF] and references therein for an introduction and some results for deterministic two-player games with infinite horizon.

Acknowlegments. This work was supported by the EU under the 7th Framework Programme FP7-PEOPLE-2010-ITN, GA number 264735-SADCO.

Setting of the problem

Let A and B be two nonempty compact subset of R m and R p respectively. For t ≥ 0, let A t := {α : (0, t) → A, measurable} and B t := {α : (0, t) → B, measurable}. We consider a dynamics f :

R d × [0, +∞[×A × B → R d such that
(H1) f is Lipschitz continuous and, for every x ∈ R d and (α, β) ∈ A t ×B t , its associated trajectory y = y α,β x defined as the (absolutely continuous) solution of

ẏ(s) = f (y(s), s, α(s), β(s)), for a.e. s ∈ [0, t], (1a) y(0) = x. (1b)
To simplify, we shall assume that f is globally Lipschitz in all its variables (although weaker assumptions can be made for part of the results).

We also assume that (K t ) t≥0 is a family of closed sets of R d (set of "constraints"), and that (C t ) t≥0 , an other family of closed sets (set of "targets").

We consider a game involving two players, starting at time t = 0. The first player wants to steer the system (initially at point x) to the target C t in some minimal time t ≥ 0, and by staying in K t (and using her input α), while the second player tries to steer the system away from C t or from K t (with her input β). More precisely we will say that the trajectory is "admissible on [0, t]" if it satisfies the constraints on the time interval [0, t]:

y α,β x (θ) ∈ K θ , ∀θ ∈ [0, t].
We define the set of non-anticipative strategies for the first player, as follows:

Γ t := a : B t → A t , ∀(β, β) ∈ B t and ∀s ∈ [0, t], β(θ) = β(θ) a.e. θ ∈ [0, s] ⇒ a[β](θ) = a[ β](θ) a.e. θ ∈ [0, s] .
Then we are interested to characterize the following capture basin for the first player: This setting includes the case of a fixed target C t ≡ C or of a fixed constraint K t ≡ K. It also contains the particular case of a one-player game (it suffices to take B = {b 0 } a fixed value, B t = {β} a fixed constant function and then any α ∈ A t represents an admissible non-anticipative strategy).

Cap f C,K (t) := x ∈ R d , ∃a ∈ Γ t , ∀β ∈ B t , (2)
y a[β],β x (t) ∈ C t , and y a[β],β x (θ) ∈ K θ , ∀θ ∈ [0, t] Thus x ∈ Cap f C,K (t)
We are also interested in computing the minimal time function T (x) defined by

T (x) := inf t ≥ 0, ∃a ∈ Γ t , ∀b ∈ B t , y a[β],β x (t) ∈ C t and y a[β],β x (θ) ∈ K θ , ∀θ ∈ [0, t] Since we have T (x) = inf{t ≥ 0, x ∈ Cap f C,K (t)}, it is sufficient to charac- terize the sets Cap f C,K (t). It is well known that the set Cap f C,K (t) is linked to a control problem. Indeed, consider some Lipschitz continuous function ϑ 0 : R d × R → R such that ϑ 0 (x, t) ≤ 0 ⇐⇒ x ∈ C t , (3)
and consider the control problem:

u(x, t) := inf a∈Γ max β∈Bt ϑ 0 (y a[β],β x (t), t), (4)
β s.t. y a[β],β x admissible on [0, t] (5)
(with value u(x, t) = +∞ if for any strategy a there is no admissible trajectory).

Let us assume that f is affine in the variable a :

(H2) f (x, t, a, b) = f 1 (x, t, b) + f 2 (x, t, b) • a
(a similar assumption, such as f (x, t, A, b) convex for all b, could also be assumed). As in [START_REF] Cardaliaguet | A differential game with two players and one target[END_REF] we have that Cap f C,K (t) is a closed set. Under (H2), the capture basin is related to the negative region of u(., t) :

Cap f C,K (t) = {x ∈ R d , u(x, t) ≤ 0}. (6)
However the characterization of u by means of an HJB equation is not easy because of the state constraints, unless some strong assumptions are satisfied. Our aim in this paper is to give a simple way to characterize the capture basin.

In section 3, we first discuss some existing and related results. In section 4 we present our results. An appendix contains the proofs.

Discussion

For the unconstrained case, several works have been devoted to the characterization of the value function u as a continuous viscosity solution of a Hamilton-Jacobi equation, see [START_REF] Falcone | Level sets of viscosity solutions : some applications to fronts and rendez-vous problems[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. In presence of state constraints (and when K t ≡ K and is different from R d ), the continuity of this value function is no longer satisfied, unless the dynamics satisfy a special controllability assumption on the boundary of the state constraints. This assumption called "inward pointing qualification (IPQ)" condition was first introduced by [START_REF] Soner | Optimal control with state-space constraint[END_REF]. It asks that at each point of K there exists a field of the system pointing inward K. Clearly this condition ensures the viability of K (from any initial condition in K, there exists an admissible trajectory which could stay for ever in K). Under the (IPQ) condition, the value function u is the unique continuous constrained viscosity solution of a HJB equation with a suitable new boundary condition.

Unfortunately, in many control problems, the (IPQ) condition is not satisfied and the value function u could be discontinuous. In this framework, Frankowska introduced in [START_REF] Frankowska | Semicontinuous solutions of Hamilton-Jacobi-Bellman equations with degenerate state constraints[END_REF] another controllability assumption, called "outward pointing condition." Under this assumption it is still possible to characterize the value function as the unique lower semi-continuous (l.s.c. for short) solution of an HJB equation.

In absence of any controllability assumption, the function u is discontinuous and its characterization becomes more complicate, see for instance [START_REF] Soravia | Optimality principles and representation formulas for viscosity solutions of Hamilton-Jacobi equations.ii. equations of control problems with state constraints[END_REF][START_REF] Bokanowski | Deterministic state constrained optimal control problems without controllability assumptions[END_REF] and the references therein. We refer also to [START_REF] Cardaliaguet | Optimal times for constrained nonlinear control problems without local controllability[END_REF][START_REF] Cardaliaguet | Pursuit differential games with state constraints[END_REF] for a characterization based on viability theory.

Several papers in the literature deal with the link between reachability and HJB equations. In the case when K = R d , we refer to [START_REF] Mitchell | A time-dependent Hamiliton-Jacobi formulation of reachable sets for continuous dynamic games[END_REF] and the references therein. The case when K is an open set in R d is investigated in [START_REF] Lygeros | On reachability and minimum cost optimal control[END_REF]. We also refer to [START_REF] Kurzhanski | Ellipsoidal techniques for reachability under state constraints[END_REF] for a short discussion linking the reachable sets under state constraints to HJB equations. The treatment in this reference assumed a C 1 value function.

In a recent work [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF], the case K t ≡ K, C t ≡ C and with no time dependency in the dynamics, has been investigated. It was shown that the capture basin Cap f C,K can be characterized by means of a control problem whose value function is continuous (even Lipschitz continuous). Let us recall here the main idea. For simplicity we consider also the one-player game (f (x, t, a, b) ≡ f (x, a)). We first consider continuous functions g : R d → R and ϑ 0 :

R d → R such that g(x) ≤ 0 ⇐⇒ x ∈ K and ϑ 0 (x) ≤ 0 ⇐⇒ x ∈ C.
Then we introduce the new control problem:

ϑ(x, t) := inf α∈At max ϑ 0 (y α x (t)), max θ∈[0,t] g(y α x (θ)) .
It is proved in [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF] that the value function ϑ is the unique continuous viscosity solution of the equation:

min ∂ t ϑ(x, t) + H(x, D x ϑ(x, t)), ϑ(x, t) -g(x)) = 0, (7a) for t ∈ [0, +∞[, x ∈ R d , ϑ(x, 0) = max(ϑ 0 (x), g(x)). ( 7b
)
where H(x, p) := max a∈A -f (x, a) • p , and furthermore we have

Cap f C,K (t) = {x, ϑ(x, t) ≤ 0}.
The main feature of ( 7) is to use a modelization with a supremum cost, in order to handle easily the state constraints and to determine the corresponding capture basins. This idea generalizes the known level-set approach usually used for unconstrained problems. Moreover, the continuous setting opens a large class of numerical schemes to be used for such problems (such as Semi-Lagrangian or finite differences schemes). We refer to [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF] for numerical results and comparison of various approaches for state-constrained problems.

We shall now consider the general problem of moving (or time-dependent) targets as well as moving obstacles.

N otations. Throughout the paper | • | is a given norm on R d (for d ≥ 1). For any closed set K ⊂ R d and any x ∈ K, we denote by d(x, K) the distance from x to K: d(x, K) := inf{|x -y|, y ∈ K}. We shall also denote by d K (x) the signed distance function to K, i.e., with d K (x) := d(x, K) for x / ∈ K, and d K (x) := -d(x, R d \K) for x ∈ K.

Main results

We assume that (H3) the set-valued applications θ K θ and θ C θ are upper semicontinuous.

We recall that if (Q t ) t≥0 denotes a family of subsets of R d+1 , then the set valued map t Q t is said to be "upper semi-continuous" if

∀ > 0, ∃α > 0, ∀θ ∈ [t -α, t + α], Q θ ⊂ Q t + B(0, 1).
Remark 1. For every t ≥ 0, the set Cap f C,K (t) contains the initial positions which can be steered to the target (exactly) at time t. Of course, we can also define the "backward reachable set", which is the set of points from which one can reach the target C t before time t. This set is also a capture basin for the dynamics f where

f (x, t, (a, λ), b) := λf (x, t, a, b), λ ∈ [0, 1]
(see [START_REF] Mitchell | A time-dependent Hamiliton-Jacobi formulation of reachable sets for continuous dynamic games[END_REF][START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF]).

We start by embedding the position y(t) and time t into a "space time" space. To do so, we set for every z = (y, t) ∈ R d × R, the set-valued map

F : R d+1 × [0, +∞) × A × B R d+1 such that : F (z, a, b) := f (y, t, a, b) × 1 ,
and we remark that F satisfies similar Lipschitz continuity and linearity assumptions as in (H1) and (H2). For a given ξ ∈ R d ×R and (α, β) ∈ A t ×B t , we can then consider z = z α,β ξ , the absolutely continuous solution of ż(s) = F (z(s), α(s), β(s)), a.e. s ∈ [0, t], z(0) = ξ.

(we shall simply denote z = z ξ if there is no ambiguity). Any solution z ξ (s) = (y(s), η(s)) of the previous system satisfies equivalently, if ξ = (x, t 0 ), ẏ(s) = f (y(s), t 0 + s, α(s), β(s)), s ≥ 0, y(0) = x, η(s) = t 0 + s, s ≥ 0 Moreover, let also introduce two subsets of R d+1 :

C := t≥0 C t × {t} and K := t≥0 K t × {t}.
We have the following elementary result : Lemma 1. Under (H3), the sets C and K are closed subsets of R d+1 .

Hence there exists Lipschitz continuous functions ϑ 0 : R d+1 → R and g : R d+1 → R such that

ϑ 0 (ξ) ≤ 0 ⇔ ξ ∈ C. ( 8
)
and

g(ξ) ≤ 0 ⇔ ξ ∈ K (9)
(for instance we may choose ϑ 0 (ξ) := d C (ξ) and g(ξ) := d K (ξ)). In particular, for any t < 0 we have ϑ 0 ((x, t)) > 0 and g((x, t)) > 0.

We can then define a capture basin associated to the new dynamics F :

Cap F C,K (τ ) := x ∈ R d+1 , ∃a ∈ Γ τ , ∀β ∈ B τ , (10) z a[β],β ξ (τ ) ∈ C, and z a[β],β ξ (θ) ∈ K, ∀θ ∈ [0, τ ]
Notice that in the case ξ = (x, 0), we have z ξ (t) = (y x (t), t) where y x is a trajectory for the dynamics f . Hence

y x (t) ∈ C t ⇔ z ξ (t) ∈ C,
and in the same way,

y x (t) ∈ K t ⇔ z ξ (t) ∈ K.
Therefore we can easily deduce the following result: Proposition 1. For all t ≥ 0, we have

x ∈ Cap f C,K (t) ⇔ (x, 0) ∈ Cap F C,K (t). (11)
Since Cap F C,K (s) has a fixed state constraint C and fixed target K and an autonomous dynamics F , we can use the results of [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF].

We consider the control problem, for ξ ∈ R d+1 and τ ≥ 0:

ϑ(ξ, τ ):= min a∈Γτ max β∈Bτ max ϑ 0 (z a[β],β ξ (τ )), max θ∈[0,τ ] g(z a[β],β ξ (θ)) (12)
where we recall that the Lipschitz function g is related to the obstacle K by (9) (we note that for t < 0 and τ ≥ 0, we have ϑ(x, t, τ ) > 0). It is the use of the supremum norm that will enable us to deal with the controllability problem, because now [START_REF] Falcone | Level sets of viscosity solutions : some applications to fronts and rendez-vous problems[END_REF] has no "explicit" state constraint. In fact, in this new setting, the term max θ∈[0,τ ] g(z ξ (θ)) plays a role of a penalization that a trajectory z ξ would pay if it violates the state-constraints. Theorem 2 will show the advantage of considering [START_REF] Falcone | Level sets of viscosity solutions : some applications to fronts and rendez-vous problems[END_REF], because ϑ will be characterized as the unique continuous solution of an HJB equation.

Theorem 1 (Characterization of the capture basin). Assume (H1)-(H3). Let ϑ 0 (resp. g) be Lipschitz continuous functions satisfying (8) (resp. (9)). Let ϑ be the value function defined by [START_REF] Falcone | Level sets of viscosity solutions : some applications to fronts and rendez-vous problems[END_REF]. For every τ ≥ 0, we have:

Cap F C,K (τ ) = ξ ∈ R d+1 , ϑ(ξ, τ ) ≤ 0 .
In particular, we have,

Cap f C,K (τ ) = x ∈ R d , ϑ(x, 0, τ ) ≤ 0 .
Now, the function ϑ can be characterized as the unique solution of a Hamilton-Jacobi equation. More precisely, considering the Hamiltonian

H F (ξ, P ) := max a∈A min b∈B -F (ξ, a, b) • P , (13)
by using [START_REF] Barron | The Bellman equation for minimizing the maximum cost[END_REF] (see also [START_REF] Quincampoix | A viability approach for optimal control with infimum cost[END_REF][START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF]), we have Theorem 2. Assume (H1), and that ϑ 0 and g are Lipschitz continuous. Then ϑ is the unique continuous viscosity solution of the variational inequation (or "obstacle" problem)

min(∂ τ ϑ + H F (ξ, D ξ ϑ), ϑ -g(ξ)) = 0, (14a) τ > 0, ξ ∈ R d+1 , ϑ(ξ, 0) = max(ϑ 0 (ξ), g(ξ)), ξ ∈ R d+1 . (14b)
For sake of completeness the notion of viscosity solution is recalled in the appendix (see Definition 1).

Application. We are thus able to compute capture basin using regular functions. There are numerous schemes that can approximate the value function ϑ of the previous HJB or HJI equations. This gives a way to compute the set Cap F C,K (t). Then in view of Theorem 1 we can find the points x that belong to Cap f C,K (t).

Remark 2. From a theoretical point of view, the choice of g is not important, and g can be any Lipschitz function satisfying [START_REF] Cardaliaguet | Pursuit differential games with state constraints[END_REF]. Of course, the value function ϑ is dependent on g, while the set {ξ ∈ R d+1 , ϑ(ξ, t) ≤ 0} does not depend on g.

There are other informations in the function ϑ. For instance, let y α,β

x,t denotes the solution of the differential equation (1a) and such that y(t) = x (instead of y(0) = x). For s ≥ t, we define

Cap f C,K (t; s) := x ∈ R d , ∃a ∈ Γ, ∀β ∈ B, (15)
y a[β],β x (s) ∈ C s , and y a[β],β x (θ) ∈ K θ , ∀θ ∈ [t, s]
Then we have Proposition 2. For any τ ≥ 0,

Cap f C,K (t; t + τ ) = {x ∈ R d , ϑ(x, t, τ ) ≤ 0}.
Indeed this comes from the fact that y α,β x (s) = y ᾱ, β x,t (s + t) (where ᾱ(s) := α(s -t) and β(s) := β(s -t)), and thus for ξ = (x, t) we can deduce that

z ξ (τ ) ∈ C ⇔ y x,t (τ ) ∈ C τ and z ξ (τ ∈ K ⇔ y x,t (τ ) ∈ K τ .
In other words, ϑ(x, t, τ ) ≤ 0 (for some τ ≥ 0) is equivalent to say that there exists some non-anticipative strategy a[•] such that (for any adverse strategy β) we can reach the target C t+τ at time t + τ starting from x at time t.

Appendix A. Proofs of the main results

For sake of simplicity these proofs are given in the one-player game setting. Proof of Theorem 1. We reproduce here the idea of the proof that can be found in [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF]. Assume that ξ ∈ Cap F C,K (τ ). Then there exists an admissible trajectory z ξ such that ϑ 0 (z ξ (t)) ≤ 0, and z ξ (θ) ∈ K for every θ ∈ [0, τ ].

Hence, max θ∈[0,τ ] g(z ξ (θ)) ≤ 0, and we have:

ϑ(ξ, τ ) ≤ max(ϑ 0 (z ξ (τ )), max θ∈[0,τ ] g(z ξ (θ))) ≤ 0.
Conversely, assume that ϑ(ξ, τ ) ≤ 0. Then there exists a trajectory z ξ for the dynamics F , such that

0 ≥ ϑ(ξ, τ ) = max(ϑ 0 (z ξ (τ )), max θ∈[0,τ ] g(z ξ (θ))).
Thus, for all θ ∈ [0, τ ], g(z ξ (θ)) ≤ 0, i.e. z ξ (θ) ∈ K, and z ξ is an admissible trajectory. Moreover, we have ϑ 0 (z ξ (τ )) ≤ 0, hence z ξ (τ ) ∈ C and we can conclude that ξ ∈ Cap F C,K (τ ). Proof of Theorem 2. It is based on a dynamic programming principle (DPP) for ϑ that we shall not reproduce here (see for instance [START_REF] Barron | The Bellman equation for minimizing the maximum cost[END_REF]Proposition 3.1]).

We recall here the definition of viscosity solution for [START_REF] Cardaliaguet | A differential game with two players and one target[END_REF] (the definition in the case of ( 14) is similar). Definition 1 (Viscosity solution). An upper semi-continuous (resp. lower semi-continuous) function ϑ : R d × R + → R is a viscosity subsolution (resp. supersolution) of (7) if ϑ(x, 0) ≤ ϑ 0 (x) in R d (resp. ϑ(x, 0) ≥ ϑ 0 (x)) and for any (x, t) ∈ R d ×(0, ∞) and any test function φ ∈ C 1 (R d ×R + ) such that ϑ-φ attains a maximum (resp. a minimum) at the point (x, t) ∈ R d ×(0, ∞), then we have min(∂ t φ + H(x, ∇φ), ϑ -g(x)) ≤ 0 (resp. min(∂ t φ + H(x, ∇φ), ϑ -g(x)) ≥ 0) .

A continuous function ϑ is a viscosity solution of (7) if ϑ is a viscosity subsolution and a viscosity supersolution of [START_REF] Cardaliaguet | A differential game with two players and one target[END_REF].

The fact that ϑ is the unique solution of ( 7) follows from the comparison principle for (7) (which is classical, see for instance [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF]), and the fact that the Hamiltonian function H satisfies for some constant C ≥ 0 and for all x i , p i , x and p in R d .

|H(x 2

 2 , p) -H(x 1 , p)| ≤ C(1 + |p|) |x 2 -x 1 |, |H(x, p 2 ) -H(x, p 1 )| ≤ C|p 2 -p 1 |,

  means that there exists a non anticipative strategy a ∈ Γ t such that for any adverse strategy β ∈ B t , we have y

	a[β],β x	(t) ∈ C t (we
	reach the target C t at time t).