B Delinchant
email: benoit.delinchant@leg.ensieg.inpg.fr

Interval Arithmetic for the Design of Electrical Engineering Systems

Keywords: Interval arithmetic, formulas stability, tolerance design, global optimization

Electrical engineering system design requires more and more light models like analytical ones. Using interval arithmetic calculation on algebraic formulas brings several benefits like stability analysis, tolerance design or global optimization. Software based on these techniques have been developed to bring these functionalities to designer hands, especially for preliminary stages of design.

I. INTRODUCTION

Analytical models in the field of electrical engineering have been developed since years. These models are more and more used especially for system design needs, multiphysics needs, quick exploration of inverse problem solution space, etc. Based on algebraic formulas, interval arithmetic computes interval instead of real number. For instance, the equation y=a. Intervals are important information for preliminary design studies as it will be shown in this paper. As detailed in [START_REF] Saxena | The use of interval mathematics in electromagnetic design[END_REF], interval analysis cannot be used for accurate modelling, but it is very useful for high level design with behavioural model.

II SOFTWARE ARCHITECTURE

Our design framework is based on software component architecture. A component represents a part of the system to design or all the system thanks to composition software [START_REF] Delinchant | A component-based framework for the composition of simulation software modeling electrical systems[END_REF]. Two kind of tools are then available, generators who create components (Fig. 2-a) and tools which bring services to components (Fig. 2-c) (like optimization framework). An interval arithmetic generator has been developed using BIAS library [START_REF] Popova | A Generalization of BIAS Specification[END_REF] (Basic Interval Arithmetic Subroutines). BIAS uses code differentiation to produce interval gradients (Fig. 2

-b), necessary for our optimization needs. The generator translates algebraic equations of the device model to a component dedicated to optimization (which computes equations propagating intervals as well as formally exact sensitivities).

A calculator is available to compute a model instance or to plot parameters. An optimizer is also available with several algorithms like SQP (sequential quadratic programming), genetic, hybrid, global, or mixed discret-continous. Three uses of this generator tool for design needs are presented in the following three parts.

III. MODEL STABILITY

First of all, out generator is helpful to detect computation rounding errors. Indeed, computer rounds number when floating point values (ANSI/IEEE Standard 754-1985, Standard for Binary Floating Point Arithmetic) are propagated through the model. When modeling a system with huge different value orders (like electromagnetic compatibility study in electromechanical converters), there can be instable formulae. Thanks to interval analysis, which adapts computer rounding method (upper for superior bound and lower for inferior bound), instability can be observed. Since it is possible to stabilize formulae by reordering numbers in the equations, interval computing is an interesting tool to debug models. For example 1 , the classical formula (1) due to Heron of Alexandria which computes the area of a needle like triangle is unstable. It can be stabilized by writing equation [START_REF] Delinchant | A component-based framework for the composition of simulation software modeling electrical systems[END_REF].

( ) )( )( ( 2 1 1 c b a s with c s b s a s s + + = - - - = ∆ (1) )) ( ))( ( ))( ( ))( ( ( 4 1 2 c b a b a c b a c c b a - + - + - - + + = ∆ (2)
a=10 14 , b=10 14 and c=1.1 gives ∆ 1 with ±1.4 % interval radius error, against ±7.10 -14 % for ∆ 2 . It highlights possible instability for engineering models, but also the fact that variables order writing is important, as well as parenthesis! Full paper will show on engineering models some practical writing formula hints to decrease overestimation of interval.

IV. TOLERANCE DESIGN

Input parameters of a model are not accurate but attached to a precision then defined as intervals. Outputs are also intervals thought the intervals propagation mechanism, and can be constrained. As an example, a transformer optimization is presented here. Analytical model [START_REF] Poloujadoff | A Procedure for illustrating the Effect of Variation of Parameters on Optimal Transformer Design[END_REF] was turned by our generator into an autonomous code which was plugged into the optimization tools and connected to an optimization algorithm (sequential quadratic programming). Several inputs were fixed with a precision (left part of Table 1) and others were optimizable (right part of Table 1). Note that winding turns was also defined as an interval of radius 1 to express its integer type. Several outputs parameters were constrained (total transformer cost, leakage reactance and total length) as show optimization results in Fig 4. Optimization iterations total length (m) Fig. 4 -interval results (inf/sup) of optimization V. GLOBAL OPTIMISATION Third use of interval arithmetic is global optimization [START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF]. A simple algorithm using interval arithmetic was implemented into our framework. It is based on branch and bound technique using gradient sign of objective interval for its upper and lower bound. Efficient algorithms are available (GlobSol2 …) and many literatures are also available3 . To illustrate rigorous global optimization of interval based algorithm, the Levy function [START_REF] Levy | Topics in Global Optimization[END_REF] as been defined as a global minimization test (780 local minima and 1 global minimum). Interval algorithm converged successfully contrary to another "global" optimizer (DIRECT [START_REF] Gablonsky | An implementation of the DIRECT algorithm[END_REF]) available in our framework, which is efficient but not rigorous as interval arithmetic.

VI. CONCLUSION

This paper shows how an interval arithmetic model generator (detailed in full paper), producing calculus component as well as sensitivity computation, can be used for the design of engineering systems. It is especially used for formula stablility debugging, toleranced design, and global optimization (more detailed and illustrated in full paper).

Fig. 1 -

 1 Fig. 1 -Example of interval arithmetic on y=a.x+b equation.

Fig. 2 -

 2 (a) Generators create (b) components / (c) Services use (b) components.

Fig. 3 -

 3 Fig. 3 -Needle like triangle.

)

TABLE 1 :

 1 FIXED INPUTS WITH THEIR INCERTITUDE AND OPTIMIZABLE INPUTS WITH THEIR BOUNDS

		value precision	unit		min	max	result	unit
	frequency	50	5%	Hz	induction	0.5	1.9	1,69	T
	apparent power	40	5%	MVA	current density	500	4500	2020	kA/mm²
	supply tension	60	5%	kV	winding turns (±1)	100	800	573	-
					transformer height	1	10	3.43	m

http://www.cs.berkeley.edu/~wkahan/Triangle.pdf

http://interval.louisiana.edu/GLOBSOL/what_is.html

http://www.mat.univie.ac.at/~neum/glopt.html


Z Badreddine

L Rifford

Measure contraction properties for two-step analytic sub-Riemannian structures and Lipschitz Carnot groups

We prove that two-step analytic sub-Riemannian structures on a compact analytic manifold equipped with a smooth measure and Lipschitz Carnot groups satisfy measure contraction properties.

Introduction

The aim of this paper is to provide new examples of sub-Riemannian structures satisfying measure contraction properties. Let M be a smooth manifold of dimension n ≥ 3 equipped with a sub-Riemannian structure (∆, g) of rank m < n, whose geodesic distance d SR is supposed to be complete. We refer the reader to Appendix A for the notations used throughout the paper. As in the previous paper of the second author on the same subject [START_REF] Rifford | Ricci curvatures in Carnot groups[END_REF], we restrict our attention to the notion of measure contraction properties in metric measured spaces with negligeable cut loci (if A ⊂ M is a Borel set then L n (A) = 0 means that A has vanishing n-dimensional Lebesgue measures in charts): Definition 1. We say that the sub-Riemannian structure (∆, g) on M has negligeable cut loci if for every x ∈ M , there is a measurable set C(x) ⊂ M with L n (C(x)) = 0, 0 Both authors are supported by the ANR project SRGI "Sub-Riemannian Geometry and Interactions", ANR-15-CE40-0018.

* Université Côte d'Azur, CNRS, Inria, Labo. J.-A. Dieudonné, UMR CNRS 7351, Parc Valrose 06108 Nice, Cedex 2, France (Zeinab.Badreddine@unice.fr)

† Université Côte d'Azur, CNRS, Inria, Labo. J.-A. Dieudonné, UMR CNRS 7351, Parc Valrose 06108 Nice, Cedex 2, France (Ludovic.Rifford@math.cnrs.fr) 1 and a measurable map γ x : (M \ C(x)) × [0, 1] -→ M such that for every y ∈ M \ C(x) the curve s ∈ [0, 1] -→ γ x (s, y)

is the unique minimizing geodesic from x to y.

Measure contraction properties consists in comparing the contraction of volumes along minimizing geodesics from a given point with what happens in classical model spaces of Riemannian geometry. We recall that for every K ∈ R, the comparison function s

K : [0, +∞) → [0, +∞) (s K : [0, π/ √ K) → [0, +∞) if K > 0) is defined by s K (t) :=        sin( √ Kt) √ K if K > 0 t if K = 0 sinh( √ -Kt) √ -K if K < 0.
In our setting, the following definition is equivalent to the notion of measure contraction property introduced by Ohta in [START_REF] Ohta | On the measure contraction property of metric measure spaces[END_REF] for more general measured metric spaces (see also [START_REF] Sturm | On the geometry of metric measure spaces[END_REF]).

Definition 2. Let (∆, g) be a sub-Riemannian structure on M with negligeable cut loci, µ a measure absolutely continuous with respect to L n and K ∈ R, N > 1 be fixed. We say that (∆, g) equipped with µ satisfies MCP(K, N ) if for every x ∈ M and every measurable set

A ⊂ M \ C(x) (provided that A ⊂ B SR (x, π N -1/K) if K > 0) with 0 < µ(A) < ∞, µ (A s ) ≥ A s s K sd SR (x, z)/ √ N -1 s K d SR (x, z)/ √ N -1 N -1 dµ(z) ∀s ∈ [0, 1],
where A s is the s-interpolation of A from x defined by

A s := γ x (s, y) | y ∈ A \ C(x) ∀s ∈ [0, 1].
In particular, (∆, g) equipped with µ satisfies MCP(0, N ) if for every x ∈ M and every measurable set A ⊂ M \ C(x) with 0 < µ(A) < ∞,

µ (A s ) ≥ s N µ(A) ∀s ∈ [0, 1].
To our knowledge, the first study of measure contraction properties in the sub-Riemannian setting has been performed by Juillet in his thesis. In [START_REF] Juillet | Geometric inequalities and generalized Ricci bounds in the Heisenberg group[END_REF], Juillet proved that the n-th Heisenberg group H n (with n ≥ 1) equipped with its sub-Riemannian distance and the Lebesgue measure L 2n+1 (in this case the ambiant space is R 2n+1 ) satisfies MCP(0, 2n + 3). This result is sharp for two reasons. First, Juillet proved that H n does not satisfy any other stronger notion of "Ricci curvature bounded from below" in metric measured spaces such as for example the so-called curvature dimension property (see [START_REF] Lott | Weak curvature conditions and functional inequalities[END_REF][START_REF] Sturm | On the geometry of metric measure spaces[END_REF][START_REF] Sturm | On the geometry of metric measure spaces[END_REF][START_REF] Villani | Optimal transport, Old and New[END_REF]). Secondly, Juillet showed that 2n + 3 is the optimal dimension for which H n satisfies MCP(0, N ), there is no N < 2n + 3 such that H n (equipped with d SR and L 2n+1 ) satisfies MCP(0, N ). The Juillet's Theorem, which settled the case of the simplest sub-Riemannian structures, paved the way to the study of measure contraction properties for more general sub-Riemannian structures. In [START_REF] Agrachev | Generalized Ricci curvature bounds for three dimensional contact subriemannian manifolds[END_REF], Agrachev and Lee investigated the case of sub-Riemannian structures associated with contact distributions in dimension 3. In [START_REF] Lee | On measure contraction property without Ricci curvature lower bound[END_REF][START_REF] Lee | Ricci curvature type lower bounds for sub-Riemannian structures on Sasakian manifolds[END_REF], Lee and Lee, Li and Zelenko studied the particular case of Sasakian manifolds. In [START_REF] Rifford | Ricci curvatures in Carnot groups[END_REF], the second author proved that any ideal Carnot group satisfy MCP(0, N ) for some N > 1 (it has been shown later by Rizzi [32] that a Carnot group is ideal if and only if it is fat). In [START_REF] Rizzi | Measure contraction properties of Carnot groups[END_REF], Rizzi showed that any co-rank 1 Carnot group of dimension k + 1 (equipped with the sub-Riemannian distance and a left-invariant measure) satisfies MCP(0, k + 3). Finally, more recently, Barilari and Rizzi [START_REF] Barilari | Sharp measure contraction property for generalized Htype Carnot groups[END_REF] proved that H-type Carnot groups of rank k and dimension n satisfy MCP(0, k + 3(n -k)). The purpose of the present paper is to pursue the qualitative approach initiated by the second author in [START_REF] Rifford | Ricci curvatures in Carnot groups[END_REF]. We aim to show that some assumptions on the sub-Riemannian structure insure that the sub-Riemannian distance enjoyes some properties which guarantee that some measure contraction property of the form MCP(0, N ) is satisfied for some N > 1 (in fact N has to be greater or equal to the geodesic dimension of the sub-Riemannian structure as introduced by Rizzi [START_REF] Rizzi | The curvature of optimal control problems with applications to sub-Riemannian geometry[END_REF]). Our approach is purely qualitative, we do not compute any curvature type quantity in order to find the best exponents. Our results are concerned with two-step analytic sub-Riemannian structures and Lipschitz Carnot groups.

Given a (real) analytic manifold M , we say that (∆, g) is analytic if both ∆ and g are analytic on M . Moreover, we recall that a distribution ∆, or a sub-Riemannian structure (∆, g), is two-step if

[∆, ∆](x) := [X, Y ](x) | X, Y smooth sections of ∆ = T x M ∀x ∈ M.
A measure on M is called smooth if it is locally defined by a positive smooth density times the Lebesgue measure L n , our first result is the following:

Theorem 3. Every two-step analytic sub-Riemannian structure on a compact analytic manifold equipped with a smooth measure satisfies MCP(0, N ) for some N > 0.

In the case of Carnot groups which are, as Lie groups equipped with left-invariant sub-Riemannian structures, analytic manifolds with analytic sub-Riemannian structures, the homogeneity allows us to extended the above result to left-invariant Lipschitz distributions.

Following [START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF], we say that a sub-Riemannian structure (∆, g) or a Carnot group whose first layer ∆ is equipped with a left-invariant metric, is Lipschitz if it is complete and the associated geodesic distance

d SR : M × M → R is locally Lipschitz outside of the diagonal D = {(x, y) ∈ M × M | x = y}.
Examples of Lipschitz sub-Riemannian structures include two-step distributions and more generally medium-fat distributions. A distribution ∆ (or a sub-Riemannian structure with distribution ∆ or a Carnot group whose first layer ∆ is equipped with a left-invariant metric) is called medium-fat if, for every x ∈ M and every smooth section X of ∆ with X(x) = 0, there holds

T x M = ∆(x) + [∆, ∆](x) + X, [∆, ∆] (x), (1.1)
where

X, [∆, ∆] (x) := X, [Y, Z] (x) | Y, Z smooth sections of ∆ .
The notion of medium-fat distribution has been introduced by Agrachev and Sarychev in [START_REF] Agrachev | Sub-Riemannian metrics: minimality of singular geodesics versus subanalyticity[END_REF]. Of course, in the case of a Carnot group the property of being medium-fat depends only on the properties of its Lie algebra. Our second result is the following:

Theorem 4. Any Lipschitz Carnot group whose first layer is equipped with a leftinvariant metric and equipped with Haar measure satisfies MCP(0, N ) for some N > 0.

The proofs of Theorem 3 and 4 are based on the fact that squared sub-Riemannian pointed distances d SR (x, •) 2 satisfy a certain property of horizontal semiconcavity. Note that for the moment, we are only able to prove this property in the analytic case under an assumption of compactness of length minimizers. The property of horizontal semiconcavity together with the lipschitzness of d SR (x, •) 2 allows us to give an upper bound for divergence of horizontal gradients of f x which implies the desired measure contraction property.

It is worth to notice that, thanks to a seminal result by Cavaletti and Huesmann [START_REF] Cavalletti | Existence and uniqueness of optimal transport maps[END_REF], measure contraction properties are strongly connected with the well-posedness of the Monge problem for quadratic geodesic distances. We refer the interested reader to [START_REF] Badreddine | Mass transportation on sub-Riemannian structures of rank two in dimension four[END_REF][START_REF] Badreddine | Mass transportation in sub-Riemannian structures admitting singular minimizing geodesics[END_REF] for further details.

We recall that all the notations used throughout the paper are listed in Appendix A. The material required for the proof of the two theorems above is worked out in Section 2. The proofs of Theorems 3 and 4 are respectively given in Sections 3 and 4. indebted to Adam Parusinski for fruitful discussions and the reference to the paper by Denef and Van den Dries [START_REF] Denef | p-adic and real subanalytic sets[END_REF].

Preliminaries

Throughout all this section, (∆, g) denotes a complete sub-Riemannian structure on M of rank m ≤ n.

The minimizing Sard conjecture

The minimizing Sard conjecture is concerned with the size of points that can be reached from a given point by singular minimizing geodesics. Following [START_REF] Rifford | Singulières minimisantes en géométrie sous-Riemannienne[END_REF], given x ∈ M , we set

S x ∆,min g := γ(1) | γ ∈ W 1,2 ∆ ([0, 1], M ), γ sing., d SR (x, γ(1)) 2 = energy g (γ) .
Note that for every x ∈ M , the set S x ∆,min g is closed and contains x (because m < n). Let us introduce the following definition. Definition 5. We say that (∆, g) satisfies the minimizing Sard conjecture at x ∈ M if the set S x ∆,min g has Lebesgue measure zero in M . We say that it satisfies the minimizing Sard conjecture if this property holds for any x ∈ M .

It is not known if all complete sub-Riemannian structures satisfy the minimizing Sard conjecture (see [START_REF] Agrachev | Some open problems. Geometric control theory and sub-Riemannian geometry[END_REF][START_REF] Rifford | Singulières minimisantes en géométrie sous-Riemannienne[END_REF]). The best general result is due to Agrachev who proved in [START_REF] Agrachev | Any sub-Riemannian metric has points of smoothness[END_REF] that all closed sets S x ∆,min g have empty interior. As the next result shows, the minimizing Sard conjecture is related to regularity properties of pointed distance functions. Following Agrachev [START_REF] Agrachev | Any sub-Riemannian metric has points of smoothness[END_REF], we call smooth point of the function y → d SR (x, y) (for a fixed x ∈ M ) any y ∈ M for which there is p ∈ T * x M which is not a critical point of the exponential mapping exp x and such that the projection γ x,p of the normal extremal ψ : [0, 1] → T * M starting at (x, p) is the unique minimizing geodesic from x to y = γ x,p [START_REF] Agrachev | Compactness for Sub-Riemannian length-minimizers and subanalyticity. Control theory and its applications (Grado[END_REF]. By Agrachev's Theorem, the set O x of smooth points is always open and dense in M . The following holds: Proposition 6. Let x ∈ M be fixed, the following properties are equivalent: Furthermore, the function y → d SR (x, y) is smooth on O x and if M and (∆, g) are analytic, then the set O x is geodesically star-shaped at x, that is

γ(s, y) ∈ O x ∀s ∈ (0, 1], ∀y ∈ O x , (2.1)
where γ x (•, y) ∈ W 1,2 ∆ ([0, 1], M ) is the unique minimizing geodesic from x to y. Proof of Proposition 6. Let x ∈ M be fixed. The part (iii) ⇒ (ii) is immediate. Let us prove that (ii) ⇒ (i). By assumption the set of differentiability D of f := d SR (x, •) has full measure in M . Recall that for every y ∈ D, there is a unique minimizing geodesic from x to y which is given by the projection of the normal extremal ψ : [START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF]Lemma 2.15 p. 54]). By Sard's Theorem, the set S of exp x (p) with p ∈ T *

[0, 1] → T * M such that ψ(1) = (y, d SR (x, y)d y f ) (see
x M critical has Lebesgue measure zero in M . Therefore, the set D \ S has full measure and for every y ∈ D \ S there is a unique minimizing geodesic from x to y and it is not singular, which shows that y does not belong to S x ∆,min g . Let us now show that (i) ⇒ (iii). By definition of S x ∆,min g , for every y / ∈ S x ∆,min g all minimizing horizontal paths between x and y are not singular. So repeating the proof of [27, Theorem 3.14 p. 98] (see also [START_REF] Cannarsa | Semiconcavity results for optimal control problems admitting no singular minimizing controls[END_REF]), we can show that the function f : y → d SR (x, y) is locally semiconcave and so locally Lipschitz on the open set U := M \ S x ∆,min g . Thus for every compact set K ⊂ U , there is a compact set P K ⊂ T *

x M such that for every y ∈ K, there is p ∈ P K with exp x (p) = y and H(x, p) = d SR (x, y) 2 /2 (in other words γ x,p : [0, 1] → M is a minimizing geodesic from x to y). By Sard's Theorem, the set S K of exp x (p) with p ∈ P K critical is a closed set of Lebesgue measure zero. For every positive integer k, set (here the diameter of the convex set d + y f is taken with respect to some geodesic distance on

T * M ) Σ k (f ) := y ∈ U | diam(d + y f ) ≥ 1/k .
By local semiconcavity of f in U , each set Σ k (f ) is a closed set in U with Lebesgue measure zero (see [14, Proposition 4.1.3 p. 79]). We claim that

S K := K ∩ k>0 Σ k (f ) ⊂ K ∩ k>0 Σ k (f ) ∪ S K . As a matter of fact, if y ∈ K belongs to ∪ k>0 Σ k (f )\∪ k>0 Σ k (f ), then d + y
f is a singleton and there is a sequence {y l } l converging to y such that all d + y l f have dimension at least one and tend to d + y f . This implies that the covector p such that exp x (p) = y and H(x, p) = d SR (x, y) 2 /2 is critical, which shows that y belongs to S K . By construction, every point in K \ S K is a smooth point. We conclude easily.

It remains to prove the second part. The smoothness of f : y → d SR (x, y) is an easy consequence of the inverse function theorem. As a matter of fact, we can show easily that for every y ∈ O x such that y = exp x (p) with H(x, p) = d SR (x, y) 2 /2 and p ∈ T *

x M non-critical, there is a neighborhood U of y in O x such that

f (z) 2 = 2H(x, exp x (z) -1 ) ∀z ∈ U,
where exp -1

x denotes a local inverse of the exponential mapping from a neighborhood of p to U . To prove (2.1), we argue by contradiction. If there are x ∈ M , y ∈ O x and s ∈ (0, 1) such that z := γ(s, y) ∈ O x then either there are two distinct minimizing geodesics from x to z or there is only one minimizing geodesic from x to y which is singular. In the first case, we infer the existence of two distinct minimizing geodesics from x to y, which contradicts the smoothness of y. In the second case, we deduce that the minimizing geodesic γ x (•, y) is the projection of a normal extremal which is regular and whose restriction to [0, s] is singular. This cannot happen under the assumption of analyticity of the datas because an abnormal extremal above γ x (•, y) over [0, s] could be extended to an abnormal extremal over [0, 1] (see [START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF]Proposition 1.11 p.21]).

Remark 7. By Proposition 6, any (complete) sub-Riemannian structure satisfying the minimizing Sard conjecture has negligeable cut loci.

Remark 8. As pointed out by the referee, the set O x could be replaced by the set A x ⊂ O x of ample points from x. This set is defined as the set of y ∈ O x where the unique minimizing geodesic from x to y is ample, that is whose growth vector saturates the tangent space. It can be shown to be open with full measure and geodesically starshaped in the smooth case, we refer the interested reader to the monograph [START_REF] Agrachev | Curvature: a variational approach[END_REF] for further details. In our case, since we need the analyticity to prove other results, we prefer to work with the simpler O x which is geodesically star-shaped in the analytic case.

Two characterizations for MCP(0, N )

The following result was implicit in the previous paper [START_REF] Rifford | Ricci curvatures in Carnot groups[END_REF] of the second author (it is also the case in [START_REF] Lee | Ricci curvature lower bounds on Sasakian manifolds[END_REF]Page 5] and [START_REF] Ohta | On the curvature and heat flow on Hamiltonian systems[END_REF]Section 6.2]). The measure contraction property MCP(0, N ) is equivalent to some upper bound on the divergence of the horizontal gradient of the squared pointed sub-Riemannian distance. This result holds at least whenever the horizontal gradient is well-defined and the sets O x are geodesically starshaped.

Proposition 9. Assume that (∆, g) satisfies the minimizing Sard conjecture and that all its sets O x are geodesically star-shaped, and let µ be a smooth measure on M and N > 0 be fixed. Then (∆, g) equipped with µ satisfies MCP(0, N ) if and only if

div µ y ∇ h f x ≤ N ∀y ∈ O x , ∀x ∈ M, (2.2)
where

f x : M → R is the function defined by f x (y) := d SR (x, y) 2 /2.
Proof. Let x ∈ M be fixed, the vector field Z := -∇ h f x is well-defined and smooth on O x . Moreover by assumption, every solution of ẏ(t) = Z(y(t)) with y(0) ∈ O x remains in O x for all t ≥ 0, we denote by {ϕ t } t≥0 the flow of Z on O x . For every y ∈ O x , the function θ

: t ∈ [0, +∞) → d SR (ϕ t (y), y) satisfies θ(0) = 0 and θ(t) = length g ϕ [0,t] (y) = t 0 |Z(ϕ s (y)| ds.
So that, for all t ≥ 0,

θ(t) = |Z(ϕ t (y)| = d SR (x, ϕ t (y)) = d SR (x, y) -d SR (y, ϕ t (y)) = d SR (x, y) -θ(t), which yields θ(t) = d SR (x, y) 1 -e -t ∀t ≥ 0.
Consequently, if A ⊂ O x is a Borel set and s ∈ (0, 1], then we have

A s = γ x (s, y) | y ∈ A = ϕ t (A) with t = -ln(s).
Let us now assume that (2.2) is satisfied. By definition of div µ Z, for every x ∈ M and any measurable set A ⊂ O x , we have for every t ≥ 0 (see for example, see [START_REF] Belotto Da Silva | The Sard conjecture on Martinet surfaces[END_REF]

, Proposition B.1]), µ (ϕ t (A)) = A exp t 0 div µ ϕs(y) (Z) ds dµ(y),
which by (2.2) implies with s = e -t ,

µ (A s ) = µ (ϕ t (A)) ≥ A exp (-N t) dµ(y) = s N µ(A).
This shows that (2.2) implies MCP(0, N ). Conversely, if (∆, g) equipped with µ satisfies MCP(0, N ) then for every x ∈ M and every small ball B δ (y) ⊂ O x (say a Riemannian ball with respect to the Riemannian extension g), we have

µ (ϕ t (B δ (y))) = B δ (y) exp t 0 div µ ϕs(y) (Z) ds dµ(y) ≥ e -N t µ (B δ (y)) ∀t ≥ 0.
For every t ≥ 0, letting δ go to 0 yields

exp t 0 div µ ϕs(y) (Z) ds ≥ e -N t .
We infer (2.2) by dividing by t and letting t go to 0.

In the case of Carnot groups, the invariance of the divergence of ∇ h f x by dilation allows us to characterize MCP(0, N ) in term of a control on the divergence over a compact set not containing the origin.

Proposition 10. Let G be a Carnot group whose first layer is equipped with a leftinvariant metric satisfying the minimizing Sard conjecture and N > 0 fixed. Then the metric space (G, d SR ) with Haar measure µ satisfies MCP(0, N ) if and only if

div µ y ∇ h f 0 ≤ N ∀y ∈ O 0 ∩ S SR (0, 1), (2.3)
where

f 0 : M → R is the function defined by f 0 (y) := d SR (0, y) 2 /2.
Proof. Since Carnot groups are indeed analytic, by the second part of Proposition 6 and Proposition 9, it is sufficient to show that (2.3) is equivalent to

div µ y ∇ h f 0 ≤ N ∀y ∈ O 0 . (2.4)
Recall that by taking a set of exponential coordinates (x 1 , . . . , x n ), we can identify G with its Lie algebra g R n and indeed consider that we work with the Lebesgue measure in R n and that the sub-Riemannian structure is globally parametrized by an orthonormal family of analytic vector fields X 1 , . . . , X n in R n satisfying

X i (δ λ (x)) = λ -1 δ λ X i (x) ∀x ∈ R n , ∀i = 1, . . . , n, (2.5)
where {δ λ } λ>0 is a family of dilations defined as (d 1 , . . . , d n are positive integers)

δ λ (x 1 , . . . , x n ) = λ d 1 x 1 , λ d 2 x 2 , . . . , λ dn x n ∀x ∈ R n .
By the homogeneity property, we have d SR (0, δ λ (x)) = λ d SR (0, x) for all x ∈ R n and λ > 0. Then we have

f 0 (δ λ (x)) = λ 2 f 0 (x) and d δ λ (x) f 0 • δ λ = λ 2 d x f 0 ∀x ∈ R n , ∀λ > 0. (2.6)
Recall that the horizontal gradient ∇ h f 0 is given by

∇ h x f 0 = m i=1 X i • f 0 (x)X i (x) ∀x ∈ R n .
Therefore, by (2.5)-(2.6), we infer that for every x ∈ R n and λ > 0,

∇ h δ λ (x) f 0 = m i=1 d δ λ (x) f 0 X i (δ λ (x)) X i (δ λ (x)) = λ -2 m i=1 d δ λ (x) f 0 δ λ X i (x) δ λ X i (x) = m i=1 d x f 0 X i (x) δ λ X i (x) = δ λ ∇ h x f 0 .
We deduce that

div µ δ λ (x) ∇ h f 0 = div µ x ∇ h f 0 ∀x ∈ R n , ∀λ > 0,
which shows that (2.3) and (2.4) are equivalent and concludes the proof.

Nearly horizontally semiconcave functions

Recall that without loss of generality, we can assume that the metric g over ∆ is the restriction of a global Riemannian metric on M . This metric allows us to define the C 2 -norms of functions from R m to M . In the following statement, (e 1 , . . . , e m ) stands for the canonical basis in R m .

Definition 11. Let C > 0 and U an open subset of M , a function f : U → R is said to be C-nearly horizontally semiconcave with respect to (∆, g) if for every y ∈ U , there are an open neighborhood

V y of 0 in R m , a function ϕ y : V y ⊂ R m → U of class C 2 and a function ψ y : V y ⊂ R m → R of class C 2 such that ϕ y (0) = y, ψ y (0) = f (y), f (ϕ y (v)) ≤ ψ y (v) ∀v ∈ V y , (2.7)
d 0 ϕ y (e 1 ), . . . , d 0 ϕ y (e m ) is an orthonormal family of vectors in ∆(y), (

and

ϕ y C 2 , ψ y C 2 ≤ C, (2.9)
where ϕ y C 2 , ψ y C 2 denote the C 2 -norms of ϕ y and ψ y . Remark 12. We refer the reader to [START_REF] Montanari | On the lack of semiconcavity of the subRiemannian distance in a class of Carnot groups[END_REF] for an other notion of horizontal semiconcavity of interest that has been investigated by Montanari and Morbidelli in the framework of Carnot groups.

If m were equal to n that is if we were in the Riemannian case, the above definition would coincide with the classical definitions of semiconcave functions (see [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF][START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF]). Here, in the case m < n, the definition tells that at each point, there is a support function from above of class C 2 which bounds the function along a C 2 submanifold which is tangent to the distribution. This type of mild horizontal semiconcavity will allows us, at least in certain cases, to bound the divergence of the horizontal gradient of squared pointed sub-Riemannian distance functions.

Before stating the main result of this section, we recall that a minimizing geodesic γ : [0, 1] → M from x to y is called normal if it is the projection of a normal extremal, that is a trajectory ψ : [0, 1] → T * M of the Hamiltonian vector field associated with (∆, g). We refer the reader to Appendix A for more details. Here is our result (|p| * = |p| *

x for every (x, p) ∈ T * M is the dual norm associated with g):

Proposition 13. Assume that M and (∆, g) are analytic and let K be a compact subset of M , U ⊂ M a relatively compact open set of M , and A > 0 satisfying the following property: For every x ∈ K, every y ∈ Ū and every minimizing geodesic γ : [0, 1] → M from x to y, there is p ∈ T * x M with |p| * ≤ A such that γ is the projection of the normal extremal ψ : [0, 1] → T * M starting at (x, p). Then there is C > 0 such that for every x ∈ K, the function y → d SR (x, y) 2 is C-nearly horizontally semiconcave in U . There is an open neighborhood U γ of γ([0, 1]) and a family F γ of m analytic vector fields X 1 γ , . . . , X m γ in M such that for every z ∈ U γ the family {X 1 γ (z), . . . , X m γ (z)} is orthonormal with respect to g and parametrize ∆ (that is Span{X 1 γ (z), . . . , X m γ (z)} = ∆(z)) and for every z ∈ M \ U γ , X 1 γ (z), . . . , X m γ (z) belongs to ∆(z). Consider the End-Point mapping from x in time 1 associated with the family

F γ = {X 1 γ , . . . , X m γ }, it is defined by E x,1 Fγ : u ∈ L 2 ([0, 1]; R m ) -→ γ Fγ u (1) ∈ M, where γ Fγ u (1) : [0, 1] → M is the solution to the Cauchy problem γ(t) = m i=1 u i (t)X i γ (γ(t)) for a.e. t ∈ [0, 1], γ(0) = x. (2.10)
Note that taking the vector fields X 1 γ , . . . , X m γ equal to zero outside of an neighborhood of U γ , we may assume without loss of generality that E x,1

Fγ is well-defined on L 2 ([0, 1]; R m ). Recall that the function E x,1

Fγ is smooth and satisfies (see [27, Proposition 1.10 p. 19])

∆ (ȳ) ⊂ Im d uγ E x,1 Fγ , (2.11)
where

u γ is the unique control u ∈ L 2 ([0, 1], R m ) such that γ Fγ u = γ. Therefore, there are v 1 γ , . . . , v m γ ∈ L 2 ([0, 1], R m ) such that d uγ E x,1 Fγ v i γ = X i γ (ȳ) ∀i = 1, . . . , m. (2.12) Define ϕ γ : R m → M by ϕ γ (v) := E x,1 Fγ u γ + m i=1 v i v i γ ∀v = (v 1 , . . . , v m ) ∈ R m .
By construction, ϕ γ is smooth and satisfies

ϕ γ (0) = E x,1 Fγ (u γ ) = γ(1) = ȳ,
and

d 0 ϕ γ (e i ) = d uγ E x,1 Fγ v i γ = X i γ (ȳ) ∀i = 1, . . . , m.
Moreover, for every v ∈ R m such that the solution to (2.10) associated with the control

u γ + m i=1 v i v i γ remains in U γ , we have d SR (x, ϕ γ (v)) 2 ≤ u γ + m i=1 v i v i γ 2 L 2 =: ψ γ (v).
By construction, ϕ γ and ψ γ are smooth, defined in a neighborhood of 0 ∈ R m and satisfy (2.7)-(2.8). It remains to show that the C 2 norms of ϕ γ and ψ γ can be taken to be uniformly bounded, this is the purpose of the next lemma.

Lemma 14. There are neighborhoods U x and U ȳ respectively of x and ȳ in M , a neighborhood U γ of u γ in L 2 ([0, 1], R m ) and C γ > 0 such that for every x ∈ K ∩ U x, y ∈ Ū ∩U ȳ and every control u x,y associated with a minimizing geodesic γ x,y : [0, 1] → M from x to y with u x,y ∈ U γ , there are

v 1 ux,y , . . . , v m ux,y ∈ L 2 ([0, 1], R m ) such that d ux,y E x,1 Fγ v i ux,y = X i γ (y) ∀i = 1, . . . , m (2.13)
and

v i ux,y L 2 ≤ C γ ∀i = 1, . . . , m. (2.

14)

Proof of Lemma 14. Note that if we prove for each i ∈ {1, . . . , m} the existence of neighborhoods U i x , U i y and U i γ such that (2.13)-(2.14) are satisfied for i, then the result follows by taking the intersections of the neighborhoods U i

x, U i ȳ and U i γ . So, let us fix i in {1, . . . , m}. By taking a chart on a neighborhood of γ([0, 1]) (that we still denote by U γ ) we may assume that the restriction of E x,1

Fγ to a neighborhood of u γ is valued in R n and doing a change of coordinates (y 1 , . . . , y n ) in a neighborhood U 0 ȳ we may also assume that there are analytic vector fields Y 1 , . . . , Y m such that

Y 1 (y) = X i (y) = ∂ y 1 , Y j (y) = ∂ y j + n l=m+1 a j l (y) ∂ y l ∀j = 2, . . . , m, ∀y ∈ U 0 ȳ and ∆(z) = Span Y 1 (z), . . . , Y m (z) ∀z ∈ U γ .
Observe that by construction, there are analytic mappings f j k on U γ for j ∈ {1, . . . , m} and k ∈ {1, . . . , m} such that

Y j (z) = m k=1 f j k (z) X k (z) ∀z ∈ U γ , ∀j = 1, . . . , m.
As a consequence, if γ : [0, 1] → U γ is solution to

γ(t) = m j=1 w j (t) Y j (γ(t)) for a.e. t ∈ [0, 1], (2.15)
for some w ∈ L 2 ([0, 1], R m ), then there holds for a.e. t ∈ [0, 1],

γ(t) = m j=1 w j (t) m k=1 f j k (γ(t)) X j (γ(t)) = m k=1   m j=1 w j (t) f j k (γ(t))   X k (γ(t)).
This shows that if we denote by F the family {Y 1 , . . . , Y m }, then any horizontal curve parametrized by F can be parametrized by F γ as above. On the other hand, any parametrization with respect to F γ leads to a parametrization with respect to F. Consequently, there is a control w ∈ L 2 ([0, 1), R m ) such that the solution of (2.15) starting at x is equal to γ and there are neighborhoods U 0 x of x and W of w together with a mapping G :

U 0 x × W → L 2 ([0, 1], R m ) of class C 1 such that E x,1 F (w) = E x,1 Fγ (G(x, w)) ∀x ∈ U 0 x , ∀w ∈ W,
where E x,1 F denotes the End-Point mapping from x in time 1 associated with the family F. Therefore, in order to prove our result it is sufficient to prove that there is C > 0 such that for every x in K close to x, every y in Ū close to ȳ and any γ minimizing geodesic from x to y close to γ associated to the control w ∈ W with respect to F, there is ω ∈ L 2 ([0, 1], R m ) such that

d w E x,1 F (ω) = X i (y) = ∂ y 1 and ω L 2 ≤ C. (2.16)
Let x ∈ U 0 x and w ∈ W be fixed, the differential of E x,1 F at w is given by (see [27, Remark 1.5 p.15])

d w E x,1 F (v) = S(1) 1 0 S(t) -1 B(t)v(t) dt ∀v ∈ L 2 ([0, 1], R m ),
where S : [0, 1] → M n (R) is solution to the Cauchy problem Ṡ(t) = A(t)S(t) for a.e. t ∈ [0, 1], S(0) = I n and where the matrices

A(t) ∈ M n (R), B(t) ∈ M n,m (R) are defined by (note that J Y 1 = J X i = 0) A(t) := m j=2 w j (t) J Y j (γ(t)) a.e. t ∈ [0, 1], B(t) = Y 1 (γ(t)) , . . . , Y m (γ(t)) ∀t ∈ [0, 1].
By construction of Y 1 , . . . , Y m , there is τ 0 ∈ (0, 1) (which does not depend upon x, y but only upon x, ȳ) such that for almost every t ∈ [1 -τ 0 , 1] the matrix A(t) has the form (2.17)

A(t) =           0 0 • • • • • • 0 0 • • • 0 0 . . . . . . 0 0 m-1,m-1 0 m-1,n-m α(t) β(t) δ(t)           , where α(t) is in R n-m , β(t) = (β(t) k,l ) is a (n -m) × (m -
and the column vector s(t) with coordinates (s m+1 (t), . . . , sn (t)) is solution of

ṡ(t) = D(t) α(1 -t) ∀t ∈ [0, τ 0 ], s(0) = 0 (2.18)
with D(t) the (n -m) square matrices satisfying

Ḋ(t) = D(t)δ(1 -t) ∀t ∈ [0, τ 0 ], D(0) = I n-m . (2.19)
Thus, a way to solve Remember now that by assumption, M and (∆, g) are analytic and for every x ∈ K and y ∈ Ū all minimizing geodesics joining x to y are normal with initial covector bounded by A. Let p ∈ T *

d w E x,1 F (ω) = ∂ y 1 is to take ω ∈ L 2 ([0, 1], R m ) of the form ω(t) = (ω 1 (1 -•), 0, . . . , 0) a.e. t ∈ [0, 1] (2.20) with Supp(ω 1 ) ⊂ [0, τ 0 ],
x M be such that γ = γ x,p where for any (x, p) ∈ T * M , γ x,p : [0, 1] → M denotes the projection of the normal extremal ψ x,p : [0, 1] → T * M starting at (x, p) (see Appendix A). Then, the desired result will follow if we show that there are a neighborhood T of (x, p) in T * M and C > 0 such that for every (x, p) ∈ T , the point x belongs to U 0

x , the curve γ x,p ([0, 1]) is contained in U γ , the associated control w = w x,p belongs to W and there is ω = ω x,p ∈ L 2 ([0, 1], R m ) satisfying (2.20)-(2.21) (with the function s = sx,p = (s m+1 (t), . . . , sn (t)) : [0, 1] → R n-m given by the first column of S(t) = Sx,p (t) := S(1)S(1 -t) -1 associated with γ x,p and w x,p which satisfies (2.17)-(2.18) with α = α x,p and δ = δ x,p associated with γ x,p and w x,p as well) such that

ω x,p L 2 ≤ C. (2.22)
We need the following lemma whose proof is postponed to Appendix B.

Lemma 15. Let K ⊂ R l be a compact set and h : [0, 1] × K → R be an analytic mapping such that h(0, κ) = 0 for all κ ∈ K. Then there are τ > 0 as small as desired and ν ∈ (0, 1) such that

τ 0 h(t, p) dt 2 ≤ ν τ τ 0 h(t, p) 2 dt ∀κ ∈ K. (2.23)
Let T be a compact neighborhood of of (x, p) in T * M such that for every (x, p) ∈ T , the point x belongs to U 0

x , the curve γ x,p ([0, 1]) is contained in U γ and the associated control w = w x,p belongs to W. Then we note that the mapping

(t, x, p, λ) ∈ [0, 1] × T × [-1, 1] n-m -→ h(t, x, p, λ) := n-m k=1 λ k sx,p k (t)
is analytic. Therefore, by Lemma 15, there are τ ∈ (0, τ 0 ) and ν ∈ (0, 1) such that

τ 0 h(t, x, p, λ) dt 2 ≤ ντ τ 0 h(t, x, p, λ) 2 dt ∀(x, p) ∈ T , ∀λ ∈ [-1, 1] n-m .(2.24) Define I ∈ L 2 ([0, 1], R) by I(t) := 1 [0,τ ] (t) for all t ∈ [0, 1] and for every (x, p) ∈ T denote by P ∈ L 2 ([0, 1], R) the orthogonal projection of I in L 2 ([0, 1], R) over the vector space V x,p := f ∈ L 2 ([0, 1], R) | f, 1 [0,τ ] sx,p l L 2 = 0, ∀l = m + 1, . . . , n .
Let (x, p) ∈ T be fixed, then there are Λ x,p 1 , . . . , Λ x,p n-m ∈ R such that

P = I + n-m k=1 Λ k 1 [0,τ ]
sx,p m+k =: I + J and there holds

P, 1 [0,τ ] sx,p l L 2 = 1 0 P (t)1 [0,τ ] (t) sx,p l (t) dt = 0 = τ 0 1 [0,τ ] (t)P (t) sx,p l (t) dt ∀l = m + 1, . . . , n -m and 1 0 1 [0,τ ] (t)P (t)(t) dt = P, I L 2 = P 2 L 2 = I 2 L 2 + J 2 L 2 + 2 I, J L 2 ≥ I 2 L 2 + J 2 L 2 -2 √ ν I L 2 J L 2 ≥ 1 - √ ν I 2 L 2 = ν 1 - √ ν ,
where we used that setting Λk = max{Λ 1 , . . . , Λ n-m } (note that Λk = 0) thanks to (2.24) yields

( I, J L 2 ) 2 = τ 0 n-m k=1 Λ k sx,p k (t) dt 2 = Λ 2 k τ 0 n-m k=1 Λ k Λk sx,p k (t) dt 2 ≤ Λ 2 k ντ τ 0 n-m k=1 Λ k Λk sx,p k (t) 2 dt = ν I 2 L 2 J 2 L 2 .
In conclusion, we deduce that for every (x, p) ∈ T , the function ω = ω x,p ∈ L 2 ([0, 1], R m ) of the form (2.20) with ω 1 = ω x,p 1 given by

ω x,p 1 (t) := 1 [0,τ ] (t)P (t) P 2 L 2 ∀t ∈ [0, 1]
satisfies (2.21) and

ω x,p 1
L 2 ≤ 1 ν (1 - √ ν) .
The proof of Lemma 14 is complete.

We conclude easily by compactness of K and Ū .

Proof of Theorem 3

Let M be an analytic compact manifold, (∆, g) a two-step analytic sub-Riemannian structure and µ a smooth measure on M . The following result, due to Agrachev and Lee [START_REF] Agrachev | Optimal transportation under nonholonomic constraints[END_REF] (see also [START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF]), is a consequence of the fact that ∆ is two-step (and the compactness of M ). We refer the reader to [START_REF] Agrachev | Optimal transportation under nonholonomic constraints[END_REF][START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF] for the proof. In fact, it is worth mentioning that Agrachev and Lee prove that a sub-Riemannian structure is two-step if and only if d SR is locally Lipschitz in charts.

Lemma 16. The function d 2 SR : M × M → R is locally Lipschitz in charts. In particular, there is L > 0 such that |∇ y f x | ≤ L for all x ∈ M and y ∈ O x . Furthermore, there is A > 0 such that for every x, y ∈ M and every minimizing geodesic γ : [0, 1] → M from x to y, there is p ∈ T *

x M with |p| * ≤ A such that γ is the projection of the normal extremal ψ : [0, 1] → T * M starting at (x, p).

By the above lemma and Proposition 13, there is C > 0 such that for every x ∈ M the function

f x : y → d SR (x, y) 2 /2 is C-nearly horizontally semiconcave in M .
Lemma 17. There is B > 0 such that for every x ∈ M the following property holds: for every y ∈ O x , there is a neighborhood U y ⊂ O x of y along with an orthonormal family of smooth vector fields X 1 , . . . , X m which parametrize ∆ in U y such that

X i C 1 ≤ B ∀i = 1, . . . , m, (3.1)
and

X i • (X i • f x ) (z) ≤ B |∇ z f x | + B ∀z ∈ U y , ∀i = 1, . . . , m, (3.2)
where ∇ z f x stands for the gradient of f x at z with respect to the global Riemannian metric g.

Proof of Lemma 17. First of all, we notice that there is A > 0 such that if v 1 , . . . , v m is an orthonormal family of tangent vectors in ∆(z) for some z ∈ M then there is an orthonormal family of smooth vector fields X 1 , . . . , X m which generates the distribution ∆ in a neighborhood of z and such that X i C 1 is bounded by A for all i = 1, . . . , m. Let x ∈ M be fixed, then by C-nearly horizontal semiconcavity of f x , for every y ∈ M , there are an open neighborhood

V y of 0 in R m , a function ϕ y : V y ⊂ R m → U of class C 2 and a function ψ y : V y ⊂ R m → R of class C 2 such that (2.7) (with f = f x ), ( 2
.8) and (2.9) are satisfied. Fix y ∈ O x and define the function F y : U y → R by F y := f x • ϕ y -ψ y , it is of class C 2 and satisfies d 0 F y = 0 and Hess 0 F y ≤ 0.

Taking a chart near y we can assume that we work in R n . Let ϕ y = (ϕ y 1 , . . . , ϕ y n ) and (x 1 , . . . , x n ) and (v 1 , . . . , v m ) the coordinates respectively in R n and R m . Then, we have

∂F y ∂v i (0) = n k=1 ∂f x ∂x k (y) ∂ϕ y k ∂v i (0) - ∂ψ y ∂v i (0) = 0 ∀i = 1, . . . , m
and for every i = 1, . . . , m,

∂ 2 F y ∂v 2 i (0) =   n k,l=1 ∂ 2 f x ∂x l ∂x k (y) ∂ϕ y k ∂v i (0) ∂ϕ y l ∂v i (0)   + n k=1 ∂f x ∂x k (y) ∂ 2 ϕ y k ∂v 2 i (0) - ∂ 2 ψ y ∂v 2 i (0) ≤ 0, which yields n k,l=1 ∂ 2 f x ∂x l ∂x k (y) ∂ϕ y k ∂v i (0) ∂ϕ y l ∂v i (0) ≤ ∂ 2 ψ y ∂v 2 i (0) - n k=1 ∂f x ∂x k (y) ∂ 2 ϕ y k ∂v 2 i (0) ≤ C + C |∇ y f x | . (3.3)
By (2.8) and the observation made at the very beginning of this proof, there is an orthonormal family of smooth vector fields X 1 , . . . , X m which generates the distribution ∆ in a neighborhood of z and such that

X i C 1 ≤ A and d 0 ϕ y (e i ) = ∂ϕ y ∂v i (0) = X i (y) ∀i = 1, . . . , m. (3.4)
Setting X i = n k=1 a i k ∂ k , we check easily that

X i • f x = n k=1 a i k ∂f x ∂x k and X i • X i • f x = n k=1 n l=1 a i l ∂a i k ∂x l ∂f x ∂x k + n k=1 a i k n l=1 a i l ∂ 2 f x ∂x l ∂x k .
The last expression at y yields, thanks to (3.3) and (3.4) (which implies a i k (y) = ∂ϕ y k ∂v i (0) for all i = 1, . . . , m and k = 1, . . . , m)

X i • (X i • f x ) (y) ≤ A 2 |∇ y f x | + C + C |∇ y f x | ∀i = 1, . . . , m.
We conclude easily by smoothness of f x in O x with U y sufficiently small and B > 0 sufficiently large.

Lemma 18. There is N > 0 such that the following property holds:

div µ y ∇ h f x ≤ N ∀y ∈ O x , ∀x ∈ M. (3.5)
Proof of Lemma 18. Let x ∈ M and y ∈ O x be fixed, by Lemma 17 there is a neighborhood U y ⊂ O x of y along with an orthonormal family of smooth vector fields X 1 , . . . , X m which parametrize ∆ in U y such that (3.2) holds. The horizontal gradient of f x in U y is given by

∇ h y f x = m i=1 X i • f x (y)X i (y).
So, we have

div µ y ∇ h f x = m i=1 X i • f x (y) div µ y X i + m i=1 X i • (X i • f ) (y).
The second term above (in the right-hand side) is bounded thanks to (3.2) and Lemma 16 and the first term is bounded by (3.1) and Lemma 16 (the quantities X i • f x (y) are indeed bounded by the fact that d SR (x, •) is solution to the horizontal eikonal equation, see [START_REF] Figalli | Mass Transportation on sub-Riemannian Manifolds[END_REF]). The proof of Lemma 18 is complete.

To conclude the proof of Theorem 3, we observe that as the functions f x = d SR (x, •) 2 /2 are Lipschitz on M the minimizing Sard conjecture is satisfied (by Proposition 6) and we note that by analyticity the the sets O x are geodesically star-shaped. Then we can apply Proposition 9 together with Lemma 18.

Proof of Theorem 4

By Proposition 10, it is sufficient to show that (2.3) holds. By assumption the function f 0 : y → d SR (0, y) 2 /2 is locally Lipschitz on G \ {0}. Thus for every relatively compact open neighrborhood U of S SR (0, 1) with Ū ⊂ G \ {0}, there is A > 0 such that for every y ∈ Ū and every minimizing geodesic γ : [0, 1] → M from 0 to y, there is p ∈ T * 0 M with |p| * ≤ A such that γ is the projection of the normal extremal ψ : [0, 1] → T * M starting at (0, p). Moreover, a Carnot group whose first layer is equipped with a leftinvariant metric is an analytic manifold equipped with an analytic sub-Riemannian structure. Consequently, by Proposition 13, the function f 0 : y → d SR (0, y) 2 /2 is C-nearly horizontally semiconcave in Ū and we can repeat the arguments used in the proof of Theorem 3 for y ∈ O 0 ∩ S SR (0, 1).

A Notations

We list below the notations used throughout this paper, we refer the reader to the monographs [START_REF] Agrachev | Introduction to Riemannian and sub-Riemannian geometry[END_REF][START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF][START_REF] Montgomery | A tour of sub-Riemannian geometries, their geodesics and applications[END_REF][START_REF] Rifford | Sub-Riemannian Geometry and Optimal Transport[END_REF] for further details:

• M is a smooth manifold of dimension n ≥ 3.

• ∆ is a smooth totally nonholonomic distribution of rank m < n.

• g is a smooth metric over ∆. Sometimes, we see g as the restriction of a global Riemannian metric g on M . We use the notation •, • instead of g x (•, •) and we denote the norm associated with g by

| • | (instead of | • | x = g x (•, •) 1/2 ). B r (x)
stands for the open geodesic ball of radius r > 0 centered at x. The dual norm associated with the Riemannian metric g on each T * x M is denoted by |p| * = |p| * x for every (x, p) ∈ T * M .

• We call horizontal path any γ : [0, 1] → M in W 1,2 which is almost everywhere tangent to ∆. We denote by W 1,2 ∆ ([0, 1], M ) the set of horizontal paths γ : [0, 1] → M endowed with the W 1,2 -topology.

• For every γ ∈ W 1,2 ∆ ([0, 1], M ), we define the length of γ (w.r.t. g) by length g (γ) = • For any x, y ∈ M , we denote by d SR (x, y) (resp. e SR (x, y)) the infimum of lengths (resp. energies) of horizontal paths joining x to y. We note that e SR = d 2 SR . We denote the open ball and the sphere centered at x with radius r > 0 respectively by B SR (x, r) and S SR (x, r). The geodesic distance d SR is said to be complete if the metric space (M, d SR ) is complete. In this case, all closed balls BSR (x, r) are compact (for any x ∈ M and any r > 0).

• We call minimizing geodesic from x to y any γ ∈ W 1,2 ∆ ([0, 1], M ) with γ(0) = x, γ(1) = y which minimizes the energy e SR (x, y) (and so the distance d SR (x, y)), that is such that energy g (γ) = e SR (γ). We note that if d SR is complete, then there exist minimizing geodesics between any pair of points.

• For every x ∈ M , we denote by W 1,2 ∆,x ([0, 1], M ) the set of paths in W 1,2 ∆ ([0, 1], M ) starting at x (that is γ(0) = x) and we define the end-point map

E x ∆ : W 1,2 ∆,x ([0, 1], M ) -→ M by E x ∆ (γ) = γ(1)
. The infinite dimensional space W 1,2 ∆,x ([0, 1], M ) has a smooth manifold structure and the end-point map E x ∆ is smooth.

• An horizontal path γ ∈ W 1,2 ∆,x ([0, 1], M ) is called singular if it is singular with respect to the end-point map E x ∆ , that is if the differential d γ E x,1
∆ is not surjective. It is convenient to rewrite the definition of singular curves in term of singular controls. If the distribution ∆ is parametrized by a family F of k smooth vector fields X

1 , . . . , X k in a open neighborhood of γ([0, 1]) and if u ∈ L 2 ([0, 1], R k ) satisfies γ(t) = k i=1 u i (t) X i (γ(t)) for a.e. t ∈ [0, 1],
then γ is singular if and only if the control u is a singular point of the smooth mapping (well-defined in an open set U)

E x,1 F : U ⊂ L 2 ([0, 1], R k ) -→ M defined by E x,1 F (v) := γ v (1) ∀v ∈ L 2 ([0, 1], R k ),
where

γ v is the curve in W 1,2 ∆,x ([0, 1], M ) solution to the Cauchy problem γv (t) = k i=1 v i (t) X i (γ v (t)) for a.e. t ∈ [0, 1], γ v (0) = x. • An horizontal path γ ∈ W 1,2 ∆,x ([0, 1], M
) is singular if and only if it is the projection of an abnormal extremal ψ : [0, 1] → T * M that never intersects the zero section of T * M , such that

ψ(t) = k i=1 u i (t) h i (ψ(t)) for a.e. t ∈ [0, 1],
where F is a family of k smooth vector fields X 1 , . . . , X k which parametrizes ∆ in a open neighborhood of γ([0, 1]) and h 1 , . . . , h k are the Hamiltonian vector fields associated canonically with h i (x, p) = p • Xi(x) in T * M . The curve ψ is called an abnormal lift of γ and γ is said to be abnormal.

• The Hamiltonian H : T * M → R associated with (∆, g) is defined by

H(x, p) := 1 2 max p(v) 2 g x (v, v) | v ∈ ∆(x) \ {0} ∀(x, p) ∈ T * M, which coincides with 1 2 m i=1 p • X i (x) 2 ,
if ∆ is parametrized locally by an orthonormal family X 1 , . . . , X m . The Hamiltonian vector field H associated with (∆, g) is the Hamiltonian vector field given by H with respect to the canonical symplectic form on T * M . In local coordinates (x, p) the trajectories ψ = (x, p) of H are solution to

ẋ = ∂H ∂p (x, p), ṗ = - ∂H ∂x (x, p),
we call them normal extremals. Any projection of a normal extremal is an horizontal path that is said to be normal.

• An horizontal path γ is called strictly abnormal if it is abnormal (singular) and not normal.

• For every x ∈ M , the exponential mapping exp x : T * x M → M associated with (∆, g) at x is defined by exp x := π(ψ x,p (1)) where ψ x,p is the trajectory of H starting at (x, p) and π : T *

x M → M is the canonical projection.

• A Carnot group (G, ) of step s is a simply connected Lie group whose Lie algebra g = T 0 G (we denote by 0 the identity element of G) admits a nilpotent stratification of step s, i.e.

g = V 1 ⊕ • • • ⊕ V s , (A.1) with V 1 , V j = V j+1 ∀1 ≤ j ≤ s, V s = {0}, V s+1 = {0}. (A.2)
By simple-connectedness of G and nilpotency of g, exp G is a smooth diffeomorphism, which allows to identify G with its Lie algebra g R n . If the first layer V 1 of G is equipped with a left-invariant metric, then there is a set of coordinates (x 1 , . . . , x n ), a one-parameter family of dilations {δ λ } λ>0 of the form δ λ (x 1 , . . . , x n ) = λ d 1 x 1 , λ d 2 x 2 , . . . , λ dn x n ∀x ∈ R n , and a orthonormal family of left-invariant vector fields generating V 1 satisfying

X i (δ λ (x)) = λ -1 δ λ X i (x) ∀λ > 0, x ∈ R n .
• A function f : U → R on a open set U ⊂ M is called locally semiconcave if for every x ∈ U there are a open neighborhood V ⊂ U of x and C > 0 such that for any y ∈ V there is a function ψ : M → R with ψ C 2 ≤ C such that f ≤ V on M and f (y) = ψ(y). For every y ∈ U , d + y f denotes the set of super-differentials of f at y, it is the set of α ∈ T *

x M for which there is a function of class C 1 , ψ : M → R such that ψ ≥ f on M , ψ(y) = f (y) and d y ψ = α.

• If f : U → M is smooth on the open set U ⊂ M , ∇ h f denotes its horizontal gradient with respect to (∆, g). For every y ∈ U , ∇ h y f is defined as the unique v ∈ ∆(y) such that d y f (w) = v, w for all w ∈ ∆(y). If ∆(y) is generated by an orthonormal family X 1 (y), . . . , X m (y), then ∇ h y f = m i=1 (X i • f )(y) X i (y).

B Proof of Lemma 15

Let K ⊂ R l be a compact set and h : [0, 1] × K → R an analytic mapping such that h(0, κ) = 0 for all κ ∈ K.

  (i) the structure (∆, g) satisfies the minimizing Sard conjecture at x ∈ M , (ii) the function y → d SR (x, y) is differentiable almost everywhere in M , (iii) the set of smooth points O x is an open set with full measure in M .

Proof of Proposition 13 .

 13 Let K be a compact set of M , U be a relatively compact open set of M and x ∈ K fixed, let us first show how to construct functions ϕ ȳ, ψ ȳ of class C 2 satisfying (2.7)-(2.8) for ȳ ∈ Ū . Pick a minimizing geodesic γ : [0, 1] → M from x to ȳ = γ(1).

  1) matrix and δ(t) is a (n -m) square matrix. Therefore, as the solution of the Cauchy problem Ṡ(t) = S(t)A(1 -t), S(0) = I n , the first column of the matrix S(t) := S(1)S(1 -t) -1 with t ∈ [0, 1] has the form ) = 0 ∀j = 2, . . . , m ∀t ∈ [0, τ 0 ]

1 0ω 1 (t) dt = 1 and 1 0ω 1

 1111 (t) sl (t) dt = 0 ∀l = m + 1, . . . , n. (2.21)

1 0 1 0

 11 | γ(t)| dt and its energy (w.r.t. g) by energy g (γ) = | γ(t)| 2 dt.

1 . 2 dt = 1 .cc

 121 Let κ ∈ K be fixed, then by [16, Lemma 4.12 p.126], there are an integer d > 0 and ρ > 0 together with analytic functions a 1 , . . . , a d on B(κ, ρ) and b 1 , . . . , b d on [0, ρ] × B(κ, ρ) such thath(t, κ) = d k=1 a k (κ) b k (t, κ) t k ∀(t, κ) ∈ [0, ρ] × B(κ, ρ) and b k (0, κ) = 1 ∀k = 1, . . . , d.Hence, by compactness of K we infer that there are an integer d > 0 and ρ > 0 such that for every κ ∈ K there are analytic functions a κ 1 , . . . , a κ d on B(κ, ρ) and bκ 1 , . . . , b κ d on [0, ρ] × B(κ, ρ) such that h(t, κ) = d k=1 a κ k (κ) b κ k (t, κ) t k ∀(t, κ) ∈ [0, ρ] × B(κ, ρ) and b k (t, κ) ∈ [1/2, 1] ∀(t, κ) ∈ [0, ρ] × B(κ, ρ), ∀k = 1, . . . , d.Let τ ∈ (0, ρ) be fixed. By Cauchy-Schwarz inequality, for every κ ∈ K we have , κ) 2 dt with equality only ifλ 1 h 2 (•, κ) = λ 2 h(•, κ) for some nonzero (λ 1 , λ 2 ) ∈ R 2 .Since h(0, κ) = 0 the case of equality may only happen whenever h(•, κ) ≡ 0 where there holdsτ 0 h(t, κ) dt 2 ≤ ντ τ 0 h(t, κ) 2dt for all ν ∈ (0, 1). In conclusion, by compactness ofK if there is no ν ∈ (0, 1) such that , κ) 2 dt ∀κ ∈ K,then there is a sequence {κ l } l in K converging to some κ ∈ K such that h (•, κ l ) ≡ 0 ∀l and lim l→∞ τ 0 h (t, κ l ) dt 2 τ τ 0 h (t, κ l ) 2 dt = Then there are analytic functions a κ 1 , . . . , a κ d , b κ 1 , . . . , b κ d and a sequence { kl } l in {1, . . . , d} such that for all l, a κ kl (κ l ) = 0 and |a κ k(κ l )| ≤ a κ kl (κ l ) ∀k ∈ 1, . . . , d , l ) a κ kl (κ l ) b κ k (t, κ l ) t k dt l ) a κ kl (κ l ) b κ k (t, κ l ) t kTherefore, since all the quantities a κ k (κ l ) /a κ kl (κ l ) belong to [-1, 1] for all k and are equal to 1 for k = kl , and since the functions b 1 (•, κ), . . . , b d(•, κ) are analytic, there are c 1 , . . . c d in [-1, 1] and k ∈ {1, . . . , d} with ck = 1 such that This means that there is a nonzero pair(λ 1 , λ 2 ) ∈ R 2 such that k b κ k (t, κ) t k = 0 ∀t ∈ [0, τ ].Let k 0 ∈ {1, . . . , d} be such that c k = 0 for all k ∈ {1, . . . , d -1} with k < k 0 . Then we havec k 0 b k 0 (t, κ) t k 0 + t k 0 k b k (t, κ) t k-k 0   = 0 ∀t ∈ [0, τ ],which is impossible because b 0 (0, κ) ∈ [1/2, 1]. The proof is complete.

Acknowledgement. The authors are grateful to the referee for useful remarks and for pointing out a gap in the initial proof of Proposition 12. The authors are also


X1

-4 -2 0 2 4 6 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 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 qq q q q q q q q q q q q q -6 -4 -2 0 2 4 6 -6

-4 -2 0 2 4 6 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 -4 -2 0 2 4 6 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 X2 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 qq 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 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 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 X3 -6 -4 -2 0 2 4 6 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 -6 -4 -2 0 2 4 6 -6

-4 -2 0 2 4 6 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 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 -6 -4 -2 0 2 4 6 q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q X4 9/24 Vincent VANDEWALLE 1,2 Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Natural extension to the heterogeneous data setting Identifiability Model identifiability is directly obtained from the identifiability of Gaussian mixture with local independence (Teicher, 1963(Teicher, , 1967)).

EM algorithm

Starting from the initial value θ [0] , iteration [r] is composed of two steps: E-step Computation of the fuzzy partitions t

[r] ibg := E[Z ibg |x i{b} , m, θ [r-1] ], hence for b = 1, . . . , B, for g = 1, . . . , G b , for i = 1, . . . , n t [r] ibg = π [r-1] bg j∈Ω b p(xij|α [r-1] jg ) G b k=1 π [r-1] bk j∈Ω b p(xij|α [r-1] jg ) ,
M-step Maximization of the expected value of the complete-data log-likelihood over the parameters,

π [r] bg = n [r] bg n and α [r] jg = arg max α jg ∈Θ j Q(αjg|xj, t [r] ω j g ),
where Q(αjg|xj, t b ) = n i=1 t ibg ln p(xij|αjg).

Remarks

Independence between the B blocks of variables permits to maximize the observed-data log-likelihood on each block separately.

Possible modification to perform the block estimation and the parameter inference simultaneously.

In practice the number of blocks, the repartition of variables into blocks, and the number of classes per block are unknown. The EM algorithm to achieve arg max (ω,θ) pen(θm|m, x), starting from (ω [0] , θ [0] ) is at iteration [r]: E-step Computation of the fuzzy partitions t

[r] ibg := E[Z ibg |xi, m, θ [r-1] ], hence for b = 1, . . . , B, for g = 1, . . . , G b , for i = 1, . . . , n t [r] ibg = π [r-1] bg j∈Ω [r-1] b p(xij|α [r-1] jg ) G b k=1 π [r-1] bk j∈Ω [r-1] b p(xij|α [r-1] jg ) , M-step1
Updating the affectation of the variables to blocks

ω [r] j = arg max ω j ∈{1,...,B}   Gω j g=1 max α jg ∈Θ j Q(αjg|xj, t [r] ω j g ) - νjGω j 2 ln n   , M-step2 Updating the model parameters π [r] bg = n [r] bg n and α [r] jg = arg max α jg ∈Θ j Q(αjg|xj, t [r] ω [r] j g
).

[r-1] {b} , z [r] b ) ≥ p(x [r-1] {b} , z [r-1] b ), where x [r-1] {b} = (xj; ω [r-1] = b).
Model step: find ω [r] such that for j = 1, . . . , d ω Extension to heterogeneous co-clustering by adding one level of latent variable to summarize the multi-partition by a single partition, while keeping the partition of the variables.

Derive some k-means type multi-partition similary to Witten & Tibshirani (2010) in the variable selection framework to deal with the very high dimensional setting.

  be decomposed in B independent blocks. The block b follows a mixture with G b components (for b = 1, . . . , B), with the assumption of class conditional independence of the variables. Notations ω = (ωj; j = 1, . . . , d) the repartition of the variables in blocks; ωj = b if variable j belongs to block b. m = (G1, . . . , GB, ω) defines the model Ω b = {j : ωj = b} the subset of variables belonging to block b θ = (π, α) model parameters p(.|αjg) the pdf of the distribution of parameters αjg. Probability distribution function of xi p(xi|m, θ) = B b=1 p(x i{b} |m, θ) with p(x i{b} |m, θ) = Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Variable selection in clustering Multiple Gaussian Mixture Proposed Multiple Partitions Mixture Properties of the model Illustration : n = 100 from a MPM with B = 2 blocks of two variables. Variable 1 and 2 belong to block 1 and variables 3 and 4 in block 2 Each block follows a bi-component Gaussian mixture (i.e., G b = 2) with equal proportions (i.e., π bg = 1/2) and µj1 = 4, µj2 = -4 and σjg = 1.

  by subsets of variables. Generalizes approaches used for variable selection in model-based clustering (if B = 2 and G1 = 1 then variables belonging to block 1 are not relevant for the clustering, while variables belonging to block 2 are relevant) MGMM permits variable selection and multiple partitions explained by subsets of variables (variables classification). Sparse model: number of parameters νm = B b=1 (G b -1) + G b j∈Ω b νj Better model search expected than in the model of Galimberti & al. (2017)

  Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Observed-data likelihood for sample x and model m (θ|m, x) = z = (z ib ; i = 1, . . . , n; b = 1, . . . , B) vectors of the component memberships z ib = (z ib1 , . . . , z ibG b ) where z ibg = 1 if observation i arose from component g for block b, and z ibg = 0 otherwise Completed-data likelihood for sample x and model m (θ|m, x, z) = B b=1 ln p(z b |π b ) + d j=1 ln p(x j |zω j , α j ), where ln p(z b |π b ) = n i=1 G b g=1 z ibg ln π bg and ln p(x j |z b , α j ) = n i=1 G b g=1 z ibg ln p(x ij |α jg ).

  Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Model collection M M = {m : ωj ≤ Bmax and G b ≤ Gmax; j = 1, . . . , d; b = 1, . . . , Bmax}, where Bmax is the maximum number of blocks and Gmax is the maximum number of components within block. Model selection Model selection often achieved by searching the model m maximizing the BIC criterion which is defined by BIC(m) = max θm pen(θm|m, x) where pen(θm|m, x) = (θm|m, x) -Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Penalised completed likelihood pen(θm|m, x, z) = (θm|m, x, z)x j |zω j , α j ) -ν j Gω j 2 ln n, Consequence For z b fixed: possibility to re-affect each variable individually to the most accurate block. Thus computationnaly attractive 15/24 Vincent VANDEWALLE 1,2 Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Combinatorial model selection through a modified the EM algorithm for B and (G 1 , . . . , G B ) fixed : choice of m ⇔ choice of ω

  Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Integrated complete-data likelihood p(x, z|m) = p(x, z|m, θ)p(θ|m)dθ. Assumptions Independence between the prior distributions Standard conjugate priors Closed form of the complete-data integrated likelihood MICL (maximum integrated complete-data likelihood) criterion MICL(m) = ln p(x, z m |m) with z m = arg max zm ln p(x, z|m). Thus (m , z m ) = arg max (m,zm) ln p(x, z|m). Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood Motivations Criteria based on the integrated complete-data likelihood are popular for model-based clustering Take into account the clustering purpose : model the data distribution and provide well-separated components Optimisation of MICL over (z, ω) for B and G1, . . . , GB fixed Starting at the initial value ω [0] , each ωj is uniformly sampled among {1, . . . , B}, the algorithm at iteration [r] is Partition step: find z [r] b such that for all b = 1, . . . , B p(x

  Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Conclusion Proposition of model-based clustering with several class variables, each one explaining the heterogeneity of a block of variables: Find groups of variables producing the same clustering of the individuals Interpret the clustering produced for each group of variables Model search performed simultaneously with parameters estimation Proposed model can be used in the heterogeneous data settings Perspectives Consider the semi-supervised setting in the multi-partition framework for new partitions discovery: i. e. z1 known and z2 unknown.

  Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering

	Vincent VANDEWALLE 1,2 Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives	Maximum likelihood inference Penalized observed-data likelihood Integrated complete-data likelihood	12/24

  Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering

	Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives Data description NBA teams for the season 2016/2017 described by 16 numerical variables 1 NBA team data Wine data total minutes played (min) field goals made rate (fgmr) field goals attempted (fga) three-pointers made rate (3pmr) three-pointer attempted (3pa) free throws made (ftm) free throw attempted (fta) offensive rebounds (or) total rebounds (tr) assists (as) steals (st) turnovers (to) blocks (bk) personal fouls (pf) technical fouls (tc) points (pts) Model selection with BIC B BIC Time (s) Block G variables 3 -1909 170 1 2 fgmr, 3pmr, pf 2 2 fga, 3pa, tr, as, st, to, pts 3 1 min, ftmr, fta, orr, bk, tc Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives NBA team data Wine data Model parameters First block: offensive vs defensive teams Three features: field goals made rate, three points made rate and personal fouls fmgr 3pmr pf offensive teams: (π11 = 0.57) mean 0.468 0.371 1628.042 high shooting ability sd 0.009 0.010 76.138 low personal fouls defensive teams: (π12 = 0.43) mean 0.446 0.342 1635.503 low shooting ability sd 0.005 0.010 173.701 high personal fouls Second block: two better statistics teams for general performances vs others Seven features: field goal attempted, 3 points attempted, total rebounds, assists, steals, turnovers and points fga 3pa tr as st to GS Warriors: (π11 = 0.07) mean 7144 2934 3642 2281 728 1186 Houston Rockets sd 3 372 3 211 59 3 Other: (π11 = 0.93) mean 6992 2162 3564 1825 626 1091 teams sd 183 262 141 131 45 105 Parameters estimation and model selection Numerical experiments Conclusion and perspectives NBA team data Wine data Data description Data 1 27 chemical and physical properties of three types of Italian wines: Barolo, Grignolino, Barbera Data collected during the time period of 1970-1979 Models selected by BIC B BIC Time Block G ARI 1 -6025.00 30 1 4 0.78 2 -5947.88 280 1 3 0.87 2 4 0.16 3 -5921.42 1590 1 4 0.74 2 4 0.20 3 2 0.02 4 -5918.06 6065 1 4 0.75 2 2 0.21 3 3 0.02 4 2 0.00 Extension of the variable selection to variable clustering Parameters estimation and model selection Numerical experiments Conclusion and perspectives NBA team data Wine data Model interpretation Block 1: the type of wines (ARI=0.75) 19 variables: Alcohol, Sugar-free Extract, Tartaric Acid, Uronic Acids, Alcalinity of Ash, Calcium, Magnesium, Phosphate, Total Phenols, Flavanoids, Non-flavanoid Phenols, Proanthocyanins, Color Intensity, Hue, OD280/OD315 of Diluted Wines, OD280/OD315 of 9481 pts 22 8600 259 Flavanoids, Glycerol, 2-3-Butanediol, Proline Barolo Grignolino Barbera Class 1 0 45 0 Class 2 0 5 48 Class 3 58 1 0 Class 4 1 20 0 Block 2: year of production 4 variables: Fixed Acidity, Malic Acid, pH, Total Nitrogen Year 1970 1971 1972 1973 1974 1975 1976 1978 Class 1 8 25 4 27 31 3 1 0 Class 2 1 3 3 2 14 6 16 29 Vincent VANDEWALLE 1,2 Extension of the variable selection to variable clustering Vincent VANDEWALLE 1,2	21/24 23/24

Third block: Six features detected as irrelevant for clustering.

1 available in the package pgmm 22/24 Vincent VANDEWALLE 1,2

http://www.dougstats.com/16-17RD.Team.Opp.txt 19/24 Vincent VANDEWALLE 1,2 Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering

http://www.dougstats.com/16-17RD.Team.Opp.txt

20/24Vincent VANDEWALLE 1,2 Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering

24/24Vincent VANDEWALLE 1,2 Joint work with Matthieu MARBAC 3 A tractable Multi-Partitions Clustering


Vincent Vandewalle

Cristian Preda

Sophie Dabo

Clustering of spatial function data: what focus and how? Focus Clustering of the times Clustering of the spatial regions Clustering of both times and spatial regions Global explicative model of the process without particular clustering structure 6/20

2,3 and Sophie DABO 2,4 Clustering spatial functional data 6/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data Clustering of spatial function data: what focus and how? Focus Clustering of the times Clustering of the spatial regions Clustering of both times and spatial regions Global explicative model of the process without particular clustering structure Method Geometric based approach: distance between curves weighted by the spatial proximity Model based approach: integrate the spatial aspect in the model 6/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data 7/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data Application on Ozone concentration data Conclusion and perspectives References Conclusion Sparse model to take into account spatial dependency Smoothing of the clustering based on spatial location 18/20

• • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Basic notations for functional data X = Xs, s ∈ R N , a measurable spatial process N ≥ 1, defined on some probability space (Ω, A, P)

X observed on some spatial region

I ⊆ R N of cardinal n, I = {s1, . . . , sn}, si ∈ R N , i = 1 . . . n.
In each location s ∈ I, the random variables Xs are valued in a metric space (E, d) of eventually infinite dimension and are locally identically distributed.

d(., .) is some measure of proximity, for instance a metric or a semi-metric u and v close ⇒ Xu and Xv have same or similar distributions, less restrictive that strict stationarity

In Functional Data Analysis (FDA): E is a space of functions, typically the space of squared integrable functions defined on some finite interval T = [0, T ], T > 0.

S the set of the n curves, S = {Xs, s ∈ I}. Let f the probability distribution of X and fg the probability distribution of X given Z = g.

The mixture model pdf:

f (x) = G g=1 πgfg(x),
where πg = P (Z = g) is the prior probability of cluster g.

Extension to the spatial setting: involving the location s into the priors

f (x|s) = G g=1 πg(s)fg(x),
the distribution of observations within the cluster g is independent of the location ⇒ spatial dependency captured by the priors πg(s). PREDA 2,3 andSophie DABO 2,4 Clustering spatial functional data Density for functional data?

X|Z = g ∼ Gaussian process
The pseudo-density defined in Delaigle and Hall ( 2010) is considered:

f (qg ) g (x; θ k ) = qg j=1 fgj (cgj (x); λgj ),
f gj is the p.d.f of j-th principal component C gj of X within the cluster g C gj (j = 1, . . . , qg) independent Gaussian zero-mean with variance equal to the eigen values λ gj of the covariance operator of X qg and d need to be defined. The EM algorithm

Proposal if X(t) =

Start with an initial random partition of data S into G clusters. Let θ (h) be the estimated parameter value at iteration h ≥ 0 of the algorithm

E step t (h+1) g (s) = E θ (h) [Zg(s)|s] = πg(s; β (h) )f (qg ) g (xs; θ (h) g ) G =1 π (s; β (h) )f (q ) (xs; θ (h) )
. Notice that β (h+1) is obtained as solution of a weighted logistic regression.

M step θ

For homoscedastic models just a simple modification of the update θ (h+1) g

. See also [START_REF] Bouveyron | High-dimensional data clustering[END_REF] for more details. q q q q q q q q q q q q q q q q q q q q qq q q q q q q q 30 33 36 Clustering spatial functional data

  accurate parametric model on fg et πg 8/20 Vincent VANDEWALLE 1,2 Joint work with Cristian

f

  gj (c gj (x), λg), λg mean of the eigen values λ gj (j = qg + 1, . . . , d) of the covariance operator of X f (qg ) g : true density if X ∈ span({φ1, . . . , φ d }) 10/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data Introduction Model based clustering clustering of spatial functional data Application on Ozone concentration data Conclusion and perspectives References Parameters estimation (2/2)

  Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data q = 18 principal components retained good separation between East cities from the West cities good separation between red and blue curves in average West cities have higher pollution than Est cities

  Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4 Clustering spatial functional data Introduction Model based clustering clustering of spatial functional data Application on Ozone concentration data Conclusion and perspectives References Conclusion Sparse model to take into account spatial dependency Smoothing of the clustering based on spatial location Perspectives Avoid prior smoothing of the data ⇒ perform the selection of the basis and other smoothing parameters based on the model Propose more flexible model for π k (s) 18/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4

  Introduction Model based clustering clustering of spatial functional data Application on Ozone concentration data Conclusion and perspectives ReferencesRunning example: US ozone data (https://www.epa.gov/outdoor-air-quality-data)

	106 monitoring stations of United
	States in 2015		0.06
	Ozone recorded hourly from July 19 at 12 am to July 20 at 11 pm Expansion of the discrete data (48	Ozone concentration (g/m3)	0.02 0.04
	data at each station) in terms of
	25 Fourier basis functions		0.00
			0	10	20	30	40
					Time

  •

  •

  Joint work withCristian PREDA 2,3 and Sophie DABO 2,4

	Introduction Introduction Introduction
	Model based clustering clustering of spatial functional data Model based clustering clustering of spatial functional data Model based clustering clustering of spatial functional data
	Application on Ozone concentration data Application on Ozone concentration data Application on Ozone concentration data
	Conclusion and perspectives Conclusion and perspectives Conclusion and perspectives
	References References References
	Clustering of spatial function data: what focus and how? Clustering of spatial function data: what focus and how? The proposed model: main idea
	Standard mixture model
	Let Z a latent categorical random variable defining G clusters
	Focus Focus
	Clustering of the times Clustering of the times
	Clustering of the spatial regions Clustering of the spatial regions
	Clustering of both times and spatial regions Clustering of both times and spatial regions
	Global explicative model of the process without particular clustering Global explicative model of the process without particular clustering
	structure structure
	Method
	Geometric based approach: distance between curves weighted by the
	spatial proximity
	Model based approach: integrate the spatial aspect in the model
	6/20
	Vincent VANDEWALLE 1,2 Clustering spatial functional data

5/20

Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4

Clustering spatial functional data

 PREDA 2,3 andSophie DABO 2,4Clustering spatial functional data

The proposed model: parametrization of the functional part Density for functional data?

The pseudo-density defined in [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF] is considered:

f gj is the p.d.f of j-th principal component C gj of X within the cluster g C gj (j = 1, . . . , qg) independent Gaussian zero-mean with variance equal to the eigen values λ gj of the covariance operator of X qg and d need to be defined.

10/20 Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 and Sophie DABO 2,4

Clustering spatial functional data

Remarks

Equivalence with parsimonious high dimensional model [START_REF] Bouveyron | High-dimensional data clustering[END_REF] d is chosen as the dimension of the basis which has been used to perform the smoothing of the data C kj can be obtained by performing PCA on the expansion coefficients of X in the metric M given by the inner product of the basis functions learning data: expansion coefficients multiplied by M 1/2 ⇔ learning a parsimonious high dimensional model (see [START_REF] Bouveyron | High-dimensional data clustering[END_REF])

Completed log-likelihood

Zg(s) the indicator random variable for the cluster g at location s.

ν G is the number of parameters of the model: spatial mixing proportions, center means, principal scores and variances and n = |I|.

Selection of G when qg (g = 1, . . . , G) are unknown Modified M step which tries to maximize the conditional expectation of the BIC criterion:

(q

, where νq g is the number of parameters required for the model with qg principal components. Vincent VANDEWALLE 1,2 Joint work with Cristian PREDA 2,3 andSophie DABO 2,4 Clustering spatial functional data

14/20

Selection of G when qg (g = 1, . . . , G) are known Maximize the Bayesian Information Criterion (BIC) criterion:

ν G is the number of parameters of the model: spatial mixing proportions, center means, principal scores and variances and n = |I|.

Selection of G when qg (g = 1, . . . , G) are unknown Modified M step which tries to maximize the conditional expectation of the BIC criterion:

(q

, where νq g is the number of parameters required for the model with qg principal components.

Remark

Also possible to consider the homoscedastic setting. we see that it is the six first hours that well separate cluster 1 (North) from cluster 3 (South). q q q q q q q q q q q q q q qq 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 30 33 36 Clustering spatial functional data


M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 9/25 1 Block Low-Rank clustering techniques M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 10/25 Clustering techniques: existing solutions • Supernode reordering techniques (supernodal) 5 4 4 3 9 2 2 3 3 9 (a) Symbolic factorization (b) First separator clustering 8 × 8 × 8 Laplacian partitioned using SCOTCH and Reordering clustering on the first separator M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 11/25

Symbolic Factorization

General approach 1. Build a partition with the nested dissection process 2. Compress information on data blocks 3. Compute the block elimination tree using the block quotient graph Adjacency graph (G). Elimination tree (T).

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Factorized matrix (L).

From full-rank to low-rank

Full-rank algorithm

Scalability of Minimal Memory strategy on Laplacian matrices

1 0 3 3 0 3 5 0 3 7 0 3 9 0 3 1 1 0 3 1 3 0 3 1 5 0 3 1 7 0 3 1 9 0 3 2 1 0 3 2 3 0 3 2 5 0 3 2 7 0 3 2 9 0 3 3 1 0 3 3 3 0 3 Size 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3

Size of the factors (GB)

Full-rank 1 0 3 3 0 3 5 0 3 7 0 3 9 0 3 1 1 0 3 1 3 0 3 1 5 0 3 1 7 0 3 1 9 0 3 2 1 0 3 2 3 0 3 2 5 0 3 2 7 0 3 2 9 0 3 3 1 0 3 3 3 0 Up to what limit?

τ =10 -12 τ =10 -8 τ =10 -4
• 96 cores machine with 3TB of memory

• 100 × 100 × 100K mesh • Full-rank would have required 11T B of memory • With τ = 10 -4 :
Only 2T B of memory for the coefficients Less than 6 hours to factorize 10 -8 solution with 15 iterations of iterative refinement

• Supernode reordering techniques (supernodal)

• k-way partitioning (dense, multifrontal)

New heuristic based on projections

ComputeProjections( nblvl, distp, width )

Lets consider a supernode and the nblvl levels below in the elimination tree, with the constraint that nodes are issued from nested dissection, and not from minimum-fill as it happens at the bottom of the elimination tree.

• Each vertex of the separator being clustered is pre-selected iff it is a neighbor to one of the children at a distance distp • The trace is then enlarged with neighbors within the separator at a distance width, similarly to BDLR.

ComputeProjections( nblvl, distp, width )

Lets consider a supernode and the nblvl levels below in the elimination tree, with the constraint that nodes are issued from nested dissection, and not from minimum-fill as it happens at the bottom of the elimination tree.

• Each vertex of the separator being clustered is pre-selected iff it is a neighbor to one of the children at a distance distp • The trace is then enlarged with neighbors within the separator at a distance width, similarly to BDLR.

How to control the number of pre-selected vertices?

• Ensure that this number if low wrt to the total number of vertices of a separator: Reordering() 10:

#selected < α √ n • Ensure that it
End If 11: End For 12: Reordering() of pre-selected vertices Projections are performed on separators of size larger than 16 × M AX BLOCKSIZE, to obtain at least four large connected subparts that will be split into four parts using kway.

Parameters

BLR solver

• Tolerance τ : absolute parameter (normalized for each block) Performance profile of the factorization time -Minimal Memory 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 q q q q q q 1e-8 1e-12 Performance profile of the memory consumption -Minimal Memory 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 q q q q q q 1e-8 1e-12 Performance profile of the factorization time -Just-In-Time 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 q q q q q q 1e-8 1e-12 Performance profile of the memory consumption -Just-In-Time 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 q q q q q q 1e-8

•
Factorization time 0% 10% 20% 30% 40% 0% 10% 20% 30%
Memory consumption 0% 10% 20% 30% 0% 10% 20% 30%

Conclusion

Ordering strategies

• Both reordering and k-way are limited to form suitable clusters reordering does not allow a good compression of separators k-way compresses well the separators, but does not take into account updates

• Pre-selection of some unknowns allows to obtain well-separated clusters and to exhibit non-compressible operations • We demonstrated it is suitable for reducing time-to-solution with only a small memory consumption increase

•

  128GB of memory • Up to 8M unknowns in full-rank • Up to 36M unknowns with τ = 10 -4 M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 7/25

  separator clustering 8 × 8 × 8 Laplacian partitioned using SCOTCH and k-way clustering on the first separator M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 11/25 How to take benefit from both algorithms? Pre-selection to detect uncompressible vertices • Pre-select unknowns that represent close interactions and will not be compressible • Idea of Boundary Distance Low-Rank (BDLR) introduced by Darve et al. • For instance, pre-select interactions between a father and its direct children • Hypothesis similar to ILU(k) factorizations

  Compression method is RRQR • Compression strategy: Minimal Memory or Just-In-Time • Blocking sizes: between 128 and 256 in following experiments Pre-selection • Graphs are reconnected using distance 1 in the Halo • Selected vertices are obtained with ComputeProjections( 3, 1, 1 ) • Projections are performed on supernodes of size larger than 16 times the blocking size parameters. • K-way is configured to get partition fitting the blocking size parameters • Connected components smaller than the minimal size for low-rank compression are merged into the pre-selected vertices set. M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 19/25

Future works •

 works Align separators to form clusters during the nested dissection process • Preliminary work using fixed vertices M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct gitlab.inria.fr/solverstack/pastix • Support shared memory with different schedulers: sequential static scheduler PARSEC/STARPU runtime systems with experimental GPU support • Low-rank support • Cholesky, LDL and LU factorizations • GMRES, CG, BiCG iterative refinements Thank you. M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 25/25

25 From full-rank to low-rank

	Just-In-Time algorithm
	• For each column-block /
	supernode:
	1. Factorize the diagonal block
	(POTRF/GETRF) • For each column-block / 2. Compress off-diagonal supernode: blocks 1. Factorize the diagonal block 3. Solve low-rank off-diagonal
	(POTRF/GETRF) blocks in the current column 2. Solve off-diagonal blocks in (TRSM)
	the current column (TRSM) 4. Update the trailing matrix 3. Update the trailing matrix (LR2FR)
	(GEMM)

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 4/

25 From full-rank to low-rank

	Minimal Memory algorithm
	• Compress all off-diagonal
	blocks
	• For each column-block /
	supernode:
	1. Factorize the diagonal block
	(POTRF/GETRF)
	2. Solve low-rank off-diagonal
	blocks in the current column
	(TRSM)
	3. Update the trailing matrix
	(LR2LR)
	• Reduce the time to solution • Reduce the memory peak of
	the solver • Memory peak is as large as full-rank • LR2LR updates may be too
	costly for small cases

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 4/M.

Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 4/25 Performance of RRQR/Just-In-Time

		1.6						τ =10 -4				τ =10 -8		τ =10 -12
		1.4
	Time BLR / Time PaStiX	0.4 0.6 0.8 1.0 1.2	1.0e-05	6.1e-10	2.3e-14	1.1e-05	6.5e-10	1.6e-14	2.6e-04	2.4e-10	2.1e-15	6.9e-05	1.4e-08	9.1e-13	4.2e-04	3.7e-10	7.4e-15	5.7e-05	1.9e-08	1.8e-12
		0.2
		0.0	lap120	atmosmodj		audi		Geo1438		Hook		Serena

M.

Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 5/25 Performance of RRQR/Minimal Memory

		4.5
		4.0						τ =10 -4				τ =10 -8		τ =10 -12
	Time BLR / Time PaStiX	1.0 1.5 2.0 2.5 3.0 3.5	7.0e-05	4.0e-08	8.9e-13	8.7e-05	1.8e-08	7.1e-13	9.3e-04	4.8e-06	4.3e-11	5.9e-04	2.1e-07	4.5e-11	1.7e-04	1.3e-07	2.4e-12	3.8e-04	4.4e-07	2.4e-11	Factorization time
		0.5
		0.0	lap120	atmosmodj		audi		Geo1438		Hook		Serena
		1.2
								τ =10 -4				τ =10 -8		τ =10 -12
	Memory BLR / Memory PaStiX	0.4 0.6 0.8 1.0	7.0e-05	4.0e-08	8.9e-13	8.7e-05	1.8e-08	7.1e-13	9.3e-04	4.8e-06	4.3e-11	5.9e-04	2.1e-07	4.5e-11	1.7e-04	1.3e-07	2.4e-12	3.8e-04	4.4e-07	2.4e-11	Memory footprint
		0.2
		0.0	lap120	atmosmodj		audi		Geo1438		Hook		Serena

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 6/25

Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 22/25

		30
	Number of matrices	10 20									Factorization time
		0
		0%	10%	20%	30%	40%	0%	10%	20%	30%	40%
					Overhead percentage wrt the optimal
					q	Projection	K-way	Reordering
		M. Faverge -

Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 23/25 Conclusion

	1e-12 10% Overhead percentage wrt the optimal 15% 20% 0% 5% q Projection K-way Reordering • Both reordering and k-way are limited to form suitable clusters 0% 5% 10% 15% 20% 0 10 20 30 Number of matrices reordering does not allow a good compression of separators k-way compresses well the separators, but does not take into account updates • Pre-selection of some unknowns allows to obtain well-separated clusters Memory consumption and to exhibit non-compressible operations • We demonstrated it is suitable for reducing time-to-solution with only a M. Faverge -Ordering strategies small memory consumption increase

M. Faverge -Supernodes

ordering to enhance Block Low-Rank compression in sparse direct solvers 24/25 Conclusion

	Ordering strategies
	• Both reordering and k-way are limited to form suitable clusters
	reordering does not allow a good compression of separators
	k-way compresses well the separators, but does not take into account
	updates
	• Pre-selection of some unknowns allows to obtain well-separated clusters
	and to exhibit non-compressible operations
	• We demonstrated it is suitable for reducing time-to-solution with only a
	small memory consumption increase
	Future works
	• Align separators to form clusters
	during the nested dissection process
	• Preliminary work using fixed vertices

M.

Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 24/25

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 12/25

M. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers 13/25

Numerical ExperimentsM. Faverge -Supernodes ordering to enhance Block Low-Rank compression in sparse direct solvers

Cost of applying a m × n update to a M × N block.