Serge Cohen

Sébastien Dejean

Sébastien

Sébastien Déjean

Sébastien Gadat

Adaptive sequential design for

Keywords: Active Learning, Optimal Designs, Stochastic Algorithms, Model Selection, Wavelets, MSC 62K05, 62L05, 93E35

Introduction

A new algorithm for building optimal designs is proposed in this paper to estimate an unknown signal η. We will use a one dimensional framework but all our results can be extended to multi-dimensional cases provided E is a finite union of finite dimensional compact and convex spaces. We discriminate competing models in a multi-resolution setting described by a set of functions (Λ r,t ) (r,t) . The signal is observed with the addition of a noise

∀x ∈ E f (x) = η(x) + σζ(x), (1)
where ζ(x) will denote a normalized white noise and σ 2 is the variance of the model. To estimate η, we allow ourself as few observations of f as possible but we can choose the locations where f is evaluated to build our design. For any set I of pairs (r, t) where r is a resolution index and t a translations factor, and any design x, we estimate η as η(x) = (r,t)∈I θr,t Λ r,t (x).

Of course, θ depends on I and x and we aim to obtain a good estimation of η with a small variance. We suppose that the cost of the evaluation of f is larger than any optimisation procedure to get I and x and our algorithm will have to select a few points among E to provide good estimations of η. A good example of such situations can be encountered in biology where experiments are often expensive and one cannot make many experiments. Another example is provided in the experimental section, dealing with the Motorcycle impact experiment introduced in [S85]. Hence, our motivation is to obtain a multi-resolution estimation with a few number of points in the design.

In the last decade, a large amount of works concerned by models selection based on stochastic or penalized approaches has been successfully developed. From a theoretical point of view, the recent ℓ 1 penalized approaches (Lasso of [T96], LARS of [START_REF] Efron | Least angle regression[END_REF] or Dantzig Selector of [START_REF] Candès | The Dantzig selector: Statistical estimation when p is much smaller than n[END_REF]) propose to introduce some ℓ 1 penalties to constraint the estimator θ to obtain sparse estimation η in the basis (Λ r,t ) (r,t) . These works reach satisfactory theoretical results for the selection of the model. But to the best of our knowledge, no exact optimization method for selecting a design dedicated to the sparse reconstruction of the model exists. Implicitly, the ℓ 1 penalized regression methods rely on properties of the design to obtain some good reconstruction but no work has been proposed to choose explicitly a design in order to satisfy such properties. The work [START_REF] Castaño | Adaptive fitting of scattered data by spline wavelets Curve and surface fitting[END_REF] uses an heuristic adaptive thresholding technique and propose a model selection procedure that deals with any (non-regular) design. One may also mention the interesting work [START_REF] Liu | Exact Wavelet Designs for Discrimination[END_REF]: their authors propose a construction of optimal and sequential designs for nonparametric regression models in the Daubechies basis with a simulated annealing strategy. This approach is also coupled with a model selection strategy. At last, the paper [START_REF] Kerkyacharian | Regression in random design and warped wavelets[END_REF] considers a warping method of the space state and then thresholds the wavelet coefficients with non-regular designs. In [START_REF] Amato | Wavelet kernel penalized estimation for non-equispaced design regression[END_REF], some methods are provided to use wavelet bases with irregular designs using wavelet kernel penalized estimation. In the experimental section, we will provide comparisons of our method with ℓ 1 penalized approach (the LARS estimator of the Adaptive Lasso [Z06]) and the thresholding methods of [START_REF] Amato | Wavelet kernel penalized estimation for non-equispaced design regression[END_REF].

Regarding now the community of optimal designs, many works have been concerned by finding designs represented as a continuous measure on the state space E. But there exist a few amount of explicit results for discrete optimal designs. Moreover, a large amount of these numerical methods yield some continuous designs although from a practical point of view, discrete designs are easier to handle. Some results can be found in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] but most of them are obtained by the optimization of numerical algorithms like simulated annealing and therefore are not explicit. Concerning the special case of multiresolution bases and optimal designs, several authors have already used a simulated annealing procedure for building robust minimax integer-valued designs (see for instance [START_REF] Fang | Integer-valued, minimax robust designs for estimation and extrapolation in heteroscedastic, approximately linear models[END_REF], [START_REF] Liu | Exact Wavelet Designs for Discrimination[END_REF] or [START_REF] Oyet | On exact minimax wavelet designs obtained by simulated annealing[END_REF]). Note that simulated annealing may be a costly task when the dimension of the problem becomes large. Among these methods, few ones propose a model selection approach except [START_REF] Liu | Exact Wavelet Designs for Discrimination[END_REF] or [START_REF] Biswas | An efficient design for model discrimination and parameter estimation in linear models[END_REF] where the model selection strategy is implemented with an hypothesis testing approach.

We want to expand the signal η in a multi-resolution basis (Λ r,t ) r,t . We successively select some points x i where f is evaluated to build a n-set x = {x 1 , . . . x n } and elements of (Λ r,t ) r,t to reach a correct estimation η of η. But for numerical reasons, it is infeasible to explore all possible subsets of functions and all n-sets x to choose among them the best model and design. In our settings, the design x = {x 1 , . . . , x n } will be adapted to the sequential ealuations of f . We propose to develop a recursive strategy: we iteratively evaluate f with a sequence (x 1 , . . . x n ) of points in E and find both a design x n = (x 1 , . . . x n ) and a family of linearly independent functions in an optimal way to estimate η.

Our contribution is twofold, we first provide an algorithm based on multi-resolution bases and adaptive strategy. This algorithm is then applied with several bases. In the very special case of triangle Schauder basis, we provide also a new theoretical result about the locations of discrete optimal designs. This result will be useful for numerical implementations. and can be tensorized to dimensions larger than one. We will use some classical ideas of optimal design theory such as D-or A-optimal designs (see e.g. [START_REF] Khuri | Response Surfaces: Designs and Analyses[END_REF][START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] for general classical reminders on the subject). The model selection strategy will use ideas coming from the adaptive regression point of view developed with Multivariate Adaptive Regression Splines (MARS) described in [F91] and classical coarse to fine multi-resolution methods (see e.g. [M90, DJ94, BG05]). Given any initialization (x 0 , I 0 ), our method will produce x n+1 from x n by the addition of a new point x n+1 ∈ E. We will not remove an element of x n at step n + 1 as it would amount to loose some available information (which is costly!), thus x n ⊂ x n+1 . The variance of our model will be controlled by the construction of a suitable design dedicated to I n . The computation of I n will rely on a stochastic approach that aims to optimize the bias of our linear reconstruction. The paper is organized as follows: next section presents our model and our algorithm. Section 3 recalls optimal design theory and the location properties of optimal designs. Section 4 precisely describes the stochastic algorithm for the model selection. At last, Section 5 provides some experimental comparisons, especially with ℓ 1 -penalized approaches which are widely used now and the kernel penalized wavelet estimator described in [START_REF] Amato | Wavelet kernel penalized estimation for non-equispaced design regression[END_REF].

Model

Integrated Mean Square Error

Assume E = [0; 1], η is supposed to be expanded in a multi-resolution basis denoted by (Λ r,t ) r≥0,0≤t≤2 r -1 :

η = r,t θ r,t Λ r,t .
(2)

Note that this assumption is true as soon as η belongs to some Besov space. More precisely, we will assume that for some unknown s > 0, η belongs to some homogeneous Besov space Ḃs,2 2 . These spaces are described for instance in the chapter VI, paragraph 9 of [M90]. We will not technically detail this point as it is not in the scope of this work, but note that such spaces contain a very large class of functions. In the sequel, we will use a generic notation I to refer to a list of pairs I = {(r 1 , t 1 ), . . . (r p , t p )}. I n will thus describe the functions used at step n among the whole set (Λ r,t ) r≥0,2 r >t≥0 .

We recall that

f (x) := η(x) + σζ(x), (3)
where ζ(x) is a white noise. In the sequel, we will suppose σ is unknown and we will provide a procedure to estimate it. Let us suppose that the design x and the set I are given. Let us denote by ηx,I the least square estimator of η with the linear model based on (Λ (r,t) ) (r,t)∈I and the observations (f (x i )) x i ∈x . One can use a standard Integrated Mean Square Error (IMSE) and define

J(x, I) = E E [η x,I (u) -η(u)] 2 du. (4)
We also denote by f (x) the column vector of the function f observed at the points of the design x and we use the notation ΛI (x) for the rectangular (p × n) matrix:

ΛI (x) =    Λ (r 1 ,t 1 ) (x 1 ) . . . Λ (r 1 ,t 1 ) (x n ) . . . . . . . . . Λ (rp,tp) (x 1 ) . . . Λ (rp,tp) (x n )    . (5)
We can expand (4) as

J(x, I) = E Var[η x,I (u)]du :=V x,I + E (E[η x,I (u)] -η(u)) 2 du :=B x,I . (6)
With standard notations, M x,I denotes the information matrix of the design x and functions Λ (r 1 ,t 1 ) , . . . Λ (rp,tp) , i.e. M x,I = ΛI (x) t ΛI (x). Let us denote by µ 1,1 (I) the integral E t ΛI (u) ΛI (u)du. Some immediate computation (see e.g. [F69] ) yields

V x,I = σ 2 T r µ 1,1 (I)M -1
x,I . The effect of the design on the bias term is not explicit and thus it will not be possible to find analytic criteria which describe precisely the influence of x on B x,I . To overcome this difficulty, it is however possible to adopt the minimax approach described for instance in [START_REF] Oyet | Exact Wavelet Designs for Discrimination[END_REF]. Given a current set of pairs I, let us denote by

I c {(r, t) | r ≥ 0, 0 ≤ t ≤ 2 r -1 , (r, t) / ∈ I}. We can decompose η in η I + η I c where η I = (r,t)∈I θ r,t Λ r,t and η I c = (r,t)∈I c θ r,t Λ r,t .
Then, it is clear that η I c will not be estimated by our estimator constructed only with the set I. If (Λ r,t ) r≥0,2 r >t≥0 is orthonormal, the bias can be decomposed as

B x,I = 1 0 [Eη x,I (u) -η I (u)] 2 + η I c (u) 2 du.
Obviously, since η I c is unknown, B x,I remains untractable but it is possible to use a worst case approach with the minimax design theory. We follow the notation used in [START_REF] Oyet | Exact Wavelet Designs for Discrimination[END_REF]. To ensure an equilibrium between the contribution of the bias and the variance to the error, we can impose a bound on the magnitude of the remainder η I c 2 ≤ τ for a constant τ > 0. We denote the least favorable case B ⋆ x,I,τ :

B ⋆ x,I,τ = sup η I c 2 ≤τ B x,I = sup η I c 2 ≤τ 1 0 [Eη x,I (u) -η I (u)] 2 + η I c (u) 2 du.
Since, J(x, I) ≤ B ⋆ x,I,τ + σ 2 V x,I , our next algorithm will detail the choice of I and x. Its dependence on the ratio ν = σ 2 /τ 2 corresponds to a bias/variance tradeoff.

The adaptive algorithm

In our adaptive framework, we need to choose successively new points of the design x and we decide to update the set I or to keep this set unchanged. As pointed in the introduction, we do not delete points of the design x. Consequently, the algorithm will necessary be of the form:

1. Step 0 • Fix any initial set of functions I 0 . For instance in a one dimensional setting with E = [0; 1],

we choose I 0 = {(0, 0); (1, 0); (1, 1)} as we do not have any prior on the support concerning η.

• For a given integer n 0 , compute the optimal design x 0 of size n 0 which minimizes the bound of the IMSE:

x 0 = arg min x B ⋆ x,I 0 ,τ + νT r µ 1,1 (I 0 )M -1 x , I 0 . (7)
2.

Step n ⋄MS⋄ Model Selection step: update the set I n to "optimize" J(x n , .). Our strategy will build I n+1 from I n and will rely on a random choice of an element (r, t) ∈ I n . We will decide if we delete (r, t) or if we add an element that does not belong to I n . We will explain the stochastic model selection in section 4.

⋄OD⋄ Optimal Design: choose x n+1 deduced from x n with the addition of one point x n+1

x n+1 = x n ∪ {x n+1 }, using the former set I n+1 . This optimization is described in section 3 and mimic equation (7).

3 Optimization of the design (⋄OD⋄ step)

In this part, we explain how to obtain the point x n+1 at step ⋄OD⋄ of the iteration n in the algorithm above. We first provide some details when ν = +∞, the case where only the variance term is considered. Next, we describe how we can deal with ν < +∞.

Generalities on ⋄OD⋄ step when ν = +∞

We will consider in this part two optimality criteria for the design based on the information matrix M x,I (for more general classical criteria, one may refer to [W92]). D-optimal designs are defined by arg min x Φ 0 (M x ) where we denote Φ 0 (M x ) := det M -1

x although A-optimal designs are based on the minimization of Φ 1,C where Φ 1,C (M x ) = T r CM -1

x . In most cases, these two criteria (and many others described in [W92]) do not yield equivalent designs but each of them aims to control the variance of the linear model. We denote by Ṽ1 (x, I) := Φ 0 (M x ) as well as Ṽ2 (x, I) = Φ 1,Id (M x ) = T r(M -1

x ), and Ṽ3 (x,

I) = Φ 1,µ 1,1 (I) (M x ) = T r(µ 1,1 (I)M -1
x ). In our work, we have mainly focused on the D-optimal designs for multi-resolution bases except in a very special case where some link with A-optimal designs is established.

Remark 1 Since we do not handle continuous measures, it is impossible to easily recover classical results for D and A optimal designs using some equivalence theorems stated in [START_REF] Kiefer | Optimum designs in regression problem[END_REF][START_REF] Fedorov | Theory of Optimal Experiments[END_REF][START_REF] Karlin | Optimal experimental designs[END_REF]. In the sequel, we claim some properties for the location of optimal designs for the Schauder and Haar bases.

Our adaptive strategy is based on solutions of

x n+1 ∈ arg min x Ṽi (x ∪ x, I n+1 ), ( 8
)
where Ṽi is defined by one of the previous criteria. Note that this step does not require any evaluation of f and the simplest way to find x n+1 seems to use some standard optimization algorithm (gradient descent or simulated annealing).

3.2 Optimal designs with Schauder/Haar basis when ν = +∞

In the sequel, the Haar mother function is H 0,0 (x) = 1 x∈[0;1/2[ -1 x∈[1/2;1] and we define the Schauder triangle mother function as S 0,0 (x) =

x 0 H 0,0 (t)dt. These functions H 0,0 and S 0,0 are examples of functions Λ 0,0 . For instance, Haar and Schauder bases are then uniquely defined by classical dilatations and translations Λ r,t (x) = 2 r/2 Λ 0,0 (2 r xt).

Some elements of the Schauder basis can be seen on Figure 1. A very important issue with Schauder basis is the fact that for every n, the argmin in (8) always occur on dyadic points x = t 2 r , whose resolutions r are bounded by the maximal resolution of the elements in I. This fact clearly speeds up the numerical resolution of equation (8). We will give a theoretical proof of this fact later. Figure 2 presents for instance the behaviour of Ṽ1 . Given any design x, we are searching the location of x ∈ E such that x ∪ {x} is D-optimal. Next theorem shows that in fact, x is necessarily a dyadic point which is a singular point of functions described by I. If Supp(Ψ) is the support of any function Ψ, we define S(I) as

S(I) = (r,t)∈I arg max x∈E Λ r,t (x) (r,t)∈I {∂Supp(Λ r,t )} .
We obtain the following result whose proof is deferred to the appendix.

Theorem 1 Let x be any fixed design and (Λ r,t ) (r,t)∈I be any finite subset of functions of the Schauder triangle basis, then arg min x Ṽ1 (x ∪ x, I) ∈ S(I).

This result is very useful for our adaptive strategy to build x n+1 from the design x n at step n. Indeed, it is sufficient to explore the small finite number of dyadic points S(I), described above and to select the point which minimizes the D-criterion. We also state a similar result for the case of ridge regession estimators.

Theorem 2 Let x be any fixed design and (Λ r,t ) (r,t)∈I be any finite subset of functions of the Schauder triangle basis, then for any θ > 0:

arg min x det (M x∪x,I + θI d ) -1 ∈ S(I).
Conjecture 1 We conjecture also that for any symmetric positive matrix C and for any t > 0, the minima of det tC + M -1 x∪x,I remain in S(I).

Even if the proof of the conjecture is still open, it is supported by many numerical experiments. Concerning the other criteria, given the last observations x, we search for the x which minimizes the criteria Ṽ2 or Ṽ3 . One can re-write the A-optimal design criterium using the useful Schur factorization. For the sake of simplicity, let µ := µ 1,1 (I). Then

T r µM -1 x∪x = T r µM -1 x - µM -1 x Λt I ΛI M -1 x 1 + t ΛI M -1 x ΛI = T r µM -1 x - t ΛI M -1 x µM -1 x ΛI 1 + t ΛI M -1 x ΛI .
Thus, the location of the optimal x at the step n + 1 is deduced from the step n by maximizing the second term of the last equation. This can be performed easily in view of next theorem.

Theorem 3 Let x be any fixed design and (Λ r,t ) (r,t)∈I be any finite subset of functions of the Schauder triangle basis. Denote by C any symmetric non negative matrix, if the Conjecture 1 is true, then arg min

x T r CM -1 x∪x,I ⊂ S(I).

We can deduce the ⋄OD⋄ step at iteration n from the last theorems for the Schauder triangle basis when the family of functions in I n+1 and the current design x n are given:

• We consider the finite set S(I n+1 ).

• Compute the criterion for each element in S(I n+1 ) and choose x n+1 in arg min

x∈S(I n+1 ) Ṽi (x n ∪x, I n+1 ).
The study of Haar basis is much simple since the function x → Ṽi (x ∪ x, I) is constant on the support of the elements of I which are intervals. Hence, it is sufficient to compute Ṽi (x ∪ x, I) at the most 2 × |I| times to obtain its maximum value and then find an optimal point into this interval. From this simple remark, one can deduce the ⋄OD⋄ step at iteration n for the Haar basis when the family of functions I n+1 and the current design x n are given:

• Rank by increasing order the elements of S(

I n+1 ) = {x 1 < x2 • • • < xl }. • Sample l + 1 elements (ξ j ) j=1...l+1 in each interval ]0; x1 [, . . . ]x j , xj+1 [, . . . ]x l ; 1[ and choose x n+1
as a solution of the optimisation problem arg min

x∈{ξ 1 ,...ξ l+1 } Ṽi (x n ∪ x, I n+1 ).

Optimal designs with general multi-resolution bases when ν = +∞

In the general case of multi-resolution wavelet bases such as Meyer basis, our approach consists in approximating solutions of several optimization problems. For D-optimal design, the algorithm uses the formula (15) provided in the appendix. The sequential D-optimal design is then equivalent to find at each step n the maximum of

x → t ΛI n+1 (x)M xn,I n+1 ΛI n+1 (x).
For A-optimal design, we use the Schur factorization formula. Hence to find the sequential A-optimal design is equivalent to find at each step n the maximum of the function

x → 1 + t ΛI n+1 (x)M xnI n+1 ΛI n+1 (x) t ΛI n+1 (x)M -1 xn,I n+1 µM -1 xn,I n+1 ΛI n+1 (x)
, where µ is either Id or µ 1,1 (I). This can be done either by using gradient descent algorithm or an exhaustive search when the set E is of small dimension (this will be the case in our simulations). Of course, the ⋄OD⋄ step in this case is much longer than the one in the Haar or Schauder case.

⋄OD⋄ step with other regressors when ν = +∞

One may naturally wonder what can be done with other regressor. Especially we think of penalized regressors. Penalized methods can be viewed as the minimization of some penalized likelihood L

θ := arg min θ {log L(x n , f (x n ), θ)) + λN (θ)} ,
where N (θ) denotes the norm of penalization. Ridge regression: For N = . 2 2 , we obtain the ridge regressor. For a fixed λ, the covariance of the ridge regression is M xn,I n+1 + λId -1 . Theorem 2 yields that optimal designs for the Schauder basis are also located in S(I n ). It is thus possible to infer some ⋄OD⋄ step for the Schauder basis using Ridge regression at step n: for each x ∈ S(I n ), compute the covariance matrix and one of the associated criterion Ṽi where M xn,I n+1 is replaced by M xn,I n+1 + λId. Then just pick the argmin x in S(I n ) and update the design x n with it. Note that this method can be obviously extended to the simpler case of Haar basis. For more general wavelet bases, the price to pay is an intensive use of optimisation algorithm to find in E the argmin points. Lasso regression: The choice of the penalization terms N = . 1 yields the lasso regressor. The effect of such a penalization on the covariance of θ is more complicated. Following the work [K01] on the soft thresholding estimator, the estimate of the covariance is approached using a Delta method by the matrix

V ar( θGsoft ) = H( θ) + λΓG( θ, σ) -1 Σ( θ) H( θ) + λΓG( θ, σ) -1
, where H (resp. Σ) is the (resp. expected) Hessian matrix of the likelihood, where G is a gaussian kernel smoothed with bandwith σ and where Γ is the diagonal matrix Γ = diag(I{θ < σ}/σ). Indeed, V ar( θGsoft ) depends on the observation x n+1 in the design which is unknown at the step n since the thresholding matrix Γ depends on f (x n+1 ). It is thus to the best of our knowledge very complicated to obtain this sequential design strategy. Choice of λ: In the penalized approach (lasso, ridge, elastic net, . . . ), an efficient way to choose λ is to use a data-driven calibration cross-validation. Indeed, such a choice increases V ar( θGsoft ) which depends on the x selected to update x n . Thus, we chose not to further investigate this question.

Minimax Optimal Design (ν < ∞)

In the general case ν = σ 2 τ 2 < ∞, we take into account the effect of the least favourable bias to choose our design. We have to find at step n the point x n that optimizes the criterion Φ

∀x ∈ E Φ(x) = B ⋆ xn∪x,In,τ + σ 2 T r µ 1,1 (I n )M -1 xn∪x,In .
The optimization of Φ is a very difficult task since B ⋆ xn∪x,In,τ is not explicit. We adopt a suboptimal approach and deal with discrete designs located on dyadic points. Given a sufficiently large maximal resolution r max ∈ N, we restrict our space of admissible x's to the set of dyadic points E max = {k2 -rmax , k = 0, . . . , 2 rmax } and we denote by 0 := ξ 1 < • • • < ξ 2 rmax +1 = 1 the points in E max . One can describe an approximation Φ(x) of Φ(x) using E max instead of E as in [START_REF] Oyet | Exact Wavelet Designs for Discrimination[END_REF]. For the sake of simplicity, denote d = 2 rmax + 1 and p = |I n+1 | the number of elements of I n+1 . Given any dyadic design x = {x 1 , . . . , x m } of length m, let us denote by P (x) the diagonal matrix of occupation measure of E max among x:

∀(k, l) ∈ {1 . . . d} × {1 . . . d} P k,l (x) = |{j|x j = ξ k }| m 1 k=l .
The application η I c ∈ L 2 (E) is replaced by its approximation vector ηI c ∈ ℓ 2 (E max ) = (η I c (ξ i )) i=1...d and we replace the constraint η I c 2 ≤ 1 by ηI c 2 2 ≤ d. If we denote by

Q I n+1 = Λ I n+1 (ξ 1 ), . . . , Λ I n+1 (ξ d ) of size p × d, the Average Squared Bias (ASB) is then approximated by ASB ≃ 1 d Q I n+1 (Q t I n+1 P (x)Q I n+1 ) -1 Q t I n+1 P (x)η I c 2 2 + ηI c 2 2 .
We assume Q I n+1 to have a full rank and we denote its the singular value decomposition by

Q I n+1 = U d×p (I n+1 )∆ p×p (I n+1 ) t V p×p (I n+1 ),
where t U (I n+1 )U (I n+1 ) = t V (I n+1 )V (I n+1 ) = I p (the identity matrix of size p). Let us define the matrix P j (x, I n+1 ) = t U (I n+1 )P j (x)U (I n+1 ), we know from [START_REF] Oyet | Exact Wavelet Designs for Discrimination[END_REF] that the least favourable approximated bias is the largest eigenvalue of P -1 1 (x, I n+1 )P 2 (x, I n+1 )P -1 1 (x, I n+1 ) denoted by λ(x, I n+1 ). The function Φ is then approximated by

Φxn,I n+1 (x) = τ 2 λ(x n ∪ x, I n+1 ) + νT r µ 1,1 (I n )M -1 xn∪x,In .
The above Φ is used to select the optimal design at step n + 1 via x n+1 = arg min x Φxn,I n+1 (x).

Stochastic Model Selection (⋄MS⋄ step)

In this section, a stochastic algorithm which updates I n is described. We provide first a method to estimate the "importance" of each point of I n with respect to the bias term in (6) and then we deduce a ranking criterion. This ranking will be plugged into the acceptation/rejection step of a stochastic Simulated Annealing (S.A.). In a sense, this strategy is a kind of forward/backward stepwise model selection strategy. We try to add points which reduce the bias term but do not increase the variance of the estimation. Thus, this idea is similar to the MARS approach in [F91], the adaptive strategy of [START_REF] Castaño | Adaptive fitting of scattered data by spline wavelets Curve and surface fitting[END_REF], or the algorithm of [START_REF] Busby | Hierarchical nonlinear approximation for experimental design and statistical data fitting[END_REF].

Ranking criterion

Recall first that if ηxn,In (x) denotes the estimation of η with model I n and design x n , then B is given by:

B xn,In = E (E[η xn,In (x)] -η(x)) 2 dx.
Let us recall that we can decompose Eη xn,In (x) in the basis

Eη xn,In (x) = (r,t)∈In θ r,t Λ r,t .
Our idea is that the importance of each function Λ r,t for the bias can be deduced from the sensitivity of the bias to a small perturbation of the coefficient θ r,t . This idea is comparable to the Recursive Feature Elimination (RFE) algorithm of [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF]: Λ r,t is important when |∂B xn,In /∂θ r,t | is large. Now, we can formally compute each partial derivative of B xn,In with respect to θ r,t to measure the importance of each Λ r,t as

∂B xn,In ∂θ r,t = 2 E Λ r,t (x) [Eη xn,In (x) -η(x)] dx .
As pointed above, the exact computation of this last term is impossible as it requires to know η. Indeed, one can first estimate the bias pointwise Eη xn,In (ξ)η(ξ) for each point of the design using a standard "leave one out" cross validation coupled with a bootstrap sampling strategy. Then we use a classical method of interpolation to estimate the bias as a function on E. Let us denote ǫ := ηxn,Inη. 3. Estimate the bias with (ǫ(ξ), ξ ∈ x n ) and a simple spline smoothing estimator using a Gaussian kernel whose bandwith parameter is chosen by the Generalized Cross Validation criterion. For more details, one may refer to [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF] for instance. This step provides a bias interpolation ǫ defined on E.

4. Estimate the derivative of the bias by a simple integration

∂ (r,t) B := 2 E Λ r,t (x)ǫ(x)dx .

Stochastic Learning of I n with a Simulated Annealing dynamic

We propose a method to update I n . This algorithm is largely inspired by Metropolis-Hastings methods.

The Simulated Annealing strategy is classically decomposed in a proposition step and an acceptation rule.

Reversible Jump proposal

Recall formally some elements of the S.A. procedure. Let Ω be a measurable set with a measure m and let µ be a measure on Ω with a density (also denoted by µ) with respect to m. We aim to minimize some cost C. The S.A. involves a simulation of a non-homogeneous Markov chain whose invariant distribution at iteration n is µ n ∝ µ -C/Tn where (T n ) n≥0 is a temperature with T n → 0. Under classical conditions (see [START_REF] Hajeck | Cooling schedules for optimal annealing[END_REF][START_REF] Geman | Stochastic Relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] for instance), µ ∞ concentrates on the set of minima of C. The S.A. method with transition distribution q(I, I ′ ) works as follow:

• from state I ∈ Ω, first propose a state I ′ with probability q(I, I ′ )

• then, accept the transition with a probability which is adjusted so that µ is invariant.

We assume the following reversible property q(I, I ′ ) > 0 ⇐⇒ q(I ′ , I) > 0.

The probability to accept the transition I to I ′ at iteration n is then defined as:

∀I ′ = I Q n (I, I ′ ) = µ n (I ′ )q(I ′ , I) µ n (I)q(I, I ′ ) ∧ 1. ( 9
)
When µ corresponds to a Gibbs field associated to a cost function C (J defined in (6) in our case), this ratio is in fact given by

∀I ′ = I Q n (I, I ′ ) = e C(I)-C(I ′ ) Tn q(I ′ , I) q(I, I ′ ) ∧ 1. ( 10
)
The main difficulty is to ensure the weak reversibility condition given in the former paragraph: q(I, I ′ ) > 0 ⇐⇒ q(I ′ , I) > 0.

Birth and Deletion transition

In our framework, we start with {(0, 0); (1, 0); (1, 1)} and we decide to use the following dynamic for the iteration I n → I n+1 :

B: Birth of any element i / ∈ I n associated to a son or a father already in I n .

R: Rebirth of the element Λ 0,0 if (0, 0) / ∈ I n .

D: Deletion of any element i ∈ I n provided that one son or father is still in I n .

Given any iteration n with a design x n and a basis I n , we fix:

• p R the probability to add Λ 0,0 , if Λ 0,0 ∈ I n , we set p R = 0 otherwise p R = 0.1 .

• p B ∈]0; 1[ the probability to add a function to I n .

• p D the probability to delete one element of I n .

We state p B = 5p D and p B + p D + p R = 1.

• In the birth case, denote by I B n the set of elements in I n such that one of their sons or father is not in I n . Then, propose the birth of the ascendant or descendant of some element Λ i , i ∈ I B n where we sample i with a discrete probability r B ∝ ∂B i .

• In the deletion case, denote I D n the set of elements in I n such that one of their descendants or ascendant is in I n and propose the deletion of one element of

I D n following the distribution r D ∝ ∂B -1 i .
The resulting transition kernel of the simulated Markov chain is then a mixture of the different transition kernels associated with the moves described above. We choose now classically

T n = C 1 C 2 +log(n)
and this yields the transition kernel q.

Remark 2 Please remark that the reachable vertices in I n are not a priori connected in the tree representation space, it is a consequence of the reversibility condition. These moves are defined by heuristic considerations, the only condition to be fulfilled is to maintain the correct invariant distribution defined by (10).

These moves are not the classical dynamic of dyadic trees for sets I n . The classical evolution would generate connected trees but such trees are not sparse. At last, the necessary reversible jump condition is fulfilled by the definition of B, R and D.

Variance estimate and summary of the ⋄MS⋄ step

In order to obtain the complete transition kernel defined in (10), it is now sufficient to estimate ∆C which is the difference of the IMSE before and after the proposed transition. It is thus necessary to estimate both the bias and the variance of the model. The quantity ∆B xn,In can be estimated using a bootstrap "leave one out" procedure as described in the paragraph 4.1. To compute ∆V xn,In , one just need to have a good estimation σ of σ defined in (1). This can be done using a maximum likelihood estimator of σ as pointed in [START_REF] Busby | Hierarchical nonlinear approximation for experimental design and statistical data fitting[END_REF]. Note that in the resulting algorithm, σ is just needed to control the acceptation / rejection rule. Thus, the practical transition rule is given by a slightly modified equation (replace σ by its estimate σ in the control term C).

To sum up, the ⋄MS⋄ step can be described as 1. Choose a jump B, D or R according to the distribution (p B , p D , p R ).

2. Compute each ranking coefficient ∂B(i), ∀i ∈ I n .

3. Compute the sampling distribution r B or r D and propose a new state I n+1 according to this discrete distribution.

4. If x n = {x 1 , . . . x n } is the design, compute σ with the M.L.E. of σ given by:

σ2 n = 1 |x n | t f (x n ) -(Λ In (x 1 ), . . . , Λ In (x n )) θn M -1 xn,In f (x n ) -(Λ In (x 1 ), . . . , Λ In (x n )) θn .

Compute the approximated differential cost ∆B (xn,In)->(xn,I n+1

) + σ2 n ∆V (xn,In)->(xn,I n+1 ) and accept the transition with the probability

Q n (I, I ′ ) = e -[∆B+σ 2 ∆V ] ∧ 1.
In the acceptation case, update I n+1 with the proposed transition. Otherwise keep I n unchanged:

I n+1 = I n .
Since f is decomposed on I n at step n, the M.L.E. σn is in general larger than the true σ. Hence, the acceptation/rejection rule Q n prevents to overfit the models I n .

Experimental results

In the sequel, we call a deterministic equispaced design a regular design whereas uniformly sampled designs will denote random designs obtained with a random uniform sampling over E. For all experiments, we set E = [0; 1] and we aim to compare:

1. Least square linear model with • Optimal design coupled with model selection for the Schauder basis.

• Optimal design coupled with model selection for the Haar basis.

• Optimal design coupled with model selection for the Meyer basis.

• Optimal minimax design coupled with model selection for the Schauder basis.

• Optimal minimax design coupled with model selection for the Meyer basis.

2. Sparse Penalized regressor Adaptive Lasso (Ada-Lasso) [Z06] solved by the LARS implementation [START_REF] Efron | Least angle regression[END_REF] (Schauder, Haar or Meyer basis).

• Ada-Lasso regressor with design uniformly sampled in E.

• Ada-Lasso regressor with regular design.

Wavelet approach of [AAP06] (Daubechies (D6) and Symmlets (S6) basis):

• Kernel penalized estimation with uniformly sampled design.

• Kernel penalized estimation with regular sampled design .

This section presents three examples where the unknown signal η must be recovered from as few observations as possible. The first synthetic example deals with the approximation of some unknown function that does not belong to finite dimensional vector spaces spanned by the triangle Schauder basis. The second example illustrates the database of Motorcycle impact experiment ( [S85]). We provide also few results in a 2-dimensional setting. We consider the evolution of the IMSE with respect to the number of evaluations of the signal. Each implementation of [START_REF] Efron | Least angle regression[END_REF] and [START_REF] Amato | Wavelet kernel penalized estimation for non-equispaced design regression[END_REF] will be used with uniformly sampled design and regular design. Note that both the Adaptive Lasso implementation and the Wavelet Penalized Kernel will be optimized using a cross-validation criterion. For the Adaptive Lasso, the bias of the penalization is removed once the model is determined using a least square estimation procedure based on the selected model. At last, the calibration of ν for minimax design is important. This parameter must be chosen sufficiently large to include the effect of the variance in our model. We have set ν = 20 in our experiments, which seems a reasonnable value according to the size of the bias and variance in our experiments.

Description of the data

We investigate first the efficiency of estimation of the methods described above. We initialize the multiresolution basis functions I 0 as ΛI 0 = {Λ 0,0 ; Λ 1,0 ; Λ 1;1 }.

Definition of the synthetic η

The synthetic function η to recover is a mixture of localized Gaussian kernel with different scaling parameters. For this, we set η to be localized around some values in E with different amplitudes and frequencies:

∀x ∈ [0; 1] η(x) = 5e -1000(x-0.25) 2 + 5e -100(x-0.75) 2 + 20e -100(x-0.5) 2 .

The signal f is thus given by f (x) = η(x) + σζ(x) and some realizations of f are shown on Figure 3.

Figure 3: Synthetic function η and some realizations of f (x) with σ = 0.5 (crossed points) or σ = 2 (dashed points).

Motorcycle impact experiment

In the Motorcycle Impact Experiment (see [S85] for a brief description of the data), the efficiency of crash helmets it studied. In [S85], the author uses 133 observations and a spline smoothing approach to estimate a curve deduced from the discrete observations. We scale the 133 observations between 0 and 1 and we compute a kernel smoothing interpolation described in [S85] to have a signal η to compare with our estimator. At last, we randomize the kernel interpolation by the addition of a white noise which yields an homoscedastic regression problem.

Remark 3 Indeed, the Motorcycle experiment belongs to the more general class of problems of heteroscedastic regressions. Note that we use this dataset as a benchmark to validate our algorithm.

Results

In order to obtain a reliable estimation of the IMSE for each algorithm, we repeat our experiments 50 times for each different method. Each experiment is reproduced also setting σ = 0.5 (low noise) and then σ = 2 (high noise). The Figure 4 shows the performance of our non-minimax sequential algorithm described in the beginning of this section on the synthetic dataset. One can see that in the low noise setting, there is no clear advantage to use Meyer basis which performs as Haar and Schauder bases. In contrast, in the high noise setting, the Meyer wavelet basis provides the best results. One can infer from this simulation that the linear model built with the Meyer wavelets is much more stable than the other ones. This is certainly due to the fact that Haar and Schauder mother functions are compactly supported and the information matrix of the least square estimator is "more singular" than the one constructed with the Meyer basis. Moreover, Figure 11 provides the detailed numerical list of IMSE performance for each of the methods listed above in the case σ = 0.5 and Figure 12 provides the same results for the large noise case σ = 2 for the synthetic dataset. Best results for each number of points in the designs are underlined in the tables.

Figure 4: Integrated Mean Square Error for the synthetic example with low (dashed curves) and high noise (continuous curves) for the Haar (red), Schauder (blue) and Meyer (green) multi-resolution families. The x-axis represents the number of iterations of our method (number of points in the optimal designs).

These numerical results confirm that the use of the Meyer wavelet basis is always better than the other multi-resolution setting in our sequential strategy and this phenomena is much more important when the variance is important. The other striking fact is the poor performance of the Adaptive Lasso with the several multi-resolution families we used. This point is mainly due to the limited size of the set of experiments and the Adaptive Lasso is completely inapropriate in this framework. Moreover, in some rare cases (depending on the basis used), some very conservative choice are made by Ada-Lasso (see for instance the case of the Schauder basis in the high noise setting).

One can also remark the very good behaviour of the thresholding procedure as soon as there are enough points in the design. Since such methods are specifically dedicated to the multi-resolution framework, it is not surprising that they behave better than Adaptive Lasso in these cases. However, note that in a very few cases of measurements, a simple linear model with the design learned by our sequential method is more efficient than non-linear thresholding procedures.

At last, the minimax sequential design strategy performs the best. But the method is costly compared to the non-minimax Schauder one and there is no clear advantage to use it for such small IMSE in the small noise setting. For the large noise case, the minimax sequential design is justified since it improves the IMSE of 50 % compared to the non-minimax Schauder sequential approach.

The Figure 8 shows the performance of the sequential algorithm for the Motorcycle dataset. We display some estimation results obtained with the Schauder basis (Figure 9) and the Meyer wavelet basis (Figure 10). Note that the conclusions for the synthetic data still hold for the Motorcycle dataset since the Meyer basis is clearly better than the other ones as pointed in the Figure 8. For the sake of simplicity, we have omitted the results of the Adaptive Lasso algorithm on Figures 11 and12. One may remark that with a very small number of points in the design, the sequential method always outperforms the wavelet thresholding procedure. But for n = 50, thresholding with Daubechies or Symlets bases are equivalent to our procedure with non-minimax design. One can remark that the minimax sequential design achieves again the best performances. This is especially the case for a very small number of experiments (n = 10): the minimax criterion tends to scatter the points of the design among E at the beginning of our strategy. It is explained by the fact that the variance term is smaller than the worst bias for small number of experiments and low resolution elements in I n . This phenomenon is reversed when the number of iterations is growing.

Figure 7 shows the good behaviour of our estimation of the bias as described in the paragraph 4.1 with a few number of functions in I n . Note that this estimation deteriorates a little bit when the number of functions in I n is growing. This last point is not very important since we mainly use this estimation as an indicator in the proposition step of our stochastic algorithm.

Figure 7: Estimation of the bias by our bootstraped procedure with 10 experiments and 5 functions in the family I on the synthetic example.

Computation Time

At last, concerning the computational costs, since the sequential design for Meyer basis is obtained by a simulated annealing algorithm, it may be interesting to switch to other multi resolution bases such as the Schauder one. Indeed, an iteration n → n + 1 with the Schauder basis requires less than 5 seconds although the same iteration with the Meyer wavelet basis needs for large n approximately 1 minute.

Although both the design and the model are estimated with our algorithm, it is as fast as the penalized approach or a thresholding method. However, the application of our method with the Meyer basis is much more costly. The minimax approach is obviously more costly than the non-minimax one when the set of functions is the Schauder one. But the cost of the minimax approach is equivalent to the one of the non-minimax approach with Meyer basis.

2-Dimensional experiments

We provide a short illustrating example in 2 dimensions with the Schauder triangle basis. Let us suppose

∀(x, y) ∈ [0; 1] 2 η(x, y) = 10xy(x -1)(y -1)e 5(x-3/4) 2 +(y-3/4) 2 -5(x-1/4) 2 +(y-1/4) 2 . ( 11
)
We present in Figure 13 the estimation of η with the tensorized Schauder basis with a small number of iterations. Such basis are deduced from the bases in the 1 dimensional case as follows:

F n = {(x, y) → Λ r 1 ,t 1 (x); (x, y) → Λ r 2 ,t 2 (y); (x, y) → Λ r 1 ,t 1 (x)Λ r 2 ,t 2 (y)} ,
where (r 1 , t 1 ) and (r 2 , t 2 ) belong to a finite set I n . Roughly speaking, one can adapt the conclusions of Theorem 1 to obtain some properties for the location of the designs. The simple adaptation of the ⋄OD⋄ step consists in:

1. Compute all singular points S x (I n ) of elements in Λ (r,t) on the x-coordinate.

2. Compute all singular points S y (I n ) of elements in Λ (r,t) on the y-coordinate. We show in Figure 13 some estimations obtained with a very small number of iterations of our algorithms. The red surface is the true response and the meshed surface is our estimation. One can remark that the algorithm still takes into account the singularities of the signal as it was the case in the one dimensional examples. Moreover, the approximation of the extrema is pretty good and this example shows that our method is promising for response surface applications.

Then S(I

n ) = S x (I n ) × S y (I n ).

Conclusion

The adaptive method developed in this paper is numerically competitive for synthetic and real examples compared to thresholding wavelet or ℓ 1 penalized methods. The iterative scheme is fast and may be very fast if it does not require any complicated optimization step as it is the case in the special case of the Schauder basis. The approximation properties of both the global signal or its maxima are satisfying. At last, the model selection ability could be of great interest for variable selection motivations and are meaningful. High resolution functions are added when needed or not used when the signal is well approximated.

But, on the theoretical side, many questions remain open. First, it would be very fruitful to generalize the result that localizes the optimal designs on dyadic points for other multi-resolution sets of smooth functions. But it is a difficult task owing to the underlying non-linear nature of the optimization problem in the case of general wavelet bases.

Second, the proof of the convergence of the stochastic coupled algorithm on (I n , x n ) remains a difficult task. To do so, it is necessary to fix a precise cooling strategy and to use the consistency result of theorem 4 provided in the appendix as a first step.

One may also consider a modified energy criterion based on a discrepancy term to bound the bias B which is unknown in our sequential framework. This would certainly enable some theoretical extensions of the consistency result (Theorem 4) in the fixed basis setting At last, it would be very fruitful to infer a sequential optimal design approach for the ℓ 1 penalized approaches. The optimization of the variance with respect to the design seems difficult, mostly because the choice of the penalty parameter as discussed in paragraph 3.4 is not explicit.

Aknowledgements We would like to gratefully acknowledge Jean-Marc Azaïs, Fabrice Gamboa and Anestis Antoniadis for their helpful comments on an early version of the manuscript and their warm support. We are very much indebted to the anonymous referees and the associate editor for their constructive comments that resulted in a major revision of the original manuscript.

Appendix

Location properties of the design

We will denote by n the number of points in a fixed design x and by p the cardinal of I. Recall that ΛI is the rectangular matrix defined by equation (5). Since I is fixed, we drop the indice I in this appendix for the notation of M x = t ΛI ΛI . We suppose that n + 1 ≥ p (which is rather trivial in our context) and denote by F the map given by F (x) = det (M x∪x ). We will show that F is a convex map on every interval where it is differentiable. Assuming x to be suitably chosen among differentiable points of ΛI , we will note Λ′ I (x) the vector composed of the differentiable maps of ΛI computed at point x and the matrix

M ′ x = ( ΛI 1 Λ′ I 2 + Λ′ I 1 ΛI 2 )(x) i 1 ,i 2 ∈I = d dx (M x∪x ) .
Using the standard Euclidean scalar product on R p , one can check immediately that

∀U ∈ R p , M ′ x U = ΛI (x); U Λ ′ I (x) + Λ′ I (x); U ΛI (x).
First, we state some classical results on matrices whose proofs are based on standard arguments on matrices of rank 1. Some details can be found in [START_REF] Meyer | The Coordinate-Exchange Algorithm for Constructing Exact Optimal Experimental Designs[END_REF] and in chapter one of [F69].

Proposition 1 If M -1
x∪x and M -1

x are non-singular,

M -1 x = M -1 x∪x + M -1 x∪x ΛI (x) t ΛI (x)M -1 x∪x 1 -t ΛI (x)M -1 x∪x ΛI (x) (12)
M -1 x∪x = M -1 x - M -1 x ΛI (x) t ΛI (x)M -1 x 1 + t ΛI (x)M -1 x ΛI (x) . ( 13
det M x = 1 1 -t ΛI M x∪x ΛI , ( 14
) det M x det M x∪x = 1 1 + t ΛI M x ΛI . ( 15
)
We now establish two technical lemma useful for our location theorem.

Lemma 1 For any symmetric matrix S, we have the relation

T r(M ′ x S) = T r(SM ′ x ) = 2 S ΛI (x); Λ ′ I (x) . (16)
Proof: Consider first the case where Λ I (x), Λ ′ I (x) are linearly independent in R p . One can show that

M ′ x S ΛI (x) = S ΛI (x); ΛI (x) Λ ′ I (x) + SΛ ′ I (x); ΛI (x) ΛI (x), and
M ′ x SΛ ′ I (x) = S ΛI (x); Λ ′ I (x) Λ ′ I (x) + SΛ ′ I (x); Λ ′ I (x) ΛI (x). Since the rank of M ′
x is 2, we can find a basis adapted to the family ( ΛI (x); Λ ′ I (x)) such that the endomorphism described by M ′

x S in the basis is  Suppose now that Λ I (x), Λ ′ I (x) are linearly dependent, we get S ΛI (x); ΛI (x) Λ ′ I (x) = SΛ ′ I (x); ΛI (x) ΛI (x), and applying the same argument as above with the endomorphism M ′ x S the rank of which is one in this case, we also obtain T r(M ′ x S) = 2 S ΛI (x); Λ ′ I (x) . If we denote by Com(M ) the matrix t cof (M ), where cof (M ) is the matrix of cofactors of A, we have the following result.

       SΛ ′ I (
Lemma 2 Assume x to be a regular point for the map ΛI , and that M x , M x∪x are non-singular, then

T r t Com(M x∪x )M ′ x = T r t Com(M x )M ′ x . Proof: Apply lemma 1 first to S = M -1 x∪x , we get T r(M -1 x∪x M ′ x ) = 2 M -1 x∪x ΛI ; Λ′ I . (17)
Moreover, lemma 1 applied now to

S = M -1 x∪x Λt I ΛI M -1 x∪x yields T r M -1 x∪x Λt I ΛI M -1 x∪x M ′ x = 2 M -1 x∪x ΛI t ΛI M -1 x∪x ΛI = ΛI ;M -1 x∪x ΛI ; Λ′ I . Thus T r M -1 x∪x Λt I ΛI M -1 x∪x M ′ x = 2 ΛI ; M -1 x∪x ΛI Λ ′ I ; M -1 x∪x ΛI . (18)
From ( 12),(17), and (18), we get

T r(M -1 x M ′ x ) = 2 M -1 x∪x ΛI ; Λ′ I + 2 ΛI ; M -1 x∪x ΛI Λ ′ I ; M -1 x∪x ΛI 1 -t ΛI M -1 x∪x ΛI = 2 M -1 x∪x ΛI ; Λ′ I × 1 + ΛI ; M -1 x∪x ΛI 1 -t ΛI M -1 x∪x ΛI T r(M -1 x M ′ x ) = T r(M -1 x∪x M ′ x ) 1 -t ΛI M -1 x∪x ΛI
.

Now, use (13) and the relation

A -1 = t (Com(A)) det(A)
to reach the conclusion of the lemma:

T r t Com(M x∪x )M ′ x = T r t Com(M x )M ′ x .
Proof of theorem 1: We will note F (x) = det (M x∪x ). Suppose first that M x is non-singular and x is not a dyadic point described by the set S . In this case, classical differentiation used with lemma 2 yields

F ′ (x) = T r t Com(M x∪x )M ′ x = T r t Com(M x )M ′ x .
Finally, since T r is a linear map, we immediately get

F "(x) = T r t Com(M x )M " x = T r t Com(M x ) Λ′ I t Λ′ I = t Λ′ I t Com(M x ) Λ′ I ≥ 0.
Thus F is a convex function on each interval outside of S Consequently, its maximum are located on some dyadic points of S. This is equivalent to the assertion of the proposition. Suppose now M x is singular, we can find a sequence M x,ǫn = M x + ǫ n Id which is non-singular such that lim

n →+∞ M x + ǫ n Id = M x .
Consider now the function F ǫn (x) defined as

F ǫn (x) = det (M x,ǫn∪x ) .
We can use the same arguments as before to conclude that arg max F ǫn ⊂ S since these arguments only rely on a slight modification of lemma 2 which becomes:

F ′ ǫn (x) = T r t Com(M x,ǫn∪x )M ′ x = T r t Com(M x,ǫn )M ′ x .
Now, remark that E is a finite set which is not varying with ǫ n and ∀x F ǫn (x) ≤ max x∈E F ǫn (x).

Taking the limit in the relation above yields the conclusion of the proof. Proof of theorem 3: Remark first that t → det tId + M -1 x∪x is a polynomial function of t whose degree p is the size of ΛI . This polynomial function is expanded in

det tId + M -1 x∪x = t p -T r M -1 x∪x 1 t p-1 + Q x (t)
where deg(Q x ) ≤ p -2. Now for x 1 , x 2 ∈ E such that

T r M -1 x∪x 1 ≥ T r M -1 x∪x 2 ,
we can immediately check that for sufficiently large t, we have

det tId + M -1 x∪x 1 ≤ det tId + M -1 x∪x 2 .
Consequently, the solutions of the trace maximization problem are the same as the one deduced from the determinant minimization problem and this remark ends the proof.

Remark 4 To extend now the proof to dimensions d higher than one with some tensorized family of Schauder functions, one just have to remark that both lemma 1 and 2 are still valid for x ∈ R d . Then a similar argument to the one used in the proof of theorem 1 shows the convexity of F except in the neighbourhood of points of the form (2 -j 1 k 1 , . . . , 2 -j d k d ).

Convergence in the case of fixed basis

We detail here the convergence of the estimation of the parameter θ with the strategy of sequential optimal design when the basis I remains fixed . As both previous criteria yield the same optimal design, we are only concerned with the study of the sequential strategy:

x n+1 = arg max

x det(M xn∪x ) and x n+1 = x n ∪ x n+1 , while θ is classically given by θn = M -1 xn ΛI (x n )f (x n ).

Theorem 4 Let f and η be defined by equations (2) and (3) with a fixed basis I, and suppose that η ∈ Span(I). Then the sequential optimal design is consistent: θn → θ a.s. Moreover, there exists a positive constant C such that θnθ ∞ ≤ C log n n .

Remark 5 The previous theorem ensures the consistency of θn if the signal η is a linear combination of the functions (Λ i ) i∈I . Note that when η / ∈ Span((Λ i ) i∈I ), the convergence to the orthogonal projection of η into Span((Λ i ) i∈I ) also holds. where λ min (M ) denotes the minimum eigenvalue of M and λ max (M ) the maximum eigenvalue of M . We establish first the condition C1. Remark that as the functions (Λ i ) i∈I are linearly independent, we can find ρ > 0 such that B(0, ρ) ⊂ Conv( ΛI (t), t ∈ E) ∪ -Conv( ΛI (t), t ∈ E),

where Conv denotes the convex hull of a set. Now, we have for any symmetric positive definite M max y∈B(0,ρ) t yM -1 y = λ min (M ) -1 ρ 2 , and since y → t yM -1 y is convex, we can state that

max x∈E t ΛI (x)M -1 ΛI (x) ≥ ρ 2 λ min (M ) . ( 19
)
Remark that all maps in ΛI are continuous and E is compact, thus

∃ L > 0 ∀t ∈ E t ΛI (t) ΛI (t) 2 ≤ L,
where A 2 := sup x∈B(0,1) Ax , where we take the Euclidean norm in the last definition. Now, the spectral radius satisfies the triangular inequality and

λ max M x k k ≤ k i=1 λ max ( ΛI (x i ) t ΛI (x i )) k ≤ L. If I k = M x k /k, the last inequality yields λ max (I k ) ≤ L. (20)
Next define ρ k = det(I k ) and d k (t) = t ΛI (t)I -1 k ΛI (t), from Proposition 1 equation (15), we have

ρ k+1 = k k + 1 p 1 + d k (x k+1 ) k ρ k ≥ ρ k k k + 1 p .
Thus, for any ǫ > 0, we can find K 1 ≥ 1 such that

∀k ≥ K 1 ρ k+1 ≥ (1 -ǫ)ρ k , (21)
and a simple induction shows that

ρ k ≥ (1 -ǫ) k-K 1 ρ K 1 . Let A k = (1 -ǫ) k-K 1 ρ K 1 , since A k → 0 as k → ∞, we can find K 2 ≥ K 1 such that ∀k ≥ K 2 ρ 2 A 1/p k > 2p and k + 1 k p ≤ 1 + 2p k . (22)
We show now by induction that ρ k is bounded from below by (1ǫ)A K 2 for sufficiently big k. This is obviously true for k = K 2 + 1.

Suppose now that ρ k ≥ (1ǫ)A K 2 . If ρ k ≥ A K 2 , in view of (21) we immediately obtain ρ k+1 ≥ (1ǫ)A K 2 . We must thus study the case Finally, the definition of x k+1 yields

A K 2 > ρ k ≥ (1 -ǫ)A K 2 .
ρ k+1 = ρ k k k + 1 p 1 + d k (x k+1 ) k = ρ k k k + 1 p 1 + max x∈E d k (x) k ≥ ρ k .
This last inequality concludes the induction and ρ k is bounded from below by a constant Γ. Now, remark that λ min (I k )λ max (I k ) p-1 ≥ det(I k ) ≥ Γ, and we obtain from equation (20) as k → +∞:

λ min (M x k ) ≥ k Γ L p-1 → +∞.
This last equation proves condition (C1).

For (C2), simple algebra yields as k → +∞:

λ min (M x k ) log (λ max (M x k )) ≥ kΓ L p-1 log (kL) → ∞,

Figure 1 :

 1 Figure 1: Several functions Λ r,t , here (r, t) belongs to {(1, 1)(3, 1)(4, 6)}.

Figure 2 :

 2 Figure 2: Evolution with respect to x of the variance Ṽ1 when I = {(0, 0); (1, 0); (1, 1); (2, 3)} for the Schauder triangle basis

1.

  Fix the number of bootstrap samplers b (for instance b = 10). 2. For each point of the design x k : (a) Build b bootstrap samples without the observation x k : (D 1 , . . . D b ). (b) Compute for each bootstrap sample the estimators ηD 1 ,In , . . . , ηD b ,In . (c) Compute the "leave one out" bootstrap bias at x k by ǫ(x k ) = 1 b b m=1 ηDm,In (x k )ηxn,In (x k ).

Figure 8 :

 8 Figure8: Integrated Mean Square Error of our sequential methods for the Motorcycle example with low (dashed curves) and high noise (continuous curves) for the Haar (red), Schauder (blue) and Meyer (green) multi-resolution families. The x-axis represents the number of iterations of our method (number of points in the optimal designs).

Figure 9 :

 9 Figure 9: Estimations on the Motorcycle example with high noise (σ = 10) at iteration 0 (a), 10 (b), 20 (c), 30 (d), 40 (e), 50 (f). Continuous curve: true signal, Dashed curve: interpolation of our sequential procedure with the Schauder wavelet basis.

Figure 10 :

 10 Figure 10: Estimations on the Motorcycle example with high noise (σ = 10) at iteration 0 (a), 10 (b), 20 (c), 30 (d), 40 (e), 50 (f). Continuous curve: true signal, Dashed curve: interpolation of our sequential procedure with the Meyer wavelet basis.

.

  x); ΛI (x) S ΛI (x); ΛI (x) 0 . . .Thus in this caseT r(M ′ x S) = 2 S ΛI (x); Λ ′ I (x) .

  Proof of theorem 4: This proof is inspired by[P00] which states the almost sure convergence of θn to θ provided the two conditionsC1 λ min [M xn ] → ∞ a.s. C2 log (λ max [M xn ]) = o (λ min [M xn ]) a.s.

.A

  From the definition of d k and (19), we have max x∈Ed k (x) ≥ k ρ 2 λ min (M x k ) ≥ k ρ 2 det(M x k ) 1/p ≥From equation (22) and our assumption on ρ k , we obtain max x∈E d k (x) ≥ ρ 2

Figure 13 :

 13 Figure 13: Estimations of the 2 dimensional signal with our adaptive sequential design method for 1, 6 (left and right top), 12, 18 (left and right middle), 24 and 30 (left and right bottom) iterations.

  Integrated Mean Square Error for the Motorcycle experiment and high noise.

	Method		IMSE (n=10) IMSE (n=30) IMSE (n=50)
	Sequential Haar		296	91	50
	Sequential Schauder		41.9	7.5	6.5
	Sequential Meyer		19.7	7.4	6.0
	Sequential minimax Schauder	23	9.5	5.5
	Sequential minimax Meyer		15.7	5.2	5.1
	Wavelet Kernel Penalized D6 Random 1549	26.4	9
	Wavelet Kernel Penalized D6 Regular 458	15	12
	Wavelet Kernel Penalized S6 Random 188	154	8.9
	Wavelet Kernel Penalized S6 Regular	28.4	11.3	9.5
	Figure 11: Integrated Mean Square Error for the Motorcycle experiment and low noise.
	Method		IMSE (n=10) IMSE (n=30) IMSE (n=50)
	Sequential Haar		477	239	232
	Sequential Schauder		171	153	152
	Sequential Meyer		104	28.5	18.3
	Sequential minimax Schauder	94	23	14
	Sequential minimax Meyer		88	21	12
	Wavelet Kernel Penalized D6 Random 1074	158	93
	Wavelet Kernel Penalized D6 Regular 556	115.7	135
	Wavelet Kernel Penalized S6 Random 180	129	30
	Wavelet Kernel Penalized S6 Regular	122	59	18
	Figure 12: Moreover,	det M x∪x
					)