A new decision making model based on Rank Centrality for GDM with fuzzy preference relations

scarcity of research dealing with the direct approaches. In this paper we present a direct approach towards aggregating several fuzzy preference relations on a set of alternatives into a single weighted ranking of the alternatives. By mapping the pairwise preferences into transitions probabilities, we are able to derive a preference ranking from the stationary distribution of a stochastic matrix. Interestingly, the ranking of the alternatives obtained with our method corresponds to the optimizer of the Maximum Likelihood Estimation of a particular Bradley-Terry-Luce model. Furthermore, we perform a theoretical sensitivity analysis of the proposed method supported by experimental results and illustrate our approach towards GDM with a concrete numerical example. This work opens avenues for solving GDM problems using elements of probability theory, and thus, provides a sound theoretical fundament as well as plausible statistical interpretation for the aggregation of expert opinions in GDM.


Introduction
Group decision making (GDM) settings involve a group of individuals (experts), where each member of the group expresses her preferences over a set of alternatives.Illustrative examples range from parliamentary groups working to converge on a political decision, to groups of friends deciding on the best choice of restaurant for a dinner.The aim of GDM is to identify the most preferred alternative for the whole group of individuals, or to derive a ranking of the alternatives that reflects the preferences of the group.
The literature proposes many different forms of expressing preferences of experts ( Capuano, Chiclana, Fujita, Herrera-Viedma, & Loia, 2017 ).Some of the most popular ones are the following: • Rankings , which are ranked lists of the alternatives from the most preferred to the least preferred one ( Seo & Sakawa, 1985 ). • Utility vectors , where each component of the vector describes the utility of the corresponding alternative, which can be seen as its ordinal strength ( Tanino, 1990 ).These are sometimes called priority vectors or weighted rankings .We use the latter expression throughout this paper.• Preference relations , where preference is expressed as a binary relation on the set of alternatives ( Kitainik, 2012 ).• Fuzzy Preference Relations (FPRs) , which relax the binary preference relations with the possibility of expressing degrees of preference among the alternatives ( Pedrycz, Ekel, & Parreiras, 2010 ).Preference degrees can be assessed using linguistic term sets which can be more natural for human expert to articulate ( Ureña, Kou, Wu, Chiclana, & Herrera-Viedma, 2019 ).
FPR is the most commonly used representation for expressing preferences in a set of alternatives, in which an expert expresses her preferences as degrees of preference assigned to each pair of alternatives.The most common way to store these pairwise preferences is in the form of a preference matrix.The main reason behind the popularity of the preference relations comes from a known fact from psychology studies that human beings are better at comparing pairs of alternatives than at coming up with a complete preference ordering of a set of alternatives ( Ureña, Chiclana, Morente-Molinera, & Herrera-Viedma, 2015 ).
Within GDM involving FPR, there are two main families of approaches: direct approaches and indirect approaches.The indirect approaches first compute the group opinion in the form of an FPR (we will call it a group or a collective FPR), usually expressed as a preference matrix, and then find a solution which is a (weighted) ranking of the different alternatives based on the collective FPR.The collective preference matrix is generally derived by applying an aggregation operator to the individual ones.On the other hand, the direct approaches do not involve constructing a collective preference matrix describing the group opinion as an intermediate step.They first compute an individual ranking for each expert based on her FPR, and then the group ranking is obtained from the individual rankings of the experts using an aggregation operator.For excellent surveys on consensus processes and preference aggregation we refer the reader to the comprehensive surveys ( Cabrerizo, Moreno, Pérez, & Herrera-Viedma, 2010;Herrera-Viedma, Cabrerizo, Kacprzyk, & Pedrycz, 2014 ) and to the book by Herrera-Viedma et al. (2011) .
While studies on indirect approaches for aggregating pairwise preferences abound, the direct approaches are not as popular, although there are a few exceptions ( Dong, Xu, & Yu, 2009;Fan, Ma, Jiang, Sun, & Ma, 2006 ).Herrera and his collaborators ( Herrera, Herrera-Viedma, & Verdegay, 1996 ) pioneered the first direct approach towards GDM based on FPR.However, the work in this direction is very scarce, although it is known that direct approaches usually possess two desirable properties, internal consistency and Pareto principle of the social choice theory ( Dong & Zhang, 2014 ).
One of the few available direct approaches in the literature was recently presented by Dong et al. in ( Dong & Zhang, 2014 ).There, the authors extended the original direct approach presented in ( Herrera et al., 1996 ) in order to support (i) different preferences representations, and (ii) a consensus process in the form of rounds where experts are required to adjust their pairwise preferences.Interestingly, in order to achieve consensus, Dong et al. resort to a form of a feedback based on measuring consensus using the individual weighted rankings of the experts.This is distinct from the main stream of research in FPR since consensus degree computation is not based on weighted rankings of individual experts but rather based on elements from the preference matrices.The approach by Dong et al. allows the experts to update their preference matrices in order to reach a consensus, defining two quantities, namely the cardinal consensus degree, based on the vector representation inspired by ( Chiclana, Herrera, Herrera-Viedma, & Martínez, 2003;Dong, Xu, Li, & Feng, 2010 ) and the ordinal consensus degree inspired by ( Herrera-Viedma, Alonso, Chiclana, & Herrera, 2007;Herrera-Viedma, Herrera, & Chiclana, 2002 ).
In this article, we take a direct approach towards group decision making given fuzzy preferences over a set of alternatives.We propose a method for aggregating the opinions of several experts, which are expressed as FPRs, into a single weighted ranking of the alternatives.Similarly to the work in ( Dopazo & Martínez-Céspedes, 2017 ), we transform the preference matrices into stochastic matrices, and then use the theory of Markov chains and random walks to compute rankings over the alternatives, as implemented in the PageRank algorithm ( Gleich, 2015 ).One main difference in this paper compared to the method in ( Dopazo & Martínez-Céspedes, 2017 ) lays in the definition of the stochastic matrix.In ( Dopazo & Martínez-Céspedes, 2017 ) the stochastic matrix is simply a column-normalization of the preference matrix so that its entries are proportional to the corresponding preferences and represent the probabilities (relative strengths) of dominance between the alternatives.In our framework we determine the entries of the stochastic matrix similarly as in ( Negahban, Oh, & Shah, 2012;2016 ) and they represent the probabilities of transiting between the corresponding alternatives in the way that the probability of transition from alternative x to alternative y is proportional to the degree of preference of y to x .The stationary vectors, however, have similar interpretation in both our approach and the approach in ( Dopazo & Martínez-Céspedes, 2017 ) as their entries represent preference strengths of the corresponding alternatives in both the cases.The difference is that the normalization we use leads to a stationary vector that satisfies the global balance property with respect to the preference matrix: the preference strength of an alternative depends on whether the alternative dominates weak or strong alternatives.This is the core idea of the Rank Centrality method ( Negahban et al., 2012;2016 ) and we discuss it in more detail in Section 3.1 .
Notice that assigning and interpreting a degree of preference is not straightforward.Using a probability value to quantify an FPR gives an intuitive interpretation of FPR itself and, moreover, enables to establish a link between probability theory and preference aggregation.Furthermore, we prove that the weighted ranking obtained as a result of the method presented in this paper corresponds to the result of Maximum Likelihood Estimation (MLE) of the Plackett-Luce model ( Plackett, 1975 ).
There is a body of literature on methods that compute weighted ranking from preference matrices based on optimization techniques such as least square method ( Gong, 2008 ), least deviation method ( Xu & Da, 2005 ), multiobjective optimization ( Fernandez & Leyva, 2004 ), new fuzzy linear programming method (FLPM) ( Zhu & Xu, 2014 ), goal programming ( Fan et al., 2006 ), etc.Although these methods are shown to provide good results, they relay on human-engineered techniques or heuristics and do not provide plausible theoretical interpretation of their computation and modelling steps.Our work is distinct from the latter works as the way we derive a ranking vector can be explained using probability theory and thus we provide a theoretical interpretation of our framework.
Indeed, as mentioned above, our work is directly inspired by a recent work on Rank Centrality algorithm ( Negahban et al., 2012;2016 ), which aggregates a set of pairwise comparisons of alternatives into a global weighted ranking.In ranking based on pairwise comparisons, the goal is to rank, for example, football teams based on results of played matches between them.This problem has an obvious analogy with ranking of alternatives based on pairwise expressed preferences, but despite the vast amount of work on ranking alternatives based on preferences, to the best of our knowledge, the ideas of Rank Centrality have not yet been adopted in the context of fuzzy preference aggregation.We argue that by properly transforming the fuzzy preferences into probabilities of transition between alternatives, probability theory can naturally be applied in preference aggregation and, consequently, we hope that our framework can inspire further research in that direction.
The remainder of this paper is structured as follows.In Section 2 we introduce the relevant background and the fundamental concepts of the state of the art in fuzzy preferences.Section 3 presents our approach for aggregating preferences based on the concept of Markov chains and discusses the conditions under which some desirable properties for the aggregation processes hold.In Section 4 we provide a theoretical sensitivity analysis of the proposed aggregation method.Concrete numerical example showing the consistency of our framework followed by an experimental sensitivity analysis is given in Section 5 .Final discussions and conclusions are provided in Section 6 .

Background and preliminary concepts
In the following we assume that E = { e 1 , . . ., e m } is a set of experts and X = { x 1 , . . ., x n } is a set of alternatives.We use the following definition of a fuzzy preference relation as provided in ( Herrera-Viedma, Herrera, Chiclana, & Luque, 2004 ).
A Fuzzy Preference Relation (FPR) P on a set of alternatives X is a fuzzy set on the product set X × X, i.e. a relation on X characterized by a membership function (1) p i j = μ P (x i , x j ) is interpreted as the preference degree of the alternative x i over the alternative x j .A usual natural assumption is that p i j + p ji = 1 , for every i, j ∈ { 1 , . . ., n } , i.e. that P is additive reciprocal .If p i j > 0 .5 , we say that x i is preferred to x j ; if p i j = 0 .5 , we say that we are indifferent between x i and x j ; and p i j = 1 indicates that x i is absolutely preferred to x j .The additive reciprocity property ensures that p ii = 0 .5 and p i j > 0 .5 iff p ji < 0 .5 .
When the set X is not too big, it is convenient to represent P as an n × n matrix of preference values, where n = | X| , i.e.P = [ p i j ] n ×n .We call this a preference matrix .For convenience, we use the same notation for both the fuzzy preference relation and the corresponding preference matrix.
We assume that each of the m experts expresses her preferences independently of each other and in the form of a fuzzy preference relation.Let us denote by P (k ) the FPR of the k -th expert and let P (k ) = [ p (k )  i j ] n ×n be the corresponding preference matrix.An indirect approach to GDM aims at reaching a collective opinion by first aggregating all the individual preference matrices into a collective FPR.A direct approach would predict the collective opinion or a GDM-"target", in general, by turning experts' FPR matrices into vectors whose entries measure the ranking of the alternatives.More formally, a weighted ranking can be defined as a function r : X → R , which maps each alternative in X into its absolute preference strength. 1Aggregating the individual ranking vectors yields a possible consensus target for the set of individuals.
Whatever approach one chooses, direct or indirect, there are two main phases in GDM based on FPR, an aggregation phase and an exploitation phase .
In the aggregation phase, the corresponding individual preference values (corresponding entries in FPR matrices or ranking vectors) are aggregated into a collective preference value using an aggregation operator .
There are many aggregation operators, such as weighted average (WA), fuzzy majority, etc.One popular example is the operator called Ordered Weighted Averaging (OWA) due to ( Yager, 1988 ).The WA and the OWA operators presume we have a list of weights w = (w 1 , . . ., w m ) , w k ∈ [0 , 1] , k = 1 , . . ., m , such that w k = 1 .Let p 1 , . . ., p m be a list of preference values to be aggregated.While the WA operator is defined as a simple weighted average of the preferences: 1 The ordering of the alternatives is implicit in each weighted ranking and follows from the linear order of the real numbers: Alternatives assigned a higher number rank higher, and alternatives assigned the same number rank the same.
the OWA operator is defined as where σ is a permutation of the set { 1 , . . ., m } such that p σ (k +1) ≥ p σ (k ) for k ∈ { 0 , . . ., m − 1 } .In the case of WA, the weights w = (w 1 , . . ., w m ) can be assumed to be corresponding to the importance of the experts in the group with respect to the particular decision making problem, or to the confidence of the experts in their opinions.In the case of OWA, assigning weights w enables weighting differently preferences of different strength, giving more value to stronger preferences, for example.By choosing different w in WA and OWA, one can implement different aggregation operators.The exploitation phase is the phase of deducing a (weighted) ranking vector based on a fuzzy preference relation (matrix).
Two relevant approaches towards defining the ranking in this phase are ( Herrera et al., 1996 ): Quantifier-Guided Dominance Degree (QGDD), where the rank of each alternative represents the dominance or importance of the alternative over the rest of the alternatives; and Quantifier-Guided Non-Dominance Degree (QGNDD), where the rank of each alternative represents the degree to which the alternative is not dominated by the rest of the alternatives.An alternative to QGDD and QGNDD is the Netflow method ( Bouyssou, 1992 ) which is also based on dominance of an alternative.2More precisely, this method defines the rank of an alternative as the difference between the inflow and the outflow of preference from it, which, under the additive reciprocity assumption ( p i j + p ji = 1 ), reduces to the following expression: Notice that the Netflow method is related to the Copeland voting rule in ( Marchant, 1996 ).An axiomatic characterization of the Copeland rule can be found in ( Henriet, 1985 ).
The indirect and the direct approaches towards GDM with FPR differ in the way the phases of aggregation and exploitation are combined: In the indirect approaches, one first applies the aggregation phase to the set of individual FPRs in order to obtain one collective FPR.Then the exploitation phase is applied to the collective FPR to give the final (weighted) ranking.In the direct approaches, the exploitation phase is applied first at each individual FPR to give the individual rankings.Then the aggregation phase is applied at the individual rankings to find the final collective ranking ( Herrera et al., 1996 ).Fig. 1 highlights these differences between the two approaches.
Note that aggregating opinions of experts into a group opinion by any of the above approaches does not necessarily amount to reaching a consensus.To guarantee agreement between the experts, one could consider a third phase after the two phases described above, in which one forces the individual opinions of the experts to get close to each other ( Palomares, Estrella, Martínez, & Herrera, 2014 ).This narrowing phase can be implemented in two different ways: (i) through automatic approaches without expert feedback, or (ii) through approaches with feedback on preferences where there is a moderator that will try to reduce the divergence of opinions between the experts.In both of these approaches consensus is considered to be a stage where all the individual opinions are sufficiently close to each other.Quantitatively, the closeness can be measured as a distance in the "space of opinions" either between each opinion and the collective (aggregated) opinion, or between pairs of individual opinions.
Having described briefly the decision framework and the main approaches for quantitative analysis in GDM, in the next section we propose a new method to address the aggregation and exploitation processes and we analyse its properties.

A rank centrality-based preference aggregation method
In this section we introduce a method for aggregating the fuzzy preference relations of the experts e 1 , . . ., e m into a collective weighted ranking of the alternatives x 1 , . . ., x n .We start by providing a way of transforming a preference matrix over the alternatives into a weighted ranking vector of the alternatives, which will be used in the exploitation phase of the general method.Then we provide the full GDM procedure.At the end we prove that, under the assumption that detailed balance relation is satisfied, the proposed method satisfies some desirable properties of aggregation processes.

From preference matrices to ranking vectors
We consider an n × n preference matrix P over the alternatives x 1 , . . ., x n ( n ≥ 2 ) which is additive reciprocal, i.e. a matrix whose elements p i j are in the unit interval and obey the condition p i j = 1 − p ji . (5) A value p i j > 0 .5 means that the alternative x i is preferred over x j , while p i j = 0 .5 means that no preference between x i and x j exists.
Inspired by recent work on ranking based on a dataset of pairwise comparisons ( Negahban et al., 2012;2016 ), our approach is based on a transformation of the given preference matrix P into a stochastic matrix S defined in the following way: ) The division by n − 1 is introduced for normalization purposes, j s i j = 1 , and to guarantee that each element s i j is proportional to the corresponding p ji and fulfills the condition 0 ≤ s i j ≤ 1 .Each row can then be seen as a probability distribution and the matrix S as the matrix of transition probabilities of a Markov chain with n states. 3The element s i j corresponds to the probability of transiting from a state x i to a state x j and, as defined in Eq. ( 6a) , this probability equals the product of the probability of choosing randomly the state x j among the n − 1 states different from x j and adopting that state with probability p ji .This ensures that the probability of transition from alternative x i to alternative x j is proportional to the preference of x j over x i .Notice that it is always possible to construct such a stochastic matrix, even if there are missing entries in the preference matrix P , i.e. if there is incomplete information ( Herrera-Viedma et al., 2007 ).
We impose that the preference matrix P fulfills the condition for every i, j ∈ { 1 , . . ., n } , and therefore also p i j = 1 , for every i, j ∈ { 1 , . . ., n } , which means that an alternative is never completely "excluded" against another one and also never "fully dominates" another one.(This may be achieved by replacing a zero preference with an arbitrarily small > 0 preference).It is easy to see that under this particular condition the matrix S is irreducible and aperiodic (see Appendix A ).Then, according to Perron-Frobenius theorem ( Horn & Johnson, 1990; MacCluer, 20 0 0 ), S being an irreducible aperiodic stochastic matrix, there is a unique stationary solution π satisfying: The stationary distribution π can be computed iteratively via a random walk on the Markov chain defined by the stochastic matrix S, or analytically via the computation of the eigenvector associ-ated with the highest eigenvalue of the matrix S. We will interpret the stationary distribution π as a weighted ranking, where each component of the vector π (each "weight") represents the importance of the corresponding alternative with respect to the whole set of alternatives, and we will refer to this method for computing rankings as a Rank Centrality (RC) method ( Negahban, Oh, & Shah, 2016 ).We call the corresponding matrix S a centrality matrix .
To justify the above interpretation of the stationary vector π , we analyse its relation with the initially given preference matrix P .
According to Eq. ( 8) , the components of the vector π read: The product π j s ji in the above sum represents the probability of transitioning from x j to x i adjusted by the weight of x j and, similarly as in ( Dopazo & Martínez-Céspedes, 2017 ), can be interpreted as the relative importance of alternative x i with respect to the alternative x j .Then, since π i , according to Eq. ( 9) , is a sum of such products, it can be interpreted as the (absolute) importance of x i .Now, according to Eq. ( 6a) , the transition probabilities s ji are proportional to preferences p i j , hence we can interpret π j s ji as relative preference strength of x i with respect to x j , and π i as (absolute) preference strength of x i .
There are two important properties of the RC method.The first property is that the ranking dynamics follows a continuous time Markov chain in global balance, and the respective balance equation can be easily derived from Eq. ( 6) and Eq. ( 8) .Namely, from Eq. ( 9) we obtain: From Eq. ( 6) we have which together with Eq. ( 11) and Eq.(6a) leads to the following equations which we call the global balance equations.Notice that the stronger condition which assumes a term-by-term equality between the sums in Eq. ( 13) , for every i ∈ { 1 , . . ., n } , is called a detailed balance .
We discuss the case of a detailed balance in a subsequent section.
The second property is given by the relation between the components of the stationary solution (ranking) and the initial preference degrees that follows directly from the second global balance equation: Note that if we did not have the terms π j on the right handside of the Eq. ( 14) , then the method would have been equivalent to the Netflow method since π i would be proportional to n j = i p i j which is the case for NF (x i ) too as seen in Eq. ( 4) .However, having π j on the right-hand side of Eq. ( 14) reflects the idea of centrality ranking.Namely, the importance of each alternative x j = x i quantified by π j is taken into account when determining π i .Thus, weak alternatives that have low π j due to being dominated by many other alternatives will not contribute much to the increase of π i even if p i j is large, because one needs to consider the product p i j π j and not only p i j as in the classical Netflow method.Informally, this means that beating weak alternatives, i.e. alternatives with low π j , does not increase much the ranking.This is one of the core ideas of the Rank Centrality method ( Negahban et al., 2016 ).

The GDM method
Our framework for preference aggregation consists of the following subsequent steps: 1. Consider a set of m experts E = { e 1 , . . ., e m } .Each expert e k , 1 ≤ k ≤ m , has a pairwise preference matrix P (k ) over the set of alternatives X = { x 1 , . . ., x n } .
2. Using Eq. ( 6) , for each matrix P (k ) we compute the corresponding stochastic matrix S (k ) .
3. We solve Eq. ( 8) for each expert e k , i.e. we solve π = π S (k ) , for k = 1 , . . ., m , and denote the unique solution by π (k ) .The vector π (k ) = [ π (k )  1 , . . ., π (k )   n ] defines a weighted ranking of the alternatives corresponding to the preferences of the expert e k over the set of n alternatives.As observed in the previous section, π (k )   i can be interpreted as the preference strength of the alternative x i according to expert e k .Then, since π (k ) is a probability distribution over X, it can be seen as representing the expert e k 's distribution of preference strengths over X. 4. We define the collective ranking vector as the arithmetic average of the individual ranking vectors, determining its components as follows: Notice that, by taking the arithmetic average of the individual ranking vectors as the aggregated state π (c) , one naturally defines the consensual stage with respect to that aggregated state: perfect consensus occurs when all the experts arrive at the same opinion given by π (c) .Moreover, such definition of perfect consensus state reflects the assumption that the perfect consensus state should not be near to some specific expert and in prejudice to the other.
From a more physical perspective, one can take π (c) as the center of mass of the set of opinions of the m experts (expressed as weighted rankings), in the space of all the possible opinions, i.e. it is the unique point from which all the experts' opinions are seen equally distributed, In other words, using the arithmetic average in the aggregation phase of the GDM process provides a safe alternative to reaching a consensus, in the absence of a special third phase for that purpose in the process.An extension to WA or OWA is straightforward and does not compromise the bulk of our framework: For any choice of weights, which should depend on the specific set of experts and/or the specific decision making context, rankings and preferences will still be treatable within our framework.

The detailed balance case
Recall that the stationary vector satisfies the detailed balance (time reversibility) property if the following equation holds: (16) for every i, j ∈ { 1 , . . ., n } .
In this section we assume that Eq. ( 16) holds for the elements of the stochastic matrix S and the components of the corresponding stationary vector obtained in steps 2. and 3. in the previous section.Then, combining Eq. ( 5) , Eq. ( 6) and Eq. ( 16) , one arrives to the following equation equivalent to Eq. ( 16) : Eq. ( 17) highlights the correspondence between pairwise expressed preferences in the form of a preference matrix and the preference strengths assigned to the alternatives in the corresponding weighted ranking vector.Note that Eq. ( 17) implies transitivity of the preference matrix P (from p i j > 0 .5 and p jk > 0 .5 follows p ik > 0 .5 ).Notice that Eq. ( 17) means that the preference relation P satisfies the Bradley-Terry-Luce (BTL) model ( van Berkum, 1997 ).It was demonstrated also by Rajkumar and Shivani ( Rajkumar & Agarwal, 2014 ) that the time-reversibility of S is equivalent to the existence of an underlying BTL model that describes the preferences in P according to their preference strengths π .
At this juncture, we shall present a theorem that provides a mathematical interpretation of the stationary distribution of the centrality matrix S corresponding to the preference matrix P , namely that the stationary distribution of S is a maximum likelihood estimate of the BTL model underlying P .
Theorem 1.Let the BTL pairwise comparison model be defined by , where p i j = P (x i > x j ) . 4 Let P = [ p i j ] be a preference matrix and let S be its corresponding centrality matrix defined by Eq. (6a) and Eq. ( 6b) .Then the Maximum Likelihood Estimate (MLE) of the BTL model satisfies the global balance equation with P .
The proof of the above theorem can be found in Appendix B .
The theorem states that the MLE of the parameter vector π of the BTL model for pairwise comparisons of n alternatives, π * , satisfies the global balance property given in Eq. ( 13) with the ground truth probabilities P of the model.This means that π * is a stationary distribution of the centrality matrix S corresponding to P , since the global balance equations in Eq. ( 13) are derived from the stationary distribution in Eq. ( 9) through a series of equivalence steps.Finally, since the stationary distribution of S is unique on the unit interval, we can conclude that it is equal to π * , the MLE of the BTL model in which P describes the probabilities of the pairwise comparisons.
In Section 3.1 we discussed two general properties of the RC method.Here, we observe that in the special case when the ranking vectors determined in step 3. of the GDM procedure satisfy detailed balance, i.e. the equations ( 16) and ( 17) , two desirable properties of the aggregation processes are satisfied: internal consistency and the Pareto principle.We interpret these properties similarly as in ( Dong & Zhang, 2014 ) and ( Chiclana, Herrera, & Herrera-Viedma, 2002 ).
Internal Consistency.In our setting, internal consistency can be understood as the consistency of the process that transforms each expert's opinion from its initial form of a preference matrix, through a stochastic matrix, to its end form of a weighted ranking vector.In other words, the property of internal consistency is satisfied if the following holds: The derived individual ranking of the k -th expert by the procedure described in Section 3.1 , reflects her initial preference relation, ranking higher (assigning higher weights to) the alternatives she prefers more: p (k )  i j ≥ p (k )   ji if and only if 4 x i > x j can be interpreted as "x i wins over x j ", "x i is preferred to x j ", etc., and ) is the probability of this event that is to be estimated from a number of pairwise comparisons in the set of alternatives X = { x 1 , . . ., x n } .
π (k )   i ≥ π (k ) j , for every i, j ∈ { 1 , . . ., n } .This property follows immediately from Eq. ( 17) .Pareto principle.The general interpretation of the Pareto principle (unanimity) in the social-choice theory is as follows: If all the experts agree upon a certain issue, then this agreement is reflected in the derived collective opinion.In our framework, it can be interpreted as the following requirement: If all the experts prefer the alternative x i over the alternative x j in their individual preference relations, then x i ranks higher than x j in the collective ranking.
More formally, if p (k )  i j ≥ p (k )   ji , for every k = 1 , . . ., m , then π c i ≥ π c j .This is a consequence of Eq. ( 17) and Eq. ( 15) : If p (k )  i j ≥ p (k )  ji , for every k = 1 , . . ., m , then, from Eq. ( 17) it follows that π (k )   i ≥ π (k )   j , for every k = 1 , . . ., m .From the last and Eq. ( 15) , it follows that π (c)   i ≥ π (c)   j .Note that the above two properties are based on qualitative comparisons between preferences where neither the magnitude of preferences and preference strengths nor the information on whether the dominated alternatives are weak or strong, is taken into account.Since these are crucial elements in the Rank Centrality method that enable deriving meaningful rankings, they justify the violation of the properties of internal consistency and Pareto optimality in the general case.This is clearly demonstrated by our example in Section 5.1 .It would be interesting, however, to examine the sensitivity of the latter properties to changes in the initial matrices.

Sensitivity Analysis (SA)
In this section we perform a theoretical sensitivity analysis of the Rank Centrality method.More specifically, we analyse the sensitivity of the output of the procedure described in Section 3.1 in face of small variations of the input, i.e. small variations in the values of the matrix parameters.The aim of this analysis is to understand how the centrality-based ranking of the alternatives is affected by small changes in the experts' opinions.As a basis for our analysis we use the derivatives of the output, and we largely apply results from Golub and Meyer ( Golub & Meyer, 1986 ).
Recall that the input of the RC ranking method is an n × n preference matrix P ( n ≥ 2 ) that is additive reciprocal, i.e. a matrix whose elements p i j are in the unit interval and obey the condition We consider an -perturbation of the matrix P resulting in a new matrix ˜ P ( ) , where the preferences associated with one particular pair of alternatives (x k , x l ) change as follows:5 with 0 < 1 , and consequently, from Eq. ( 5) , Let S and ˜ S ( ) be the centrality matrices of P and ˜ P ( ) respectively.Note that ˜ P (0) = P and ˜ S (0) = S. Whenever there is no confusion, we will omit the dependency on from the notation.According to Eq. ( 6) , the corresponding entries of ˜ S ( ) are then given by: As for the diagonal terms, the only terms affected by the perturbation in the preference matrix are ˜ s kk and ˜ s ll : ˜ All these variations can be written in a compact form for the centrality matrix as where t ji denotes the mean first passage time from x i to x j , i.e. the expected number of steps to reach state x j for the first time starting from state x i , in the unperturbed centrality matrix ˜ The proof of Theorem 2 is given in Appendix C .Note that Eq. ( 24) implies The perturbation of the preference associated with the pair (x k , x l ) affects every component ˜ π i in the stationary distribution ˜ π ( ) of the centrality matrix ˜ S ( ) .For sufficiently small these components can be written as6 Equation ( 25) can be interpreted as follows: For similar first passage times between alternatives, the perturbation of preferences in a single alternative-pair has an effect in the stationary distribution that is proportional to the amplitude of its components: dominant alternatives (higher rank) are more affected than other alternatives.One important consequence of Eq. ( 25) is that one can estimate first passage times for all possible transitions between alternatives.Indeed, for k = i and l = i , since t kk = t ll = 0 by definition of first passage time, Eq. ( 25) yields respectively ) Equation ( 25) can be used to estimate the global deviation of ˜ π ( ) from the "unperturbed" stationary distribution ˜ π (0) : where (28) This equation can be interpreted as follows.The global deviation of ˜ π ( ) from the unperturbed stationary distribution ˜ π (0) increases with a "weighted" second moment of the components of the unperturbed stationary distribution.The weights for the component π i are given by the difference between the first passage times from the corresponding alternative x i to the "perturbed altrenatives" x k and x l .This observation has two main consequences.First, the result in Eq. ( 27) uncovers another intuitive consequence: the components of the stationary distribution which remain unchanged by perturbing the preference for a pair of alternatives (x k , x l ) are those for which the mean first passage times from the corresponding alternative x i to each alternative in the perturbed pair is equal, i.e. t ki = t li .In particular, perturbing the preference for pairs of alternatives (x k , x l ) , which are evenly chosen in front of all other alternatives x i , i.e. t ki = t li for all i = k and i = l, will have no impact on the global preference strengths of the alternatives.Their impact is reduced to local changes of the amplitudes of each alternative in the perturbed pair, namely ˜ Second, since deviations between perturbed and unperturbed stationary distribution are easy to measure in practice, Eq. ( 25) and ( 27) provide new insight for establishing a framework to assess, at least at a qualitative level, the impact of local perturbations in the global dynamics towards consensus.Namely, the results of this section have some implications when it comes to reaching consensus worth pursuing as future work.In fact, we can formalize the consensus problem as a gradient descent optimization where experts need to do small adjustment to their preference matrices.The sensitivity of Markov chains to their transition probabilities can be used for computing the gradient in order to make the individual stationary distribution of each expert move towards collective stationary distribution, by only adjusting the corresponding preference matrices.

Testing GDM with FPR: numerical experiments
In this section we provide a concrete example of a GDM with FPR using our method.We first compare the results of using our method against the results of using the popular Netflow method.Then, we perform an experimental sensitivity analyses to explore what happens to the resulting ranking vector if we vary the entries of the initial preference matrices 1. by several values of > 0 applied to the same pair of alternatives; and 2. by a small fixed > 0 applied to different pairs of alternatives in the matrix.
We start by a numerical example highlighting the similarities and differences between our algorithm and an algorithm that uses the standard Netflow (NF) method.

A numerical example
We consider a first scenario with two experts ( m = 2 ) and four alternatives ( n = 4 ) 7 The experts' opinions are expressed in the fol-lowing preference matrices:

(29b)
The entries of the preferences matrices P 1 and P 2 differ by 5% or more, in a scale from 0 to 1.They can be provided by experts directly as numerical values or obtained, for example, by qualitative preference modelling. 8 Assuming, for the sake of illustration, that these preferences represent relationships between pairs of four different sport teams in two different leagues, the relationships in both leagues share two common features: (i) the third team is moderately better than the first and the second team, while the fourth team is much better than the first and the second team, and (ii) the third team seems, against intuition, to be better than the fourth team.We observe what happens with the ranking vectors derived from these two matrices when we apply our method.
First of all, we normalize the preference matrices into stochastic matrices applying the Eq. ( 6) from Section 3.1 .This produces the following centrality matrices: Then we apply two direct methods to derive a group ranking of the alternatives: • (NF+WA) : exploitation is performed by computing a per-expert ranking using Netflow, and then aggregation is performed by computing the final collective ranking using WA with uniform weights.• (RC+WA) : exploitation is performed by computing a per-expert ranking using Rank Centrality, and then aggregation is performed by computing the final collective ranking using WA with uniform weights.
The output weighted rankings produced by these methods are the following: output NF+WA = −0 .383 −0 .517 0 .36 0 .54 output RC+WA = 0 .087 0 .069 0 .388 0 .456 , 8 One can, for example, use linguistic preference modelling ( Herrera, Alonso, Chiclana, & Herrera-Viedma, 2009 ), and even give each individual expert the possibility to use different preference domains to express their respective preferences.This issue is studied in ( Delgado, Herrera, Herrera-Viedma, & Martinez, 1998 ) where it is shown that in such heterogeneous decision contexts, it is possible to achieve a solution by first making the preferences uniform by converting them into FPRs over [0,1] by means of transformation functions.
corresponding to the following qualitative rankings: rank NF+WA = 4 3 1 2 rank RC+WA = 4 3 1 2 , where the i-th element is the index of the alternative with the ith best preference strength (in this case, the 4-th alternative is the most preferred, the 3-rd alternative is the second most preferred, and so on).The results are identical, and they rank the fourth alternative at the top, which captures the intuition in the football teams example.Note that in the above example, we have violation of both Internal Consistency and Pareto Optimality at the alternatives 3 and 4, but we still obtain a meaningful ranking.
We now consider a second scenario, again with two experts ( m = 2 ) and four alternatives ( n = 4 ), but with the following preference matrices: Notice that this new set of preference matrices preserves the relationships of dominance described in i) and ii) above; however the degrees by which the third alternative dominates the first and the second one have been increased, although they remain lower than the corresponding ones of the fourth alternative.
After normalization according to Eq. ( 6) from Section 3.1 we obtain the following stochastic centrality matrices: This time, the results differ in the ordering of the third and the forth alternative.The difference in the results between the first and the second scenario shows that the approach based on rank centrality has a distinctive sensitivity to the magnitude of the relations of preference among the alternatives.In particular, the approach based on RC+WA may rank higher alternatives that marginally dominate weak alternatives but also marginally dominate strong alternatives.This is consistent with the interpretation of rank centrality provided in ( Negahban et al., 2012;2016 ).29a) and (30) .

Experimental sensitivity analysis
In this section we perform experimental sensitivity analysis, investigating how the result of the group decision making changes under a perturbation of amplitude in the initial preference matrices.We restrict to the case when only one of the initial preference matrices changes in only one pair of alternatives, and we observe how various degrees of change affect the result.We consider again the first scenario we presented in the previous simulation, but this time we instantiate a parametric version of the matrix P 1 for the first expert: with the parameter assuming values in the interval [0,0.3) in order to satisfy the requirement that, for every entry, 0 < p i j < 1 .The parameter allows us to increase the margin of the preference of the third alternative over the first, and, consequently, to narrow the gap between the third and the fourth alternative.The preference matrix of the second expert is taken to be the same fixed matrix P 2 used in the previous section.We then apply the two direct methods we considered before (NF+WA and RC+WA) to the matrices ˜ P 1 ( ) and P 2 , while changing the value of the parameter .
Fig. 2 shows the variation in the output of the two methods as a function of the parameter .The results show that the two direct methods respond differently to similar changes.As the parameter increases, the gap between the third and the fourth alternative narrows more significantly when using RC+WA instead of NF+WA.This change also leads to a decrease in the value of the first alternative that is more marked for the RC+WA method; indeed, for values of around 0.20 the first alternative becomes less preferable than the second one; for the range of that we considered, we do not notice a similar change in the ordering of the first and the second alternative in the NF+WA method.
When applying RC+WA in a scenario with many alternatives where we only perturb one pair of them, it is not always obvious how this perturbation will reflect on the resulting preference strengths of the alternatives: which of them will have their pref- erence strengths increased or decreased and for how much.As we will see, our theoretical perturbation analysis given by Eq. ( 24) can predict and interpret the changes based on the mean first passage times and the magnitudes of the preference strengths.The true mean first passage times are given in Table 1 .In our case, the pair that is changed is (x k , x l ) = (3 , 1) .As per our theoretical results, the preference strength of the alternative 3 for expert 1 computed by RC will increase proportionally to t 13 = 12 .7326 moderated by the strength π (1)   3 of the alternative itself. 9Similarly the preference strength of the alternative 1 for expert 1 computed by Deviation (%) by the strength the alternative itself π (1) 2 .Since π (1)   2 is small, the magnitude of the latter changes is small too as can be seen in the aggregated weighted ranking given in Fig. 2 .
To end the sensitivity analysis for this particular numerical experiment, we fix the perturbation amplitude at a small value, namely = 0 .001 , and apply it to each pair of alternatives (x k , x l ) in the matrix P 1 .The results are shown in Fig. 3 .While the larger deviations are typically observed for the components of the stationary distribution associated with the perturbed pair, namely π k and π l , one also observes significant changes in the other components.For example, when perturbing the pair (4,2) one observes also a significant change in the amplitude of the first component, and when perturbing the pair (3,4), all components are significantly affected.Note that, in all cases the deviations for (x k , x l ) and (x l , x k ) are symmetric, as expected from Eq. ( 25) .
Finally, in Fig. 4 we show the difference of the first passage times, t li − t ki , for each perturbated pair of alternatives (x k , x l ) considered in Fig. 3 , as computed directly from Eq. ( 25) .Similarly, one also observes cases where the perturbation of one particular pair of alternatives induces a change in the transitions from another al-  ternative to the alternatives in the perturbed pair.Moreover, applying Eq. ( 26) we estimate the first passage times for all pairs of alternatives.Notice that, repeating the numerical experiment interchanging the role of k and l enables to make two estimates for each first passage time.As shown in Tab. 1 , in all the cases both estimates are close to each other, showing the ability of our procedure for estimating this dynamical property of consensus processes.

Conclusions and future work
In this paper we have presented a direct approach to aggregating fuzzy preference relations proposing a GDM method based on rank centrality.The method has the advantage of providing a natural interpretation of the preference degrees as transition probabilities in a Markov chain and obtaining the corresponding weighted rankings by well-established computational methods.Moreover, as we show with our numerical examples, our approach shows more sensitivity to small variations in the preference values compared to other similar approaches.
The natural next step is to design an experiment to test our framework and compare it with other GDM with FPR frameworks in the literature.We have implemented an online platform for collecting data during an iterative process towards consensus, which will enable to investigate the distances in the opinion space, either to the collective opinion, or pairwise distances between experts' opinions.In this way we can test our framework for modelling processes towards reaching a consensus, and examine which initial preference matrices lead to a consensual opinion.Finally, such experimental setup will also enable to investigate the time interval needed for achieving consensus in various real-life applications, and determine which framework enables the fastest converging iterative processes.Moreover, investigating the inconsistency ( Kou, Ergu, & Shang, 2014;Kou & Lin, 2014;Lin, Kou, Peng, & Alsaadi, 2020 ) directly from a centrality matrix is a future research direction worth investigating.
We hope that the current work can fuel more research interest in bridging the gap between the GDM community and researchers in probability theory and its applications.
any state i , there is a non-zero probability of returning to state i in k steps, for k = { 1 , 2 , 3 , . . .} .The greatest common divisor of this number of steps is then 1. Thus every state i is aperiodic, and the Markov chain with its associated centrality matrix is aperiodic.

Appendix B. Proof of Theorem 1
To prove Theorem 1 , let us consider a BTL model with parameter vector π such that p i j = π i π i + π j .Let M i j be the number of sample comparisons between x i and x j and m i j the samples in which alternative x i wins over x j .Hence, there is a real number, M, such that M i j = m i j + m ji ≡ M(π i + π j ) .(B.1) where .denotes the operator that rounds to nearest integer.In what follows we will only need the limit of large values of M i j (i.e.M i j → ∞ ) for which M → ∞ and M(π i + π j ) → M(π i + π j ) .
Under the assumption of independent and identically distributed samples, the likelihood function of the BTL model is given by: where we disregard the case π i = 0 , which only occurs in the "pathological" case where for some j = i , p i j = 0 .We take the loglikelihood function of (B.
We determine the MLE π * as the value of π at which the partial derivatives are zero, that is: By applying the law of large numbers, p i j can also be defined as: * i be the i-th column of A * .Then, we can apply Theorem 3.2 from Golub and Meyer ( Golub & Meyer, 1986 ) to study the sensitivity of the ranking to changes of , namely: where δ i j is the Kronecker-delta, δ i j = 1 if i = j and zero otherwise.

Fig. 2 .
Fig. 2.Output of the two GDM strategies (NF+WA and RC+WA) when applied to ˜ P 1 ( ) and P 2 as a function of the parameter .See Eq. (29a) and (30) .
RC will decrease proportionally to t 31 = 4 .8463 moderated by the strength π(1)   1 of the alternative itself.But what about the rest of alternatives?The changes in the preference strength of the alternative 4 for expert 1 is proportional to t 34 − t 14 = 5 .3184 − 15 .1259 = −9 .8075 moderated by π(1)

Fig. 3 .
Fig. 3. Deviation of each component of the perturbed stationary distribution ˜ π ( ) in percentage of the respective component of the unperturbed stationary distribution π (0) , when the preference matrix is changed for one single pair of alternatives (x k , x l ) .In all cases = 0 .001 .See Eq. (25) .

Fig. 4 .
Fig. 4. For each perturbation of the preference matrix shown in Fig. 3 one plots the difference of first passage times, t li − t ki , from each alternative x i to each alternative of the perturbed pair (x k , x l ) .The values were obtained by solving Eq. (25) for t li − t ki .
Maximum Likelihood Estimate (MLE) of the BTL model we take the partial derivative of the likelihood function with respect to the parameters of the model: Please note that we used the argument that π * i + π * j ≈ π i + π j by virtue of the MLE.Similarly,p ji = lim M→∞ m ji M(π * i + π * j ) Hence in the limit of M → ∞ m i j = p i j M(π * i + π * j )Let A * be the group inverse of A = I − S( Golub & Meyer, 1986 ), i.e.A * is the unique matrix 10 satisfying the three equations AA * A = A , A * AA * = A * and A * A = AA * .Let a * i j denote the entries of A * , and let A * e i = (0 , . . ., 0 , 1 , 0 , . . ., 0) T be the unit vector with entries δ ik for k = 1 , . . ., n .Then, e k − e k e l − e l e l + e l e k ] ,where denotes the tensor product between unit vectors.Substituting this expression in Eq. (C.1) yields e k e k ) rs − (e k e l ) rs −(e l e l ) rs + (e l e k ) rs a * si δ kr δ ks − δ kr δ ls −δ lr δ ls + δ lr δ ks ) a * 50 0 .60 0 .20 0 .10 0 .40 0 .50 0 .15 0 .05 0 .80 0 .85 0 .50 0 .55 0 .90 0 .95 0 .45 0 .50

Table 1
Estimates of all first passages times ˆ t i j according to Eq. (26) for matrix S 1 , compared with the values computed directly from the respective Markov chain simulation.