Arthur Fétiveau
email: arthur.fetiveau@univ-ubs.fr

Permutations based model for business performance

Keywords: Computational statistics, Generalized Mallows models, Kemeny consensus, Kendall distance, Machine learning, Mallows models, Permutations

This article is devoted to the development of a machine learning statistical framework to drive company's objectives. To this end, their sales data are used to target efficiently the issues or opportunities by a ranking. We implement a permutation based model using generalized Mallows models dealing with quantitative values, considering that a ranking is a permutation. The advantage of the generalized version is the possibility to differentiate the cost to move each element in the permutation. In our model, we differentiate the cost of an inversion in the permutation by using the gap value between the two elements. We propose model parameters estimators and we illustrate our estimation procedure on simulated data and on a real application.

Introduction

In this paper, we are interested in the development of a machine learning statistical framework to drive company's objectives. The defined approach helps the decision makers choices using the company's sales data. We have to learn what is important to know from the decision makers point of view to efficiently help them to target the issues or opportunities.

To explain the objectives and for the sake of example, we study a company selling fishing products worldwide. This company sells a large number of different products like fish hooks, fishing nets, lures, fishing rods and more. We may have variations of each of those products in size, weight, material, color or even package, multiplying the total number of unique references to evaluate and analyze by the company. As this company makes profits worldwide, it sells each of those references in a bunch of stores, located in many countries with different climates. For example, this company can sell fishing products in France, Iceland, Mexico and Australia. In every country, we can also make a distinction between each selling region. For instance in France, we can sell in Brittany, New Aquitaine, Normandy, which are coastal regions but also in Great East as you can also fish in lakes or rivers. Depending on your responsibilities in the company, you are interested in a set of references in a geographical area. For example, the sales director of Brittany is interested in products sold in Brittany and is not interested if the company makes profit selling fishing rods in Australia.

For every product in stores, we can have information on a lot of indicators. For example, the revenue of the product, the selling price, the expected revenue leading to information on seasonality and many other information. Each criterion has to be taken into account to target the best interest of the user. For instance, the Normandy sales director is not interested to analyze the lure sales for salmons in the Channel during the winter as we can only fish them in spring and summer. Even with the deletion of the uninterested products, there is still a wide range of potentially interesting products. It would be really expensive, in money or time, to analyze every one of them at once. So, it is important to rank all of the references in a customized importance order. It is mandatory to prioritize the elements with the best potential for our Normandy sales director, but also every other sales director. For example, in September, it is more important for the Normandy sales director to analyze an important gap between the expected and the real revenue in the lures for cuttlefish than a huge gap between the expected and the real revenue in the lures for salmons. Indeed, as it is the beginning of the fishing season of cuttlefish and the end of the fishing season of salmons, a smaller gap may lead to a more important future profit. Once every sales director is provided with the indicator of his own importance order, he can study the elements in the decreasing order, and makes a plan for each one to increase their benefits. If our Normandy sales director has the time after analyzing the cuttlefish's lures sales, he may then study the salmon's lures sales.

Our objective is then to create a statistical and computational framework or-dering every reference to match the user's preferences. This method uses the indicators given for all the products at every geographical level. However, this is not an easy supervised learning as we do not have any truth about the potential of the product. To help us with this, we use an information given by the user by reinforcement learning. First, our Normandy sales director, the user, have a ranking of his references; the ranking begins with the cuttlefish's lures, followed by the salmon's lures, and so on. He then analyzes the first references, the most as he can in the resources allowed to this task. Then, he returns an evaluation of his interest for every analyzed reference. Based on the user returns, we adjust the statistical model to fit the user preferences. This process iterates, month after month, to find the perfect ranking for the user.

As it is impossible for the user to give the perfect ranking of several hundreds of references or even to quantify the difference of two references, we will only ask for a value on a one to five ladder based on his interest. We then need to use a method able to deal with criteria composed of numerical continuous values, like the revenue where a difference of 10 euros is not the same as a difference of 100 euros, as well as ranking only criterion, like the partial ranking of the user preferences. To this end, we develop a permutation based procedure using generalized Mallows models. We consider here that a ranking is a permutation. The advantage of the generalized version it to be able to differentiate the cost to move each element in the permutation. However, we need to be able to differentiate the cost of an inversion in the permutation, for example by using the gap value between the two elements.

As we have seen in our example, the sales director of New Aquitaine where the fishermen fish the mackerel in the Bay of Biscay does not have the same kind of interest in the different criteria as the sales director of Great East. The seasonality of the mackerel is not as important as the salmon's seasonality in the large mountain rivers. Maybe, if you are the world sales director, the large amount of different climates and latitudes will smooth all the seasonality. It is then important to have a method able to be generic, to fit to every possible desire of the user. Each profile must be able to have its own set of parameters. We must also be able to compare each set of parameters between the different user to understand the differences between them.

This paper is organized as followed. In Section 2, we define the Mallows models, its generalization and limitations in our case. Section 3 deals with the parametric estimation of the model parameters. Section 4 introduces our modifications to integrate the cost of an inversion. Section 5 deals with the study of our procedure on simulated data while section 6 is devoted to real data application.

Mallows model and its generalization 2.1 Mallows model

The Mallows model [START_REF] Mallows | Non-null ranking models[END_REF] is based on the distances between the permutations of a set of n items to order. A permutations is a one to one application from {1, ..., n} to itself. The set of all the permutations of n items is denoted by S n . We write π(i) the position associated to the item i in the permutation named π for i ∈ {1, ..., n} and π -1 (i) the item associated to the i th position of the permutation named π for i ∈ {1, ..., n}. The cardinal of S n is equal to n!.

In this framework, we consider the permutation π as the perfect theoretical ranking for the user. The Mallows model is then defined as the probability to select randomly a permutation σ, in the set of all the permutations S n , based on the permutation π with a dispersion parameter θ. We have

P θ (σ) = exp(-θ d(σ, π)) Z(θ) (1)
with θ ≥ 0, d(., .) a right invariant distance between two permutations of S n and

Z(θ) = σ∈Sn exp(-θ d(σ, π)), (2)
the normalization term which is not influenced by π, thanks to the right invariant property. A distance d(., .) is right invariant [START_REF] Diaconis | Group representations in probability and statistics[END_REF] if d(π, σ) = d(πτ, στ ) for every π and σ in S n . The normalization term Z(θ) only varies with θ. This normalization term is then easy to compute in a time of order O(n). Moreover, in the Mallows we tries as much as possible to avoid computing this value. As we can see in Section 3, we do not need it to estimate our parameters. When θ = 0, then each permutation σ of S n has the same probability. The larger θ, the larger the density of σ around the modal permutation π. The distance is right invariant, meaning that a change of the labels in the permutations does not change the distance. To continue with our fishing example, if our permutation σ is first the mackerels fishing nets, then the salmon's lures and finally the fishing rods, and our permutation π is first the fishing rods, then the mackerels fishing nets and finally the salmon's lures, then renaming the mackerels fishing nets to 1, the salmon's lures to 2 and the fishing rods to 3 will gave the same distance as renaming the fishing rods to 1, the mackerels fishing nets to 2 and the salmon's lures to 3.

In the Figure 1, every node represents a permutation and every link represents a distance of 1 between 2 permutations. Here, we take the black node as π. Then, the probability of the black permutation, π which is at a distance of 0 is exp(0) Z(θ) . The probabilities of the red permutations are exp(-θ)

Z(θ)
and the probabilities of the orange permutations are exp(-2θ) Z(θ) .

Generalized Mallows model

The generalized Mallows model [START_REF] Michael | Distance based ranking models[END_REF] is an extension of the Mallows model based on a decomposition of a right invariant distance. If we have

d(σ, π) = p k=1 S k (σπ -1 , e), (3)
with e the identity permutation obtained by ππ -1 , then, the parameter θ can be replaced by θ k for k = 1, . . . , p. We obtain

P θ (σ) = 1 Z(θ) exp p k=1 -θ k S k (σπ -1 , e) (4)
where

Z(θ) = π∈Sn exp p k=1 -θ k S k (σπ -1 , e) (5)
corresponds to the normalization term not influenced by π. The Mallows model is a particular case of the generalized Mallows model with θ = θ k , k = 1, . . . , p when all θ k are the same.

We decide here to use the Kendall distance [START_REF] Kendall | A new measure of rank correlation[END_REF]. The Kendall distance represents the minimal number of adjacent transpositions to transform a permutation σ into a permutation π, both from the set of permutations S n :

d k (π, σ) = n-1 i=1 j>i 1 ((π(i) > π(j) ∩ σ(i) < σ(j)) ∪ (π(i) < π(j) ∩ σ(i) > σ(j))) (6)
where 1(P) takes the value 1 or 0 depending on whether the condition P is satisfied or not. This distance is right invariant.

Using our fishing example, the Kendall distance between our permutation σ, the mackerels fishing nets (A), the salmon's lures (B) and the fishing rods (C), and our permutation π, the fishing rods (C), the mackerels fishing nets (A) and the salmon's lures (B), is 2 as the mackerels fishing nets (A) is inverted with the fishing rods (C) and the salmon's lures (B) is inverted with the fishing rods (C).

With the Kendall discrepancy function, we can determine for k = 1, . . . , n -1, S k (σπ -1 , e). The minimum number of adjacent inversion to assign the first item of e in the first position of σπ -1 is S 1 . We then remove this item from both permutations to compute S 2 which will be the minimum number of adjacent inversion to assign the remaining first item in the remaining first position. We repeat this process until there is not any remaining item. We have S k (σπ -1 , e) = l>k 1((σπ -1 ) -1 (k) > (σπ -1 ) -1 (l)). [START_REF] Mallows | Non-null ranking models[END_REF] With this distance, we assume that each item of π has its own parameter θ k , depending on its position in π, to move and does not depend on the elements inverted with this item. However, in our problem, for some criteria, it may be important to make a difference between each inversion, and we want to determine a new parameter associated to each element inversion.

Estimation of the model parameters

We assume that the permutations σ j are known and represent the ranking of the J criteria, each of them picked according to a Mallows model with the same modal permutation π but with different dispersion parameters θ j for j ∈ {1, . . . , J}. We want to estimate π and every θ j . Each θ j is invariant through time, meaning that each month the criteria keep the same importance for the user. The permutation π changes each month however, due to the variation of important products according to the criteria. So, each month, we have J permutations picked from the Mallows models, with still the same θ j but a different π. In this Section, the Kendall distance is considered.

Estimations of π

Knowing, for j ∈ {1, . . . , J}, all the σ j and their associated dispersion parameter θ j for a month, we can estimate the permutation π associated to this month. We use the maximum likelihood estimator method. The likelihood function is given using (1) by

L(π|σ, θ) = J j=1 exp(-θ j d(σ j , π)) Z(θ j ) . ( 8
)
As Z(θ j ) does not depend on π and the exponential function is monotonous, the maximum likelihood estimator π is given by

π = argmax π∈Sn J j=1 -θ j d(σ j , π) (9)
where S n is the set of every possible n items permutations. We use a method based on Blin idea [START_REF] Blin | Median of an odd number of permutations[END_REF], to find the median of J permutations. First, we need to initialize π. We use for this initialization the permutation σ j associated with the biggest θ j . Then, we compute the matrix of discrepancy between π and σ j , that we name G π,σ j . The matrix G π,σ j is a n × n symmetric matrix with 0 on the diagonal since an item cannot be inverted with itself. Each item of this discrepancy matrix between π and σ j , G π,σ j ,a,b equal 1 if items a and b are in a different ordered between π and σ j and 0 otherwise. We create each discrepancy matrix in the same order as the initialization of π in order to sum them. Finally, this matrix is multiplied by θ j .

Each item of our discrepancy matrix G = J j=1 G π,σ j takes its value x between 0 and J j=1 θ j . Furthermore, if we invert two elements of our permutation π, then the new value of the element in our discrepancy function x new , is

x new = J j=1 θ j -x. ( 10
)
While an inversion of two items can diminish the sum of elements of G, we move items in π. When there is no more possible moves, we stop and return π, the estimation of π.

With this method, we face a few issues. First, we can have a lot of permutations maximizing the likelihood function. But the main issue appears if our parameters θ j are too different from each other. If there is j such as J i=1,i̸ =j θ i -θ j < 0 then, the maximum likelihood is given for the permutation π = σ j . So, we cannot have any compromise. To solve this issue, we can introduce, to estimate π the following minimization problem

π = argmin π∈Sn J j=1 (θ j d(σ j , π)) 2 . ( 11
)
Using (11), we need to modify our algorithm. To this end, we use the matrix G π,σ j and modify π while we can reduce J j=1 (θ j d(σ j , π)) 2 where d(σ j , π) is the sum of each element of G π,σ j , corresponding to the distance between π and σ j .

We can also use the Borda [START_REF] De Borda | Mémoire sur les élections au scrutin: Histoire de l'académie royale des sciences[END_REF] method to find the permutation π. In this method, we compute a weighted sum of each element position and sort the results in the ascending order. For our example with the same weight to each permutation, the best compromise between, the mackerels fishing nets, the salmon's lures and the fishing rods permutation and the fishing rods, the mackerels fishing nets and the salmon's lures permutation, is the mackerels fishing nets (1+2), the fishing rods (3+1) and finally the salmon's lures (2+3) permutation. This method is significantly faster than the previous ones. A lot of different algorithms were presented and tested in [START_REF] Ali | Experiments with kemeny ranking: What works when?[END_REF].

Algorithm 1 Estimate π

Require: σ 1 , . . . , σ J ∈ S J n and θ 1 , . . . , θ J ∈ R J π ← σ j where j = argmax j∈1,...,J θ j for {j ∈ 1, . . . , J} do for {i ∈ 1, . . . , n} do for {k ∈ i, . . . , n} do if

(π(i) > π(k) ∩ σ j (i) < σ j (k)) ∪ (π(i) < π(k) ∩ σ j (i) > σ j (k)) then G π,σ j ,i,k ← θ j else G π,σ j ,i,k ← 0 end if end for end for end for G ← j∈J G π,σ j G new is an actualization of G with one element moved in π,
G new ≤ G We stop the algorithm when there is no more G new ≤ G

Estimation of θ

The estimator of θ is discussed in Fligner and Verducci [START_REF] Michael | Distance based ranking models[END_REF], among others. To estimate the θ parameter in the Mallows model, we use the moment generative function of the random variable D(π, σ), the distance between π and σ, in the uniform case, when θ = 0.

Theorem 1. Let P θ (σ) defined by (1) be the distribution of σ, with π the modal permutation parameter and θ the dispersion parameter, and let D(π, σ) the random variable of the right invariant distance d(π, σ) between π and σ. Then the expected mean of the random variable D(π, σ) for t = -θ is given by

E θ (D) = ∂ ∂t log(M D,0 (t)), (12)
and the theoretical variance of the random variable D(π, σ) for t = -θ is

V ar θ (D) = ∂ 2 ∂t 2 log(M D,0 (t)). ( 13
)
Proof of Theorem 1. In the sequel, we use D instead of D(π, σ). The moment generative function of D for θ = 0 is given by

M D,0 (t) = σ∈Sn exp(t d(σ, π))P 0 (σ), (14)
Where P θ (σ) is defined by (1) and corresponds to the probability to select σ. As given in Section 2.1, every permutation has the same probability when θ = 0. We obtain

M D,0 (t) = σ∈Sn exp(t d(σ, π)) 1 n! . ( 15
)
The normalization term Z(θ) given in (2) can be written using the moment generative function of D for θ = 0 and t = -θ:

Z(θ) = n! M D,0 (-θ). ( 16
)
Finally, we express the moment generative function of D in θ by

M D,θ (t) = σ∈Sn exp(t d(σ, π))P θ (σ). ( 17
)
The probability of σ knowing θ is given by (1) and using (16) on

M D,θ (t) = σ∈Sn exp(t d(σ, π)) exp(-θ d(σ, π)) Z(θ) = σ∈Sn exp((t -θ) d(σ, π)) Z(θ) , (18)
we obtain

M D,θ (t) = σ∈Sn exp((t -θ) d(σ, π)) n! M D,0 (-θ) = 1 M D,0 (-θ) σ∈Sn exp((t -θ) d(σ, π)) n! (19)
where we recognize the expression (15) with t = t -θ. Finally, we have

M D,θ (t) = M D,0 (t -θ) M D,0 (-θ) . ( 20
)
The moment generative function properties allow us to calculate the mathematical expectation and the theoretical variance of D by a simple derivative. The mathematical expectation knowing θ of D(π, σ) is given for t = 0 by

E θ (D) = ∂ ∂t M D,θ (t). ( 21
)
We know with (20) the link between the moment generative function with θ ̸ = 0 and the moment generative function where θ = 0, the case where every permutation σ has the same probability. We can then replace M D,θ (t) by its equivalent with

M D,0 (t) in t = -θ E θ (D) = 1 M D,0 (-θ) ∂ ∂t M D,0 (t). (22)
We recognize in the mathematical expectation E θ (D) that the numerator is the derivative of the denominator and thus for t = -θ we obtain (12)

E θ (D) = ∂ ∂t log(M D,0 (t)).
The variance is given in t = 0 by

V ar θ (D) = ∂ 2 ∂t 2 M D,θ (t) - ∂ ∂t M D,θ (t) 2 . (23)
Using (20), we have in t = -θ

V ar θ (D) = ∂ 2 ∂t 2 M D,0 (t) M D,0 (-θ) -    ∂ ∂t M D,0 (t) M D,0 (-θ)    2 = M D,0 (-θ) ∂ 2 ∂t 2 M D,0 (t) - ∂ ∂t M D,0 (t) 2 M D,0 (-θ) 2 .
(24) We recognize here the second derivative of the logarithm. Then for t = -θ, we obtain (13)

V ar θ (D) = ∂ 2 ∂t 2 log(M D,0 (t)).
However, even with those simplifications with θ = 0, the calculations of the mathematical expectation and the theoretical variance are still time consuming due to the sum of n! terms. As previously, we use the Kendall distance. This distance helps us to simplify the calculations of the mathematical expectation and the theoretical variance of D(π, σ).

Theorem 2. Under the same assumptions of Theorem 1 with d(π, σ) the Kendall distance. Then the expected mean of the random variable D(π, σ) is given by

E θ (D) = n exp(-θ) 1 -exp(-θ) - n j=1 j exp(-jθ) 1 -exp(-jθ) , ( 25
)
and the theoretical variance of the random variable D(π, σ) is

V ar θ (D) = n exp(-θ) (1 -exp(-θ)) 2 - n j=1 j 2 exp(-jθ) (1 -exp(-jθ)) 2 . ( 26
)
Proof of Theorem 2. If we can rewrite the random variable D(σ, π) as

p j=1 S j (σπ -1 , e), (27)
where each S j (σπ -1 , e) is independent for θ = 0, then we have the following property between M S j ,0 (t), the moment generative functions of S j (σπ -1 , e), and M D,0 (t), the moment generative function of D(σ, π):

M D,0 (t) = j M S j ,0 (t). ( 28
)
The Kendall distance (6), as shown in Section 2.2, can be rewritten by d(σ, π) = n-1 j=1 S j (σπ -1 , e) where each S j (σπ -1 , e) is independent from the others. As the Kendall distance is right invariant, we can use π = e where e is the identity permutation. The product of the moment generative functions S j (σ, e) is given by

n-1 j=1 M S j ,0 (t) = n-1 j=1 1 n! σ∈Sn exp(t S j (σ, e)). (29)
No matter the permutation σ ∈ S n , S j (σ, e), the minimum number of inversion to assign the j th element in the first remaining position of σ, takes its values in {0, . . . , n -j}. As we are in the case where each permutation has the same probability, we need to count all the possibilities for each value of S j (σ, e). Each value of S j (σ, e) is obtained by the same number of permutations n! n-j+1 . Then, we can write

n-1 j=1 M S j ,0 (t) = n-1 j=1 1 n -j + 1 n-j i=0 exp(t i).
(30)

As we can calculate easily the sum of the geometric sequence n-j i=0 exp(t i) and using (30), the derivative of E θ (D) given in (12

) in t = -θ is ∂ ∂t log(M D,0 (t)) = ∂ ∂t log n-1 j=1 1 n -j + 1 1 -exp((n -j + 1) t) 1 -exp(t) . ( 31
)
Calculating our derivative, we have with t = -θ

∂ ∂t log(M D,0 (t)) = (n -1) exp(t) 1 -exp(t) - n-1 j=1 (n -j + 1) exp((n -j + 1) t) 1 -exp((n -j + 1) t) . ( 32
)
We can rewrite (32) as

∂ ∂t log(M D,0 (t)) = (n -1) exp(t) 1 -exp(t) - n j=2 j exp(jt) 1 -exp(jt) . (33)
Finally, for j = 1, both parts of the equation are the same and we obtain (25)

E θ (D) = n exp(-θ) 1 -exp(-θ) - n j=1 j exp(-jθ) 1 -exp(-jθ)
However, this last equation is not defined for θ = 0, which can induce some difficulties during the optimization process. To stop those issues, we can keep the sum of exponential. In this case, we have

E θ (D) = n-1 j=1 n-j i=0 i exp(-iθ) n-j i=0 exp(-iθ) . ( 34
)
To calculate the variance of D(σ, π) in θ, we determine the derivative of (25

) for t = -θ ∂ ∂t E -t (D) = ∂ 2 ∂t 2 n exp(t) 1 -exp(t) - n j=1 j exp(jt) 1 -exp(jt) . ( 35
)
So with

∂ 2 ∂t 2 n exp(t) 1 -exp(t) = n exp(t)(1 -exp(t)) -exp(t)(-exp(t)) (1 -exp(t)) 2 = n exp(t) (1 -exp(t)) 2 (36) and ∂ 2 ∂t 2 j exp(jt) 1 -exp(jt) = j 2 exp(jt)(1 -exp(jt)) -j exp(jt)(-j exp(jt)) (1 -exp(jt)) 2 = j 2 exp(jt) (1 -exp(jt)) 2 ,
(37) we obtain (26)

V ar θ (D) = n exp(-θ) (1 -exp(-θ)) 2 - n j=1 j 2 exp(-jθ) (1 -exp(-jθ)) 2 .
For a number of months M , we have all the σ m , for m ∈ {1, . . . , M }, each one of them selected according to a Mallows model with a known modal permutation π m , not necessarily the same each month, and an unknown dispersion parameter θ, always the same, that we want to find. We then estimate θ by the maximum likelihood method. We have the following log-likelihood function for θ

log(L(σ|θ, π)) = M m=1 log(P θ (σ m )) = M m=1 log exp(-θ d(σ m , π)) Z(θ) (38)
and so

log(L(σ|θ, π)) = M m=1 (-θ d(σ m , π) -log(Z(θ))) . (39)
Using (16), we obtain

log(L(σ|θ, π)) = M m=1 -θ d(σ m , π) - n i=1 log(i) -log(M D,0 (-θ)) . ( 40
)
To maximize the likelihood over the parameter θ, we differentiate the loglikelihood with respect to θ

∂ ∂θ log(L(σ|θ, π)) = M m=1 d(σ m , π) - ∂ ∂θ log(M D,0 (-θ)) . (41)
Then, the maximum likelihood estimator θ of θ is solution of

∂ ∂θ log(M D,0 (-θ)) = 1 M M m=1 d(σ m , π). ( 42
)
Using the Kendall distance, with (12) and (25), we have the maximum likelihood estimator θ of θ the solution of the following equation :

n exp(-θ) 1 -exp(-θ) - n j=1 j exp(-j θ) 1 -exp(-j θ) = 1 M M m=1 d(σ m , π) (43)
The Fisher information of the θ parameter is obtained differentiating (42) with respect to θ

I(θ) = E ∂ 2 ∂θ 2 log(L(σ|θ, π)) = E ∂ 2 ∂θ 2 M D,0 (-θ) . ( 44
)
With the Kendall distance, we use ( 13) and ( 26), and obtain

I(θ) = n exp(-θ) (1 -exp(-θ)) 2 - n j=1 j 2 exp(-jθ) (1 -exp(-jθ)) 2 . ( 45
)
This second derivative does not depend on the observations σ m for m ∈ {1, . . . , M }. The asymptotic behavior of the maximum likelihood estimator θ of θ is given by the central limit theorem. Theorem 3. Under regularity conditions, we have when M tend to infinity :

√ M θ -θ D -→ N 0, I(θ) -1 (46)
where I(θ) is the Fisher information of θ given in (45).

As said in the introduction however, the real permutation π is really difficult to obtain. Our best option can give us a partial permutation on a smaller set of element. The formulas (25) and (26) only stands in the case π is a permutation, a ranking without any tie. In a case of partial rankings, one needs to adapt its distance to take the ties into account.

New method

When we use a criterion based on a numerical continuous value, we need to use the difference between two values to determine the cost of an inversion so a transformation of the difference is applied. For example, if our Normandy sales director has a reference of lures for the trout with a difference between the expected and the real revenue of 1 euro and a reference for the salmon's lures with a difference between the expected and the real revenue of 101 euros, then our Normandy sales director will clearly be more interested to analyze the salmon's lures reference than the trout's lures reference. However, if the difference between the expected revenue and the real revenue for the trout's lures is 1000 euros and the difference for the salmon's lures reference is 1100 euros, the Normandy sales director will only be slightly more interested to analyze the salmon's lures reference. So in the first case, the cost, for this criterion, of the inversion between the two references must be higher than the cost of the inversion of the two references in the second case.

As it is really difficult to have labeled data, non parametric statistical method is difficult to use. Moreover, all the information we have are on the final result : "Is this reference interesting or not?" based on the user's opinion. This function's parameters and the Mallows model dispersion parameters θ need to be computed at the same moment. As a consequence, a logit transformation is used and we assign the weight of an inversion at the value of the difference between two references after this transformation. With this transformation, the Mallows model corresponds to the special case where each weight is the same.

We still use the same formalism as before but we do not assume that σ j for j ∈ 1, . . . , J, with J the number of criteria used, is selected randomly based on the the modal permutation π, the theoretical best ranking for the user. We determine the probability that π is the best ranking for the user. We have the values of each σ j , the ranking obtained for the criterion j and its associated vector x j ∈ R n where n is the number of items to order. Then, we give the probability of a permutation π based on a modal permutation σ j and the matrix of inversions costs θ j associated to the criterion j with

P σ j ,θ j (π) = exp(-d p (σ, π)) Z(σ j , θ j ) (47)
where d p (., .) is a weighted distance between two permutations of S n and

Z(σ j , θ j ) = τ ∈Sn exp(-d p (σ j , τ )), (48)
the normalization term influenced by σ j due to a change of values for each inversion and θ j but does not depend on π.

We use as d p (., .) a variation of the Kendall distance given by

d pk (π, σ) = n-1 i=1 k>i θ jik 1 ((π(i) > π(k) ∩ σ(i) < σ(k)) ∪ (π(i) < π(k) ∩ σ(i) > σ(k))) (49)
where θ jik is the difference of the transformations

θ jik = f (x j (i)) -f (x j (k)). ( 50
)
The value of the item i for the criterion j is given by x j (i). Here, we use the logit function for the transformation

f (x j (i)) = α 1 + exp(r(a -x j (i))) + β, (51)
with β the minimum value after the transformation, α + β the maximum value after the transformation, r the velocity parameter to go from the minimum to the maximum and a the symmetrical point. These 4 parameters α, β, r and a fully determine all the θ jik . Learning them give the matrix θ j .

Estimation of π

To choose π, assuming that all the criteria are independent. We maximize the function

Ψ(σ, θ j , π) = J j=1 P σ j ,θ j (π), (52)
where P σ j ,θ j (π) is defined in (47). Then, taking the logarithm of Ψ(σ, θ j , π), we have log(Ψ(σ, θ j , π)) =

As the second part of the equation (53) does not depend on π, we can ignore it and maximize log(Ψ(σ, θ j , π))

= J j=1 -d p (π, σ j ). ( 54
)
To maximize this value, we use the method described in Section 3.1 to maximize (9). In this variation, we can reinitialize G π,σ j ,ik with

G π,σ j ,ik = θ jik 1 ((π(i) > π(k) ∩ σ(i) < σ(k)) ∪ (π(i) < π(k) ∩ σ(i) > σ(k))) .
(55)

Estimation of θ

We assume independence between all the criteria. With this independence, maximizing the probability P σ j ,θ j (π), for any j ∈ 1, . . . , J, maximize the function Ψ(σ, θ j , π). Figure 2 gives a graphical representation of what we want to do. In this representation, we want to minimize the cost of the red links and maximize the cost of the black links. Thus, we can estimate the parameter θ j for any criterion j independently of the others. We have, for a number of months M , π m , the wanted permutation considering the criteria. This permutation is the best for the user for the month m. We also know the real permutations σ jm , associated to their real values x jm for each month m ∈ {1, . . . , M }. Also, θ j is the same for each month. As θ j is defined with f , we need to obtain the parameters α, β, r and a. To determine those values, we use the maximum likelihood estimators. The likelihood function is given by :

L(x j , σ j , π; a, α, β, r) = M m=1 exp(-d pk (σ jm , π m )) Z(a, α, β, r, x j , σ jm ) , (56)
where d pk (., .) is defined in (49). Thus, the log likelihood function is given by log(L(x j , σ j , π; a, α, β, r))

= M m=1 -d pk (σ jm , π m ) - M m=1
log(Z(a, α, β, r, x j , σ jm )).

(57) In this case, not as before, the time complexity of Z(a, α, β, r, x j , σ jm ) cannot be reduced in O(n). It is really expensive to compute as the time complexity is in O(n!). For this reason, it is necessary to approximate Z(a, α, β, r, x j , σ jm ). We know how to determine this function with the generalized Mallows model, when all the inversions of the first element with its followings are the same cost, the second with its followings are the same cost as well and so on. To approximate, we can bound the value of Z(a, α, β, r, x j , σ jm ) by using the minimum and maximum value of the transposition cost of an element with its followings. A central value of this bounding, as the mean, is often a satisfying approximation of Z(a, α, β, r, x j , σ jm ).

In the case of the Kendall distance, we have the possibility to separate the distance in a sum of independent distances d k (π, σ) = n-1 i=1 S k (π, σ) as described in [START_REF] Mallows | Non-null ranking models[END_REF].

Using this decomposition, we have for an approximation of Z(θ j , σ jm ) using a unique value of the inversion cost of an element with its followings where θ k is the cost to invert the k th element with its followings. As the S k (τ, σ jm ) are independent, we can, for every permutation of S n , combine every possible number of inversion of the first element, with every possible number of inversion of the second element, with every possible number of inversion of every other element. We obtain

Z approx (θ j , σ jm ) = τ ∈Sn exp n-1 k=1 -θ k S k (τ, σ jm ) , (58)
Z approx (θ j , σ jm ) = n-1 k=1 n-k i=0 exp(-i θ k ). ( 59
)
Finally, the logarithm of the approximation of the normalization term is given by log

(Z approx (θ j , σ jm )) = n-1 k=1 log(1 -exp(-(n -k) θ k )) -log(1 -exp(-θ k )). (60)

Simulation studies

The estimators of π in the Mallows model are studied through simulation. The estimators considered are described in Section 3.1. The real parameter θ j for each permutation σ j , for j ∈ 1, . . . , J is fixed. Each σ j is selected following the Mallows model with the parameters π and θ j . We evaluate π, the estimator of π, using two different criteria: the distance between π and π must be the lowest as possible and the estimator π must be a compromise of every σ j , for j ∈ 1, . . . , J. We consider different combinations of θ j to describe the behavior of the estimators. The real value of π is not important to know as it doesn't affect the results due to the right invariant property of the distance used. In this simulation, the Kendall distance ( 6) is used.

For the results, we denote the estimator of π by MLE the Maximum Likelihood Estimator method given in (9), by QUAD the QUADratic sum of the distances method given in (11) and by BORDA the BORDA method [START_REF] De Borda | Mémoire sur les élections au scrutin: Histoire de l'académie royale des sciences[END_REF].

The simulation results given in Table 1 compare the three methods for 50 σ j , permutations with n = 30 items, each one associated to the same θ over 3000 realizations. The Kendall distance distance between π and π is considered and the means and standard deviations results on the distance are given. This distance is between 0, the minimum value, and n(n-1)

2

, the maximum value. The best results are for a distance of 0 meaning that the two permutations on n items are the same, and a distance of n(n-1) 2 as π is the unique exact inverse permutation of π. The worst distance value is n(n-1) 4 corresponding to an uninformative random noise. For n = 30, the maximum value for the Kendall distance is 435 and the worst value is 217.5.

In Table 1, we observe that the BORDA method gives the best estimation results for π considering the distance to π criterion for the lowest θ parameters. The lowest θ parameters values, the less information contains the criteria associated to the permutations σ j . The BORDA method gives in our simulations, the best results in the sense of the minimum variance comparing to the MLE and QUAD methods for the lowest parameters θ. The maximum likelihood estimator method is better than the BORDA method when θ ≥ 0.1 and the difference with the BORDA method for the distance to π is more important when θ increases. For the standard deviation, the maximum likelihood estimator has a smaller value than the BORDA method when θ ≥ 0.25. Finally, the QUAD estimator is only slightly worst, less than 1 adjacent inversion, than the MLE one but the gap between the two never really change too much.

In Table 2, we consider the same simulation data but we study a compromise of every σ j , for j ∈ 1, . . . , J. To evaluate this information, the Kendall distance between each σ j and the estimator π is considered. The standard deviations of those distances reflect the harmony between the use of all σ j in the decision of a compromise π. The lower this value, the best our estimators for compromises.

In Table 2, we observe that the QUAD method is clearly the best for this simulation. The QUAD method is better than the maximum likelihood estimator no matter the value of θ and better than the BORDA method until θ = 0. the average values are close in values. However, as seen previously in Table 1 for θ ≥ 0.5, the QUAD method give better results than the BORDA method with a lower distance between π and π.

For this criterion however, the results of the maximum likelihood estimator, even though they are the worst, are the most similar to the results obtained with the real value of π instead of any of its estimators (results not shown).

We show in Figure 3 the proximity to π criterion, when the number of criteria used is reduce to 5, still keeping the same value for all θ. We do not observe a significant difference in the distances with the 50 criteria simulation. The BORDA method is better than the maximum likelihood estimator method when θ ≤ 0.35 and the differences between every estimator are lower than 3 adjacent inversions. It is really interesting however to study the differences on the compromise criterion.

In Table 3, we observe the value of the compromise criterion for the different estimators. We observe that the maximum likelihood estimator is clearly outplayed by the others and that the QUAD estimator is clearly better than the others. When we reduce the number of permutations σ j , for j ∈ {1, . . . , J}, taken into account, it is obvious that the maximum likelihood estimator cannot match the QUAD estimator performance on this criterion as we approach the tyranny described in the end of Section 3.1.

In Table 4 and Figure 4, the mean and the standard deviations of the Kendall distance between π and π are given for different values of J, the number of cri- 2: Mean ± standard deviation of the standard deviation of the Kendall distances between σ and π for 50 σ j , permutations of n = 30 items, with the same θ and the same π over 3000 realizations. 3: Mean ± standard deviation of the standard deviation of the Kendall distances between σ and π for 5 σ j , permutations of n = 30 items, with the same θ and the same π over 3000 realizations.

teria used. All the criteria have the same importance. The value of θ and n are respectively 0.3 and 30. The BORDA method is better than the two others in the simulation, for less than 7 permutations to use. We observe that the maximum likelihood estimator is the worst when we only have two permutations. In this particular case, we are in the tyranny due to the limitations of the method. However, the maximum likelihood estimator become the best as the number of permutations used increases. Concerning the QUAD estimator, the results are only slightly worst than the maximum likelihood estimator on the proximity to π criterion. The Table 5 represents the compromise criterion for the same set of parameters than the Table 4. The tyranny of the maximum likelihood estimator when we only have two permutations to take a decision, is also shown by this criterion. Here, we observe also the perfect compromise with the QUAD estimator, which does not necessarily lead to the best result as it is 12 adjacent inversions away in average from the BORDA method for the proximity to π criterion. We also notice that all the methods tend to the same level of compromise when the number of permutations used to find π is increased. For less than 10 permutations however, the QUAD estimator is better, especially when J is even. When J is odd, the compromise is not as good.

We are now interested to evaluate the performances of the three methods with different associations of parameters to find witch one is the more robust to every possibilities. We determine the effect of adding low informative criteria into the decision. We are still comparing on two different objectives. The distance between π and π does not need any change to adapt for this kind of tests but the compromise criterion cannot be taken the same way as before. Obviously, the distance between π and σ strongly depends on θ. As θ increases, we have lower distances between π and σ. Then, we adapt the criterion using the expected distance between π and σ knowing θ given for the Kendall distance with (25). We then have, for the compromise criterion, to study

C 2 = V ar(d(π, σ) -E θ (σ)). (61)
We want to study what happen when we use a different mix of very informative criteria, with θ = 0.6, and low informative criteria with θ = 0.05. In Table 6, we study the Kendall distance between π and π for those cases, with n = 30. We notice that for two parameters θ = 0.6, the set of parameters P 3 , the BORDA method is clearly the best and the maximum likelihood estimator is the worst due to the tyranny problem. We observe particularly for the set of parameters P 5 that adding a lot of low informative criteria strengthen the results for the maximum likelihood estimator. However, the BORDA and QUAD method decrease the quality of the estimation while adding low informative criteria. The set of parameters P 4 shows that adding only a small amount of low informative criteria has the same effect as adding a lot of low informative criteria, like in the set of parameters P 5 , for the BORDA method and the maximum likelihood estimator. This set of parameters P 4 also shows better results for the QUAD estimator than only very informative criteria in the set of parameters P 3 . A mix with a lot of low informative criteria can have better results for the maximum likelihood estimator than the others methods but using only very informative criteria with the BORDA method can show a better precision in the estimation of π.

We need to be careful when we have a lot of very informative criteria. With the set of parameters P 6 and P 7 , we witness that with 25 criteria having a parameter θ = 0.6, adding 25 criteria with a parameter θ = 0.05 does not increase the results, even with the maximum likelihood estimator. However, for the maximum likelihood estimator, the loss in precision is not really a problem as it is only 0.02 more adjacent inversions in average to reach π and this estimator is better than the BORDA and QUAD methods on this criterion anyway.

We use for the Table 7 the same sets of parameters than in the Table 6 but for the compromise criterion using (61). With this criterion, we cannot compare the different lines as different values of θ will also change the variance around the expected value of the distance, as it is well shown in Figure 5. This Figure, for the set of parameter P 4 shows also the main difference between the QUAD method and the two others. The QUAD method tends to favor low informative criteria more than high informative criteria in comparison with the BORDA method and obviously the maximum likelihood estimator. When coming to the compromise criterion, the QUAD method gives the best results. However, in the sets of parameters tested here and excepting the special case P 3 , the difference between the three methods is not statistically significant. The behavior of our different estimators is studied in two cases: all different values of θ and with a big value of θ equalized by a sum of small values of θ. Those two cases are shown in Table 8, for the criterion of the distance between π and π, and are denoted respectively by P 8 and P 9 . We can see that for this criterion and the set of parameters P 8 , the maximum likelihood estimator, as expected, is the best. The surprising part of this Table is with the set of parameters P 9 . We would have expected the maximum likelihood estimator to be the worst as with the set of parameters P 3 due to the tyranny. The results for the maximum likelihood estimator are the same as P 3 . However, it is by far the best, meaning that the confusion induced by the low informative criteria with θ = 0.1 confuses the two other methods and more significantly the QUAD method. The Figures 6 and7 represent, for respectively the sets of parameters P 8 and P 9 , the distances between π and σ. We can see in Figure 6 that for the set of parameters P 8 , the maximum likelihood estimator over-valuate the criterion with the highest value and under-valuate the others. This is the exact opposite as the QUAD method which tends to under-valuate the criterion with the highest value and over-valuate the others. Here, only the BORDA method seems to respect the importance of every criterion. For the Figure 7, it is essentially the same problems than Figure 6 at a higher level. Here, the maximum likelihood estimator uses the ranking of the criterion associated to the θ with a value of 0.6 as π without any compromise. This is a tyranny.

Application on real data

Now, we want to compare the different models on real data. The data set is composed of items to order for 4 consecutive months. Each month has a different number of items. The first month has 325 items, the second 328 items, the third 331 items and the last 340 items. Those items have to be ordered with respect to 7 different criteria. The 7 criteria are not developed here as this is not important to understand the results. Due to the complexity to get the real value of π, in this experiment, we first use a Mallows model with values of θ given by the user as a first intuition. Considering the results, the user evaluates the 40 best items obtained by this Mallows model considering the maximum likelihood estimator (9) used to find π. The user gives to the 40 items a notation between 0 and 9, where 0 correspond to the best notation, a very interesting product to study and 9 correspond to the worst notation, a useless product to study. Based on these notations, we make a partial ranking of π for those 40 items. So, we learn the parameters θ j for the Mallows Model (1), denoted as MM, and a, α, β, r for the model given in (47) which uses the real values of each criteria, considering different approximations of log(Z(a, α, β, r, x j , σ jm )) in the likelihood (57).

As an approximation of log(Z(a, α, β, r, x j , σ jm )), we use 3 possibilities. The MEDian, denoted as MED, the MEAN denoted as MEAN, and MINMAX, the mean of the MINimum and MAXimum values of the logarithm. For the 4 models, we also use the MLE and QUAD estimators of π to compare the results.

All these models are compared over the 40 evaluated items of each month with the π obtained and the one given by the user. We compare also the order of the 40 items when we rank only those 40 items and when we rank every item of the month. For 40 items, the maximum Kendall distance is 780. In the Table 9, we compare for every model with the 2 different estimators of π the Kendall distance between π and π. Here, we observe that using the real values, the results are better. The MEAN estimator of the logarithm of the normalization term is the best. Concerning the estimator of π for the MEAN model, both QUAD estimator and Maximum Likelihood estimator give similar results.

Month

An interesting property for a prioritization is the indifference of irrelevant alternative, meaning that if we add an element to rank, then the order found without 29 this element has to be the same than the order of the same items found with this element. This property does not stand in the case of Mallows models. However, it can be interesting to have a proximity in the orders given for the full ranking and the orders given with the 40 items. Table 10 Kendall distance between π ranked over every items and π ranked over the 40 evaluated items, on the 40 evaluated items.

In Table 10, we observe for the 4 most interesting models, the Kendall distance between π ranked over every items and π ranked over the 40 evaluated items, on the 40 evaluated items. A distance of 0 may seems interesting but with the Maximum Likelihood Estimator, we have the possibility of a tyranny, which in the case of the Mallows Model case is the tyranny of the third criterion. For the MEAN model with the Maximum Likelihood Estimator, it is not exactly the case as we have a close proximity to the first and third criteria but we do not have the exact same permutation. Having a so close proximity between the 2 different π with the MEAN model and the QUAD estimator is good.

Conclusions

In this article, we focus on a permutation based method to find the best ranking of items to study. We choose the permutation space to use information given in rankings and explain our method to include other types of information given in the value of the variable.

An estimation of the dispersion parameters of the model is given and we propose different algorithms to find the best compromise for the ranking to study, according to the estimation of those parameters. For each estimation of the compromise ranking, the expected advantages and limitations of our choices is explained.

The limitations and the benefits of the three algorithms for the Mallows model is studied on simulated data. The behavior of the estimation of the compromise is studied with respect to the number of ranking to take into account and their dispersion parameter.

Finally, we apply our new model on a real data set and describe our procedure idea to find the dispersion parameters and the best compromise for the user with

Figure 1 :

 1 Figure 1: Permutohedron of order 4 for the Kendall Distance.

  d p (τ, σ j ) .

Figure 2 :

 2 Figure 2: Permutohedron of order 4 for the Kendall Distance.

Figure 3 :

 3 Figure 3: Boxplot of the Kendall distances between π and π for 5 σ j , permutations of n = 30 items, with the same θ and the same π over 3000 realizations.

Figure 4 :

 4 Figure 4: Kendall Distance between π and π for different numbers of σ j , permutations of n = 30 items, with the same θ = 0.3 and the same π over 3000 realizations.

Figure 5 :

 5 Figure 5: Boxplot of the Kendall distance between π and σ for n = 30, in the set of parameters P 4 .

Figure 6 :

 6 Figure 6: Boxplot of the Kendall distance between π and σ for n = 30, in the set of parameters P 8 .

Figure 7 :

 7 Figure 7: Boxplot of the Kendall distance between π and σ for n = 30, in the set of parameters P 9 .

  Month MM + MLE MM + QUAD MEAN + MLE MEAN +

Table 1 :

 1 Means ± standard deviations of Kendall distances between π and π for 50 σ j , permutations with n = 30 items, with the same θ and the same π over 3000 realizations.

	5 when

  ± 2.785 27.86 ± 2.801 25.21 ± 2.571 0.05 26.17 ± 2.616 26.64 ± 2.646 25.13 ± 2.492 0.1 23.95 ± 2.449 24.38 ± 2.491 23.46 ± 2.388 0.15 21.15 ± 2.180 21.55 ± 2.210 20.90 ± 2.158 0.2 18.46 ± 1.929 18.79 ± 1.972 18.31 ± 1.918 0.25 16.17 ± 1.686 16.42 ± 1.707 16.07 ± 1.668 0.3 14.21 ± 1.509 14.41 ± 1.532 14.15 ± 1.505 0.35 12.67 ± 1.331 12.83 ± 1.350 12.63 ± 1.329 0.4 11.38 ± 1.196 11.51 ± 1.210 11.35 ± 1.195 0.45 10.31 ± 1.116 10.42 ± 1.118 10.30 ± 1.108

	Choice of θ	BORDA	MLE	QUAD
	0.01 27.40 0.5 9.43 ± 0.999	9.52 ± 1.005	9.42 ± 0.994
	0.55	8.65 ± 0.920	8.72 ± 0.927	8.65 ± 0.920
	0.6	8.02 ± 0.847	8.07 ± 0.848	8.02 ± 0.844
	Table

Table 4 :

 4 Means ± standard deviations of Kendall distances between π and π for 2 to 50 σ j , permutations of n = 30 items, with the same θ = 0.3 and the same π over 3000 realizations.

		BORDA	MLE	QUAD
	2	51.52 ± 10.300 69.13 ± 14.736 63.59 ± 12.554
	3	42.51 ± 8.530	43.76 ± 9.118	43.83 ± 9.074
	4	37.00 ± 7.409	39.49 ± 8.238	38.81 ± 7.908
	5	32.73 ± 6.766	32.91 ± 7.222	33.10 ± 7.117
	6	29.97 ± 6.166	30.65 ± 6.768	30.40 ± 6.661
	7	27.42 ± 5.642	26.66 ± 5.921	26.94 ± 5.929
	8	25.37 ± 5.320	25.20 ± 5.752	25.27 ± 5.797
	9	23.77 ± 5.037	22.79 ± 5.371	23.12 ± 5.307
	10	22.07 ± 4.728	21.45 ± 5.157	21.53 ± 5.137
	15	17.33 ± 3.992	15.79 ± 4.178	16.25 ± 4.244
	20	14.29 ± 3.588	12.78 ± 3.635	13.02 ± 3.648
	25	12.23 ± 3.255	10.46 ± 3.398	10.82 ± 3.314
	30	10.66 ± 3.055	8.95 ± 2.976	9.15 ± 3.044
	40	8.57 ± 2.655	6.82 ± 2.504	7.05 ± 2.570
	50	7.07 ± 2.431	5.29 ± 2.204	5.49 ± 2.242

Table 5 :

 5 Mean ± standard deviations of the standard deviation of the Kendall distances between σ and π for 2 to 50 σ j , permutations of n = 30 items, with the same θ and the same π over 3000 realizations.

		BORDA	MLE	QUAD
	2	2.40 ± 1.845 50.24 ± 8.921 0.25 ± 0.250
	3	7.01 ± 3.675 11.74 ± 6.181 7.55 ± 4.217
	4	8.81 ± 3.824 13.38 ± 5.848 5.09 ± 3.924
	5	9.92 ± 3.670 12.63 ± 4.659 9.66 ± 3.640
	6	10.73 ± 3.554 13.64 ± 4.446 8.89 ± 3.733
	7	11.25 ± 3.440 13.10 ± 3.899 10.93 ± 3.208
	8	11.74 ± 3.224 13.47 ± 3.670 10.56 ± 3.262
	9	12.16 ± 3.221 13.57 ± 3.581 11.77 ± 2.949
	10	12.22 ± 3.083 13.57 ± 3.408 11.38 ± 3.116
	15	13.16 ± 2.575 13.95 ± 2.724 12.83 ± 2.448
	20	13.58 ± 2.233 14.16 ± 2.331 13.28 ± 2.204
	25	13.77 ± 2.054 14.23 ± 2.127 13.57 ± 2.016
	30	13.84 ± 1.867 14.21 ± 1.930 13.69 ± 1.855
	40	14.07 ± 1.618 14.33 ± 1.658 13.97 ± 1.609
	50	14.22 ± 1.482 14.41 ± 1.498 14.14 ± 1.480
	Choice of θ	BORDA	MLE	QUAD
	P 1	63.08 ± 10.843 64.61 ± 11.530 65.09 ± 11.590
	P 2	1.23 ± 1.060	0.54 ± 0.712	0.62 ± 0.779
	P 3	22.83 ± 5.811	32.67 ± 8.107	30.90 ± 7.436
	P 4	25.20 ± 5.661	27.38 ± 6.716	29.76 ± 6.395
	P 5	26.94 ± 5.456	22.74 ± 5.631	33.13 ± 6.380
	P 6	3.32 ± 1.625	2.04 ± 1.374	2.25 ± 1.437
	P 7	3.51 ± 1.664	2.06 ± 1.372	2.70 ± 1.564

Table 6 :

 6 Mean ± standard deviations of Kendall distances between π and π for a set of parameters θ : P 1 : 50 parameters at 0.05 ; P 2 : 50 parameters at 0.6 ; P 3 : 2 parameters at 0.6 ; P 4 : 2 parameters at 0.6 and 5 parameters at 0.05 ; P 5 : 2 parameters at 0.6 and 48 parameters at 0.05 ; P 6 : 25 parameters at 0.6 ; P 7 : 25 parameters at 0.6 and 25 parameters at 0.05. The permutations are of n = 30 items.

Table 7 :

 7

	Parameters θ	BORDA	MLE	QUAD
	P 1	26.17 ± 2.616 26.64 ± 2.646 25.13 ± 2.492
	P 2	8.02 ± 0.847	8.07 ± 0.848	8.02 ± 0.844
	P 3	1.73 ± 1.344 25.12 ± 5.053	0.25 ± 0.25
	P 4	20.36 ± 6.735 20.97 ± 6.689 19.53 ± 6.405
	P 5	25.99 ± 2.639 26.20 ± 2.657 25.45 ± 2.570
	P 6	7.76 ± 1.187	7.91 ± 1.209	7.70 ± 1.177
	P 7	19.62 ± 2.511 19.66 ± 2.507 19.58 ± 2.505

Mean ± standard deviations of the standard deviation of the Kendall distances between σ and π for a set of parameters θ : P 1 : 50 parameters at 0.05 ; P 2 : 50 parameters at 0.6 ; P 3 : 2 parameters at 0.6 ; P 4 : 2 parameters at 0.6 and 5 parameters at 0.05 ; P 5 : 2 parameters at 0.6 and 48 parameters at 0.05 ; P 6 : 25 parameters at 0.6 ; P 7 : 25 parameters at 0.6 and 25 parameters at 0.05. The permutations are of n = 30 items.

Table 8 :

 8 Mean ± standard deviations of the standard deviation of the Kendall distances between σ and π for a set of parameters θ : P 8 : 8 different parameters, θ 1 = 0.1, θ 2 = 0.15, θ 3 = 0.2, θ 4 = 0.25, θ 5 = 0.3, θ 6 = 0.35, θ 7 = 0.4 and θ 8 = 0.6 ; P

	Choice of θ	BORDA	MLE	QUAD
	P 8	23.95 ± 5.075 21.78 ± 5.368 24.74 ± 5.531
	P 9	38.01 ± 7.309 32.66 ± 7.921 45.98 ± 8.356

9 : 1 parameters at 0.6 and 6 parameters at 0.1. The permutations are of n = 30 items.

Table 9

 9 Kendall distance between π and π when we rank only the 40 evaluated items.

		MM + MLE	MM + QUAD MINMAX + MLE MINMAX + QUAD
	1	160	181	189	183
	2	91	147	121	137
	3	128	135	168	168
	4	137	173	182	198
	Month MED + MLE MED + QUAD	MEAN + MLE	MEAN + QUAD
	1	135	136	128	133
	2	106	106	97	97
	3	125	132	123	125
	4	132	132	135	134


Christophe Negre
email: christophe.negre@univ-perp.fr

Address-Bit Differential Power Analysis on Boolean Split Exponent Counter-Measure

 preventing differential power analysis on address bits. We show that one of their proposed protections has a flaw. We derive an attack exploiting this flaw and we successfully apply it on a simulated power consumption of an RSA modular exponentiation.

INTRODUCTION

Side channel analysis is a serious threat for devices performing cryptographic computation. Specifically, it can threaten exponentiation x κ involved in currently used cryptosystems (RSA, DSA and ECDSA). In 1996, in a seminal work [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF] Kocher showed that a statistical analysis of the computation time can leak out the secret key used in RSA cryptosystem. Afterwards, it was shown that the power consumption can be used to extract secret information: with simple power analysis (SPA) [START_REF] Kocher | Differential Power Analysis[END_REF] one can read on the power trace the sequence of operations (squaring/multiplication) performed during an exponentiation and then deduce the secret exponent. In [START_REF] Kocher | Differential Power Analysis[END_REF] the authors introduce the differential power analysis (DPA) which computes differences on the power traces to determine the bits of the key.

To counter-act these side channel analyses, the basic approach is to randomly mask sensitive data: we can blind the point x by multiplying it by a random element [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF] or use randomized representation [START_REF] Clavier | Horizontal Correlation Analysis on Exponentiation[END_REF], and we can use additive or multiplicative mask of the exponent [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF][START_REF] Tunstall | Coordinate Blinding over Large Prime Fields[END_REF]. But this is not always sufficient, since it was shown [START_REF] Itoh | Address-Bit Differential Power Analysis of Cryptographic Schemes OK-ECDH and OK-ECDSA[END_REF] that the loads and stores performed during the exponentiation induce a power consumption correlated to the corresponding address bits. This can be exploited to mount a differential power analysis on the address bits (ADPA). Recently some strategies [START_REF] Izumi | Improved countermeasure against address-bit DPA for ECC scalar multiplication[END_REF][START_REF] Tunstall | Boolean Exponent Splitting[END_REF] were proposed to counter-act this attack by randomizing the address bits used in the loads and stores of the exponentiation algorithm.

Contribution. In this paper we analyse a potential flaw in the signed version of the boolean exponent splitting approach of [START_REF] Tunstall | Boolean Exponent Splitting[END_REF]. We show that the sequence of operations (squarings and multiplications) are not fully regular: we show that if we can distinguish a squaring from a multiplication we can deduce the bit used to randomize the address bits. Then a classical differential power analysis on the address bits can be performed to determine the secret exponent. We validate this attack by simulating power traces based on the Hamming weight model. We provide a method for distinguishing a squaring from a multiplication which has a high level of confidence. Then we show that an ADPA attack can be performed and successfully deduce the secret exponent with a few thousand power traces.

Organization of the paper. In Section 2 we review power analyses on exponentiation algorithm and related counter-measures. In Section 3 we show that there is a flaw in the counter-measures "boolean splitting exponent" of [START_REF] Tunstall | Boolean Exponent Splitting[END_REF]. In Section 4 we present our approach for power consumption simulation of modular exponentiation, and we apply an attack on the "boolean splitting exponent" exploiting the flaw shown in Section 3. We ends the paper in Section 5, with a few concluding remarks.

REVIEW OF POWER ANALYSES

In currently used public key cryptosystems (e.g. RSA, DSA, ECDSA), the sensitive operation is an exponentiation x κ for x ∈ G a finite group or ring and κ a secret exponent. With the square-and-multiply algorithm one can compute x κ as a sequence of squarings R ← R 2 followed by a multiplication g ← g × x when the i-th bit κ i = 1. This basic approach is not secure when considering side channel analysis. For example, if we monitor the power consumption of an exponentiation, then, a simple power analysis [START_REF] Kocher | Differential Power Analysis[END_REF] can identify on the trace the operations performed (square or multiplication) on the device and then deduce the secret exponent κ. To counteract this attack modified versions of the square-and-multiply exponentiation are recommended for which the sequence of computed operations is not related to bits of the exponent κ. The most popular method is the Montgomery ladder (Joye and [START_REF] Yen | The Montgomery Powering Ladder[END_REF] which involves two variables R 0 and R 1 , for which, during the exponentiation, we always have R 1 = xR 0 and the loop iteration is always a multiplication followed by a squaring (cf. Algorithm 1).

Algorithm 1 Montgomery-ladder

Require: x ∈ G, one n-bit integers κ = ∑ n-1 i=0 κ i 2 i Ensure: x κ 1: R 0 ← 1 G ; R 1 ← x; 2: for i = n -1 down to 0 do 3: R 1-κ i ← R κ i • R ¬κ i ; 4: R κ i ← R 2 κ i ; 5: return R 0
The Montgomery ladder is a good protection against SPA, but it is not robust against attack like differential power analysis (DPA) [START_REF] Kocher | Differential Power Analysis[END_REF], collision attack [START_REF] Fouque | The Doubling Attack -Why Upwards Is Better than Downwards[END_REF][START_REF] Yen | Relative Doubling Attack Against Montgomery Ladder[END_REF]. Indeed these attacks guess one or several bits of the key and they predict the power consumption at some iteration of the exponentation. If this prediction is correct this means that the guessed bits are also correct otherwise one can try another guess. To counteract such attacks at the algorithmic level, the best approaches randomize data in order to render power consumption unpredictable:

• Point blinding. The strategy here is to hide the data x and R 0 , R 1 by either multiplying x by a random element [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF], adding a random mask or by randomizing the representation of x [START_REF] Clavier | Horizontal Correlation Analysis on Exponentiation[END_REF].

• Exponent masking. In this case we randomly modify the exponent. [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF] proposed to add to κ a random mutiple r × N where N is the order of G. One can also [START_REF] Tunstall | Coordinate Blinding over Large Prime Fields[END_REF] randomize the exponent as κ = β -1 κ mod N and compute x κ as (x β ) κ .

Point blinding and Exponent masking counter-measures induce an overhead which, for the later, is important since the level of randomization have to be larger than the longest run of 0 or 1 in κ [START_REF] Smart | Randomised representations[END_REF]. Point blinding approaches used alone are not sufficient to counter-act DPA attack on the address bits (ADPA [START_REF] Itoh | Address-Bit Differential Power Analysis of Cryptographic Schemes OK-ECDH and OK-ECDSA[END_REF]). Indeed guessing a few key bits leads to a prediction of the address bits involved in the loads and stores performed during the exponentiation. Predicting these address bits, leads to predicting the behavior of the power consumption during the loads or stores and performing a DPA on several power traces validates the guessed key bit or not. Since point blinding counter-measures do not modify the sequence of operations done during the exponentiation, they do not protect the implementation from an ADPA.

Consequently, to counter-act such DPA on address bits, the authors in [START_REF] Izumi | Improved countermeasure against address-bit DPA for ECC scalar multiplication[END_REF] proposed a first version of the Montgomery-ladder which randomizes the loads and stores. This work was subsequently improved in [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] which provides two kinds of randomized Montgomery-ladder. In their first approach they use one random bit a i to split the key bit κ i = a i ⊕ b i and which decides the order of loads R 0 and R 1 for the multiplication R 0 × R 1 in Step 3 of Algorithm 1. This bit a i is also used to randomly store the results of the squaring in Step 4 and the multiplication in Step 3 of Algorithm 1 done in one iteration of the Montgomery ladder. This requires a bit b which memorizes how the value are stored in the two registers R 0 and R 1 at the end of a loop iteration. This approach is shown in Algorithm 2.

The authors in [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] propose a second approach (Algorithm 3) which randomizes the signed version of the Montgomery ladder. This algorithm perform the exponentiation as a sequence of two multiplications one between R 0 and R 1 and between R 0 and either U 0 = x or U 1 = x -1 . In Algorithm 3 the addresses used in the Algorithm 2 Montgomery-ladder with XOR split exponent I [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] Require

: x ∈ G, two n-bit integers A = ∑ n-1 i=0 a i 2 i and B = ∑ n-1 i=0 a i 2 i Ensure: x κ where κ = A ⊕ B 1: R 0 ← 1 G ; R 1 ← 1 G ; R 2 ← 1 G ; b ← R {0, 1}; R ¬b ← x; 2: for i = n -1 down to 0 do 3: R 2 ← R a i • R ¬a i ; 4: R a i ← R 2 (b i ⊕b )⊕a i ; 5:
R ¬a i ← R 2 ; 6:

b ← b i ; 7: return R b store instructions are fixed and then they are not correlated to the bits of the exponent. The random bit a i (and thus b i = κ i ⊕ a i ) is used to randomize the load instructions in Step 4 and Step 5. The bit b is used to memorize how the data are placed in the registers R 0 and R 1 at the end of a loop iteration.

Algorithm 3 Montgomery-ladder with XOR split exponent II [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] Require:

x ∈ G, two n-bit integers A = ∑ n-1 i=0 a i 2 i and B = ∑ n-1 i=0 a i 2 i Ensure: x κ where κ = A ⊕ B 1: R 0 ← 1 G ; R 1 ← 1 G ; R 2 ← 1 G ; b ← R {0, 1}; R ¬b ← x; 2: for i = n -1 down to 0 do 3: R 0 ← R b i ⊕b • R (b i ⊕b )⊕a i ; 4: R 1 ← R 0 •U b i ; 5: b ← b i ; 6: return R b

WEAKNESSES OF BOOLEAN SPLIT RANDOMIZATION

In [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] the authors claim that their algorithm combined with a point blinding technique has the same security level but with a lower cost as the method based on randomizing the exponent [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF][START_REF] Tunstall | Coordinate Blinding over Large Prime Fields[END_REF]. For Algorithm 2 and 3, this is not entirely true:

• Randomizing the exponent alone prevents from the three attacks : DPA, collision and ADPA. This approach ensure that the sequence of points R 0,i and R 1,i in G for i = n -1, . . . , 0 computed during the exponentiation are always different. This renders impossible to predict the power consumption of operation involving R 0,i and R 1,i and then prevents DPA and collision attack. Randomizing the exponent also randomly changes the bits of the exponent, which implies that the stores and loads are also randomized and thus unpredictable.

• Randomizing the loads and stores like it is done in Algorithm 2 and 3 does not modify the sequence of computed element R 0,i and R 1,i of G. So a correct guess of the bits of the exponent would lead to a correct guess of R 0 and R 1 at the considered iteration. So if we want to hide these values we have to inject a high level of randomization otherwise a small amount of leakage could be exploited and would lead to a successful DPA. The only approach which produces a randomized and unpredictable sequence of computed points R 0,i and R 1,i in the exponentiation is the point blinding of [START_REF] Coron | Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems[END_REF] which multiply x by a random value.

The table below summarizes the strength of the counter measures reviewed in Section 2 when they are used alone. This table shows that exponent randomizations are the most robust methods. Algorithm 3 has a more important flaw. The main problem of Algorithm 3 is that in Step 3, a multiplication of two different data R 0 and R 1 is done when a i = 1 but when a i = 0 this multiplication is either R 0 × R 0 or R 1 × R 1 , which are squarings done using a multiplication routine. But distinguishig such a squaring done with multiplication routine from a genuine multiplication can be detected by power analysis. In [START_REF] Hanley | Using templates to distinguish multiplications from squaring operations[END_REF] the authors showed that the power consumption of these operations can be distinguished using template methodology.

If we successfully distinguish the power trace of R 0 × R 1 from either R 1 × R 1 or R 0 × R 0 , for unknown R 0 and R 1 , we can deduce the value of a i for i = n -1, . . . , 0. Then the randomization of loads and stores is broken and DPA can be conducted on the address bits of Algorithm 3.

EXPERIMENTATION

In this section we present experimental results of the proposed attack on address randomization of Algorithm 3. The attack is performed on a simulated power traces using the hamming weight leakage model.

Simulation of power consumption

We target an RSA modular exponentiation with an RSA modulus N of bit length 1024. We consider an implementation of the modular exponentiation RSA based on a word-level Montgomery modular multiplication algorithm [START_REF] Bosselaers | Comparison of Three Modular Reduction Functions[END_REF]. Then in order to simulate the power consumption of a full modular multiplication we only need to generate the power consumption of w-bit additions and multiplications. for i = 0 to 32 do r ← r + 2 i (b i × a) The i-th bit r i is output after the i-th clock-cycle. The additions in Fig. 1 are done through carry save adder in order to reduce the critical path delay. This means that the carries are not propagated and but are saved in the flip-flops c i for i = 0, . . . , 31. After the first 32 iterations the first 32 bits r 0 , . . . , r 31 of r are generated and output. But it remains to perform 32 more iterations, with input 0 in place of b i ,to propate the carries and generate r 32 , . . . , r 63 .

We do not provide the circuit of the 32-bit adder, since it is a classical 32 bit adder, which can be easily found in the literature.

To get the power consumption of one clock-cycle of the multiplier we compute the hamming weight of signal flowing in the wires. We split the clock cycle into four parts P 1 , P 2 , P 3 and P 4 . We assume that at beginning of a cycle there are only 0 on all wires. Then we propagate data from the flip-flop and deduce the simulated consumption:

• Each wire containing a signal 1 and connecting a flip-flop and a first gate contributes to 1 on the consumption of P 1 , P 2 , P 3 and P 4 .

• Each wire containing a 1 and connecting a first gate and a second gate contributes to to the consumption of P 2 , P 3 and P 4 .

• Each wire containing a 1 and connecting a second gate and a third gate contributes to 1 to the consumption of P 3 and P 4 .

• Each wire containing a 1 and connecting a third gate and the flip-flop contributes to 1 to the consumption of P 4 . The above simulation of the consumption is highly simplified but we believe that it is sufficient to validate the potential threat of the proposed attack on Algorithm 3. We show in Fig. 2 an example of the simulated power consumption of one w-bit multiplication.

Word level Montgomery modular multiplication. For a modular multiplication we use the word level version [START_REF] Bosselaers | Comparison of Three Modular Reduction Functions[END_REF] of the Montgomery modular multiplication with word size w. Given two integers A and B consisting of s words the modular multiplication are performed as follows 1: R ← (0, . . . , 0) 2 w 2: for i = 0 to s -1 do

3: q ← (R + A[i] × B) × N mod 2 w 4: R ← (R + q × N + A[i] × B)/2 w
Step 3 consists of two w-bit multiplications and one w-bit addition. Step 4 involves two products of a ws-bit integer by a w-bit integer and two additions of two ws-bit integers. Each of these operations are computed through a sequence of w-bit multiplications and/or additions.

Distinguishing a square from a multiplication

We would like to distinguish a multiplication (A × B) mod N from a squaring (A × A) mod N. The main idea to get such distinguisher is that during the multiplication of (A × A) mod N for each i = j the same product is done twice as

A[i] × A[ j] and in reverse operand A[ j] × A[i] in
Step 3 the word level Montgomery multiplication. But since these multiplications involve the same data their power trace must be correlated.

Then we tried to compute the covariance of the power traces of the 32-bit multiplication in order to determine if they are correlated or not. But this strategy was not successful. Probably this is due to the non symmetric form of the considered 32-bit multiplier. We tried another strategy where we evaluate the difference of the mean of the power traces:

σ i, j (A, B)=| Tr(Mult(A[i], B[ j])| -| Tr(Mult(B[i], A[ j])| .
This computed value might be low if A = B, since, in this case, many computed bits appearing in the multiplication Mult(A[i], A[ j]) also appear in Mult(A[ j], A[i]) and they cancel out in the difference of σ i, j (A, A). If A = B these bits do not cancel out resulting in a higher value for σ i, j (A, B). This effect is amplified if we add up σ i, j (A, B) over all i < j:

σ(A, B) = ∑ 0≤i< j≤s σ i, j (A, B).
In order to check this fact we performed the following experiment: for a sufficiently large number of time ( ∼ = 10000) we chose randomly A and B and computed σ(A, B) and σ(A, A). We obtained the resulting density distribution function shown in Fig. 3. We chose k = 3 leading to a level of confidence in our distinguisher of 1 -1/9.This is validated by an experiment for 5000 samples for which we get 99.5% of success.

FINAL STEP : ADPA

Now, since we have a reliable distinguisher, we can proceed to the next step: finding the bits of the exponent with a DPA attack on the address bits. Given a power trace of an exponentiation x κ mod N we can find the random bit a i used in Algorithm 3 for all i = 0, . . . , n -1. Our goal now is to find difference of the power trace related to the address bits involved in loads or stores. Let us analyze the loads and stores in order to select interisting power traces and compute a difference producing a peak if the guess is incorrect and a flat trace if the guess is correct. The table below shows the relation between the (i + 1)-th iteration and i-th iteration for the operation done in Step 5 in Algorithm 3 depending on the value of b , k i and a i . If we focus on the cases a i = 0 (the row in bold in

= 0 R 1 ← R 0 U 0 k i = 0, a i = 0 R 1 ← R 0 U 0 b = 0 R 1 ← R 0 U 0 k i = 0, a i = 1 R 1 ← R 0 U 1 b = 1 R 1 ← R 0 U 1 k i = 0, a i = 0 R 1 ← R 0 U 0 b = 1 R 1 ← R 0 U 1 k i = 0, a i = 1 R 1 ← R 0 U 1 b = 0 R 1 ← R 0 U 0 k i = 1, a i = 0 R 1 ← R 0 U 1 b = 0 R 1 ← R 0 U 0 k i = 1, a i = 1 R 1 ← R 0 U 0 b = 1 R 1 ← R 0 U 1 k i = 1, a i = 0 R 1 ← R 0 U 1 b = 1 R 1 ← R 0 U 1 k i = 1, a i = 1 R 1 ← R 0 U 0 the table)
, we can notice that:

• Case 1: b ⊕ k i = 0. The operations in Step 5 of the two iterations are the same. In particular the address for U i is the same.

• Case 2. b ⊕ k i = 1. The operations in Step 5 of the two iterations are different. In particular the address bits for U i are different.

Consequently, we proceed by guessing the value of g = k i ⊕ k i+1 . Then we select the power traces such that a i = 0 and such that a i+1 ⊕ g = 0. If our guess is correct we will have:

0 = a i+1 ⊕ g = a i+1 ⊕ k i+1 ⊕ k i = b i+1 ⊕ k i = b ⊕ k i .
This means that if the guess g is correct, then Case 1 applies, which means that computing the difference of the power traces of loop i and i + 1 would lead to a zero difference for the address of U i , and the difference will be flat. If the guess is not correct, we are in Case 2, and the difference would be equal to the consumption of the address of U 0 minus the one for U 1 . Adding a sufficient number of such differences would lead to a flat difference if the guess is correct and a peak if it is not correct. In Fig. 4 and Fig. 5 we provide the ADPA obtained for an RSA of size 1024 bits and using 4000 traces. We can see the peak in Fig. 5 showing that the guess g is not correct. There is no peak in Fig. 4 which means that the guess is correct.

This experimentation shows that the proposed approach is effective to extract the whole key. This means that Algorithm 3 does not provide the claimed protection from ADPA. This attack works even if the elements are blinded at the beginning of the exponentiation, by either a randomized representation or a multiplication with a random element.

CONCLUSION

In this paper we considered two exponentiation algorithms (Algorithm 2 and 3) proposed in [START_REF] Tunstall | Boolean Exponent Splitting[END_REF] with randomized store and load in order to counter-act address bit differential power analysis. We analyzed the security of these approaches, and we showed that Algorithm 3 has a an important flaw. Indeed, the operation done in Step 3 of Algorithm 3 is a square or multiplication depending on the bit used for load and store randomization. With a simulated power consumption we showed that we can distinguish a square from a multiplication. This means that the randomization of loads and stores in Algorithm 3 is not effective anymore and an ADPA can be conducted to recover the whole secret key with a few thousand power traces.

Figure

  Figure 1: 32-bit sequential multiplier

Figure 2 :

 2 Figure 2: Power trace of a 32-bit multiplication

Fig. 3

 3 clearly shows that we can easily distinguish a squaring from a genuine multiplication. We estimated the mean value m and the standard deviation d for the two cases: squaring (S) and genuine multiplication (M):m S = 310000, d S = 28000 m M = 710038, d M = 82000Given a power trace for a multiplication A × B mod N, we compute σ(A, B) and we deduce that A = B if σ(A, B) is close to m S and A = B if it is close to m M . To a get an estimation of the probability of success of this approach,

Figure 3 :

 3 Figure 3: Square and non square repartition of power consumption .

Figure 4 :Figure 5 :

 45 Figure 4: Difference obtained for a correct guess on k i+1 ⊕ k i