Alain Gerard

Olivier Cahuc

Miron Zapciu

SELF-EXCITED VIBRATIONS IN TURNING: EFFORTS TORSOR AND AVERAGE FRICTION COEFFICIENT ANALYSIS

Keywords: Self-excited vibrations, torsor measurement, friction in turning, torsor central axis

An experimental device in turning including, in particular, a six-component dynamometer is exploited to measure the complete torque of cutting forces in a case of self-excited vibrations. For the tests, the tool used is type TNMA 16 04 12 carbide not covered (nuance carbide-SUMITIMO ELECTRIC), without chip breeze. The machined material is an alloy of chrome molybdenum type 42 CrMo24. The test-tubes are cylindrical with a diameter of 120 mm and a length of 30 mm. The effort of analysis relates to the moments. In particular, to the tool tip point, when feed rate increases the friction coefficient of swivelling is increasing while that of bearing is decreasing. Innovative remarks on the evolution of the average friction coefficient for various depths of cut and feed rate are presented.

Rezumat:

In lucrare se utilizeaza o instalatie de achizitie de date experimentale in cazul strunjirii, in particular fiind utilizat un dinamometru cu sase componente, pentru masurarea fortelor si momentelor in cazul in care sunt prezente autovibratii in cadrul procesului. Pentru aceste teste scula aschietoare este de tip TNMA 16 04 12 din aliaj fara durificare superficiala (structura -SUMITIMO ELECTRIC), fara spargator de aschii. Materialul prelucrat este un aliaj cu crom si molibden de tip 42 CrMo24. Piesele de proba sunt cilindrice, cu diametrul de 120 mm si lungimea de 30 mm. Efortul de analiza se refera in mod special la momente. Sunt prezentate concluzii inovatoare relative la coeficientul mediu de frecare pentru diferite cazuri de variatie pentru adancimea de aschiere si avansul utilizat. Rezultatele prezentate corespund unui numar de aproximativ 7590 masuratori.

Introduction

In the three-dimensional cutting case, the mechanical actions torsor (forces and moments) is often truncated: because the torsor moment part is probably a neglected fault of access to an adapted metrology [START_REF] Mehdi | Dynamic behavior of thin wall cylindrical workpiece during the turning process, Part 2: Experimental approach and validation[END_REF][START_REF] Yaldiz | Design, development and testing of a turning dynamometer for cutting force measurement[END_REF][START_REF] Yaldiz | A dynamometer design for measurement the cutting forces on turning[END_REF]. Unfortunately, until now, the results on the cutting forces are almost still validated using platforms of forces (dynamometers) measuring those three components [START_REF] Marui | Chatter vibration of the lathe tools. Part 2 : on the mechanism of exciting energy supply[END_REF][START_REF] Lian | Self-organizing fuzzy control of constant cutting force in turning[END_REF]. However, forces and pure moments (or torque) can be measured [START_REF] Couétard | Capteurs de forces à deux voies et application à la mesure d'un torseur de force[END_REF]. Recently, an application consisting in six component measurements of the actions torsor in cutting process was carried out for the case of high speed milling [START_REF] Couétard | Mesure des 6 actions de coupe en fraisage grande vitesse[END_REF], drilling [START_REF] Laporte | Phenomenological model for drilling operation[END_REF][START_REF] Yaldiz | Design, development and testing of a fourcomponent milling dynamometer for the measurement of cutting force and torque[END_REF], etc. Cahuc et al., [START_REF] Cahuc | Experimental and analytical balance sheet in turning applications[END_REF], present another use of this six-component dynamometer in an experimental study: the taking into account of the cut moments allows a better machine tool power consumption evaluation. This led to a better cut approach [START_REF] Laporte | Phenomenological model for drilling operation[END_REF][START_REF] Couétard | Caractérisation et étalonnage des dynamomètres à six composantes pour torseur associé à un système de forces[END_REF][START_REF] Bisu | Nouvelle analyse des phénomènes vibratoires en tournage[END_REF] and should enable us to reach new properties of the vibrations of the system piece-tool-matter.

Moreover, the tool torsor has the advantage of being transportable in any space point and, in particular, at the tool tip in O point. This study follows the cut torsor and is carried out in several stages, including two major one. The first is related to the forces and the moments friction coefficient analysis at the tool tip during the cut. The second is dedicated to the torsor central axis determination during the cut. The central axes beams deduced from the multiple tests strongly confirm the moments presence to the tool tip. In a precise way, to section 2 we present primarily the experimental device used. Section 3 is devoted to measure the cut actions torsor reduced to the center of the six-component dynamometer and then transported to the tool tip point. An analysis of the forces exerted during the cut actions is carried out. It allows to establish in experiments certain properties of the resultant of the cut actions. The evolution of the average friction coefficient associate to a depth of cut for several feed rate is also obtained. The case of the moments to the tool tip point is examined. The torsor central axis is required (Section 4) and the central axes beams deduced from multiple measurements confirm the presence of moments to the tool tip point. In Section 5 we more particularly carry out the analysis of the moments to the central axis by looking at the case most sensitive to the vibrations (ap = 5 mm and f = 0.1 mm/rev). This study gives a certain number of properties. Before concluding, some innovative reflections are obtained, in the case of turning presented here.

Experimental device

The experimental device presented Figure 1 is a conventional lathe (Ernault HN400) for which the spindle speed does not exceed 3,500 rpm. The machining system behavior is analyzed through one three direction accelerometer fixed on the tool and using two unidirectional accelerometers positioned on the lathe, on the front bearing of the spindle, to identify the influence of this one on the cutting process. Moreover, a six-component dynamometer [START_REF] Couétard | Capteurs de forces à deux voies et application à la mesure d'un torseur de force[END_REF], being used as toolholder [START_REF] Bisu | New method to characterize a machining system: application in turning[END_REF], is positioned on the lathe to measure all the cutting forces (forces and torques). On the test machining system, the instantaneous speed of rotation is controlled by a rotary encoder directly related to the workpiece. The connection is carried out by a rigid steel wire, which allows a better transmission of the behavior (Figure 1). During the cutting process, the number of revolutions is controlled constantly at nearly 690 rpm, and a negligible variation cutting speed of about 1% is detected. The cylindrical test-tubes have a diameter of 120 mm and a length of 30 mm (Figure 2). The dimensions retained for these test tubes were selected using the finite elements method coupled to an optimization method using SAMCEF software presented in [START_REF] Bisu | Etude des vibrations auto-entretenues en coupe tridimensionnelle : nouvelle modélisation appliquée au tournage[END_REF]. This process describes in [START_REF] Zapciu | Experimental study of machining system: dynamic characterization[END_REF] allows to ensure that the natural frequencies band of the machining system is well outer from that of the machine natural frequencies. Thus under a load P = 1,000 N, the material having a Young modulus E = 21.10 5 N/mm² the workpiece dimensions selected are: D 1 diameter = 60 mm, L 1 length = 180 mm for a bending stiffness of 7.10 7 N/m (Figure 2). This value is including in the higher part of the interval of the acceptable rigidity values for conventional lathe [START_REF] Ispas | Vibrations des systèmes technologiques[END_REF][START_REF] Koenigsberger | Machine Tools Structures[END_REF][START_REF] Konig | Zweikomponenten schnittkraftmesser[END_REF]. Geometry of holding fixture / workpiece During the tests, the used tool was a noncoated carbide tool (TNMA 160412) without chip break geometry. The cutting material is a chrome molybdenum alloy type 42CrMo24. Moreover, the geometry of the tool (Figure 3) [START_REF] Laheurte | Application de la théorie du second gradient à la coupe des matériaux[END_REF] is characterized by the cutting angle γ, the clearance angle α, the edge angle of inclination λ s , the direct angle κ r , the nozzle radius r  and the sharpness radius R. The tool insert is examined after each test and is changed if necessary of wear along the cutting face (Vb ≤ 0.2 mm ISO 3685), that may disturb the studied phenomenon. The tool characteristics used are presented in the Table 1.

    r r  R -6° 6° -6° 91°
1,2 mm 0,02 mm Figure 3 The tool geometrical details Table 1 The tool geometrical characteristics

Cutting action torsor

Tests

The experiments are performed within a framework similar to that exposed in Cahuc and al., [START_REF] Cahuc | Experimental and analytical balance sheet in turning applications[END_REF]. For each test, the mechanical actions are measured according to the feed rate using the six-component dynamometer [START_REF] Couétard | Mesure des 6 actions de coupe en fraisage grande vitesse[END_REF] following the method initiated in Toulouse [START_REF] Toulouse | Contribution à la modélisation et à la métrologie de la coupe dans le cas d'un usinage tridimensionnel[END_REF], developed and finalized by Couétard [START_REF] Couétard | Caractérisation et étalonnage des dynamomètres à six composantes pour torseur associé à un système de forces[END_REF]. These mechanical actions are evaluated for four depths of cut ap (= 1 mm; 2 mm; 3.5 mm; 5 mm) and according to four feed rate f (= 0.1 mm/rev; 0.075 mm/rev; 0.0625 mm/rev; 0.05 mm/rev).

Measurements are taken in the six-component dynamometer transducer O' center, and then transported to the tool point O via the moment transport traditional relations. Uncertainties of measurement with the six-component dynamometer used are about ± 4% for the forces components and ± 6% for the moments components. To check the repeatability and accuracy of identifications all tests and measurements are carried out four times and the average is selected for each configuration specified above. For each series of measure, we analyzed the results over 15 seconds of recording. We have 44 points of measurements recorded per full piece rotation. Thus, the average values discussed thereafter correspond finally to the average values of 7,590 measures.

Resultant of the cutting forces analysis

Among the four values feed rate f indicated above, two examples of the cutting forces resultant measurements applied to the tool tip point are presented in the stable case (quasi non-existent vibrations) ap = 2 mm (Figure 4), and in the unstable case (self-excited vibrations) ap = 5 mm (Figure 5). In the stable case (ap = 2 mm, f = 0.1 mm/rev), we see that the amplitude of all forces components remain almost constant in the course of time. Only is detected a weak amplitudes variation about 1 to 2 N around their nominal values. Let us take for reference the absolute value of the component (F x ) of the forces. In the case of turning here studied, we observe that, whatever the depth of cut ap, there exists the following order relation between the modulus of the average cutting forces components:

|F x |  |F z |  |F y |. ( 1
)
The six-component dynamometer used gives the instantaneous values of all torque components of cutting action in the reference system of machine tool ( x  , y  , z  ) (Figure 1). The average values of the forces components (F x , F y , F z ) in this reference system are presented to Table 2 in the instability case (ap= 5 mm, feed rate variable f). To determine the average values of friction coefficients of the tool/piece it is initially necessary to have the forces components (F x , F ycop , F zcop ) in the tool reference system ( x  , cop y  , cop z  ). This is acquired by making the change of reference system [START_REF] Bakkal | Chip formation, cutting forces, and tool wear in turning of Zr-based bulk metallic glass[END_REF] : In Eq. 2 the F ycop component is carried by the normal with the tool while the components F x and F zcop are in the plan of this one. Thus, it is proceeded to the determination of the average values tangential component modulus of the cutting forces in the tool plan ||T cop ||. This one is given classically by the relation:

F ycop = F y cos  -F z sin ; F zcop = F y sin  + F z cos   = 6° (Table 1) (2)
||T cop || = 2 2 zcop x F F  . ( 3
)
The average normal component modulus is: ||N|| = ||F ycop| |. These elements allow to draw Table 2 the forces components average values applied to the tool in the reference system of machine tool, then in the reference system of the tool for a fixed depth of cut (ap = 5 mm) and a variable feed rate (f = 0.1 mm/rev; 0.075 mm/rev; 0.0625 mm/rev).

f (mm/tr) F x (N) F y (N) F z (N) F ycop (N) F zcop (N) ||T cop || (N) 0,1 200 -1512 1073 -1616 909 946
Table 2. Forces components average values in the reference system of machine tool and tool for the unstable case ap = 5mm, f = 0.1 mm/rev and N = 690 rpm These results lead to the average friction coefficient of Coulomb type represented in Figure 6. We see that the average values coefficient friction is a decreasing function of the feed rate. These results are coherent with those of the literature. However, in our case of three-dimensional cut we cannot neglect the F x component as in [START_REF] Seker | The effect of feed rate on the cutting forces when machining with linear motion[END_REF] for example. As for the torsor resultant, we give two examples of the moment components measurement to the tool tip point in the reference system of the machine tool. As comparison we present initially the stable case (without vibrations ap = 2 mm) Figure 7, then the unstable case (with vibrations ap = 5 mm) Figure 8. As for the resultant in the stable case (ap = 2 mm, f = 0.1 mm/rev), the average moment components to the tool tip point are slightly disturbed while in the unstable case (ap = 5 mm, f = 0.1 mm/rev) the average values of the moment components of the tool tip point are definitely more chaotic.

However, one can note that the average moment components to the tool tip point are more disturbed than their equivalents on the level of the resultant. Thus, the moment components seem more sensitive to the self-sustained vibrations than the resultant components of the forces applied. The follow-up of the evolution of these could, for example, being a good means of detecting the existence of regenerative vibrations precociously.

As for the forces, on the moment components level we have, always in the case of turning examined here, the following order relation between the average values modulus of the moment components to the tool tip point (on 7,590 points of measurements).

|M oy |  |M oz |  |M ox | (4)
According to the confrontation of the Eq. ( 1) and the Eq. ( 4), it thus proves that on the order relations level between the modulus of the forces and the moment components average values to the tool tip point, the role of the x  and y  axes is reversed. This can be allotted to transport the moments from the tool tip point to the central axis.

In order to evaluate the average values friction coefficients of swivelling and bearing, we must now place ourselves in the reference system of the tool. Relations similar to Eq. ( 2) are then applied to the moment components transported to the tool tip point O. For the unstable case (ap =5 mm, f = 0.1 mm/rev) the moment average values are consigned to Table 3 Tableau 3. Moment components average values of the forces in the reference system of the machine tool and this one of the tool for the unstable case ap = 5mm, f = 0.1 mm/rev and N = 690 rpm.

From this table of values, we deduce the average friction coefficients from bearing (M o , Roul ) and from swivelling (M o , piv ) to the tool tip point starting from the following traditional definitions:

||M otcop || = 2 2 ozcop ox M M  , m o,Roul = N M otcop ; m o,piv = N M oycop ( 5
)
Figure 9 shows the evolution of these average friction coefficients according to feed rate. It is noted that the average friction coefficient of swivelling to the tool tip point is nearly constant even very slightly growing with f. On the other hand the average friction coefficient of bearing has a rather different behavior; this last is, all things considered, strongly decreasing when feed rate increases.

It is also noted that the average friction coefficients of bearing and swivelling are confused for the strong values feed rate (f = 0.075 mm/rev and f = 0.1 mm/rev). On the other hand for the low values feed rate (f= 0.05 mm/rev and f= 0.0625 mm/rev) these points are well dissociated. A thorough study relating to one more a large number of feed rate is planned in order to consolidate these results.

Central axis determination

It is well-known that, with any torsor, it is possible to associate a central axis (except the torsor of pure moment), which is the single object calculated starting from the six components torsor [START_REF] Brousse | Cours de mécanique (1 er cycle et classe prépa[END_REF]. A torsor [A] O in a point O is composed of a resultant forces R  and the resulting moment o M  :

[A] o =    o M R (6)
The central axis is the line defined classically by:

OA = 2 R M R O  +  R . (7)
In the Eq. ( 7) O is the point where the mechanical actions torsor was moved (here, the tool tip point ) and A is the current point describing the central axis.

Thus, OA is the vector associated with the bipoint [O, A] (Figure 10). This line (Figure 10a) corresponds to geometric points where the mechanical actions moment torsor is minimal. The central axis calculation consists in determining the points assembly (a line) where the torsor can be expressed according to a slide block (straight line direction) and the pur moment (or torque). The central axis is also the point where the resultant cutting force is colinear with the minimum mechanical moment (pure torque). The test results enable us to check for each point of measurement where the colinearity between the resultant cutting force R  and moment A M  calculated is related to the central axis (Figure 10b). The meticulous examination of the six mechanical action torsor components shows that the forces and the moments average values are not null. For each measure point, the central axis is calculated, in the stable (Figure 11a) and unstable modes (Figure 11b). In any rigor, the case ap = 2 mm and f = 0.1 mm/rev should be described as "quasistable" movement because the vibrations exist but, as we already noticed, their amplitudes are very low, about the micrometer. Thus, they are quasi null compared to the other studied cases. Indeed, for example for ap = 5 mm and f = 0.0625 mm/rev, the recorded amplitude was 10 times more important. In the presence of vibrations (ap=5 mm), it is possible to observe for 68 revolutions of the machined piece the dispersive character of the central axes beam (fig. 11b), compared to the stable mode, where this same beam (fig. 11a) is tightened more or less tilted compared to the normal axis on the plan (x,y). This central axes dispersion can be explained by the self-sustained vibrations which cause the variable forces and moments generation.

Moments friction coefficient analysis at central axis

While transporting the moment from tool tip point to the central axis, the minimum moment (or pure torque) A M  is obtained. From the moment values to the central axis it is possible to deduce the constant and variable part of this one, as for the forces [START_REF] Bisu | Self-excited vibrations in turning: Forces torsor analysis[END_REF]. The analysis shows that the vibrations generate rotations of the tool, causes variations of contact and thus generates variable moments. This representation allows to express the moments along the three axes of the machine tool: swivelling moment in the y  direction and the two bearing moments following the x  and the z  direction, then in the ( x  , cop y  , cop z  ) reference system of the tool ( With the central axis the experimental data which preceding show, always in the case of turning studied here, the following order relation:

||M Ax ||  ||M Az ||  ||M Ay ||. (8)
Thus, according to Eq. 1 the modulus of the average moments components to the central axis are in the same order as those of the average values force components. This is not the case to the tool tip point where the role of the axes x  and y  is reversed according to Eq. 4. Thus, the effect of the moment transport is sensitive. Lastly, by using the central axis the definitions given in Eq. 5, we obtain the swivelling and the bearing average values coefficients of the pure torques. The evolution of those according to the feed rate is represented Figure 12. It is noted that the swivelling friction coefficient is constant to the central axis, while that of the bearing is slightly increasing when feed rate increases. This last tendency is opposite of that which we found to the tool tip, and in less accentuated. For the swivelling friction coefficient if it is constant on the central axis (Figure 12), on the other hand, to the tool tip (Figure 9) it slightly believes with feed rate (in the case of turning studied here). Thus one notes that with the central axis, the average values friction coefficients are distinct (with a common average value in swivelling and bearing for f = 0.075 min/rev). This is obviously a tangible consequence of more than influence of the moments transport and need for analyzing the pure torques in machining what is the principal originality of this work.

Conclusion

The experimental procedures installation allow to determine the elements necessary to a rigorous analysis of the average friction coefficients from the mechanical actions torsor. For the first time, it seems to us, is presented within the framework of the self-excited vibrations in turning a brief study of the resulting moments of the mechanical action torsor to the tool tip point and to the central axis. If the force components are compared with those of the moment, it is possible to note, that in absolute value, the average values are classified in the same order if the average values of the moment components are taken to the central axis. On the other hand, the roles of these same average values on the axes x  and y  are reversed if we look at the moments to the tool tip. The z  axis however keeps a role of "pivot" (in the case of turning examined here).

Bearing and swivelling average values friction coefficients analysis allow to establish a rather clear difference between what occurs to tool tip and to central axis. In particular to the central axis, the average values friction coefficient of swivelling is constant what is not completely the case to the tool tip where a light growth is observed according to feed rate selected here. In the same way a notable difference exists between the average values friction coefficient of bearing to the central axis and to the tool tip; the evolutions are reversed. The decrease of the average friction coefficient of bearing is more important to the tool tip than the growth of this one to the central axis when the feed rate increases.

Figure 1

 1 Figure 1 Experimental machining system and metrology environment in turning process

Figure 2

 2 Figure 2Geometry of holding fixture / workpiece

Figure 4 Figure 5

 45 Figure 4 Signals of the forces components exerted to the tool tip point according to the three directions of space machine in the stable case ap = 2 mm, f = 0.1 mm/rev, N = 690 rpm Figure 5 Signals of the forces components exerted to the tool tip point according to the three directions of space machine in the unstable case ap = 5 mm, f = 0.1 mm/rev, N = 690 rpm

Figure 6 .

 6 Figure 6. Average values friction coefficient evolution (||T cop || / ||N||) according to the feed rate for the unstable case (ap = 5 mm, f = 0.1 mm/rev and N = 690 rpm)

Figure 7 .Figure 8 .

 78 Figure 7. Signals of the moment components exerted to the tool tip point according to the three directions of the machine tool in the stable case ap = 2 mm, f = 0.1 mm/rev, N = 690 rpm Figure 8. Signals of the moment components exerted to the tool tip point according to the three directions of the machine tool in the unstable case ap = 5 mm, f == 0,1 mm/rev, N = 690 rpm

Figure 9 .

 9 Figure 9. Swivelling and bearing friction coefficient evolution to the tool tip point.

Figure 10 .

 10 Figure 10. Central axis representation (a) and of the colinearity between vector sum R  and

Figure 11

 11 Figure 11 Central axes representation obtained for 68 rpm the workpiece speed and feed rate f = 0.0625 mm/rev (a) stable process ap = 2 mm; (b) unstable process ap = 5 mm

Figure 12 .

 12 Figure 12. Friction coefficient evolution of swivelling and bearing moments to the central axis.

  where ||M otcop || indicate the module of the average moment in the tool plan (cf. Eq. 3).

	f (mm/tr)	M ox (dNm)	M oy (dNm)	M oz (dNm)	M ycop (dNm)	M zcop (dNm)	||M otcop || (dNm)
	0,1	168	-101	-135	-86,4	-144,8	225,8

Table 4

 4

		).
	f (mm/tr) M Ax (dNm)	M Ay (dNm)	M Az (dNm)	M Aycop (dNm)	M Azcop (dNm)	||M Atcop || (dNm)
	0,1	-4,1	-33,5	8,6	-34,2	5	41,7

Table 4 .

 4 Average values of the moment components to central axis in the machine tool reference

	system (x, y, z) and in the tool reference system ( x  , cop y  , cop z  ) for the unstable case (ap = 5 mm,
	f = 0.1 mm/rev and N = 690 rpm)

Prof., Dr, Université Bordeaux 1 (Sciences et Technologies) et CNRS Laboratoire de Mécanique Physique UMR 5496, Chair de Mécanique, 351 cours de la Libération, 33405 Talence, France, (email: alain.gerard@u-bordeaux1.fr).

[START_REF] Yaldiz | Design, development and testing of a turning dynamometer for cutting force measurement[END_REF] Prof., Dr, Eng., Universitate Politehnica from

Bucharest, 313 splaiul Independentei, MSP, 60042 Bucharest, Romania.


Nicolas Verzelen

Data-driven neighborhood selection of a Gaussian field

Keywords: Gaussian field, Gaussian Markov random field, Data-driven calibration, model selection, pseudolikelihood

We study the nonparametric covariance estimation of a stationary Gaussian field X observed on a lattice. To tackle this issue, a neighborhood selection procedure has been recently introduced. This procedure amounts to selecting a neighborhood m by a penalization method and estimating the covariance of X in the space of Gaussian Markov random fields (GMRFs) with neighborhood m. Such a strategy is shown to satisfy oracle inequalities as well as minimax adaptive properties. However, it suffers several drawbacks which make the method difficult to apply in practice. On the one hand, the penalty depends on some unknown quantities. On the other hand, the procedure is only defined for toroidal lattices. The present contribution is threefold. A data-driven algorithm is proposed for tuning the penalty function. Moreover, the procedure is extended to non-toroidal lattices. Finally, numerical study illustrate the performances of the method on simulated examples. These simulations suggest that Gaussian Markov random field selection is often a good alternative to variogram estimation.

Sélection automatique de voisinage d'un champ gaussien

Résumé : Nous étudions l'estimation non-paramétrique d'un champ gaussien stationnaire X observé sur un réseau régulier. Dans ce cadre, nous avons précédemment introduit une procédure de sélection de modèle [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. Cette procédure revient à sélectionner un voisinage m grâce une technique de pénalisation puis à estimer la covariance du champ X dans l'espace des champs de Markov gaussiens de voisinage m. Une telle stratégie satisfait des inégalités oracles et des propriétés d'apdaptation au sens minimax. En pratique, elle présente néanmoins quelques inconvénients. D'une part, la pénalité dépend de quantités inconnues. D'autre part, la procédure est uniquement définie pour des réseaux toriques. La contribution de cet article est triple. Nous proposons un algorithme automatique pour calibrer la pénalité. De plus, nous introduisons une extension à des réseaux non-toriques. Enfin, nous étudions les performances pratiques de la procédure sur des données simulées. Ces simulations suggèrent que la sélection de champs de Markov gaussiens est souvent une bonne alternative à l'estimation de variogramme.

Mots-clés : Champ gaussien, champ de Markov gaussien, calibration automatique, sélection de modèle, pseudo-vraisemblance.

Introduction

We study the estimation of the distribution of a stationary Gaussian field (X[i,j]) (i,j)∈Λ indexed by the nodes of a rectangular lattice Λ of size p 1 ×p 2 . This problem is often encountered in spatial statistics or in image analysis. Classical statistical procedures allow to estimate and subtract the trend. Henceforth, we assume that the field X is centered. Given a n-sample of the field X, the challenge is to infer the correlation. In practice, the number n of observations often equals one. Different methods have been proposed to tackle this problem.

A traditional approach amounts to computing an empirical variogram and then fitting a suitable parametric variogram model such as the exponential or Matérn model (see [START_REF] Cressie | Statistics for spatial data[END_REF] Ch.2 or [START_REF] Stein | Interpolation of spatial data[END_REF]). The main disadvantage with this method is that the practitioner is required to select a good variogram model. When the field exhibits long range dependence, specific procedures have been introduced (e.g. Frías et al. [START_REF] Frías | Semiparametric estimation of spatial long-range dependence[END_REF]). In the sequel, we focus on small range dependences. Most of the nonparametric (Hall et al. [START_REF] Hall | On the nonparametric estimation of covariance functions[END_REF]) and semiparametric (Im et al. [START_REF] Im | Semiparametric estimation of spectral density with irregular observations[END_REF]) methods are based on the spectral representation of the field. To our knowledge, these procedures have not yet been shown to achieve adaptiveness, i.e. their rate of convergence does not adapt to the complexity of the correlation functions.

In this paper, we define and study a nonparametric estimation procedure relying on Gaussian Markov random fields (GMRF). This procedure is computationally fast and satisfies adaptive properties. Let us fix a node (0, 0) at the center of Λ and let m be a subset of Λ \ {(0, 0)}. The field X is a GMRF with respect to the neighborhood m if conditionally to (X[k,l]) (k,l)∈m , the variable X[0,0] is independent from all the remaining variables in Λ. We refer to Rue and Held [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF] for a comprehensive introduction on GMRFs. If we know that X is a GMRF with respect to the neighborhood m, then we can estimate the covariance by applying likelihood or pseudolikelihood maximization. Such parametric procedures are well understood, at least from an asymptotic point of view (see for instance [START_REF] Guyon | Random fields on a network[END_REF] Sect.4). However, we do not know in practice what is the "good" neighborhood m. For instance, choosing the empty neighborhood amounts to assuming that all the components of X are independent. Alternatively, if we choose the complete neighborhood, which contains all the nodes of Λ except (0, 0), then the number of parameters is huge and estimation performances are poor.

We tackle in this paper the problem of neighborhood selection from a practical point of view. The purpose is to define a data-driven procedure that picks a suitable neighborhood m and then estimates the distribution of X in the space of GMRFs with neighborhood m. This procedure neither requires any knowledge on the correlation of X, nor assumes that the field X satisfies a Markov condition. Indeed, the procedure selects a neighborhood m that achieves a trade-off between an approximation error (distance between the true correlation and GMRFs with neighborhood m) and an estimation error (variance of the estimator). If X is a GMRF with respect to a small neighborhood, then the procedure achieves a parametric rate of convergence. Alternatively, if X is not a GMRF then the rate of convergence of the procedure depends on the rate of approximation of the true covariance by GMRFs with growing neighborhood. In short, the procedure is nonparametric and adaptive.

Besag and Kooperberg [START_REF] Besag | On conditional and intrinsic autoregressions[END_REF], Rue and Tjelmeland [START_REF] Rue | Fitting Gaussian Markov random fields to Gaussian fields[END_REF], Song et al. [START_REF] Song | A comparative study of gaussian geostatistical models and gaussian markov random field models[END_REF], and Cressie and Verzelen [START_REF] Cressie | Conditional-mean least-squares of Gaussian Markov random fields to Gaussian fields[END_REF] have considered the problem of approximating the correlation of a Gaussian field by a GMRF, but this approach requires the knowledge of the true distribution. Guyon and Yao have stated in [START_REF] Guyon | On the underfitting and overfitting sets of models chosen by order selection criteria[END_REF] necessary conditions and sufficient conditions for a model selection procedure to choose asymptotically the true neighborhood of a GMRF with probability one. Our point of view is slightly different. We do not assume that the field X is a GMRF with respect to a sparse neighborhood. We do not aim at estimating the true neighborhood, we rather want to select a neighborhood that allows to estimate well the distribution of X (i.e. to minimize a risk). The distinction between these two points of view has been nicely described in the first chapter of MacQuarrie and Tsai [START_REF] Mcquarrie | Regression and time series model selection[END_REF].

In [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF], we have introduced a neighborhood selection procedure based on pseudolikelihood maximization and penalization. Under mild assumptions, the procedure achieves optimal neighborhood selection. More precisely, it satisfies an oracle inequality and it is minimax adaptive to the sparsity of the neighborhood. To our knowledge, these are the first results of neighborhood selection in this spatial setting.

If the procedure exhibits appealing theoretical properties, it suffers several drawbacks from a practical perspective. First, the method constrains the largest eigenvalue of the estimated covariance to be smaller than some parameter ρ. In practice, it is difficult to choose ρ since we do not know the largest eigenvalue of the true covariance. Second, the penalty function pen(.) introduced in Sect.3 of the previous paper depends on the largest eigenvalue of the covariance of the field X. Hence, we need a practical method for tuning the penalty. Third, the procedure has only been defined when the lattice Λ is a square torus.

Our contribution is twofold. On the one hand, we propose practical versions of our neighborhood selection procedure that overcome the previously-mentioned drawbacks:

The procedure is extended to rectangular lattices. We do not constrain anymore the largest eigenvalue of the covariance. We provide an algorithm based on the so-called slope heuristics of Birgé and Massart [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] for tuning the penalty. Theoretical justifications for its use are also given. Finally, we extend the procedure to the case where the lattice Λ is not a torus.

On the other hand, we illustrate the performances of this new procedure on numerical examples. When Λ is a torus, we compare it with likelihood-based methods like AIC [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] and BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF], even if they were not studied in this setting. When Λ is not toroidal, likelihood methods become intractable. Nevertheless, our procedure still applies and often outperforms variogrambased methods.

The paper is organized as follows. In Section 2, we define a new version of the estimation procedure of [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] that does not require anymore the choice of the constant ρ. We also discuss the computational complexity of the procedure. In Section 3, we connect this new procedure to the original method and we recall some theoretical results. We provide an algorithm for tuning the penalty in practice in Section 4. In Section 5, we extend our procedure for handling nontoroidal lattices. The simulation studies are provided in Section 6. Section 7 summarizes our findings, while the proofs are postponed to Section 8.

Let us introduce some notations. In the sequel, X v refers to the vectorialized version of X with the convention X[i,j] = X v [(i-1)×p2+j] for any 1 ≤ i ≤ p 1 and 1 ≤ j ≤ p 2 . Using this new notation amounts to "forgetting" the spatial structure of X and allows to get into a more classical statistical framework. We note X 1 , X 2 , . . . , X n the n observations of the field X. The matrix Σ stands for the covariance matrix of X v . For any matrix A, ϕ max (A) and ϕ min (A) respectively refer the largest eigenvalue and the smallest eigenvalues of A. Finally, I r denotes the identity matrix of size r.

INRIA 2 Neighborhood selection on a torus

In this section, we introduce the main concepts and notations for GMRFs on a torus. Afterwards, we describe our procedure based on pseudolikelihood maximization. Finally, we discuss some computational aspects. Throughout this section and the two following sections, the lattice Λ is assumed to be toroidal. Consequently, the components of the matrices X are taken modulo p 1 and p 2 .

GMRFs on the torus

The notion of conditional distribution is underlying the definition of GMRFs. By standard Gaussian derivations (see for instance [START_REF] Lauritzen | Graphical models, volume 17 of Oxford Statistical Science Series[END_REF] App.C), there exists a unique p 1 × p 2 matrix θ such that θ[0,0] = 0 and

X[0,0] = (i,j)∈Λ\{(0,0)} θ[i,j]X[i,j] + ǫ[0,0] , (1)
where the random variable ǫ[0,0] follows a zero-mean normal distribution and is independent from the covariates (X[i,j]) (i,j)∈Λ\{(0,0)} . The linear combination (i,j)∈Λ\{(0,0)} θ[i,j]X[i,j] is the kriging predictor of X[0,0] given the remaining variables. In the sequel, we note σ 2 the variance of ǫ[0,0] and we call it the conditional variance of X[0,0]. Equation (1) describes the conditional distribution of X[0,0] given the remaining variables. By stationarity of the field X, it holds that that θ[i,j] = θ[-i,-j]. The covariance matrix Σ is closely related to θ through the following equation:

Σ = σ 2 [I p1p2 -C(θ)] -1 , (2)
where the

p 1 p 2 × p 1 p 2 matrix C(θ) is defined by C(θ)[(i1-1)p2+j1,(i2-1)p2+j2] := θ[i2-i1,j2-j1]
for any 1 ≤ i 1 , i 2 ≤ p 1 and 1 ≤ j 1 , j 2 ≤ p 2 . The matrix (I p1p2 -C(θ)) is called the partial correlation matrix of the field X. The so-defined matrix C(θ) is symmetric block circulant with p 2 × p 2 blocks. We refer to [RH05] Sect.2.6 or the book of Gray [START_REF] Gray | Toeplitz and Circulant Matrices: A Review[END_REF] for definitions and main properties on circulant and block circulant matrices.

Identities (1) and (2) have two main consequences. First, estimating the p 1 × p 2 matrix θ amounts to estimating the covariance matrix Σ up to a multiplicative constant. We shall therefore focus on θ. Second, by Equation (1), the field X is a GMRF with respect to the neighborhood defined by the support θ. The adaptive estimation issue of the distribution of X by neighborhood selection therefore reformulates as an adaptive estimation problem of the matrix θ via support selection.

Let us now precise the set of possible values for θ. The set Θ denotes the vector space of the p 1 × p 2 matrices that satisfy θ[0,0] = 0 and θ[i,j] = θ[-i,-j], for any (i, j) ∈ Λ. Hence, a matrix θ ∈ Θ corresponds to the distribution of a stationary Gaussian field if and only if the p 1 p 2 × p 1 p 2 matrix (I p1p2 -C(θ)) is positive definite. This is why we define the convex subset Θ + of Θ by

Θ + := {θ ∈ Θ s.t. [I p1p2 -C(θ)] is positive definite} . ( 3
)
The set of covariance matrices of stationary Gaussian fields on Λ with unit conditional variance is in one to one correspondence with the set Θ + . We sometimes assume that the field X is isotropic. The corresponding sets Θ iso and Θ +,iso for isotropic fields are introduced as:

Θ iso := {θ ∈ Θ , θ[i,j] = θ[-i,j] = θ[j,i] ,
∀(i, j) ∈ Λ} and Θ +,iso := Θ + ∩ Θ iso .

RR n 6798

Description of the procedure

Let |(i, j)| t refer to the toroidal norm defined by

|(i, j)| 2 t := [i ∧ (p 1 -i)] 2 + [j ∧ (p 2 -j)] 2 ,
for any node (i, j) ∈ Λ.

In the sequel, a model m stands for a subset of Λ \ {(0, 0)}. It is also called a neighborhood. For the sake of simplicity, we shall only use the collection of models M 1 defined below.

Definition 2.1. A subset m ⊂ Λ \ {(0, 0)} belongs to M 1 if and only if there exists a number

r m > 1 such that m = {(i, j) ∈ Λ \ {(0, 0)} s.t. |(i, j)| t ≤ r m } . (4)
In other words, the neighborhoods m in M 1 are sets of nodes lying in a disc centered at (0, 0). Obviously, M 1 is totally ordered with respect to the inclusion. Consequently, we order the models m 0 ⊂ m 1 ⊂ . . . ⊂ m i . . .. For instance, m 0 corresponds to the empty neighborhood, m 1 stands for the neighborhood of size 4, and m 2 refers to the neighborhood with 8 neighbours. See Figure 1 for an illustration. For any model m ∈ M 1 , the vector space Θ m is the subset of matrices Θ whose support is included in m. Similarly Θ iso m is the subset of Θ iso whose support is included in m. The dimensions of Θ m and Θ iso m are respectively noted d m and d iso m . Since we aim at estimating the positive matrix (I p1p2 -C(θ)), we also consider the convex subsets of Θ + m and Θ +,iso m which correspond to non-negative precision matrices.

Θ + m := Θ m ∩ Θ + and Θ +,iso m := Θ iso m ∩ Θ +,iso . (5)
For any θ ′ ∈ Θ + , the conditional least-squares (CLS) criterion γ n,p1,p2 (θ ′ ) [START_REF] Guyon | Estimation d'un champ par pseudo-vraisemblance conditionnelle: étude asymptotique et application au cas markovien[END_REF] is defined by

γ n,p1,p2 (θ ′ ) := 1 np 1 p 2 n i=1 (j1,j2)∈Λ X i [j1,j2] - (l1,l2)∈Λ\{(0,0)} θ ′ [l1,l2]X i [j1+l1,j2+l2] 2 . (6)

INRIA

The function γ n,p1,p2 (.) is a least-squares criterion that allows us to perform the simultaneous linear regression of all X i [j1,j2] with respect to the covariates (X i [l1,l2]) (l1,l2) =(k1,k2) . This criterion is closely connected with the pseudolikelihood introduced by Besag [START_REF] Besag | Statistical Analysis of Non-Lattice Data[END_REF]. The associated estimator is slightly less efficient estimator than maximum likelihood estimation ([Guy95] Sect.4.3). Nevertheless, its computation is much faster since it does not involve determinants as for the likelihood. See [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] Sect. 7.1, for a more complete comparison between CLS and maximum likelihood estimators in this setting. For any model m ∈ M 1 , the estimators are defined as the unique minimizers of γ n,p1,p2 (.) on the sets Θ + m and Θ +,iso

Since the lattice Λ is a torus, the computation of the estimators θ m is performed efficiently thanks to the following lemma.

Lemma 2.1. For any p × p matrix A and for any 1 ≤ i ≤ p 1 and 1 ≤ j ≤ p 2 , let λ[i,j](A) be the (i, j)-th term of two-dimensional discrete Fourier transform of the matrix A, i.e.

λ[i,j](A) := p1 k=1 p2 l=1 A[i,j] exp 2ιπ ki p 1 + jl p 2 , ( 9
)
where ι 2 = -1. The conditional least-squares criterion γ n,p1,p2 (θ ′ ) simplifies as

γ n,p1,p2 (θ ′ ) = 1 np 2 1 p 2 2 p1 i=1 p2 j=1 [1 -λ[i,j](θ)] 2 n k=1 λ[i,j] (X k ) λ[i,j] (X k ) .
A proof is given in Section 8. Optimization of γ n,p1,p2 (.) over the set Θ + m is performed fastly using the fast Fourier transform (FFT). Nevertheless, this is not the privilege of CLS estimators, since maximum likelihood estimators are also computed fastly by FFT when Λ is a torus.

In Section 5, we mention that the computation of the CLS estimators θ m remains quite easy when Λ is not a torus whereas likelihood maximization becomes intractable.

Theoretical results

Throughout this section, Λ is assumed to be a toroidal square lattice and we note p its size. Let us mention that the restriction to square lattices made in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] allows to simplify the proofs but is not necessary so that the theoretical results hold. In this section, we first recall the original procedure and we emphasize the differences with the one defined in the previous section. We also mention a result of optimality. This will provide some insights for calibrating the penalty pen(.) in Section 4.

Given ρ > 2 be a positive constant, we define the subsets Θ + m,ρ and Θ +,iso m,ρ by

Θ + m,ρ := θ ∈ Θ + m , ϕ max [I p1p2 -C(θ)] < ρ (10) Θ +,iso m,ρ := θ ∈ Θ +,iso m , ϕ max [I p1p2 -C(θ)] < ρ .
Then, the corresponding estimators θ m,ρ and θ iso m,ρ are defined as in ( 7), except that we now consider Θ + m,ρ instead of Θ + m . Let us mention that the estimator θ m corresponds to the estimator

θ m,ρ1 defined in [Ver09] Sect.2.2 with ρ 1 = +∞. θ m,ρ := arg min θ ′ ∈Θ + m,ρ γ n,p,p (θ ′ )
and θ iso m,ρ := arg min

θ ′ ∈Θ +,iso m,ρ γ n,p,p (θ ′ ) .
Given a subcollection M of M 1 and a penalty function pen(.), we select the models m ρ and m iso ρ as in (8) except that we use θ m,ρ and θ iso m,ρ instead of θ m and θ iso m . We also note θ ρ and θ iso ρ for θ b mρ,ρ and θ iso b m iso ρ ,ρ . The only difference between the estimators θ and θ ρ is that the largest eigenvalue of the precision matrix (I p 2 -C( θ)) is restricted to be smaller than ρ. We make this restriction in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] to facilitate the analysis.

In order to assess the performance of the penalized estimator θ ρ and θ iso ρ , we use the prediction loss function l(θ 1 , θ 2 ) defined by

l(θ 1 , θ 2 ) := 1 p 2 tr [(C(θ 1 ) -C(θ 2 ))Σ(C(θ 1 ) -C(θ 2 ))] . (11)
As explained in [Ver09] Sect.1.3, the loss l(θ 1 , θ 2 ) expresses in terms of conditional expectation

l(θ 1 , θ 2 ) = E θ E θ1 X[0,0]|X Λ\{0,0} -E θ2 X[0,0]|X Λ\{0,0} 2 , ( 12
)
where E θ (.) stands for the expectation with respect to the distribution N (0, σ 2 (I p1p2 -C(θ)) -1 ). Hence, l( θ, θ) corresponds the mean squared prediction loss of X[0,0] given the other covariates.

A similar loss function is also used by Song et al. [START_REF] Song | A comparative study of gaussian geostatistical models and gaussian markov random field models[END_REF], when approximation Gaussian fields by GMRFs. For any neighborhood m ∈ M, we define the projection θ m,ρ as the closest element of θ in Θ + m,ρ with respect to the loss l(., .). θ m,ρ := arg min

θ ′ ∈Θ + m,ρ l(θ ′ , θ) and θ iso m,ρ := arg min θ ′ ∈Θ +,iso m,ρ l(θ ′ , θ) .
We call the loss l(θ m,ρ , θ) the bias of the set Θ + m,ρ . This implies that θ m,ρ cannot perform better than this loss.

Theorem 3.1. Let ρ > 2, K be a positive number larger than an universal constant K 0 and M be a subcollection of

M 1 . If for every model m ∈ M, it holds that pen(m) ≥ Kρ 2 ϕ max (Σ) d m + 1 np 2 , (13)
INRIA then for any θ ∈ Θ + , the estimator θ ρ satisfies

E θ [l( θ ρ , θ)] ≤ L(K) inf m∈M [l(θ m,ρ , θ) + pen(m)] , (14)
where L(K) only depends on K. A similar bound holds if one replaces θ ρ by θ iso ρ , Θ + by Θ +,iso , θ m,ρ by θ iso m,ρ , and d m by d iso m .

Although we have assumed the correlation is non-singular, the theorem still holds if the spatial field is constant. The nonasymptotic bound is provided in a slightly different version in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. It states that θ ρ achieves a trade-off between the bias and a variance term if the penalty is suitable chosen. In Theorem 3.1, we use the penalty Kρ 2 ϕ max (Σ)(d m + 1)/(np 2 ) instead of the penalty Kρ 2 ϕ max (Σ)d m /(np 2 ) stated in the previous paper. This makes the bound (14) simpler. Observe that these two penalties yield the same model selection since they only differ by a constant. Let us further discuss two points.

In this paper, we use the estimator θ rather than θ ρ . Given a collection of models M, there exists some finite ρ > 2, such that these two estimators coincide. Take for instance ρ = sup m∈M sup θ∈Θ + m ϕ max (I p1p2 -C(θ)). Admittedly, the so-obtained ρ may be large, especially if there are large models in M. The upper bound (14) on the risk therefore becomes worse. Nevertheless, we do not think that the dependency of (14) on ρ is sharp. Indeed , we illustrate in Section 6 that the risk of θ exhibits good statistical performances. Theorem 3.1 provides a suitable form of the penalty for obtaining oracle inequalities. However, this penalty depends on ϕ max (Σ) which is not known in practice. This is why we develop a data-driven penalization method in the next section.

Slope Heuristics

Let us introduce a data-driven method for calibrating the penalty function pen(.). It is based on the so-called slope heuristic introduced by Birgé and Massart [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] in the fixed design Gaussian regression framework (see also [START_REF] Massart | Concentration inequalities and model selection[END_REF] Sect.8.5.2). This heuristic relies on the notion of minimal penalty. In short, assume that one knows that a good penalty has a form pen(m) = N F (d m ) (where d m is the dimension of the model and N is a tuning parameter). Let us define m(N ) the selected model as a function of N . There exists a quantity N min satisfying the following property: If N > N min , the dimension of the selected model d b m(N ) is reasonable and if N < N min , the dimension of the selected model is huge. The function pen min (.) := N min F (.) is called the minimal penalty. In fact, a dimension jump occurs for d b m(N ) at the point N min . Thus, the quantity N min is clearly observable for real data sets. In their Gaussian framework, Birgé and Massart have shown that twice the minimal penalty is nearly the optimal penalty. In other words, the model m := m(2 N min ) yields an efficient estimator.

The slope heuristic method has been successfully applied for multiple change-point detection [START_REF] Lebarbier | Detecting multiple change-points in the mean of a gaussian process by model selection[END_REF]. Applications are also being developed in other frameworks such as mixture models [START_REF] Maugis | Slope heuristics for variable selection and clustering via gaussian mixtures[END_REF], clustering [START_REF] Baudry | Selecting models focussing the modeller's purpose[END_REF], estimation of oil reserves [START_REF] Lepez | Some estimation problems related to oil reserves[END_REF], and genomic [START_REF] Villers | Tests et sélection de modèles pour l'analyse de données protéomiques et transcriptomiques[END_REF].

If this method was originally introduced for fixed design Gaussian regression, Arlot and Massart [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] have proved more recently that a similar phenomenon occurs in the heteroscedastic random-design case. In the GMRF setting, we are only able to partially justify this heuristic. For the sake of simplicity, let us assume in the next proposition that the lattice Λ is a square of size p.

Verzelen

Proposition 4.1. Consider ρ > 2, and η < 1 and suppose that p is larger than some numerical constant p 0 . Let m ′ be the largest model in M 1 that satisfies d m ′ ≤ np 2 . For any model m ∈ M 1 , we assume that

pen(m ′ ) -pen(m) ≤ K 1 (1 -η)σ 2 ϕ min I p 2 -C(θ) ∧ ρ -ϕ max I p 2 -C(θ) d m ′ -d m np 2 , ( 15
)
where K 1 is a universal (constant defined in the proof ). Then, for any θ ∈ Θ + m ′ ,ρ , it holds that

P d b mρ > L np 2 ∧ p 2 ≥ 1 2 ,
where L only depends on η, ρ, ϕ min I p 2 -C(θ) , and

ϕ max I p 2 -C(θ) .
The proof is postponed to Section 8. Let us define

N 1 := K 1 σ 2 {ϕ min (I p1p2 -C(θ)) ∧ [ρ -ϕ max (I p1p2 -C(θ))]} ,
and let us consider penalty functions pen(m) = N dm np1p2 for some N > 0. The proposition states that if N is smaller than N 1 , then the procedure selects a model of huge dimension with large probability, i.e d b m(N ) is huge. Alternatively, let us define

N 2 := K 0 σ 2 ρ 2 ϕ min (I p1p2 -C(θ)) d m np 1 p 2 ,
where the numerical constant K 0 is introduced in Theorem 3.1 in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. By Theorem 3.1, choosing N > N 2 ensures that the risk of θ ρ achieves a type-oracle inequality and the dimension d b mρ(N ) is reasonable. The quantities N 1 and N 2 are different especially when the eigenvalues of (I p1p2 -C(θ)) are far from 1. Since we do not know the behavior of the selected model m ρ (N ) when N is between N 1 and N 2 , we are not able to really prove a dimension jump as the fixed design Gaussian regression framework. Besides, we have mentioned in the preceding section that we are more interested in the estimator θ than θ ρ . Nevertheless, we clearly observe in simulation studies a dimension jump for some N between N 1 and N 2 even if we use the estimators θ m instead of θ m,ρ . This suggests that the slope heuristic is still valid in the GMRF framework.

m(N ) ∈ arg min m∈M γ n,p1,p2 θ m + N d m np 1 p 2 . 2. Find N min > 0 such that the jump d b m " [ b Nmin] - " -d b m " [ b Nmin] + " is maximal. 3. Select the model m = m(2 N min ).
The difference f (x -)f (x + ) measures the discontinuity of a function f at the point x.

Step 2 may need to introduce huge models in the collection M all the other ones being considered as "reasonably small". As the function m(.) is piecewise linear with at most Card(M) jumps, so that steps 1-2 have a complexity O (Card(M)) 2 . We refer to App.A.1 of [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] for more details on the computational aspects of steps 1 and 2. Let us mention that there are other ways INRIA of estimating N min than choosing the largest jump as described in [AM09] App.A.2. Finally, the methodology described in this section straightforwardly extends to the case of isotropic GMRFs estimation by replacing m(N ) by m iso (N ) and d m by d iso m .

In conclusion, the neighborhood selection procedure described in Algorithm 4.1 is completely data-driven and does not require any prior knowledge on the matrix Σ. Moreover, its computational burden remains small. We illustrate its efficiency in Section 6.

Extension to non-toroidal lattices

It is often artificial to consider the field X as stationary on a torus. However, we needed this hypothesis for deriving nonasymptotic properties of the estimator θ in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. In many applications, it is more realistic to assume that we observe a small window of a Gaussian field defined on the plane Z 2 . If we are unable to prove nonasymptotic risk bounds in this new setting. Nevertheless, Lakshman and Derin have shown in [START_REF] Lakshmanan | Valid parameter space for 2-D Gaussian Markov random fields[END_REF] that there is no phase transition within the valid parameter space for GMRFs defined on the plane Z 2 . Let us briefly explain what this means: consider a GMRF defined on a square lattice of size p, but only observed on a square lattice of size p ′ . The absence of phase transition implies the distribution of this field observed on this fixed window of size p ′ does not asymptotically depend on the bound conditions when p goes to infinity. Consequently, it is reasonable to think that our estimation procedure still performs well to the price of slight modifications. In the sequel, we assume that the field X is defined on Z 2 , but the data X still correspond to n independent observations of the field X on the window Λ of size p 1 ×p 2 . The conditional distribution of X[0,0] given the remaining covariates now decomposes as

X[0,0] = (i,j)∈Z 2 \{(0,0)} θ[i,j]X[i,j] + ǫ[0,0] , (16)
where θ[.,.] is an "infinite" matrix defined on Z 2 and where ǫ[0,0] is a centered Gaussian variable of variance σ 2 independent of (X[i,j]) (i,j)∈Λ\{(0,0)} . The distribution of the field X is uniquely defined by the function θ and positive number σ 2 . The set Θ +,∞ of valid parameter for θ is now defined using the spectral density function. We refer to Rue and Held [RH05] Sect.2.7 for more details.

Definition 5.1. A function θ : Z 2 → R belongs to the set Θ +,∞ if it satisfies the three following conditions:

1. θ[0,0] = 0.

For any

(i, j) ∈ Z 2 , θ[i,j] = θ[-i,-j].

For any

(ω 1 , ω 2 ) ∈ [0, 2π) 2 , 1 -(i,j)∈Z 2 θ[i,j] cos (iω 1 + jω 2 ) > 0.
Similarly, we define the set Θ +,∞,iso for the isotropic GMRFs on the lattices. As done in Section 2 for toroidal lattices, we now introduce the parametric parameter sets. For any model m ∈ M 1 , the set Θ +,∞ m refers to the subset of matrices θ in Θ +,∞ whose support is included in m. Analogously, we define the parameter set Θ +,∞,iso m corresponding to isotropic GMRFs.

We cannot directly extend the CLS empirical contrast γ n,p1,p2 (.) defined in (6) in this new setting because we have to take the edge effect into account. Indeed, if we want to compute the RR n 6798 conditional regression of X i [j1,j2], we have to observe all its neighbors with respect to m, i.e. {X i [j1+l1,j2+l2], (l 1 , l 2 ) ∈ m}. In this regard, we define the sublattice Λ m for any model m ∈ M 1 .

Λ m := {(i 1 , i 2 ) ∈ Λ , (m + (i 1 , i 2 )) ⊂ Λ} ,
where (m + (i, j)) denotes the set m of nodes translated by (i, j). For instance, if we consider the model m 1 with four nearest neighbors, the edge effect size is one and Λ m contains all the nodes that do not lie on the border. The model m 3 with 12 nearest neighbors yields an edge effect of size 2 and Λ m contains all the nodes in Λ, except those which are at a (euclidean) distance strictly smaller than 2 from the border.

For any model m ∈ M 1 , any θ ′ ∈ Θ +,∞ m , and any sublattice Λ ′ ⊂ Λ m , we define γ Λ ′ n,p1,p2 (.) as an analogous of γ n,p1,p2 (.) except that it only relies on the conditional regression of the nodes in Λ ′ .

γ Λ ′ n,p1,p2 (θ ′ ) := 1 nCard(Λ ′ ) n i=1 (j1,j2)∈Λ ′ X i [j1,j2] - (l1,l2)∈m θ ′ [l1,l2]X i [j1+l1,j2+l2] 2 .
Then, the CLS estimators θ Λ ′ m and θ Λ ′ ,iso m are defined by

θ Λ ′ m ∈ arg min θ ′ ∈Θ +,∞ m γ Λ ′ n,p1,p2 (θ ′ ) and θ Λ ′ ,iso m ∈ arg min θ ′ ∈Θ +,∞,iso m γ Λ ′ n,p1,p2 (θ ′ ) .
Contrary to θ m , the estimator θ Λm m is not necessarily unique especially if the size of Λ m is smaller than d m . Let us mention that it is quite classical in the literature to remove nodes to take edge effects or missing data into account (see e.g. [START_REF] Guyon | Random fields on a network[END_REF] Sect.4.3). We cannot use anymore fast Fourier transform for computing the parametric estimator. Nevertheless, the estimators θ Λ ′ m are still computationally amenable, since they minimizes a quadratic function on the closed convex set Θ +,∞ m . Suppose we are given a subcollection M of M 1 . We note Λ M the smallest sublattice among the collection of lattices Λ m with m ∈ M. In order to select the neighborhood m, we compute the estimators θ ΛM m and minimize the criteria γ ΛM n,p1,p2 ( θ ΛM m ) penalized by a quantity of the order d m /(nCard(Λ M )). We compute the quantities γ ΛM n,p1,p2 ( θ ΛM m ) instead of γ Λm n,p1,p2 ( θ Λm m ) since we want to compare the adequation of the models using the same data set.

We now describe a data-driven model selection procedure for choosing the neighborhood. It is based on the slope heuristic developed in the previous section.

Algorithm 5.1. (Data-driven penalization for non-toroidal lattice).

Compute the selected model m(N ) as a function of

N > 0 m(N ) ∈ arg min m∈M γ ΛM n,p1,p2 ( θ ΛM m ) + N d m nCard(Λ M ) . 2. Find N min > 0 such that the jump d b m " [ b Nmin] - " -d b m " [ b

Nmin]

+ " is maximal.

3. Select the model m = m(2 N min ). ). As for Algorithm 4.1, it is advised to introduce huge models in the collection M in order to better detect the dimension jump. However, when the dimension of the models increases the size of Λ m decreases and the estimator θ Λm m may become unreliable. The method therefore requires a reasonable number of data. In practice, Λ should not contain less than 100 nodes.

Simulation study

In the first simulation experiment, we compare the efficiency of our procedure with penalized maximum likelihood methods when the field is a torus. In the second and third studies, we consider the estimation of a Gaussian field observed on a rectangle. The calculations are made with R [R D08]. Throughout these simulations, we only consider isotropic estimators.

Isotropic GMRF on a torus

First, we consider X an isotropic GMRF on the torus Λ of size p = p 1 = p 2 = 20. There are therefore 400 points in the lattice. The number of observations n equals one and the conditional variance σ 2 is one. We introduce a radius r := √ 17. Then, for any number φ > 0, we define the p × p matrix θ φ as:

   θ φ [0,0] := 0 , θ φ [i,j] := φ if |(i, j)| t ≤ r and (i, j) = (0, 0) , θ φ [i,j] := 0 if |(i, j)| t > r .
In practice, we set φ to 0, 0.0125, 0.015, and 0.0175. Observe that these choices constrain θ φ 1 < 1. The matrix θ φ therefore belongs to the set Θ +,iso m10 of dimension 10 introduced in Definition 2.1.

First simulation experiment. In Section 3, we have advocated the use of the estimator θ instead of θ ρ , although theoretical results are only available for θ ρ with ρ < ∞. We recall that θ = θ ρ with ρ = ∞. We check in this simulation study that the performances of θ and θ ρ with different values of ρ are similar.

We consider the collection of neighborhoods M := {m 0 , m 1 , . . . , m 20 } whose maximal dimension d iso m20 is 21. The estimator θ iso is built using the CLS model selection procedure introduced in Algorithm 4.1. The estimators θ iso ρ are computed similarly, except that they are based on the parametric estimators θ iso m,ρ (Sect. 3) instead of θ iso m . The Gaussian field X with φ = 0.015 is simulated by using the fast Fourier transform. The quality of the estimations is assessed by the prediction loss function l(., .) defined in (11). The experiments are repeated 1000 times. For ρ = 2, 4, 8, we evaluate the risks E θ φ [l( θ iso , θ φ )] and E θ φ [l( θ iso ρ , θ φ )] as well as the corresponding empirical 95% confidence intervals by a Monte-Carlo method. We also estimate the risks of θ iso m and θ iso m,ρ for each model m ∈ M. It then allows to evaluate the oracle risks E θ φ [l( θ iso m * ,ρ , θ φ )] and the risk ratios

E θ φ [l( θ iso ρ , θ φ )]/E θ φ [l( θ iso m * ,ρ , θ φ )].
The risk ratio measures how well the selected model m iso performs in comparison to the "best" model m * . Moreover, the risk ratio roughly illustrates the oracle type inequality presented in Theorem 3.1. Indeed, the infimum inf m∈M [l(θ m,ρ , θ) + pen(m)] in ( 14) is a good measure of the risk E θ φ [l( θ iso m * ,ρ , θ φ )] as explained in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] Sect.4. The results are given in Table 1.

RR n 6798

Verzelen

They corroborate that the estimators θ iso and θ iso ρ perform similarly. Moreover, the risk ratios

E θ φ [l( θ iso ρ , θ φ )]/E θ φ [l( θ iso m * ,ρ , θ φ )] correspond to the ratios ρ 2 4 8 ∞ E θ φ [l( θ iso ρ , θ φ )] × 10 2 4.1 ± 0.1 4.2 ± 0.2 4.2 ± 0.1 4.2 ± 0.3 E θ φ [l( θ iso ρ , θ φ )]/E θ φ [l( θ iso m * ,ρ , θ φ )] 1.3 ± 0.1 1.3 ± 0.1 1.3 ± 0.1 1.3 ± 0.2
Table 1: First simulation study.

Estimates and 95% confidence intervals of the risks

E θ φ [l( θ iso , θ φ )], E θ φ [l( θ iso ρ , θ φ )]
, and of the ratios

E θ φ [l( θ iso , θ φ )]/E θ φ [l( θ iso m * , θ φ )] and E θ φ [l( θ iso ρ , θ φ )]/E θ φ [l( θ iso m * ,ρ , θ φ )
] with φ = 0.015 and ρ = 2, 4, 8.

Second simulation experiment. We compare the efficiency of the method with two alternative model selection procedures. For each of them, we use the collection M as in the previous experiment. The two alternative procedures are based on likelihood maximization. In this regard, we first define the parametric maximum likelihood estimator θ mle m for any model m ∈ M,

θ mle m , σ mle m := arg min θ ′ ∈Θ +,iso m ,σ ′ -L p (θ ′ , σ ′ , X) ,
where L p (θ ′ , X) stands for the log-likelihood at the parameter θ ′ . We then select a model m applying either an AIC-type criterion [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] or a BIC-type criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF]:

m AIC := arg min m∈M -2L p ( θ mle m , σ mle m , X) + 2d iso m , m BIC := arg min m∈M -2L p ( θ mle m , σ mle m , X) + log(p 2 )d iso m .
For short, we write θ AIC and θ BIC for the two obtained estimators θ mle b m AIC and θ mle b m BIC . Although AIC and BIC procedures are not justified in this setting, we still apply them as they are widely used in many frameworks. Their computation is performed efficiently using the fast Fourier transform described in Section 2.3.

The experiments are repeated 1000 times. The Gaussian field is simulated using the fast Fourier transform. The quality of the estimations is assessed by the prediction loss function l(., .). For any φ and any of these three estimators, we evaluate the risks

E θ φ [l( θ AIC , θ φ )], E θ φ [l( θ BIC , θ φ )]
, and E θ φ [l( θ iso , θ φ )] as well as the corresponding empirical 95% confidence intervals by a Monte-Carlo method. We also estimate the risk ratios E θ φ [l( θ iso , θ φ )]/E θ φ [l( θ iso m * , θ φ )] The results are given in Table 2.

The BIC criterion outperforms the other procedures when φ = 0, 0.0125, or 0.015 but behaves bad for a large φ. Indeed, the BIC criterion has a tendency to overpenalize the models. For the two first values of φ the oracle model in M is m 0 . Hence, overpenalizing increases the performance of estimation in this case. However, when φ increases, the dimension of the oracle model is larger and BIC therefore selects too small models.

In contrast, AIC and the CLS estimator exhibit similar behaviors. If we forget the case φ = 0 for which the oracle risk is 0, the risk of θ iso is close to the risk of the oracle model (the ratio is close to one). Hence, the neighborhood choice for θ iso is almost optimal.

INRIA φ × 10 2 0 1.25 1.5 1.75 E θ φ [l( θ AIC , θ φ )] × 10 2 1.2 ± 0.2 3.1 ± 0.2 4.3 ± 0.2 6.4 ± 0.2 E θ φ [l( θ BIC , θ φ )] × 10 2
0.01 ± 0.01 1.9 ± 0.1 3.7 ± 0.1 9.7 ± 0.3

E θ φ [l( θ iso , θ φ )] × 10 2 1.6 ± 0.2 3.2 ± 0.2 4.2 ± 0.1 7.2 ± 0.3 E θ φ [l( θ iso , θ φ )]/E θ φ [l( θ iso m * , θ φ )] +∞ 1.9 ± 0.7 1.3 ± 0.2 1.5 ± 0.3
Table 2: Second simulation study.

Estimates and 95% confidence intervals of the risks

E θ φ [l( θ AIC , θ φ )], E θ φ [l( θ BIC , θ φ )]
, and E θ φ [l( θ iso , θ φ )] and of the ratio

E θ φ [l( θ iso , θ φ )]/E θ φ [l( θ iso m * , θ φ )].
In conclusion, θ iso or θ AIC both exhibit good performances for estimating the distribution of a regular Gaussian field on a torus. The strength of our neighborhood selection procedure lies in the fact it easily generalizes to non-toroidal lattices as illustrated in the next section.

Isotropic Gaussian fields on Z 2

First simulation experiment. We now consider X an isotropic Gaussian field defined on Z 2 but only observed on a square Λ of sizes p = p 1 = p 2 = 20 or p = p 1 = p 2 = 100. This corresponds to the setting described in Section 5. The variance of X[0,0] is set to one and the distribution of the field is therefore uniquely defined by its correlation function ρ(k, l) := corr(X[k,l], X[0,0]). Again, the number of replications n is chosen to be one. In the first experiment, we use four classical correlation functions: exponential, spherical, circular, and Matérn (e.g. [Cre93] Sect.2.3.1 and [START_REF] Matérn | Spatial variation[END_REF]).

Exponential: ρ(k, l) = exp - d(k, l) r Circular: ρ(k, l) =      1 -2 π d(k,l) r 1 -d(k,l) r 2 + sin -1 d(k,l) r if d(k, l) ≤ r 0 else Spherical: ρ(k, l) = 1 -1.5 d(k,l) r + 0.5 d(k,l) r 3 if d(k, l) ≤ r 0 else Matérn: ρ(k, l) = 1 2 κ-1 Γ(κ) d(k, l) r κ K κ d(k, l) r ,
where d(k, l) denotes the euclidean distance from (k, l) to (0, 0) and K κ (.) is the modified Bessel function of order κ. In a nutshell, the parameter r represents the range of correlation, whereas κ may be regarded as a smoothness parameter for the Matérn function. In this simulation experiment, we set r to 3. When considering the Matérn model, we take κ equal to 0.05, 0.25, 0.5, 1, 2, and 4. The Gaussian fields are simulated using the function GaussRF in the library Random-Fields [START_REF] Schlather | RandomFields: Simulation and Analysis of Random Fields[END_REF]. For each of experiments, we compute the estimator θ iso based on Algorithm 5.1 with the collection M := {m ∈ M 1 , d iso m ≤ 18}. Since the lattice Λ is not a torus, methods based on likelihood maximization exhibit a prohibitive computational burden. Consequently, we do not use MLE in this experiment. We shall compare the efficiency of θ iso with a variogrambased estimation method.

We recall that the linear combination (i,j)∈Λ\{(0,0)} θ[i,j]X[i,j] is the kriging predictor of X[0,0] given the remaining variables (Equation (1)). A natural method to estimate θ in this spatial setting amounts to estimating the variogram of the observed Gaussian field and then performing ordinary kriging at the node (0, 0). More precisely, we first estimate the empirical variogram by applying the modulus estimator of Hawkes and Cressie (e.g. [START_REF] Cressie | Statistics for spatial data[END_REF] Eq.(2.2.8)) to the observed field of 400 points. Afterwards, we fit this empirical variogram to a variogram model using the reweighted least-squares suggested by Cressie [Cre85]. This procedure therefore requires the choice of a particular variogram model. In the first simulation study, we choose the model that has generated the data. Observe that this method is not adaptive since it requires the knowledge of the variogram model. In practice, we use Library geoR [START_REF] Ribeiro | geoR: a package for geostatistical analysis[END_REF] implemented in R [R D08] to estimate the parameters r, var(X[0,0]) and eventually κ of the variogram model. Then, we compute the estimator θ V by performing ordinary kriging at the center node of Λ. For each of these estimations, we assume that the variogram model is known. For computational reasons, we use a kriging neighborhood of size 11 × 11 that contains 120 points. Previous simulations have indicated that this neighborhood choice does not decrease the precision of the estimation. For the Matèrn model with κ = 2 and 4, the covariance is almost singular. There are sometimes inversion difficulties and we therefore use kriging neighborhood of respective size 7 × 7 and 3 × 3.

We again assess the performances of the procedures using the loss l(., .). Even if this loss is defined in (11) for a torus, the alternative definition (12) clearly extends to this non-toroidal setting. Consequently, the loss l( θ, θ) measures the difference between the prediction error of X[0,0] when using (i,j)∈Λ\{(0,0)} θ[i,j]X[i,j] and the prediction error of X[0,0] when using the best pre-

dictor E[X[0,0]|(X[i,j]) (i,j)∈Λ\{(0,0)} ].
In other words, l( θ, θ) is the difference of the kriging error made with the estimated parameters θ and the kriging error made with the true parameter θ.

The experiments are repeated 1000 times. For any of the four correlation models previously mentioned, we evaluate the risks E θ [l( θ iso , θ)] and E θ [l( θ V , θ)] by Monte-Carlo. In order to assess the efficiency of the selection procedure, we also evaluate the risk ratio

Risk.ratio = E θ [l( θ ΛM,iso b m , θ)] E θ [l( θ ΛM,iso m * , θ)] .
As in Section 6.1, the oracle risk E[l( θ ΛM,iso m * , θ)] is evaluated by taking the minimum of the evaluations of the risks E[l( θ ΛM,iso m , θ)] over all models m ∈ M. Results of the simulation experiment are given in Table 3 and4.

Observe that none of the fields considered in this study are GMRFs. Here, the GMRF models should only be viewed as a collection of approximation sets of the true distribution. This simulation experiment is in the spirit of Rue and Tjelmeland's study [START_REF] Rue | Fitting Gaussian Markov random fields to Gaussian fields[END_REF]. However, there are some major differences. Contrary to them, we perform estimation and not only approximation. Moreover, our lattice is not a torus. Finally, we use our prediction loss l(., .) to assess the performance, whereas they compare the correlation functions.

Comments on Tables 3 and4. In both tables, the ratio

E θ [l( θ ΛM,iso b m , θ)]/E θ [l( θ ΛM,iso m *
, θ)] stays close to one. Hence, the model selection is almost optimal from an efficiency point of view. In most of the cases, the estimator θ iso outperforms the estimator θ V based on geostatistical methods. This is particularly striking for the Matérn correlation model because in that case the computation of θ V requires the estimation of the additional parameter κ. Indeed, let us recall that the exponential model and the Matérn model with κ = 0.5 are equivalent. For κ = 0.5, the risk of θ V is 100 times higher when κ has to be estimated than when κ is known.

INRIA

Model

Exponential Circular Spherical E θ [l( θ V , θ] × 10 2 0.08 ± 0.01 9.1 ± 0.5 2.9 ± 0.1 E θ [l( θ iso , θ)] × 10 2 1.08 ± 0.01 6.5 ± 0.1 3.4 ± 0.1 Risk.ratio 3.6 ± 0.4 1.4 ± 0.1 1.6 ± 0.1 Table 3: Estimates and 95% confidence intervals of the risks E θ [l( θ V , θ)] and E θ [l( θ iso , θ)] and of Risk.ratio for the exponential, circular and spherical models with p = 20.

κ 0.05 0.25 0.5 1 E θ [l( θ V , θ)] × 10 3
91.8 ± 0.7 80.0 ± 0.2 18.0 ± 0.1 2.5 ± 0.1 E θ [l( θ iso , θ)] × 10 3 2.24 ± 0.01 0.62 ± 0.01 0.33 ± 0.01 0.08 ± 0.01

Risk.ratio

1.3 ± 0.1 1.7 ± 0.2 1.5 ± 0.2 1.3 ± 0.1 κ 2 4 E θ [l( θ V , θ)] × 10 4 6.3 ± 1.1 0.011 ± 0.001 E θ [l( θ iso , θ)] × 10 4 1.9 ± 0.1 0.17 ± 0.01 Risk.ratio 2.6 ± 0.2 1.1 ± 0.1
Table 4: Estimates and 95% confidence intervals of the risks E θ [l( θ V , θ)] and E θ [l( θ iso , θ)] and of Risk.ratio for Matérn model with p = 100.

Second simulation experiment. The kriging estimator θ V requires the knowledge or the choice of a correlation model. In the second simulation experiment, the correlation of X is the Matèrn function with range r = 3 and κ = 0.05. The size p of the lattice is chosen to be 100. We now estimate θ using different variogram models, namely the exponential, the circular, the spherical and the Matèrn model. The estimator θ iso for such a field was already considered in Table 4. The experiment is repeated 1000 times.

Model

Exponential Circular Spherical Matèrn E θ [l( θ V , θ)] × 10 3 48.3 ± 0.4 461 ± 16 293 ± 7 91.8 ± 0.7

Table 5: Estimates and 95% confidence intervals of the risks E θ [l( θ V , θ)] for Matérn model with κ = 0.05 when using the exponential, circular, spherical, and Matèrn models with p = 100.

Comments on Table 5. One observes that circular and spherical models yield worse performances than Matèrn model. In contrast, the exponential model behaves better. The choice of the variogram model therefore seems critical to get good performances. The model selection estimator θ iso (Table 4) exhibits a smaller risk than the exponential model.

Anisotropic Gaussian fields on Z 2

We still consider X a Gaussian field observed on a square Λ of size 100 × 100. Contrary to the previous study, the field is not assumed to be isotropic. To model the geometric anisotropy, we suppose that X is an isotropic field on a deformed lattice Λ ′ . The transformation consists in multiplying the original coordinates by a rotation R and a shrinking matrix T . For the sake of simplicity, we take the identity for R. The shrinking matrix T is defined by the anisotropy ratio (Ani.ratio). It corresponds to the ratio between the directions with smaller and greater continuity in the field X, i.e the ratio between maximum and minimum ranges. In this experiment, X follows a Matèrn correlation with range r = 3, κ = 0.05, 0.25, 0.5, 1, 2, and 4 and Ani.ratio=2 or 5. We compute the anisotropic estimator θ based on Algorithm 5.1 with the collection M := {m ∈ M 1 , d m ≤ 28}. As a benchmark, we also compute the variogram-based estimator θ V based on the Matèrn model. In order to compute θ V , we assume that we know the anisotropy ratio and the anisotropy directions. Observe that the estimator θ does not require any assumption on the form of anisotropy, while θ V uses the geometric parameters of the anisotropy.

The experiments are repeated 1000 times. We evaluate the risks E θ [l( θ V , θ)] and E θ [l( θ, θ)] and the risk ratio defined by

Risk.ratio = E θ [l( θ ΛM b m , θ)] E θ [l( θ ΛM m * , θ)] . κ 0.05 0.25 0.5 1 E θ [l( θ V , θ)] × 10 2 15.8 ± 0.1 13.9 ± 0.1 3.3 ± 0.1 0.30 ± 0.01 E θ [l( θ, θ)] × 10 2
0.65 ± 0.01 0.20 ± 0.01 0.089 ± 0.001 0.17 ± 0.01 Risk.ratio

1.2 ± 0.1 1.1 ± 0.1 1.1 ± 0.1 1.7 ± 0.2 κ 2 4 E θ [l( θ V , θ)] × 10 4
9.8 ± 0.1 0.020 ± 0.001 E θ [l( θ iso , θ)] × 10 4 45.0 ± 0.1 4.3 ± 0.1 Risk.ratio 2.9 ± 0.2 22.3 ± 1.7

Table 6: Estimates and 95% confidence intervals of the risks E θ [l( θ V , θ)] and E θ [l( θ, θ)] and of Risk.ratio for Matérn model and Ani.ratio= 2. κ 0.05 0.25 0.5 1 E θ [l( θ V , θ)] × 10 2 11.2 ± 0.1 14.9 ± 0.1 3.7 ± 0.1 2.9 ± 0.1 E θ [l( θ, θ)] × 10 2 0.66 ± 0.1 0.40 ± 0.01 0.081 ± 0.001 0.14 ± 0.01 Risk.ratio

1.1 ± 0.1 1.1 ± 0.1 1.2 ± 0.1 3.4 ± 0.8 κ 2 4 E θ [l( θ V , θ)] × 10 4 30.6 ± 0.1 0.22 ± 0.01 E θ [l( θ iso , θ)] × 10 4 38.0 ± 0.1 39.6 ± 0.1 Risk.ratio 2.1 ± 0.1 9.0 ± 1.4
Table 7: Estimates and 95% confidence intervals of the risks E θ [l( θ V , θ)] and E θ [l( θ, θ)] and of Risk.ratio for Matérn model and Ani.ratio= 5. 6 and7. Except for the cases κ = 2, 4, the estimator θ performs better than the variogram-based estimator θ V , although θ V uses the true anisotropy parameters. For κ = 4, the neighborhood selection is no performed efficiently (the risk ratio is large).

Comments on Tables

INRIA 7 Discussion

In this paper, we have extended a neighborhood selection procedure introduced in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. On the one hand, an algorithm is provided for tuning the penalty in practice. On the other hand, the new method also handles non-toroidal lattices. The computational complexity remains reasonable even when the size of the lattice is large.

In the case of stationary fields on a torus, our neighborhood selection procedure exhibits a computational burden and statistical performances analogous to the AIC procedure. Even if AIC has not been analyzed from an efficiency point of view, this suggests that AIC may achieve an oracle inequality in this setting. Moreover, we have empirically checked that θ performs almost as well as the oracle model since the oracle ratio E[l( θ, θ)]/E[l( θ m * , θ)] remains close to one.

The strength of this neighborhood selection procedure lies in the fact it easily extends to non-toroidal lattices. We have illustrated that our method often outperforms variogram-based estimation methods in terms of the mean-squared prediction error. Moreover, the procedure behaves almost as well as the oracle. In contrast, variogram-based procedures may perform well for some covariances structure but also yield terrible results for other covariance structures. These results illustrate the adaptivity of the neighborhood selection procedure.

In many statistical applications, Gaussian fields (or Gaussian Markov random fields) are not directly observed. For instance, Aykroyd [START_REF] Aykroyd | Bayesian estimation for homogeneous and inhomogeneous gaussian random fields[END_REF] or Dass and Nair [START_REF] Dass | Edge detection, spatial smoothing, and image reconstruction with partially observed multivariate data[END_REF] use compound Gaussian Markov random fields to account for non stationarity and steep variations. The wavelet transform has emerged as a powerful tool in image analysis. The wavelet coefficients of an image are sometimes modeled using hidden Markov models [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF][START_REF] Portilla | Image denoising using scale mixtures of Gaussians in the wavelet domain[END_REF]. More generally, the success of the GMRFs is mainly due to the use of hierarchical models involving latent GMRFs [START_REF] Rue | Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations[END_REF]. The study and the implementation of our penalization strategy for selecting the complexity of latent Markov models is an interesting direction of research.

Proofs

Let us introduce some notations that shall be used throughout the proofs. For any 1 ≤ k ≤ n, the vector X v k denotes the vectorialized version of the k-th sample of X. Moreover, X v is the matrix of size p 1 p 2 × n of the n realisations of the vector X v k . Throughout these proofs, L, L 1 , L 2 denote constants that may vary from line to line. The notation L(.) specifies the dependency on some quantities. Finally, the γ(.) function stands for an infinite sampled version of the CLS criterion γ n,p1,p2 (.): γ(.) := E[γ n,p1,p2 (.)].

Proof of Lemma 2.1

Let us provide an alternative expression of γ n,p1,p2 (θ ′ ) in term of the factor C(θ ′ ) and the empirical covariance matrix X v X v * .

γ n,p1,p2 (θ ′ ) = 1 np 1 p 2 tr (I p1p2 -C(θ ′ ))X v X v * (I p1p2 -C(θ ′ )) . (17)
for any 1 ≤ i ≤ p 1 and 1 ≤ j ≤ p 2 .

This lemma is proved as in [RH05] Sect.2.6.2 to the price of a slight modification that takes into account the fact that P is orthogonal and not unitary. The difference comes from the fact that contrary to Rue and Held we also assume that C(θ) is symmetric. Lemma 8.1 states that all symmetric block circulant matrices are simultaneously diagonalizable. Observe that for any 1 ≤ i ≤ p 1 and 1 ≤ j ≤ p 2 , it holds that D(θ)[(i-1)p2+j,(i-

1)p2+j] = λ[i,j](θ) since θ[k,l] = θ[p1-k,p2-l]. Hence, Expression (17) becomes γ n,p1,p2 (θ ′ ) = 1 np 1 p 2 p1 i=1 p2 j=1 [1 -λ[i,j](θ)] 2 n k=1 [P * X v k (X v k ) * P ] [(i-1)p2+j,(i-1)p2+j] ,
where X v k is the vectorialized version of the k-th observation of the field X. Straightforward computations allow us to prove that the quantities

(P * X v k (X v k ) * P ) [(i-1)p2+j,(i-1)p2+j] + (P * X v k (X v k ) * P ) [(p1-i-1)p2+p2-j,(p1-i-1)p2+p2-j] and 1 √ p 1 p 2 λ[i,j](X v k )λ[i,j](X v k ) + 1 √ p 1 p 2 λ[p1-i,p2-j](X v k )λ[p1-i,p2-j](X v k )
are equal for any 1 ≤ i ≤ p 1 and 1 ≤ j ≤ p 2 . Here, the entries of the matrix λ(.) are taken modulo p 1 and p 2 and the entries of [P * X v k (X v k ) * P ] are taken modulo p 1 p 2 . The result of Lemma 2.1 follows.

Proof of Proposition 4.1

Proof of Proposition 4.1. We only consider the anisotropic case, since the proof for isotropic estimation is analogous. For any model m ∈ M 1 , we define ∆(m, m ′ ) := γ n,p,p θ m,ρ + pen(m)γ n,p,p θ m ′ ,ρpen(m ′ ) .

We aim at showing that with large probability, the quantity ∆(m, m ′ ) is positive for all small dimensional models m. Hence, we would conclude that the dimension of m is large. In this regard, we bound the deviations of the differences

γ n,p,p θ m,ρ -γ n,p,p θ m ′ ,ρ = γ n,p,p θ m,ρ -γ n,p,p (θ m,ρ ) + γ n,p,p (θ m,ρ ) -γ n,p,p (θ) + γ n,p,p (θ) -γ n,p,p θ m ′ ,ρ .
Lemma 8.2. Let K 2 be some universal constant that we shall define in the proof. With probability larger than 3/4,

γ n,p,p (θ) -γ n,p,p (θ m,ρ ) ≤ K 2 2 ρ 2 ϕ max (Σ) d m ∨ 1 np 2 and γ n,p,p (θ m,ρ ) -γ n,p,p θ m,ρ ≤ K 2 2 ρ 2 ϕ max (Σ) d m np 2 for all models m ∈ M 1 .
INRIA Lemma 8.3. Assume that p is larger than some numerical constant p 0 . With probability larger than 3/4, it holds that

γ n,p,p (θ) -γ n,p,p ( θ m ′ ,ρ ) ≥ K 3 σ 2 ϕ min I p 2 -C(θ) ∧ ρ -ϕ max I p 2 -C(θ) d m ′ np 2 ,
where K 3 is a universal constant defined in the proof.

Let us take K 1 to be exactly K 3 . Gathering the two last lemma with Assumption (15), there exists an event Ω of probability larger than 1/2 such that

∆(m, m ′ ) ≥ σ 2 np 2 K 1 ηd m ′ ϕ min I p 2 -C(θ) ∧ ρ -ϕ max I p 2 -C(θ) -K 2 (dm∨1)ρ 2 ϕmin(I p 2 -C(θ))
, for all models m ∈ M 1 . Thus, conditionally to Ω, ∆(m, m ′ ) is positive for all models m ∈ M 1 that satisfy

d m ∨ 1 d m ′ ≤ K 3 η K 2 ρ 2 ϕ min I p 2 -C(θ) ϕ min I p 2 -C(θ) ∧ ρ -ϕ max I p 2 -C(θ) .
By Lemma 8.7 in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF], the dimension d m ′ is larger than 0.5[ np 2 ∧ (p 2 -1)]. We conclude that

d b mρ ∨1 ≥ np 2 ∧ p 2 -1 K 3 η K 2 ρ 2 ϕ min I p 2 -C(θ) ϕ min I p 2 -C(θ) ∧ ρ -ϕ max I p 2 -C(θ) ,
with probability larger than 1/2.

Proof of Lemma 8.2. In the sequel, γ n,p,p (.) denotes the difference γ n,p,p (.) -γ(.). Given a model m, we consider the difference γ n,p,p (θ)γ n,p,p (θ m,ρ ) = γ n,p,p (θ)γ n,p,p (θ m,ρ )l(θ m,ρ , θ) .

Upper bounding the difference of γ n,p,p therefore amounts to bounding the difference of γ n,p,p . By definition of γ n,p,p and γ, it expresses as

γ n,p,p (θ) -γ n,p,p (θ m,ρ ) = 1 p 2 tr (I p 2 -C(θ)) 2 -(I p 2 -C(θ m,ρ )) 2 X v X v * -Σ .
γ n,p,p (θ) -γ n,p,p (θ m,ρ ) = 1 p 2 tr (D(θ m,ρ ) -D(θ)) 2I p 2 -D(θ) -D(θ m,ρ ) D Σ YY * -I p 2 , (19)
where the matrix Y is defined as P √ Σ -1 X v P * . Its components follow independent standard Gaussian distributions. Since the matrices involved in (19) are diagonal, Expression (19) is a linear combination of centered χ 2 random variables. We apply the following lemma to bound its deviations.

Verzelen

T := D i=1 a i Y 2 i -1 .
Then, the following deviation inequality holds for any positive x

P T ≥ 2 a 2 √ x + 2 a ∞ x ≤ e -x .
This result is very close to Lemma 1 of Laurent and Massart in [START_REF] Laurent | Adaptive estimation of a quadratic functional by model selection[END_REF]. The only difference lies in the fact that they constrain the coefficients a i to be non-negative. Nevertheless, their proof easily extends to our situation. Let us define the matrix a of size n × p 2 as

a i [j] := D Σ [i,i] (D(θ m,ρ )[i,i] -D(θ)[i,i]) (2 -D(θ[i, i] -D(θ m,ρ )[i, i]) np 2 ,
for any 1 ≤ i ≤ n and any 1 ≤ j ≤ p 2 . Since the matrices I -C(θ) and I -C(θ m,ρ ) belong to the set Θ + ρ , their largest eigenvalue is smaller than ρ. By Definition (11) of the loss function l(., .), a 2 ≤ 2ρ ϕ max (Σ)l(θ m,ρ , θ)/(np 2 ) and a ∞ ≤ 4ρ 2 ϕ max (Σ)/(np 2 ). By Applying Lemma 8.4 to Expression (19), we conclude that

P γ n,p,p (θ) -γ n,p,p (θ m,ρ ) ≥ l(θ m,ρ , θ) + 12ρ 2 ϕ max (Σ) np 2 x ≤ e -x ,
for any x > 0. Consequently, for any K > 0, the difference of γ n,p,p (.) satisfies

γ n,p,p (θ) -γ n,p,p (θ m,ρ ) ≤ K 2 ρ 2 ϕ max (Σ) d m ∨ 1 np 2 ,
simultaneously for all models m ∈ M 1 with probability larger than 1 -m∈M1\∅ e -K(dm∨1)/24 . If K is chosen large enough, the previous upper bound holds on an event of probability larger than 7/8. Let us call K ′ 2 such a value.

Let us now turn to the second part of the result. As previously, we decompose the difference of empirical contrasts γ n,p,p (θ m,ρ )γ n,p,p θ m,ρ = γ n,p,p (θ m,ρ )γ n,p,p θ m,ρl θ m,ρ , θ m,ρ Arguing as in the proof of Theorem 3.1 in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF], we obtain an upper bound analogous to Eq.( 49) in [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] γ n,p,p (θ m,ρ )γ n,p,p θ m,ρ ≤ l( θ m,ρ , θ m,ρ ) + ρ 2 sup ∧ α 2 1+α/2 )]. Thus, there exists some numerical constant α 0 such that the upper bound (20) with α = α 0 holds simultaneously for all models m ∈ M 1 \ ∅ with probability larger than 7/8. Choosing K 2 to be the supremum of K ′ Let us consider Ψ i1,j1 , . . . , Ψ i d m ′ ,j d m ′ a basis of the space Θ m ′ defined in Eq.( 14) of [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. Let α be a positive number that we shall define later. We then introduce θ ′ as

R∈B H ′ m 2 ,
θ ′ := ϕ min I p 2 -C(θ) α p 2 d m ′ k=1 tr C (Ψ i k ,j k ) ZZ * Ψ i k ,j k .
Since θ is assumed to belong to Θ + m ′ ,ρ , the parameter θ ′ belongs to (Θ + m ′ ,ρθ) if Let us now bound the deviations of the two random variables involved in ( 21) and ( 23) by applying Markov's and Tchebychev's inequality. For the sake of simplicity, we assume that d m ′ is smaller than (p 2 -2p)/2. In such a case, all the nodes in m ′ are different from their symmetric in Λ. We omit the proof for d m ′ larger than (p 2 -2p)/2 because the approach is analogous but the computations are slightly more involved. Straightforwardly, we get

INRIA

In the sequel, we assume that p is larger than some universal constant p 0 , which ensures the dimension d m ′ to be larger than 4L ′ 1 . Gathering (21) with the upper bound on tr C(θ ′ )ZZ * yields where K 3 is an universal constant.

θ ′ 1 ≤ 2 √ 2αϕ min I p 2 -C(θ) d m ′

Figure 1 :

 1 Figure 1: (a) Model m 1 with first order neighbors. (b) Model m 2 with second order neighbors. (c) Model m 3 with third order neighbors.

  Algorithm 4.1. (Data-driven penalization with slope heuristic). Let M be a subcollection of M 1 . 1. Compute the selected model m(N ) as a function of N > 0

4. Compute the estimator θ Λ c m b m . INRIA

 . This procedure straightforwardly extends to the case of isotropic GMRFs estimation by replacing m(N ) by m iso (N ) and d m by d iso m . For short, we write θ (resp.θ iso ) for θ Λ c m b m (resp. θ Λ c m ,iso b m

  The matrices Σ, (I p 2 -C(θ)), and (I p 2 -C(θ m,ρ )) are symmetric block circulant. By Lemma 8.1, they are jointly diagonalizable in the same orthogonal basis. If we note P an orthogonal matrix associated to this basis, then C(θ m,ρ ), C(θ), and Σ respectively decompose in C(θ m,ρ ) = P * D(θ m,ρ )P , C(θ) = P * D(θ)P and Σ = P * D(Σ)P , where the matrices D(θ m,ρ ), D(θ), and D(Σ) are diagonal.

2

  and 2L 1 (1 + α 0 /2) allows to conclude.Proof of Lemma 8.3. Thanks to the definition (17) of γ n,p,p (.) we obtainγ n,p,p (θ)γ n,p,p ( θ m ′ ,ρ ) = 1 p 2 sup θ ′ ∈Θ + m ′ ,ρ tr (C(θ ′ ) -C(θ)) 2I p 2 -C(θ) -C(θ ′ ) ΣZZ * ,where the p 2 × n matrix Z is defined byZ := √ Σ -1 X v .We recall that the matrices Σ, C(θ) and C(θ ′ ) commute since they are jointly diagonalizable by Lemma 8.1. Let (Θ + m ′ ,ρθ) be the set Θ + m ′ ,ρ translated by θ. Since C(θ) + C(θ ′ ) = C(θ + θ ′ ), we lower bound the difference of γ n,p,p (.) as followsγ n,p,p (θ)γ n,p,p ( θ m ′ ,ρ ) = 1 p 2 sup θ ′ ∈ " Θ + m ′ ,ρ -θ " 2σ 2 tr C(θ ′ )ZZ *tr C(θ ′ ) 2 ΣZZ * ≥ σ 2 p 2 sup θ ′ ∈ " Θ + m ′ ,ρ -θ " 2tr C(θ ′ )ZZ *ϕ -1 min I p 2 -C(θ) tr C(θ ′ ) 2 ZZ * .

ϕVerzelen.

  max [C(θ ′ )] ≤ ϕ min I p 2 -C(θ) and ϕ min [C(θ ′ )] ≥ -ρ + ϕ max I p 2 -C(θ) . The largest eigenvalue of C(θ ′ ) is smaller than θ ′1 whereas its smallest eigenvalue is larger thanθ ′ 1 . Let us upper bound the l 1 norm of θ ′ :θ ′ 1 = 2ϕ min I p 2 -C(θ) α p 2 d m ′ k=1 tr C (Ψ i k ,j k ) ZZ * ≤ 2 α p 2 ϕ min I p 2 -C(θ) d m ′ tr C(θ ′ )ZZ * . (21) Hence, θ ′ belongs to (Θ + m ′ ,ρθ) if θ ′ 1 ≤ ϕ min I p 2 -C(θ) ∧ ρϕ max I p 2 -C(θ) .(22)Thus, we get the lower boundγ n,p,p (θ)γ n,p,p ( θ m ′ ,ρ ) ≥ σ 2 p 2 2tr C(θ ′ )ZZ *ϕ -1 min I p 2 -C(θ) tr C(θ ′ ) 2 ZZ * ,(23)as soon as Condition (22) is satisfied.

E

  tr C(θ ′ )ZZ * = 4αϕ min I p 2 -C(θ) d m ′ n ,since the neighborhood m ′ only contains points (i, j) whose symmetric (-i, -j) is different. A cumbersome but pedestrian computation leads to the upper boundvar tr C(θ ′ )ZZ * ≤ L 1 α 2 ϕ 2 min I p 2 -C(θ) d m ′ n 2 ,where L 1 is a numerical constant. Similarly, we upper bound the expectation of tr C(θ′ ) 2 ZZ * E tr C(θ ′ ) 2 ZZ * ≤ L 2 α 2 ϕ 2 min I p 2 -C(θ) d m ′ n .Let us respectively apply Tchebychev's inequality and Markov's inequality to the variables tr C(θ ′ )ZZ * and tr C(θ ′ ) 2 ZZ * . Hence, there exists an event Ω of probability larger than 3/4 such that2tr C(θ ′ )ZZ *ϕ -1 min I p 2 -C(θ) tr C(θ ′ ) 2 ZZ * ≥ ϕ min I p 2 -C(θ) θ ′ )ZZ * ≤ 4αϕ min I p 2 -C(θ)

  2αϕ min I p 2 -C(θ) , since d m ′ ≤ p √ n. If 2 √ 2α is smaller than 1 ∧ ρϕ max I p 2 -C(θ) ϕ -1 min I p 2 -C(θ) , then Condition (22) is fulfilled on the event Ω and it follows from (23) thatP γ n,p,p (θ)γ n,p,p ( θ m ′ ,ρ ) ≥ 4σ 2 ϕ min I p 2 -C(θ) d m ′ np 2 αα 2 L ′ 2 (I p 2 -C(θ))4ϕmin(I p 2 -C(θ)) , we getP γ n,p,p (θ)γ n,p,p ( θ m ′ ,ρ ) ≥ K 3 σ 2 ϕ min I p 2 -C(θ) ∧ ρϕ max I p 2 -C(θ) d m ′ np 2 ≥ 3 4 ,

  There exists an orthogonal matrix P which simultaneously diagonalizes every p 1 p 2 × p 1 p 2 symmetric block circulant matrices with p 2 × p 2 blocks. Let θ be a matrix of size p 1 × p 2

	This is justified in [Ver09] Sect.2.2. such that C(θ) is symmetric. The matrix D(θ) = P * C(θ)P is diagonal and satisfies p1 k=1 p2 l=1 Lemma 8.1. Verzelen D(θ)[(i-1)p2+j,(i-1)p2+j] = θ[k,l] cos [2π(ki/p 1 + lj/p 2 )] ,

  Lemma 8.4. Let (Y 1 , . . . , Y D ) be i.i.d. standard Gaussian variables. Let a 1 , . . . , a D be fixed numbers. We set

				D
	a ∞ := sup i=1,...,D	|a i |,	a 2 2 :=	i=1	a 2 i
	Let T be the random variable defined by

  m 2 1 p 2 tr RD Σ YY * -I p 2 H ′ m 2 ,m 2 is defined in the proof of Lemma 8.2 in[START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. Its precise definition is not really of interest in this proof. Coming back to the difference of γ n,p,p (.), we getγ n,p,p (θ m,ρ )γ n,p,p θ m,ρ ≤ ρ 2We consecutively apply Lemma 8.3 and 8.4 in[START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF] to bound the deviation of this supremum. Hence, for any positive number α, γ n,p,p (θ m,ρ )γ n,p,p θ m,ρ ≤ L 1 (1 + α/2)ρ 2 ϕ max (Σ)

	  	sup R∈B H ′ m 2 ,m 2	1 p 2 tr RD Σ YY * -I p 2	 2  	.
				d m np 2 .		(20)
	with probability larger than 1 -exp[-L 2	√ d m ( α √ 1+α/2

2

.

INRIA

The set B

RR n 6798

Acknowledgements

I am grateful to Pascal Massart and Liliane Bel for many fruitful discussions. I also thank the referees and the associate editor for their suggestions that led to an improvement of the manuscript.