Arnaud Lejeune

Hakim Boudaoud

Michel Potier-Ferry

Isabelle Charpentier

Hamid Zahrou

Automatic solver for non-linear partial differential equations with implicit local laws. Application to unilateral contact

Keywords: Asymptotic Numerical Method (ANM), contact mechanics, shell element, implicit surface

In general, nonlinear continuum mechanics combine global balance equations and local constitutive laws. In this work, frictionless contact between a rigid tool and a thin elastic shell is considered. This class of boundary value problems involves two nonlinear algebraic laws: the first one gives explicitly the stress field as a function of the strain throughout the continuum part, while the second one is a nonlinear equation relating the contact forces and the displacement at the boundary.Given the fact that classical computational approaches sometimes require significant effort in implementation of complex nonlinear problems, a computation technique based on Automatic Differentiation of constitutive laws is presented in this paper. The procedure enables to compute automatically the higher-order derivatives of these constitutive laws and thereafter to define Taylor series that are the basis of the continuation technique called Asymptotic Numerical Method. The algorithm is about the same with an explicit or an implicit constitutive relation. In the modeling of forming processes, many tool shapes can be encountered. The presented computational technique permits an easy implementation of these complex surfaces, for instance in a finite element code: the user is only required to define the tool geometry and the computer is able to get the higher-order derivatives.

Introduction

Many industrial processes require the mathematical description of complex surfaces. This work is a part of a research project dealing with the simulation of a sheet metal conveying, with application to steel fabrication plants, especially continuous annealing, levelling or galvanization. In the conveying process, the rolls are imperfect cylinders with a non smooth surface due to wear or to improve the control of process. As a consequence of this non cylindrical shape, wrinkling phenomena can occur, leading to difficulties in the management of the plant [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF] [START_REF] Fischer | Buckling phenomena related to rolling and levelling of sheet metal[END_REF]. Unilateral contact involves a non smooth relation between the contact forces and the distance between the two bodies. The numerical techniques to solve unilateral contact are often based on Lagrange multipliers [START_REF] Papadopoulos | A lagrange multiplier method for the finite element solution of frictionless contact problems[END_REF], augmented Lagrangian [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solutions of Boundary Value-Problems[END_REF] [START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF] and regularization of the non smooth function like in penalty method [START_REF] Kikuchi | A smoothing technique for reduced integration penalty methods in contact problems[END_REF] [7] [START_REF] Wriggers | Finite element algorithms for contact problems[END_REF]. The tool shape may have significant impact on the deformation and the stability of the sheet [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF]. Hence, it is important to describe the complex surfaces with accuracy. Bézier, B-spline and NURBS are among the most popular in the field of simulation of forming processes, see for instance [START_REF] Oliveira | Algorithms and strategies for treatment of large deformation frictional contact in the numerical simulation of deep drawing process[END_REF].

This paper is devoted to a generic modelling approach of unilateral contact, so that the results of the paper could be applied to any structural and processing cases involving contact. A new procedure involving Automatic Differentiation (AD) and Asymptotic Numerical Method (ANM), to account for complex tool shape, is presented. The Asymptotic Numerical Method is an efficient path following technique for the solution of strongly non-linear partial differential equations. For instance it has been applied to problems of unilateral contact [START_REF] Elhage-Hussein | An asymptotic numerical algorithm for frictionless contact problems[END_REF] [START_REF] Potier-Ferry | Traitement des fortes non-linéarités par la méthode asymptotique-numérique[END_REF] [12] [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF], plasticity [START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF] [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF], viscoplasticity [START_REF] Zahrouni | Asymptotic numerical method for strong nonlinearities[END_REF] and to problem coupling several non-linearities [START_REF] Zahrouni | Asymptotic numerical method for strong nonlinearities[END_REF] [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF].

The appearances of wrinkles in strip conveying [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF] or of flatness defects after rolling process of thin sheets [START_REF] Fischer | Buckling phenomena related to rolling and levelling of sheet metal[END_REF] are instability phenomena induced by residual stresses due to fabrication processes. Because ANM has proved to be a very efficient continuation algorithm in order to predict all sorts of buckling phenomena, see for instance [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF], it is a good candidate to design software dedicated to the simulation of strip conveying. Up to recently, ANM was not well adapted to the description of generic tool shapes and to account for plasticity and friction. The present work focuses on the first point with the help of Automatic Differentiation while the treatment of plasticity and friction is presented in companion papers [START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF] [START_REF] Béchet | Taylor series to solve friction problems[END_REF].

The ANM algorithm does not differ significantly from the one presented in [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF]. The contact law is regularized by a smooth function and the main unknowns are the displacement and the contact forces as in many contact algorithms. Up to now within ANM, the tool surfaces have been represented by simple surfaces as planes or cylinders [START_REF] Elhage-Hussein | An asymptotic numerical algorithm for frictionless contact problems[END_REF] [START_REF] Aggoune | High order predictioncorrection algorithms for unilateral contact problems[END_REF] [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF]. The main goal of this paper is to establish an efficient and generic methodology to account for complex surfaces. Then, AD is used to compute the higher-order derivatives of constitutive laws. This allows for the automation of the implementation of various complex tool shapes, tool motions or gap-pressure relations. There have been many works about Automatic Differentiation [START_REF] Griewank | On the efficient generation of taylor expansions for dae solutions by automatic differentiation[END_REF] [START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF] and about its applications to partial differential equations using or not finite element method ( see for instance [START_REF] Korelc | Automation of primal and sensitivity analysis of transient coupled problems[END_REF], [START_REF] Tijskens | FastDer++, efficient automatic differentiation for non-linear PDE solvers[END_REF], [START_REF] Tijskens | Automatic differentiation for solving nonlinear partial differential equations: An efficient operator overloading approach[END_REF], [START_REF] Bartlett | Automatic differentiation of C++ codes for large-scale scientific computing[END_REF], [START_REF] Bartlett | Hybrid differentiation strategies for simulation and analysis of applications in C++[END_REF],... ), but only few papers have described coupling of AD and ANM [27] [28] [21] [START_REF] Lejeune | Object-oriented design to automate a high order nonlinear solver based on asymptotic numerical method[END_REF].

The layout of the paper is as follows: Section 2 presents the higher-order differentiation of some local constitutive laws for nonlinear PDE problems and Section 3 introduces the treatment of the global problem with unilateral contact. Section 4 illustrates the resolution of a global contact problem with complex contact surfaces. Some perspectives are proposed as a conclusion.

Higher-order differentiation of local constitutive laws

In Asymptotic Numerical Method, the unknowns U of the non-linear problem are written in the form of truncated Taylor series in the vicinity of a starting point U 0 :

a ∈ R → U (a) = U 0 + K k=1 a k U k . (1)
Hence the computation of a solution branch is equivalent to the computation of the coefficients of the expansion that are obviously related to the higher-order derivatives of the function :

U k = 1 k! ∂ k U ∂a k (0) . (2)
Suppose now that two variables of the problem are related by an algebraic relation such as the constitutive law of non-linear elasticity that connects the strain and the stress. It is well known within ANM [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF] [START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF] that such a law leads to a constitutive law at order k: the higher-order derivative of the stress is a function of derivatives of the strain. The goal of this part is to show how this constitutive law at order k could be implemented automatically using Automatic Differentiation, even if the starting law is complex. As usual in continuum mechanics, the full problem associates global balance equations and local constitutive laws. In this part, we focus only on the treatment of constiitutive laws. The full problem will be discussed in Part 3.

Algebraic constitutive laws

Material modeling

Let us recall that the deformation gradient tensor F is related to the displacement gradient tensor θ = ∇ x u by F = I 2 + θ where I 2 is the second order identity tensor. The Green-Lagrange strain tensor γ is defined by

γ = 1 2 T F F -I 2 = 1 2 T θ + θ + T θθ .
(3) noted T and S related through T = (I 2 + θ) S.

In elasticity framework, S and T are related to the Green Lagrange second order strain tensor by the fourth order elastic tensor L:

T = (I 2 + θ) S, (4)
S = L : γ, (5)
γ = 1 2 T θ + θ + T θθ . ( 6
)
The computation of the derivatives of the constitutive laws ( 4)-( 6) is well known and is accounted by straightforward formulae :

T k = (I 2 + θ 0 ) S k + k i=1 θ i S k-i , (7)
S k = L : γ k , (8)
γ k = 1 2 ( T θ k + θ k ) + k i=0 T θ i θ k-i . (9)
AD avoids the implementation of the derivative formulae ( 7)-( 9) : the user has only to implement the basic laws ( 4)-( 6), the derivatives being evaluated by the use of operator overloading technique. This is a direct mode of automatic differentiation.

Contact modeling : the geometrical part

Projection algorithm consists in computing the distance h and the normal vector n between the slave point and the rigid master surface. The normal vector is defined on the rigid master surface. By denoting X 0 the initial coordinate of the observed slave point and u its displacement field, the slave node is located at the coordinate X = X 0 +u. Forthcoming subsections briefly describe different kinds of functions used to express the projection function linking the displacement field u to the projection variables h and n. This function depends on the master surface geometric model. Papers on sheet metal forming are related to Bezier, NURBs or analytic description [START_REF] Oliveira | Algorithms and strategies for treatment of large deformation frictional contact in the numerical simulation of deep drawing process[END_REF].

Explicit Projection This case is the most simple. The projection variables h, n are explicitly related to the displacement field u:

proj (h, n, u) = T {h, n} -f (u) = 0. ( 10
)
with f a given function depending on the surface definition. Higher-order differentiation is of no problem and treated as the elastic constitutive equation.

Implicit Projection with non-parametric surface The master surface is described by geometrical variables (for example radius and center coordinates for a sphere). The relationship between the normal vector, the distance and the displacement field is given by the following implicit function :

proj (h, n, u) = 0. ( 11
)
For example, let us consider a perfect cylinder, given by its center point C, its radius value r with an axis parallel to (Oy). The projection on this cylinder is given by:

proj (u, h, n) =        (h + r)n x -x 0 + C x -u x = 0 n y = 0 (h + r)n z -z 0 + C z -u z = 0 n, n -1 = 0 ( 12
)
Implicit Projection with parametric surface Another way to describe an obstacle is to use parametric surfaces such as Bezier surfaces for instance.

Let us denote by S (s) the parametric surface where we suppose that s is in a bounded domain of dimension L.

To compute the variables h and n, one has to know the coordinate of the point projected from the continuum media (deformable solid) onto the parametric surface. Then finding the distance h and the normal n is equivalent to solving the following problem: Find (h, n, s) such that:

proj(h, n, s, u) =    X 0 + u -S (s) -hn = 0 n, n -1 = 0 n, ∂S ∂s (s) = 0 (13)
So considering parametric surfaces leads to increase of the dimension of the local projection problem by the dimension of parametric variables.

Contact modeling : the mechanical part

In this paper, we only discuss frictionless contact. The contact pressure is assumed to be normal to the two bodies:

R = R N n. (14)
The Asymptotic Numerical Method (ANM) requires smooth constitutive laws. As proposed in [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF] the hard contact law (R N ≥ 0, h ≥ 0, R N h = 0) may be replaced by a smooth relation between the normal pressure R N and the gap h, following:

R N ( R N α + h) = η(∆ -h). (15)
In [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF], ∆ is an initial gap (∆ > 0) and α and η are two parameters. The given number η has the dimension of a pressure and it is assumed to be small. For a very small η, the relation [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF] tends to a piecewise linear function

R N ≈ 0 if h > 0, R N = α(-h) if h < 0.
Hence the parameter α is a kind of contact stiffness parameter and it has to be large. If 1/α = 0 the penetration is prevented and the relation [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF] leads to a barrier function as in [START_REF] Zavarise | Real contact mechanisms and finite element formulation: a coupled thermomechanical approach[END_REF]. The barrier contact law (1/α = 0) has been used in many papers about ANM, for instance [START_REF] Elhage-Hussein | A numerical continuation method based on padé approximants[END_REF] [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF]. By combining [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF] with the projection equations of subsection 2.1.2, the contact normal pressure R N can be considered as an implicit function of the displacement u. In the same way, the pressure vector R is considered as a function of R N and u. In the next part, we present a general numerical procedure to compute the higher-order derivatives of these functions u → R N and (R N , u) → R.

DIAMANT approach applied to local laws

Continuum mechanics deal with non-linear problems implying complex local non-linear laws. Such problems can be solved using the ANM through the computation of higher-order derivatives in a particular direction. The approach, named DIAMANT, coupling Automatic Differentiation (AD) and ANM has been proposed in [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche diamant[END_REF].

DIAMANT approach

As noted previously in the introduction of part 2, unknowns are expanded in Taylor series versus a scalar path parameter. ANM algorithm consists in finding each Taylor coefficient. These terms are related to the higher-order derivative with respect to the path parameter. The non-linear problem is also expanded in Taylor series leading to a system of linear equations at each order. Expressing higher-order equations with respect to Taylor coefficients could be difficult and error-prone. Then coupling ANM with AD (namely the DIAMANT approach) allows to automate Taylor series terms in a specific direction using well-known recurrence formula.

Recently, DIAMANT approach has been proposed and applied to several kinds of academic problems [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche diamant[END_REF], [START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF], [START_REF] Koutsawa | A generic approach for the solution of nonlinear residual equations. part i: The diamant toolbox[END_REF]. Right-hand side terms and tangent matrix occurring in ANM are obtained by application of AD technique. This technique has been originally proposed in [START_REF] Griewank | On the efficient generation of taylor expansions for dae solutions by automatic differentiation[END_REF], in which the package ADOL-C for automatic differentiation [START_REF] Walther | A Package for the Automatic Differentiation of Algorithms Written in C/C++ -Version 2[END_REF] is employed. It is quite related to several works about solving ODE problems (see for example [START_REF] Jorba | A software package for the numerical integration of odes by means of high-order taylor methods[END_REF]) or to Tadiff software [START_REF] Bendtsen | TADIFF, a flexible C++ package for automatic differentiation using taylor series expansion[END_REF].

Thus we retrieve the same differences between our approach and the use of ADOL-C as noted in [START_REF] Jorba | A software package for the numerical integration of odes by means of high-order taylor methods[END_REF]. Moreover, [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche diamant[END_REF] recalls that Taylor series computation does not need higher-order derivatives in all directions. In this sense, it is more simple than classical higher-order AD tools.

To illustrate this approach, let us recall the automatic computation of a variable defined by an explicit function itself expressed as a sequence of elementary operations. Let y be a variable defined as:

y = f (x) . ( 16
)
Using Faà di Bruno's formula as in [START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF], one could demonstrate that the latter relation is linear with respect to x k and therefore can be written as:

y k = ∂f ∂x (x 0 )x k + f nl k (x k-1 , ..., x 1 , x 0 ) , f nl k = f k (0, x k-1 , ..., x 1 , x 0 ) .
In the latter formula, two main quantities appear that are fundamental within ANM; the tangent matrix ∂f ∂x (x 0 ) and the right-hand side f nl k that can be seen as a higher-order derivative of the function f (x). Hence, implementing ANM is equivalent to computing higher-order derivatives of function. Derivatives up to the order k of y could be computed with AD by using operator overloading. A sequence of computation could be:

y 0 = f (x 0 ) , y 1 = f 1 (x 1 , x 0 ) , . . . y K = f K (x K , . . . , x 1 , x 0 ) .
Principle of computation of k th derivative of f function by AD is based on the differentiation chain rule. The analytical function f depends on x from a compound of elementary arithmetic operators. As the higher-order differentiation of these elementary operators is well known, applying the function f onto differentiation of variable would lead to the computation of higher-order derivatives f k of the function f . As an example, consider this arbitrary function:

f (x) = x 2 + e 3x . ( 17
)
Naïve successive computation operations could be described by pseudo-code given in Algorithm 1. AD by operator overloading described in subsection 2.2.3 would lead to the computation of derivatives using recurrence formulae. For instance a differentiation code specific of the function described by the Algorithm 1 is given in Algorithm 2. The generic AD code combines elementary operations that corresponds to each line of Algorithm 2.

Algorithm 1 Successive operations to compute

f (x) = x 2 + e 3x 1: u = x * x 2: v = 3 * x 3: w = e v 4: f = u + w Algorithm 2 Successive operations to compute derivatives of f (x) = x 2 + e 3x 1: u k = k l=0 x l * x k-l 2: v k = 3 * x k 3: w k = 1 k k-1 l=0 (k -l) w l v k-l 4: f k = u k + w k 2.2.

DIAMANT approach extended to local implicit functions

Until now, problems treated with DIAMANT approach only deal with local explicit functions. A difficulty arises when considering local implicit equations for the contact projection especially for complex contact surfaces. Then we would like to extend DIAMANT approach to treat more complex industrial problems modelled with continuum mechanics and complex local laws. This extension is turned into a generic approach to solve analytical non-linear problems with complex local laws. The originality of our work consists in mixing non-linear global equations and computation of local variables from local implicit and explicit non-linear laws. Compared to explicit function, dealing with an implicit function is slightly different. Explicit function is computed with a formula that is composed of its variables and elementary operations. On the contrary, an implicit function needs the use of an algorithm to express one variable with respect to the other one. Assume y is related to x by such a function (with all the needed properties for the sequence):

g (x, y) = 0. ( 18
)
One starts from a known solution (x 0 ,y 0 ): in practice, this starting point will be the endpoint of the previous step. The derivative expression of the order k > 0 of equation ( 18) is given by:

g k = ∂g ∂x x k + ∂g ∂y y k + g nl k (x k-1 , . . . , x 1 , x 0 ; y k-1 , . . . , y 1 , y 0 ) = 0. ( 19
)
As the former derivatives are linear with respect to x k , y k , the relation [START_REF] Béchet | Taylor series to solve friction problems[END_REF] Algorithm 3 Computation of a term defined by an implicit function Require: x 0 , y 0 / g (x 0 , y 0 ) = 0

1: Compute k T = -∂g ∂y -1 2: ORDER 1 3: Compute g1 = g 1 | y1=0 4: Compute y 1 = k T g1 5: ORDER k > 1 6: for k = 2 to K do 7: Compute gk = g k | y k =0 8:
Compute y k = k T gk 9: end for 10: return y 0 , y 1 , . . . , y k can be inverted by classical algebra:

y k = - ∂g ∂y -1 ∂g ∂x x k + g nl k (x k-1 , . . . , x 1 , x 0 ; y k-1 , . . . , y 1 , y 0 ) . (20)
To automate the process, we note that the term on the right hand side corresponds to the computation of g k restricted to y k = 0. Consequently, a driver has been proposed to compute implicit function derivatives up to the order k following Algorithm 3. When the implicit function is a local one included in a more general process, computation time is reduced by storing each result at a given order and just computing the given order. Note that this computation of the higher-order derivative can be seen as a specific higher-order extension of the Post-Differentiation defined in [START_REF] Bartholomew-Biggs | Using forward accumulation for automatic differentiation of implicitly-defined functions[END_REF]. They conclude that, in the first order derivative case, the Post-Differentiation approach seems to be less expensive than a Full Differentiation one in terms of cost of iterations.

AD with Operator Overloading technique

This section is a brief recall of the AD technique in forward mode applied to the computation of higher-order Taylor series terms.

Following the definition of [START_REF] Bischof | Computing derivatives of computer programs[END_REF], AD could be viewed as the way to generate a code that computes the derivative of a function itself given by a code. A natural approach is to use the forward mode by applying operator overloading. This is notably well-suited to C++ language [START_REF] Jerrell | Function minimization and automatic differentiation using c++[END_REF]. This language allows to overload classical operators in a simple way [START_REF] Stroustrup | The C++ Programming Language[END_REF]. Notice that AD by operator overloading has already been implemented since a long time in other classical languages [START_REF] Jerrell | Automatic differentiation using almost any language[END_REF]. Interested readers can refer to [START_REF] Bischof | On the implementation of automatic differentiation tools[END_REF] for a good review of AD techniques. But in most cases, works are done for low order derivatives (first or second order). We are here concerned with higher-order derivatives in direct mode as in [START_REF] Rall | Automatic Differentiation: Techniques and Applications[END_REF] but we only need the computation of the derivatives in a specific direction.

From a programing point of view, higher-order differentiation by operator overloading may be realized using template facilities of C++ language. First, the function to be derived is implemented and tested in the continuous domain. Then the template parameter is changed into the higher-order derivative one. Consequently, one single implementation allows the computation of continuous residuals and derivatives. A particuliar procedure is applied in the case of local implicit functions due to the different nature between the continuous function (implicit) and the derivative function (explicit).

As an illustration a naïve templated implementation of the example function f (x) = x 2 + e 3x presented in subsection 2. ✡ ✝ ✆ One could note that switching between scalar and derivative computation results in changing variable type. As recurrence formulae are described in the T Diff class with operator overloading ability, computing the scalar function would then apply the Algorithm 1 successive operations while computing the derivative function woud refer to Algorithm 2 operations.

Higher-order derivative of the projection equations

Explicit Projection

This is the simplest case. Following AD well-known technique, derivative of the local law at the order k is given by:

T {h k , n k } = f k (u) . ( 21
)
with f k the derivative at the order k of the given function f depending on the surface definition. This derivative could be expressed in term of linear and non-linear part (f nl k ) as:

f k (u) = ∂f ∂u u k + f nl k (u 0 , u 1 , . . . , u k-1 ) . (22)
Using AD technique results in implementing function f and letting AD compute the Taylor series terms (linear and non linear).

Implicit Projection with non-parametric surface

Dealing with an implicit local projection law leads to solving the following problem Knowing u, find h (u) and n (u) such as proj (h, n, u) = 0.

Using subsection 2.2.2 and denoting Π = (h, n), one could apply the Algorithm 3 so that the non-linear local laws are replaced by linear ones:

Π k = -∂proj ∂Π (Π 0 , u 0 ) -1 ∂proj ∂u (Π 0 , u 0 ) u k + proj nl k (u 0 , . . . , u k-1 ; Π 0 , . . . , Π k-1 ) . (24)

Implicit Projection with parametric surface

This is the same as implicit projection with non-parametric surface except that the problem dimension is increased by the number of parameters defining the parametric surface:

Knowing u, find h (u), n (u) and s (u) such as proj (h, n, s, u) = 0.

Denoting Π = (h, n, s) yields a formula similar to [START_REF] Tijskens | Automatic differentiation for solving nonlinear partial differential equations: An efficient operator overloading approach[END_REF] :

Πk = - ∂proj ∂ Π Π0 , u 0 -1 ∂proj ∂u Π0 , u 0 u k + proj nl k u 0 , . . . , u k-1 ; Π0 , . . . , Πk-1 . ( 26
)
3 Higher-order differentiation of global problems with unilateral contact

The library described in part 2 is applied in solving the finite element modeling of unilateral contact problems between a rigid tool and a deformable elastic body. The constitutive assumptions are those presented in section 2.1: the geometrical quantities (strain, gap between the two bodies) are defined exactly. The stress-strain relation is linear, the contact pressure is defined by a regularized version of the hard frictionless contact law. These constitutive assumptions are accounted by three local functions. The elasticity laws leads to an explicit function θ = ∇u → T , the smooth contact law requires two functions: the first one defines the normal pressure as a function of the displacement (u → R N ) and the second one defines the pressure vector as a function of the displacement and of the normal pressure (R N , u) → R. To implement the ANM in a finite element code, one needs only linear relations deduced from the higher-order derivatives of those functions: these higher-order constitutive laws are computed according to the automatic differentiation procedures presented in section 2.2. The constitutive functions of the present paper have been chosen for simplicity. The contact law is more or less generic in the frictionless framework, but the stress-strain relation is very specific. The present approach could be extended for instance to hyperelasticity, which has been solved by ANM in few papers [START_REF] Galliet | Une version parallèle des méthodes asymptotiques numériques: Application à des structures complexes à base d'élastoméres[END_REF] [44] [START_REF] Nezamabadi | Solving hyperelastic material problems by asymptotic numerical method[END_REF], because hyperelasticity leads also to smooth algebraic functions θ → T (θ). Alternative constitutive laws like plasticity [START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF] or friction [START_REF] Béchet | Taylor series to solve friction problems[END_REF] implies Differential Algebraic Equations (DAE) that requires a slightly different treatment, but the DIAMANT approach can also be used in case of DAEs.

Unilateral contact with an elastic 3D body

The deformable body is discretized with 3D elements. A mixed formulation in displacement/contact pressure is proposed. The balance of momentum equation and the boundary conditions are accounted to define a well posed boundary value problem. The equivalent weak form is the virtual work equation. Most of the computer codes involving contact mechanics are based on a mixed approach: the discrete unknowns are the displacement and the contact pressure. Within the ANM, this Lagrange multiplier method seems better than a pure displacement method [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF]. Hence the full elastic-contact problems can be summarized as follows:

Find u, R N , λ such as :

       Ω T .∇δudΩ - ∂Ω R(R N , u).δudΓ -λP e (δu) = 0, R N = RN (u). ( 27
)
where the constitutive functions T = T (∇u), R = R(R N , u) and R N = RN (u) have been defined in part 2. To compute the series u (a), R N (a) and λ (a), one needs the derivatives at any order k of their constitutive functions. They are in the form :

T k = ∂ T ∂θ θ k + T nl k = ∂ T ∂θ ∇u k + T nl k , (28)
R k = n 0 .R N k + ∂ R ∂u .u k + R nl k , (29)
R N k = ∂ RN ∂u .u k + R N nl k . ( 30
)
In the formula (28), ∂ T ∂θ is nothing but the usual tangent stiffness matrix at the starting point of the series. The details of the functions T k (θ k , θ k-1 , . . . , θ 1 , θ 0 ) have been given in formulas ( 7), ( 8) and ( 9). The latter formulae ( 29) and ( 30) combine the order k derivatives of the scalar functions R N (h) given in [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF] and the k th -derivatives of the functions h (u) and n (u) computed like in subsection 2.3.3.

Finally, combination of ( 28), ( 29) and ( 30) and the obvious k th -derivative of the problem ( 27) leads, at a given order k, to the linear set of equations satisfied by the mixed unknown u k , R N k and the k th component λ k of the load parameter:

Find u k , R N k , λ k such as :            Ω ∇δu. ∂ T ∂θ .∇u k dΩ -∂Ω δu. ∂ R ∂u .u k dΓ -∂Ω c δu.n 0 R N k dΓ c -λ k P e (δu) = ∂Ω c δu.R nl k dΓ c -Ω ∇δu.T nl k dΩ -∂ RN ∂u u k + RN k = R N nl k ( 31
)
In this form, the linear problem [START_REF] Elhage-Hussein | A numerical continuation method based on padé approximants[END_REF] does not differ significantly from the formulations presented in [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF]. So the discretization is quite the same. The main difference lies in the computation method of the right hand sides R nl k , T nl k and R N nl k and in the application to any tool shape described with explicit or implicit functions.

Unilateral contact with an elastic shell

In view of application to sheet metal conveying process simulation, the body has to be discretized by a shell element. As in [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF], a three-dimensional shell formulation proposed by Buchter-Ramm-Roehl in [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF] is used. An advantage of this element is its strong similarity with volumetric elements and moreover the three-dimensional constitutive law could be implemented directly. Hence the mathematical formulation and the discretized technique would be quite similar to those of the three dimensional case. The main difference is the introduction of a new strain variable, denoted by γ to account for a linear variation of γ 33 through the thickness. This extra variable is incorporated via the Enhanced Assumed Strain (EAS) concept introduced in [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF]. The elastic constitutive law (4), ( 5), ( 6) is modified. It is an explicit function (θ = ∇u, γ) → (T, S) that is written as:

T = (I 2 + θ) S (32)
S = L : (γ c (θ) + γ) , (33)
γ c (θ) = 1 2 T θ + θ + T θθ . ( 34
)
where [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche diamant[END_REF], γ c (θ) is the conventional Green-Lagrange strain.

As a counterpart, a new equation is introduced that is an orthogonality condition between the second Piola-Kirchhoff stress tensor and the enhanced part of the strain γ.

The full contact problem ( 27) is modified as follows:

Find u, R N , λ such as :

         Ω T (u, γ) : ∇δudΩ - ∂Ω R(R N , u).δudΓ -λP e (δu) = 0 R N = RN (u) Ω Ŝ (u, γ) : δγdΩ = 0 (35)
The higher-order differentiation of ( 35) is straightforward and quite similar to the volumetric case, the only difference being the account of the last equation leading to an affine relationship between ∇u k and γk .

Finite element discretization, ANM algorithm

The key point of this work is the automatic computation of the linear problem defining the derivatives at any order which has been presented up to now. The remainder of the method is very classical: spatial discretization using the Finite Element Method and continuation technique within the ANM.

In the application discussed in the following section, we limit ourselves to the contact between a thin shell and a rigid tool. The deformable body is discretized using 8-nodes quadrilateral shell elements with quadratic shape functions introduced in [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF] and based on the EAS concept [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF].

The implementation in a finite element framework is the same as in [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. The contact pressures are represented by forces at each node as in [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF].

The Finite Element discretization and the assembly step are not impacted by the AD technique that is only applied at the level of the elements. For instance, the four degrees of freedom accounting for the additional strain γ are eliminated at this level using a condensation technique. Note also that the continuation algorithm is really independent of differentiation method. Then it is not affected by the AD technique.

Once the series are computed, the length of the solution branch is defined by an estimate of the convergence radius of the series, as discussed for instance in [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]:

a max = δ u 1 u K 1 K-1 . ( 36
)
where δ is a user-defined accuracy parameter (typically δ = 10 -3 or δ = 10 -6 ). Next another computation of series is performed from the step end defined using a max and so on. Two numerical tests of contact between a thin plate and a rigid roll are considered. The aim is to show that complex tool shapes can be easily and efficiently modelled with the help of the presented numerical tool software that couples ANM and AD. Two roll shapes are considered. In the first test, it is a perfect cylinder (Figure 1) and in the second case a tapered roll that combines a perfect cylinder and a cone (Figure 2). The main geometric data are the dimensions of the plate (length = 400 mm, width = 20 mm, thickness = 0.8 mm) and the radius of the roll (r = 50 mm) used in the two cases. According to Figure 2, additional data defining the tapered roll are the flat length (length of the cylindrical part) l 1 = 15 mm, the length l 2 = 5 mm and the angle between the cone and the cylinder p = 0.1 rad. The implementation of the cylindrical roll is done by the implicit representation given in equation ( 13). In the case of the tapered roll, the geometric equation is given by:

proj (u, h, n) =      (h + r p ) n -X + C -u = 0 r p = r 1 + tan(p) 2 (y + u y -l 1 ) 1 + tanh y+uy-l1 ǫ n 2 -1 = 0 (37)
where r, p, l 1 are the parameters described above, C is the center point of the roll and ǫ = 23 a parameter allowing a smooth approximation of the tanh function. With the equation [START_REF] Bischof | Computing derivatives of computer programs[END_REF], the non smooth surface of the tapered roll is replaced by a smooth one that differs slightly to the previous one near the vertex. Let us underline that the AD method makes easier modelling of the tapered surface, that requires only the implementation of the three formulae [START_REF] Bischof | Computing derivatives of computer programs[END_REF]. The computation of the Taylor series is done automatically as explained previously.

As the aim of this section is to validate the approach, only 20 elements discretize the plate (Figures 1 and2).

Concerning the boundary conditions applied on the deformable body (i.e. the plate), the long edges are free. One of the short edges is clamped and the second one is submitted to an uniformly distributed load λF , where λ is the load parameter mentioned in equation [START_REF] Bendtsen | TADIFF, a flexible C++ package for automatic differentiation using taylor series expansion[END_REF] and the nominal force F is 100 N. The material is elastic, linear and isotropic with data E = 220 MPa and ν = 0.3. The contact law is the regularized law (15) that depends on two numerical parameters: a stiffness parameter α = 10 3 MPa.m and a non dimensional regularisation parameter η = 0.1. They have been chosen in such a way that the regularized contact law is sufficiently close to the hard contact law. The initial distance between the plate and the tools is 4 mm.

Numerical results

The numerical results obtained by the present method are compared with a well established commercial code (Abaqus) using a hard contact law. Let us mention that the efficiency of Asymptotic Numerical Method for contact problem has been assessed in several previous papers, see for instance [START_REF] Aggoune | High order predictioncorrection algorithms for unilateral contact problems[END_REF] [17] and this will not be repeated here. Few numerical results are presented in Figures 3,4, 5 and 6. In the case of a cylindrical roll (Figure 3), the vertical displacement on a point submitted to the applied load is plotted versus the applied load λ. One sees many slope discontinuities, due to appearance or disappearance of new nodes in contact. To get the response curve of Figure 3 with an ANM order equals 20 and the accuracy parameter δ = 10 -6 , 95 ANM steps have been necessary, see [START_REF] Aggoune | High order predictioncorrection algorithms for unilateral contact problems[END_REF] or [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF]. Note that the number of steps could be reduced by using Padé approximants [START_REF] Elhage-Hussein | A numerical continuation method based on padé approximants[END_REF] or by changing the accuracy parameter δ or the regularization parameters, η or α. In Figure 4, the evaluation of three contact forces according to λ is plotted and compared to Abaqus results (reference curves). At the beginning (λ = 0), there is no contact forces, next the contact occurs at each of these three points, then the nodes take off from the cylinders. The comparison with Abaqus results are in good agreement. Figures 5 and6 illustrate the effect of the conical part on the displacement curve between a cylindrical roll and a tapered roll on three points positioned on Figure 2. Comparison between ANM solver and Abaqus (reference curves) again shows good agreement. Looking at the point 3 curve on Figure 6, it is important to notice that the regularisation with a tanh function for the ANM solver allows to obtain a good description of the tapered roll surface, notably the intersection between conical and cylindrical parts. At our knowledge, it is the first time that an unilateral contact problem with this kind of non regular obstacle is achieved using the ANM. The rolls are imperfect cylinders with a convex shape. We consider here the process on one tapered roll, as described in Figure 7, using the implicit geometric representation equation [START_REF] Bischof | Computing derivatives of computer programs[END_REF] and the equilibrium defined by the equation [START_REF] Bendtsen | TADIFF, a flexible C++ package for automatic differentiation using taylor series expansion[END_REF]. This example is quite similar to the previous one excepted that a tensile load is applied, only half of the plate is modelled and geometric data are similar to those typically encountered in industrial conveying process. The main geometric data on Figure 7 rad. Note that for the projection equation ( 37), the parameter l 1 equals lt 2 . The parameter ǫ equals 23 as in the previous case.

The upper long edge of the strip is free from any boundary condition. The lower long edge uses symmetric boundary conditions: the model is constrained in the Y-translation. One of the short edges uses symmetric boundary conditions: the model is constrained in the X-translation. To avoid strip wandering during the conveying process, the second short edge is submitted to an uniformly distributed tensile load λF , where λ is the load parameter and the nominal force F equals 1 MPa. The material is elastic, linear and isotropic with data E = 210 MPa, ν = 0.3.

The contact law is the same used in the previous validation cases: regularized law [START_REF] Zahrouni | Solving plasticity problems by a perturbation technique[END_REF] with a stiffness parameter α = 10 6 MPa.m and an adimensional regularisation parameter η = 0.1.

The deformable body is discretized using 256 8-noded quadrilateral shell elements with quadratic shape functions mentioned in section 3.3.

The non convex profile roll leads, as explained in [START_REF] Gueydan | Modélisation numérique du flambage de bandes en acier dans les recuits continus[END_REF], to buckling which is sensitive to the profile geometric parameters: the table length l t and the angle p. This problem is simulated with 28 ANM steps using the ANM order equals 20 and the accuracy parameters δ = 10 -6 .

Figure 8 shows the deformed strip for a tensile load λF equal to 40 MPa. On this figure, the vertical displacement on the free part of the strip was increased fortyfold to make the buckling phenomenon visible on the figure.

In Figure 9 one sees the vertical displacement at the position x = 470 mm for different values of λF . For weak value of λF (i.e 8 MPa) the strip is slightly curved. As λF exceeds 13.1 MPa the buckling appears. For larger values of λF (i.e 44 MPa and 120 MPa), the buckling is well established but the distortion remains small. The strip profile and the start of the buckling are in good agreement with the results obtained by Gueydan [START_REF] Gueydan | Modélisation numérique du flambage de bandes en acier dans les recuits continus[END_REF] where the buckling initiates at λF =12.9 MPa.

Towards the simulation of irreversible wrinkling in industrial processes

In [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF], sheet wrinkling in continuous annealing furnace was predicted in two steps by using a standard finite element code. In the first step, the combination of a traction load and of a roll profile led to an elastic buckling phenomenon. In a second step, a strip displacement was prescribed and the friction between roll and strip leads to an increase of the compressive stress and finally a wrinkle formation due to a plastic cumulative mechanism. We have established in the section 4.2.1 that the present solver is completely able to simulate the first step governed by geometric non-linearity and frictionless unilateral contact. The use of Automatic Differentiation was a crucial tool in a view to describe easily a general roll shape. In the same way, a general nonlinear constitutive law should be introduced, but we simply considered the linear Saint-Venant constitutive law (3) that is sufficient for a small strain buckling phenomenon. The relevance of the Saint-Venant constitutive law could be questioned, but not for such small strains. The framework of the present paper can be easily extended to hyperelasticity as in [START_REF] Galliet | Une version parallèle des méthodes asymptotiques numériques: Application à des structures complexes à base d'élastoméres[END_REF] [44] [START_REF] Nezamabadi | Solving hyperelastic material problems by asymptotic numerical method[END_REF], but this will be insufficient for the aforementioned second step of [START_REF] Jacques | Buckling and wrinkling during strip conveying in processing lines[END_REF]. For this second step, the plastic constitutive law and the friction law are not algebraic, but Differential Algebraic Equations (DAE). It is clear that the ideas of the present paper could also be applied to DAEs, but for simplicity we do not address this issue in this paper and we refer to [START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF] [START_REF] Béchet | Taylor series to solve friction problems[END_REF] for the treatment of plasticity and friction within ANM.

The small strain framework is sufficient in the study of strip conveying, but this is not the case for most of the forming processes that require large strain plasticity. Let us mention that ANM has been used with standard large strain plastic models [START_REF] Assidi | Méthode asymptotique numérique pour la plasticité[END_REF]. So an extension of the present paper to DAEs would lead to a non-linear solver for many forming processes. The key point would be the automatic treatment of DAEs large strain plasticity.

Concluding remarks

In this paper, we presented an automatic approach to implement non-linear problems arising from continuum mechanics and other fields of engineering sciences. We focussed on problems including non-linear laws expressed with implicit and explicit functions. To achieve an automatic solver, an extension of the DIAMANT approach has been proposed. Then a library solving a wide range of non-linear analytic problems is available.

Moreover coupling our library with an automated Finite Element tool such as the Fenics project [START_REF] Logg | Automating the finite element method[END_REF] could result in a real generic finite element solver and could be of great interest.

As applications, we solved non-linear unilateral contact problems modelled with Finite Element Method employing classical 3D or shell elements. This work is the first step to achieve the simulation of the strip conveying process. The next steps consist in accounting friction phenomena and including plasticity models.

This work has been supported by the French "Agence Nationale de la Recherche" (ANR) in the framework of the Project INSTABANDE.

  It is then possible to compute the function value in the continuous domain ✞ ☎ here suppose that T_Diff is the class describing derivative types T_Diff computeDerivatives(T_Diff x)

4. 2

 2 Ability of the solver to simulate strip conveying and other forming processes 4.2.1 Simulation of an elastic buckling problem under traction and unilateral contact with a tapered roll

  are the dimensions of the strip (horizontal length L = 1000 mm, width l b 2 = 200 mm, thickness = 0.2 mm, table length

lt 2 =Figure 1 :Figure 2 :Figure 3 : 20 Figure 4 :Figure 5 :

 21232045 Figure 1: Cylinder/sheet geometry

Figure 6 :

 6 Figure 6: Vertical displacement profile along a conical cylinder. ANM versus Abaqus results.

Figure 7 :

 7 Figure 7: Industrial geometry

Figure 8 :

 8 Figure 8: Zoom x40 on the strip buckling

Figure 9 :

 9 Figure 9: Vertical displacement observed at x=470 mm for different tensile load


Convex set detection

Keywords: change point, detection, hypothesis testing, minimax, sepration rate, set estimation

come

Introduction

Consider the statistical model:

Y i = 1(X i ∈ G) + ξ i , i = 1, . . . , n, (1)
where 1(•) is the indicator function. The set made of the points X i is called the design. The unknown set G is a segment on [0, 1], the noise terms ξ i are i.i.d. random variables, independent of the design. We distinguish two types of design:

(DD) Deterministic, and regular design: X i = i/n, i = 1, . . . , n;

(RD) Random, uniform design: the variables X i , i = 1, . . . , n, are i.i.d. uniform on [0, 1].

In the sequel, the design will be denoted by X . This model can be interpreted as a partial and noisy observation of an image, [0, 1], in which there is an unknown object G. From this observation, one would like to determine whether it is true that there is an object -the unknown set G might be empty -, and/or to recover that object, i.e. to estimate G. Another framework for this model is the noisy observation of some signal, here 1(X i ∈ G), and one would like to determine if the observation comes from a pure noise, or if there actually is some signal, and/or to recover that signal.

As it was mentioned in [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF], what makes difficult estimation of G is the complexity of the class of possible candidates, and detectability of G. The unknown set belongs to the class of all segments [a, b] of [0, 1]. This is a parametric class, with two parameters 0 ≤ a ≤ b ≤ 1. Therefore, one may expect to be able to estimate G at the parametric speed 1/n, up to a positive multiplicative constant. However, it turns out that G may be too hard to detect and thus, the speed of estimation would be deteriorated. For instance, if the segment is too small, and if no point of the design falls inside G, no procedure would see G. We would like to understand under which assumptions detection is not an obstacle for estimating G. In particular, the two following assumptions will be of interest for us : Assumption 1. The set G is of the form [0, θ], for some unknown number 0 ≤ θ ≤ 1.

Assumption 2. |G| ≥ µ.

Here and in the rest of the paper, µ ∈ (0, 1) is a given positive number. The brackets | • | stand for the Lebesgue measure.

Let us discuss these two assumptions in order to understand, intuitively, why they facilitate detection of G. The first one gives some information on the location of G. In other words, it tells that the set G starts from the left side of the frame. The second one tells that G is not too small, and thus should not be unnoticed by the statistician. Model (1), together with Assumption 1, is well known under the name of the change point problem. It can be rewritten as :

Y i = 1(X i ≤ θ) + ξ i , i = 1, . . . , n, (2)
for some number θ ∈ [0, 1]. The change point problem was studied in [KT93, Sec. 1.9], and a continuous-time version of this model is addressed in [START_REF] Korostelev | On minimax estimation of a discontinuous signal (translation of a paper in russian[END_REF].

The aim is to estimate the breakpoint θ. In the continuous time version, Korostelev [START_REF] Korostelev | On minimax estimation of a discontinuous signal (translation of a paper in russian[END_REF] proposed a more general framework. Instead of the indicator function in the regression equation (2), that is a function with a jump at the point θ, and satisfying a Lipschitz condition both on the left and on the right sides of θ. In these two works [START_REF] Korostelev | On minimax estimation of a discontinuous signal (translation of a paper in russian[END_REF], [START_REF] Korostelev | Minimax Theory of Image Reconstruction[END_REF], the change point θ is estimated with a precision, in expectation, of the order of 1/n, if θ is assumed to be separated from 0 and 1 : h ≤ θ ≤ 1-h, for some h ∈ (0, 1/2). Ibragimov and Khasminskii [START_REF] Ibragimov | Statistical Estimation: Asymptotic Theory[END_REF], under, among others, the same assumption of separation from 0 and 1, also proposed a consistent estimator of the discontinuity point of a regression function, with precision of order 1/n as well. This separation hypothesis is common in this kind of estimation problems. For instance, it is made in [START_REF] Korostelev | Asymptotically minimax image reconstruction problems[END_REF] and [KT93, Chap 3], where Tsybakov and Korostelev propose an estimator for boundary fragments. A boundary fragment, in dimension d ≥ 1, is a subset G of R d which can be described as the subgraph of a positive function g

: [0, 1] d-1 -→ [0, 1]: G = (x 1 , . . . , x d ) ∈ [0, 1] d : 0 ≤ x d ≤ g(x 1 , . . . , x d-1 ) . (3)
The authors study a similar model to Model (1), where G is a boundary fragment instead of a segment. This is one possible generalization in higher dimensions of our problem. They assume that the true g belongs to some Hölder class, and is separated from 0 and 1:

h ≤ g(x) ≤ 1 -h, ∀x ∈ [0, 1] d-1 ,
for a given parameter h ∈ (0, 1/2). In [START_REF] Korostelev | Asymptotically minimax image reconstruction problems[END_REF], the same authors estimate the support of a uniform density, assuming it is a boundary fragment. Again, they make the hypothesis that the underlying function g is separated from 0 and 1. In both models, they build a piecewise polynomial estimator of the function g, and prove that it is optimal in a minimax sense -see more details bellow -. One goal of the present work is to understand whether this separation hypothesis is necessary, in the simple case of dimension 1. Two other extensions in higher dimension have been addressed in [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF]. A segment of [0, 1] can be interpreted both as a one dimensional convex body, or as a one dimensional convex polytope. If the dimension is greater than 1, then the class of convex bodies is much bigger than that of convex polytopes. This is what explains that the optimal speed of estimation is better for polytopes than for general convex bodies -details in [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF] -. In this previous work, we did not study the impact of such assumptions as Assumptions 1 and 2. This question is beyond the scope of the present paper, but we believe that the results would be similar in the case of convex polytopes. Indeed, that case allows to keep the parametric property and, as we mentioned in [Bru13, Section 5], detectability becomes the main obstacle for estimation of the unknown set. In [START_REF] Gayraud | Minimax hypothesis testing about the density support[END_REF], a testing problem is addressed, in a different model: the set of interest, as in [START_REF] Korostelev | Asymptotically minimax image reconstruction problems[END_REF], is the support of a uniform density, assumed to be a boundary fragment (see (3)). Of course, detection is not pertinent in that model, since, as long as there are observations, there must be a nonempty support. However, the author studied the separation rate for distinguishing hypotheses of the type

H 0 : G = G 0 and H 1 : d(G, G 0 ) ≥ h,
where G 0 is a given boundary fragment, d is the Nykodim distance between sets, defined below, and h is a given positive number which may depend on n.

In [START_REF] Chan | Detection with the scan and the average likelihood ratio[END_REF], the detection question is addressed in a slightly different framework. In Model (1), we assume that the strength of the signal is given, equal to 1. It is also interesting to deal with the case of a signal of unknown strength, i.e.

Y i = δ1(X i ∈ G) + ξ i , i = 1, . . . , n, (4)
where δ is a positive number. For the signal to be detectable, there should be a tradeoff between its length |G| and its strength δ. Naturally, if δ is small, then the set G should be big enough and conversely, if δ is large, the set G is allowed to be small, so that the signal can be detected. Testing the presence of a signal, i.e. whether δ = 0 or not, is considered in [START_REF] Chan | Detection with the scan and the average likelihood ratio[END_REF]. This work is concerned with the power of two tests: the scan -or maximumlikelihood ratio, and the average likelihood ratio.

δ n|G| ≥ 2 ln 1 |G| + b n , for some sequence b n such that b n -→ ∞, then
there is a test whose power is asymptotically 1. Note that here, δ |G| is exactly the L 2 -norm of the signal, δ1(• ∈ G) 2 = δ |G|. In [START_REF] Lepski | Asymptotically exact nonparametric hypothesis testing in sup-norm and at a fixed point[END_REF], signals of unknown shape but known smoothness were considered. Exact minimax separation rates, in terms of the L 2 -norm of the signal, for distinguishing the null hypothesis, under which observations are pure noise, and the alternative one, under which there is a signal, are given. Detection is harder in that framework, because unlike in Model (1) or (4), where the shape of the signal is known -it is piecewise constant -, only its smoothness is known, and the separation rates are larger than those of models (1) and (4), in the sense that they allow less freedom for the size of the signal. However, this problem is different from ours, since we are concerned with the location of the signal, not the signal itself. Model 1 deals with change points in the mean of the observations, conditionally to the design. Under Assumption 1, there is only one change point, and under Assumption 2, there are two. A problem of interest in time series analysis is that of detecting change points in the mean, or in the covari-ance function, or other characteristics of the series. We refer to [START_REF] Shao | Testing for change points in time series[END_REF] and the references therein. In [FMSnt], a sample of n independent observations Y 1 , . . . , Y n is given, and one assumes that Y i admits a density f (•, ϑ( i n )) with respect to a given measure, for i = 1, . . . , n, where f belong to an exponential parametric class of densities. The real valued function ϑ is assumed to be piecewise constant on [0, 1], with a finite number K of jumps, not necessarily known. Under a similar condition on δ and |G| to that of [START_REF] Chan | Detection with the scan and the average likelihood ratio[END_REF], it is shown that at least one change point is consistently detectable: the proposed estimator of K is positive with probability that goes to 1, as n → ∞. When f (•, θ) is the density of a Gaussian distribution with mean θ and given variance σ 2 > 0, the problem was addressed in [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF]. This model includes models (1) and (4), when the design is (DD). However, Assumption 2 is not considered in that work, and optimality of the estimator of the change point is not treated, in the case when it is assumed to be unique -which corresponds to Assumption 1 -or when it is known that there are only two of them.

Notation

If G is a segment of [0, 1], we denote, respectively, by P G , E G and V G the joint probability measure of the observations ((X 1 , Y 1 ), . . . , (X n , Y n )) that come from Model (1), and the corresponding expectation and variance operators. If the event inside the probability, expectation or variance sign does not depend on G, we omit the subscript and write only P, E and V.

If G 1 and G 2 are two segments of [0, 1], we denote by G 1 △G 2 their symmetric difference. The Lebesgue measure of the symmetric difference of G 1 and G 2 , |G 1 △G 2 |, is also called the Nykodim distance between G 1 and G 2 .

In the present work, we do two types of inference on the underlying set G: hypothesis testing and estimation. The tests consist of deciding whether G is empty or not. Estimation consists of giving an estimator of G. A subset Ĝn of [0, 1], whose construction depends on the sample is called a set estimator or, more simply, an estimator. Given an estimator Ĝn , we measure its accuracy in a minimax framework. The risk of Ĝn on a class C of Borel subsets of R is defined as

R n ( Ĝn ; C) = sup G∈C E G [|G△ Ĝn |]. ( * )
The rate (a sequence depending on n) of an estimator on a class C is the speed at which its risk converges to zero when the number n of available observations tends to infinity. For all the estimators defined in the sequel, we are interested in upper bounds on their risk, in order to get information about their speed of convergence. The minimax risk on a class X , when n observations are available, is defined as

R n (C) = inf Ĝn R n ( Ĝn ; C), ( * * )
where the infimum is taken over all set estimators depending on n observations. If R n (C) converges to zero, we call the minimax rate of convergence on the class C the speed at which R n (C) tends to zero. It is interesting to provide a lower bound for R n (C). By definition, no estimator can achieve a better rate on C than that of the lower bound. This bound gives also information on how close the risk of a given estimator is to the minimax risk. If the rate of the upper bound on the risk of an estimator matches the rate of the lower bound on the minimax risk on the class C, then the estimator is said to have the minimax rate of convergence, or to be optimal in the minimax sense -up to constants -on this class. A test consists of deciding whether to reject or not a given hypothesis, called the null hypothesis, when it is compared to an alternative one. Let h ∈ (0, 1). In the whole paper, we will consider the following null hypothesis:

H 0 : G = ∅,
and the alternative hypothesis:

H 1 : |G| ≥ h.
Testing H 0 against H 1 is equivalent to deciding whether there is a set G or not. A test is a random variable τ n , which is built from the data set, and whose possible values are 0 and 1. The decision associated to the test τ n is to reject H 0 if and only if τ n = 1. We measure the performance of a test τ n on a class C using

γ n (τ n , C) = P ∅ [τ n = 1] + sup G∈C,|G|≥h P G [τ n = 0] .
This quantity is the sum of the errors of the first and the second kinds of the test τ n . We say that τ n is consistent on the class C if and only if γ n (τ n , C) -→ 0, when n → ∞. Let us allow the number h to depend on n. We call the separation rate on the class C any sequence of positive numbers r n such that:

• if h r n -→ n→∞ ∞
, then there exists a consistent test on C, and

• if h r n -→ n→∞ 0, then no test is consistent test on C.
With regard to both Assumptions 1 and 1, we focus on three different classes of sets, which are defined as bellow:

-S = {[a, b] : 0 ≤ a ≤ b ≤ 1} is the class of all segments on [0, 1], -S 0 = {[0, θ] : 0 ≤ θ ≤ 1} is the class of all segments on [0, 1], satisfying Assumption 1, -S(µ) = {G ∈ S : |G| ≥ µ} is the class of all segments on [0, 1], satisfying Assumption 2.
In the whole paper, we assume that the noise terms ξ 1 , . . . , ξ n from Model (1) satisfy

E e uξ i ≤ e σ 2 u 2 /2 , ∀u ∈ R, (5)
for some positive constant σ > 0. This constant need not be known. By satisfying this inequality, the noise is said to be subgaussian. An important case of such subgaussian variables is zero mean Gaussian variables, for which σ 2 turns to be the variance.

Since the design and the noise are assumed to be independent, reordering the X i 's does not modify the model. Indeed, there exists a reordering {i 1 , . . . , i n } of {1, . . . , n}, such that X i 1 ≤ . . . ≤ X in . The random indexes i 1 , . . . , i n are independent of the noise, and therefore the new noise vector (ξ i 1 ≤ . . . ≤ ξ in ) has the same distribution as (ξ 1 ≤ . . . ≤ ξ n ). With regard to this remark, we assume from now on that X is the reordering of a preliminary design, and therefore X 1 ≤ . . . ≤ X n almost surely, without loss of generality.

For two real valued sequences A n and B n , and a parameter ϑ, which may be multidimensional, we will write A n ≍ ϑ B n to say that there exist positive constants c(ϑ) and C(ϑ) which depend on ϑ, such that c(ϑ)B n ≤ A n ≤ C(ϑ)B n , for n large enough. If we put no subscript under the sign ≍, this means that the involved constants are universal, i.e. do not depend on any parameter.

Contributions of the paper

In this paper, we propose three new results.

1. The first result is a preliminary to the rest of the paper, and has been already mostly covered in the existing literature. Model (1) is a simple change point problem. If Assumption 1 is satisfied, the change point is unique. Otherwise, it is known that there are two change points, and the size of the jump is also known. In a hypothesis testing setup, we give the the minimal size -as a function of n -of the unknown set G, in Model (1), must satisfy, for that set to be detectable, depending on Assumption 1 being satisfied or not. We prove our results under both designs (DD) and (RD).

2. The second result concerns the change point problem is not necessary. We show that the speed 1/n of convergence in estimating the breakpoint θ in (2) can be achieved without assuming a ≤ θ ≤ b for some positive numbers 0 < a ≤ b < 1.

3. The last result shows that recovering a set requires that set to be sufficiently big, in order to achieve the parametric rate 1/n uniformly. We show that without this assumption, a logarithmic factor appears in the minimax rate.

Our work is made in the minimax setup. We provide exact rates of convergence of the minimax risks on the classes which we study. These rates are obtained by precise deviation inequalities for all the estimators that we propose. In the second and third results, the minimax rates correspond exactly to the separation rates given by the first result, for testing emptiness of the set G.

Organization of the paper

This paper is organized as follows. Next section deals with the question of deciding whether there is a set or not. We give the separation rates for the corresponding detection test, for the change point problem -i.e. under Assumption 1 -and on the class S. The third section deals with estimating the unknown set, and we give the minimax rates of estimation on the classes S, S 0 and S(µ). The forth section is dedicated to the conclusion, followed by a discussion about possible extensions. The fifth and last section is devoted to the proofs.

Detection of a set

We test the null hypothesis H 0 : G = ∅, against the alternative hypothesis H 1 : |G| ≥ h. We find the asymptotic minimal value of h, as a function of n, so the two hypotheses are well separated and there exists a consistent test. This minimal value depends on the class of possible sets G. Intuitively, the a priori knowledge that the unknown set G belongs to the class S 0 gives an important piece of information about the location of this set, and therefore it should be easier to detect it. Actually, the following theorem confirms this intuition, by showing that the separation rate is smaller -by a logarithmic factor -for the subclass S 0 than for the whole class S.

The idea, for the class S 0 , is the following. Under H 1 , [0, h] ⊆ G. Therefore, we check among those observations (X i , Y i ) for which X i ≤ h if there is a sufficiently large number of Y i 's that are large, e.g. larger than 1/2. Let N = max{i = 1, . . . , n :

X i ≤ h} = # (X ∩ [0, h]).
Let S be the following test statistics:

S = #{i = 1, . . . , N : Y i ≤ 1 2 }.
If the alternative hypothesis holds, i.e. if |G| ≥ h, all the X i , i = 1, . . . , N fall inside the set G, and the corresponding Y i should not be too small. The test statistic S counts how many of these Y i 's are suspiciously small. This is how is built the test T 0 n :

T 0 n = 1(S ≤ cN ),
where c

∈ (P[ξ 1 ≤ -1/2]), P[ξ 1 ≤ 1/2]) .
For the class S, we propose a scan test, i.e., a procedure which scans the whole frame [0, 1] and seeks for a large enough quantity of successive observations for which

Y i is large. If G ∈ S, let R(G) = n i=1 Y i 1(X i ∈ G) -#(X ∩G)

2

, and R = sup |G|≥h R(G). Under the alternative hypothesis, R should be quite large, and we define the test

T 1 n = 1(R ≥ 0).
Theorem 1. Let Model (1) hold.

1. Assume that the design is (DD) or (RD), and that the noise satisfies:

P[ξ 1 ≤ -1/2] < P[ξ 1 ≤ 1/2].
Then, if nh -→ ∞, the test T 0 n is consistent, i.e. γ n (T 0 n , S 0 ) -→ 0. If, in addition, the noise is Gaussian, then a separation rate on the class S 0 is r n = 1/n.

2. Assume that the design is (DD) or (RD). Then, if nh/ ln n -→ ∞, the test T 1 n is consistent, i.e. γ n (T 1 n , S) -→ 0. If, in addition, the noise is Gaussian, then a separation rate on the class S is r n = ln(n)/n.

In the next section, we show that the separation rates given in Theorem 1 are the minimax rates of convergence on the corresponding classes.

3 Estimation of a set 3.1 Least square estimators

Let Model (1) hold. For G ′ ∈ S, let A 0 (G ′ ) = n i=1 (Y i -1(X i ∈ G ′ ))
2 be the sum of squared errors. A way to estimate G is to find a random set Ĝn which minimizes A 0 (G ′ ), among all possible candidates G ′ . Note that minimizing A 0 (G ′ ) is equivalent to maximizing

A(G ′ ) = n i=1 (2Y i -1)1(X i ∈ G ′ ). (6)
Denote by S = {i = 1, . . . , n : X i ∈ G} and by

S ′ = {i = 1, . . . , n : X i ∈ G ′ },
for some G ′ ∈ S. Denote by #(•) the cardinality, for finite sets. The criterion A(G ′ ) becomes, if denoted as a function of S ′ ,

A(S ′ ) = i∈S ′ (2Y i -1) = i∈S ′ (21(X i ∈ G) + 2ξ i -1) = 2#(S ∩ S ′ ) -#S ′ + 2 i∈S ′ ξ i , so, A(S ′ ) -A(S) = -#(S△S ′ ) + 2   i∈S ′ \S ξ i - i∈S\S ′ ξ i   . (7)
A subset S ′ of {1, . . . , n} is called convex if and only if it is of the form {i, . . . , j}, for some 1

≤ i ≤ j ≤ n. It is clear that if a convex subset S ′ of {1, . . . , n} maximizes A(S ′ ) -A(S), then by defining G ′ = [X min S ′ , X max S ′ ], the segment G ′ maximizes A(G ′ ) (cf. (6)).

Estimation of the change point

P G | Ĝn △G| ≥ x n ≤ C 0 e -x/(8σ 2 ) , ∀x > 0,
where C 0 is a positive constant which depends on σ only.

A simple application of Fubini's theorem leads to the following result.

Corollary 1. Let the assumptions of Theorem 2 be satisfied. Then, for all q > 0, there exists a positive constant A q which depends on q and σ only, such that sup

G∈S 0 E G | Ĝn △G| q ≤ A q n q .
This corollary shows that the minimax risk on the class S 0 is bounded from above by 1/n, up to multiplicative constants. Next theorem proves that up to constants, 1/n is also a lower bound on the minimax risk, if the noise is Gaussian.

Theorem 3. Consider Model (1), with design (DD). Then for all integer n ≥ 1,

R n (S 0 ) ≥ 1 2n .
Combining Theorems 2 and 3 yields Theorem 4. Consider Model (1), with design (DD). Then, the minimax risk on the class S satisfies

R n (S 0 ) ≍ σ 1 n .

Recovering any set

Let us now assume that the unknown set G does not necessarily contain 0.

We shall prove that whether to assume that G belongs to the class S(µ) or not does not lead to the same minimax rate. As we saw in Section 3.1, an estimator of G in Model (1) can be obtained by maximizing the Gaussian process (7) over all segments of {1, . . . , n}. This is not the track that we will borrow, but it would be interesting to work precisely on this process. This would probably be the first step to extensions of our results in higher dimensions. However, this problem remains open for now. The methods that we develop in this section are quite different. If G is only assumed to belong to the biggest class S, the proposed estimator is the LSE, which was already detailed in [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF] for convex polytopes, in higher dimension. If |G| is a priori known to be greater or equal to µ, then we first build a preliminary estimator of G -the LSE -, using one half of the observed sample. This estimator is not optimal, but it is close to G with high probability. We show that the middle point mn of this estimator is in G with high probability. This brings us back to the change point problem, where 0 is now replaced by mn , and we use the second half of the observed sample to estimate two change points. Let us first state the following theorem, which is, for the design (RD), a particular case of [Bru13, Theorem 1], for d = 1.

Theorem 5. Let n ≥ 2. Let Model Model hold, with design (DD) or (RD). Let Ĝn ∈ ArgMax G ′ ∈S A(G ′ ) be a LSE estimator of G. Then, there exist two positive constants C 1 and C 2 which depend on σ only, such that

sup G∈S P G n | Ĝn △G| - 4 ln n C 2 n ≥ x ≤ C 1 e -C 2 x , ∀x > 0.
The expressions of C 1 and C 2 are given in the proof of [Bru13, Theorem 1], for the design (RD). For the design (DD), we do not give a proof of this theorem here, but it can be easily adapted from that of the case of the design (RD). The next corollary is immediate.

Corollary 2. Let the assumptions of Theorem 5 be satisfied. Then, for all q > 0, there exists a positive constant B q which depends on q and σ only, such that

sup G∈S E G | Ĝn △G| q ≤ B q ln n n q .
This corollary shows that the minimax risk on the class S is bounded from above by ln(n)/n, up to a multiplicative constant. The following theorem establishes a lower bound, if the noise is supposed to be Gaussian. Theorem 6. Consider Model (1), with design (DD) or (RD). Assume that the noise terms ξ i are i.i.d. Gaussian random variables, with variance σ 2 > 0. For any large enough n,

R n (S) ≥ α 2 σ 2 ln n n ,
where α is a universal positive constant.

This lower bound comes from [Bru13, Theorem 2] in the case of the design (RD), and the proof is easily adapted for the design (DD). Eventually, the minimax risk on the class S is of the order ln(n)/n: Theorem 7. Consider Model (1), with design (DD) or (RD). Assume that the noise terms ξ i are i.i.d. Gaussian random variables, with variance σ 2 > 0. The minimax risk on the class S satisfies, asymptotically: R n (S) ≍ σ ln n n .

For the design (DD), we combine both Theorems 2 and 5 to find the minimax rate on the class S(µ). Let Model (1) hold, and let G ∈ S(µ). First, we split the sample into two equal parts. Let D 0 be the set of sample points with even indexes, and D 1 the set of sample points with odd indexes. Note that D 0 ∪D 1 is exactly the initial sample, that these two subsample are independent, and that each of them is made of at least (n -1)/2 data. Let Ĝn be the LSE estimator of G given in Theorem 5, built from the subsample D 0 . Let mn be the middle of Ĝn . As it will be shown in the proof of the next theorem, mn satisfies both following properties, with high probability:

1. mn ∈ G, 2. µ/2 ≤ mn ≤ 1 -µ/2.
This brings us to estimating two change points -the endpoints of G -, using the second subsample D 1 . From Theorem 2, we know that this can be done at the speed 1/n, up to multiplicative constants.

Theorem 8. Consider Model (1), with design (DD). There exists an estimator Gn of G, such that sup G∈S(µ)

P G | Gn △G| ≥ x n ≤ 2C 0 e -µx/(256σ 2 ) + C 1 n 4 e -C 2 µn/2 , ∀x > 0,
for n large enough. The positive constants C 0 and C 2 appeared in Theorems 2 and 5.

Naturally, Theorem 8 leads to the next corollary.

Corollary 3. Let the assumptions of Theorem 8 be satisfied. Then, for all q > 0, there exists a positive constant B ′ q which depends on q, µ and σ only, such that sup

G∈S(µ) E G | Gn △G| q ≤ B ′ q n q .
This corollary, for q = 1, shows that the minimax risk on the class S(µ) is bounded from above by 1/n, up to a multiplicative constant. A very similar proof to that of Theorem 3 yields a lower bound for this minimax risk, which leads to the next theorem.

Theorem 9. Consider Model (1), with design (DD). The minimax risk on the class S(µ) satisfies, asymptotically,

R n (S(µ)) ≍ µ,σ 1 n .
Remark 1. Note that, in Theorem 8, the upper bound contains one residual term which does not depend on x. This term, in order to be sufficiently small, requires that µ -if allowed to depend on n -is of larger order than ln(n)/n. This in an echo to Theorem 1, in which we showed that the smallest set which can be detected has measure of this order exactly. In addition, if µ is of the order of ln(n)/n, then the proof of the lower bound of Theorem 6 can be applied, and the minimax risk on the class S(µ) will satisfy R n (S(µ)) ≍ σ ln n n .

Conclusion and discussion

We summarize our results in Table 1. The rates that are written in this table hold for Gaussian noise, which is the most important case. In each case we indicate the design for which the rate holds.

S S 0 S(µ) Minimax rate ln(n)/n (DD,RD)

1/n (DD) 1/n (DD) Separation rate ln(n)/n (DD,RD) 1/n (DD,RD)

•

Table 1: Minimax risks and separation rates for the classes S, S 0 and S(µ).

Note that in two cases, only the design (DD) has been considered. This is mainly for technical reasons, and we believe that the rates are still the same for the design (RD).

It comes out that asymptotically, a segment can be estimated infinitely faster when it is a priori supposed either to contain a given point (here, 0), or to be large enough. The main question that remains is: does this phenomenon still hold for two -or higher -dimensional sets ? An important -if not essential -assumption which has been done all over this paper is the convexity of the unknown set. In dimension 1, the class of convex subsets of [0, 1] is simple, and parametric. In any higher dimension d, the class C d of convex subsets of [0, 1] d is much more complex. In particular, its metric entropy is much larger than that of a parametric family [START_REF] Bronshtein | ǫ-entropy of convex sets and functions[END_REF], and it seems to us that it is the complexity of this class that makes it harder to estimate a set, than detectability.

In Model (1), if G belongs to C d , estimation of G can be done at the minimax rate n -2/(d+1) [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF]. However, an adaptation of the proof of [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF]Theorem 2] shows the following: on any subclass of C d invariant under translations and invertible affine transformations -which keep a set G inside the frame [0, 1] d -, the minimax rate is at least of order ln(n)/n. In addition, we believe this is the separation rate for the detection problem, on any such subclass of C d .

For the whole class C d , since ln(n)/n is much smaller than n -2/(d+1) , the minimax rate is of the order of n -2/(d+1) . However, for a parametric subclass, such as that of all convex polytopes with a given number of vertices, we believe that the complexity of the class leads to a term of order 1/n in the minimax risk, although detectability leads to a term of order ln(n)/n. In our opinion, this is what explains that, as shown in [START_REF] Brunel | Adaptive estimation of convex polytopes and convex sets from noisy data[END_REF], the minimax risk on the class of all convex polytopes of [0, 1] d , with a given number of vertices, is of order ln(n)/n.

Motivated by Theorem 8, we also conjecture that 1/n is the minimax rate on the class of all convex polytopes of [0, 1] d , of volume greater than a given µ > 0. However, it is not possible to extend Theorem 8 to higher dimensions. An adaptation of the proof of [Bru13, Theorem 2], by taking, as the M hypotheses used in the proof, sets which contain the origin and have pairwise zero measure intersections, would show that (ln n)/n remains a lower bound under that assumption. Yet, we believe that if the unknown set is assumed to contain a given section of positive d -1 dimensional Lebesgue measure, of a given hyperplane in R d , then an analog of Theorem 8 should hold, and the minimax rate should be of order 1/n.

To end this discussion, we note that Model (4) can be rewritten in the vector form:

y = β + ξ, (8)
where

y = (Y 1 , . . . , Y n ) ⊤ , β = (δ1(X 1 ∈ G), . . . , δ1(X n ∈ G)) ⊤ ∈ {0, δ} n and ξ = (ξ 1 , . . . , ξ n ) ⊤ .
If we denote by A the n × n lower triangular matrix with coefficients on and under the main diagonal equal to 1, equation ( 8) can be rewritten as:

y = Aγ + ξ, (9)
where γ = (β 1 , β 2β 1 , . . . , β nβ n-1 ) ⊤ is the vector of differences of β. The vector γ is sparse: it has exactly two nonzero coefficients, corresponding to the change points locations. This formulation of our initial problem leads to a high dimensional regression setup, under sparsity constraints. If δ is known, equal to 1 -this is Model (1)-, we are not interested in the estimation of the vector γ itself, but only of its support, which indicates the location of G, if the design is (DD). If δ is unknown, the two nonzero coefficients of γ take values δ and -δ, and it is also of interest to estimate γ itself, in addition to its support. Penalized regression methods, with Lasso-type penalizations (see [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] for details about Lasso), have been addressed in [START_REF] Mammen | Locally adaptive regression splines[END_REF],

[BKL + 09], among others. In these two works, the penalization is written in terms of the number of jumps in the vector β in (8), which is equal, for Model (1), to the number of nonzero coefficients of γ in (9). Therefore, the penalty is equivalent to a L 0 or BIC-type penalty, see for example [START_REF] Bunea | Aggregation for gaussian regression[END_REF]. Estimation of the support is usually a secondary question, which is addressed in terms of consistency, and not of rate of convergence. The main focus is about estimation of γ, or prediction, i.e., estimation of β. We believe that the parameter µ, if both µ and G are unknown, can be estimated at the speed (ln n)/n. In that case, some results about estimation rates of G, for the design (DD), are given in [BKL + 09], but the risks, in expectation, are not computed.

It would be interesting to understand how Lasso or BIC-type estimators could be adapted in order to achieve the rate 1/n on the classes S 0 and S(µ). The penalties of the classical estimators do not use the information that is contained in those classes, and we believe that these estimators are suboptimal. However, one should try to modify the penalties, according to that information, in order to improve the accuracy of these estimators.

Proofs

Proof of Theorem 1

On the class S 0 Upper bound Let us first prove the upper bound, i.e. assume that nh → ∞, and prove that there exists a consistent test. Recall that N = max{i = 1, . . . , n :

X i ≤ h} = # (X ∩ [0, h]).
If the design is (DD), then N is just equal to the integer part of nh. If the design (RD), then N is a binomial random variable, with parameters n and h. Let us show first that the error of the first kind of the test T 0 n goes to zero, when n → ∞.

P ∅ [S ≤ cN ] = P ∅ # i = 1, . . . , N : Y i > 1 2 ≥ (1 -c)N ≤ E P ∅ # i = 1, . . . , N : ξ i > 1 2 ≥ (1 -c)N |X .
Since the ξ i 's are independent of X , the distribution of #{i = 1, . . . , N : ξ i > 1 2 } conditionally to X is binomial, with parameters N and β, where β = P [ξ 1 > 1/2] ∈ [0, 1). Thus, by Bernstein's inequality for binomial random variables, by defining γ =

(1 -c -β) 2 2β(1 -β) + (1 -c -β)/3 > 0, P ∅ [S ≤ cN ] ≤ E [exp (-γN )] .
If X satisfies (DD), then N ≥ nh -1 and it is clear that

P ∅ [S ≤ cN ] -→ 0. If X satisfies (RD), then E [exp (-γN )] = exp -nh 1 -e -γ , so P ∅ [S ≤ cN ] -→ 0.
Let us show, now, that the error of the second kind goes to zero as well. Let G ∈ S 0 satisfying the alternative hypothesis, i.e. |G| ≥ h. Denote by

β ′ = P[ξ 1 ≤ -1/2], and by γ ′ = (c -β ′ ) 2 2β ′ (1 -β ′ ) + (c -β ′ )/3 > 0 P G [S > cN ] = P G # i = 1, . . . , N : Y i ≤ 1 2 > cN ≤ E P ∅ # i = 1, . . . , N : ξ i ≤ - 1 2 > cN |X ≤ E exp -γ ′ N ,
by a similar computation to that for the error of the first kind. Since the right-side of the last inequality does not depend on G,

sup |G|≥h P G [S > cN ] ≤ E exp -γ ′ N
and therefore, by the same argument as for the error of the first kind, goes to zero when n → ∞, for both designs (DD) and (RD).

Lower bound Assume, now, that nh → 0. Let τ n be any test. Let G 1 = [0, h]. We denote by H the Hellinger distance between probability measures. The following computation uses properties of this distance, which can be found in [START_REF] Tsybakov | Introduction to nonparametric estimation[END_REF].

γ n (τ n , C) ≥ E ∅ [τ n ] + E G 1 [1 -τ n ] ≥ min (dP ∅ , dP G 1 ) ≥ 1 2 1 - H(P ∅ , P G 1 ) 2 2 . ( 10
)
Let G, G ′ ∈ S. A simple computation shows that, for the design (DD),

1 - H(P G , P G ′ ) 2 = exp - # (X ∩ (G△G ′ )) 8σ 2 , ( 11
)
and for the design (RD),

1 - H(P G , P G ′ ) 2 = 1 -1 -e -1 8σ 2 |G△G ′ | n . (12)
In particular, for the design (DD),

1 - H(P ∅ , P G 1 ) 2 ≥ exp - nh 8σ 2 ,
and for the design (RD),

1 - H(P ∅ , P ⊗n G 1 ) 2 = 1 -1 -e -1 8σ 2 h n .
In both cases, we showed that the right side of (10) tends to 1/2, when n → ∞. Therefore the test τ n is not consistent.

On the class S

Upper bound Assume that nh ln n -→ ∞. Let us first show that the error of the first kind of T 1 n goes to zero, when n → ∞. Recall that

T 1 n = 1(R ≥ 0), where R = sup |G|≥h R(G) and R(G) = n i=1 Y i 1(X i ∈ G) - #(X ∩G) 2
, for all G ∈ S. Note that R(G) is piecewise constant, and can only take a finite number of values. It is clear that

{R(G) : G ∈ S, |G| ≥ h} = {R([X k , X l )) : 1 ≤ k < l ≤ n, X l -X k > h}. Recall that for 1 ≤ k < l ≤ n, R([X k , X l )) = 1 2 l-1 i=k (2Y i -1). Therefore, P ∅ [R ≥ 0] = P ∅   max 1≤k<l≤n X l -X k >h R([X k , X l )) > 0   ≤ P ∅   1≤k<l≤n {R([X k , X l )) > 0} ∩ {X l -X k > h}   ≤ 1≤k<l≤n P ∅ [R([X k , X l )) > 0, X l -X k > h] ≤ 1≤k<l≤n P ∅ l-1 i=k (2ξ i -1) > 0 P[X l -X k > h]. (13)
For 1 ≤ k < l ≤ n,

P ∅ l-1 i=k (2ξ i -1) > 0 ≤ exp - (l -k)σ 2 8 , (14)
using Markov's inequality and (5).

If the design is (DD), then P[X l -X k > h] is 1 if and only if lk > nh, 0 otherwise, so from ( 13) and ( 14) we get that

P ∅ [R ≥ 0] ≤ l-k>nh exp - (l -k)σ 2 8 ≤ l-k>nh exp (-nh)σ 2 8 ≤ n 2 2 exp (-nh)σ 2 8 -→ 0,
when n → ∞, which proves that the error of the first kind goes to zero. If the design is (RD), let us use the following Lemma.

Lemma 1. Let X 1 , . . . , X n be the (RD) design. Then, for any 1 ≤ k < l ≤ n, and h > 0,

P[X l -X k > h] ≤ n exp -nh(1 -e -u ) + u(l -k) , ∀u > 0. Proof of Lemma 1 Note that the event {X l -X k > h} is equivalent to {#(X ∩ (X k , X k + h)) < l -k}.
Let us denote by X ′ 1 , . . . , X ′ n the preliminary design, from which X 1 , . . . , X n is the reordered version. The random variables X ′ 1 , . . . , X ′ n are then i.i.d., with uniform distribution on [0, 1]. Hence,

P[X l -X k > h] = n j=1 P #(X ∩ (X k , X k + h)) < l -k, X k = X ′ j ≤ n j=1 P #(X ∩ (X ′ j , X ′ j + h)) < l -k ≤ n j=1 E P #(X ∩ (X ′ j , X ′ j + h)) < l -k|X ′ j ≤ n j=1 E P n -#(X ∩ (X ′ j , X ′ j + h)) ≥ n -l + k|X ′ j ≤ n j=1 E f (X ′ j ) , (15)
where f

(x) = P [n -#(X ∩ (x, x + h)) ≥ n -l + k], for all x ∈ [0, 1]. The random variable n -#(X ∩ (x, x + h))
is binomial with parameters n and 1h, and by Markov's inequality, for all u > 0,

f (x) ≤ e nu 1 -h(1 -e -u ) n e -u(n-l+k) , (16)
and ( 15) and ( 16) yield the lemma.

Therefore, by ( 13), ( 14) and (15), and Lemma 1 with u = σ 2 /8,

P ∅ [R ≥ 0] ≤ 1≤k<l≤n exp - (l -k)σ 2 8 -nh(1 -e -u ) + u(l -k) ≤ n 2 2 exp -nh(1 -e -σ 2 /8 ) -→ 0,
when n → ∞, which proves that the error of the first kind goes to zero. Let us bound, now, the error of the second kind. Let G ∈ S satisfying |G| ≥ h. For this G, denote by N G = #(X ∩ G). Then,

P G [R < 0] ≤ P G [R(G) ≤ N G /2] ≤ P n i=1 ξ i 1(X i ∈ G) ≤ -N G /2 . (17)
For the design (DD), N G is the integer part of n|G|, so N G ≥ nh. Therefore, by Markov's inequality, and by ( 5), (17) becomes

P G [R < 0] ≤ exp - N G 8σ 2 ≤ exp - nh 8σ 2 . ( 18
)
For the design (RD), N G is a random binomial variable, with parameters n and |G|. By conditioning to the design and using Markov's inequality, (17) becomes

P G [R < 0] ≤ P n i=1 -ξ i 1(X i ∈ G) ≥ N G /2 ≤ E exp - N G 8σ 2 ≤ exp (-Cn|G|) ≤ exp (-Cnh) , (19)
where C = 1e -1 8σ 2 . In both cases (18) and (19), the right side does not depend on G, and goes to zero as n → ∞. We conclude that, for both designs (DD) and (RD),

sup |G|≥h P G [R < 0] -→ 0,
which ends the proof of the upper bound.

Lower bound We more or less reproduce the proof of [START_REF] Gayraud | Minimax hypothesis testing about the density support[END_REF], Theorem 3.1. Here, the noise is supposed to be Gaussian, with variance σ 2 . Let us assume that nh ln n -→ 0. Let M = 1/h, assumed to be an integer, without loss of generality. For q = 0, . . . , M , let G q = [qh, (q+1)h]. For q = 1, . . . , M , let Z q = dP Gq dP ∅ (X 1 , Y 1 , . . . , X n , Y n ), and denote by Z = 1 M M q=1 Z q . Let τ n be any test. Then,

γ n (τ n , S) ≥ P ∅ [τ n = 1] + 1 M M q=1 P Gq [τ n = 0] ≥ 1 M M q=1 P ∅ [τ n = 1] + P Gq [τ n = 0] ≥ 1 M M q=1 E ∅ [τ n ] + E Gq [1 -τ n ] ≥ 1 M M q=1 E ∅ [τ n + (1 -τ n )Z q ] ≥ E ∅ τ n + (1 -τ n ) Z ≥ E ∅ τ n + (1 -τ n ) Z 1( Z ≥ 1/2) ≥ 1 2 P ∅ Z ≥ 1/2 . ( 20
)
Let us prove that E ∅ [ Z] = 1, and that V ∅ [ Z] -→ 0. This will imply that the right side term of (20) goes to zero, when n → ∞. For q = 1, . . . , M , under the null hypothesis,

Z q = exp - 1 2σ 2 n i=1 (Y i -1(X i ∈ G q )) 2 -Y 2 i = exp 1 2σ 2 n i=1 (2ξ i -1)1(X i ∈ G q ) . (21)
By its definition, Z q has expectation 1 under P ∅ :

E ∅ [ Z] = 1. (22)
Since, almost surely, no design point falls in two G q 's at the time, a simple computation shows that the random variables Z q , q = 1, . . . , M , are not correlated. Thus,

V ∅ [ Z] = 1 M 2 M q=1 V ∅ [Z q ].
Let us bound from above V ∅ [Z q ], for q = 1, . . . , M :

V ∅ [Z q ] ≤ E ∅ [Z 2 q ] = E exp - #(X ∩ G q ) σ 2 E ∅ exp 2 σ 2 n i=1 ξ1(X i ∈ G q ) |X = E exp #(X ∩ G q ) σ 2 . ( 23
)
If the design is (DD), then we get that

V ∅ [Z q ] ≤ exp nh + 1 σ 2 ,
and the variance of Z is then bounded from above:

V ∅ [ Z] ≤ h exp nh + 1 σ 2 . ( 24
)
If the design is (RD), then #(X ∩ G q ) is a binomial random variable with parameters n and h, so from (23), we get that

V ∅ [Z q ] ≤ 1 + e 1/σ 2 -1 h n ≤ exp (Cnh) ,
where C = e 1/σ 2 -1, and the variance of Z is then bounded from above:

V ∅ [ Z] ≤ h exp (Cnh) . ( 25
)
Since we assumed that nh/ ln n -→ 0, the right side terms of ( 24) and (25) go to zero, and therefore, for both designs (DD) and (RD),

V ∅ [ Z] -→ 0. ( 26
)
Finally, we get from (20), ( 22) and (26), that liminf n→∞ γ n (τ n , S) ≥ 1 2 .

This concludes the proof.

Proof of Theorem 2

The beginning of this proof holds for any design {X 1 , . . . , X n }, independent of the noise ξ i , i = 1, . . . , n. Let G ∈ S 0 . Let M = max{i = 1, . . . , n : X i ∈ G} -set M = 0 if the set is empty -. Then, Mn ∈ ArgMax M ′ =1,...,n

(F (M ′ ) -F (M )), and, by (7),

F (M ′ ) -F (M ) = -|M ′ -M | +          2 M i=M ′ +1 ξ i if M > M ′ , 0 if M ′ = M, -2 M ′ i=M +1 ξ i if M < M ′ .
Let us complete the i.i.d. sequence ξ 1 , . . . , ξ n to obtain an infinite double sided i.i.d. sequence (ξ i ) i∈Z , independent of the design. Let k ∈ N * be any positive integer. Define, for i ∈ Z, ξi = ξ i+M . Since M depends on the design only, it is independent of the ξ i , i ∈ Z, and therefore, the ξi , i ∈ Z are i.i.d., with same distribution as ξ 1 . Let E k be the event { Mn ≥ M + k}. If and, by choosing u = 1/(4σ 2 ),

P G [E k ] ≤ Ce -k/(8σ 2 ) ,
where C = 1e -1/(8σ 2 ) -1 is a positive constant. By symmetry, we obtain that :

P G [| Mn -M | ≥ k] ≤ 2Ce -k/(8σ 2 ) . ( 27
)
If the design is (DD), the conclusion is straightforward, since for all i, j = 1, . . . , n, |X i -X j | = |i-j| n , and Theorem 2 is proved.

If the design is (RD), it is not clear how to go from (27) to an upper bound for the probability P G [| θn -θ| ≥ ǫ], for ǫ > 0. This question should be addressed in a future work.

Proof of Theorem 3

The proof is straightforward. Let G 1 = [0, 0] and G 2 = [0, 1/(2n)]. Then P G 1 = P G 2 , since no point of the design falls in G 1 △G 2 , and for any estima-tor Ĝn , sup

G∈S 0 E G | Ĝn △G| ≥ E G 1 | Ĝn △G 1 | + E G 2 | Ĝn △G 2 | ≥ E G 1 | Ĝn △G 1 | + | Ĝn △G 2 | ≥ E G 1 [|G 1 △G 2 |]

Proof of Theorem 8

Let I 0 be the set of even positive integers less or equal to n, and I 1 the set of odd such integers. Note that {X i : i ∈ I 0 } is a deterministic and regular design, with step 2/n. Let Ĝn ∈ ArgMax G ′ ∈S i∈I 0 (2Y i -1)1(X i ∈ G ′ ) be the LSE estimator given in Theorem 5, built using only the subsample {X i : i ∈ I 0 }. Let x > 0, whose value will be specified in the course of the proof. Consider the event E x = {| Ĝn △G| < x ln n n }. By Theorem 5, this event holds with probability at least 1-C 1 e -(C 2 x-4) ln n . Choose x such that

Let

  Model (1) hold, with design (DD). Assume that G belongs to S 0 . This is the change point problem. For some θ ∈ [0, 1], G can be written as G = [0, θ]. Let us make one preliminary remark. For any estimator Ĝn of G, the random segment Gn = [0, sup Ĝn ] performs better than Ĝn , since | Gn △G| ≤ | Ĝn △G| almost surely. Therefore, it is sufficient to consider only estimators of the form Ĝn = [0, θn ], where θn is a random variable. Then, | Ĝn △G| = | θn -θ|, and the performance of the estimator Ĝn of G is that of the estimator θn of θ. Let us build a least square estimator (LSE) of θ. For M = 1, . . . , n, let F (M ) = A({1, . . . , M } = M i=1 (2Y i -1). Let Mn ∈ ArgMax M =1,...,n F (M ), and θn = X Mn . The following theorem follows. Theorem 2. Let n ≥ 1. Let Model (1) hold, with design (DD). Let Ĝn = [0, θn ]. Then, sup G∈S 0

  E k holds, then :0 ≤ F ( Mn ) -F (M ) = M -Mnall u > 0, P G [E k ] ≤ P G max k≤j u+2σ 2 u 2 )j , by (5)

  x ln n n ≤ µ/2. This choice implies that on the event E x , | Ĝn △G| < µ ≤ |G|, so necessary, Ĝn and G must intersect. Thus, still on the event E x ,| Ĝn △G| = | bn -b| + |â n -a|,where we denoted by G = [a, b] and Ĝn = [â n , bn ]. Let m = a+b 2 and mn = ân+ bn 2 be, respectively, the middle points of G and Ĝn . From now on, let us assume that E x holds. Then, | mn -m| ≤ 1 2 (| bn -b| + |â n -a|) ∈ G and, combining (28) with the fact that |G| ≥ µ, min( mn , 1mn ) ≥ µ 4 . (29)

Let us define

 andI - 1 = {i ∈ I 1 : X i ≤ mn }. By (29), #I ǫ 1 ≥ µn 8 -1 ≥ µn 16 for n large enough, and for ǫ ∈ {+, -}. Note that {X i :

) is a deterministic and regular design of the segment [ mn , 1] (resp. [0, mn ]), of cardinality greater or equal to µn 16 , as we saw just before. Then, since we have both

, the change points a and b can be estimated as in Theorem 2, using the subsamples {(X i , Y i ) : i ∈ I + 1 } and {(X i , Y i ) : i ∈ I - 1 } respectively, and we get two estimators ãn and bn which satisfy:

and

for all y > 0. Set Gn = [ã n , bn ], on the event E x , and Gn = ∅ on its complementary Ēx . By setting x = µn 2 ln n , which is the maximal value authorized in this proof,

≤ 2C 0 e -µy/(256σ 2 ) + C 1 n 4 e -C 2 µn/2 , which ends the proof of Theorem 8.