Iosif Bena
email: iosif.bena@cea.fr

Hagen Triendl
email: hagen.triendl@cea.fr

Bert Vercnocke
email: bert.vercnocke@cea.fr

Camouflaged Supersymmetry

We establish a relation between certain classes of flux compactifications and certain families of black hole microstate solutions. This connection reveals a rather unexpected result: there exist supersymmetric solutions of N = 8 supergravity that live inside many N = 2 truncations, but are not supersymmetric inside any of them. If this phenomenon is generic, it indicates the possible existence of much larger families of supersymmetric black rings and black hole microstates than previously thought.

Introduction

There is an extensive body of work on obtaining supersymmetric and non-supersymmetric vacua for flux compactifications of string theory and studying their phenomenology, and a parallel extensive body of work on constructing supersymmetric and non-supersymmetric black hole microstate solutions to understand black hole physics in string theory. While the physical motivations are different, the technical tools are rather close. In particular, the equations underlying supersymmetric solutions are well-understood and classified: On the flux compactification side (see for example [START_REF] Grana | Supersymmetric three form flux perturbations on AdS(5)[END_REF][START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF][START_REF] Grana | Supersymmetric backgrounds from generalized Calabi-Yau manifolds[END_REF]) in ten dimensions, on the black hole microstate side for the underlying supergravity in five dimensions [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. Furthermore, some of the methods for constructing non-supersymmetric solutions from supersymmetric ones are strikingly similar. These methods include slightly deforming the supersymmetric solution by additional fluxes [START_REF] Grana | Supersymmetric three form flux perturbations on AdS(5)[END_REF][START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF], flipping some signs [START_REF] Goldstein | Almost BPS black holes[END_REF], or writing some effective Lagrangian as a sum of squares for black holes [START_REF] Ferrara | Black holes and critical points in moduli space[END_REF][START_REF] Denef | Supergravity flows and D-brane stability[END_REF][START_REF] Ceresole | Flow equations for non-BPS extremal black holes[END_REF][START_REF] Andrianopoli | First Order Description of Black Holes in Moduli Space[END_REF][START_REF] Cardoso | First-order flow equations for extremal black holes in very special geometry[END_REF] or flux backgrounds [START_REF] Lust | Generalized non-supersymmetric flux vacua[END_REF][START_REF] Held | DWSB in heterotic flux compactifications[END_REF].

It is therefore not surprising that one can find a relation between certain types of solutions on the two sides. Indeed, as we will show below, certain supersymmetric flux backgrounds of the type [START_REF] Becker | M theory on eight manifolds[END_REF] where the "internal" (non-compact) manifold contains a hyper-Kähler factor can be interpreted as certain non-rotating solutions in the classification of [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. (One can similarly relate non-supersymmetric solutions. The story is more intriguing and will be alluded to in this letter, but we leave the details for a companion publication [16].) The main purpose of this letter is to show that there are other supersymmetric solutions of the same class of flux compactifications which, when interpreted as black hole microstates in N = 2 supergravity, do not fall into the classification of supersymmetric solutions [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. 1 Hence, from the point of view of N = 2 supergravity, these solutions should be non-supersymmetric. However, they are supersymmetric inside N = 8 supergravity! As we will explain below, these solutions have the right field content to fit into many possible N = 2 truncations, and hence they will always be solutions of these N = 2 theories. However, the unbroken supercharges are projected out in all possible N = 2 truncations and hence from the point of view of N = 2 supergravity none of these solutions are supersymmetric2 . A simple way to understand this is to recall that all N = 2 supersymmetric solutions in the class [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF] have (in our conventions) anti-self-dual fields on a hyper-Kähler base, while our solutions have both anti-self-dual and self-dual fields.

Our results have quite a few unexpected implications. First, it is widely believed that all supersymmetric microstate geometries of three-charge black holes in five dimensions are described by the equations of [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. Our results indicate that many solutions that are not described by these equations are also supersymmetric in the parent N = 8 theory. This implies that beside the classes of microstate solutions constructed so far there may exist many more supersymmetric microstates, which would contribute to the entropy count.

Second, it has been conjectured [START_REF] Ferrara | Magic supergravities, N= 8 and black hole composites[END_REF] and argued [START_REF] Bossard | 1/8 BPS Black Hole Composites[END_REF] that all multicenter supersymmetric solutions of N = 8 supergravity must live inside an N = 2 truncation, and one may believe that this implies that the solutions of [START_REF] Denef | Supergravity flows and D-brane stability[END_REF][START_REF] Bates | Exact solutions for supersymmetric stationary black hole composites[END_REF] capture all supersymmetric multicenter N = 8 solutions. Our results show that this is not so.

Third, it is well-known that the supersymmetric black ring in five dimensions [START_REF] Elvang | A supersymmetric black ring[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF][START_REF] Elvang | Supersymmetric black rings and three-charge supertubes[END_REF][START_REF] Gauntlett | General concentric black rings[END_REF]] is part of a truncation to N = 2 ungauged supergravity and belongs to the class of solutions [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. Our results indicate that there may exist a new, more general supersymmetric black ring with more dipole charges (coming from the extra self-dual fluxes). Besides its interest as a new solution, if this black ring existed, it may also help to account for the missing entropy between the D1-D5 CFT and the dual bulk in the moulting black hole phase [START_REF] Bena | Moulting Black Holes[END_REF].

More generally, the relation between black hole microstates and flux compactifications that we outline will likely prove fruitful in both directions. There exists a whole methodology for constructing flux compactifications by writing the effective Lagrangian governing these compactifications as a sum of squares of calibrations [START_REF] Lust | Generalized non-supersymmetric flux vacua[END_REF][START_REF] Held | DWSB in heterotic flux compactifications[END_REF]. Under the guise of "floating branes", calibrations have also been used to find non-supersymmetric black hole microstates [START_REF] Bena | Supergravity Solutions from Floating Branes[END_REF], and relating the two approaches is likely to yield novel classes of solutions on both sides. We plan to report on this relation in an upcoming companion paper [16]. Furthermore, it has been recently discovered that even some non-extremal cohomogeneity-two black holes, black rings and microstates are calibrated [START_REF] Bena | Imaginary Soaring Branes: A Hidden Feature of Non-Extremal Solutions[END_REF]. If one could use this to write down a new decomposition of the effective Lagrangian (similar to the one of non-extremal cohomogeneity-one solutions [START_REF] Gibbons | Antigravitating Black Hole Solitons with Scalar Hair in N=4 Supergravity[END_REF][START_REF] Miller | Nonextremal black holes are BPS[END_REF][START_REF] Perz | First-order flow equations for extremal and non-extremal black holes[END_REF][START_REF] Galli | Non-extremal black holes of N=2, d=4 supergravity[END_REF]) one would obtain a systematic method to construct new highly-non-trivial and physically-interesting solutions.

The Solution

We focus on a class of solutions to five-dimensional N = 8 supergravity that arises as the lowenergy limit of a T 6 compactification of eleven-dimensional supergravity. The spatial part of the five-dimensional spacetime is given by a hyper-Kähler space M 4 , and the warp factor A depends only on the M 4 coordinates. The full eleven-dimensional metric is

ds 2 11 = -e -2A dt 2 + e A ds 2 (M 4 ) + e A (dx 2 5 + dx 2 6 + dx 2 7 + dx 2 8 ) + e -2A (dx 2 9 + dx 2 10 ) , (1)
with coordinates x 5 . . . x 10 on T 6 . The four-form field strength is

F mag 4 = d(e -3A ) ∧ dt ∧ dx 9 ∧ dx 10 + [Θ + -Θ -] ∧ dx 5 ∧ dx 8 + [Θ + + Θ -] ∧ dx 6 ∧ dx 7 + Θ+ ∧ (dx 6 ∧ dx 8 -dx 5 ∧ dx 7 ) (2)
where Θ + , Θ+ are self-dual two-forms on M 4 and Θ -is an anti-self-dual one. With hindsight, we focus on a solution whose self-dual forms obey the relation

(Θ + + i Θ+ ) ∧ (Θ + + i Θ+ ) = 0 , (3)
which implies that Θ + + i Θ+ defines a complex structure on M 4 under which it is a holomorphic two-form. As we will see in Section 2.1, this ensures that the solution is supersymmetric. Finally, the warp factor is determined by

∆ 4 e 3A = (Θ 2 + + Θ2 + + Θ 2 -) + ρ M 2 , (4)
where ∆ 4 is the Laplacian on M 4 and ρ M 2 the M2 brane density. This solution has the electric charge of a set of M2 branes extended along the x 9 and x 10 directions and smeared on the other compact directions of T 6 . The magnetic component of the four-form can be thought of as being sourced by four types of M5 branes on the corresponding Poincaré dual cycles. We summarize that in Table 1.

0 9 10 5 6 7 8 M 4 M2 × × × M5 × × × × × γ 1 M5 × × × × × γ 2 M5 × × × × × γ 3 M5 × × × × × γ 4
Table 1: The brane charges for our configurations along the T 6 directions x 5 . . . x 10 . A brane is localized in directions marked "×" and smeared in the other ones. The M5 branes each wrap a 1-cycle γ i in the hyper-Kähler space M 4 , determined by the (anti)-selfdual fields Θ ± , Θ+ .

Interpretation as a flux compactification

We now argue that this solution is a supersymmetric solution of 11-dimensional supergravity. By swapping the roles of M 4 and T 2 9,10 as external and internal spaces, we see that the above solution is actually an eight-dimensional Calabi-Yau 'compactification' of M-theory, of the type discussed first in [START_REF] Becker | M theory on eight manifolds[END_REF]. The eleven-dimensional spacetime has the form M 1,10 = M 1,2 × X 8 , where X 8 = M 4 × T 4 5,6,7,8 . The metric and the gauge field preserve three-dimensional Poincaré invariance, as can be seen by rewriting (1) and (2) as ds 2 11 = e -2A (-dt 2 + dx 2 9 + dx 2 10 ) + e A ds 2 (X 8 ) ,

F 4 = d(e -3A vol 3 ) + Im [(Θ + -i Θ+ ) ∧ dz ∧ dw + Θ -∧ dz ∧ d w] , (5)
where vol 3 = dt∧dx 9 ∧dx 10 is the volume form of three-dimensional spacetime and A only depends on the coordinates of the internal manifold X 8 . Furthermore, we defined the holomorphic oneforms dz = dx 5 + i dx 6 , dw = dx 7 + i dx 8 .

The supersymmetry conditions require ds 2 (X 8 ) to be a Calabi-Yau metric for X 8 and the internal components of F 4 to be a primitive (2, 2)-form. The first two requirements are fulfilled since (1) and ( 2) give a Calabi-Yau metric ds 2 (X 8 ) = ds 2 (M 4 ) + dzdz + dwd w .

Since the anti-self-dual two-forms on hyper-Kähler manifolds are (1, 1), eq. ( 5) implies that the internal components of F 4 indeed make up a primitive (2, 2)-form if (Θ + + i Θ+ ) ∧ dz ∧ dw is the holomorphic four-form of X 8 (such that (Θ + -i Θ+ ) is antiholomorphic on M 4 ). This in turn can only be realized if condition (3) holds. The equation of motion for the gauge field then determines the warp factor in general as

d * 8 d A = 1 6 F mag 4 ∧ F mag 4 , (8)
which reduces to (4) when X 8 = M 4 × T 4 5,6,7,8 . Note that the described background is dual to a supersymmetric flux background of IIB string theory in the GKP class [START_REF] Grana | Supersymmetric three form flux perturbations on AdS(5)[END_REF][START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF].

Relation to five-dimensional STU solutions

Finally, we can interpret our supersymmetric solution in eleven-dimensional supergravity compactified on a six-torus (T 6 (5,6,7,8,9,10) ) which descends to five-dimensional N = 8 supergravity. There exists a very large class of solutions to this theory, that fit inside an N = 2 truncation with two vector multiplets: they describe black rings, black holes as well as microstate solutions that have the same charges as these objects but no horizon.

All supersymmetric solutions of this truncation are known [START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF][START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF], and are given by:

ds 2 11 = -Z -2 (dt + k) 2 + Z ds 2 4 + Z 3 I=1 ds 2 I Z I , F 4 = dA (I) ∧ ω I = 3 I=1 -d dt + k Z I + Θ (I) ∧ ω I , (9)
where Z ≡ (Z 1 Z 2 Z 3 ) 1/3 , ds 2 I and ω I are respectively a unit metric and a unit volume form on the three T 2 's inside T 6 and ds 2 4 is a four-dimensional hyper-Kähler metric. When this metric has a translational U(1) isometry it becomes a Gibbons-Hawking metric; if one then compactifies along the Gibbons-Hawking fiber, one obtains a solution of the four-dimensional STU model. Note that we work in a convention in which the three curvature two-forms of the hyper-Kähler base are self-dual, and hence the Θ (I) of a supersymmetric solution are anti-self-dual.

The metric and the timelike (electric) components of the four-form of our solution (1,2) are of the form (9) with Z 1 = Z 2 = 1 and k = 0. However, the spacelike (magnetic) four-form field strengths have more components, and only reduce to the N = 2 truncation above when Θ + = Θ+ = 0. Hence, despite having the right electric charges, the supersymmetric N = 8 solution we found does not fit into the standard "STU" N = 2 truncation. In the next section we discuss the supersymmetry of this solution, and how it fits into a larger N = 2 truncation.

3 Supersymmetry in N = 8 and N = 2

We have shown already in Section 2.1 that the solution (1, 2) is a Calabi-Yau four-fold flux background, and hence preserves at least four supercharges [START_REF] Becker | M theory on eight manifolds[END_REF]. We first analyze the supersymmetry in detail and then discuss whether the solution and its supercharges fit inside the largest N = 2 truncation of the N = 8 theory.

1/8 BPS solutions in N = 8 supergravity

Clearly, the hyper-Kähler background breaks half of the supersymmetry, as it admits only a covariant spinor of (say) positive chirality. This corresponds to the projection Γ 1234 η = -η, where η is a spinor on the internal eight-dimensional manifold. Furthermore, the flux F 4 breaks more supersymmetry. Its electric component (corresponding to an M2-brane charge along the 9, 10 directions) breaks another half of supersymmetry, by the projection Γ 12345678 η = η.

To understand how the magnetic components of F 4 affect the supersymmetry, it is best to choose an appropriate vierbein e i , i = 1, . . . , 4, on the hyper-Kähler space M 4 , such that ( 3) is fulfilled and we can identify the self-dual two-forms of (2) as

Θ + = θ + (e 1 ∧ e 3 + e 4 ∧ e 2 ) , Θ+ = θ + (e 1 ∧ e 4 + e 2 ∧ e 3 ) . (10)
The supersymmetry conditions / F η = 0 and / F m η = 0 [START_REF] Becker | M theory on eight manifolds[END_REF] will now contain an additional projector, which further halves the amount of supersymmetry. More precisely:

0 = 1 4! F ijkl Γ ijkl η = 1 4 [(Θ + ) ij Γ ij58 + ( Θ+ ) ij Γ ij68 ](1 -Γ 5678 )(1 -Γ 1234 )η -1 4 (Θ -) ij Γ ij58 (1 + Γ 5678 )(1 + Γ 1234 )η , (11)
where we have inserted the projectors 1 2 (1 ± Γ 1234 ) by making use of the (anti-)self-duality of Θ ∓ . The term containing the anti-self-dual flux Θ -vanishes on the Killing spinors annihilated by the two earlier projectors 1 2 (1 + Γ 1234 ) and 1 2 (1 -Γ 12345678 ), and this agrees with the known structure of BPS three-charge solutions, in which turning on an anti-self-dual field strength on the base does not affect supersymmetry.

For arbitrary self-dual forms Θ + , Θ+ , the first line is not zero and supersymmetry is broken. However, for the specific choice [START_REF] Ceresole | Flow equations for non-BPS extremal black holes[END_REF] this term contains a new projector:

0 = 2θ + Γ 1358 (1 + Γ 3456 )η , (12)
which is compatible with the first two. More generally, under the condition (3) we always find such a projector and the solution has four supercharges. It is not hard to see that the equations / F m η = 0 do not impose any extra conditions on the remaining Killing spinors, essentially because the flux pieces that are self-dual on the hyper-Kähler manifold always combine into the projector 1 2 (1 + Γ 3456 ), while the anti-self-dual components always give either 1 2 (1 + Γ 1234 ) or 1 2 (1 + Γ 5678 ), depending on the index m. Therefore, the solution is 1/8 BPS, and its 4 Killing spinors are annihilated by the projectors:

1 2 (1 + Γ 1234 ) , 1 2 (1 + Γ 3456 ) and 1 2 (1 + Γ 5678 ) . ( 13
)

A puzzle

The 1/8 BPS solution we gave in (1,2) has not been found in the literature. Moreover, its magnetic field strength (2) has both self-dual and anti-self-dual components on the hyper-Kähler space. This is surprising since all 1/2 BPS solutions in N = 2 supergravity in five dimensions have only anti-self-dual fluxes on the hyper-Kähler space, as shown in [START_REF] Gauntlett | All supersymmetric solutions of minimal supergravity in five dimensions[END_REF][START_REF] Gutowski | General supersymmetric AdS(5) black holes[END_REF]. This indicates that our solution cannot be a 1/2 BPS solution of N = 2 supergravity. In the following we want to discuss what happens to the 1/8 BPS solution (1,2) when mapped to the maximal N = 2 truncation of N = 8 supergravity.

N = 2 truncations and supersymmetry

In order to find a supergravity with eight supercharges in five dimensions, we have to perform a truncation of N = 8 supergravity. The field content of these truncated theories (also called 'magical supergravities') has been discussed for instance in [START_REF] Gunaydin | Exceptional Supergravity Theories and the MAGIC Square[END_REF][START_REF] Gunaydin | The Geometry of N=2 Maxwell-Einstein Supergravity and Jordan Algebras[END_REF]. The N = 2 truncation with the maximal field content (and only vector multiplets) is the magical supergravity related to the Jordan algebra over the quaternions and it admits the global symmetry group SU * (6). It has the same bosonic field content as five-dimensional N = 6 supergravity. As we show in a more detailed work [16], the projection to this N = 2 supergravity in five dimensions corresponds to fixing a complex structure I on T 6 and projecting out some representations of the related SL(3, C). The surviving vector fields of the N = 2 projection contain all gauge fields coming from the eleven-dimensional three-form potential with two legs on T 6 that are (1, 1) with respect to I. Note that I does not have to be related to the complex structure under which dz and dw are holomorphic, as long as the metric given in (1) respects it. If we choose a complex structure I on T 6 such that dz 1 = dx 8 + i dx 5 , dz 2 = dx 6 + i dx 7 and dz 3 = dx 9 + i dx 10

are holomorphic one-forms under I, then the flux given in ( 2) is (1, 1) on T 6 , and we see that our solution indeed gives a solution to N = 2 supergravity. Now let us understand the amount of supersymmetry of the solution in N = 2 supergravity. The complex structure above is different from the complex structure chosen in [START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF], and under the new complex structure the flux F 4 (5) has a piece that is (3, 1) ⊕ (1, 3) and therefore the configuration is not supersymmetric in N = 2 supergravity. More precisely, the projection to N = 2 breaks the N = 8 R-symmetry group USp(8) to USp(6)×SU [START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF], where the latter factor is the R-symmetry of the N = 2 theory. The action of USp(6) on the spinors defines the projection to N = 2. The generator C ≡ 1 2 (Γ 85 -Γ 67 ) [START_REF] Becker | M theory on eight manifolds[END_REF] commutes with the complex structure I, the Cartan generator of SU(2), and hence is a generator of USp [START_REF] Bena | One ring to rule them all ... and in the darkness bind them?[END_REF]. In particular, the requirement Cη = 0 implies

1 2 (1 -Γ 5678 )η = 0 . ( 16
)
This projects out all four Killing spinors of the 1/8 BPS solution, cf. [START_REF] Lust | Generalized non-supersymmetric flux vacua[END_REF]. Hence, when we projected to the N = 2 SU * (6) supergravity, we projected out all supercharges which remain unbroken in the solution [START_REF] Grana | Supersymmetric three form flux perturbations on AdS(5)[END_REF][START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF]. Therefore, the solution is non-BPS in N = 2 supergravity.

We use four-dimensional supersymmetry conventions. For instance, all N =

theories, regardless of dimension, have 8 supercharges.[START_REF] Giddings | Hierarchies from fluxes in string compactifications[END_REF] The fact that a non-supersymmetric solution of an N = 2 or an N = 4 theory can become supersymmetric when embedded in N = 8 has been know for quite a while[START_REF] Ferrara | N=8 non-BPS Attractors, Fixed Scalars and Magic Supergravities[END_REF][START_REF] Khuri | A Nonsupersymmetric dyonic extreme Reissner-Nordstrom black hole[END_REF]. However, in all these examples, there always exists an N = 2 or N = 4 truncation in which the N = 8 solution is supersymmetric. In our example no such truncation exists, and the supersymmetry of the solution cannot be captured in any daughter N = 2 theory.

Acknowledgments

We are very grateful for useful discussions with S. Ferrara, S. Giusto, M. Graña, M. Gunaydin and M. Shigemori. We thank the Aspen Center for Physics and the Centro de Ciencias de Benasque


M Corneli

E Erosheva

A Bayesian approach for clustering and exact finite-sample model selection in longitudinal data mixtures

A Bayesian approach for clustering and exact finite-sample model selection in longitudinal data mixtures

Marco Corneli, Elena

Introduction and background

The question of clustering longitudinal trajectories arises in a number of social sciences concerned with describing human behavior over time (see Erosheva et al., 2014, and references therein). Here, we focus on longitudinal data involving N individuals. Each individual is associated with one (and only one) hidden group. Latent random variable z i labels the group (cluster) of the i-th individual. Assuming that there is a variable of interest that is measured for all individuals, a vector y i ∈ R D keeps track of the measurements for the i-th individual, where D is the number of measurements. In the following, it is assumed that exactly D measurements are taken at the same times across N individuals and y ij denotes the j-th measurement for the i-th individual at time t j , with j ≤ D. Assume the following generative model

y ij = φ(t ij ) T β z i + σ ij , ∀i ≤ N, j ≤ D, (1)
where φ(•) : R → R K is a user defined function (e.g., identity, polynomial, etc.), β 1 , . . . , β Q are parameters in R K with Q the number of mixing components, and σ is the scalar standard deviation of the noise term. The residuals ij are all i.i.d Gaussian N (0, 1) distributed. The following more compact notation will also be employed

y i = Φβ z i + σ i , ( 2
)
where Φ ∈ R D×K is the design matrix, whose j-th row is φ(t ij ) T , y i and i are vectors in R D and i follows a multivariate isotropic Gaussian distribution N (0, σ 2 I D ), with I D denoting the identity matrix of order D. In the following, y 1 , . . . , y N are assumed independent, given z 1 , . . . , z N . Finally, the model can be extended straightforward by assuming that the standard deviation σ also depends on z i . All the results reported in this paper are still valid.

A frequentist approach and some related drawbacks

For a given Q, assuming that z 1 , . . . , z N are i.i.d. random variables and such that P{z i = q} = π q , ∀q ∈ {1, . . . , Q}, the observed data likelihood of the generative model described in the above section is

p(y 1 , . . . , y N |θ) = N i=1 Q q=1 π q g(y i ; Φβ q , σ 2 I D ), (3)
where θ := {β q , π q , σ 2 } q≤Q is the set of the model parameters and g(•; µ, Σ) denotes the probability density function of a multivariate Gaussian distribution with mean µ and covariance matrix Σ. Clearly, this generative model can be seen as a constrained Gaussian mixture model (see Appendix A for more details and a proof that time matters in this model).

A standard approach to estimate the model parameters (see, for instance, [START_REF] Nagin | Group-based modeling of development[END_REF] would: 1. set φ(•) to some polynomial of reasonable order (possibly one!), depending on the the data, 2. maximize the log-likelihood N i=1 log p(y i |θ, Q) with respect to θ to obtain θML . Posterior distributions p(z i |y i , θML , Q) could then be computed for all i and they might used to cluster the observations. Repeating steps 1 and 2 for different Q, the number of clustering components Q might be estimated via some model selection criterion (e.g. AIC or BIC) (see e.g. [START_REF] Nagin | Group-based modeling of development[END_REF][START_REF] Muthén | Growth mixture modeling: Analysis with nongaussian random effects[END_REF].

However, in practice, AIC and BIC are often not able to provide reasonable choice for choosing the number of longitudinal trajectory clusters. It is common that these goodness-of-fit criteria keep improving as the number of groups Q increases until models encounter convergence problems [START_REF] Erosheva | Breaking bad: Two decades of life-course data analysis in criminology, developmental psychology, and beyond[END_REF]. In addition, when the choice of Q is perceived as large because the identified clusters do not represent substantively distinct trajectory patterns, the practical advice is to reduce Q to a lower number that is more meaningful in the applied context (e.g., [START_REF] Nagin | Group-based modeling of development[END_REF].

We now describe some know drawbacks to using AIC/BIC for selecting the number of clusters in standard mixture data analysis. We conjecture that these drawbacks could provide a possible explanation as to why AIB/BIC often fail to produce meaningful number of clusters with longitudinal trajectories.

First, we note that BIC (as well as AIC) is an asymptotic criterion needing a large number of observations to be given and some relevant assumptions to be fulfilled by the data that may not be straightforward in case of longitudinal observations. In addition, BIC was shown to possibly overestimate of the number of clusters in mixture models (see for instance [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Baudry | Combining mixture components for clustering[END_REF]. This point is illustrated with an example in Figure 1. One thousand points were independently sampled from a mixture of Q = 4 Gaussian distributions in dimension two. The actual number of mixing components can be safely recovered by BIC model selection by some software in standard statistical libraries (e.g. mclust on R). However, Q is not necessarily the most meaningful number of clusters. Indeed, in Figure 1 there are 2 disjoint sets of points. In our simulations in Section 4, we show that the over-estimation problem illustrated in Figure 1 might exist also in longitudinal data mixtures.

Second, a more parsimonious alternative to AIC/BIC for model selection in mixture models is the Integrated Classification Likelihood (ICL, [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. See Appendix B for a formal definition of ICL and the proof that ICL is a lower bound of BIC. Although ICL might select a smaller number of components than BIC, it still remains an asymptotic criterion. Moreover, whereas ICL was shown to solve the model selection issue described in Figure 1 (see for instance [START_REF] Baudry | Combining mixture components for clustering[END_REF], it does not solve the overestimation problem for longitudinal data mixtures, as shown in Section 4.

In an attempt to address these issues, this paper focuses on a Bayesian cluster-

-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Figure 1: One thousand data points independently sampled from a mixture of Q = 4 bi-variate Gaussian distributions (identified by the red ellipses), forming two separated groups.

ing framework for the generative model detailed so far allowing us to compute an exact (i.e. a finite sample) version of the ICL by adopting prior distributions on the model parameters and integrating them out. A similar approach was adopted for Gaussian mixture models by [START_REF] Bertoletti | Choosing the number of clusters in a finite mixture model using an exact integrated completed likelihood criterion[END_REF] or by [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]; [START_REF] Corneli | Exact icl maximization in a nonstationary temporal extension of the stochastic block model for dynamic networks[END_REF], in the context of graph data clustering.

We should note that non-parametric Bayesian approaches such as a Bayesian extension of the generative model in Eq. 3 to a Dirichlet process (DP) mixture model (see e.g. [START_REF] Maceachern | Estimating mixture of dirichlet process models[END_REF] would represent an alternative way to analyze longitudinal mixtures. Bayesian non-parametric techniques as well as their extensions to grouped data sharing cluster memberships (Hierarchical DP, [START_REF] Teh | Sharing clusters among related groups: Hierarchical dirichlet processes[END_REF] have been widely used in the context of mixture models and topic models [START_REF] Griffiths | Finding scientific topics[END_REF]. Apart from their flexibility, one of the main advantages of the DP mixture models is their ability to perform clustering and selection of the number of components simultaneously in one shot. In addition, DP mixture models can be used with a reduced computational cost thanks to variational inference [START_REF] Blei | Variational inference for dirichlet process mixtures[END_REF]. However, we do not pursue Bayesian nonparametric methods here because one of the main aims of this paper is to detail a model selection method that could be used on a given clustering separately from the clustering process. In this way, practitioners will have a choice of using their preferred clustering tool and just adopt our model selection criterion or using our complete approach of clustering and model selection in longitudinal data. This objective is difficult (and even unnatural) to be reached via DPMMs.

We organize the remainder of our paper as follows. Section 2 develops our Bayesian framework to clustering longitudinal trajectories. Section 3 develops the associated estimation routine whose computational cost remains reasonable. Section 4 presents a simulation study comparing our method to some state-of-theart alternatives, in order to highlight the main features of our approach. A final discussion, in Section 5, concludes the paper.

A Bayesian perspective

In this section, we detail a Bayesian extension of the model in Eq. ( 2) allowing us to i) cluster the observations y 1 , . . . , y N in Q groups and ii) select Q in a non asymptotic framework. The target probability distribution, that we would like to maximize with respect to the pair (Z, Q) is

p(Z, Q|Y ) = p(Z, Q|Y, θ)p(θ)dθ (4)
where Z := (z 1 , . . . , z N ), Y := (y 1 , . . . , y N ) and the model parameters θ are seen as random variables and integrated out. Note that, from a full Bayesian perspective, the number of clusters Q is also viewed as a random variable in the above equation.

In order to develop an estimation algorithm that is reasonably fast, we choose not to implement MCMC algorithms to simulate the above posterior distribution and present the following alternative strategy that relies on the Bayes rule

p(Z, Q|Y ) = p(Y, Z|Q)p(Q) p(Y ) .
where p(Q) is a prior distribution over Q. Since the denominator does not depend on (Z, Q), it holds that arg max

(Z,Q) p(Z, Q|Y ) = arg max (Z,Q) p(Y, Z|Q)p(Q) = arg max Q arg max Z|Q p(Y, Z|Q) p(Q) .
Note that, because max Z|Q p(Y, Z|Q) does not depend on the functional form of p(Q), the results shown in next sections, only involving max Z|Q p(Y, Z|Q), remain valid for any choice of p(Q). However, to simplify the exposition, we assume that Q is uniform distributed (p(Q) ∝ 1) and the above equation then reduces to arg max

(Z,Q) p(Z, Q|Y ) = arg max (Z,Q) p(Y, Z|Q).
(5) Thus, computing the maximum posterior (MAP) estimates of Z and Q reduces to maximize the complete data integrated log-likelihood.

A closed form complete data integrated log-likelihood

The complete data integrated log-likelihood on the right hand side of Eq. ( 5) is

p(Y, Z|Q) = p(Y, Z|θ, Q)p(θ|Q)dθ, ( 6
)
where we recall that θ := {β q , π q , σ 2 } q . We stress that, in general, this quantity is not tractable, since the integral on the right hand side of the above equality cannot be computed. However, some further assumptions allow us to obtain p(Y, Z|Q) in a closed form, which is essential in order to develop the clustering/model selection strategy detailed in the next section. First, we assume that the prior distribution factorizes over the model parameters, namely

p(β, σ, π) = p(β, σ)p(π),
where, in order to keep the notation uncluttered, we denote β := {β q } q and π := {π q } q . Thus, the integrated log-likelihood in Eq. ( 6) factorizes

p(Y, Z|Q) = p(Y |Z, Q)p(Z|Q). ( 7
)
Then, we assume prior conjugate distributions for β, σ and π. Specifically, conditionally on σ 2 , assume that β 1 , . . . , β Q follow independent Gaussian prior distributions

β q ∼ N 0, σ 2 η q I K , (8)
where η 1 , . . . , η Q are positive parameters and I K is the identity matrix of order K. The β q s are further assumed to be independent from i , for all i. Likewise, assume σ 2 follows an Inverse Gamma prior distribution

σ 2 ∼ IΓ(a, b), ( 9
)
where a, b > 0. Finally, assume π follows a Dirichlet prior distribution

π ∼ Dir(α 1 , . . . , α Q ), (10)
where α q > 0 for all q.

The integration with respect to β is detailed in the next section, since it allows us to highlight a central feature of our Bayesian model. The details of the integration with respect to σ 2 and π are postponed in Appendices C and D.

Integrating with respect to β

By integrating out β q in Eq. (2), we obtain the marginal conditional density of y i

y i |z i , σ 2 ∼ N 0, σ 2 η z i ΦΦ T + I D . ( 11
)
By the law of iterated expectations it follows that Cov(y i , y

j |z i , z j , σ 2 ) = E y i y T j |z i , z j , σ 2 = E Φβ z i β T z j Φ T |z i , z j , σ 2 = σ 2 η z i ΦΦ T 1 {z i =z j } ,
where 1 A (•) denotes the indicator function over a set A. The above equation has an important consequence: after β q is integrated out, the random vectors y i sharing the same cluster are no longer independent, as they were in the frequentist model. Moreover, let us denote by

y (q) := {y i |i ≤ N, z i = q}, ( 12
)
the set of trajectories forming the q-th cluster, whose cardinality is denoted by C q . Since all vectors in y (q) are Gaussian distributed, their joint conditional density can be specified as

Y q |Z, σ 2 ∼ N 0, σ 2 G q , ( 13
)
where Y q ∈ R DCq is a column vector obtained by concatenating all the observations in cluster y (q) and G q ∈ R DCq×DCq is a block matrix. The blocks on the main diagonal are of the form η q ΦΦ T + I D , whereas the blocks outside the main diagonal look like η q ΦΦ T .

Remark. Equation (11) shows that the generative model detailed so far can entirely be expressed in terms of kernels. In other terms, we can directly work with ΦΦ T without previous specifying the feature map Φ(•). By doing that, the generative approach described so far embodies 1. mixtures of linear regression models, with [ΦΦ T ] jl := t j t l , 2. mixtures of polynomial regression models, with [ΦΦ T ] jl := (1 + t j t l ) r being a polynomial kernel of degree r,

mixtures of non-parametric regression models, with [ΦΦ

T ] jl := exp -(t j -t l ) 2 2γ
being a radial basis function (RBF) kernel, for some positive scale parameter γ, 4. any other mixture of regression models depending on the choice of a symmetric, positive definite kernel matrix ΦΦ T .

Inference

The closed form expression of the integrated log-likelihood is reported in Eqs ( 20)-( 21) in Appendices C-D. This quantity depends on the value of the hyper-parameters ι := {η, a, b, α}, that must be set based on the prior knowledge about the the data.

If the number of groups Q is assumed to vary in a range {1, . . . , Q max }, for a given Q, we aim at estimating

Z * Q := arg max Z|Q log p(Y, Z|Q).
Then, the Q leading to the highest value of log p(Y, Z * Q |Q) will be the estimated number of clusters (see Section 2)1 .

The estimation of Z * Q is challenging. The strategy that we propose relies on a greedy maximization of log p(Y, Z|Q). In simple terms, we compute changes in the integrated log-likelihood obtained by switching y i from its current cluster to other proposal clusters, and retain the proposal leading to the highest increase in the loglikelihood. We stress that having a close form expression for log p(Y, Z|Q) is crucial to adopt such strategy. All vectors y 1 , . . . , y N are switched once. Notice that, since Q is given, if y i is alone in its cluster, no movement is allowed. This classification step (CS) is repeated recursively until no further increase of log p(Y, Z|Q) is possible. We stress that the CS described so far is not guaranteed to converge to a global optimum. Indeed, CS is a greedy step. In order to reduce the risk of being trapped into local optima, one should either run the algorithm several times (for each Q), with different initialization, or run the algorithm with a "clever" initialization of Z, for instance the one obtained by k-means clustering. We adopt this latter solution in experiments in Section 4. However, notice that also in case of smart initialisation, it might be useful to run the CS several times, randomizing over the switch order of the observations.

The CS is described in more details in Appendix E, where some original results are reported in order to reduce the computational complexity of this step. Instead, we conclude this section with a couple of important remarks. First, we notice that the algorithm described in this section is reminiscent of the Classification EM algorithm (C-EM, [START_REF] Celeux | A classification em algorithm for clustering and two stochastic versions[END_REF] that replaces the expectation step of the standard EM algorithm with a classification step similar to the one described here. However, while EM algorithms maximize the complete data log-likelihood with respect to the model parameters, here, those parameters are integrated away and the maximization is performed with respect to the cluster assignments alone. Second, it has been shown in the literature that, in latent variable models, likelihood maximization strategies involving a greedy search over the group labels can outperform alternative variational approaches in terms of classification accuracy (see for instance [START_REF] Bertoletti | Choosing the number of clusters in a finite mixture model using an exact integrated completed likelihood criterion[END_REF][START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF][START_REF] Corneli | Exact icl maximization in a nonstationary temporal extension of the stochastic block model for dynamic networks[END_REF].

Experiments

Based on simulation studies, this section compares the Bayesian approach previously described2 with some alternative methods, both in terms of pure clustering (Q fixed) and model selection. The methods considered for comparison are:

1. flexmix3 : in the eponymous R package, it is a function fitting mixtures of generalized linear models. Here, it is used to fit the generative model in Eq. ( 3). Since it is the frequentist version of the generative model introduced in Section 2, we consider flexmix as our major competitor.

2. funHDDC4 : it performs (in R) clustering and model selection of functional longitudinal data based on a latent mixture model which fits high dimensional data in group-specific sub-spaces [START_REF] Bouveyron | Model-based clustering of time series in groupspecific functional subspaces[END_REF][START_REF] Schmutz | Clustering multivariate functional data in group-specific functional subspaces[END_REF].

3. BayesianGaussianMixture5 : from the Python library scikit-learn, it is used in order to compare our approach with a very popular model based clustering technique. Here, longitudinal data are seen as unconstrained multivariate mixtures of Gaussian distributions (despite the discussion in Appendix A), possibly with non diagonal covariance matrices. A prior Dirichlet process, allowing for a potentially infinite number of components is considered and the cluster assignments are estimated via variational inference.

4. mclust6 : in the eponymous R package, it fits mixtures of Gaussian distributions to the data, in a frequentist framework. Definitely a gold standard in model based clustering.

Setup I

Aim. One of the major claims of this paper is that, in longitudinal data mixture analysis, our Bayesian approach might outperform asymptotic model selection cri-teria (BIC and ICL) in some scenarios. This might happen either because BIC tends to overestimate the number of components and/or the number of observations is not large enough to ensure the convergence of the asymptotic criteria to the corresponding posterior probabilities. This section focuses on such scenarios.

Settings. An increasing number of simulated trajectories, around Q = 6 mean signals, is considered. The mean signals are fixed (coloured curves in Figure 2b). Each trajectory is assigned to a group uniformly at random (π q = 1/6 for all q) and the number of simulated trajectories increases from N = 20 to N = 190. Each simulated trajectory is obtained via perturbations around the corresponding signal. We set D = 12 points per trajectory, not equally spaced. Figure 2a shows 40 (linearly interpolated) sampled trajectories over the time horizon [5,60]. The same number of sampled trajectories and the underlying signals are shown in Figure 2b.

For each value of N , we simulated 10 independent data sets (around the same 6 signals) and each model was fitted on each data set testing a value of Q from 1 to 8. In order to avoid convergence to local maxima, each algorithm run 10 times, one for each data set and for each value of Q, provided with a k-means initialization of Z. The estimates leading to the highest value of the objective function were finally retained. Our model and flexmix were equipped with a polynomial kernel of order 4, the choice of the hyper-parameters for our model is discussed later in this section. Finally, we tested two different families of random perturbations, Gaussian and Gumbel, respectively, both centered and with standard deviation σ = 0.5. The results for model selection are reported in Figure 3. Note that the legend reports our method as "Bayes" and BayasianGaussianMixture as "NPBGMM". For each value of N we report the mean value of the selected Q by each model (over 10 runs). For all models, except Bayes and NPBGMM, the reported Q was selected via ICL, whose calculation was supported by different packages. Interestingly, the differences with respect to BIC were absolutely negligible, knowing that Q BIC ≥ Q ICL . In case of flexmix we always found that

Q BIC = Q ICL .
Discussion. Although in a context of mixture density estimation it is of course desirable to be able to select Q = 6 mixture densities, in a clustering perspective it is fully reasonable to desire a model selection criterion selecting three main groups of trajectories, especially for small N , as it can be understood when looking at Figure 2a (three groups of well separated trajectories and not six). This point was raised and discussed in more general terms in Section 1.1. Thus, when looking for "clusters", both flexmix and mclust might be overestimating the number of groups for small values of N both in case of Gaussian and Gumbel perturbations (in the Gaussian case the overestimation is more striking). Not surprisingly, a Gaussian mixture model is unfit to model longitudinal data with Gumbel residuals and this 3a) and Gumbel residuals (Figure 3b). mclust, in Figure 3b, even for large N . flexmix always selects more groups than our model and, interestingly, BIC and ICL produce the very same selected Q. This suggests that also the over penalization induced by ICL (see Appendix B) fails to produce a more conservative model selection criterion than BIC. In our model (Bayes), the role of the hyper-parameters is crucial. We will come onto this point in a while. The other competitor approaches (funHDDC and NPBGMM) seem to suffer less overestimation issues for small N . However, the low values of Q estimated by funHDDC are certainly dues to a poor clustering more then to a more penalizing model selection design. This can be seen in Figure 4. It reports the mean Adjusted Rand indexes (ARIs, [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF]) (over 10 runs) obtained by each method when fixing Q = 6. The first four values for funHDDC are not reported since, for Q = 6, the log-likelihood went to minus infinity and the model produced NA. Starting from N > 50 the ARIs can always be computed but they are lower then those obtained by other methods. We think that funHDDC might need more observations per trajectory to work properly. On the opposite, our clustering method outperforms the competitors most of the time. We stress that, for values of N around 50 trajectories, the ARIs are almost one, meaning that the model can perfectly recover the Q = 6 mixture components. Still, if it is allowed to choose the number of clusters, it selects Q = 3. NPBGMM has good clustering performances for high N . The fact that it needs more observations to recover the true clustering is consistent with the discussion in Appendix A. Note, however, that NPBGMM is not capable of performing clustering with a fixed number of groups. This might penalize it for smaller values of N when it selects Q < 6. Several values of the hyper-parameters were tested both for the mean precision and the weight concentration priors in NPBGMM but the best results were always obtained with the default parameters provided by the function.

The experiments reported so far were performed with the simple settings showed in Table 1. To simplify the hyper-parameter choice we assumed a symmetric Dirchlet distribution in Eq. 10 and fixed η q = η for all q. In this setup, the value of

Parameter Value η 1.0 a 1.0 b 1.0 α 10
Table 1: Values of the prior hyper-parameters.

α was not relevant when the algorithm was provided with a "smarter" k-means initialization. On the contrary, we observed that with a random initialization of Z, the choice of the Dirichlet parameter α had an impact on the final clustering (not reported). In more detail, we noticed that the Classification Step detailed in Appendix E tended to switch too many observations at once, thus emptying some clusters. A similar remark was previously made by [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]; [START_REF] Corneli | Exact icl maximization in a nonstationary temporal extension of the stochastic block model for dynamic networks[END_REF] in the context of graph clustering. However, such a tendency to over-switching can be resisted by picking a large value for the initial α (see Ch.4 of Fruhwirth-Schnatter et al., 2019, for a detailed discussion about the Dirichlet prior hyper-parameters in mixture analysis.). A value of α = 10 or even higher can make the job in this framework. We noticed that the value of the hyper-parameters a and b did not impact on the final result. Instead, the hyper-parameter η has an important role. We performed again model selection in the same settings described at the beginning of this section, with Gaussian residuals, but testing several values of η. As it can be observed in Figure 5, our model selection criterion could even be more parsimonious for small values of η, thus needing higher values of N to select highers Q's. On the opposite, higher values of η allow us to select a higher number of clustering components for smaller values of N . However, the "convergence" of our model selection criterion toward its asymptotic counterpart (ICL) is never fulfilled when N < 30: the selected Q is still 3.

Setup II

Aim. As stated in Section 2.1, the Bayesian formulation of the generative model in Eq. 2 allows us to directly fit mixtures of linear, polynomial or non-parametric regression models to the data. A non-parametric scenario is considered in this section. with and without the underlying signals, respectively, for a value of σ = 0.5. Unlike the previous setup, here trajectories are not assigned to clusters in equal proportions, but according to the following scheme: 4 trajectories around s 1 (•), 8 trajectories around s 2 (•), 12 trajectories around s 3 (•) and, finally, 16 trajectories around s 4 (•). Note that, such a setting is more challenging with respect to the equally balanced case, since most of the trajectories are affected to s 3 (•) and s 4 (•) who are the nearest signals. As long as σ increases, it becomes impossible to separate these two clusters. For different values of σ, 10 independent samples of N trajectories each simulated. Our model was equipped with a RBF kernel

Settings

[ΦΦ T ] jl = exp - (t j -t l ) 2 2γ
with γ = 0.04, selected by cross validation. The hyper-parameters are set as in the previous setting, except for η, whose role is discussed at the end of this section. We stress that, once fitted to the data, our approach implicitly provides estimates for the mean signals in Figure 6b. Indeed, the estimated parameters (as well as the estimated Z) can be used to compute the mean and variance of the (Gaussian) predictive distribution p(y i (t)|y i , ẑi , σ2 , ι), where t is a new time point corresponding to a potential measurement for the i-th individual. In Figure 7 the estimated mean signals (the mean values of the predictive distributions over time) are plotted together with the original simulated data clustered by our algorithm. Dashed lines delimit 95% empirical confidence regions, obtained via sampling according to the estimated predictive distributions. All models were fitted 5 independent times on each simulated dataset, provided with a k-means initialisation. The run leading to the highest value of the objective function was finally retained. The results for clustering and model selection are reported in Figure 8.

Discussion.

As it can clearly be seen in Figure 8a, our algorithm and flexmix outperform other approaches. In case of mclust and NPBGMM this is not surprising: in this setup, the number of points for trajectory (D = 100) is higher than the number of observations (N = 40). Thus, both these approaches work in very high dimension (see also Appendix A). For NPBGMM we tried several options for the maximum number of iterations, the mean precision prior and the weight concentration prior. The best results are reported in Figure 8 and correspond to the default values of the hyper-parameters provided by the Python function.

We also report that, since flexmix does not support mixture of non-parametric regressions, we needed to perform progressively higher order polynomial regressions (up to order 6) in order to obtain the best results showed in Figure 8a. Still, our method slightly outperforms flexmix. For completeness, results of model selection are shown in Figure 8b. As in the previous section, the model selection via ICL criterion was reported for flexmix, Mclust and funHDDC: the differences between ICL and BIC were negligible. Being the number of mixing components smaller in this setup, model selection was less challenging. Basically all methods tend to select Q = 4 groups for the smallest value of σ = 0.5. Then, trajectories around s 3 (•) and s 4 (•) are less and less separated and models tend to select Q = 3 groups and finally Q = 2 (our approach still selects Q = 3 for some samples, for σ = 1.1). Finally, Figure 9 reports the effect of the hyper-parameter η on the selection of Q. Similarly to the previous experiment, smaller values of η can be used to have a more parsimonious model selection. Notice that the value η = 1.2 is the one used to obtain the curve in Figure 8b. Moreover, as in the previous scenario, different values of η affect the model selection but not the clustering: the four values of η always lead to the Bayes ARIs reported in Figure 8a (not reported).

Discussion and perspectives

The experiments in the previous section suggest that our approach can be used as a viable alternative as a pure clustering and model selection tool in applications. Indeed, it outperforms the state-of-the art when performing clustering (with a fixed Q) of longitudinal data and provides the user with a model selection criterion which is by far more flexible than BIC and ICL. Our criterion can be very parsimonious depending on the choice of the hyper-parameter η. As usual in Bayesian statistics, this parameter must be set by the statistician based on the prior knowledge he/she has about the phenomenon under observation. Smaller prior variances should be preferred when the aim is to detect very well separated clusters of trajectories, higher prior variances when one is interested in the number of Gaussian mixing components needed to fit the data distribution. The higher the number of observations the less crucial the choice.

Although both the proposed Bayesian modeling (based on conjugated prior distributions on the model parameters) and the inference strategy (based on a greedy search over trajectory labels) can easily be extended to multivariate regression models not involving time or to longitudinal data mixtures whose measurements are not taken at the same times, we recall that, on those cases, the results presented in Appendix E would no longer be valid. Future researches might focus on strategies to keep under control the complexity of the estimation procedure under less restrictive assumptions. Also, it would be interesting to extend the generative model considered here in order to account for individual effects and/or autoregressive components.

A Relation with multivariate Gaussian mixture models

To simplify the exposition, let us consider here the simple case where φ(•) is the identity function. Then, Eq. ( 2) reduces to

y i = tβ z i + σ i ,
where β z i is a scalar and t = (t 1 , . . . , t D ) T ∈ R D . Thus

y i |z i ∼ N tβ z i , σ 2 I D .
First, we stress that, since a Gaussian distributions is fully characterized by its mean and variance matrix, a permutation of the times inside t induces a different probability distribution on y i . Thus, unless

β 1 = • • • = β Q = 0,
p(y i |π, β, σ 2 ) = Q q=1 π q g(y i ; tβ q , σ 2 I D ), (14)
where π := (π 1 , . . . , π Q ) T are the mixing proportions, g(•; µ, Σ) denotes the pdf of a multivariate Gaussian distribution with mean µ and variance Σ and z i was integrated out. Since, y i follows a multivariate Gaussian mixture distribution, would it be reasonable to attack the problem of clustering y 1 , . . . , y N via standard unconstrained Gaussian mixtures machinery (thus ignoring the data dependence on time)? Not really. It is worth noticing that the generative model induced by the above equation has 2Q free parameters to estimate. On the contrary, if also we saw that model as a not constrained Gaussian mixture model with spherical shared variance (the most parsimonious version), namely

p(y i |π, β, σ 2 ) = Q q=1 π q g(y i ; µ q , σ 2 I D ),
with µ q ∈ R D , unknown, the above equation would count Q(D+1) free parameters to estimate. Thus, as long as D ≥ 2 the longitudinal formulation is more parsimonious. In particular, in a high dimensional framework, the difference between the two approaches would be dramatic.

B ICL and BIC

We consider a standard mixture model where the N observed variables are denoted by X := (x 1 , . . . , x N ) and the corresponding latent variables by Z = (z 1 , . . . , z N ).

The unknown number of mixing components is K. The (asymptotic) BIC and ICL criteria are defined as follows

BIC K := max θ log p(X|θ, K) - ν(K) 2 log N (15)
and

ICL K := max θ log p(X, Z|θ, K) - ν(K) 2 log N, ( 16
)
where θ denotes the set of the model parameters, ν(K) is the number of model parameters and log p(•) denotes the log density of the observations.

Remark. The ICL criterion in Eq. ( 16) is an approximation of the marginal loglikelihood log p(X, Z|K) = log

θ p(X, Z|θ, K)p(θ|K)dθ, ( 17
)
where the model parameters are integrated out [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. If the prior distribution p(θ|K) is conjugated, the quantity on the left hand side of the above equation can be computed explicitly. We sometimes call it exact ICL.

The following notations are adopted θ = arg max

θ log p(X|θ, K), ( 18
) θ = arg max θ log p(X, Z|θ, K) (19)
and we stress that, in general, θ = θ.

The following Proposition formally shows that ICL K in Eq. ( 16) is a lower bound of BIC K . This result was mentioned in [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Baudry | Combining mixture components for clustering[END_REF] but not formally proven.

Proposition 1.

ICL K ≤ BIC K .

Proof.

ICL K -BIC K = log p(X, Z|θ) -log p(X| θ) = log p(X, Z|θ) p(X| θ) = log p(X, Z|θ)p(X|θ) p(X|θ)p(X| θ) log p(Z|X, θ) + log p(X|θ) p(X| θ) ≤ 0,
where the dependence on K was omitted for simplicity and the last inequality comes from the discrete nature of the random variables z 1 , . . . , z N and the definitions of θ and θ.

C Integrating with respect to σ 2

In Section 2.1.1 we saw that the marginal conditional density p

(Y |Z, σ 2 , Q) is p(Y |Z, σ 2 , Q) = Q q=1 p(Y q |Z, σ 2 ) = 1 (2πσ 2 ) DN 2 Q q=1 det(G q ) exp   - 1 2σ 2 Q q=1 Y T q (G q ) -1 Y q   ,
where G q ∈ R DCq×DCq is the block matrix introduced in Eq. ( 13). Since, Eq. ( 9) states that

p(σ 2 |a, b) = b a Γ(a) 1 σ 2 a-1 exp - b σ 2 1(σ 2 ) ]0,∞[ ,
when looking at the joint conditional density p(Y, σ 2 |Z, Q) as a function of σ 2 we recognize the pdf of an Inverse Gamma distribution IΓ a

+ DN 2 , b + 1 2 Q q=1 Y T q G -1 q Y q . Therefore, σ 2 can be integrated out to obtain p(Y |Z, Q) = 1 (2π) N D 2 q det(G q ) b a Γ(a) Γ DN 2 + a b + 1 2 q Y T q (G -1 q )Y q DN 2 +a .
(20)

D Integrating with respect to π

The second integral on the right hand side of Eq. ( 7) can be computed in a similar fashion. Due to Eq. ( 10), the posterior distribution p(π|Z, Q) is still a Dirichlet. It can easily be seen that

p(Z|Q) = Γ Q q=1 α q Q q=1 Γ(α q ) Q q=1 Γ(C q + α q ) Γ N + Q q=1 α q , ( 21
)
where we recall that C q is the number of trajectories in cluster q.

E Classification step in depth

Consider the observation y i ∈ R D and assume it currently belongs to the q-th cluster, namely z i = q. The change in the integrated log-likelihood due to a switch of y i from the q-th cluster to the l-th cluster can be computed as

∆ i:q→l ll : = log p(Y, Z a |Q) -log p(Y, Z b |Q) = log p(Y |Z a , Q) -log p(Y |Z b , Q) + log p(Z a |Q) -log p(Z a |Q), (22)
where Z a and Z b denote the configurations after and before the switch, respectively. As it can be seen by looking at Eqs. ( 20)-( 21), the calculation of ∆ i:q→l ll basically boils down to compute i) the determinant of G q and ii) the quadratic term Y T q G -1 q Y q , for all q ≤ Q. We report in the following some results allowing us to speed up the calculation of these two terms.

First term. First, recall that

G q =       B q A q . . . A q A q B q . . . A q A q A q . . . A q A q . . . . . . B q       (23)
where A q = η q ΦΦ T ∈ R D×D and B q = A q + I D . Since the size of G q changes whenever an observation is switched from one cluster to another, we need a fast way to compute det(G q )7 . Theorem 1 in [START_REF] Silvester | Determinants of block matrices[END_REF], whose precise formulation is reported in Appendix F, can help us. This theorem basically allows us to compute det(G q ) in two steps:

1. a first intermediate determinant (ID q ) is computed as if A q , B q in Eq. ( 23) where numbers

ID q := det(G q ) ∈ R D×D , (24)
where the over line is used to differentiate this determinant from the real one, that we are actually trying to compute. Then, 2. since ID q is itself a matrix, det(G q ) is computed as det(G q ) = det(ID q ). ( 25)

According to Theorem 1 in [START_REF] Silvester | Determinants of block matrices[END_REF], the above equality holds as long as all the possible matrix products between two blocks in G q are commutative. This condition is fulfilled as stated in the following Proposition 2. The product between each pair of blocks of G q is commutative.

Proof. For simplicity, in this proof the subscript q is removed from A q and B q , since not needed. Clearly the products AA and BB are commutative. Moreover,

AB = A(A + I D ) = AA + A = (A + I D )A = BA
and the proposition is proven.

Now we can state the following Theorem

Theorem 1. The determinant of G q can be computed in O(D) as

det(G q ) = D j=1 (1 + C q λ (q) j ), ( 26
)
where λ (q) 1 , . . . , λ

D are the eigenvalues of A q and C q is the cardinality of the q-th cluster (see Eq. ( 12)).

Proof. The proof of this Theorem relies on Lemma 1 in Appendix F, stating that ID q in Eq. ( 24) is

ID q = I D + C q A q ,
which is a matrix in R D×D . Now, since A q is symmetric it admits a diagonal representation A q = Q q Λ q Q T q , where Λ q ∈ R D×D is a diagonal matrix whose non null entries are the eigenvalues of A q and Q q is an orthogonal matrix whose columns are the corresponding eigenvectors. Thus

det(ID q ) = det(I D + C q Q q Λ q Q T q ) = det Q q (I D + C q Λ q )Q T q = det (I D + C q Λ q ) = D j=1 (1 + C q λ (q) j ).

Second term. The quadratic form

Y T q G -1 q Y q (27)
is now considered. Adopting the notation in Eq. ( 23), Theorem 3 in Appendix F states that the inverse matrix G -1 q is also a diagonal block matrix

G -1 q =       V q W q . . . W q W q V q . . . W q W q W q . . . W q W q . . . . . . V q      
where

W q : = -(I D + C q A q ) -1 A q , V q : = W q + I D . ( 28
)
Thus, to compute the quadratic form in Eq. ( 27) it is not required to invert the whole G q but only the matrix (I D + C q A q ). Notice that the computational cost of this operation is independent of the number of observations N8 . Moreover

Y T q G -1 q Y q = Cq i=1 y (q) i T V q y (q) i + Cq i=1 Cq j=1 j =i y (q) i T W q y (q) j = Cq i=1 y (q) i T I d y (q) j + Cq i=1 Cq j=1 y (q) i T W q y (q) j = Y q 2 2 +   Cq i=1 y (q) i T   W q   Cq j=1 y (q) j   ,
where • 2 denotes the Euclidean norm and we denoted by y (q) i ∈ R D the i-th column vector in cluster y (q) . Since the sum of all observations in cluster y (q) (namely Cq j=1 y (q) j ) can be pre-computed before and after each switch and W q does not depend on i, the last term on the right hand side of the above equation can be computed in O(D 2 ) that can be sensibly smaller than O(N 2 D 2 ) needed for a direct calculation of Y T q G -1 q Y q .

F Linear algebra results

Theorem 2 [START_REF] Silvester | Determinants of block matrices[END_REF]). Let R be a commutative subring of F n×n , where F is a field and F n×n denotes the set of matrices n × n over F . Let M ∈ R m×m . Then det

F M = det F det R (M ) .
Lemma 1. Consider a R N ×N square matrix A such that

A ij =    a + if i = j a otherwise
where a, are two real constants. Then det(A) = N -1 ( + N a).

(29)

Proof. We proceed by recurrence. For N = 1, A = (a + ) and Eq. ( 29) is verified. Now, let us assume that Eq. ( 29) holds for all i ≤ N . The case where N is an even number is considered at first. Thus det M N +1 = det where the sign of the last term on the right hand side (r.h.s.) of the equality is due to the assumption of an even N . Thanks to the inductive assumption, the first term on the r.h.s of the equality is (a + ) det(M N ) = (a + ) N -1 ( + N a)

= a N -1 ( + N a) + N ( + N a)

= N ( + (N + 1)a) + a 2 N N -1 .

Similarly, the inductive assumption can bu used to replace det(M N -i-1 ) = N -i-2 ( -(N -i -1)a) in the second term on the r.h.s. of Eq. ( 30). Thus, by developing the sum over i we obtain det(M N -1 ) = N ( + (N + 1)a)

+ a 2 N N -1 - a 2 N N -1 + • • • + a N N ! - a N N ! = N ( + (N + 1)a).
The case where N is odd is analogous and the lemma is proven.

Lemma 2. Consider a R N ×N square matrix A such that

A ij =    a + if i = j a otherwise
where a, are two real constants. Then the inverse A -1 is Proof. We need to prove that M M -1 is a block matrix whose non-diagonal blocks are matrices in R D×D having zero everywhere (henceforth denoted by 0 R D×D ) and whose diagonal blocks are I D . We observe that (I D +N A) -1 = det(M )

A -1 ij =    +(N -
-1

, where det(•) is defined in Eq. ( 24) and (I D + N A) is invertible thanks to the assumption of invertible M combined with Lemma 1 and Theorem 1 in [START_REF] Silvester | Determinants of block matrices[END_REF]. Now, the following notation is introduced

D M := (I D + N A) -1
to simplify the exposition. Moreover, with a slight abuse of notation we denote by (M M -1 ) ij the block (and not the real entry!) at position (i, j) in M M -1 . Then, for j = i

Figure 2 :

 2 Figure 2: Figure 2a shows N = 40 random trajectories distributed around Q = 6 signals (not showed). Figure 2b shows another sample of N = 40 trajectories as well as the underlying signals. In both cases, the random trajectories are obtained by means of centered Gaussian perturbations around the mean signals (σ = 0.05).

Figure 3 :

 3 Figure3: The mean Q selected by each model over ten runs in case of Gaussian residuals (Figure3a) and Gumbel residuals (Figure3b).

Figure 4 :

 4 Figure 4: Adjusted Rand indexes obtained by the clustering algorithms for Q = 6. Gaussian residuals in Figure 4a, Gumbel residuals in Figure 4b.

Figure 5 :

 5 Figure 5: Selection of Q performed by our model, with Gaussian perturbations, increasing number of trajectories and for different values of the prior hyper-parameter η.

Figure 6 :

 6 Figure 6: Figure 6a shows N = 40 random trajectories distributed around Q = 4 signals (not showed). Figure 6b show the same sample with the underlying signals. The random trajectories are obtained by means of centered Gumbel perturbations around the mean signals (σ = 0.5).

Figure 7 :

 7 Figure 7: Trajectories clustered via our algorithm (σ = .5). The solid dark lines are the means of the predictive distributions on time. The dashed lines delimit 95% empirical confidence regions.

Figure 8 :Figure 9 :

 89 Figure 8: Adjusted Rand indexes obtained by the clustering algorithms for Q = 4 and selected number of clusters.

  + ) det(M N ) -a 2 N det(M N -1 ) + a 2 N (N -1) det

  It suffices to verify that, given A -1 in Eq. (31), AA -1 = I N .Theorem 3. Consider an invertible square block matrix M ∈ R DN ×DN such that diagonal blocks A ∈ R D×D are symmetric matrices and the diagonal blocks are B = A + I D . Then, the inverse matrix M -1 is still a block matrixM the non-diagonal blocks W are W = -(I D + N A) -1 A, (33)and the diagonal blocks V areV = (I D + N A) -1 (I D + (N -1)A) = W + I D .(34)

(

  M M -1 ) ij = -BD M A + AD M (I D + (N -1)A) -(N -2)AD M A = -D M A + AD M = 0 R D×D ,where the last equality comes fromA(I D + N A) = (I D + N A)A ⇒ A = (I D + N A)A(I D + N A) -1 ⇒ (I D + N A) -1 A = A(I D + N A) -1 . (35) Similarly, for i = j (M M -1 ) ii = BD M (I d + (N -1)A) -(N -1)AD M Aand some easy calculations show that the above quantity is equal to I D .

  time matters. Second, recalling that y 1 , . . . y N are assumed independent, the above equation defines a constrained version of a multivariate Gaussian mixture model, since

More generally, when Q in not (a priori) uniformly distributed, the selected Q is the one maximizing log p(Y, Z * Q |Q) + log p(Q), where the last term is independent of Z * Q .

The Python code implementing the clustering/model selection method described in this paper is publicly available at https://github.com/marcogenni/BayesianLongitudinal.

https://cran.r-project.org/web/packages/flexmix/index.html.

https://cran.r-project.org/web/packages/funHDDC/index.html.

https://scikit-learn.org/stable/modules/mixture.html# variational-bayesian-gaussian-mixture.

https://cran.r-project.org/web/packages/mclust/index.html.

For instance, the cost of computing det(G q ) via an LU decomposition is O(D 3 N 3 ) and this approach is used by most linear algebra libraries.

For instance, relying on the Gauss-Jordan elimination, the computational cost of the inversion would be O(D 3 ), which is smaller than O(N D 3 ) required to invert G q .