MAP5 F Comte
email: comte@biomedicale.univ-paris5.fr.

Y Rozenholc

A NEW ALGORITHM FOR FIXED DESIGN REGRESSION AND DENOISING

In this paper, we present a new algorithm to estimate a regression function in a fixed design regression model, by piecewise (standard and trigonometric) polynomials computed with an automatic choice of the knots of the subdivision and of the degrees of the polynomials on each sub-interval. First we give the theoretical background underlying the method: the theoretical performances of our penalized least-squares estimator are based on non-asymptotic evaluations of a mean-square type risk. Then we explain how the algorithm is built and possibly accelerated (to face the case when the number of observations is great), how the penalty term is chosen and why it contains some constants requiring an empirical calibration. Lastly, a comparison with some well-known or recent wavelet methods is made: this brings out that our algorithm behaves in a very competitive way in term of denoising and of compression.

Introduction

We consider in this paper the problem of estimating an unknown function f from [0, 1] into IR when we observe the sequence Y i , i = 1, . . . , n, satisfying [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] Y i = f (x i ) + σε i , for fixed x i , i = 1, . . . , n in [0, 1] with 0 ≤ x 1 < x 2 < • • • < x n ≤ 1. Most of the theoretical part of the work concerns any type of design but only the equispaced design x i = i/n is computationally considered and implemented. Here ε i , 1 ≤ i ≤ n is a sequence of independent and identically distributed random variables with mean 0 and variance 1.

The positive constant σ is first assumed to be known. Extensions to the case where it is unknown are proposed. We aim at estimating the function f with a data driven procedure. In fact, we want to estimate f by piecewise standard and trigonometric polynomials in a spirit analogous but more general than e.g. [START_REF] Denison | Automatic Bayesian curve[END_REF]. We also want to choose among "all possible subsets of a large collection of pre-specified candidates knot sites" as well as among various degrees on each subinterval defined by two consecutive knots.

Our method is based on recent theoretical results obtained by [START_REF] Baraud | Model selection for regression on a fixed design[END_REF][START_REF] Baraud | Model selection for regression on a random design[END_REF], Baraud et al. (2001a, b) who adapted to the regression problem general methods of model selection and adaptive estimation initiated by [START_REF] Barron | Minimum complexity density estimation[END_REF] and developed by [START_REF] Birgé | Minimum contrast estimators on sieves: exponential bounds and rates of convergence[END_REF], [START_REF] Barron | Risks bounds for model selection via penalization[END_REF]. Results on Gaussian regression can also be found in [START_REF] Birgé | Gaussian model selection[END_REF]. It is worth mentioning that a similar (theoretical) solution to our regression problem, in a context of regression with random design, is studied by [START_REF] Kohler | Nonparametric estimation of piecewise smooth regression functions[END_REF]: he proposes also piecewise smooth regression functions to estimate the regression function, and he uses a penalized least-squares criterion as well. The approach is similar to [START_REF] Baraud | Model selection for regression on a fixed design[END_REF] and he uses Vapnis-Chervonenkis theory in place of Talagrand's or deviation inequalities.

All the results we have in mind about fixed design regression have the specificity of giving non asymptotic risk bounds and of dealing with adaptive estimators. The first results about adaptation in the minimax sense in that context were given by [START_REF] Efromovich | Learning algorithm for nonparametric filtering[END_REF]. Some asymptotic results have been also proved by [START_REF] Shibata | An optimal selection of regression variables[END_REF], [START_REF] Li | Asymptotic optimality for c p , c l cross-validation and generalized cross-validation: discrete index set[END_REF], [START_REF] Polyak | Asymptotic normality of the c p test for the orthogonal series estimation of regression[END_REF]. An overview of most nonparametric techniques is also given by [START_REF] Hastie | Generalized additive models[END_REF]. Lastly, it is worth mentioning that most available algorithms deal with equally spaced design; results and proposals concerning the more general case of a non necessarily equi-spaced design are quite recent. Some of them can be found for instance in [START_REF] Antoniadis | Wavelet regression for random or irregular design[END_REF], see also the survey in [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF].

An attractive feature of the method which is developed here is that, once a calibration step is done, everything is completely automatic and reasonably fast. [START_REF] Friedman | Flexible parsimonious smoothing and additive modeling[END_REF] already gave an algorithm for optimizing the number and location of the knots of the partition in an adaptive way: this algorithm is used by [START_REF] Denison | Automatic Bayesian curve[END_REF] but the latter calibrate a piecewise cubic fit. In other words, all their polynomials have the same degree a priori fixed to be 3. Ours have variable degrees between 0 and r max (which is r max = 75 in most experiments) and those degrees are also automatically chosen. This provides an important flexibility for denoising step signals for instance. Note that we did not deal with splines (and associated constraints) which were studied for instance by Lindstrom (1999) using a penalized estimation method as well. Moreover, the calibration operation being done once and for all, the only input of the algorithm is the maximal number of knots to be considered and the maximal degree r max . Contrary to many MCMC methods, we do not have any complicated nor arbitrary stopping criterion to deal with, nor do we have any problems of initialization either. A great number of wavelet methods have also been recently proposed in the literature. For an exhaustive presentation and test of these methods, the reader is referred to [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF]. Therefore, we compared our method with standard toolboxes implemented by [START_REF] Donoho | Ideal space adaptation by wavelet shrinkage[END_REF], as well as with some more recent methods tested in [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF]. The performances of our algorithm prove that our method is very good, for any sample size, any type of function f , and whether σ is known or not. Let us mention also that our method seems to globally behave in a very competitive way, in term of L 2 -error performances as well as in term of sparseness of the representation of the signal. In addition, we deal with much more general frameworks. Our main drawback until now is in term of the complexity of our algorithm,which is of the order O(n 2 ) linear operations or O(n 3 ) elementary operations (+, ×, <) when the one of wavelet algorithms is of the order O(n log 2 (n)) elementary operations. Actually we propose a quick but approximated version with complexity of the order O(n) linear operations or O(n 2 ) elementary operations (+, ×, <). As a counterpart the analysis provided by wavelet methods includes 2 n/2 bases which are constructed on dyadic partitions whereas ours includes about (2r max ) n different bases which are constructed on general partitions.

Section 2 gives some more details on the theoretical part of the procedure. We present first its formal principle and a theoretical result is stated. Finally, the general form of the penalty we are working with is written. In Section 3, details about how the estimate is computed are given, two relevant bases are described (one for the space of standard polynomials, the other for the space of trigonometric polynomials) and the reason for the choice of the form of the penalty term involved in the computation of the estimate is explained. Section 4 presents the algorithm: the two main ideas, namely localization and dynamical programming are developed. The scheme for accelerating the algorithm without losing its good properties is introduced. Section 5 presents the empirical results for both the complete algorithm and the accelerated algorithm. The calibration procedure is led by the complete algorithm. Then both methods are compared (in term of L 2error and of compression performances) with wavelet denoising developed by [START_REF] Donoho | Ideal space adaptation by wavelet shrinkage[END_REF] and [START_REF] Donoho | Wavelet shrinkage: asymptopia (with discussion)[END_REF] whose toolbox is available on the Internet together with the test functions we used. Comparison results with 8 other recent methods are also provided.

The general method

2.1. General framework. We aim at estimating the function f of model (1) using a data driven procedure. For that purpose we consider families of linear spaces generated by piecewise (standard or trigonometric) polynomials bases and we compute for each space (basis) the associated least-squares estimator. Our procedure chooses among the resulting collection of estimators the "best" one, in a sense that will be later specified. The procedure is the following. Let D max and r max be two fixed integers and D an integer such that 0 ≤ D ≤ D max . For any D, we choose a partition of [0, 1], that is a sequence Given some m in M n , we define the standard least-squares estimator fm of f in S m by

a 1 , . . . , a D-1 of D -1 real numbers in [0, 1] such that 0 = a 0 < a 1 < • • • < a D-1 < a D =
(3) n i=1 (Y i -fm (x i )) 2 = min g∈Sm n i=1 (Y i -g(x i )) 2 .
In other words, we compute the minimizer fm for all g in S m of the contrast γ(g) where ( 4)

γ(g) = 1 n n i=1 [Y i -g(x i )] 2 .
Each model m being associated with an estimator fm , we have a collection of estimators { fm , m ∈ M n } and we look for a data driven procedure m = m(Y i , i = 1, . . . , n) which selects automatically among the set of estimators the one that is defined as the estimator of

f : f = f m.
m is a vector "number of bins, partition, degree of the polynomials on each piece, type of polynomial on each piece" with values in M n , ( D, ( â1 , . . . , â D-1 ), (r 1 , . . . , r D), ( B1 , . . . , B D)) based solely on the data and not on any a priori assumption on f . Let us precise what the "best" estimator is and how to select it. We measure the risk of an estimator via the Expected Average Squared Error. Namely, if f is some estimator of f , the risk of f is defined by

d 2 n (f, f ) = E 1 n n i=1 (f (x i ) -f (x i )) 2 .
The risk of fm , where fm is an estimator built as in relation [START_REF] Antoniadis | Regularization of wavelets approximations (with discussion)[END_REF], can in fact be proved to be equal (see equation [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF] in [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]) to

d 2 n (f, S m ) + dim(S m ) n σ 2
where d n (f, S m ) = inf t∈Sm d n (f, t) and dim(S m ) denotes the dimension of S m . Therefore an ideal selection procedure choosing m should look for an optimal trade-off between d 2 n (f, S m ), the so-called bias term and σ 2 dim(S m )/n, the so-called variance term. In other words, we look for a model selection procedure m such that the risk of the resulting estimator f m is almost as good as the risk of the best least-squares estimator in the family. More precisely, our aim is to find m such that (5)

d 2 n (f, f m) ≤ C min m∈Mn d 2 n (f, S m ) + σ 2 L m dim(S m ) n ,
where the L m 's are some weights related to the collection of models {S m , m ∈ M n }. This inequality means that, up to a constant C (which has to be not too far from one for the result to be of some interest) our procedure chooses an optimal model and inside that model an optimal estimator in the sense that it realizes a L m -trade-off between the bias and the variance terms. We consider the selection procedure based on a penalized criterion of the following form

m = arg min m∈Mn [γ( fm ) + pen n (m) n ]
where pen n (m) is a penalty function mapping M n into R + . We will precise this penalty later on and just mention that it is closely related to the classical C p criterion of [START_REF] Mallows | Some comments on C p[END_REF].

The procedure is then as follows: for each model m we compute the normalized residual sum of squares, γ( fm ), where γ is defined by (4), we choose m in order to minimize among all models m ∈ M n the penalized residual sum of squares γ( fm ) + pen n (m)/n and we compute the resulting estimator, f m. Mallows' C p criterion corresponds to pen n (m) = 2σ 2 dim(S m )/n where σ2 denotes a suitable estimator of the unknown variance of the ε i 's. Our penalty term is similar but uses an unknown universal constant instead of 2 and the factor L m allowing for very rich collections of models (see the further discussion on the choice of the L m 's). When σ 2 is unknown we replace it by an estimator.

Theoretical results.

From the theoretical point of view, [START_REF] Kohler | Nonparametric estimation of piecewise smooth regression functions[END_REF], [START_REF] Baraud | Model selection for regression on a fixed design[END_REF][START_REF] Baraud | Model selection for regression on a random design[END_REF], Baraud et al. (2001a, b) obtained several results. We formulate in detail a result corresponding to model (1) satisfying the following condition:

(H ǫ ) The ε i 's are i.i.d.centered variables and satisfy, ∀u ∈ R E(exp uε 1 ) ≤ exp (u 2 s 2 /2) for some positive s. This assumption allows the variables ε i 's to be Gaussian with variance s 2 or to be bounded by s. The particular case of Gaussian variables is given in [START_REF] Baraud | Model selection for regression on a random design[END_REF] and the following result is a simplified version of Theorem 1 in Baraud et al. (2001).

Theorem 1. Consider model [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] where f is an unknown function. Assume that the ε i 's satisfy Assumption (H ε ) and that the family of piecewise polynomials described in section 2.1 has dimensions D m such that [START_REF] Baraud | Model selection for regression on a random design[END_REF] m∈Mn e -LmDm = Σ < +∞ where the L m 's are nonnegative numbers (to be chosen) and Σ is a constant independent of n. Then there exists a universal constant θ > 0 such that if the penalty function is chosen to satisfy

pen n (m) ≥ θs 2 D m (1 + L m )
then the estimator f = f m satisfies

(7) d 2 n (f, f ) ≤ C inf m∈Mn d 2 n (f, S m ) + pen n (m) n + C ′ s 2 Σ n
where C and C ′ are universal constants.

This kind of result can be extended to variables ε i 's admitting only moments of order p, provided that p > 2 (see [START_REF] Baraud | Model selection for regression on a fixed design[END_REF]) for regular collections of models only (i.e. collections with one model by dimension, as in example (RP) defined below). It is worth mentioning that [START_REF] Baraud | Adaptive estimation in an autoregressive and a geometrical β-mixing regression framework[END_REF] allows to compute the rate of the estimator f as soon as f is assumed to belong to some class of regularity authorizing an evaluation of the bias term in function of the dimension of the projection space S m . 2.3. Collections of models and choice of the weights. Let us now illustrate condition [START_REF] Baraud | Model selection for regression on a random design[END_REF] in order to better see the role of the L m 's. Roughly speaking, the final rate for estimating a function of smoothness α is the minimax rate n -2α/(2α+1) when the L m 's can be chosen constant. In most other cases, the L m 's are required to be of order ln(n) and the rate falls to (n/ ln(n)) -2α/(2α+1) . Let us give some (standard) examples for the choice of the spaces when the design is equispaced, namely when x i = i/n: (RP) Regular piecewise polynomials (and regular S m 's). This is typically what is meant when talking about regular collections of models. We work with constant degrees We keep all the degrees constant equal to r -1. We choose the D -1 values of

r 1 = • • • = r D =
a 1 < • • • < a D-1 in the set {j/n, j ∈ {1, . . . , n -1}} for D = 1, . . . , D max = [n/r].
We then have for

L m = L n m∈Mn e -LmDm = [n/r] D=1 n -1 D -1 e -rDLn .
Therefore, if we choose

L m = L n = ln(n)/r it implies that m∈Mn e -LmDm ≤ n-1 k=0 n -1 k e -r(k+1)Ln ≤ n-1 k=0 n -1 k 1 n k+1 = 1 n 1 + 1 n n-1 ≤ 1 + 1 n n ≤ e
and condition ( 6) is satisfied with logarithmic weights. (ITC) Irregular trigonometric polynomials with constant degrees.

The partitions and the a j 's are selected just as before. The degree in this example (but not in practice) is fixed to 2r + 1 in the sense that, on an interval I of weight ℓ we consider Trig ℓ 0

(x) = 1/ℓ1l I (x), Trig ℓ 2p (x) = 2/ℓ cos 2nπ ℓ px 1l I (x), Trig ℓ 2p+1 (x) = 2/ℓ sin 2nπ ℓ px 1l I (x)
, for p = 1, . . . , r. Let us mention that the Trig polynomials would have to be multiplied by √ n to be normalized in L 2 . For the same reason as above, this would lead to weights L m 's of order ln(n), in order for ( 6) to be fulfilled. Note that the first example is meaningless for a non equally spaced design, but the second and third ones can be extended to the general fixed design case by simply choosing the knots

a 1 < • • • < a D-1 in the set {x 1 , x 2 , . . . , x n } with still a 0 = 0 and a D = 1.
Moreover the degrees of the polynomials are supposed to be fixed (to r or 2r + 1) in the previous examples for the sake of simplicity but are variable in the set {0, . . . , r max } in the algorithm developed below. In such case, the dimensional constraint

D max = [n/r] becomes D d=1 (r d + 1)
≤ n where r d ∈ {0, . . . , r max } is the degree of the polynomial on I d . This implies a greater number of models.

2.4.

The aim of the calibration study. The order of the penalty as given in the theoretical results above is only a crude approximation which technically works. One of the aims of the empirical work is to find a more precise development for the choice of the penalty and to calibrate empirically the universal constants involved. For instance, if we think of a penalty: [START_REF] Baraud | Model Selection for (auto-)regression with dependent data[END_REF] pen

n (m) = s 2 c 1 ln n -1 D -1 + c 2 (ln(D)) c 3 + c 4 D d=1 (r d + 1) + c 5 D d=1 [ln(r d + 1)] c 6
for m defined by (2), we need to check that it satisfies [START_REF] Baraud | Model selection for regression on a random design[END_REF]. Then since we believe the constants c i , i = 1, . . . , 6 to be universal constants, we want to compute them using simulation experiments. Note that complementary terms in a penalty function have been studied in a theoretical framework (but for another problem and with a penalty having a different form) by [START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike. (Histograms selection with an Akaike type criterion)[END_REF]. On the other hand, empirical experiments for calibrating a penalty have already been performed for density estimation with regular histograms by [START_REF] Birgé | How many bins should be put in a regular histogram[END_REF]. For all degrees set to zero and regular partitions, they proposed pen n (D) = [ln(D)] 2.5 + D -1. Here, we take c 1 = c 4 = 2 and s 2 = σ 2 . We look for c 2 , c 3 , c 5 and c 6 .

The choice s 2 = σ 2 appears in many theoretical results. If this variance is known, we keep it as the multiplicative factor. Otherwise it can be estimated by the least-squares residuals:

σ2 = 1 n n i=1 (Y i -fm (x i )) 2
for a fm computed on a well chosen S m . For instance, for equispaced design regression, we can take the space generated by a d = d/D with D = n/ ln(n). This has been proved to allow an extension of the theoretical results in the case of regular subdivisions in Comte and Rozenholc (2002).

γ( fm ) = 1 n D d=1 min P ∈P B r d x k ∈I d (Y k -P (x k )) 2 .
In other words, for some given m, we replace the global minimization of the contrast γ in S m by D minimizations of local contrasts denoted

(9) γ I d (g) = 1 n {k/x k ∈I d } (Y k -g(x k )) 2
for g varying in P B r d . Then we have to compute for any degree r, any basis B and any interval I, the polynomial P I,B r ∈ P B r such that:

γB I (r) def = γ I (P I,B r ) = min P ∈P B r {k/x k ∈I} (Y k -P (x k )) 2 /n = 1 n   {k/x k ∈I} Y 2 k - {k/x k ∈I} (P I,B r (x k )) 2   .
This defines the contribution of the interval I to the global contrast. This local contrast is defined only by the points x k and Y k for the indexes k such that x k ∈ I; therefore, any interval I containing the same x k leads to the same minimization procedure and to the same polynomial P I . So there is no loss of generality to consider intervals with bounds chosen among the x k 's.

It is well known that, for any basis B r = (B 0 , . . . , B r ) of a linear space P B r , the contrast minimizer

P I r = α 0 B 0 + α 1 B 1 + . . . + α r B r
is the solution of the system of equations C I r A I r = D I r where (denoting by X ′ the transpose of the vector X), A I r = (α 0 , . . . , α r ) ′ , (10)

C I r = (c s,t ) 1≤s,t≤r , c s,t = k/x k ∈I B s (x k )B t (x k ), (11)
D I r = (d 0 , d 1 , . . . , d r ) , d s = k/x k ∈I Y k B s (x k ). (12)
Let us denote by X I r the matrix (B s (x k )), s = 0, . . . , r, k ∈ {j/x j ∈ I} with r + 1 rows and with #{k/x k ∈ I} = ℓ(I) columns, and by Y I the vector of the Y k 's for x k falling in I. The minimum of the contrast satisfies

n γB I (r) = (Y I ) ′ Y I -(A I r ) ′ X I r (X I r ) ′ A I r = (Y I ) ′ Y I -(D I r ) ′ (C I r ) -1 C I r (C I r ) -1 D I r .
and therefore [START_REF] Birgé | How many bins should be put in a regular histogram[END_REF] γB

I (r) = 1 n [(Y I ) ′ Y I -(D I r ) ′ (C I r ) -1 D I r ].
3.2. Choice of a relevant basis. Since ( 13) is valid for any basis B r of P B r , we look for a relevant choice of the basis of polynomials B r = (B 0 , . . . , B r ) on the interval I. In other words, we aim at choosing the basis such that C I r = I r (I r is the r × r identity matrix), that is

c s,t = k/x k ∈I B s (x k )B t (x k ) = δ s,t
where δ s,t is the Kronecker symbol such that δ s,t = 1 is s = t and δ s,t = 0 otherwise.

In the case of a general design (x i ) 1≤i≤n (not necessarily equispaced), for each interval I, we can build by Gram-Schmidt orthonormalization and using a Q-R decomposition of X, an orthonormal basis of polynomials of any degree r, with respect to the discrete scalar product associated to the x k 's in I. The problem here is that for each possible interval I, and degree r max a specific orthonormalized basis must be computed, which is feasible but costing a lot of time from a computational point of view. Consequently, some other ideas for accelerating the method have to be found. The ideas below are relevant in the equispaced case only.

3.3.

Choice of a relevant basis in the case x i = i/n.

Polynomial basis.

We use the discrete Chebyshev polynomials defined as follows (see [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]). The discrete Chebyshev polynomial on {0, 1, . . . , ℓ-1} with degree r is

(14) Cheb ℓ r (x) = 1 ℓ-1 i=0 [C ℓ r (i)] 2 1/2 C ℓ r (x)
where C ℓ 0 (x) = 1 and

C ℓ r (x) = 1 (r!) 2 ∆ r r s=0 g ℓ (x -s) , where g ℓ (x) = x(x -ℓ) and ∆f (x) = f (x + 1) -f (x), ∆ r = ∆ r-1 o∆. Those polynomials satisfy ℓ-1 k=0 Cheb ℓ r (k)Cheb ℓ s (k) = δ s,t , for 0 ≤ s, t ≤ r.
Therefore, choosing on the intervals I = [i/n, . . . , (i + ℓ + 1)/n[, the basis B I s (x) = Cheb ℓ s (nx -i), will do the job. This leads to [START_REF] Buckheit | [END_REF] γCheb

I (r) = 1 n (Y I ) ′ Y I -(D I r ) ′ D I r ,
where D I r is the vector with components ( 16)

d s = k/n∈I Y k B I s (k/n) = ℓ(I)-1 k=0 Y k+i Cheb ℓ s (k).
3.3.2. Trigonometric basis. The case of piecewise trigonometric bases is even simpler since the basis described in (ITC) is naturally orthonormal with respect to the discrete scalar product considered with a regular design, for ℓ = ℓ(I):

x k ∈I Trig ℓ s (x k )Trig ℓ t (x k ) = δ s,t .
Therefore, choosing on the intervals I = [i/n, . . . , (i + ℓ + 1)/n[, the basis

B I s (x) = Trig ℓ s (x)
, will do the job. This leads to [START_REF] Cai | Incorporating information on neighboring coefficients into wavelet estimation[END_REF] γTrig

I (r) = 1 n (Y I ) ′ Y I -(D I r ) ′ D I r ,
where D I r is the vector with components ( 18)

d s = k/n∈I Y k B I s (k/n) = x k ∈I Y k Trig ℓ s (x k ).
3.4. The choice of the penalty. Equation ( 8) is our choice for the global form of the penalty. For the results given in Theorem 1 to hold, we must prove that m∈Mn e -LmDm < +∞ with pen n (m) = s 2 (1 + L m )D m and m defined by [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF].

2 D exp -c 1 ln n -1 D -1 + c 2 [ln(D)] c 3 +c 4 D d=1 (r d + 1) + c 5 D d=1 [ln(r d + 1)] c 6 + D = Dmax D=1 2 D n -1 D -1 exp -c 1 ln n -1 D -1 + c 2 [ln(D)] c 3 + D × rmax r 1 =0 • • • rmax r D =0 exp -c 4 D d=1 (r d + 1) -c 5 D d=1 [ln(r d + 1)] c 6 = Dmax D=1 2 D exp ln n -1 D -1 -c 1 ln n -1 D -1 + c 2 [ln(D)] c 3 + D × rmax r=0 e -c 4 (r+1)-c 5 [ln(r d +1)] c 6 D . Therefore, if c 1 ≥ 1, c 2 ≥
pen n (m) = s 2 c 1 ln n -1 D -1 + c 2 (ln(D)) c 3 + c 4 D d=1 (r d + 1) + c 5 D d=1 [ln(r d + 1)] c 6
is such that m∈Mn e -LmDm = m∈Mn e -pen n (m)/s 2 +Dm converges with exponential rate, provided that c 1 ≥ 1, c 2 ≥ 0, c 4 ≥ 1.87 and c 5 ≥ 0.

It should be noted that the total number of visited bases is asymptotically (for great values of n and fixed r max ) of the order

n D=1 n -1 D -1 (2r max ) D = (2r max ) (2r max + 1) n-1 = O((2r max ) n ).

Description of the algorithm

In the sequel, both for the description of the algorithm and for the empirical results, we consider only the regular design defined by x i = i/n. Then we find a localized decomposition of the penalized contrast:

nγ( fm ) + pen n,c (m) = pen n,c 1 ,c 2 (D) + D d=1 nγ B r d I d (r d ) + pen n,c 4 ,c 5 ,c 6 (r d ) ,
where the first part of the penalty pen n,c 1 ,c 2 (D) is the global penalization concerning the number of sub-intervals and the second part pen n,c 4 ,c 5 ,c 6 (r d ) is the local part concerning the degree on each sub-interval. We recall that γB I (r) is defined by ( 13) for a basis B r , and more precisely given by ( 15)-( 16) or ( 15)- [START_REF] Cai | Incorporating information on neighboring coefficients into wavelet estimation[END_REF]. This can also be written:

nγ( fm ) + pen n,c (m) = n k=1 Y 2 k + pen n,c 1 ,c 2 (D) + D d=1 pen n,c 4 ,c 5 ,c 6 (r d ) -p B r d I d (r d )
where on the interval I = [i/n, . . . , (j + 1)/n[

p B s I (s) = s t=0 j-i k=0 Y k+i B I t (k) 2 def = p B s (i, j).
The quantity p B s (i, j) represents precisely the weight of the contrast when going from i to j (j included), so that p B s (i, i) is defined. Note that, for 1 ≤ ℓ ≤ n, those quantities are systematically computed by setting Y ℓ = (Y i+k-1 ) 1≤i≤ℓ,1≤k≤n-ℓ and B ℓ = (B I i (k)) 0≤i≤rmax,0≤k≤ℓ-1 and by computing and storing (B ℓ Y ℓ ) •2 where A •2 = (a 2 i,k ) 1≤i≤p,1≤k≤q for A = (a i,k ) 1≤i≤p,1≤k≤q . Then considering different values of ℓ amounts to take into account intervals of any weight ℓ = 1, . . . , n. For j ≥ i, the procedure of minimization first computes:

p(i, j) def = min 1≤s≤rmax min B s σ 2 (c 4 (1 + s) + c 5 ln(1 + s) c 6 ) -p B s (i, j) ,
so that the best basis and the best degree are chosen.

Dynamical programming.

Here we reach the point where we need to use dynamical programming (see [START_REF] Kanazawa | An optimal variable cell histogram based on the sample spacings[END_REF]). The fundamental idea of dynamical programming is that to go to point j with d steps (i.e. pieces), we must first go to some k < j with d -1 steps and then go from k to j in one step. A very similar idea is developed from a theoretical point of view by [START_REF] Kohler | Nonparametric estimation of piecewise smooth regression functions[END_REF] in order to define a procedure in the same type of regression function estimation problem. Let q(d, k) be the minimum of the contrast -penalized in degree with basis selectionto go from 1 to k with d pieces; this value is thus associated to a best partition, d best bases and a choice of d best degrees which fulfill the localization constraints. First note that q(1, k) = p(1, k) which gives an initialization; then

(19) q(d + 1, j) = min d≤k<j [q(d, k) + p(k + 1, j)]
which represents 2j operations. Then a Q matrix can be filled in, with two possible strategies:

(1) "Off line" method: Compute the q(1, j) for j = 1, . . . , n and then do a recursion on d using [START_REF] Coifman | Translation-invariant de-noising[END_REF]. The drawback of the method is that the actualization (i.e. if some more observations are available and so n increases), everything must be done again whereas it is clear that only changes in the last column are useful. (2) "In line" method: Assume that you have built (q(d, j)) 1≤d≤j≤n and you want to increase n and compute the q(d, n + 1), d = 1, . . . , n + 1. Then as q(1, n + 1) = p(1, n + 1) and

q(d + 1, n + 1) = min d≤k<n+1 [q(d, k) + p(k + 1, n + 1)] ,
you only need to compute the p(k + 1, n + 1), 1 ≤ k ≤ n, the q(d, k) being already known. The first part of the work, namely the computation of the coefficients p(i, ℓ) requires O(n 3 r max ) elementary operations, and the dynamical programming part requires O(n 2 D max ) operations. The global complexity of the algorithm is therefore of the order

n 3 r max + n 2 D max .
The implemented method is the former (off line), but for an actualization purpose, the latter method is preferable. Now, on the last column of Q, there are the q(d, n)'s, 1 ≤ d ≤ n, which are the minima of the contrast penalized in degree, to go from 1 to n with d pieces. Thus the last thing to do is to choose

D = arg min d=1,...,n q(d, n) + c 1 ln n -1 d -1 + c 2 ln(d) c 3 .
Of course, the involved partitions must be stored, and not only their number of pieces. As a summary, let us give the steps of the algorithm:

Proposition 2. A model is selected by the algorithm following the steps:

1. On any interval I = [i/n, (j + 1)/n[, compute p B s (i, j) = s t=0 j-i k=0 Y k+i B I t (k) 2 for 1 ≤ i ≤ j ≤ n, 0 ≤ s ≤ r max ,
(i, j) = min 1≤s≤rmax σ 2 (c 4 s + c 5 ln(s) c 6 ) -p B s (i, j) for 1 ≤ i ≤ j ≤ n, 3. Compute p(i, j) = min B∈{Cheb,Trig} p B (i, j), 4. Initialize q(1, k) = p(1, k) for 1 ≤ k ≤ n, and compute recursively for 1 ≤ d ≤ n -1, q(d + 1, n) = min d≤k<n [q(d, k) + p(k + 1, n)] ,
5. Then choose D = arg min d=1,...,n

q(d, n) + c 1 ln n -1 d -1 + c 2 ln(d) c 3 .
The positions of the knots of the involved partitions as well as the selected degrees in step 2 and bases in step 3 must be stored.

4.3.

A fast version of the algorithm. We have also implemented a quick but approximated version of the algorithm, with complexity of order D opt r max n 2 . We must admit that other algorithms, like those based on wavelets for example, have a lower complexity, of the order O(n log 2 (n)). This is a drawback of our procedure. Let us now describe our fast procedure. We construct iteratively a sequence of partitions defined by 0 = a t 0 < a t 1 < • • • < a t Dt-1 < a t Dt = 1 using the following scheme.

• We start at time t = 0 with D 0 = 1, i.e. with the partition with a unique interval.

• At time t of the procedure, -we first check, when D t > 3, if the penalized contrast decreases when removing one a t j for j = 1, . . . , D t -1 from the partition at time t. If some such points exist, we consider the one implying the most significant decrease.

-Next, we check if the penalized contrast decreases when adding one point a in between the a t j 's for j = 0, . . . , D t of the partition at time t. If some such points exist, we consider the one implying the most important decrease.

-If only one such points exist, we follow the associated strategy.

-If two such points exist, we choose the strategy which leads to the most important decrease.

-If no such point exists, the procedure stops. During the procedure, on each involved interval, the best basis (Chebyshev or trigonometric polynomials), the best degrees and the best coefficients are chosen. The interest of such a procedure is that it requires only comparisons between ). Without the potential "removing step", this procedure is like the expansion phase of a CART procedure, see Breiman et al. (1984). In that way, this procedure is related to the CART-type procedure proposed by [START_REF] Kohler | Nonparametric estimation of piecewise smooth regression functions[END_REF], except that he does not consider either the possibility of removing some points from the partition.

Empirical results

Risks and calibration.

First we take for the penalty as defined in Proposition 1,

s 2 = σ 2 , c 1 = c 4 = 2.
Concerning the determination of c 1 , . . . , c 6 , we determined them as follows, using the test functions (and some others, to avoid over-fitting) described below. With the notations of formula (8), we took c 4 = 2 to be in accordance with standard Akaike criterion from an asymptotic point of view. For c 5 and c 6 , we followed some ideas developed by [START_REF] Birgé | How many bins should be put in a regular histogram[END_REF] for the choice of the number of bins of an histogram in density estimation. This question is similar to the choice of the degree in polynomial regression. For c 1 and c 2 , the problem amounts to the choice of an irregular partition when the degree is fixed, for example to zero. We used a rough optimization on the square (c 1 , c 2 ) ∈ [0, 5] × [0, 5] to fix our choice. In any case, we claim that our global choice is satisfactory even if probably not optimal. This choice is the following:

c 1 = c 2 = c 4 = c 5 = 2, c 3 = c 6 = 2.5.
Secondly, we use a set of 16 test functions with very different shapes and regularity. The test functions are given in Figure 1. Functions 1 to 4 and 7 to 14 are the same as the ones used by [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF], functions 1 to 14 come from the Wavelab toolbox developed by Donoho (see [START_REF] Buckheit | [END_REF]) and functions 15 and 16 have been added in order to test also the estimators for some regular functions. Third, we consider different levels (namely 3, 5, 7 , 10) of noise which are evaluated in terms of a signal to noise ratio, denoted by s2n, and computed as

s2n = 1 n n i=1 (f (i/n) -f ) 2 σ 2 , f = 1 n n i=1 f (i/n).
Lastly, the performances are usually compared to a reference value called an oracle. This oracle is the lowest value of the risk. It is computed by using the fact that, in the simulation study, the true function is known, and that, if all the estimates are computed, the one with the smallest risk can be found, as well as its associated risk. In other words in this case, we would have to compute all f -f 2 n , for some known f and all possible f , for all the sample paths. This can be done for regular models but would imply much too heavy computations for the general models considered here. Therefore, another reference must be found to evaluate our estimator. Before giving the details about the wavelet methods, let us explain the reference we use. The index w j denotes the wavelet method number j where 48 wavelet methods are considered and fw j denotes the estimate of f obtained by using the method w j explained below. We generate K (K = 100) samples with length n (n = 128, 512) in the regression model, and denote by f (k) an estimate of f (computed with any method, f (k) w j with the method w j ) based on the k th sample. Then

ℓ 2 2 (f, f (k) ) = 1 n n i=1 (f -f (k) ) 2 (i/n),
and

E * [ℓ 2 2 (f, f )] = 1 K K k=1 ℓ 2 2 (f, f (k) ).
We use the following ratios

(20) R 2 (f ) = E * min j=1,...,48 ℓ 2 2 (f, fw j ) E * ℓ 2 2 (f, f )
.

The ratios are compared to one: the higher over 1 the ratio, the better our method. If our test functions f 1 , . . . f 16 lead to values of R 2 (f i ) for i = 1, . . . , 16 such that ∀i ∈ {1, . . . , 16}, R 2 (f i ) ≥ a, then this means that, for any f ∈ {f 1 , . . . , f 16 },

E * [ℓ 2 2 (f, f )] ≤ 1 a E * min j=1,...,48 ℓ 2 2 (f, fw j ) .
We must emphasize that we chose diadic values of n (n = 128 = 2 7 or n = 512 = 2 9 ) in order to be able to apply all wavelet methods, but our method does not require diadic samples and can be used for any n without any change. We present in Figure 2 an example of data set and estimated signal as performed by our algorithm. The signal has been built with three pieces using functions 15, 8 and 6. The fourth picture gives the variation of the residuals and shows that the algorithm has found an estimator with four pieces, the first one is a standard polynomial with degree 6 corresponding to the estimation of the first function, the second one is a trigonometric polynomial with degree 26 corresponding to the estimation of the second function, the last two pieces correspond to the estimation of the third function, and are polynomials of degree 59 and 54.

5.2.

Comparison with wavelet methods.

Comparison with standard wavelet methods.

Let us be more precise about the wavelets. We use both the MathWorks toolbox developed by [START_REF] Misiti | The Wavelet Toolbox[END_REF] and the WaveLab toolbox developed by [START_REF] Buckheit | [END_REF], following the theoretical works by [START_REF] Donoho | Denoising by soft-thresholding[END_REF], [START_REF] Donoho | Ideal space adaptation by wavelet shrinkage[END_REF], [START_REF] Donoho | Wavelet shrinkage: asymptopia (with discussion)[END_REF]. The abbreviations below refer to the MathWorks toolbox. We use the 6 following basic wavelets: the Haar wavelet (well suited for step signals), two Daubechies DB4 and DB15 wavelets (well suited for smooth signals), two symmetric wavelets abbreviated as Symmlets, sym2 and sym8, the bi-orthogonal wavelet bior3.1 (well suited for signals with rupture). The wavelets are associated with 4 types of threshold: the threshold 2 log(n) called "sqtwolog", the minimax threshold called "minimaxi", the SURE (Stein's Unbiased Risk Estimate) threshold called "Rigsure", an heuristical version of SURE threshold using a correcting term for small values of n, called "Heursure". Lastly, we use the two standard types of threshold, hard and soft thresholding1 . This explains the 6 * 4 * 2 = 48 indexes for the wavelets methods. This is what we call in the following "standard wavelet methods", as opposed to some more recent methods described further. We report in Figure 3 the performances of our estimate when the maximal degree is set to r max = 74 and for s2n= 3, 5, 7, 10, the functions f i being as given in Figure 1, K = 100 and n = 128. More precisely, Figure 3 plots the values obtained for R 2 (f ) relative to the test functions given in Figure 1. We emphasize that the ratio we compute is very unfavorable to our method, because for each sample, we compare our risk to the one of the best (unknown in practice) wavelet method. The ordinate of the lower point is anyway greater than 0.65, which is relatively good since we recall that it means that, for any

f ∈ {f 1 , . . . , f 16 }, E * [ℓ 2 2 (f, f )] ≤ 1.54E * min j=1,...,48 ℓ 2 2 (f, fw j ) .
We have also computed the risks in term of ℓ 1 -type error (squares are replaced by absolute values) where we obtain the same type of results except for the functions f 8 and f 16 where the results are even better: this can be explained by the fact that the ℓ 1 -distance reduces the weights of the discontinuities which are inherent in our method. In addition, the errors are centered Gaussian with known variance, but we also considered centered uniform and Cauchy errors and the results were similar.

5.2.2.

Comparison with recent wavelet methods. We also implemented for comparison some more recent wavelet methods, already studied in [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF] and therefore quite reproducible for such a test. More precisely we considered the following methods, implemented using the Gaussian Wavelet Denoising Library built by [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF] (see http://www.jstatsoft.org/v06/i06/), using either a Haar or a Symmlet8 filter:

W1 Coifman and Donoho (1995)'s translation invariant method using soft thresholding (TI-soft), coded with the function "recTI" in the library, W2 Coifman and Donoho (1995)'s translation invariant method using hard thresholding (TI-hard), coded with the function "recTI" in the library, W3 Cai (1999)'s method using a block non-overlapping thresholding estimator, reusing the first few empirical coefficients to fill the last block, coded with the function "recblockJS" in the library W4 Cai (1999)'s previous method, the last few remaining empirical coefficients being unused, coded with the function "recblockJS" in the library, W5 Huang and Lu (2000)'s method based on nonparametric mixed-effect models, coded with the function "recmixed" in the library, W6 Cai and Silverman (2001)'s method using an overlapping block thresholding estimator, coded with the function "recneighblock" in the library, W7 Antoniadis and Fan (2001)'s hybrid method using a "keep", "shrink" or "kill" rule (SCAD), W8 [START_REF] Vidakovic | BAMS method: theory and simulations[END_REF]'s bayesian adaptive multiresolution method coded with the function "recbams" in the library. For a more precise description of those methods, we refer to [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF]. There are a selection of recent methods that [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF] describe and test, namely methods number [START_REF] Baraud | Model selection for regression on a fixed design[END_REF][START_REF] Baraud | Model selection for regression on a random design[END_REF][START_REF] Birgé | Gaussian model selection[END_REF][START_REF] Birgé | How many bins should be put in a regular histogram[END_REF][START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF][START_REF] Birgé | Minimum contrast estimators on sieves: exponential bounds and rates of convergence[END_REF][START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike. (Histograms selection with an Akaike type criterion)[END_REF] and 34 respectively in their Table 3. We work first with σ = 1 assumed to be known. Moreover, in all the following of the subsection, we use the quick version of our algorithm.

Table 1 below gives the L 2 -errors for the 16 test functions and signal to noise ratio s2n=5 obtained with the quick version of our method (when using Chebyshev polynomials (CP) or both Chebyshev and trigonometric polynomials (CP/TP)) and with the other methods W1 to W8. We must say that we did not succeed in making W8 work, but this may be an error of ours. Besides we found out that the method of [START_REF] Coifman | Translation-invariant de-noising[END_REF] with hard thresholding (W2 or TI-hard) seems to be almost always better than all the other wavelet methods. Our method behaves very well and is in general better than all the other methods. Even when we do not have the lowest errors, we are close to it. Globally, the CP/PT method seems to be preferable: the losses are never very important but the gains are sometimes decisive, when compared to the wavelet methods in competition. 1. L 2 -errors for s2n= 5, n = 512, CP is our method when considering piecewise Chebyshev polynomials only, CP/TP is our method when considering both Chebyshev and trigonometric (piecewise) polynomials, W1 to W8 are the wavelet methods described above with Symmlet8 filter. σ = 1 is known. • gives the best wavelet method, * gives the best method between CP/TP and W1-W8.

Since we found the method of [START_REF] Coifman | Translation-invariant de-noising[END_REF] with hard thresholding (TI-hard) to be the better one, we present a more precise comparison of our results with theirs in Figure 4, in order to illustrate the influence of the choice of the filter (either the Symmlet8 filter or the Haar filter) in the wavelet methods. Our method does not require such a choice, which seems to be sometimes decisive (Signals 7 and 8).

Lastly, we compared our method with others when σ is unknown. We implemented our method using a preliminary estimator of σ 2 based on the mean square residuals obtained with an estimator of the regression function computed on a regular model with D = [n/ ln(n)] intervals in the subdivision and degree r = 3. This estimator is used for the penalization procedure. The estimate of f is then used to re-evaluate σ and to initialize a second penalization. For the test functions 5 and 6, it appears clearly that almost no method makes the job in this case, neither wavelets, nor ours. The only good wavelet method is Huang and Lu (2000)'s method W5 which is never better than the other wavelet methods for the other signals. Note that the test functions 5 and 6 are not used by [START_REF] Antoniadis | Wavelet Estimators in nonparametric regression: a comparative simulation study[END_REF] in their experiments. In the other cases again and as shown by the results given in Table 2, one of the better wavelet methods remains Coifman and Donoho (1995)'s method, contrary to Antoniadis et al. ( 2002)'s conclusion that the best method highly depends on the type of the signal function. Note that we gave for this method the results using both the Symmlet8 and the Haar filter. We give in Table 3 above the estimation performances in term of R 2 (f ) (with respect to the 48 standard wavelet methods) and in term of CPU time (for the same samples) of the accelerated algorithm compared to the standard one.

It appears that except for the first signal, which is better identified by the complete algorithm, the quick algorithm performs very well both in term of risk (which was expected) and time (which was the aim). More precisely and if we except Signal 1, there is essentially no loss in term of risk when using the quick algorithm, but it is between five and twenty times faster for a sample with size n = 512. This effect naturally increases with the sample size. As a conclusion, it is clear that both the standard and the accelerated algorithm work very well. 5.4. Compression performances. We already computed the complexity of our algorithm so that it is clear that even in its quick version, it remains slower than wavelet methods. But it has two decisive advantages with respect to those methods, in addition to its completely automatic feature: first, it performs very well whatever the type of signal, and particularly when discontinuities arise, and secondly, its compression properties are reasonably good, and in particular much better in many cases than wavelets, which was not expected.

Therefore, we also made a simple comparison of the standard wavelet methods and of our algorithm in terms of their compression performances. For each estimated function,

Error ratio

Error ratio

Error ratio

Error ratio

Length ratio

Length ratio Length ratio Length ratio Error ratio in function of the compression ratio. The ratios are our method over all wavelet methods (dotted lines correspond to levels 1). K = 100 samples with length n = 128 and s2n=5.

we compute three types of code lengths: an "integer length" which is a number of integers, namely the number of nonzero wavelet coefficients for the wavelet methods, N intw, and twice the chosen number D of intervals in the Chebyshev piecewise polynomials method N intCP (1 integer for D, D integers for each degree r d , D -1 integers for the length of the intervals), a "real length" which is the number of real coefficients of the developments for each method, N realw and N realCP , and a global length, N w and N CP , defined in both cases as (N int/4) + N real to take into account the fact that an integer is four times smaller than a real number in terms of code length. Figure 5 above plots the error ratios R (j) (f ) = E * ℓ 2 2 (f, fw j ) /E * ℓ 2 2 (f, f ) in function of the global length ratios N w j /N CP where j is the index of the wavelet method, for f taken as each signal of Figure 1 and for s2n= 5. We can see that both the estimation and the compression performances of our algorithm are most of the time better than wavelet methods since both ratios are higher than 1. Signals 5 and 6 are the only one for which the wavelets are better but it appears that they are better either in term of risk or in term of compression, but not both; on the contrary for the other signals, our method is better in terms of both risk and compression performance.

N

  1, a sequence of degrees, that is integers r 1 , . . . , r D , such that for any d, 1 ≤ d ≤ D, 0 ≤ r d ≤ r max , and a sequence C 1 , . . . , C D of binary variables such that if C d = P standard polynomials are considered on [a d-1 , a d [ and if C d = T trigonometric polynomials are considered on [a d-1 , a d [, for d = 1, . . . , D. Then, denoting by (2) m = (D, a 1 , . . . , a D-1 , r 1 , . . . , r D , C 1 , . . . , C D ) we define a linear space S m as the set of functions g defined on [0, 1] that admit the following kind of decomposition: let I d = [a d-1 , a d [ for d = 1, . . . , D -1, and I D = [a D-1 , a D ], then g(x) = D d=1 P d (x)1I I d (x), P d polynomials with degree r d , d = 1, . . . , D, where the polynomial P d is standard if C d = P and trigonometric if C d = T . Note that we may consider standard polynomials exclusively and in such a case C 1 , . . . , C D is not necessary in the definition of m. We define ℓ(I) as the number of x k falling in the subinterval I and we call it "weight of I". The space S m generated in this way has the dimension D m = D i=1 (r i + 1). If we call M n ⊂ {1, ..., D max }× Dmax D=1 [0, 1] D-1 ×{0, . . . , r max } D ×{P, T } D a finite set of all possible choices for m, the family of linear spaces of interest is then {S m , m ∈ M n }.

  r -1 and we choose a j = j/D for j = 0, . . . , D (regular partition of [0, 1]). Then m = (D, a 1 , . . . , a D-1 , r-1, . . . , r-1, P, . . . , P), dim(S m ) = rD, we take D = 1, . . . , D max and we simply impose that rD max ≤ n, i.e. D max = [n/r]. Then we look for L m 's such that m∈Mn e -LmDm = [n/r] D=1 e -LmD ≤ Σ < ∞. Therefore L m = 1 (or L m = 2 ln(D)/D) suits. (IPC) Irregular piecewise polynomials with constant degrees. This illustrates by comparison the extension from regular to general collections of models.

3. Computation of the estimate 3 . 1 .

 31 The general formula. Given a basis B = (B 0 , . . . , B r , . . . ) of polynomials, with degree of B r = r, we denote by B r = (B 0 , . . . , B r ) and by P B r the linear space spanned by B 0 , . . . , B r . The first step for the computation of f = f m is the computation of the fm 's for m varying in M n among which we choose it. Let m = (D, a 1 , . . . , a D-1 , r 1 , . . . , r D , C 1 , . . . , C D ) be given and recall that I d = [a d-1 , a d [ for d = 1, . . . , D -1, and I D = [a D-1 , a D ], a 0 = 0 and a D = 1. Then fm satisfies

  m∈Mn e -LmDm = m∈Mn e -pen n (m)/s 2 +Dm = 1≤D≤Dmax,0≤r d ≤rmax

0 and c 5 ≥ 1 -e -c 4 D,Proposition 1 .

 5141 0, we can give the following bound m∈Mn e -LmDm ≤ Dmax D=1 (2e) D e -c 4 D 1 -e -c 4 (rmax+1)1 -e -c 4 and this last term is bounded provided that2e 1-c 4 1 -e -c 4 < 1 that is, if c 4 > ln(1 + 2e) ≃ 1.862.Thus in the general case, the chosen penalty is of the form: The following choice of the penalty:

4. 1 .

 1 Localization. Here we should emphasize the two basic ideas of our procedure. The first one is based on a localization of the problem. With the results and notations of Section 3.1 and the subsections following, the global value of the contrast is γ( fm ) = D d=1 γB r d I d (r d ) and we look for m = arg min m∈Mn γm + pen n (m) n where we must consider here that pen n (m) = pen n,c (m) with c = (c 1 , c 2 , c 3 , c 4 , c 5 , c 6 ) and pen n,c (m) = σ 2 c 1 ln n -1 D -1 + c 2 ln(D) c 3 + D d=1 σ 2 [c 4 (1 + r d ) + c 5 ln(1 + r d ) c 6 ] def = pen n,c 1 ,c 2 (D) + D d=1 pen n,c 4 ,c 5 ,c 6 (r d ).

  γB [a,b[ (r) + pen n,c 1 ,c 2 (D t ) + pen n,c 4 ,c 5 ,c 6 (r) and γB ′ [a,c[ (r ′ ) + γB ′′ [c,b[ (r ′′ ) + pen n,c 1 ,c 2 (D t + 1) + pen n,c 4 ,c 5 ,c 6 (r ′ ) + pen n,c 4 ,c 5 ,c 6 (r ′′

Figure 1 .

 1 Figure 1. Test functions.

Figure 2 .

 2 Figure 2. An example of decomposition of a signal by the complete algorithm.

16 EFigure 3 .

 163 Figure 3. Performance ratios R 2 (f ) relative to the 16 test functions for K = 100 and n = 128. The greater than one, the better our algorithm.

Figure 4 .

 4 Figure 4. Comparison of the L 2 -errors for the 16 test functions and the four s2n, of our method, using both Chebyshev and trigonometric polynomials (thick curve) and Coifman and Donoho (1995)'s method for the Haar (dashed curve) and the Symmlet8 filters (thin curve). K = 100, n = 512, σ = 1 known.

Figure 5 .

 5 Figure 5. Error ratio in function of the compression ratio. The ratios are our method over all wavelet methods (dotted lines correspond to levels 1). K = 100 samples with length n = 128 and s2n=5.

Table 2 .

 2 L 2 -errors for s2n= 5, n = 512, and σ = 1 is unknown. CP/TP is our method when considering both Chebyshev and trigonometric piecewise polynomials, W1 to W8 are the wavelet methods described above with Symmlet8 filter, W2H is the method W2 when using the Haar filter. Comparison of the complete and fast algorithms. We have described in Section 4.3 an accelerated version of our complete algorithm, but we had to test if the performances of this method were indeed of the same order as the standard one, and nevertheless appreciably faster.

	Signal CP/TP W1	W2	W2H	W3	W4	W5	W6	W7	W8
	1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16	0.092 * 2.08 0.494 0.111 • 0.914 0.934 0.470 1.070 0.976 0.871 0.367 * 2.990 0.461 0.385 • 0.946 0.960 0.510 1.030 1.200 0.871 0.063 * 0.209 0.111 • 0.120 0.201 0.230 0.136 0.217 0.149 0.856 0.206 0.873 0.187 * 0.461 0.274 0.305 0.240 0.215 0.365 0.857 • 4.690 11.4 2.57 16.1 2.38 2.38 0.694 * 2.4 5.98 0.884 • 23.1 23.7 23.5 23.6 23.6 22.8 0.986 * 23.6 23.1 4.05 • 0.013 * 0.571 0.140 0.026 • 0.404 0.420 0.192 0.437 0.295 0.86 0.060 * 0.610 0.073 • 0.290 0.243 0.260 0.275 0.259 0.402 0.864 0.052 * 0.390 0.081 0.064 • 0.167 0.184 0.144 0.208 0.164 0.859 0.090 0.319 0.083 * 0.104 0.128 0.154 0.120 0.136 0.182 0.859 • 0.088 0.204 0.053 * 0.066 0.110 0.142 0.089 0.083 0.092 0.856 • 0.088 0.240 0.050 * 0.123 0.085 0.117 0.094 0.080 0.107 0.858 • 0.150 0.913 0.136 * 0.181 0.265 0.274 0.233 0.250 0.373 0.867 • 0.053 * 0.179 0.063 • 0.072 0.098 0.097 0.074 0.096 0.084 0.857 0.028 * 0.206 0.073 • 0.084 0.158 0.177 0.138 0.169 0.124 0.856 0.080 * 0.751 0.118 • 0.427 0.239 0.258 0.256 0.234 0.362 0.868

• : best wavelet method, * : best method.

5.3.

Table 3 .

 3 Quick and complete algorithm comparison : Risk ratio and CPU Time ratio, ratio = Quick / complete, n = 512 and K = 100.

Table 4 .

 4 ,05 2,09 2,22 1,61 1,20 1,24 1,26 2,30 1,79 1,84 1,93 2 2,92 2,89 1,52 1,47 4,32 4,20 2,35 2,44 3,13 3,08 1,64 1,59 3 2,85 3,57 3,33 3,24 5,25 6,46 6,31 7,30 3,14 3,92 3,68 3,64 4 2,23 2,12 2,11 1,20 7,61 8,72 9,74 6,86 2,59 2,50 2,50 1,44 5 1,47 1,01 1,17 1,06 16,74 12,58 14,98 13,64 1,80 1,23 1,43 1,29 6 1,13 1,16 1,19 1,22 28,89 33,98 35,01 36,30 1,40 1,44 1,47 1,52 7 2,35 2,26 2,40 2,36 1,32 1,33 1,35 1,34 2,03 1,98 2,08 2,05 8 1,64 2,33 2,55 2,43 9,06 18,00 21,97 22,10 1,96 2,83 3,10 2,95 9 1,35 3,71 3,55 2,80 1,80 5,56 5,69 4,94 1,42 3,97 3,84 3,06 10 3,00 4,97 5,19 5,75 4,80 6,87 6,26 6,27 3,24 5,26 5,37 5,84 11 2,18 4,72 4,85 4,68 3,54 8,07 8,89 8,84 2,36 5,15 5,33 5,17 12 4,43 3,15 4,67 4,50 8,79 7,48 11,41 11,52 4,92 3,57 5,30 5,12 13 3,50 3,67 3,68 1,85 5,10 5,40 5,63 2,93 3,73 3,92 3,95 1,99 14 3,40 2,49 2,46 6,11 5,49 4,53 4,68 12,79 3,68 2,73 2,72 6,82 15 2,09 3,45 3,14 3,21 3,42 6,41 6,04 6,85 2,26 3,80 3,47 3,59 16 1,40 1,61 1,83 1,97 9,95 13,01 15,99 18,50 1,69 1,95 2,23 2,40 Compression performances for the three ratios in function of the s2n and of the signal for K = 100 samples with length n = 128.

		realw j /N realCP		N intw j /N intCP			N w j /N CP
			s2n				s2n				s2n
	Signal	3	5	7	10	3	5	7	10	3	5	7	10
	1	2,57 2

Note that the estimation is improved by using the function "wmaxlev" to select the maximum level of the wavelets instead of the standard level round[log

(n)] and we therefore use this MathWorks function as well.

We give also in Table 4 the means (on the K samples) of the ratios: mean of N realw j divided by mean of N realCP , mean of N intw j divided by mean of N intCP and mean of N w j divided by mean of N CP , for the index j corresponding to the best wavelet in term of approximation performance for each sample path; we also distinguish between the values s2n of the signal to noise ratios. We can see that all ratios are greater than 1, which means that our method is better in term of compression that all standard wavelets. In addition, the compression improvement of our algorithm increases when the s2n increases.