F Valdivieso
email: fvaldivieso@emse.fr

M Tupin
email: tupin@emse.fr

M Pijolat

M Soustelle

Oxidation kinetics of ZrNbO in steam: differences between the pre-and post-transition stages

Keywords: ZrNbO, oxidation, kinetics, water vapour, steady state, rate-limiting step

Zirconium based alloys (and particularly Zircaloy-4) are widely used as cladding material of fuel rods in water-cooled nuclear reactors. However, advanced alloys are now used for longer operation time in more severe conditions. One of these new alloys is a ZrNbO, which contains 1% of Nb as alloying element (instead of Sn in the case of Zircaloy-4).

This ZrNbO alloy, like many zirconium alloys, undergoes a kinetic transition, which is a sharp increase in the oxidation rate when the oxide thickness exceeds a critical value. In the pre-transition stage, the kinetic oxidation curves are approximately parabolic (which is generally interpreted by a diffusion controlling step), and tend towards a quasi-linear behaviour after the transition. Recent studies have shown that the oxygen pressure has an accelerating effect on the oxidation rate of ZrNbO; a platinum layer deposited on the oxide surface has the same effect, both in oxygen and in water vapour. Thus it has been suggested that the oxidation in the pre-transition stage could be described by a mixed reaction-diffusion kinetics. Due to the lack of consistent data on the effect of water vapour on the oxidation rate of ZrNbO before and after the transition, and to confirm (or not) these interpretations, we have studied the oxidation of ZrNbO by isothermal thermogravimetry at 550°C, under controlled partial pressures of water vapour (13 to 80 hPa) and hydrogen (10 hPa). The aim is to verify if the oxidation proceeds in a steady state and if there is a rate-limiting step in one (or both) stages, and to clearly evidence the differences between the pre-and post transition stages from the kinetic point of view.

TG-DSC experiments have shown that the system is in a steady state from the beginning of the oxidation, in the pre-and post-transition stages. Then, the existence of a ratelimiting step was verified using an experimental method based on temperature or pressure jumps: it has been concluded that this assumption is valid in both stages, which is not compatible with the assumption of a mixed reaction-diffusion controlled rate). Moreover, it comes out that the kinetic behaviour is different before and after the transition: the influence of temperature jumps is greater in pre-transition, whereas the effect of water vapour partial pressure is more pronounced in post-transition (nevertheless, an accelerating effect is also observed before the transition). No influence of hydrogen partial pressure has been observed. Besides, the higher the Nb content in the alloy, the higher the oxidation rate (in pre-transition). A mechanism has been proposed to account for the results obtained in pretransition, involving the diffusion of adsorbed species in the porous part of the oxide layer as rate-determining step.

Introduction

Despite numerous works dedicated to the oxidation kinetics of zirconium alloys by pressurized water, water vapour or oxygen, the mechanisms and the rate-limiting steps for zirconia growth are not already established. It is well known that after a pseudo-parabolic first period [1][2][3], a kinetic transition occurs when the thickness of the oxide layer exceeds a critical value, which corresponds to an increase in the oxidation rate, associated with the appearance of cracks and pores in the oxide layer [2,4,5]. In the post-transition stage, the oxidation curves are considered to be either linear [4] or to result from cyclic periods of increasing and decreasing rate [6].

In a previous work dedicated to the oxidation of Zircaloy-4 by water vapour [7], we have shown that great differences exist between the pre and post-transition stages, from the kinetic point of view: before the transition, the oxidation rate is controlled by the diffusion of oxygen vacancies in the oxide layer, as usually suggested in the literature data, whereas in the post-transition stage, the assumption of the rate-limiting step is no more valid.

As far as the oxidation of Nb-containing alloys is concerned, the shape of the kinetic curves (weight gain versus time) is qualitatively similar to that obtained with zircaloys. Nevertheless, Nb-containing alloys are oxidised faster than recrystallised Zircaloy in oxygen and in water vapour [4,[8][9][10][11][12], at the beginning of the oxidation, but the kinetic transition is not so sharp in Nb-containing alloys and it occurs after a longer time. In autoclave, different results are obtained, these alloys showing in most cases a better resistance to corrosion [START_REF] Sabol | Zirconium in the Nuclear Industry, 8th International Symposium[END_REF][START_REF] Mardon | Proceedings of the International Meeting on Light Water Reactor Fuel Performance[END_REF].

Another difference between the two alloys is related to the effect of platinum on the oxidation rate: if platinum is deposited on the surface of an oxide layer grown on ZrNbO, in pre-transition, the oxidation rate increases, in oxygen [10] or in water vapour [11], whereas Pt has no effect on the oxidation of Zy4 [11,12].

Moreover, it has been shown that the oxidation rate depends on the niobium concentration in the zirconium alloy [START_REF] Choo | [END_REF], and on the oxygen pressure [17].

In order to account for these results, the pre-transition stage of the oxidation of ZrNbO is suggested to be controlled by a mixed diffusion-reaction regime, and therefore in steady conditions [10,[START_REF] Frichet | 12 th Symposium on Zirconium in the Nuclear Industry[END_REF].

Finally, due to the lack of consistent data on the effect of water vapour on the oxidation rate before and after the transition, and in order to confirm (or not) these interpretations, we have studied the oxidation of a Nb-containing zirconium alloy, ZrNbO alloy, at 530°C, following the same methodology as in the case of Zircaloy-4 [7]. This approach is principally based on the verification of kinetic assumptions that are usually made a priori. The oxide scale growth is assumed to be described by a series of elementary steps, one of them being rate-determining. This assumption implies that (i) the system is in a steady state and (ii) the oxidation rate can be written as:

( ) ( ) t E . P , T dt d i Φ = α (1)
in which α is the fractional conversion, Φ is a rate per unit area (mol.s -1 .m -2 ), it depends on the nature of the rate-limiting step (diffusion, interface reaction), it is independent on time but may be a function of temperature T and partial pressures of the gases P i . E(t) (m 2 .mol -1 ) is a function related to the extent of the reaction zone where the rate-limiting step is located. In this paper, the assumptions (i) and (ii) will be verified experimentally, then the variations of Φ with P H2O and P H2 will be obtained, using a method based on pressure jumps. These results will allow to clearly identify the differences between the pre and post-transition stages, and a mechanism will be proposed for the pre-transition.

Experimental results

Experimental. A 0.4 mm thick sheet of recrystallised ZrNbO alloy was used, and samples were cut to 10 x 10 mm for thermogravimetric experiments, and 15 x 15 mm for simultaneous thermogravimetry and calorimetry experiments. The alloy composition is indicated in Table I. The sample surface was cleaned with an equimolar solution of ethanol and acetone in ultrasonic waves, then with pure ethanol and dried with compressed air.

Fe (ppm)

O (ppm) Nb (%) 354 1303 1.03 Table I : Composition of the ZrNbO alloy The oxidation curves were obtained in isothermal and isobaric conditions at 500 and 530°C with a symmetrical thermoanalyser SETARAM TAG16, under a flowing mixture of water vapour and hydrogen in helium. The flowrate of the gases was controlled by massflowmeters (Brooks 5850S), the total flowrate being 2.3 l.h -1 , and the partial pressure of hydrogen being usually equal to 10 hPa. The water vapour partial pressure was fixed in the range 13-80 hPa, using thermoregulated baths. It was controlled using humidity sensors (Transmicor 241-242 Coreci), placed at the inlet and outlet of the thermobalance furnace.

Jumps in hydrogen pressure were carried out by changing the flowmeter setpoint, whereas jumps in water vapour pressure were performed by switching the gaseous flow from one water bath to another one maintained at a different temperature.

The simultaneous thermogravimetry and calorimetry experiments were performed using a Setaram TG/DSC 111, under flowing mixtures of water vapour and hydrogen in helium.

The morphology of the oxide layers was observed by scanning electron microscopy (SEM DSM 960A Zeiss). Shape of the kinetic curves. Fig. 1 represents a kinetic curve giving the mass gain per unit area m/S(t) versus time, and its derivative ( dt d(m/S) ), obtained at 530°C (after an initial temperature rise of 30°C/min from room temperature to 530°C), under 10 hPa in hydrogen and 13 hPa in water vapour. The minimum of the rate corresponds to an oxide thickness (calculated from the weight gain) of about 10 µm. The curve appears to be approximately parabolic in the pre-transition stage, whereas the rate is constant after the transition. (13 hPa) and hydrogen (10 hPa), showing the pre-and post-transition stages. scaling factor allowing to superimpose the two curves could be found, which shows that the reacting system is in a steady state before the transition. Several experiments were performed up to the transition, confirming the results of Fig. 2. Unfortunately, only the very beginning of the transition was reached in TG/DSC experiments, the linear regime characteristic of the post-transition stage could not be reached, even for experiments lasting more than one month. Consequently, the steady state approximation could not be verified after the transition.

Fig. 1 : Weight gain (⎯) and its derivative (---) versus time for ZrNbO at 530°C in water vapour

Rate-determining step assumption.

When a rate-limiting step can be assumed, Eq. ( 1) gives the variations of the oxidation rate with the intensive variables (T, P i …) and the time. In isobaric and isothermal conditions, Φ(T,P i ) is constant, then a sudden change (jump) in temperature or partial pressure during an experiment will then lead to a change in Φ, while E(t) will remain approximately constant (provided that the time necessary for the T or P change is short enough).

Thus the ratio of the rates measured after (at the right side) and before (at the left side) the jump is equal to l r Φ Φ (E(t) being eliminated in the ratio). Performing a series of similar jumps at various reaction times will lead to a series of l r Φ Φ ratios, which must be identical if Eq. ( 1) can be applied. This method, that we have named the "ΦE test", has been successfully used in previous works [7,[17][18][19][20][21][22]. The results are indicated on Fig. 3, for sudden changes in temperature from 500 to 530°C (P H 2 O = 13 hPa, P H 2 = 10 hPa) in pre-transition stage (Fig. 3a), and from 530 to 500°C in post-transition (Fig. 3b). Considering the experimental errors bars, the ratio keeps a constant value (2.17 ± 0.1 for temperature jumps from 500 to 530°C (Fig. 3a) and 2.13 ± 0.1 for temperature jumps from 530 to 500°C (Fig. 3b)) during the pre-transition stage, up to about 6 µm. Then the ratio decreases between 6 µm and 12 µm, during the transition, and becomes again constant after 17 µm, in the post-transition stage (characterized by a linear regime with a constant rate).

Consequently, it can be concluded that the "ΦE test" is validated in the pre-transition domain; as the steady state assumption is also verified, there exists a rate-limiting step in this stage (moreover, a mixed diffusion-reaction controlling rate is not in agreement with the results of the "ΦE test", because in such a case the ratios of the rates would not be constant [11,20]). The "ΦE test" is also verified in the post-transition stage, but the value of the ratio is clearly lower than its value in the pre transition stage (1.55 ± 0.08 instead of 2.13 ± 0.1). Obviously, a change in the rate-limiting step occurs between the pre and post-transition stages, since the shape of the curves and the values of the rate ratios are different.

Effect of H 2 O and H 2 partial pressures.

The variations of Φ with a gas partial pressure P can easily be obtained by performing sudden changes, from P 0 to P, at a given thickness X, the ratio of the rates being then equal to ( ) ( ) 17,19,21]. The experiments were carried out at 500°C, with 10 hPa in hydrogen and a water vapour pressure varying from 13 hPa (P 0 ) to 80 hPa after the jump. For the hydrogen pressure jumps, the water pressure was fixed at 13 hPa and the hydrogen pressures changed from 10 hPa (P 0 ) to 40 hPa. The variations of Φ with P H 2 O in pre-transition are represented on Fig. 4a: the water vapour pressure has a slight accelerating effect. In the post-transition stage, jumps in water vapour pressure have been performed at an oxide thickness equal to 25 µm (in the linear part of the curves), from 13 hPa to 80 hPa: the Φ function increases with P H 2 O , the accelerating effect being more pronounced than in the pre-transition stage (Fig. 4b). On the contrary, the hydrogen pressure has no effect on the oxidation rate [22]. Effect of Nb content. An alloy containing 0.4%wt in niobium (Zr0.4%NbO) has been used, in order to compare its oxidation behaviour to ZrNbO one. This alloy does not contain β-Nb precipitates but contains Zr(Nb, Fe) 2 phases [23]. The oxidation curves giving the oxide thickness X versus time, obtained for the two alloys at 520°C under 13 hPa in hydrogen and 67 hPa in water vapour are represented on Fig. 5a. It can be seen that the higher the Nb content, the higher the oxidation rate in pre-transition (which is in agreement with the literature data) [START_REF] Choo | [END_REF][24][25][26]. The post-transition stage was not reached, in the experiment, with Zr0.4%NbO alloy. Fig. 5b shows the rate curves dt dX versus X for both alloys. As a scaling factor allowing to superimpose the two curves could be found, the oxidation rates of these two alloys remain proportional during all the reaction time, with a ratio equal to 1.5 (which is lower than the ratio of the Nb concentration in the alloys (equal to 2.5)).

0 P P Φ Φ [7,

4%Nb (⎯) and ZrNbO (---).

Characterisation of the samples. The cross sectional views of the oxide scale grown during the pre-transition stage present a continuous and uniform layer adherent to the substrate. Fig. 6 shows micrographs of samples oxidised at 530°C in 13 hPa and 10 hPa of water vapour and hydrogen respectively. The thickness calculated from the weight gain is 3.5 µm for Fig. 6a and 6 µm for Fig. 6b. The metal/oxide interface is more or less undulated, and short cracks parallel to the interface appear regularly inside the layer [9]. A delayed oxidation can also be observed straight above some of these cracks.

Similar cracks are observed in samples oxidised after the kinetic transition, up to an oxide thickness of 12 μm and more. Moreover, new large cracks are observed which are perpendicular to the metal/oxide interface and connected to the gaseous atmosphere.

Discussion.

Pre-transition. Since the steady-state and the "ΦE test" have been validated experimentally, the assumption of a rate determining step can be made and the oxidation rate is given by (X is the oxide thickness):

) t ( E ) P , T ( S 2 n V dt dX i o o ox Φ = (2)
where V ox is the molar volume of the oxide (m 3 .mol -1 ), n 0 is the initial amount of metal (mol), S 0 is the surface of the samples (m 2 ). It has been observed that the Φ function increases with the partial pressure of water vapour. A simple mechanism involving oxygen vacancies can be written to describe zirconia growth [7], which provides a good interpretation of the oxidation of Zircaloy-4. From this mechanism, it is well-known that only a rate-limiting step located at the external interface (gas/oxide) can account for an accelerating effect of the partial pressure on the Φ function [7].

If such an interface step was rate-limiting, the oxidation rate should be constant with time and the kinetic curves should be linear, which is not the case. Some authors [10,[START_REF] Frichet | 12 th Symposium on Zirconium in the Nuclear Industry[END_REF] have interpreted the curves with a mixed reaction-diffusion rate, but this assumption is not in agreement with our experimental results [20,22].

Consequently, a new mechanism has to be proposed to describe zirconia growth on ZrNbO alloy. The dissociation of water vapour at the surface of the oxide probably leads to the formation of adsorbed hydroxyl groups (OH) and then, the diffusion of such species through the oxide layer may predominate. Considering the SEM micrographs of our samples (Fig. 6), the oxide layer appears to be homogeneous (excepted a few cracks). Nevertheless, it has been found by impedance spectroscopy [10,[START_REF] Garzarolli | Zirconium in the Nuclear Industry, 9th International Symposium[END_REF] that the oxide layer grown on ZrNbO could be divided into two sub-layers. For the kinetic modelling, two possibilities can be imagined, involving one or two oxide layers: -either one oxide layer with micropores and diffusion of adsorbed species via this porosity, -or two oxide layers: the microporous layer and a very thin and dense layer near the metal/oxide interface (for example, the native oxide layer), in which the diffusion of oxygen vacancies is very rapid. The thickness of this dense layer is assumed to remain very small during all the pre-transition stage.

In both cases, the weight gain is considered to be due only to the growth of the microporous oxide layer. The oxygen transport through this layer is supposed to occur via adsorbed hydroxyl groups (as seen in the mechanism below) ; thus, a part of the hydrogen produced by the oxidation is released near the metal/oxide interface. Adsorbed hydrogen can migrate through the micropores towards the external interface.

The various steps involved in the growth mechanism are detailed in Table II

O s H ' e V s OH O + - + → • • ← -+ - (7)
+ + • • → ← Table II: Oxide growth mechanism in pre-transition.
The hydrogen amount in the metal can be neglected in the total weight gain (this is in agreement with the very low amount of hydrogen measured by Cezus in the metal in the pretransition stage (about 30 ppm)).

In order to account for the decrease in the rate versus time in pre-transition, we have to consider a diffusion step as rate-limiting. The kinetic curves being not strictly parabolic, we have successfully fitted them using the following equation (as previously done in our study of Zircaloy-4 oxidation [7]):

( X k exp X k dt dX 2 1 - = ) (3)
in which X represents the thickness of the microporous layer.

The curves giving the rate dt dX versus X were fitted using several laws (parabolic, cubic, power…), the comparison of these fits are given in Fig. 7. Eq. ( 3) always leads to the best agreement with all our experiments. From Eq. ( 3), by comparison to Eq. ( 2), the Φ(T,P i ) and E(t) functions can be easily calculated [22]:

( ) . , i O D C T P X Φ = and ( ) ( ) X X k exp n X S C . D t E 2 0 0 0 - = (4)
X 0 is a characteristic length, for example the initial thickness of the sample. The variations with temperature and partial pressures of the Φ(T,P i ) function are deduced from the growth mechanism detailed above. Assuming that the rate-limiting step is the diffusion of adsorbed hydroxyl groups OH --s from the external interface into the porous layer via the surface of the pores, and neglecting the concentration of OH --s groups at the interface of arrival, the Φ function can be calculated :

( ) ( ) 4 H 4 H 4 H 1 3 2 1 O H O H 3 2 1 0 0 i K P K P K P K K 1 K K P P K K K X D P , T 2 2 2 2 2 + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + β + β θ = Φ (5)
in which K i is the equilibrium constant of the step (i), θ 0 is the concentration of adsorption sites per unit of surface area, and β is the electrons activity in the oxide layer, which is supposed to be fixed by the Nb content in solid solution in the oxide ( [ ]

' Zr e Nb • ⎡ = ⎣ ⎤ ⎦ =β).
The variations of Φ(T,P i ) with P H 2 O follow an hyperbolic law, which is in a good agreement with the experimental data, as shown in Fig. 4a (continuous line).

From Eq. ( 5), Φ(P) should vary with hydrogen pressure like , but no influence of P H 2 was observed experimentally in the range (2.5 -40hPa). However, this law exhibits a constant behaviour when P H2 increases, so experiments at lower hydrogen pressure would have been useful to validate it. Besides, Eq. (5) shows that if the value of the equilibrium constant K 4 is high compared to the values of P H2 , the term 4 H K P 2 may become very small and may be neglected besides the other terms. In that case, no effect of hydrogen pressure would be observed on the oxidation rate. Unfortunately, information on the values of the various equilibrium constants are not available in order to support one of these suggestions. Moreover, the Φ(T,P i ) function increases when the Nb content in solid solution in the oxide, β, increases, which is in agreement with our experimental results [22].

It has also been observed that platinum has an accelerating effect on the oxidation rate of ZrNbO [10,11], in water vapour or oxygen. As the platinum layer leads to an increase of the equilibrium constant K 3 [22], this model allows to propose an interpretation of the effect of platinum on the oxidation rate of ZrNbO.

Finally, we have shown that the pre-transition stage in ZrNbO oxidation cannot be accounted for by the mechanism usually proposed for the oxidation of zirconium alloys, which involves the diffusion of oxygen vacancies through the dense oxide layer as the rate controlling step. Thus, we have proposed another mechanism allowing to account for the accelerating effect of water vapour and Nb content. The important feature is to demonstrate that it is possible to explain the influence of the gases without assuming necessarily an interface rate-limiting step. In the case of ZrNbO, the oxidation is assumed to be controlled by a diffusion step, in which the concentration gradient of the diffusing species is fixed by a step located at the external interface (and not by a step located at the metal/oxide interface, as in the mechanism proposed usually).

Transition and post-transition. It is clear from the results of the jumps method (Fig. 3b and4) that the diffusion step controlling the beginning of the oxidation is no longer valid, even before the rate reaches its lowest value (since the ratios of the rates obtained with the jumps method decreases from about 7 µm, whereas the oxidation rate is minimum at about 10 µm, see Fig. 3b). Large cracks connecting a part of the oxide layer to the gaseous atmosphere appear from this point, which could be related to the increase in the rate after the transition.

Then, when the oxide thickness exceeds 17 µm, the shape of the kinetic curve is linear. Moreover, the ratio of the rates remains constant (with a value different from the pretransition one, see Fig. 3b). Consequently, it can be deduced that the post-transition stage is controlled by an interface step, largely influenced by the water vapour pressure. Another possibility could be the diffusion through a layer with a constant thickness, but this is less probable due to the influence of the water vapour pressure.

The pre-transition mechanism could explain the post-transition kinetics if the rate of one of the interface steps of this mechanism became far lower in post-transition than it was in pretransition [22]. Another way to explain the post-transition kinetics would be a change in the mechanism, involving completely different steps. Unfortunately, we have no experimental evidence to support one of these suggestions. More investigations would be necessary to give a valuable interpretation of the post-transition kinetics.

Conclusions

The oxidation of a ZrNbO alloy in a mixture of hydrogen and water vapour at 530°C exhibits clear differences between the pre-and post-transition stages.

In pre-transition, the oxidation proceeds in a steady state, and is influenced by water vapour and Nb content, whereas hydrogen pressure has no effect in the studied range. The rate decreases with time, according to a sub-parabolic law.

A mechanism has been proposed to account for these results, which involves the diffusion of adsorbed OH groups in a microporous part of the oxide layer, as rate-determining step.

Then, after a transition period beginning at about 7 µm at 530°C, a rate limiting step can again be considered when the oxide thickness is higher than 17 µm. As the curves are linear, the oxidation is probably controlled by an interface step. This step has not been determined, but it is probably different from the steps proposed in the pre-transition mechanism.

Consequently, the transition is accompanied by a change in the oxidation mechanism, that could be linked to the change in the morphology of the oxide layer.

  /S)/dt (mg/cm 2 /s) time (h) Steady state approximation. In Fig. 2, the variations of the rate of weight gain ( dt dm ) and the heat flow dt dQ versus time are represented in the pre-transition stage, at 550°C. A

Fig. 2 :

 2 Fig. 2 : Rate of weight gain (---dm/dt) and heat flow (⎯) versus time for ZrNbO at 550°C in water vapour(13 hPa) and hydrogen(10 hPa), for the pre-transition stage.

Fig. 3 :

 3 Fig.3: Rate of weight gain (530°C, P H2O = 13 hPa, P H2 = 10 hPa) and ratios of the rate before and after the temperature jumps, in the pre-(a) and post-transition (b) stages.

Fig. 4 :

 4 Fig. 4 : Variations of the reactivity of growth Φ with the partial pressure of water vapour, (a)in the pre-transition stage (jumps from P H2O = 13 hPa to P, for X = 2.5μm, T = 500°C, P H2 = 10 hPa), (b) in post-transition (jumps from P H2O = 13 hPa to P, for X = 25μm, T = 530°C, P H2 = 10 hPa).

Fig. 5 :

 5 (a) Oxide thickness versus time for ZrNbO and Zr0.4%Nb (T = 520°C, P H2O = 13hPa, P H2 = 10 hPa), and (b) oxidation rate (dX/dt) versus time for Zr0.

Fig. 6 :

 6 Cross-sectional views of oxide layers grown on ZrNbO at 530°C (P H2O = 13 hPa, P H2 = 10 hPa) : X ≈ 3.5 μm (a) and X ≈ 6 μm (b).

Fig. 7 :

 7 Fig.7 : Rate of oxidation of ZrNbO as a function of the oxide thickness before the kinetic transition -comparison with various rate laws.

  .

	Interfaces External interface	Elementary steps (1) adsorption of water on adsorption sites s: (2) dissociation of water into OH groups: (3) reduction of hydroxyl groups: (4) desorption of hydrogen:	s H + -s O -2 s -OH H OH + → s ← s → ← ' e -+ -→ s 2 H s 2 + ← -→ O s 2 + s -H 2 H -OH O H 2 ←	s
		(5) diffusion of adsorbed OH -groups in the
		porous layer.
		(6) interface step :

O

Acknowledgements

This work has been carried out with the financial support of FRAMATOME-ANP. We would like thank P. Barberis (Cezus) and A. Frichet (Framatome-ANP) for helpful discussions.


Julien Chiquet

Pierre Gutierrez

Guillem Rigaill

Fast tree inference with weighted fusion penalties

Given a data set with many features observed in a large number of conditions, it is desirable to fuse and aggregate conditions which are similar to ease the interpretation and extract the main characteristics of the data. This paper presents a multidimensional fusion penalty framework to address this question when the number of conditions is large. If the fusion penalty is encoded by an q -norm, we prove for uniform weights that the path of solutions is a tree which is suitable for interpretability. For the 1 and ∞-norms, the path is piecewise linear and we derive a homotopy algorithm to recover exactly the whole tree structure. For weighted 1-fusion penalties, we demonstrate that distance-decreasing weights lead to balanced tree structures. For a subclass of these weights that we call "exponentially adaptive", we derive an O(n log(n)) homotopy algorithm and we prove an asymptotic oracle property. This guarantees that we recover the underlying structure of the data efficiently both from a statistical and a computational point of view. We provide a fast implementation of the homotopy algorithm for the single feature case, as well as an efficient embedded cross-validation procedure that takes advantage of the tree structure of the path of solutions. Our proposal outperforms its competing procedures on simulations both in terms of timings and prediction accuracy. As an example we consider phenotypic data: given one or several traits, we reconstruct a balanced tree structure and assess its agreement with the known taxonomy.

Introduction

As data floods in, it is now possible to compare many features across a very large number of conditions in various fields of science. To cite but a few, we encounter this setting in genomics where high-throughput technologies allow us to monitor the expression level of many genes (the features) at various stages of a given biological process (the conditions); this also occurs in phylogenetics where several quantitative traits (the features) are available for many species (the conditions). Beyond biological sciences, sets of data gathered in astronomy are now routinely composed of millions of conditions for hundreds of features. An interesting question is to group together -or fuse -these conditions across the feature space, arguing that these conditions should not really be considered as different. In other words, we aim at recovering an interpretable clustering of those conditions.

There are basically two cases: either a prior group structure between the conditions is known, or it is not. In the first case, one typically applies one-way ANOVAor MANOVA for multiple features -to test for any significant difference between all pairs of groups. The final structure between the conditions then depends on the level of significance used to test for differences. However, when the number K of groups is large, which typically occurs for a large number n of conditions, this leads to a multiple-testing issue and algorithmic problems since the number of pairwise tests is in O(K 2 ). Furthermore, each test is performed independently and the resulting structure is not necessarily simple and easily interpretable.

In the second case, when no prior group structure is available, we basically face a clustering problem over the multidimensional space of the features. A popular heuristic to solve this problem is agglomerative clustering, which defines a hierarchical structure between the conditions. Hierarchies are very appealing for interpretability. A serious bottleneck of agglomerative clustering when analyzing large data sets is its complexity in O(n 3 ), which can be reduced to O(n 2 ) using single-linkage clustering.

There are two major issues for large values of n: i) the need for an interpretable structure between the conditions and ii) the need for a computationally and statistically efficient estimation procedure. These two goals cannot be reached simultaneously, neither by MANOVA nor by agglomerative clustering algorithms, due to restrictions either on the interpretability of the inferred structure or the computational burden of the procedure. This paper presents a unifying approach to tackle these two problems simultaneously by means of a weighted fusion penalty that constructs a hierarchical structure on the conditions at a low computational cost and reaching the two aforementioned goals. Section 2 presents our proposal in detail and puts it in perspective with existing methods. Then we use the optimality conditions detailed in Section 3 to characterize the regularization path (Section 4). In Section 5, we propose weights for which the path is provably without splits. For the 1 -norm some of those weights lead to a desirable balanced tree structure. In Section 6 we present a homotopy algorithm which is in O(K log K) for well chosen weights. We also provide an efficient embedded cross-validation procedure to tune up the level of aggregation -or fusion -between groups in the ANOVA settings. Numerical experiments illustrate the extremely competitive performance of our algorithm in terms of timings. Section 7 presents consistency results that bring statistical guarantees for our approach. We illustrate our theorem on a simulation study that shows that our weights are more efficient than those of its competing procedures. Finally, Section 8 is dedicated to a complete example in phylogenetics where our method is applied to the reconstruction of a balanced tree structure across several phylogenetic features between many species. We assess its relevance by comparison with the known phylogeny.

A penalized framework for tree inference

To bring MANOVA and hierarchical clustering together in the same unifying penalized framework, note that the latter can be considered as a particular case of the former when there is only one condition per group, i.e, when K = n. This can be thought of as a non-informative prior on the clustering between the conditions.

To be more specific, we set y ij the observation of a continuous random variable that describes the intensity of the jth feature in condition i, with i ∈ {1, . . . , n} and j ∈ {1, . . . , p}. The p-dimensional vector y i = (y i1 , . . . , y ip ) encompasses the data related to condition i across the p features. We are given a partition with K groups as prior knowledge that is depicted by the indexing function κ : {1, . . . , n} → {1, . . . , K}. In words, κ indicates the group to which condition i is allocated a priori. The number of elements in group k is denoted by n k = card {i : κ(i) = k}, such that k n k = n. One-way MANOVA is a multivariate linear regression problem whose parameters are fitted by minimizing the residual sum of squares, i.e., , where β kj is the coefficient for the jth feature in the kth group, such that β k = (β k1 , . . . , β kp ) ∈ R p . The final structure between the conditions is obtained by testing for significant differences between all pairs of estimated means ( βkj , β j ) using Fisher statistics.

Compared to MANOVA, hierarchical clustering assumes one individual per group, that is K = n or equivalently κ(i) = i for all i = 1, . . . , n. It performs agglomeration by recursively joining the closest points. As suggested by [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF], hierarchical clustering aims at solving the following optimization problem:

minimize β∈R np n i=1 y i -β i 2 2 , s.t. i>i 1 β i =β i ≤ t. (1)
The complete hierarchy between the conditions is recovered by starting from t = n(n -1)/2, where no constraint applies, then by decreasing t until all points agglomerate. This immediately suggests a corresponding scheme for agglomerating groups of conditions in MANOVA just by using the prior grouping knowledge encoded by κ in the square loss. However, Problem (1) and its MANOVA counterpart are difficult combinatorial problems in general. To overcome this restriction, we consider the following convexified Lagrangian formulation which includes the whole family of optimization problems discussed throughout this paper:

minimize β∈R Kp 1 2 n i=1 y i -β κ(i) 2 2 + λ k, :k = w k Ω(β k -β ). (2)
In general, Ω is a norm and w k are positive, symmetric weights over all pairs of groups in {1, . . . , K} such that w k > 0 and w k = w k . The penalty term and the choice of Ω is designed to encourage elements of β to "fuse" by enforcing similarity between pairs of vectors (β k , β ) as in the fused-Lasso signal approximator [START_REF] Friedman | Pathwise coordinate optimization[END_REF], which is an 1 -based method designed to aggregate pairs of elements. As such, we refer to the penalty term in (2) as a "fusion" penalty. In the multidimensional case though, other choices are possible for Ω that induce a fusion effect. The level of fusion is tuned by two parameters: the global level of penalty λ and the group specific weights w k , the choice of which is of the highest importance. It conditions both i) the ability of the method to infer an interpretable structure between the conditions, ii) the existence of fast algorithms to fit the parameters β for various values of λ and iii) the existence of statistical guarantees for the estimator. The main objective of this paper is to study classes of weights that reach these three goals simultaneously.

Links to existing works. Problem (2) is a generalization of two interesting existing procedures related to ours. The first one arose in the clustering framework and is known as the Clusterpath [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]. The Clusterpath covers cases in (2) where K = n and Ω(•) = • q with q ∈ {1, 2, ∞}. Still, for general weights, the complexity of the associated algorithms does not improve over the agglomerative clustering, and the inferred structure is not a tree. However, when q = 1, the path of solutions is linear with respect to λ and a homotopy algorithm is used by Hocking et al. to recover the solutions over all the values of λ that either correspond to events of fusion or split between a couple (β k , β ). Moreover, if w k = 1 and q = 1, they showed that no split event can occur and that a homotopy algorithm can be implemented in O(n log(n)). In other words, the reconstructed structure is a tree in this case. However the unitary weights typically lead to unbalanced hierarchies which are not fully satisfactory.

A second close cousin to our approach is the Cas-ANOVA of Bondell and Reich (2008b). Cas-ANOVA is a 1 -penalized version of the ANOVA which corresponds to (2) in the univariate setting where p = 1 and Ω(•) = • 1 . The main contribution of this proposal is statistical: Bondell and Reich introduce adaptive weights w kl ∝ √ n k + n /(ȳ k -ȳ ), where n k is the number of conditions in group k and ȳk = i:κ(i)=k y i /n k is the corresponding empirical mean. Similar weights have been proposed in [START_REF] Gertheiss | Sparse modeling of categorial explanatory variables[END_REF] to cope with ordered categorical variables. These weights have an adaptive property such that the corresponding estimator of β enjoys asymptotic consistency, in the manner of the adaptive Lasso (Zou, 2006). Still, Cas-ANOVA weights do not lead to a tree when the number of individuals per condition is unbalanced, i.e., n k = n for any couple (k, ). Moreover, the optimization procedure is in O(K 2 ) and only provides the solution for a given λ. We also experienced numerical instability using Cas-ANOVA weights.

Contributions. Compared to these two works, our contributions are the following:

• We prove that no split can occur along the path of solutions in (2) when w k = n k • n and Ω(•) is an q -norm. As a consequence, this proves that the Clusterpath does not split for unitary weights, whatever the choice of the norm (as conjectured by Hocking et al. for the 2 -norm).

• When Problem (2) is separable across the features (e.g., when Ω is the 1norm), we introduce distance-decreasing weights for which we prove that the path is a tree. From an interpretation point of view, this family of weights is particularly interesting as it leads to balanced tree structures.

• For the 1 -norm, we introduce exponentially adaptive weights that enter the family of distance-decreasing weights. They enjoy asymptotic oracle properties that guarantee selection of the true underlying structure for a large scale of possible λ. This shows that our estimator shares the same asymptotic properties as Cas-ANOVA, but for a larger range of λ and at a much lower computational cost.

• We provide a general homotopy algorithm for (2) when Ω(•) = • 1 , whatever the choice of w k . On a single feature, the initialization for unspecified weights is in O(K 2 ) and the homotopy itself is in O(K log(K)). However, we propose a faster initialization procedure for exponentially adaptive weights such that the whole complexity for p features is in O(pK log(K)) -or O(pn log(n)) in the clustering framework.

• When the number K of prior groups is smaller than n (e.g., in the ANOVA settings, when there are some replicates per condition/group), a natural crossvalidation (CV) error can be defined. In this case, we develop a fast procedure that takes advantage of the DAG (directed acyclic graph) structure of the path of solutions along λ. This approach has a lower complexity than a standard CV procedure.

In short, we propose choices for weights in (2) that induce a balanced tree structure between the conditions such that the associated estimation procedure enjoys the good computational properties of the 1 -Clusterpath with unitary weights, with stronger statistical guarantees than Cas-ANOVA.

Motivating example in phylogeny. As a simple motivating example, we consider a univariate problem in phylogeny. We want to reconstruct a tree between many species based on some simple features (like the height, or the weight of individuals). Ideally, this tree should resemble the known phylogeny. We illustrate this task on the "Animal Ageing Longevity Database"1 , which provides various features for many animal species. Here, we consider classifying bird species based on their birth weight. The known phylogeny groups these n = 184 individuals into 40 bird families, themselves grouped into 15 orders. We reconstruct the tree based on the weights and check whether it matches the orders and the family classification. Recovered solution paths of (2) are plotted in Figure 1 for a) the Cas-ANOVA weights (Bondell and Reich, 2008b) ; b) the "default" Clusterpath weights [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]; and c) our own weights that we call "fused-ANOVA" weights. On the left panel, the Cas-ANOVA path includes many splits which make interpretation rather difficult. On the middle panel, default Clusterpath weights, as expected, provide a tree structure. Still, the structure of this tree is unbalanced and thus not fully satisfactory in the sense that small groups often fuse with very large ones. Specifically, the Clusterpath tree does not capture the simple fact that there are visibly three groups corresponding to light, medium or heavy birds. Conversely, the fused-ANOVA tree in the right panel is more balanced and clearly exhibits these three groups. Furthermore, it is in better agreement with the known phylogenetic classification, improving the rand index by 5% compared to ClusterPath. estimated coefficients βλ Multidimensional 1 Clusterpath and fused-ANOVA. In the previous example, we consider only one feature. In practice, one often has to consider multiple features at the same time. This is possible with our proposed weighted 1penalty. Indeed, as noted by [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF], Problem (2) is separable on dimensions when considering the 1 -penalty, which is also the case for our weighted fused-ANOVA scheme. Thus, Clusterpath and fused-ANOVA algorithms solve the multidimensional problem in two steps:

1. First, they recover p independent trees (one per dimension). This task can be easily executed in parallel.

2. Second they aggregate those p trees in a consensus tree. This is done by considering the same penalty value (λ) corresponding to a given height in those trees. Two individuals k and are in the same multidimensional cluster if they have been fused on every dimension.

This multidimensional classification is recovered on a grid of λ in the clusterpath package and in the fusedanova package.

Note however that the classification recovered over all the dimensions is not necessarily better than those recovered on single, well-chosen features. We illustrate this point at the end in Section 8 on phylogenetic data: in a number of cases, the best agreement with the known phylogeny is obtained by a single-feature-based tree.

Optimality conditions and consequences

We start by characterizing Problem (2), giving elementary facts which are at the basis of most of our results. Note that the objective function in (2) is a nonsmooth function which is strictly convex in β and thus admits a unique solution when λ ≥ 0. This solution can be characterized by the KKT (Karush-Kuhn-Tucker) conditions that may be derived thanks to subgradient calculus (see, e.g., [START_REF] Boyd | Convex optimization[END_REF]). In the case at hand, β is optimal if, for all k ∈ {1, . . . , K}, β k verifies the following subgradient equations:

0 p = -n k (ȳ k -β k ) + λ : =k β k =β w k τ k + λ : =k β k =β w k ∂Ω(β k -β ) ∂β k , (3)
where ȳk = i:κ(i)=k y i /n k is the vector of empirical means for the kth group across every feature. The p-dimensional vectors τ k are such that, for any k, there exists = k with β k = β such that τ k = -τ k and Ω(τ k ) ≤ 1. We omit the proof as it is a straightforward adaptation of the fused-Lasso subgradient equations [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] to the multidimensional case, with a general norm Ω.

Interesting consequences arise when summing the subgradient equations (3) for all β k which are "fused" in the same cluster, as stated in the following Lemma.

Lemma 1. Consider a cluster C = {k : β k = β C } formed by some β k , where β is the solution to (2). Then we have

β C = ȳC - λ n C / ∈C w C ∂Ω(β C -β ) ∂β C , (4)
where n C = k∈C n k , ȳC = k∈C ȳk /n C and w C = k∈C w k .

Proof. By summing (3) for all k ∈ C, we have

0 p = -n C ȳC + n C β C + λ k, ∈C:k = w k τ k + λ k∈C, / ∈C w k ∂Ω(β k -β ) ∂β k .
Then, by the KKT conditions, we must have τ k = -τ k for some k, ∈ C. Thus the third term on the left-hand side of the above expression vanishes by symmetry of the weights w k . Also notice that ∂Ω(

β k -β )/∂β k = ∂Ω(β k -β )/∂β k for any k, k ∈ C, /
∈ C, and we easily get the desired result.

Regularization path and tree structure

Characterization of the minimization Problem (2) in terms of its optimality conditions is essential in many ways. In particular, Lemma 1 allows us to characterize the regularization path of solutions {β(λ), λ > 0} depending on the choices of Ω and w k . This is important for our problem since the shape of the path is actually the structure recovered between the conditions. This is also important since it may induce some computational properties that guarantee a low complexity of the associated fitting procedure. This section investigates which conditions must be imposed on the regularization path to ensure a structure that is fully satisfactory both in terms of algorithmic complexity and interpretability, namely, a balanced tree structure.

The mildest condition which is required is continuity of the regularization path, that is to say, of the function {β(λ), λ > 0}: without continuity, interpretability of the recovered structure is obviously out of reach. This property is straightforward for solutions of problems of form (2) which is strictly convex. However, continuity of the path is not enough to provide an interpretable structure, and we shall investigate conditions ensuring that the inferred structure is a tree. In terms of regularization path, it requires that any couple of parameters which have fused at a certain time λ 0 such that β k (λ 0 ) = β (λ 0 ) = β C cannot "split" anymore in the future, that is, for any value λ > λ 0 that would correspond to a higher level in the hierarchy of the tree. Insights on this remark can be found in Figure 2, where various regularization paths are plotted in the univariate case. Paths on the top and bottom left panels contain splits, while the remainders do not and lead to trees with different shapes the properties of which are discussed later in this section.

Piecewise quadratic path

Piecewise linear path not a tree (splits) unbalanced tree balanced tree

Figure 2: Various typologies of the regularization paths in the single feature case that lead to more or less interpretable structures.

Though highly desirable, guaranteeing a tree is complicated as the absence of splits in the path of (2) depends jointly on the choice of the weights w k and on the fusing norm Ω(•). In the following Theorem, we provide a simple generic choice for the weights that ensures the absence of splits in the general formulation with q -norms.

Theorem 1. If Ω is an q -norm with q ∈ {1, . . . , ∞} and w kl = n k • n , the path of solutions {β(λ) : λ > 0} of (2) contains no splits.

The proof is postponed to Appendix A.1. Schematically, it investigates the subgradient equations of ( 2) and shows that given a solution at λ 0 , we can always explicitly construct for λ > λ 0 a valid subgradient not involving any split. Theorem 1 generalizes the results of [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF] obtained for Ω(•) = • 1 in the clustering case when w k = n k = n = 1 to any q -norm Ω.

A consequence of Theorem 1 is that the existing implementation of the 2 Clusterpath -or any other q solver -can be simplified by no longer considering the eventuality of splits with default weights (see Algorithm 1 in [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF].

As said before, guaranteeing a tree-structure is the first step towards interpretability. As such, Theorem 1 characterizes an interesting family of problems. Still, the scope of arbitrary norms with uniform weights is not fully satisfactory because, even when the structure is a tree,

• the path is not a linear function of λ in general, as illustrated on the first row of Figure 2. In this situation, detecting the events of fusion may be expensive.

It might be impossible to provide an efficient algorithm to infer the structure at a low computational cost.

• the inferred structure may be highly unbalanced. By unbalanced, we mean a tree where two parameters initially close to one another at λ = 0 fuse relatively late in the path of solutions. Such situations are depicted on the second column of Figure 2. It is obvious that disequilibrium may significantly narrow the potential for interpretability of the tree.

First, equilibrium of the inferred structure is a property that is mainly controlled by the w k . We cannot limit ourselves to w k = n k • n and must exhibit weights sharing both the equilibrium and the non-split property. This will lead us to the distance-decreasing weights described in the next section.

Second, piecewise-linearity -and thus existence of a fast path-following algorithm -is a property of the norm Ω. A solution path which is piecewise linear can be computed efficiently (and exactly) with a homotopy algorithm like the LARS for the LASSO [START_REF] Efron | Least angle regression[END_REF]. More generally, [START_REF] Rosset | Piecewise linear regularized solution paths[END_REF] give conditions for the existence of such a property in a broad penalized framework. These results are easily adapted to the case at hand, where we roughly have to differentiate Expression (4) over λ to conclude: for Ω(•) = • q any q-norm with q ≥ 1, then

∂β C ∂λ = 1 n C / ∈C w C signs(β -β C ) • |β -β C | q-1 β -β C q-1 q , ( 5
)
where | • | and signs(•) apply element-wise and • is the element-wise product. Application of Proposition 1 of Rosset and Zhu to these expressions implies that the path is piecewise linear only for q ∈ {1, ∞}. In other words, there must exist a homotopy algorithm to infer the structure between the conditions for the 1 and ∞ -norms.

More generally, we could use any norm Ω that builds on 1 and ∞ such as the OS-CAR (Bondell and Reich, 2008a). Note, however, that there is no guarantee that the number of steps will be small in the homotopy algorithm for general weights. In fact, [START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF] exhibits pathological cases for the LARS algorithm where the number of kinks in the piecewise linear path of solutions grows exponentially with the number of variables. Such cases can be transposed to the weighted fusion penalty with Ω(•) = • 1 , which corresponds to situations where there is a large number of splits along the path. To overcome this restriction and guarantee that the number of iterations required to fit the whole path of solutions will be small, we introduce in the next section a family of weights that ensures no split along the path of solutions for the particular case of the 1 -norm.

5 Distance-decreasing weights guaranteeing no split

In this section, we focus on the 1 -norm and generalize Theorem 1 to a larger class of weights that we call distance-decreasing weights, defined in Theorem 2. Indeed, although uniform weights ensure the absence of split, the recovered tree structure is often unbalanced. Intuitively, distance-decreasing weights should ensure that close neighbors fuse quickly. Here, we demonstrate that for such weights there is no split. Thus, the algorithm proposed by [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] for the generalized fused-Lasso is considerably simplified since there is no need to check for possible split events, and thus there is no need to solve potentially numerically unstable maximum flow problems.

Remark 1. Note that the absence of splits does not ensure a fast algorithm. Indeed, the initialization of the generalized fused-Lasso algorithm is for most weights in K 2 . We exhibit in Section 6 a subset of distance-decreasing weights for which initialization is linear and for which we can guarantee good statistical properties in Section 7.

Another advantage of the 1 -norm is that it brings separability across the p features in (2), that is to say, that the p-dimensional problem splits into p univariate problems. To recover a consensus classification, we first infer p independent trees (one per dimension) and then aggregate those p trees by considering the same penalty value λ. Thus, without loss of generality, we restrict the discussion to the following 1 univariate problem which is a weighted generalized fused-Lasso problem:

minimize β∈R K 1 2 K k=1 n k (ȳ k -β k ) 2 + λ k, :k = w k |β k -β |. (6)
For this problem, we get the following result:

Theorem 2. The path of solutions does not contain splits when weights are chosen such that

w k = n k n f (|ȳ k -ȳ |),
where f (•) is a decreasing positive function.

Schematically, the proof is based on two ingredients:

1. first, using geometrical arguments, it is possible to show that absence of splits is equivalent to preservation of the order along the path, that is to say, ȳk ≤ ȳ ⇔ βk (λ) ≤ β (λ); 2. second, by considering a problem that is dual to (6) as in [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] for the generalized Lasso, we show that distance-decreasing weights preserve the order.

The proof is detailed in Appendix A.2.

6 Fast homotopy algorithm for 1 weighted penalties

In this section, we consider algorithmic issues when Ω is the 1 -norm. As in Section 5, we restrict the discussion to univariate Problem (6) and thus give the numerical complexity in the case p = 1. For a p-dimensional problem, we aggregate the p univariate trees by considering the same values of λ for all trees.

An algorithm for general weights and its limitations. Optimization problem (6) can be solved for general weights w k by the homotopy algorithm proposed in [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] for the generalized fused-Lasso. This is also the solution retained in the clustering framework by [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF]. A schematic view of this algorithm adapted to ( 6) is depicted in 1.

Algorithm 1: Homotopy algorithm for the generalized fused-Lasso

Input: data, weights and initial groups {y i , w k , κ} Initialization for λ = 0 Initialize β k parameters (equal to the empirical means ȳk ) Initialize the list of possible next events (only fusion at this stage) while all groups are not fused do Find the next event (having the smallest λ), it can be a split or a fusion Update β k parameters accordingly Update the list of possible next events (fusion and split) end Output: Directed acyclic graph (DAG) of fusion and split events and associated values of the parameters This procedure for general weights has two major flaws that may have detrimental effects on its computational performance:

• By piecewise-linearity of the solution path, the total number of iterations (that is, the total number of events before all the groups have fused) is bounded. However, by rewriting (6) as a Lasso problem -which only requires straightforward algebra -we may construct pathological cases where there are (3 K + 1)/2 linear segments in the path of solutions (see [START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF], a complexity that we cannot afford even for a moderate number of conditions K.

• While detecting fusion events in Algorithm 1 may be cheap since it roughly only requires calculation of the slopes ∂β k (λ)/∂λ, checking for the possibility of split events boils down to maximum-flow problems the resolution of which at large scale may clearly be a bottleneck (see [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF].

To circumvent these limitations, we shall consider weights that prevent split events. Although the choice w k = n k n has been shown to prevent splits in Theorem 1, it will typically lead to fusion events occurring very late (that is, for large λ), even between groups having close empirical means. This corresponds to an unbalanced tree structure between the conditions, which is hardly interpretable. On the contrary, using the family of distance-decreasing weights, introduced in Section 5, prevents split events and leads to a balanced tree structure. In this case the total number of events is exactly K -1, which is the number of iterations required to fuse K groups into 1, assuming that there cannot be a fusion of more than two groups at once. As for the maximum-flow problems, they are completely eluded from the algorithm with these weights. Still, we have to take into account the cost of detecting successive fusion events and of updating the coefficients β k (λ) along the K -1 steps. In the next paragraph, we propose a solution inducing a global complexity of O(K log(K)) for a given choice of weights belonging to the family of distance-decreasing weights.

Weights with an O(K log(K)) implementation. First we need to define the next time a fusion event is going to happen. We proceed mainly as in [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] for the one-dimensional fused-Lasso signal approximator, except that the initial ordering is not defined by the neighborhood between the coefficients, but by the ordering of the empirical means ȳk . And thanks to the property of the distancedecreasing weights, this ordering remains the same throughout the algorithm, which allows us to compute the path in O(K log K) operations. Here are some details.

At the initialization step, one has λ 0 = 0, and the next time a fusion occurs is

t(λ) = arg min t k (λ)>λ0 t k , t k (λ) = λ 0 -(β k (λ 0 )-β (λ 0 )) ∂β k ∂λ (λ 0 ) - ∂β ∂λ (λ 0 ) -1 . (7)
In words, it is the smallest value of λ among all the values such that two coefficients fuse. The main cost in ( 7) is due to the calculation of the slopes ∂β k /∂λ at λ 0 = 0. Note that β k (0) = ȳk , and by Lemma 1 and ( 5), one has

∂β k ∂λ (0) = - 1 n k =k w k signs(ȳ k -ȳ ). ( 8
)
For general weights w k , computing these slopes for all k requires O(K 2 ) operations and is the limiting factor of the algorithm. However, we provide a O(K log(K)) procedure for a special case of our distance-decreasing weights that we call "exponentially adaptive weights" because of their statistical properties (see Section 7). They are defined by

w k = n k n exp{-α √ n|ȳ k -ȳ |}, α > 0, (9)
for α a positive constant. The key idea to achieve O(K log(K)) complexity with these weights is that each slope can be computed as the sum of two terms, for which there exists a simple recurrence formula: first, we order the ȳk in decreasing order, which can be done in O(K log(K)) operations. Assuming this is done, we obtain

∂β k ∂λ (0) = - =k n signs(ȳ k -ȳ ) exp -α √ n|ȳ k -ȳ | = <k n exp -α √ n(ȳ -ȳk ) - >k n exp -α √ n(ȳ k -ȳ ) = exp α √ nȳ k <k n exp -α √ nȳ L k -exp -α √ nȳ k >k n exp α √ nȳ R k . The recurrence formulae are R k+1 = R k +n k exp {-αȳ k } and L k-1 = L k +n k exp {αȳ k }.
By this means, the initial slopes (8) and thus the first fusion time can be computed in O(K log(K)).

Then, for each of the K -1 steps of the algorithm, we only need to update the two slopes and the two coefficients which are currently fusing. This only requires a constant number of operations. Concerning the next fusion time, however, the new minimum among the updated t k (λ + 0 ) is found in log(K) if stored in an appropriate structure. This way we can reach O(K log(K)) for the global complexity.

As a final remark, note that we use the same storage solution -namely a binary tree -as did [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] for the one-dimensional fused-Lasso. By this means, we maintain the memory requirement at a low level that only grows linearly in K.

An embedded cross-validation procedure. Providing the whole path of solutions is clearly interesting for interpretability, since we force it to be a tree by means of an appropriate weighting scheme coupled with the 1 -norm for fusion. Still, it is always necessary to provide a practical way to choose the tuning parameter, which corresponds in the case at hand to choosing the level at which to cut the tree. This also gives a fixed classification between the initial conditions.

When the number K of prior groups is smaller than n (e.g., in the ANOVA settings), a natural cross-validation (CV) error can be defined. Although CV is often incriminated for being time-consuming, it is possible in this case to rely on the tree structure of the solution -or DAG in the case where split is allowed in the algorithm -to enhance the performance. Indeed, we can first build a tree using a training set (in which all prior groups are present) and then assess its performance by measuring its ability to predict the remaining individuals of the test set for any given value of λ. Here, we perform the CV on a predefined grid of L values of λ because the fusion times will be different for every new training set and it would be memory intensive to store the CV-error for every of those fusion time.

To be more specific, we consider a split of the data in a train set D and a test set T such that each prior group is represented in the train set. Using D, we recover a fused-ANOVA tree and an estimator βD κ(i) (λ). The test error on T is

CV err (D, T , λ) = i∈T y i -βD κ(i) (λ) 2 . (10)
A naive approach to computing (10) is to consider each prior group at a time on a grid of λ. Computing the prediction based on a given fitted regularization path requires O(log(K)) operations to search through the tree of solutions. This has to be done for the K prior groups and for the L values in the grid of λ. Hence, computing the sum (10) naively has a total complexity of O(LK log(K)) (which dominates the complexity in O(K log(K)) of the fit itself!).

On the contrary, our embedded cross-validation procedure takes advantage of the tree structure of the fit in the computations whenever possible. Indeed, along the branch of a cluster C, the estimator βD κ(i) (λ) is a piecewise linear function of λ and thus the error (10) is a piecewise quadratic function of λ. The coefficients of this quadratic function are easily updated when constructing the tree, and the error along this branch is computed in O(1) for any λ rather than O(|C| log(K)). More precisely the error in (10) of cluster C decomposes thanks to the Huygens formula as

i∈T :κ(i)∈C y i -βD κ(i) (λ) 2 = i∈T :κ(i)∈C (y i -ȳT C ) 2 + n T C ȳT C -βD κ(i) (λ) 2 ,
where n T C = card({i ∈ T : κ(i) ∈ C}) and ȳT C is the empirical mean of individuals of cluster C, i.e.,

ȳT C = 1 n T C i∈T :κ(i)∈C y i .
It is difficult to assess exactly the gain brought by using the tree structure for computing the CV error in general. Indeed, it depends on the tree itself, the length of its branches, its height and so on. Assuming a binary balanced tree of height log(K), with branches of equal length and an equally spaced grid of λ, we can show that the complexity is in O(LK/ log(K)). If some groups fused rapidly (as with the fused-ANOVA weights), the gain could be even greater. In practice (see Figure 3.c), we often see a ten-fold difference between our CV procedure and a naive implementation.

Timings. We implemented both the general and the without-split version of Algorithm 1 in C++ embedded in an R-package called fusedanova distributed on R-forge. It contains a wide family of weights which are not mentioned in this paper due to space requirements. Figure 3 illustrates the rather good performance of our algorithm and implementation through three numerical experiments: a) In the left panel, we illustrate the capability of our method to treat large scale problems extremely fast: we generate a size-n vector y such that y i ∼ N (0, 1) and assume n = K, meaning one condition per group2 . We vary n from 10 2 to 10 8 and record the corresponding timing in seconds. We apply our method with the exponentially adaptive weights and average over 10 trials. As can be seen, we can reconstruct a tree on n = 10 6 observations in about 10 seconds.

b) The middle panel illustrates the gain in runtime due to the fact that we no longer have to check for splits in the homotopy algorithm using a maximum-flow solver.

We generate data as in the preceding experiment but with K conditions each containing n k = 20 replicates. When K = 10 3 , the gain in seconds brought by not checking for the possibility of splits is of almost 2 orders of magnitude.

c) The right panel illustrates the performance of our embedded CV procedure compared to the naive implementation. We used the same settings as in the previous experiment.

timings in seconds (log) q q q q q q q 10 -2 10 -1 10 0 10 1 10 2 10 3 10 2 10 3 10 4 10 5 10 6 10 7 10 8 q q q q q q q q q q q q q q q q q q q q q q q 10 -2 10 -1.5 10 -1 10 -0.5 10 0 10 0.5 10 1 10 1 10 2 10 3 algorithm q q general no split q q q q q q q q q q q q q q 10 -2 10 -1.5 10 -1 10 -0.5 10 0 10 0.5 10 1 10 1 10 1.5 10 2

cross-validation

q q embedded classical number of conditions K (log) We tried other implementations to solve (6) such as the Clusterpath package by [START_REF] Hocking | Clusterpath: an algorithm for clustering using convex fusion penalties[END_REF], the flsa package by [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF] or the genlasso package by [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF]. These implementations do not fully exploit the structure of the problem and have runtimes considerably longer than ours, even for moderate K. Thus, we do not report their timings here.

Statistical guarantees

Asymptotic settings. To discuss the asymptotic properties of our exponentially adaptive weights (9), we shall consider the following univariate3 ANOVA model

y i = β κ(i) + ε i , s.t. E(ε i ) = 0, Var(ε i ) = σ 2 , i = 1, . . . , n, (11)
where β = (β 1 , . . . , β K ) is the true vector of parameters and ε i are iid residuals.

The correct structure between the coefficients -or classification -in β is denoted by A = {(k, l) :

β k = β }.
A usual technical assumption is to consider designs the associated gram matrices of which converge to positive definite matrices. In the oneway ANOVA settings, we just need to assume that when n → ∞, then n k /n → ρ k < ∞ for all k = 1, . . . , K. We denote by D the corresponding asymptotic covariance matrix which is a K-diagonal matrix with diagonal entries equal to ρ 1 , . . . , ρ K .

In the univariate case like in ( 11), the estimator associated with Problem (2) using the 1 -norm for fusion is

β(n) = arg min β∈R K 1 2 K k=1 n k (ȳ k -β k ) 2 + λ n k = w k |β k -β |, (12)
which is just a rewriting of ( 6) where the dependency on n of the estimator and the tuning parameter is stated explicitly for the purpose of asymptotic analysis. Similarly, we denote by Ân = (k, ) :

β(n) k = β(n) the estimated group structure.
Exponentially adaptive weights and the fused-ANOVA. In this paragraph, we study the exponentially adaptive weights, which we recall here:

w FA k = n k n exp{-α √ n|ȳ k -ȳ | γ }, α, γ > 0.
We show that they enjoy some "oracle properties" in the sense of [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties[END_REF], that is, both i) right model identification (recovering the true classification A ) and ii) optimal estimation rate √ n. In the context of the penalized ANOVA problem (12), we denote these weights by w FA k and call the associated estimator the fused-ANOVA. These weights are adaptive as in the adaptive-Lasso of Zou (2006): it is known that raw 1 methods like the Lasso do not enjoy the aforementioned oracle properties, yet this can be fixed by choosing judicious weights that depend on an estimator of β which is asymptotically √ n-consistent -like the ordinary least squares, which equals (ȳ 1 , . . . , ȳK ) in the case at hand. Here we are interested in the differences between the entries of β; thus the quantity √ n|ȳ k -ȳ | seems quite natural in (9). While studying the asymptotic of our estimator, we came across the proposal of Bondell and Reich (2008b) for adaptive weights: they consider Problem (12) with additional constraints on the β k 's -that must sum to zero -and the following weights, which we refer to as the Cas-ANOVA weights:

w CA k = √ n k + n |ȳ k -ȳ | . ( 13
)
As we shall see, though quite interesting, Cas-ANOVA weights are adaptive on a smaller range of λ n than are fused-ANOVA weights. Moreover, they lead to splits. Thus, we believe that fused-ANOVA is computationally and statistically more efficient for solving Problem (12). We now proceed to the Theorem stating the required conditions on λ n for the fused-ANOVA to enjoy the oracle properties.

Theorem 3 (Oracle properties). Suppose that λ n n 3/2 exp {-α √ n} → 0 and λ n n 3/2 → ∞ when n → ∞. Then the fused-ANOVA enjoys asymptotic normality and consistency for recovering the true classification, i.e.,

√ n β(n) -β → d N (0, σ 2 D -1 ) and P( Ân = A ) → 1 when n → ∞.
The proof is postponed to Appendix A.3 and roughly follows that of Zou. We have, however, some comments related to this Theorem.

Remark 2 (On the exponentially adaptive weights). The key idea behind this theorem is that when n goes to infinity, then w FA k /

√ n goes to infinity if (k, ) ∈ A and to zero exponentially fast if (k, ) / ∈ A . This is due to the joint effect of the √ n-consistency of the ȳk and of the exponential. This is to be compared with Cas-ANOVA weights, where, when n → ∞,

w CA k / √ n goes to infinity if (k, ) ∈ A , but only to a constant if (k, ) / ∈ A .
Remark 3 (On the range of λ n ). Theorem 3 is true for a large range of λ n values. In particular it is true for a constant λ n . Asymptotically all groups belonging to the same class fuse almost immediately (i.e., for small values of λ of the order n 3/2 exp {-α √ n}) and the groups belonging to different classes fuse for very large λ, i.e., of the order n 3/2 . Numerical illustration in the univariate case. We generate data from model (11) as follows, for K the number of prior groups and n being fixed: the true vector of parameters β is composed of K entries picked up randomly among {1, 2, 3}, such that the correct structure A is always composed of 3 groups. Then, the initial group sizes n k are drawn from a multinomial distribution M(n, (p 1 , . . . , p K )) with p k = 1/K for all k = 1, . . . , K, such that the n k are approximately balanced. Finally, we let ε i ∼ N (0, 1) to generate the vector of data y = (y 1 , . . . , y n ).

We compare the capability of three weighting schemes to recover the true grouping A , namely the fused-ANOVA weights, the Cas-ANOVA weights, and the so-called default weights corresponding to w k = n k n , which are not adaptive but produce a path of solutions that contains no split. Such weights correspond to the Clusterpath weights adapted to the ANOVA setup. We use our own code for each method. Typically, the computational burden required by Cas-ANOVA is huge, compared to the other two procedures as the path of solutions may contain splits. Qualitatively, the difference would be as in Figure 3, middle panel. Thus, we typically force the algorithm not to split when using the Cas-ANOVA weights.

We generate data as specified below, and for each procedure we check whether there exists at least one λ for which the correct structure is identified along the path of solutions. The probability of true support recovery is evaluated by replicating this experiment a large number of times (8096 times4 ). To investigate the asymptotic behavior of each method, we vary n from 50 to 1,000 and consider two scenarios for the initial number of groups K. First, K is fixed at 10 such that the number of elements in each group grows with n. In the second scenario, K grows with n through the relationship K = 2.5 • log(n). The results are reported on Figure 4, with the first (resp. the second) scenario on the left (resp. the right) panel. The results confirm Theorem 3. The two adaptive procedures, Cas-ANOVA, and to a greater extent, fused-ANOVA, dominate the non-adaptive weights. As expected from Section 7, fused-ANOVA always dominate Cas-ANOVA, as experienced in other scenarios (e.g., K = C • √ n) not reported here to save space.

a) K = cst. b) K = C log(n)
P( Ân = A ) q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 250 500 750 1000 n q q q q q q q q q q q q q q q q q q q 0.00 Numerical illustration in the bivariate case. Theorem 3 characterizes the asymptotic of the fused-ANOVA estimators when considering one dimension at a time. Concerning the multidimensional setting, there are two situations. On the first hand, there exists a dimension such that all the true groups are different, i.e β kj = β j . In this case, our theorem guarantees that, using this particular dimension, the recovered classification will converge to the true one. On the second hand, there exists no dimension such that the true groups are all different. In that case, we have no theoretical guarantee to support the fused-ANOVA weights. It is nonetheless possible to aggregate the classification obtained in each dimension to a consensus classification. For a given λ, two individuals k and are in the same multidimensional cluster if they have been fused on every dimension.

In order to evaluate empirically the performance of the aggregation step, we consider a two dimensional classification problem with three classes and two scenarii. Each prior group is drawn from one of three classes. In the first scenario, the three classes have different means on the first dimension and the same mean on the second dimension. The mean vectors are (1, 1.5); (2, 1.5); (3, 1.5), as in top left panel of Figure 5. In the second scenario, both dimensions are informative: the first dimension separates classes {1, 2} from {3} while the second dimension separates classes {1, 3} from {2}. The mean vectors are (1, 1); (1, 2); (2, 1), as in top right panel of Figure 5). We increase the difficulty in each scenario by adding a Gaussian noise with increasing standard deviation σ. Results in Figure 5 corresponds to the estimated probability of true classification recovery along the path, averaged over 2,000 runs.

In both scenarii, the fused-ANOVA weights with aggregation outperform the multidimensional 2 -Clusterpath as well as the single linkage hierarchical clustering. The Ward hierarchical clustering shows better performance but at a much higher computational cost.

In this simple multidimensional numerical study, we always aggregate the classification over the dimensions. However, this aggregation is not necessarily better than performing classification on single, well-chosen feature. We illustrate this point in the following section on phylogenetic data.

ȳk,1 σ= 0.25 σ= 1.25 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 1 2 3 1 2 3 σ= 0.25 σ= 1.25 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 0 1 2 0 1 2 ȳk,2 ȳk,2 P( Ân = A ) q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 σ q q q q q q q q q q 0.00

A complete example in phylogeny

Evolutionary trees -sometimes referred to as "trees of life" -are one of the most emblematic hierarchical representations in computational biology. They are typically used in phylogenetics to compare biological species based on their similarities regarding one or several features. These features could be phenotypic traits or genetic characteristics. In these tree structures, each node corresponds to a taxonomic unit, the root node being the most recent common ancestor to all leaves on the tree. All other intermediate nodes between root and leaves represent the taxonomic knowledge between the species of interest. The study depicted in [START_REF] Vetrovsky | The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses[END_REF] enters this framework by more specifically considering features associated with bacterial genomes to determine the phylogenetic relationships between the taxa. The data set consists of various genetic features associated with n = 1, 690 complete bacterial genomes classified in K = 903 known bacterial species. We apply our method on this data set to assess its capability of capturing the true underlying taxonomic structure. To do so, we consider the following genetic features to construct the hierarchy: the number of known genes, the number of known proteins and the genome size (measured by the number of bases in millions). We apply the univariate model ( 6) on each feature to reconstruct a tree structure. The indexing function κ is built from the lowest level of classification available that splits the genomes into K = 903 bacterial species. We use the default weights and the fused-ANOVA weights (9) with α chosen specifically for each feature (see below). We also apply hierarchical clustering using Ward's criterion and starting from the known classification in bacterial species. Hierarchical clustering is applied individually on each feature, as well as across the three features using the Euclidean distance to build the similarity matrix. To assess the relevance of the inferred trees, we compare them with various levels of the known taxonomic classification above the species level, namely genus (470 groups), family (216 groups), order (100 groups), class (46 groups) and phylum (27 groups). To this end, we compute the best adjusted randindex between the respective reference classifications and the classifications obtained by cutting an inferred tree at all the possible levels of the hierarchy. As an example, we report in Figure 6 a subset of the tree inferred by the fusion penalty with fused-ANOVA weights and the cutting level that leads to the best performance in terms of adequacy with the true phylum taxonomy. More quantitative results are reported in Figure 7, with the adjusted rand indexes for the taxonomic classifications in terms of phylum, order and family, using either the number of genes, the number of proteins or the genome size as the feature variable for classification. We also represent the consensus/multidimensional classifications either obtained by aggregating the three univariate fused-ANOVA trees or by considering the three features together for Ward hierarchical clustering. Note that for the fused-ANOVA weights, we apply our method on a grid of α and report the results obtained for the best α in terms of adjusted rand-index.

First, we notice that the fused-ANOVA weights always outperform the default weights. This is expected since the former weights are a special case of the latter when α → 0. Second, we note that the consensus classification -or the one obtained by multivariate hierarchical clustering -is not always the best choice. This is particularly obvious for the phylum classification, where the "size" feature leads to very good results in terms of adjusted rand-index. These results considerably deteriorate for the consensus classification, due to the relatively poor results obtained from the "genes" and "proteins" features. Finally, the most striking result in Figure 7 is that the fusion penalty approaches clearly outperform the Ward hierarchical clustering. At first glance, one might argue that the weighting scheme used in fused-ANOVA is responsible for such good performance. However, the fusion penalty with default weights remains competitive in a few cases. This supports the fact that the regular-Adjusted Rand-Index izing virtue of the fusion penalty is of great help when the problem size is high.

Vivian Viallon, Sophie Lambert-Lacroix, Hölger Hoefling, and Franck Picard. On the robustness of the generalized fused lasso to prior specifications.

A Proofs

A.1 Theorem 1 (absence of splits with norms)

For the sake of brevity the proof is detailed only in the clustering framework, i.e., when κ(i) = i and n k = 1 for all k = 1, . . . , K. The generalization to groups with more than one individual is straightforward and follows the exact same line.

Consider the objective function in (2) and a time λ 0 at which we have a valid set of clusters. It is obvious that clusters containing only one individual cannot split. We will thus consider clusters grouping together more than one element. We denote by C = {k : β k (λ 0 ) = β C (λ 0 )} such a cluster, with β C the current estimated mean. For unitary weights, Lemma 1 implies that

0 p = -ȳ C + β C + λ 0 i / ∈C ∂Ω(β C -β i ) ∂β C (λ 0 ).
Subtracting the above equation from the subgradient equation ( 3) for i ∈ C, one has ȳC -

y i + λ 0 j∈C τ ij (λ 0 ) = 0 p . (14)
We now consider any time λ ≥ λ 0 such that no fusion has occurred. Let us show that for τ ij (λ) = λ0 λ τ ij (λ 0 ), it is possible to solve the KKT conditions, and thus show that no split occurs.

First, the proposed τ ij (λ) are valid subgradients as Ω(τ ij (λ)) ≤ 1 since Ω(τ ij (λ 0 )) ≤ 1 and λ > λ 0 . Second, for this particular choice of subgradient and thanks to (14), the KKT conditions for all C and all i ∈ C simplify as follows:

β C -y i + λ j∈C λ 0 λ τ ij (λ 0 ) + λ C =C |C | ∂Ω(β C -β C ) ∂β C (λ) = β C -ȳC + λ C =C |C | ∂Ω(β C -β C ) ∂β C (λ).
It now remains to check that we can find a β which zeroes this subgradient equation. Note that for all C = C, the differential ∂Ω(β C -β C )/∂β C (λ) is well defined. Then, by multiplying the above expression by |C|, we obtain the gradient of the following objective function

1 2 n i=1 y i -β i 2 2 + λ C,C :C =C |C| • |C | Ω(β C -β C ).
This is a strictly convex problem admitting one unique solution which is solved by zeroing its gradient. Thus we necessarily have

β C -ȳC + λ C =C |C | ∂Ω(β C -β C ) ∂β C (λ) = 0 p ,
which ends the proof.

A.2 Theorem 2: absence of splits with distance-decreasing weights in 1-d

For the sake of brevity the proof is detailed only in the case where κ(i) = i and n k = 1 for all k = 1, . . . , K. The generalization to groups with more than one individual is straightforward, seeing that we can replace a group κ(i) by n κ(i) individuals with value j y j /n κ(i) . Also, when Ω ≡ 1 , the proof remains valid but should be done separately on each dimension.

Throughout the proof, we may thus consider the estimator defined by

β(λ) = arg min β∈R n 1 2 n i=1 (y i -β i ) 2 + λ i,j:i =j w ij |β i -β j |. (15)
The proof proceeds in two steps detailed hereafter:

1. in subsection A.2.1, we show that absence of splits is equivalent to preservation of the order along the path;

2. in subsection A.2.2, we show that distance decreasing weights preserve the order, by considering a dual formulation of Problem (15).

For simplicity, we consider that the data vector y is initially ordered such that

y 1 ≥ . . . , y i ≥ y i+1 ≥ • • • ≥ y n .

A.2.1 Preserving the order

We say that the loss is order-preserving, if y i ≤ y j implies that βi (λ) ≤ βj (λ), for all λ ≥ 0.

Lemma 2. The absence of splits is equivalent to preservation of the order along the path for Problem (15).

Proof. First of all, in the absence of splits in the path, it is clear that the order is preserved.

Conversely, assume that there is an event at λ 0 that splits a group C into C down and C up , where βdown (λ) < βup (λ) for all λ ≥ λ 0 . By means of Equation (4), we necessarily have ȳdown ≥ ȳup as illustrated on Figure 8. However, if the order is preserved, for all (i, j) ∈ C down × C up , we must have y i < y j and ȳdown < ȳup , which leads to a contradiction.

A.2.2 The dual problem

We follow arguments developed by [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] for the generalized Lasso. Indeed, Problem (15) can be recast as

β(λ) = arg min β∈R n 1 2 y -Xβ 2 2 + λ W Dβ 1 , (16)
a generalized Lasso problem with X = I nn , W a diagonal matrix the diagonal of which is the n(n -1)/2 vector given by diag(W ) = (w 11 , . . . , w 1n , w 23 , . . . , w 2n , w 34 , . . . , w (n-1)n )

and D is a n(n -1)/2 × n matrix that performs the pairwise differences such that

D =                  couple (i,j) (1,1) 1 -1 (1,2) 1 -1 . . . . . . . . . (1,n) 1 -1 (2,3) 1 -1 . . . . . . . . . (2,n) 1 -1 . . . (n-1,n) 1 -1                  . ( 17
)
In what follows, it will be convenient to index rows of the matrix D in terms of the couple (i, j), as is done in Expression ( 17).

We then rely on the Lagrangian dual of the primal problem ( 16) studied in Tibshirani and Taylor (2011), which is û(λ) = arg min u∈R (n(n-1)/2)

1 2 y -D T W u 2 2 subject to u ∞ ≤ λ, (18)
and where the correspondence between the primal and dual variables is

β = y -D T W û.
The dual solution must satisfies ûij ∈

     {+λ} if (W D β) ij > 0, {-λ} if (W D β) ij < 0, [-λ, +λ] if (W D β) ij = 0,
where we use the indexing in terms of (i, j) for the vector u. We also define B, the set of (i, j) such that |u ij | = λ, that is, the ones reaching the boundary in the dual.

The key point is to note that the order is not preserved if and only if, at some point of the path, there exists some (i, j) and λ such that ûij (λ) = -λ, meaning that (W Dβ) ij < 0. The rest of the proof will show that this event is not possible for distance decreasing weights and the matrix D given by (17). To this end, we proceed by contradiction, by supposing that the order is not preserved along the path. We thus consider the first split event that will disrupt the order, which occurs at λ 0 . At this point, the order is preserved and there is an ε > 0 such that on ]λ 0 , λ 0 + ε], the order is not preserved. We note that λ 0 > 0 since the order is necessarily preserved up to the first fusion event that fuses data points with different values. At λ 0 , we must have a couple (i 0 , j 0 ) such that ûi 0 j 0 (λ 0 ) = -λ 0 that reaches the boundary. Moreover, the left derivative ∂ -ûi 0 j 0 (λ) must be less than -1 because the path is continuous (see [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] and because we consider the first event disrupting the order. We provide geometrical insight into this point on Figure 9. We now show that we necessarily have ∂ -ûi 0 j 0 (λ) > -1, leading to a contradiction. To this end, we consider the set C of indices which are fused with i 0 and j 0 just before λ 0 , that is, C = i : βi = βi 0 = βj 0 . We denote by Just before λ 0 and for all (i, j) ∈ D out , we have ûij (λ) = λ and so sign(û ij ) = 1, because the order is preserved at this point. Then, following [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF], the KKT conditions of ( 18 where A + denotes the Moore-Penrose pseudo-inverse of A. Note that such a choice is important since it guarantees that û(λ) is a continuous function of λ.

Expression of (19) greatly simplifies by exploiting an explicit formula for the pseudo-inverse, which we derive in the next paragraph.

Analytic form of the pseudo-inverse. In this paragraph, we consider only the D Din×C matrix and the W D 2 in matrix, which correspond to the set of intra C differences and their weights. For simplicity, we just denote them D and W here and call n the group size. We have Let us consider the size-n C vector V , which includes the differences between elements in C and elements outside C. Note that the ith column of D Dout×C is zero everywhere, except for the elements of D out containing i. In the last case, it is equal to 1 if y i ≥ y j and to -1 otherwise. Hence,

V i = (D Dout×C ) T W D 2 out 1 Dout i = j∈ C w ij sign(y i -y j ).
Also recall that the matrix D Din×C encodes the pairwise positive differences. Then for y i > y i , the (i, i ) element of D Din×C V equals (D Din×C V ) ii = V i -V i = j∈ C w ij sign(y i -y j ) -w i j sign(y i -y j ).

There are two possibilities: either y j > y i ≥ y i or y i ≥ y i > y j . We thus split the summation in the above equation into two parts:

(D Din×C V ) ii = j∈ C yj >yi≥y i (w i j -w ij ) + j∈ C yi≥y i >yj (w ij -w i j ).
which is in contradiction with the rest of the subgradient equation of β k since we recall that the left-hand side is O P (1). Therefore we must have P (k, ) ∈ Ân → 1 for all (k, ) ∈ A , which completes the proof of the consistency part in Theorem 3.

Figure 1 :

 1 Figure 1: Reconstructed phylogenetic trees for various weighting schemes. Families classified in the same order share the same color.

Figure 3 :

 3 Figure 3: timing experiments: a) time in seconds as a function of the number of conditions K; b) timing comparison for general/without-split algorithm; and c) timing comparison for naive/embedded cross-validation.

Figure 4 :

 4 Figure 4: Univariate case: estimated probability of consistency as a function of the sample size n, for various weights and in two scenarios: the number of initial groups K is either a) fixed to a constant (10) or b) increases in C log(n) with C = 2.5. The true number of groups in A is 3.

Figure 5 :

 5 Figure 5: Bivariate example: estimated probability of consistency as a function of the noise standard deviation σ, for various clustering methods. The initial number of groups K fixed to a constant (10). The true number of groups in A is 3.

Figure 6 :

 6 Figure6: Tree reconstructed with fused-ANOVA from the "Size.Mb" trait (Adjusted Rand-Index=0.71). Projected colors correspond to the true taxonomy (phylum).

Figure 7 :

 7 Figure 7: Adequacy of various tree-based clustering methods to different levels of phylogenetic classification.

Figure 8 :

 8 Figure 8: Equivalence between preserving the order and absence of splits relies on a simple geometrical argument.

Figure 9 :

 9 Figure 9: Geometrical insight into a event in the dual.

D

  in = {(i, j) ∈ C × C : i < j} the set of intra C differences andD out = {(i, j) ∈ B : i < j, i ∈ C or j ∈ C} ,the set of differences between C and other groups. Finally we denote by R the set of all other indices which are not in D out and D in . Given those sets we reindex the matrix D and W as follows. for all (i, j) ∈ R, i and j do not belong to C and thus D R×C = 0. By simple matrix algebra, the restriction of D T W to the rows in C can be written(D T W ) C = (D Din×C ) T W D 2 in + (D Dout×C ) T W D 2 out .

  ) restricted to D in imply that ûDin (λ) = W D 2 in D Din×C D T Din×C W D 2 in + W D 2 in D Din×C y -λ(W D 2 out D Dout×C ) T 1 Dout ,(19)

D

  T D = n I n n -1 n 1 T n , D1 n = 0 n , and from this we get DD T D = n D and thus, D + = D T /n . Finally (DD T ) + = 1 n 2 DD T , and (DD T ) + D = D n .If we now consider the weighted version of Problem (16), one has(W DD T W ) + W D = W -1 (DD T ) + W -1 W D = W

  Statistics and Computing, pages 1-17, 2014. H. Zou. The adaptive lasso and its oracle properties. JASA, 101(476):1418-1429, 2006.

publicly available at http://genomics.senescence.info/species/

With this simulation setting, there is no underlying clustering since our point is to compare run times here.

We numerically study the multidimensional case at the end of this section.

this number arises from the manifold computer cores available.

Acknowledgments

We would like to thank Mahendra Mariadassou for useful discussions about the bacterial genomes data set and for sharing his knowledge of phylogeny.

And from this we see that if the weights are positive and distance decreasing, all the (D Din×C V ) ii are negative. To conclude, the slopes in Expression (20), that is,

are positive, which is in contradiction with ∂ -ûi 0 j 0 (λ) ≤ -1.

A.3 Theorem 3: consistency for exponentially adaptive weights

We essentially follow the same line as for the adaptive Lasso in Zou ( 2006), yet adapted to the fusion penalty as in Viallon et al. (2014); Bondell and Reich (2008b).

The main difference comes from the use of the exponentially adaptive weights w FA k .

We start by asymptotics in the vein of [START_REF] Fu | Asymptotics for lasso-type estimators[END_REF] for Lasso-type procedures: Lemma 3 below gives the limiting distribution of the fused-ANOVA estimator (12) on the range of interest for the penalty λ n which essentially proves the asymptotic normality part of the Theorem.

where, for W ∼ N (0, σ 2 D),

) -where β is the true vector of parameters and

Note that u n is also the minimizer of

Let us study the limiting behavior of V n . The basic assumptions for our fused-ANOVA Problem ( 12) are having a design such that lim n→∞ n k /n = ρ k and having i.i.d residuals with zero mean and common variance σ 2 . Thus, the first two terms in V n respectively converge to a constant u T Du and to a Gaussian W = N (0, σ 2 D), where D is a K-diagonal matrix such as D kk = ρ k . For the third term, there are two possibilities: either β k = β or β k = β . In other words, (k, ) belongs to A or does not. First note that

In words, this part of the third term converges to a finite constant in both situations which is null as soon as u k = u . Second, consider the remaining part of this third term which involves the weights w FA k . It suffices to use the √ n-consistency of the OLS estimators (ȳ 1 , . . . , ȳK ) coupled with assumptions made on the limiting behavior of λ n to see that

Application of Slutsky's Lemma gives the limiting behavior of the third term in V n and we finally get V n (u) → V (u) with V defined as in Lemma 3. The final convergence of u n → d arg min u V (u) is obtained by applying epiconvergence results of [START_REF] Geyer | Simultaneous factor selection and collapsing levels in ANOVA[END_REF].

Turning back to the proof of Theorem 3, just note that the unique minimizer of the convex function

) such that u k = u for all (k, ) ∈ A and the asymptotic normality part is proved.

We now proceed to the consistency in terms of support recovery. First, concerning the elements of β that should not fuse according to the true A , Lemma 3 indicates that

Second, regarding elements of β that must fuse, we need to prove that

To this end, we proceed as in Viallon et al. ( 2014) to get a contradiction by considering the largest βk such that βk = β even though β k = β . This can be done by inspecting the KKT conditions asymptotically. In the univariate case and for Ω the 1 -norm, an optimal β verifies the following subgradient equation for all k = 1, . . . , K. This is written

Now, in the first term of the right-hand side, suppose that there exists at least one such that (k, ) ∈ A and βk = β simultaneously; consider βk with k = arg max :(k, )∈A { β }, the largest coefficients verifying these conditions: we must have τ k = 1 for all such that β = βk and β = β k . Now if we sum equation ( 21) for all that are fused with k we obtain:

By Lemma 3 and asymptotic normality, the left-hand side in (22) converges to a O P (1). Then, the second term on the right-hand side (that is, elements that should not fuse) tends to 0 since λ n w FA k / √ n → 0 when (k, ) / ∈ A , as seen previously. Finally we have :