Keywords: Rater model, Kendall’s tau, minimax optimality, ranking, missing data, U-empirical process

1 sec:Introduction↩︎

Rank aggregation, the process of synthesizing multiple orderings or ratings into a single global ranking, is a fundamental problem in information retrieval and decision science. Models and methods for rank aggregation have broad applications in practice, ranging from ranking players or teams in sports to aggregating customers’ preferences in marketing surveys. Traditionally, this task is typically solved using parametric models, including the famous Thurstonian and Bradley-Terry (BT) models [1], [2] for pairwise comparison data, and the Plackett–Luce (PL) model [3], [4] for multi-way comparison data. Under parametric model assumptions, these models transform a rank aggregation problem, which is discrete in nature, into a parameter estimation problem that can be solved via continuous optimization. The statistical properties of the BT model have been established under both dense [5] and sparse settings [6]–[8], as well as sequential settings where comparison pairs are adaptively chosen [9]. More recently, theoretical results have also been established for the PL model [10]–[12].

In many real-world applications, the parametric assumptions of the Thurstonian, BT, and PL models can be overly simplified, resulting in biased rankings. Moving beyond restrictive parametric assumptions, nonparametric models and methods have been developed for pairwise comparison data, including rank aggregation based on Borda counting [13], estimation of pairwise comparison probabilities under nonparametric BT models based on universal singular value thresholding [14] (USVT) and monotone constraints [15], and noisy sorting [16], [17]. To our knowledge, however, no nonparametric models or methods have been developed under general multi-way comparison settings.

This paper proposes a MAximum Score esTimator for aggEgating Ranks (MASTER), a model-free approach for rank aggregation based on multi-way comparison data derived from raters’ scores. More specifically, we consider a system of \(n\) items and \(p\) raters, where each rater ranks a subset of the items, with the objective of recovering the latent global ranking. In the limiting case where each rater evaluates only two items, our framework reduces to the classical pairwise comparison setting. The MASTER analyzes observed multi-way comparison data under a highly flexible setting, in which the data distribution is allowed to be highly heterogeneous across raters. This setting covers rank aggregation based on pairwise comparisons under the Weak Stochastic Transitivity (WST) assumption, as studied in [17]. We establish the theoretical foundations of the proposed estimator by proving its consistency, showing that the proportion of discordant pairs (Kendall’s tau) converges to zero in probability as the number of raters diverges. Furthermore, we derive error bounds for a performance metric based on Kendall’s tau, covering both dense and sparse regimes where the number of ratings per rater can be either a fixed constant or diverge with the total number of items. Under certain asymptotic regimes, these error bounds imply that the proposed estimator is nearly minimax optimal. Simulation studies in both pairwise and multi-way comparison settings show that MASTER substantially improves over existing methods under heterogeneous preference distributions and non-uniform sampling, while remaining competitive when the BT or PL model is correctly specified. We further apply MASTER to a professional tennis dataset based on pairwise match outcomes and a sushi preference dataset based on multi-way partial rankings. The real-data analyses show that MASTER yields interpretable aggregate rankings and reveals meaningful departures from parametric rankings in the presence of preference heterogeneity.

The proposed method is closely related to the noisy sorting framework investigated by [17], yet there are fundamental differences. While [17] addressed rank aggregation for pairwise comparisons under the WST and uniform-sampling assumptions, our work considers a significantly more complex setting involving multi-way comparisons and rater heterogeneity. In our framework, the model space is considerably larger and more complex; notably, we allow a general sampling scheme, moving beyond the uniform-sampling requirements of prior work. Consequently, the theoretical machinery employed in [17] to establish minimax optimality is insufficient for analyzing MASTER. In the special case of pairwise comparisons, our method coincides with an estimator used in [17] for minimax analysis. By using more general proof techniques, our theoretical results not only recover the minimax rates established therein but also provide theoretical guarantees to broader, non-uniform, and multi-way comparison regimes.

In summary, the contributions of this paper are threefold. First, we establish a general framework for rank aggregation under significantly weaker assumptions and in more general settings than those in the existing literature. Our model-free approach encompasses traditional parametric models, such as the Bradley-Terry (BT) and Plackett-Luce (PL) models, as well as nonparametric BT variants based on Strong Stochastic Transitivity (SST). Unlike previous methods, such as noisy sorting, our framework does not rely on uniform sampling, offering greater flexibility in accounting for heterogeneity in rater behavior. Second, we provide rigorous error bounds that illuminate the information-theoretic fundamental limits and convergence rates of this rank aggregation problem, which, to the best of our knowledge, is the first model-free theory for rank aggregation from heterogeneous multi-way comparisons. The key technical contribution is a uniform concentration analysis for a nonstandard weighted U-empirical process indexed by permutations. The resulting estimator is nearly minimax optimal in the bounded-comparison-size regime and remains consistent under substantially more general sampling and preference heterogeneity than existing approaches. Finally, we develop a practical estimation algorithm for the MASTER objective. Although maximizing the MASTER objective function is inherently NP-hard due to its combinatorial structure, we demonstrate through extensive simulations and real-data applications that our algorithm is empirically robust and yields high-quality rankings in practice.

The following notation is used in the rest. For an integer \(n\), denote \([n] = \{1,...,n\}\) and define \(S_n\) as the permutation set of \([n]\). We write \(1\{\cdot\}\) for the indicator function and \(|A|\) for the cardinality of a finite set \(A\). Further define \(S_n^+ = \{\boldsymbol{\pi}= (\pi_1,...,\pi_n): \pi_j = 1 + |\{i\in[n]: y_i<y_j\}|~~\text{for some } \mathbf{y}\in\mathbb{R}^n\}\) as all rank vector representations of weak orderings that allow ties. For any vectors \(\mathbf{y}= (y_1,...,y_n)\in\mathbb{R}^n\), define the rank vector for \(\mathbf{y}\) as \(\text{rank}(\mathbf{y}) = (\pi_1,...,\pi_n)\in S_n^+\), where \(\pi_j = 1 + |\{i\in[n]: y_i<y_j\}|\) represents the rank of \(y_j\) in an increasing order for \(j\in[n]\). For example, for \(\mathbf{y}= (0.1, 0.2, 0.2, 0)\), we have \(\text{rank}(\mathbf{y}) = (2, 3, 3, 1)\). For two rankings \(\boldsymbol{\pi},\boldsymbol{\pi}'\in S_n\), define Kendall’s tau distance as \[\tau(\boldsymbol{\pi},\boldsymbol{\pi}') = \left| \{(i,j):i<j,\;(\pi_i-\pi_j)(\pi'_i-\pi'_j)<0\} \right|.\] For positive sequences \(\{a_n\}\) and \(\{b_n\}\), we denote \(a_n \lesssim b_n\) or, equivalently, \(b_n \gtrsim a_n\), if there exists a constant \(c>0\) that \(a_n \leq c b_n\) for all \(n\). We also write \(a_n \asymp b_n\) if \(a_n \lesssim b_n\) and \(a_n \gtrsim b_n\) hold.

The remainder of the paper is organized as follows. Section 2 introduces the MASTER framework and its underlying assumptions, followed by the derivation of upper and lower convergence bounds and an algorithm for its computation. In Section 3, we evaluate the performance of the proposed method against several state-of-the-art competitors through extensive simulations and two real-world applications: ranking professional tennis players using pairwise comparisons and aggregating sushi preferences from multi-way comparison data. Finally, Section 4 concludes the paper with a summary of our contributions and a discussion of potential future directions.

2 Proposed Method↩︎

2.1 Proposed MASTER method↩︎

Consider \(n\) items ranked by \(p\) raters, where each rater \(k\in[p]\) ranks a subset of \(n\) items, denoted as \(\mathcal{C}_k\subseteq[n]\) with \(L_k = |\mathcal{C}_k|\). Let \(y_{ik}\) be the score on item \(i\) given by rater \(k\), allowed to be of different data types, such as continuous, ordinal, or count. A higher score on an item than another indicates that the rater ranks it higher. Our goal is to derive a global ranking based on the observed data \(y_{ik}\), \(k=1, ..., p\), \(i \in \mathcal{C}_k\). We define MASTER as follows: \[\label{eq:mse} \widehat\boldsymbol{\pi}= (\hat{\pi}_1, ..., \hat{\pi}_n)^\top \in \mathop{\mathrm{arg\,max}}_{\boldsymbol{\pi}\in S_n} L(\boldsymbol{\pi}),\tag{1}\] where \(\boldsymbol{\pi}= (\pi_1, ..., \pi_n) \in S_n\) is a permutation of \((1, ..., n)\) and \[L(\boldsymbol{\pi}) = \sum_{k=1}^p \sum_{i<j\in \mathcal{C}_k} (1_{\{y_{ik} < y_{jk}\}}-1_{\{y_{jk}< y_{ik}\}}) 1_{\{\pi_i < \pi_j\}}\] is our objective function. The output \(\widehat \boldsymbol{\pi}\) gives the estimated global ranking for the \(n\) items, where \(\pi_i\) gives the ranking of item \(i\) among the \(n\) items. Specifically, \(\hat{\pi}_i = n\) means item \(i\) is ranked the highest among all the items, while \(\hat{\pi}_i = 1\) means it is ranked the lowest. We provide a few remarks about this estimator.

2.2 Theoretical properties of MASTER↩︎

In what follows, we first introduce the assumptions underlying MASTER and then establish its theoretical properties.

Before presenting our main theoretical results, we provide some intuition about the proposed estimator. Consider the expectation of our objective function \(l(\boldsymbol{\pi})\), \[\begin{align} \mathbb{E}L(\boldsymbol{\pi}) =& \sum_{k=1}^p\sum_{i<j\in\mathcal{C}_k}\left(\Pr(y_{ik}<y_{jk}) - \Pr(y_{jk}< y_{ik})\right)1_{\{\pi_i < \pi_j\}}. \end{align}\] According to Assumption 1, the multiplicative factor in each term of \(\mathbb{E}L(\boldsymbol{\pi})\), \(\Pr(y_{ik}<y_{jk}) - \Pr(y_{jk}< y_{ik})\), is positive when \(\pi_i^* < \pi_j^*\) and negative otherwise. Thus, it is obvious that the expected objective function is maximized when \(\boldsymbol{\pi}= \boldsymbol{\pi}^*\). In our theoretical analysis, we establish the concentration of \(L(\boldsymbol{\pi})\) to \(\mathbb{E}L(\boldsymbol{\pi})\), based on which we show the convergence of \(\hat{\pi}\) towards \(\pi^*\).

For \(k\in[p]\), define \(\mathcal{P}_k = \{(i,j):i<j\in\mathcal{C}_k\}\) as the set of item pairs compared by rater \(k\). For any \(\boldsymbol{\pi},\boldsymbol{\pi}'\in S_n\), define the discordant pair set between them as \[B(\boldsymbol{\pi},\boldsymbol{\pi}') = \{(i,j):~i<j ~\text{and}~ (\pi_i-\pi_j)(\pi_i'-\pi_j')<0\}.\] For \(1\leq i<j\leq n\) and \(k\in[p]\), define the probability gap, which quantifies the difficulty for rater \(k\) to distinguish between items \(i\) and \(j\), as \(q_{ijk} = |\Pr(y_{ik}<y_{jk}) - \Pr(y_{jk}< y_{ik})|.\) We assess the accuracy of \(\hat{\pi}\) by a modified Kendall’s tau distance \[\tau_{m}(\widehat\boldsymbol{\pi},\boldsymbol{\pi}^*) = \frac{(\sum_{k=1}^p\sum_{(i,j)\in B(\widehat\boldsymbol{\pi},\boldsymbol{\pi}^*)\cap\mathcal{P}_k}q_{ijk})^2}{\sum_{k=1}^p|B(\widehat\boldsymbol{\pi},\boldsymbol{\pi}^*)\cap \mathcal{P}_k|}.\] The metric \(\tau_{m}\) is the intrinsic loss induced by the fixed sampling design and heterogeneous pairwise signal strengths. The Proposition II.1 below shows that, under random sampling, it controls the standard Kendall’s tau distance up to the sampling intensity \(\bar N\) and the minimum signal strength \(q_{\min}\).

In Theorem 1 below, we establish an upper bound for the modified Kendall’s tau metric.

We relate the modified Kendall’s tau distance to the standard Kendall’s tau distance, under a flexible non-uniform random sampling design as specified by the following Assumption 2.

Define \(N_{ij} = \sum_{k=1}^p 1_{\{\{i,j\}\subseteq \mathcal{C}_k\}}\) as the number of raters that rate both items \(i\) and \(j\), and \(\bar{N} = \sum_{k=1}^p L_k(L_k-1)/(n(n-1))\) as the average number of total comparisons for different item pairs. Let \(q_{\min} = \min_{k, i<j} q_{ijk} > 0\). The following Proposition 1 establishes a non-asymptotic upper bound for the Kendall’s tau distance under Assumption 2.

We now specify certain scenarios under which the proportion of discordant pairs converges to zero.

As mentioned in Section 2.1, the discrete optimization 1 contains \(n\)! elements in its feasible space, making it a challenge to find the global maximum. However, we point out that we do not have to find the global maximum; instead, an approximate solution \(\widetilde{\boldsymbol{\pi}}\) satisfying \(L(\widetilde{\boldsymbol{\pi}})\geq L(\boldsymbol{\pi}^*)\) with high probability is sufficient. We formalize this in Proposition 3 below, the proof of which is given in the Appendix.

In the following, we turn to establish the minimax lower bound for the rank estimator with respect to the modified Kendall’s tau distance. To this end, suppose that for each rater \(k\in[p]\), we only observe the ranking of his ratings on items in \(\mathcal{C}_k\), e.g. we observe \(\mathcal{D}= \{(\mathcal{C}_k,\boldsymbol{\pi}_{\text{obs},k})\}_{k=1}^p\), where \(\boldsymbol{\pi}_{\text{obs},k} = \text{rank}(\mathbf{y}_{\mathcal{C}_k,k})\in S_{L_k}^+\). Consider estimators based on \(\mathcal{D}\), e.g. \(\check\boldsymbol{\pi}= T(\mathcal{D})\). We establish the minimax lower bound under the random sampling design of Assumption 2. Define the parameter space for the distribution of \(\mathcal{D}\) as \(\mathcal{G}\) which consists of all joint distributions for \(\mathcal{D}\) satisfying Assumptions 1 and 2. For any \(G\in \mathcal{G}\), due to the data generating process for \(\mathcal{D}\), there exists a unique ranking, denoted as \(\boldsymbol{\pi}(G)\in S_n\), such that \(\Pr(y_{ik}<y_{jk}) > \Pr(y_{jk}< y_{ik})\) for all \(i,j\in \mathcal{C}_k\) and \(k\in[p]\). In the following, we consider to derive a lower bound for \(\min_{\check\boldsymbol{\pi}= T(\mathcal{D}) }\max_{G\in\mathcal{G}}\tau_{m}(\check\boldsymbol{\pi},\boldsymbol{\pi}(G))\).

2.3 Computation↩︎

The objective function for MASTER in 1 involves an NP-hard discrete optimization problem. To find a good solution to 1 , we develop a greedy local search algorithm and combined it with an initial ranking obtained via smoothing. Our initialization is motivated by a smoothed maximum score estimator for binary regression [19]. Let \(w_{ij}\) and \(l_{ij}\) be the numbers of raters that rank item \(i\) higher and lower than item \(j\), respectively. Then the objective function for MASTER can be written as \(L(\boldsymbol{\pi}) = \sum_{i<j} (l_{ij} - w_{ij})1_{\{\pi_i<\pi_j\}}\). We consider the following smooth optimization \[\label{eq:surrogate} {\boldsymbol{\beta}}_{\text{init}} = \mathop{\mathrm{arg\,max}}_{{\boldsymbol{\beta}}\in\mathbb{R}^n} \sum_{i<j}(l_{ij} - w_{ij})\frac{e^{(\beta_j - \beta_i)/h}}{1 + e^{(\beta_j - \beta_i)/h}},\quad\text{subject to}\quad \sum_{i=1}^n\beta_i=0, \quad \frac{1}{n}\sum_{i=1}^n\beta_i^2=1.\tag{2}\] where we replace the indicator function \(1_{\{\pi_j-\pi_i > 0\}}\) in \(L(\boldsymbol{\pi})\) by the function \(\phi_h(x)=(1+\exp(-x/h))^{-1}\) for some small \(h>0\) to facilitate smooth optimization. Note that as \(h\to0\), \(\phi_h(x)\to 1\) if \(x>0\), and \(\phi_h(x)\to 0\) if \(x<0\). Therefore we use the objective function in 2 to approximate \(L(\cdot)\). The constraints fix the scale of \({\boldsymbol{\beta}}\), so that the smoothing parameter \(h\) cannot be absorbed by rescaling. We solve this problem sequentially for \(h=1, 0.5, 0.25, 0.1,0.05,0.01\) using the solution from the previous value of \(h\) as the warm start for the next one. After each gradient ascent step, \({\boldsymbol{\beta}}\) is re-centered and re-scaled by \[{\boldsymbol{\beta}} \leftarrow \frac{\sqrt n({\boldsymbol{\beta}}-\bar\beta{\boldsymbol{1}})}{\|{\boldsymbol{\beta}}-\bar\beta{\boldsymbol{1}}\|_2}, \qquad \bar\beta=\frac{1}{n}\sum_{i=1}^n\beta_i.\] Our initial ranking \(\boldsymbol{\pi}_{\text{init}}\) is then obtained by \(\boldsymbol{\pi}_{\text{init}} = \text{rank}({\boldsymbol{\beta}}_{\text{init}})\). After getting \(\boldsymbol{\pi}_{\text{init}}\), we perform a consecutive \(K\)-tuple searching algorithm to obtain the final ranking, as described in Algorithm 1 below.

The feasible space of this \(n\)-dimensional discrete optimization problem 1 consists of \(n!\) elements. There is no guarantee that Algorithm 1 finds the global optimal solution. However, as shown in Proposition 3, we only need \(\boldsymbol{\pi}_{\text{curr}}\) to satisfy \(L(\boldsymbol{\pi}_{\text{curr}})>L(\boldsymbol{\pi}^*)\) with high probability. According to our simulation study with \(K=3\), the estimates given by Algorithm 1 satisfy this condition in more than \(95\%\) of experiments.

3 Numerical Experiments↩︎

3.1 Simulation study↩︎

Pairwise comparison. We first consider a pairwise comparison setting, which corresponds to the special case \(L_k=2\) in our general framework. For each simulation setting, we generate \(p=n^2\) raters. To allow non-uniform comparison frequencies, we first generate pair-specific sampling weights \(\xi_{ij}\sim {\rm Uniform}(0.3,0.5)\) for \(1\le i<j\le n\). Then each rater independently selects one item pair according to \[\Pr(\mathcal{C}_k=\{i,j\}\mid \{\xi_{ab}\}) = \frac{\xi_{ij}}{\sum_{1\le a<b\le n}\xi_{ab}}, \qquad k=1,\ldots,p.\] This design is a pairwise instance of Assumption 2, since the sampling probability of each pair is of order \({n\choose 2}^{-1}\). Let \(N_{ij}=\sum_{k=1}^p 1_{\{\mathcal{C}_k=\{i,j\}\}}\) denote the realized number of comparisons between items \(i\) and \(j\). Given the selected pair \(\mathcal{C}_k=\{i,j\}\) with \(i<j\), we generate the comparison outcome \[Z_k=1_{\{\text{rater }k\text{ ranks item }j\text{ above item }i\}}\quad \text{from}\quad Z_k\sim {\rm Bernoulli}(p_{jik}),\] where the rater-specific comparison probabilities \(p_{jik}\) are generated according to one of the following scenarios. Throughout the three scenarios, we set the true global ranking as \(\pi_i^*=i\).

1. Random heterogeneity. For each observed pair \(\mathcal{C}_k=\{i,j\}\) with \(i<j\), generate \(p_{jik}\sim {\rm Uniform}(0.55,0.95)\).

2. Group-structured heterogeneity. Divide the items into two groups, \(\{1,\ldots,n/2\}\) and \(\{n/2+1,\ldots,n\}\). For each observed pair \(C_k=\{i,j\}\) with \(i<j\), generate \[p_{jik}\sim \begin{cases} {\rm Uniform}(0.75,0.85), & \text{if } i,j \text{ are in the same group},\\ {\rm Uniform}(0.65,0.75), & \text{otherwise}. \end{cases}\]

3. BT model. Generate latent scores \(x_1,\ldots,x_n\) independently from \(N(0,2)\) and sort them in increasing order as \(x_{(1)}<\cdots<x_{(n)}\). Set \[p_{jik} = \frac{1}{1+\exp[-(x_{(j)}-x_{(i)})]}, \qquad i<j.\]

It is not difficult to see that all three scenarios satisfy the pairwise version of the EWST assumption. In the first two scenarios, for every observed pair \(C_k=\{i,j\}\) with \(i<j\), the rater-specific comparison probability satisfies \(p_{jik}>1/2\), and hence item \(j\) is more likely to be ranked above item \(i\) by rater \(k\). However, these two scenarios do not satisfy SST in general. In the first scenario, this is due to the randomness and heterogeneity in generating the rater-specific pairwise comparison probabilities. In the second scenario, the violation is more systematic: comparison probabilities tend to be closer to \(1/2\) when the rank difference between two items is larger, i.e., when they are not in the same group, which is opposite to the monotonicity required by SST. In contrast, the third scenario follows a BT model. The comparison probability \(p_{jik}\) is monotone increasing in \(x_{(j)}-x_{(i)}\). Therefore, the third scenario satisfies SST and reduces to the standard homogeneous BT model.

We compare the proposed MASTER with competitors designed for pairwise comparison data including the maximum likelihood estimator based on the BT model in [6], the Borda counting algorithm proposed in [13], and universal singular value thresholding (USVT) [14]. For the USVT method, we apply the algorithm in [14] to the standardized skew-symmetric matrix \(\mathbf{X}= (x_{ij})_{n\times n}\) with \[x_{ij} = \left\{ \begin{align} &\frac{2w_{ij}}{w_{ij}+l_{ij}} - 1,\quad &&\text{if~~} w_{ij} + l_{ij} > 0,\\ &0,\quad&& \text{if~~} w_{ij} + l_{ij} = 0, \end{align} \right.\] to estimate the pairwise comparison probabilities, and then rank the players by the row sums of the estimated pairwise comparison probability matrix.

By combining the choices of item size with the three scenarios for comparison probabilities, we obtain nine simulation settings. For each setting, we generate 100 independent datasets. For each dataset, the methods mentioned previously are applied. For the MASTER, we run Algorithm 1 with search size \(K =3\). The accuracy of the aggregated ranking is measured by the standard Kendall’s tau ranking error, which is calculated as the proportion of discordant pairs, i.e., \(2\tau(\hat{\pi},\pi^*) / (n(n-1))\). The results are summarized in Table 1.

Table 1 shows that the proposed method substantially outperforms the competitors, achieving an error reduction by factors ranging from 2.8 to 16.8 compared to the second-best method. The advantage is more obvious in scenario 2, where the SST assumption is violated more severely. These results are due to the fact that all competitors assume SST. When the SST assumption is violated, the aggregated ranks from these methods are biased, resulting in higher errors compared to the proposed method. On the other hand, in scenario 3, where SST holds, all three competitors perform slightly better than the MASTER, and unsurprisingly, the BT model, under which the data are generated, performs the best. However, we should note that the MASTER is only slightly inferior to its competitors, due to the efficiency loss caused by having a larger parameter space. We should also note that the MASTER remains consistent in scenario 3, as SST implies WST.

Multi-way comparison. We next consider a multi-way comparison setting. For each simulation setting, we generate \(p=n^2\) raters. Let \(m=n/2\), and divide the items into two groups \(G_1=\{1,\ldots,m\}\) and \(G_2=\{m+1,\ldots,n\}\). We generate item-specific sampling weights \(\xi_i\sim\text{Uniform}(0.5,1.5)\), \(i\in[n]\). The number of items evaluated by rater \(k\) is generated independently from \(L_k\sim\text{Uniform}\{3,4,5\}\). Conditional on \(L_k\), with probability \(1-\eta\), where \(\eta=0.2\), we first draw \(g_k\sim\text{Uniform}\{1,2\}\) and then sample \(\mathcal{C}_k\subset G_{g_k}\) with \(|\mathcal{C}_k|=L_k\) according to \[\Pr(\mathcal{C}_k=\mathcal{C}\mid g_k,L_k,\{\xi_i\}_{i=1}^n) = \frac{\prod_{i\in\mathcal{C}}\xi_i}{\sum_{\mathcal{C}'\subset G_{g_k}:|\mathcal{C}'|=L_k}\prod_{i\in\mathcal{C}'}\xi_i}.\] With probability \(\eta\), we sample \(\mathcal{C}_k\subset[n]\) with \(|\mathcal{C}_k|=L_k\) according to \[\Pr(\mathcal{C}_k=\mathcal{C}\mid L_k,\{\xi_i\}_{i=1}^n) = \frac{\prod_{i\in\mathcal{C}}\xi_i}{\sum_{\mathcal{C}'\subset[n]:|\mathcal{C}'|=L_k}\prod_{i\in\mathcal{C}'}\xi_i}.\] This design makes within-group comparisons more frequent while still allowing cross-group comparisons.

Given \(C_k\), we generate latent scores according to one of the following scenarios and observe only the induced ranking \(\boldsymbol{\pi}_{\text{obs},k}=\text{rank}(\mathbf{y}_{\mathcal{C}_k,k})\). We compare MASTER with the maximum likelihood estimator under the Plackett–Luce model [11] and a Borda counting method [24]. Throughout the three scenarios, the true global ranking is \(\pi_i^*=i\) and \(\theta_i=3i/n\).

1. Group-structured heteroscedastic scores. Given \(\mathcal{C}_k\), generate \(y_{ik}=\theta_i+\sigma_{g(i)}\epsilon_{ik}\) for \(i\in\mathcal{C}_k\), where \(\epsilon_{ik}\sim N(0,1)\) independently, \(g(i)\in\{1,2\}\) is the group label of item \(i\), and \(\sigma_1=0.3,\sigma_2=3\).

2. Bottom-block heavy-tail score model. Let \(B=\{1,\ldots,\lfloor bn\rfloor\}\) with \(b=0.25\). Given \(\mathcal{C}_k\), generate \[y_{ik}=\theta_i+M1_{\{i\in B\}}\eta_{ik}+ J1_{\{i\notin B\}}+\sigma\epsilon_{ik},\] where \(\eta_{ik}\sim\text{Bernoulli}(q)\), \(\epsilon_{ik}\sim N(0,1)\), independently across \(i\in\mathcal{C}_k\) and \(k\in[p]\). We use \(q=0.3\), \(M=10\), \(J = 0.05\) and \(\sigma=0.03\).

3. Plackett–Luce model. Given \(\mathcal{C}_k\), generate \(y_{ik}=\theta_i+\epsilon_{ik}\) for \(i\in\mathcal{C}_k\), where \(\epsilon_{ik}\sim\text{Gumbel}(0,1)\) independently across \(i\) and \(k\). Then the induced ranking on \(\mathcal{C}_k\) follows a Plackett–Luce model with item parameters \(\{\exp(\theta_i):i\in\mathcal{C}_k\}\).

All three scenarios satisfy the EWST assumption. Scenario 1 introduces group-dependent heteroscedasticity and therefore violates the homogeneous PL structure. Scenario 2 further creates a top-choice bias: items in the bottom block occasionally receive large shocks and appear near the top of some observed rankings. The additional shift \(J\) for items outside the bottom block ensures that the pairwise majority direction remains consistent with the true ranking \(\boldsymbol{\pi}^*\), so the EWST condition is still satisfied. Scenario 3 is correctly specified by the Plackett–Luce model and serves as a favorable benchmark for the PL MLE.

Table 2 shows that the Borda method performs poorly in all three scenarios. This is mainly due to the non-uniform sampling design: stronger items are more likely to be compared with stronger items, and weaker items are more likely to be compared with weaker items. As a result, Borda-type scores mainly reflect local winning rates rather than global ranking positions. In the first two scenarios, MASTER substantially outperforms the PL-MLE, with error reduction factors ranging from about \(2.3\) to \(62.1\). This is because the PL structure is violated by the heterogeneous score distributions in Scenario 1 and by the bottom-block heavy-tail effect in Scenario 2, while the EWST assumption still holds. In contrast, Scenario 3 is correctly specified by the PL model, and hence the PL-MLE performs the best. Nevertheless, MASTER remains competitive and substantially improves over the Borda-type methods, reflecting the efficiency loss of using a model-free estimator compared with a correctly specified parametric model.

3.2 Real data analysis↩︎

Tennis data (pairwise comparison). We apply the same methods considered in the simulation in Section 3.1 to a professional tennis match dataset, which contains the results of all men’s matches organized by the Association of Tennis Professionals (ATP) from 2000 to 2018, collected from http://www.tennis-data.co.uk. The data covers Grand Slams, ATP Masters 1000, and other professional events. The same data has been considered in [6] and [25]. After removing players who lost all games, the set contains 875 players, with 7.13% of player pairs having competed at least once. Here, each match is treated as a pairwise comparison event with \(L_k=2\); thus, the proposed method is directly applicable even though the data do not involve human raters in the usual sense.

Table 3 gives the top 10 players ranked by each method. The top 10 lists provided by the counting and USVT methods do not align well with people’s impression about the top male tennis players. For example, the USVT method ranks Roddick A. third, ahead of, for example, Djokovic, N., and the counting algorithm ranks Ferrer D. fourth, ahead of, for example, Murray, A., which does not seem sensible. In fact, we do not expect the counting and USVT methods to work well for the current data, because these methods are based on the assumption that the missingness probabilities for the competing pairs are uniform. This assumption is unlikely to hold here as the likelihood of two professional tennis players meeting in a match depends on their ranking. The BT model and the MASTER do not need this assumption.

The rankings produced by the MASTER and the BT model are positively associated with substantial differences. Kendall’s tau and Spearman’s rho rank correlations between the two rankings are 0.47 and 0.65, respectively. We provide some discussions about the differences between the top 10 player lists from the two methods. First, the rivalry between Nadal, R. and Federer, R. is one of the most significant in tennis history, which also appears here. The BT model ranks Federer as the second and Nadal as the third, while the MASTER reverses this ordering. Considering that Nadal has won 24 out of 40 head-to-head competitions against Federer during their careers, and 10 out of 14 Grand Slam competitions, we believe that ranking Nadal ahead of Federer makes more sense when SST is not assumed. Second, it is interesting to see that the BT model ranks Roddick, A., as the fifth and Hewitt, L., as the eleventh, while the MASTER reverses this ordering. We believe that Roddick’s high ranking under the BT model is boosted by his strong head-to-head record against Djokovic, N., who is ranked first by both methods. Roddick and Djokovic competed nine times, out of which Roddick won five times. This strong record carries significant weight in the likelihood of the BT model, resulting in a high latent score estimate for Roddick. In contrast, Hewitt competed against Djokovic seven times, winning only once. To better understand the rankings, we check some head-to-head records. The record of Roddick against Hewitt is tied at 7:7. However, the records of Roddick and Hewitt against Agassi, A. and del Potro, J.M., who rank sixth and seventh by both methods, show that having Hewitt in the fifth place makes more sense. Specifically, the records are 1:5, 1:4, 4:4, and 3:2, for Roddick against Agassi, Roddick against del Potro, Hewitt against Agassi, and Hewitt against del Potro, respectively, which suggests that ranking Roddick ahead of Agassi and del Potro is not sensible.

Sushi data (multi-way comparison). We further apply MASTER to the Sushi Preference Dataset², which contains the preference records of 5,000 users regarding a set of 100 different types of sushi, ranging from common seafood to regional specialties. The data were collected through a large-scale questionnaire survey, in which each of the 5,000 respondents ranked 10 types of sushi from most to least preferred. We compare MASTER with the estimate from the PL model.

The Kendall’s tau correlation between the estimated rankings from the MASTER and PL models is 0.81, indicating that the two methods exhibit a strong overall consensus on the dataset’s broader preference structure. We further compare the top ten items given by the two methods, as shown in Table 4, where the sushi types are ordered from the first to the tenth for each method. As a reference, we also show the pairwise winning rates between the top ten sushi under MASTER in Table 5, where each entry shows the proportion that the sushi type listed in the corresponding row is preferred over that in the corresponding column, among raters who have rated both types of sushi. For example, it shows that “toro" is preferred over”chu_toro" 49% of the raters who have rated both. As shown in Table 4, the top two rankings for “chu_toro" and”toro" are consistent for both methods. In addition, the ranking of the sushi types in the third to the tenth are also largely consistent. However, there are also notable differences. For example, the sushi type “maguro" is ranked third by the PL model but ranked eighth by the MASTER. Its third rank by the PL model may be boosted by the PL model’s parametric assumption. As we can see from Table 5, the eighth ranking of”maguro" perfectly aligns with the WST assumption, in the sense that its winning rate is below 0.5 when compared with the first seven sushi types and above 0.5 when compared with the last two. Another notable difference is the ranking of “kurumaebi”, which is ranked fifth by the PL model but to the tenth by MASTER. The higher ranking under the PL model may be boosted by its high winning rate (0.69) against “negi_roro", which is ranked high by both methods, under the restricted parametric assumption of the PL model. On the other hand, it is consistently less preferred over the types”tarabagani", “maguro", and”toro_samon" in the seventh to the ninth ranks under the MASTER, which leads to its tenth aggregated rank under the MASTER.

4 Discussions↩︎

In this paper, we propose a novel framework for rank aggregation that operates under substantially weaker assumptions than existing methodologies. We establish the theoretical properties of the proposed estimator, including its statistical consistency and error bounds. To facilitate implementation, we develop an efficient algorithm for computing the estimator and demonstrate its advantages over state-of-the-art methods through extensive simulation studies. Furthermore, the practical utility of our approach is validated by applying it to rank professional tennis players using data from the Association of Tennis Professionals (ATP) and to analyze sushi preference data. In both cases, the proposed method yields meaningful results.

One limitation of the MASTER framework is the computational challenge posed by its discrete optimization objective. While the proposed search algorithm performs effectively within our simulation settings, its performance in higher-dimensional or more general configurations warrants further scrutiny. To address this, future work should explore more sophisticated stochastic search techniques, such as simulated annealing. By incorporating controlled randomness, such algorithms could better navigate the complex discrete landscape to escape local optima and more reliably identify the global maximizer.

A distinct advantage of our approach, which remains underexploited, is its natural compatibility with privacy-preserving randomization. Protecting rater privacy is paramount when handling preference data, as individual rankings often reveal sensitive personal information. As our assumptions are highly nonparametric, they are likely to hold even after data perturbation. Consequently, the MASTER approach may maintain its theoretical consistency and error bounds under privatization without the need for the complex debiasing procedures required by parametric models (see, e.g., [26]), whereas randomization often leads to violation of the model’s structural assumptions. We leave the formal development of these privacy-preserving extensions for future investigation.

5 Proof of Theoretical Results↩︎

References↩︎

---

1. Haoran Zhang is with Department of Statistics and Data Science, Southern University of Science and Technology, Shenzhen, China (e-mail: zhanghr@sustech.edu.cn). Yunxiao Chen is with Department of Statistics, London School of Economics and Political Science, London, U.K. (e-mail:y.chen186@lse.ac.uk).↩︎

2. Data are available through the link: https://www.kamishima.net/sushi/↩︎