Nicolas Meyer
email: nicolas.meyer@umontpellier.fr

Olivier Wintenberger
email: olivier.wintenberger@sorbonne-universite.fr

Multivariate sparse clustering for extremes

Keywords: Euclidean projection onto the simplex, model selection, multivariate extremes, regular variation

Identifying directions where extreme events occur is a major challenge in multivariate extreme value analysis. In this paper, we use the concept of sparse regular variation introduced by Meyer and Wintenberger (2021) to infer the tail dependence of a random vector X. This approach relies on the Euclidean projection onto the simplex which better exhibits the sparsity structure of the tail of X than the standard methods. Our procedure based on a rigorous methodology aims at capturing clusters of extremal coordinates of X. It also includes the identification of the threshold above which the values taken by X are considered as extreme. We provide an efficient and scalable algorithm called MUSCLE and apply it on numerical examples to highlight the relevance of our findings. Finally we illustrate our approach with financial return data.

Introduction

The aim of this article is to study the tail dependence of a random vector X ∈ R d + with continuous marginals. In this context it is customary to assume that X is regularly varying (see e.g. [START_REF] Resnick | Extreme Values, Regular Variation and Point Processes[END_REF], [START_REF] Resnick | Heavy-Tail Phenomena: Probabilistic and Statistical Modeling[END_REF], Hult and Lindskog (2006)), i.e. that there exist a n → ∞ and a non-zero Radon measure µ on the Borel σ-field of R d + \ {0} such that

nP(a -1 n X ∈ •) v → µ(•) , n → ∞ , (1.1)
where v → denotes the vague convergence in the space of nonnegative Radon measures on [0, ∞] d \ {0}. The limit measure µ is called the tail measure of the regularly varying vector X. It satisfies the homogeneity property µ(tB) = t -α µ(B), for any set B in R d + \ {0} and any t > 0. The parameter α is called the tail index of X. It highlights the intensity of the extremes. The smaller this index is, the heaviest the tail of X is likely to be.

It is often more convenient to decompose the former convergence into a radial and an angular part (see for instance [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF], Section 8.2.3): the regular variation property is equivalent to the convergence

+ := {x ∈ [0, ∞) d : x 1 + • • • + x d = 1}.
In order to study the support of the spectral measure we partition the simplex in terms of the nullity of some coordinates [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF], [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF], [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF]). For β ⊂ {1, . . . , d} the subspace C β is defined as

C β = {x ∈ S d-1 + : x i > 0 for i ∈ β, x i = 0 for i / ∈ β} . (1.3)
This partition highlights the extremal structure of X. For a given β ⊂ {1, . . . , d} the inequality P(Θ ∈ C β ) > 0 implies that the marginals X j , j ∈ β, are likely to take simultaneously large values while the ones for j ∈ β c are of smaller order. Hence the identification of clusters of directions β which concentrate the mass of the spectral measure brings out groups of coordinates which can be large together.

Highlighting such groups is at the core of several recent papers on multivariate extremes, all of them relying on some hyperparameters [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF], [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF], [START_REF] Chiapino | Identifying groups of variables with the potential of being large simultaneously[END_REF], [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF]). This approach faces a crucial issue, namely the difference of support between Θ and X/|X|. Indeed, the spectral measure is likely to place mass on low-dimensional subspaces C β , β = {1, . . . , d}. We say that this measure is sparse when the number of coordinates in the associated clusters β is small.

Conversely, the distribution of the self-normalized vector X/|X| only concentrates on the subset C {1,...,d} since X has continuous marginals.

All the existing approaches proposed in the literature rely on nonstandard regular variation for which α = 1 and all marginals are tail equivalent, possibly after a standardization.

However, sparsity arises all the more for standard regular variation (1.2). In this case, it is possible that the marginals of X are not tail equivalent so that the support of the spectral measure is included in S r-1 + for r d. This is the approach we use in this article. For a comparison of standard and nonstandard regular variation we refer to [START_REF] Resnick | Heavy-Tail Phenomena: Probabilistic and Statistical Modeling[END_REF], Section 6.5.6.

In this article we provide a method which highlights the sparsity of the tail structure by exhibiting sparse clusters of extremal directions. By sparse clusters we mean groups of coordinates β which contain a reduced number of directions compared to d. We refer to this method as sparse clustering. The statistical procedure we propose to achieve this clustering relies on the framework of [START_REF] Meyer | Sparse regular variation[END_REF] which allows to circumvent the estimation's issue that arises with the spectral measure. The angular component X/|X| in (1.2) is replaced by π(X/t), where π denotes the Euclidean projection onto S d-1 + [START_REF] Duchi | Efficient projections onto the 1 -ball for learning in high dimensions[END_REF], [START_REF] Kyrillidis | Sparse projections onto the simplex[END_REF], [START_REF] Condat | Fast projection onto the simplex and the 1 ball[END_REF]). This substitution leads to the concept of sparse regular variation. A random vector X is said to be sparsely regularly varying if

P((|X|/t, π(X/t)) ∈ • | |X|> t) w → P((Y, Z) ∈ •) , t → ∞ , (1.4)
where Z is a random vector on the simplex S d-1 [START_REF] Meyer | Sparse regular variation[END_REF] proved that under mild assumptions both concepts of regular variation (1.2) and (1.4) are equivalent (see Theorem 1 in their article). In particular, the

+ and P(Y > y) = y -α , y > 1.
relation Z = π(Y Θ) holds.
Similarly to the existing approaches with Θ, we are willing to capture the tail dependence of X via the identification of the clusters β which satisfy P(Z ∈ C β ) > 0. We call such β's the extremal clusters. They can be identified via the study of π(X/t) since the convergence P(π(X/t) ∈ C β | |X|> t) → P(Z ∈ C β ) holds for any β ⊂ {1, . . . , d} (see [START_REF] Meyer | Sparse regular variation[END_REF], Proposition 2). This encourages to consider for any β the quantity

T n,k (β) = k j=1 1{π(X (j) /|X (k+1) |) ∈ C β } , (1.5)
where X 1 , . . . , X n is a sample of iid sparsely regularly varying random vectors, k = k n is an intermediate sequence called level which satisfies k → ∞ and k/n → 0, and X (j) denotes the observation with j-th largest norm:

|X (1) |≥ • • • ≥ |X (n) |.
It turns out that the number of positive T n,k (β) often overestimates the total number of extremal clusters. We call the clusters which satisfy T n,k (β) > 0 and P(Z ∈ C β ) = 0 the biased clusters. The approach we propose to reduce this bias relies on model selection. It consists in fitting a multinomial model to the data and to compare the Kullback-Leibler divergence between the data and this theoretical model. We obtain a minimization criterion based on a penalized likelihood similarly to Akaike's criterion [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]). This approach provides a way to select the appropriate number of extremal clusters for a given level k. This is the first step of our procedure, which we call the bias selection.

The second step then consists in extending the procedure in order to automatically select an appropriate level k. We call this step the level selection. Several authors have pointed out that choosing a reasonable level, or equivalently a reasonable threshold above which the data are considered as extreme, is a challenging task in practice. This issue is tackled in a few articles [START_REF] Stȃricȃ | Multivariate extremes for models with constant conditional correlations[END_REF], [START_REF] Abdous | Non-parametric estimators of multivariate extreme dependence functions[END_REF], [START_REF] Kiriliouk | Peaks over thresholds modeling with multivariate generalized Pareto distributions[END_REF], [START_REF] Wan | Threshold selection for multivariate heavy-tailed data[END_REF], see also the review on marginals threshold selection by [START_REF] Caiero | Threshold selection in extreme value analysis[END_REF]). It turns out that in the sparse regular variation framework the choice of such a level and the identification of the extremal clusters are closely related. Therefore our approach consists in extending the bias selection by including k as a parameter to tune. Since Akaike's procedure only holds for a constant sample size we have to adapt the standard approach to an extreme setting where the number of extremes varies. Therefore we include the non-extreme values in the model and separate the data into an extreme group and a non-extreme one. The procedure then provides a level k for which this separation is reasonable. To the best of our knowledge, our work is the first one which simultaneously tackles this issue with the study of tail dependence.

Outline of the paper The paper is organized as follows. Section 2 introduces the theoretical background on sparse regular variation and level selection that is needed throughout the paper. In Section 3 we introduce the statistical framework of our method and establish asymptotic results for the estimators of the probabilities P(Z ∈ C β ). Section 4 details the methodology of our approach. We develop the two steps of the model selection, the bias selection and the level selection. In Section 5 we illustrate our findings on numerical results and compare our approach with the existing procedures proposed by [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] and [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF]. Finally we illustrate our approach on financial data in Section 6.

The proofs are given in the Supplementary Material.

Sparse regular variation

We consider a sparsely regularly varying random vector X ∈ R d + as defined in (1.4) and focus on its angular component π(X/t):

P (π(X/t) ∈ • | |X|> t) w → P(Z ∈ •) , t → ∞ .
(2.1)

The orthogonal projection on the simplex enjoys many sparsity properties which justifies its use to study high-dimensional data. The vector π(X/t) may put mass in every subspace C β even if X is almost surely positive. This is a key difference with the self-normalized vector X/|X| which shares the same sparsity properties as X, and therefore always concentrates on the interior C {1,...,d} of the simplex.

Remark 1. Our statistical methodology exhibits the choice of a level k which corresponds to the number of vectors among a sample X 1 , . . . , X n which are considered as extreme. This is achieved by studying also the n-k non-extreme vectors. In terms of the convergence (2.1) the latter vectors correspond to vectors whose norm is below the threshold t = |X (k+1) |. In order to propose a consistent methodology based on these non-extreme vectors we need to slightly modify the projection and to consider π as the Euclidean projection onto the unit

positive 1 -ball B d + = {x ∈ R d + : x 1 + . . . + x d ≤ 1}.
It does not change the theory of sparse regular variation since projecting onto the sphere or the ball is equivalent for vectors with norm larger than 1. The only difference is that a vector v such that |v|< 1 now satisfies

π(v) = v.
Our aim is to infer the distribution of the angular vector Z in order to identify the extremal directions of X. This is achieved by focusing on the probabilities p * (β) :

= P(Z ∈ C β ) for β ∈ P * d .
We define the set of extremal clusters

S * (Z) := {β : p * (β) > 0} , (2.2)
and denote by s * its cardinality. [START_REF] Meyer | Sparse regular variation[END_REF] proved that for any β we have the convergence

P(π(X/t) ∈ C β | |X|> t) → p * (β) , t → ∞ . (2.3)
This convergence allows one to study the behavior of Z on the subsets C β via the one of π(X/t). The aim of this paper is to build a statistical procedure to identify the extremal clusters β ∈ S * (Z).

Example 1 (Discrete spectral measure). For β ∈ P * d , we denote by e(β) the sum j∈β e j so that the vector e(β)/|β| belongs to the simplex S d-1 + (recall that |β| corresponds to the length of the cluster β). We consider the following family of discrete distributions on the simplex:

β∈P * d c(β) δ e(β)/|β| , (2.4)
where (c(β)) β is a probability vector on R If we choose c(β) = 0 for all β's except the ones of length 1 then the spectral measure becomes d j=1 c j δ e j , (c j ) 1≤j≤d ∈ S d-1 + . This corresponds to asymptotic independence (see e.g. [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF]Tawn (2004), de Haan and[START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Section 6.2).

If Θ places mass on a subset C β then so does Z, but the converse is not true. Thus the set of clusters we identify with our method includes the usual ones on which Θ puts mass. However, the notion of maximal cluster (an extremal cluster which is not included in another extremal one) defined by [START_REF] Meyer | Sparse regular variation[END_REF] coincide for Θ and Z and links both types of clusters.

Example 2. Consider a spectral measure in dimension 2 with Θ 1 ∼ U (0, 1). Then the distribution of Z is given by Z 1 = 1 4 δ 0 + 1 2 U (0, 1)+ 1 4 δ 1 , see [START_REF] Meyer | Sparse regular variation[END_REF], Example 1. In this case the clusters {1} and {2} are extremal clusters for Z but not for Θ. The only maximal cluster for Z and Θ is C {1,2} .

Impact of the level on the sparsity structure

We briefly explain in this section how the choice of a threshold t > 0 influences the sparsity of the projected vector π(x/t) for x ∈ R d + . For t > 0, let us denote by π t the Euclidean projection onto the positive sphere {x ∈ R d + :

x 1 + • • • + x d = t}.
The relation π t (x) = tπ(x/t) implies that the sparsity structures of π t (x) and π(x/t) are the same. The number of null coordinates of the projected vector π t (x) strongly depends on the choice of t. Indeed, if t is close to |x|, then π t (x) has only non-null coordinates (as soon as x itself has non-null coordinates). On the contrary, the vector π t (x) is sparse if t |x|.

Moving on to a statistical framework, we consider a sample X 1 , . . . , X n of iid random

vectors in R d + .
It is common in extreme value theory to define a level k = k n satisfying k → ∞ and k/n → 0 (see e.g. de [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF][START_REF] Resnick | Heavy-Tail Phenomena: Probabilistic and Statistical Modeling[END_REF]). It leads to the choice of a threshold t = u n → ∞ such that

n k P(|X|> u n ) → 1 , n → ∞ . (2.5)
The level k must be seen as the number of extreme vectors used for the statistical analysis.

It is therefore natural to consider the k-largest vectors in terms of their norm, i.e. the vectors

X (1) , • • • , X (k) where |X (j) | denotes the j-th largest norm |X (1) |≥ • • • ≥ |X (n) |.
Note that since we assumed the marginals of X to be continuous, these inequalities are strict almost surely. This encourages to work with the random threshold |X (k+1) |. By Vervaat's Lemma (see Lemma 1.0.2 in de Haan and Ferreira ( 2006)), the assumption (2.5) implies that |X (k+1) |/u n converges to 1 in probability as n → ∞.

0 x 1 x 2 1 1 × × × × × × × × × × × × • u n 0 x 1 x 2 1 1 × × × × × × × × • u n
{y ∈ R d + : y 1 + • • • + y d = u n }.
Figure 1: Influence of the level k on the sparsity structure of the data. The threshold u n corresponds to the norm of the vector X (k+1) which is represented by a bullet.

A small k corresponds to a large threshold u n and vice versa. In this case only a few extreme vectors are kept for the statistical analysis and they are close to the threshold.

Thus, these vectors are projected on subsets C β with large |β|'s which means that the projected vectors are not very sparse. On the other hand, choosing a large k means choosing a low threshold u n so that we move away from the extreme region. In this case the largest vectors are projected on subsets C β with small |β|'s, i.e. the projected vectors are sparse.

We refer to Figure 1 for an illustration of these two cases. Following these remarks we have to make a balanced choice between providing a sparse structure for the data and staying in the extreme region.

Asymptotic analysis of the extremal clusters

We consider a sequence of iid sparsely regularly varying random vectors X 1 , . . . , X n with generic distribution X and angular limit vector Z. We also consider a level k satisfying k → ∞ and k/n → 0 and a threshold u n such that (2.5) is satisfied. In order to identify the set S * (Z) defined in (2.2) we provide suitable estimators for the probabilities p * (β),

β ∈ P * d .
We define the estimators

T n (x, β) = n j=1 1 {X/un∈A(x,β)} , β ∈ P * d , x > 0 , where A(x, β) = {y ∈ R d + : x|y|> 1, π(xy) ∈ C β } so that the estimator T n,k (β) defined in (1.5) satisfies T n,k (β) = T n (u n /|X (k+1) |, β). An empirical version of S * (Z) is then given by S n := {β ∈ P * d : T n,k (β) > 0} . (3.1)
We denote by ŝn the cardinality of this set. Finally, we define

p n (β) = P(π(X/u n ) ∈ C β | |X|> u n ) , β ∈ P * d .
which converges to p * (β) ∈ P * d for any β, see Equation (2.3).

3.1

The bias between S n and S * (Z)

In this section we compare the set S * (Z) with its empirical counterpart S n . We first establish the consistency of our estimator.

Proposition 1. For any

β ∈ P * d , T n,k (β) k = 1 k k j=1 1 {π(X (j) /|X (k+1) |)∈C β } → p * (β) , n → ∞ , (3.2)
in probability.

The proof relies on Proposition 2.2 of de Haan and [START_REF] De Haan | Estimating the limit distribution of multivariate extremes[END_REF]. It suffices to prove that Z does not put any mass on the boundary of {|x|> 1} ∩ π -1 (C β ), which has already been established in Proposition 2 of [START_REF] Meyer | Sparse regular variation[END_REF].

Proposition 1 implies that if p * (β) = 0, i.e. if Z does not place mass on the subset C β , then T n,k (β)/k becomes smaller and smaller as n increases. Actually as soon as the dimension d is large a lot of T n,k (β)'s are even equal to 0 since the number of extreme vectors, that is k, is far below the number of clusters, that is 2 d -1.

In order to study the bias between S n and S * (Z) we focus on the speed of convergence of P(T n,k (β) = 0), and thus on the one of [START_REF] Meyer | Sparse regular variation[END_REF] established the equivalence

P(X/u n ∈ A(x, β)) = P({y ∈ R d + : x|y|> 1, π(xy) ∈ C β }).
π(x) ∈ C β if and only if ∀i ∈ β c , ∀j ∈ β, x i ≤ |x β | |β| < x j .
In other words, all x j , j ∈ β, should be of the same order, while the x i , i ∈ β c , should be of smaller order. We set x β, i = j∈β (x j -x i ) for any i and define

C β = {x ∈ R d + : min i∈β c x β, i ≥ 0} = x ∈ R d + : j∈β (x j -max i∈β c x i ) ≥ 0 , which forms a cone of R d + .
Studying the convergence of P(X/u n ∈ A(x, β)) then boils down to studying the asymptotic behavior of X on the cone C β . Based on the theory of hidden regular variation (HRV) by [START_REF] Lindskog | Regularly varying measures on metric spaces: Hidden regular variation and hidden jumps[END_REF], we make the following assumption on X.

Assumption (HRV). For every β ∈ P * d the vector X is regularly varying on R d + \ C β with tail index α(β) and exponent measure µ β satisfying

µ β ({x ∈ R d + : max i∈β x β, i < 1, min i∈β c x β, i ≥ 1}) > 0 .
This assumption allows one to deal with the asymptotic behavior of P(T n,k (β) = 0) even when p * (β) = 0, as stated in the following lemma.

Lemma 1. Under Assumption (HRV) we have for every

β ∈ P * d , log P(T n,k (β) = 0) -kp n (β) → 1 , n → ∞ .
Lemma 1 encourages to focus on the quantity kp n (β) and to consider the set

S ∞ = {β ∈ P * d : kp n (β) → ∞ as n → ∞} . (3.3)
We denote by s ∞ its cardinality. This set contains S * (Z) so that we have the inequality

s * ≤ s ∞ . Subsequently, Lemma 1 implies that P(S ∞ ⊂ S n ) = 1 -P(∃β ∈ S ∞ , β / ∈ S n ) ≥ 1 - β∈S∞ P(T n,k (β) = 0) → 1 ,
as n → ∞. This leads to the following proposition.

Proposition 2. Under Assumption (HRV) the inclusions

S * (Z) ⊂ S ∞ ⊂ S n
hold true with probability converging to 1.

These inclusions highlight the fact that the observations T n,k (β) tend to overestimate the number of clusters β in S * (Z). They imply that we only have a "one-side bias" composed of clusters that appear empirically but which theoretically do not contain any mass. One of the main challenge of our study is the derivation of the asymptotic properties of T n,k (β)

for biased clusters β ∈ S n \ S * (Z).

By Lemma 1 the inclusion S n ⊂ S ∞ means that we only observe faces C β for which P(T n,k (β) = 0) does not decrease super-exponentially fast with n. We define the sets of admissible sequences (k n ) by

K = {(k n ) : k n → ∞, k n /n → 0, S ∞ = S n a.
s. for all n large enough} .

That K is non empty is a strong assumption equivalent to the fact that S n converges a.s.

to S ∞ for some sequence of levels (k n ). Combining the definition of K with Assumption (HRV) we can rely on the statistics T n,k (β), (k n ) ∈ K, which are non-null sufficiently often even when p n (β) → p * (β) = 0, in order to quantify the bias.

Asymptotic normality

We now establish a convergence result for the joint distribution of T n,k (β) for β ∈ S ∞ .

This is achieved via the study of the joint distribution of T n (x, β) for x ∈ [ 1 1+τ , 1 + τ ], τ > 0. Having in mind the model selection proposed in Section 4 we consider for any 0 ≤ s < r ≤ s ∞ and any disjoint clusters β 1 , . . . , β r ∈ S ∞ the vectors

T s,r n (x) = T n (x, β 1 ), . . . , T n (x, β s ), r j=s+1
T n (x, β j ) ∈ R s+1 , and

P s,r n (x) = x α(β 1 ) p n (β 1 ), . . . , x α(βs) p n (β s ), r j=s+1 x α(β j ) p n (β j ) ∈ R s+1 .
For τ > 0 we denote by ∞ ([ 1 1+τ , 1 + τ ]) the set of functions defined and bounded on

[ 1 1+τ , 1 + τ ].
Theorem 1. Let Assumption (HRV) hold. Assume that there exists (k n ) ∈ K and choose

u n such that k ∼ nP(|X|> u n ) as n → ∞.

The following convergence holds in

∞ ([ 1 1+τ , 1 + τ ]) as n → ∞: √ k Diag(P s,r n (x)) -1/2 T s,r n (x) k -E T s,r n (x) k ; (1 + τ ) -1 ≤ x ≤ 1 + τ s<r d → (N s,r ) s<r , (3.4)
where the constant limit process is identified to N s,r , a standard centered multivariate Gaussian vector in R s+1 .

2. If we assume moreover that for any β ∈ S ∞ ,

sup x∈[ 1 1+τ ,1+τ ] k p n (β) n k P(X/u n ∈ A(x, β)) -x α(β) p n (β) → 0 , n → ∞ , (3.5)
then we have

√ k Diag(P s,r n (x)) -1/2 T s,r n (x) k -P s,r n (x) ; (1 + τ ) -1 ≤ x ≤ 1 + τ s<r d → (N s,r ) s<r , (3.6) in ∞ ([ 1 1+τ , 1 + τ ]) as n → ∞.
Based on Theorem 1, we establish the asymptotic behavior of the estimators T n,k (β).

We define

T s,r n,k = T s,r n (u n /|X| (k+1) ) = T n,k (β 1 ), . . . , T n,k (β s ), r j=s+1 T n,k (β j ) ∈ R s+1 .
Proposition 3. Under the assumptions of Theorem 1, under (3.5), and under the bias

assumption √ k(p n (β) -p * (β)) → 0 , n → ∞ , β ∈ S ∞ , (3.7)
we have the convergence

√ k Diag(P s,r n (1)) -1/2 T s,r n,k k -P s,r n (1) d → (Id s+1 - √ P s,r • √ P s,r )N , n → ∞ , (3.8)
where N ∈ R s+1 is a standard centered multivariate Gaussian vector, and where P s,r is the limit vector of P s,r n (1): p n (β j ) .

P s,r = p * (β 1 ), . . . , p * (β s ),
Remark 2. The bias assumption (3.7) holds for

β ∈ S ∞ \ S * (Z) if k = o(n κ ) as n → ∞
where κ > 2(α(β)-α) 2α(β)-α for every β ∈ S ∞ \ S * (Z). We refer to the Supplementary Material for a proof.

Remark 3. For r = s ∞ the matrix Id s+1 -√ P s,r • √ P s,r is symmetric and satisfies

(Id s+1 - √ P s,r • √ P s,r ) 2 = Id s+1 -2 √ P s,r • √ P s,r + ( √ P s,r • √ P s,r ) 2 = Id s+1 - √ P s,r • √ P s,r , since √ P s,r • √ P s,r = r j=1 p * (β j ) = 1.
Therefore it corresponds to an orthogonal projection with rank s. Cochran's theorem then ensures that the 2 -norm of the vector (Id s+1 -√ P s,r √ P s,r )N follows a chi-squared distribution with s degrees of freedom.

Going back to Proposition 3 we obtain the following convergence:

k s j=1 (T n,k (β j )/k -p n (β j )) 2 p n (β j ) + k [ r j=s+1 (T n,k (β j )/k -p n (β j ))] 2 r j=s+1 p n (β j ) d → ψ(s) , (3.9)
where ψ(s) follows a chi-squared distribution with s degrees of freedom. This convergence is useful to identify the parameter s in the bias selection, see Lemma 4 in the Supplementary Material.

Methodology

We develop in this section our methodology to estimate the set S * (Z). We use the same notation as in Section 3.

Bias selection

We consider the vector T n,k ∈ R 2 d -1 with components T n,k (β) whose distribution P k is multinomial with probability weights (p n (β)) β∈P * d , and adding up to k. We propose a bias selection which consists in comparing the distribution P k with the theoretical multinomial model M k with 2 d -1 outcomes adding up to k and a probability vector (p 1 , . . . , p s , p, . . . , p, 0, . . . , 0)

∈ [0, 1] 2 d -1 , with p 1 ≥ • • • ≥ p s > p and r -s components p satisfying p 1 + • • • + p s + (r -s)p = 1.
The parameters p j model the probability that Z belongs to the associated subsets C β while the parameter p models the probability that a biased cluster appears. We denote by p the vector (p 1 , . . . , p s ) ∈ B s + (0, 1). The likelihood L M k of the model M k is given by The identification of the extremal clusters β in S * (Z) is achieved by choosing the model M k which best fits the sample T n,k . Following the AIC approach of [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF], we select the multinomial model which minimizes the expectation of the Kullback-Leibler (KL) divergence (see [START_REF] Kullback | On information and sufficiency[END_REF]) between the true distribution P k and the model M k evaluated at p, where p denotes the maximum-likelihood estimator of p.

L M k (p; y) = k! 2 d -1 i=1 y i ! s i=1 p y i i r i=s+1 1 -s j=1 p j r -s y i 1 {y r+1 =•••=y 2 d -1 =0} , ( 4
Hence we consider the quantity

E[KL(P k M k )| p=p ] = E[ log L P k (T n,k )] -E[E[ log L M k (p; T n,k )]| p=p ] , (4.2)
where L P k denotes the likelihood of the distribution P k . Theorem 2 below provides an asymptotic expansion of this quantity.

Before stating this result we compute the maximum likelihood estimator of the model M k . The first components of the model M k being associated to the extremal clusters, we reorder the coordinates of the vector T n,k so that its components are ordered in the decreasing order. Hence we define T n,k,1 = max β T n,k (β) and

T n,k,j = max {T n,k (β), β ∈ P * d } \ {T n,k,1 , . . . , T n,k,j-1 } , j = 2, . . . , 2 d -1 .
The expression in (4.1) is also maximal when r corresponds to the number ŝn of clusters that appear empirically. This leads to the following expression of the log-likelihood

log L M k (p; T n,k ): log(k! ) - 2 d -1 i=1 log(T n,k,i ! ) + s i=1 T n,k,i log(p i ) + r i=s+1 T n,k,i log 1 -s j=1 p j r -s . (4.3)
The optimization of this quantity then provides the maximum likelihood estimator p ∈ R s with components pj := T n,k,j /k for 1 ≤ j ≤ s.

Theorem 2. Under the assumptions of Proposition 3 the following convergence holds:

E[KL(P k M k )| p=p ] -E[ log L P k (T n,k )] + E[log L M k (p; T n,k )] → s , n → ∞ .
Based on Theorem 2 we choose the model M k which minimizes the quantity

-log L M k (p; T n,k ) + s . (4.4)
Therefore for a given sequence (k n ) ∈ K the bias selection procedure consists in choosing the parameter ŝ(k) which minimizes this penalized log-likelihood.

(B2) For n sufficiently large, k ∈ K, there exist c, C > 0 such that cnq n ≤ k ≤ Cnq n .

Assumptions (B1) and (B2) allow one to control the bias between T n,n-T n,2 d ,j | T n,2 d and T n,k,j , and between k and nq n respectively.

The following theorem provides an asymptotic expansion of the expectation of this estimator.

Theorem 3. Under (B1), (B2) and the assumptions of Proposition 3 we have

E[KL(P n M n )| p ] = nq n E[ -log L M k (p; T n )] + s k + log(k/n) + O(nq n ) , n → ∞ .
Theorem 3 encourages to choose a level k which minimizes the penalized log-likelihood

-log L M k (p; T n ) + s k + log k n .
It turns out that the additive penalization term log(k/n) leads to numerical instability as k/n is small. To cope with this issue we upper bound it by k/n -1. The level k which minimizes the criterion is then smaller than the one that appears with log(k/n). Thus it satisfies more likely the bias assumptions (B1), (B2). So in practice we choose a level k which minimizes the following penalized log-likelihood

-log L M k (p; T n ) + s k + k n . (4.6)
Note that the two steps of our procedure are clearly identified in the penalized loglikelihood. The term -log L M k (p; T n ) + s corresponds to the bias selection and the multiplicative factor and the additional term to the level one.

Algorithm: MUltivariate Sparse CLustering for Extremes (MUSCLE)

In practice we choose a large range K of k (often between 0.5% and 10% of n) and we compute the value of (4.6) for these k and for s = 1, . . . , ŝn , where ŝn depends on the chosen level k. We choose k which minimizes the penalized log-likelihood (4.6) and then choose ŝ( k) which minimizes (4.6) for k = k. Then we define S * as the set gathering the ŝ( k) clusters corresponding to the largest T n, k(β)'s. Finally we consider the probability vector ζ defined by

ζ(β) := T n, k(β) γ∈ S * T n, k (γ)
, for β ∈ S * and 0 elsewhere, as an estimator of p * . Our procedure entails the following parameter-free algorithm called MUSCLE for MUltivariate Sparse CLustering for Extremes.

Remark 4. While our procedure leads to the choice of a unique k, we expect that this approach is not too sensitive to this choice. Therefore, it is relevant to plot the function k → ŝ(k) which provides the chosen value of s for every k ∈ K. We expect that this function is approximately constant around the chosen value k.

5 Numerical results

Overview

The aim of the numerical results is to compare the extremal clusters given by MUSCLE with the theoretical ones in S * (Z). To this end we compare the estimated probability

h(p * , ζ) = 1 √ 2 β∈P * d (p * (β) 1/2 -ζ(β) 1/2 ) 2 1/2 . ( 5.1)
The closer h(p * , ζ) is to 0, the better ζ estimates p * . In order to compare our method with some existing ones, we also compute the Hellinger distance between the true probabilities P(Θ ∈ C β ) and the estimated ones given by the algorithm called DAMEX of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] and the two methods of [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF]. We represent the mean Hellinger distance over N = 100 simulations. The parameters in the method of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] are chosen to be = 0.1, k = √ n, and p = 0.1, see the notation in their paper. Regarding the methods of [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF] we use the parameters given by the authors in Section 4.2 of their paper, i.e. we set π = 0.01, and p = 0.5 and u β to be the 0.75 quantile of observed Q values in region C β for the first method, and δ = 0.5 and u β to be the 0.85 quantile of observed Q values in region C β for the first method. We refer to [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF] for more insights on these parameters.

Remark 5. Contrary to the aforementioned methods, we recall that MUSCLE does not require any hyperparameter. This is a main advantage from a statistical and computational point of view. On the contrary, for the other methods these values could be tuned via crossvalidation. For the numerical results we do not choose this approach and keep the values fixed by the authors of the cited papers.

In the following section we develop the example of a max-mixture distribution. The code related to this article can be found at https://drive.google.com/drive/folders/ 11TvhbVMPXcSkxmdnnAySvZt64lpMKZqL?usp=sharing. Another example related to asymptotic independence is given in the Supplementary Material.

Max-mixture distribution

For any β ∈ P * d , let A β ∈ R |β| + be a random vector with standard Fréchet marginal distributions and with dependence structure given below, and let {A β : β ∈ P * d } be independent random vectors. Then the vector X = (X 1 , . . . , X d ) whose components are defined via X j = max β∈P * d :i∈β λ i, β X j ,β , with λ i, β ∈ [0, 1] and β∈P * d :i∈β λ i, β = 1, has also standard Fréchet marginal distributions and is regularly varying.

For our simulations we consider the five-dimensional example introduced by Simpson et al. ( 2020) which we recall for completeness. We consider two bivariate Gaussian copulas with correlation parameter ρ and Fréchet marginals A {1,2} and A {4,5} , and three extremevalue logistic copulas with dependence parameter α and Fréchet marginals A {1,2,3} , A {3,4,5} , and A {1,2,3,4,5} . For ρ < 1, the Gaussian copula is asymptotically independent (see Example 1 for more insights on this notion) and thus the spectral measure defined in (1.2) concentrates on the subsets C {1} , C {2} , C {4} , and C {5} . For α ∈ (0, 1) the logistic distribution is asymptotically dependent so that the spectral measure also places mass on the subsets C {1,2,3} , C {3,4,5} , and C {1,2,3,4,5} . Following [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF], we set λ {1,2} = (5, 5)/7 , λ {4,5} = (5, 5)/7 λ {1,2,3} = (1, 1, 3)/7 , λ {3,4,5} = (3, 1, 1)/7 , λ {1,2,3,4,5} = (1, 1, 1, 1, 1)/7 , so that equal mass is assigned to each of the seven aforementioned subsets. In order to compute the mass the distribution of Z assigns to every subset C β we start from the distribution of Θ and use Monte-Carlo simulation. We then compare these probabilities with their estimated ones ζ given by MUSCLE.

We run our algorithm for different values of ρ ∈ {0, 0.25, 0.5, 0.75} and α ∈ {0.1, 0.2, . . . , 0.9}.

Figures 2 and3 shows the average mean Hellinger distance for our method, the one of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF], and the two of [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF] over 100 simulations. Our method provides a mean Hellinger distance which stabilizes between 0.2 and 0.3 for all values of ρ and α. For α ≤ 0.7 the distance slightly decreases with alpha, while it increases for α ≥ 0.8. The standard deviation is quite small for α ≤ 0.7 and then increases with α. Regarding the approach of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF], the mean Hellinger distance tends to increase with α and with ρ. The smallest values is obtained for ρ ∈ {0, 0.25} and for small α. distance of 0.2. This distance is the smallest one for ρ = 0.5 and ρ = 0.75 for all α compared to the one of the three other methods. For small ρ, MUSCLE better performs for large α.

For small α the second method of [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF] provides better results than our approach, while its standard deviation is larger. It turns out that except for small α with ρ = 0 and ρ = 0.25 our algorithm better detects the extremal clusters.

6 Application to real-world data 6.1 Preprocessing for real-world data

In Section 5 we considered an example with standard Fréchet marginal distributions so that the tail index (see Equation (1.2)) of the considered vectors is equal to 1. The influence of this index on the extremal clusters has been studied on some numerical results by [START_REF] Meyer | Sparse regular variation[END_REF]. It turns out that a large tail index does not provide accurate results while a small one highlights one-dimensional clusters, see Remark 11 in their article.

A tail index of α = 1 seems to provide the best results.

For real-world data the estimation of the tail index of a sample x 1 , . . . , x n is achieved with a Hill plot [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]). It consists in plotting and to choose α as the value around which the plot stabilizes. Then, we consider the power transform x j = (x j ) α. This transformation highlights the tail structure of the data without modifying the support of the spectral measure, see [START_REF] Meyer | Sparse regular variation[END_REF], We conclude that the aforementioned 14 clusters given by MUSCLE correspond to subsets C β which gather the mass of Z. Among them, eight gather some mass of Z and are not included in larger subsets on which Z places mass. Following [START_REF] Meyer | Sparse regular variation[END_REF], Theorem 2, these maximal subsets also concentrate the mass of the spectral measure. We refer to [START_REF] Meyer | Sparse regular variation[END_REF], Section 3.2, for a discussion on maximal and non-maximal subsets. Standard approaches which hold for low-dimensional extremes can then be applied to these subsets, see [START_REF] Einmahl | Estimating a multidimensional extreme-value distribution[END_REF], [START_REF] Einmahl | Estimating the spectral measure of an extreme value distribution[END_REF], [START_REF] Einmahl | Maximum empirical likelihood estimation of the spectral measure of an extreme-value distribution[END_REF] After removing the 12 extremal components we reapply MUSCLE to obtain the dependence structure of the non-extremal portfolios. The algorithm provides a unique cluster with all 37 remaining portfolios. Hence these portfolios tend to have a dependent tail structure: their extreme variability is strongly correlated.

α(k) = 1 k k j=1 log(|x| (j) ) -log(|x| (k) ) -1 , k = 2, . . . ,

Conclusion

The statistical analysis introduced in this article provides a new approach to detect the extremal directions of a multivariate random vector X. This method relies on the notion of sparse regular variation which better highlights the tail dependence of X. Several convergence results are established in Section 3 and are used to build a rigorous statistical method based on model selection. This approach provides not only the clusters of directions on which the extremes of X gather but also a reasonable threshold above which the data are considered as extreme values. The latter issue has always been challenging and no theoretical-based procedure has been provided in a multivariate setting yet, even if it has been the subject of much attention in the literature. The choice of the directions is achieved with an AIC-type minimization whose penalization allows to reduce the number of selected subsets. Including the choice of the level k of the random threshold |X| (k) then entails multiplicative and additive penalization terms. This approach leads to the parameter-free algorithm MUSCLE whose purpose is to recover the extremal clusters of a sample of iid sparsely regularly varying random vectors X 1 , . . . , X n .

The absence of any hyperparameter is a main difference with the existing methods [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF], [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF], [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF], [START_REF] Chiapino | Identifying groups of variables with the potential of being large simultaneously[END_REF]). Another main advantage of our procedure is that it is still efficient for large d.

This follows from the expected linear-time algorithm introduced by [START_REF] Duchi | Efficient projections onto the 1 -ball for learning in high dimensions[END_REF] to compute the Euclidean projection.

The numerical experiments on max-mixture distributions provide promising results.

Our algorithm provides better results than the ones of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] and [START_REF] Simpson | Determining the dependence structure of multivariate extremes[END_REF] for ρ close to 1, or small ρ and α close to 1. Moreover the results do not vary a lot with ρ and α. Finally, the application of our algorithm on financial data highlights sparse clusters and thus reduces the dimension of the study. We obtain a sparse tail dependence structure for for the extreme variability of several industry portfolios. This reinforces the relevance of our approach for reducing the dimension in Extreme Value Theory.

→

  P((|X|/t, X/|X|) ∈ • | |X|> t) w → P((Y, Θ) ∈ •) , t → ∞ , denotes weak convergence, and where Θ is a random vector on the positive unit sphere {x ∈ [0, ∞) d : |x| = 1} independent of the random variable Y which satisfies P(Y > y) = y -α , y > 1. The random vector Θ is called the spectral vector and its distribution P(Θ ∈ •) the spectral measure. Its support indicates the directions supported by large events. The subspaces of the positive unit sphere on which the spectral vector puts mass correspond to the directions where large events are likely to appear. Note that the choice of the norm |•| in Equation (1.2) is arbitrary. In this article we choose the 1 -norm and thus focus on the simplex S d-1

  Symbols in bold such as x ∈ R d are column vectors with components denoted by x j , j ∈ {1, . . . , d}. Operations and relationships involving such vectors are meant componentwise. If x = (x 1 , . . . , x d ) ∈ R d , then Diag(x) or Diag(x 1 , . . . , x d ) denotes the diagonal matrix whose diagonal is x. We denote by Id s the identity matrix of R s . We define R d + := {x ∈ R d : x 1 ≥ 0, . . . , x d ≥ 0}, 0 := (0, . . . , 0) ∈ R d , and 1 := (1, . . . , 1) ∈ R d . For j = 1, . . . , d, e j denotes the j-th vector of the canonical basis of R d . In all the paper we denote the 1 -norm by |•|. For d ≥ 1 we denote by P d the power set of {1, . . . , d} and by P * d the set P d \ {∅}. If β ∈ P d we denote by |β| the number of coordinates in β.

  2 d -1 (see[START_REF] Segers | Max-stable models for multivariate extremes[END_REF], Example 3.3).[START_REF] Meyer | Sparse regular variation[END_REF] proved that in this case we have Z = Θ a.s. and that the family of distribution in (2.4) is the only possible discrete distributions for Z. For this type of distributions we have S * (Z) = {β : c(β) > 0}.

  (a) For k = 12 the points in the blue area are projected on the interior of the positive sphere {y ∈ R d + : y 1 + • • • + y d = u n } while the ones in the red area are projected on the edges of this sphere. (b) For k = 8 all points in the blue area are projected on the interior of the positive sphere

  β j ) = lim n→∞ p n (β 1 ), . . . , p n (β s ), r j=s+1

  .1) for any vector p ∈ B s + (0, 1) = {u ∈ R s + : u 1 + • • • + u s ≤ 1} and any y ∈ N 2 d -1 0 adding up to k, where N 0 denotes the sets of non-negative integers.

Algorithm 1 :

 1 MUltivariate Sparse CLustering for Extremes (MUSCLE) Data: A sample X 1 , . . . , X n ∈ R d + and a range of values K for the level Result: A list S * of clusters β and the associated probability vector ζ. for k ∈ K do Compute u n = |X| (k+1) the (k + 1)-th largest norm; Assign each π(X j /u n ) the subsets C β it belongs to; Compute T n,k (β) for each β ∈ P * d ; Compute the minimizer ŝ(k) which minimizes the criterion given in Equation (4.4); end Choose k which minimizes (4.6) plugging in the minimal value in (4.4); Output: S * = {the clusters β associated to the T n, k,1 , . . . , T n, k,ŝ( k) } and ζ as above. vector ζ with the theoretical one p * via the Hellinger distance

  The estimation particularly deteriorates for ρ = 0.75. Finally both methods proposed by Simpson et al. (2020) provide a mean Hellinger distance which increases with α and ρ. The second one seems to provide almost always better results than the first one. While all methods provided by Goix et al. (2017) and Simpson et al. (2020) deteriorate when ρ or α increase, our procedure provides results which stabilize around a mean Hellinger

Figure 2 :

 2 Figure 2: Mean Hellinger distance over 100 simulations for ρ = 0 (top) and ρ = 0.25 (bottom). The abbreviation SWT1 (resp. SWT2) refers to the first (resp. second) method of Simpson et al. (2020).

Figure 3 :

 3 Figure 3: Mean Hellinger distance over 100 simulations for ρ = 0.5 (top) and ρ = 0.75 (bottom). The abbreviation SWT1 (resp. SWT2) refers to the first (resp. second) method of Simpson et al. (2020).

Figure 4 :

 4 Figure 4: Evolution with respect to k of the penalized log-likelihood given in (4.6) (left) and of ŝn (k) (right) for the financial data.

Figure 5 :

 5 Figure 5: Representation of the 14 clusters and their inclusions. The abbreviations are the following ones: Softw = Computer Software, Txtls = Textiles, Hlth = Healthcare, RlEst = Real Estate, Hardw = Hardware, Chips = Electronic Equipment, Fin = Finance.

  .

		Coal, Banks, Fin
	Hlth, Softw		Coal, Banks	Steel, Coal	Gold, Coal		Hardw, Chips
	Softw	Hlth	Banks	Coal	Txtls	Gold	RlEst	Smoke

Acknowledgments We are grateful to two referees for careful reading of the paper and for useful suggestions.

Level selection

The second step of the model selection consists in considering k as a parameter which has to be estimated and tuned. It is therefore necessary to consider all observations X 1 , . . . , X n and not only the extreme ones. We consider a vector T n ∈ R 2 d such that

The last component T n,2 d corresponds to the number of non-extreme values of the sample. We assume that this vector follows a multinomial distribution P n with parameter n and probability vector

Similarly to Section 4.1 we consider a multinomial model M n with probability vector given by (q p 1 , . . . , q p s , q p , . . . , q p , 0, . . . , 0, 1 -q ) ∈ R 2 d with p 1 ≥ . . . ≥ p s > p and r -s components q p satisfying the relation p 1 +. . .+p s +(r -s )p = 1. Here q models the proportion of extreme vectors. We denote by p the vector (p 1 , . . . , p s , q ) ∈ B s + × (0, 1). We consider the Kullback-Leibler divergence between P n and M n given by

where L P n (resp. L M n ) denotes the likelihood of the distribution P n (resp. M n ). Following the same ideas as in Section 4.1 and similarly to an AIC procedure we estimate the Kullback-Leibler divergence in Equation (4.5) by the estimator KL(P n M n )| p where p denotes the maximum likelihood estimator of p .

We make the following assumptions.

(B1) For k ∈ K and β j ∈ S ∞ we have

Remark 8. It differs from the standardization discussed in the introduction for which the vectors are normalized via a rank transform.

In the following section we apply MUSCLE to financial data. An application on wind speed data can be found in the Supplementary Material.

Extreme variability for financial data

The data set we use corresponds to the value-average daily returns of 49 industry portfolios compiled and posted as part of the Kenneth French Data Library. They are available at https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html. A related study on a similar dataset has been conducted by [START_REF] Cooley | Decompositions of dependence for highdimensional extremes[END_REF]. We restrict our study to the period 1970-2019 which provides n = 12 613 observations denoted by x obs 1 , . . . x obs n ∈ R 49 . Our goal is to study the variability of these returns so that we take the componentwise absolute value x j = |x obs j | of the data. Thus, we study the non-negative vectors x 1 , . . . , x n in R d + with n = 12 613 and d = 49. Following Section 6.1, we consider the vectors x j = (x j ) α, where α = 2.99 is the Hill estimator of the sample |x 1 |, . . . , |x n |.

Following Remark 4, we plot the evolution of the estimator of the Kullback-Leibler divergence in (4.6) as a function of k. We see on Figure 4


Adrien Loiseau

Thibault Marsan

Pablo Navarro

Bruno Watier

Yann Landon

Optimizing racing wheelchair design through coupled biomechanical-mechanical simulation

Keywords: Wheelchair design, Racing, Athlete positioning, Para athletics

The purpose of this study is to optimize the design of racing wheelchairs to improve the performances of the athletes. The design of manual wheelchair allows athletes to express their full potential. Two models have then been created. The first one to compute the optimal position of the shoulder of the athlete relatively to the wheelchair to obtain the maximal wheelchair speed for long distance races. The second one was designed to represent the 100 m race and to optimize the pelvis position of the athlete on the wheelchair to reduce the time to reach 100 m. Our model quantified the maximal speed reached by the wheelchair to 32 km/h and the optimal time to 14.35 s. To obtain these performances, the athlete would be in a lying position, with the vertical position of the pelvis centre close to the vertical position of the shoulder. The second program also returned the optimal speed curve of the wheelchair during the 100 m race. The coaches could then use the optimal acceleration curve found in this study to match the acceleration of the wheelchair of their athlete.

Introduction

Manual wheelchair is a key element in the performance of wheelchair racing athletes. For this reason, the design of such wheelchairs is a particular issue for high performance athletes. Many advances have been made since the first racing wheelchair [START_REF] Bundon | Expert Users' Perceptions of Racing Wheelchair Design and Setup: The Knowns, Unknowns, and Next Steps[END_REF] but improvements remained to be made. The purpose of this study was to optimize the design of racing wheelchairs to improve the performance of the athlete. As a design requirement, the optimal position of the athlete in his wheelchair must be determined for each sport and each athlete, as well as the global mass repartition [START_REF] Mason | The ergonomics of wheelchair configuration for optimal performance in the wheelchair court sports[END_REF]. Indeed, in wheelchair racing the athletes present very different pathologies which necessitate customized settings. The work presented in this paper concerns athletics, especially wheelchair racing. The manual wheelchair designed for racing is specific to this sport. At first, it was just a modified wheelchair [START_REF] Macleish | Design of a composite monocoque frame racing wheelchair[END_REF]. One of the major changes in the manual racing wheelchair was the evolution of a four-wheeled wheelchair to a threewheeled wheelchair [START_REF] Cooper | Wheelchair racing sports science: a review[END_REF], reducing the total weight of the wheelchair. Then, the materials evolved [START_REF] Cooper | Engineering and Technology in Wheelchair Sport[END_REF], reducing the weight once more to make the wheelchair easier to propel and maneuver. It would then be interesting to study the optimal position of the athlete on the wheelchair, that would allow him to express his full potential.

It has already been demonstrated that the position on the wheelchair would influence both the pushing pattern and the kinematics of the subject [START_REF] Mâsse | Biomechanical analysis of wheelchair propulsion for various seating positions[END_REF]. Different methods have been used in the literature. Experimental data were collected using inertia rollers, inertial measurement units, or motion capture [START_REF] Hybois | Approche numérique pour l'optimisation personnalisée des réglages d'un fauteuil roulant manuel[END_REF]. Propulsion models have also been developed [START_REF] Morrow | A 2-D model of wheelchair propulsion[END_REF]. It has been stated that, since most of the motion occurred in the sagittal plane, a 2D model of the wheelchair would be sufficient to study the wheelchairs' propulsion. Most of the developed models in the literature were used to compute optimal parameters of the wheelchair. However, none of the models of the literature have been used to obtain the optimal position of the athlete on racing wheelchair.

Para-athletes compete in different racing length, from 100 m to marathons, and the athletes' kinematics changes depending on the race. For instance, T53 athletes, which is a category of para-athletes with abdominal deficiencies who have little to no control over their trunk and their legs, use their shoulder to propel the wheelchair during the first acceleration but not once they reached their maximal speed [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF]. However, T54 athletes, who also have low control over their legs but can control their trunk, use their shoulders both during acceleration and while maintaining their maximal speed. During the 100 m race, athletes keep accelerating the wheelchair since they reach their maximal speed at the end of the race [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF]. It would then be interesting to develop two models, one with a fixed shoulder, which would be useful for long-distance races, and another one with a moving shoulder, for the short-distance races.

In this context, the desired optimal position of the athlete in his wheelchair aims at: reducing the aerodynamic drag, maximizing the mechanical work exerted on the pushing rim at each pushing cycle, and minimizing the loss of energy especially for long races. Thus, this study focused on obtaining the optimal position of the athlete on the wheelchair for two scenarios, using dynamics models. The first model aimed at maximizing the wheelchair speed during the permanent regime, for long-distance races. The second one focused on minimizing the time to reach 100 m.

Material and methods

Two models have been designed, with different objectives and hypotheses. The main difference between the two models was the motion of the shoulder, which was fixed for the long-distance race model but not for the short distance model. This was due to the fact that athletes had a higher torso amplitude, therefore, a high shoulder motion, during the first accelerations, than when they reached their cruising speed [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF].

Long distance race model

This model aimed at finding the best position of the athlete in its wheelchair to reach a maximal permanent speed.

For that objective, the energy provided by the athlete to the wheelchair was computed, depending on the position of the athlete on the wheelchair and its velocity. The optimal position was defined when this energy was maximized. This long-distance race model (Fig. 1 on the left) was implemented with the following assumptions:

-The athlete's trunk was fixed.

-The shoulder and elbow joints were assimilated as pivots.

-The hand was assimilated as a point on the pushing rim.

-The model only considered the rear wheels and the athlete.

-The speed during the pushing cycle was considered constant. The model varied the position of the athlete's shoulder relatively to the centre of the wheel. For each position, the total energy of the pushing cycle needed for a specific speed of the wheelchair was calculated using an incremental method. A pushing cycle was defined as the elementary cycle of action of the athlete on the pushing rim, between the first contact with the pushing rim until the athlete removed his hands. The rotation was split into 360 points, each point corresponding to one degree (angle a). Between two points on the pushing rim, the point representing the wrist made a circular rotation of one degree.

The rolling resistance was modelled using the model developed by Masson [START_REF] Masson | Amélioration des performances d'athlètes en fauteuil roulant d'athlétisme à partir d'une configuration optimale[END_REF] and the aerodynamic drag was adapted from the Samozino's model [START_REF] Samozino | A simple method for measuring power, force, velocity properties, and mechanical effectiveness in sprint running[END_REF].

For each point on the pushing rim, the position of each joint centre and the angles of the joints were computed. Then, the joint speed was also computed using the information of a second point on the pushing rim close to the first one. Since, we did not have access to data about the forces generated by the athletes at their shoulder and torque during a propulsion, the values have been based on the literature. It was then possible to obtain the shoulder [START_REF] Mayer | Normal values of isokinetic maximum strength, the strength/velocity curve, and the angle at peak torque of an degrees of freedom in the shoulder[END_REF] and elbow torques [START_REF] Provins | Maximum torque exerted about the elbow joint[END_REF][START_REF] Colson | Isokinetic elbow flexion and coactivation following eccentric training[END_REF] based on the joint speeds. The torque value used in our study was the maximal torque exerted. The torques developed at the joints by the muscles created a force on the pushing rim, and the tangential part of this pushing force was considered.

The scalar product of the tangential force and the motion of one degree of the hand on the pushing rim corresponded to the energy supplied by the athlete to the wheelchair between two points. The energy computed for all the points of the push cycle was summed up to find the total energy. If the total energy (W) was greater than the energy required to counteract the resistive forces (En), the algorithm increased the wheelchair speed and restarted the calculation for the same position, but with the new wheelchair speed. Otherwise, the maximum speed was reached.

Sprint race model

The long-distance model used the maximal speed to obtain the best position. However, for a 100 m race, the maximal speed was reached only on the last few meters. Indeed, since the athletes keep accelerating all along the 100 m, the maximal speed reachable by the athlete in the wheelchair might not be the best criterion for performance. Consequently, another model was developed to determine the best position in the wheelchair in order to minimize the time needed to reach 100 m.

When starting, athletes tended to be mobile with their trunk whereas their trunk was fixed when they reached their cruise speed [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF]. It was therefore considered in this second model since the angle between the trunk and the wheelchair influenced the shoulder and elbow angles.

This sprint race model (Fig. 1 on the right) was implemented with the following assumptions:

-The shoulder and elbow joints were assimilated as pivots.

-The speed of the wheelchair was assumed to be constant during one pushing cycle.

-The athlete considered is in the category T53 which meant that he did not have control over his abdominal muscles to maintain his posture, and that his back was bound to the wheelchair. The same methodology than the first model, in terms of computation of the tangential force, was applied to this second one with the difference that the shoulder position was not fixed. Also, the parameter computed to stop the iterating process was the distance reached by the athlete. The model stopped when the distance travelled (Length) reached 100 m. The time needed to reach these 100 m was then computed.

Dataset considered to test the models

To test the different models, a dataset including the entry parameters of the models was defined. First, the anthropometric data such as torso, arm, and forearm length were needed. As an example, to try out the models, test-data were taken from one of the authors. The arm length was set at 250 mm, the forearm length at 300 mm and the torso length at 550 mm. The parameters of the wheelchair needed were only the rear wheel diameter and the pushing rim diameter. The rear wheel diameter was set at 700 mm and the pushing rim diameter at 350 mm, which were based on the observation of an existing racing wheelchair. Then, constraints had to be implemented. First, the degrees of freedom of the joints were set between a minimal and a maximal value: the torso angle was set between 𝛽 "#$ = 20° and 𝛽 ")* = 45°, the shoulder joint angle between 𝛿 "#$ = -30° and 𝛿 ")* = 170°, and the elbow joint angle between 𝜃 "#$ = 20°and 𝜃 ")* = 175°. These values were also used as an example to test the model. These angles (Fig. 2) were used to determine when the hand is in contact with the pushing rim. Then, for each position of the hand on the pushing rim, the position of the other joints was computed as well as their instantaneous speed. The torques generated by the joints were then taken from the literature based on their angle and speed [START_REF] Mayer | Normal values of isokinetic maximum strength, the strength/velocity curve, and the angle at peak torque of an degrees of freedom in the shoulder[END_REF][START_REF] Provins | Maximum torque exerted about the elbow joint[END_REF][START_REF] Colson | Isokinetic elbow flexion and coactivation following eccentric training[END_REF]. The torque considered was assumed to be the maximal torque obtained in those studies. This assumption could be discussed as the subjects of these studies were able-bodied whereas our study concerned para-athletes.

Results

Results of the long-distance race model

Given the dataset considered, the long-distance race model was first tested. The first iteration reached a maximal speed of 8.33 m/s with a position of the shoulder relative to the centre of the wheel of 𝑋 56789:;< = 100𝑚𝑚 and 𝑌 56789:;< = 350𝑚𝑚 using a step of 𝛥 = 20𝑚𝑚 in both axes. Then, with a second optimisation using a step of 𝛥 = 2𝑚𝑚, the optimal position of the shoulder relative to the wheel centre was 𝑋 56789:;< = 80𝑚𝑚, and 𝑌 56789:;< = 332𝑚𝑚. With this position, the maximal velocity reached by the athlete on his wheelchair was 8.88 m/s. This first result seemed to be in accordance with the actual position used by the athletes. These iterations allowed us to map the optimal position of the shoulder of the athlete as a function of the total energy generated by athlete on his wheelchair (Fig. 3). It was observed that, for the athlete with the anthropometry and the wheelchair geometry used in this study, the optimal position of his shoulder would be between 𝑋 56789:;< = 0 to 90mm, and 𝑌 56789:;< = 300to 400mm. The evolution of the total energy as a function of the 𝛼 angle is presented on Fig. 4. It can be observed that there was no evolution of the total energy between 𝛼 = 0° to 70° and between 𝛼 = 260° to 360°. This was because the athlete was not able to hold on the pushing rim for these angles. This incapacity to grab the pushing rim was linked to the anthropometry of the athlete and his position on the wheelchair.

Results of the sprint race model

Using the same dataset, the sprint race model was then tested. The first iteration of the second model used the following position of the pelvis, with respect to the wheel centre: 𝑋 >;9?#5 = 100 mm and 𝑌 >;9?#5 = 50 mm. The time needed to reach 100 m was then of 𝑡 FGG" = 21.21𝑠 and the maximal wheelchair velocity reached was 𝑣 ")* = 5.8 𝑚 𝑠 ⁄ . After the optimisation of the position of the pelvis on the wheelchair, the time returned was 14.35s to complete the race. The maximal velocity of the wheelchair reached was then 9.08 m/s. The optimal position of the hip was 𝑋 >;9?#5 = 386𝑚𝑚 and 𝑌 >;9?#5 = 386mm.

The results of the first model suggested that the shoulder should be between 𝑌 56789:;< = 300 and 400 mm to reach an optimal position. The second model returned about the same vertical position of the pelvis. This means that the athlete should have different settings depending on the race he's participating in.

The wheelchair speed as a function of the time was computed (Fig. 5) and it can be observed that it followed an exponential curve. Also, the speed always increased, and it can be suspected that there would be a plateau when the maximal speed was reached but this maximal speed was not reached during the 100 m race.

Discussion

The values of maximal wheelchair speed obtained with the Long-distance race model was 𝑣 ")* = 8.88 m/s. This value was close to the one found in the study by Sauret et al. [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF] where experimental data have been captured using IMU. Indeed, they found a maximal speed of 8.86 m/s for the first athlete and 8.09 m/s for the second athlete. However, the position of the shoulder used by these athletes was more towards the front of the chair rather than close to the wheel centre, which was suggested by our optimization.

Concerning sprint races, in the literature [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF], two athletes performed a 100 m race and the best of them had a time of 15.38s and reached his maximal speed of 8.86 m/s at the end of the race. Since it was a short race, the speed of the athlete kept increasing, until the 100 m. However, for longer races, a plateau could be seen. This phenomenon was also observed in our results with the sprint race model (Fig. 5) where the speed kept increasing during the first hundred meters.

Therefore, the use of both models can be useful to coaches and athletes, as it may give the optimal positions of the pelvis and the shoulder to obtain either the best acceleration pattern and/or the highest maximal speed during both short and long-distance races. Nevertheless, the comfort of the athlete should also be considered when selecting the optimal position on the wheelchair. Indeed, if the athlete would not be able to maintain this optimal position due do discomfort, then this position would not be optimal for the performance of the athlete.

The dataset used in this study was not taken from a Paralympic athlete, which could explain the differences found in this study. Indeed, the parameters used in the models were based on the measurements made on one of the authors. Differences could exist between the actual measurements of a Paralympic athlete and the model could gain in accuracy with a dataset coming from an actual athlete.

One of the hypotheses made for the first model was that the shoulder was not moving during the propulsion when the athlete reached his maximal speed. This was the case for the T53 athletes since they do not have control over their abdominal muscles but T54 athletes do have control over those muscles. Therefore, they had a non-negligible motion of their shoulders with an amplitude of almost 40° [START_REF] Sauret | Performance analysis of wheelchair racing using IMU sensors[END_REF]. The optimal position of the shoulder found with the first model could also be used for the athletes with control over their abdominal muscles as a mean position of the shoulder during propulsion, and the athletes would be able to generate even more mechanical work during the cycle.

The speed of the wheelchair was deemed constant throughout the pushing cycle; however, as it can be observed on the Fig. 5, the hand was not in contact with the pushing rim during the whole cycle. Neglecting this allowed us to compute the optimal acceleration curve of the wheelchair. For coaches, this might be seen as the optimal mean acceleration curve that the athletes should aim to achieve.

As it can be observed on the Fig. 5, the sprint race model over-estimated the speed of actual athletes. Those differences could be due to the fact that the athletes may not have the optimal position in their wheelchairs, or the wheelchair might be poorly designed. It might also be due to the different hypotheses made in this study. Indeed, the torques computed in this study were based on the literature, and they were computed based on the angles and rotational speeds of the joints. The torque applied was then the maximal torque the joint could supply. Even if the athletes were the best in their field, they were not able to provide maximal torque during the whole cycle. Also, only the in-plane components of the force applied by the hand on the pushing rim was considered. However, the force applied by the hand is three-dimensional [START_REF] Chénier | A high sample rate, wireless instrumented wheel for measuring 3D pushrim kinetics of a racing wheelchair[END_REF]. Therefore, the out-of-plane part of the force was neglected, and this force may be detrimental to the propulsion, so neglecting it may result in an increased wheelchair speed. This could be one of the improvements to be made in the future. Considering the three components of the force on the pushing rim dynamics should result in a diminution of the wheelchair speed, and the simulated data would probably match better the experimental ones. The model would also be more realistic if the wheelchair speed was not supposed constant during the pushing cycle. Indeed, considering the freewheel phase of the pushing cycle would lead to a more realistic propulsion phase. Also, going from a 2D model to a 3D model could help considering more than the tangential component of the force of the hand on the pushing rim [START_REF] Chénier | A high sample rate, wireless instrumented wheel for measuring 3D pushrim kinetics of a racing wheelchair[END_REF]. Another possible evolution of the models would be to consider the shoulder and elbow as actual joints with multiple degrees of freedom [START_REF] Bourgain | Effect of shoulder model complexity in upper-body kinematics analysis of the golf swing[END_REF], which would result in a more complex, but more realistic model of the propulsion.

Finally, after the optimization proposed, a last model was developed to determine the optimal position of the centre of mass of the athlete and his wheelchair set in order to avoid the wheelie effect during the acceleration phase [START_REF] Kirby | The manual wheelchair wheelie: A review of our current understanding of an important motor skill[END_REF].

Conclusion

The models developed in this study allowed to compute the optimal position of both the shoulder and the pelvis of the athlete to generate the maximal mechanical work on the pushing rim, leading to a maximal wheelchair speed. This could help athletes, professionals as beginners, to set up and/or optimize their wheelchair and to be seated at the optimal position that would allow the maximal speed. The second model might also be helpful for coaches, who could set the optimal acceleration curve as an objective to their athletes.

To improve the models, the 3D force applied by the hand on the pushing rim might be considered. Also, considering the freewheel phase of the pushing cycle could provide information about a more realistic pushing cycle.

Annex

The formulas used to compute the kinematic parameters of the model are presented in this annex. This allowed us to find the shoulders and elbow angles.

The first step to calculate the kinematics parameters of the Long-distance model was based on the current position of the shoulder centre (𝑋 56789:;< and 𝑌 56789:;< respectively the horizontal and vertical position): For the sprint model, the first step to compute the kinematic parameters was based on the dimension of the torso (𝐿 i7<57 ), the horizontal and vertical position of the pelvis (respectively 𝑋 >;9?#5 and 𝑌 >;9?#5 ), and 𝛽, the angle between the torso and a horizontal line (figure attached).

Fig. 1 .

 1 Fig. 1. Workflow of the long-distance race model (left) and sprint race model (right)

Fig. 2 .

 2 Fig.2. Definition of the angles used in the model. 𝛼 is the angle between the hand and the shoulder, 𝛽 is the torso angle, 𝛿 the shoulder joint angle, and 𝜃 the elbow joint angle. 𝑋 56789:;< , 𝑌 56789:;< , 𝑋 >;9?#5 , and 𝑌 >;9?#5 are respectively the horizontal and vertical position of the shoulder and the pelvis, with respect to the wheel centre.

Fig. 3 .

 3 Fig. 3. Map of the total energy supplied by the athlete to the wheelchair as a function of the horizontal and vertical position of the shoulder.

Fig. 4 .

 4 Fig. 4. Evolution of the simulated total energy as a function of the angle 𝛼. The hatched area represents the angles for which the athlete was not able to grab the pushing rim.

Fig. 5 .

 5 Fig. 5. Superposition of our simulated data (black) with the sprint race model and the experimental data of the literature [9] (with authorization) for different distances.

  𝑋 56789:;< = 𝐿 i7<57 cos(β) -𝑋 >;9?#5 ; 𝑌 56789:;< = 𝐿 i7<57 sin(𝛽) + 𝑌 >;9?#5 𝑖𝑛𝑐 = arctan n 𝑌 56789:;< 𝑋 56789:;< o ; 𝐷 = p(𝑋 56789:;< ² + 𝑌 56789:;< ²)Then, the next step was the same for both models. The parameters used were: 𝛼, the

Acknowledgments

This work has benefited funds managed by the National Research Agency (ANR) under the "Future Investment 3" plan within the framework of the "Priority Sport Research Program of Very High Performance" bearing the reference ANR--19-STHP-0005.