Parameter Estimation for Thurstone Choice Models

TL;DR

This paper analyzes the accuracy of parameter estimation in Thurstone choice models, including the Luce and Bradley-Terry models, providing theoretical bounds for maximum likelihood and rank-breaking estimators.

Contribution

It offers a tight theoretical characterization of the mean squared error for maximum likelihood and rank-breaking estimators in Thurstone models, extending understanding of estimation accuracy.

Findings

Derived tight bounds for MLE estimation error

Analyzed rank-breaking method under independence assumptions

Connected estimation accuracy to classification error in binary parameters

Abstract

We consider the estimation accuracy of individual strength parameters of a Thurstone choice model when each input observation consists of a choice of one item from a set of two or more items (so called top-1 lists). This model accommodates the well-known choice models such as the Luce choice model for comparison sets of two or more items and the Bradley-Terry model for pair comparisons. We provide a tight characterization of the mean squared error of the maximum likelihood parameter estimator. We also provide similar characterizations for parameter estimators defined by a rank-breaking method, which amounts to deducing one or more pair comparisons from a comparison of two or more items, assuming independence of these pair comparisons, and maximizing a likelihood function derived under these assumptions. We also consider a related binary classification problem where each individual parameter takes value from a set of two possible values and the goal is to correctly classify all items within a prescribed classification error.

Tables

Table 1. Table 1: The values of parameters for our examples of 𝒯 F subscript 𝒯 𝐹 {\mathcal{T}}_{F} .

$F$	$\frac{\partial p_k\left(\mathrm{𝟎}\right)}{\partial x_1}$	${\gamma }_{F,k}$
Gaussian	$O\left(\frac{1}{k^{2-ϵ}}\right)$	$\Omega \left(\frac{1}{k^{2ϵ}}\right)$
Double-exponential	$\frac{1}{\beta k^2}$	${\beta }^2\frac{k}{k-1}$
Laplace	$\frac{1-1/2^{k-1}}{\beta k\left(k-1\right)}$	${\beta }^2\frac{k-1}{k{\left(1-1/2^{k-1}\right)}^2}$
Uniform	$\frac{1}{2a\left(k-1\right)}$	$4a^2\frac{k-1}{k^3}$

Equations

$$p_{i,S}(\theta )=p_k({\mathbf{x}}_{i,S}(\theta )),fori\in S$$

$$p_k(\mathbf{x})={\int }_{\mathbf{R}}\left(\prod _{v=1}^{k-1}F(x_v+z)\right)f(z)dz,for\mathbf{x}\in {\mathbf{R}}^{k-1}.$$

$$p_{i,j}(\theta )=p_2({\theta }_i-{\theta }_j)$$

$$p_2(x)={\int }_{\mathbf{R}}F(x+z)f(z)dz,forx\in \mathbf{R}.$$

$$p_k(\mathbf{x})=\frac{1}{1+{\sum }_{i=1}^{k-1}{e}^{-x_i/\beta }},for\mathbf{x}\in {\mathbf{R}}^{k-1}.$$

$$\ell (\theta )=\sum _{t=1}^m\log (p_{y_t,S_t}(\theta )).$$

$$\ell (\theta )=\sum _{S\subseteq N}\sum _{i\in S}{w}_{i,S}\log (p_{|S|}({\mathbf{x}}_{i,S}(\theta )))+const$$

$$\ell (\theta )=\sum _{i=1}^n\sum _{j=1}^n{w}_{i,j}\log \left(p_2({\theta }_i-{\theta }_j)\right)+const.$$

$${\ell }_{k-1}(\theta )=\sum _{t=1}^m\sum _{v\in S_t\setminus \{y_t\}}\log (p_{y_t,v}(\theta )).$$

$${\ell }_1(\theta )=\sum _{t=1}^m\log (p_{y_t,z_t}(\theta )).$$

$$MSE(\hat{\theta },{\theta }^{\star })=\frac{1}{n}\Vert \hat{\theta }-{\theta }^{\star }{\Vert }_2^2.$$

$$L_{\mathbf{A}}=\mathrm{diag}(A1)-\mathbf{A}$$

$$m_{i,j}=\frac{n}{m}\sum _{k\ge 2}w(k)m_{i,j}(k).$$

$$m_{i,j}=\frac{1}{n-1}\sum _{k\ge 2}w(k)k(k-1)\mu (k).$$

$${\lambda }_2(L_{{\mathbf{M}}_w})=\left(1-\frac{1}{n}\right)\sum _{k\ge 2}w(k)k(k-1)\mu (k)$$

$$\Vert \hat{\theta }-{\theta }^{\star }{\Vert }_2\le \frac{2{\Vert \nabla g({\theta }^{\star })\Vert }_2}{{\min }_{\theta \in \Theta }{\lambda }_2({\nabla }^{2}g(\theta ))}.$$

$$\mathbf{P}\left[{\Vert S_m\Vert }_2\ge x\right]\le 2e^2{e}^{-\frac{x^2}{2m{\sigma }^2}}.$$

$$\mathbf{P}\left[{\lambda }_1(S_m)\le (1-ϵ){\lambda }_1(\mathbf{E}[S_m])\right]\le n{e}^{-\frac{{ϵ}^2{\lambda }_1(\mathbf{E}[S_m])}{2\sigma }}.$$

$$\frac{d^2\log (p_2(x))}{d{x}^2}\le -Afor\; allx\in [-2b,2b].$$

$$\frac{d\log (p_2(x))}{dx}\le Bfor\; allx\in [-2b,2b].$$

$$MSE(\hat{\theta },{\theta }^{\star })\le D^2\frac{n(\log (n)+2)}{{\lambda }_2(L_{{\mathbf{M}}_{1/4}}{)}^2}\frac{1}{m}$$

$$m\ge \frac{1}{{ϵ}^2}{D}^2\frac{1}{{\lambda }_2(L_{{\mathbf{M}}_{1/4}}{)}^2}n(\log (n)+2).$$

$$m\ge \frac{4}{{ϵ}^2}{D}^{2}n(\log (n)+2).$$

$$\ell (\theta )=\sum _{t=1}^m\log (p_2({\theta }_{y_t}-{\theta }_{z_t}))$$

$$\underset{\theta \in \Theta }{\min }{\lambda }_2({\nabla }^2(-\ell (\theta )))\ge 4A\frac{m}{n}{\lambda }_2(L_{{\mathbf{M}}_{1/4}}).$$

$$\frac{d^2\log (p_2({\theta }_i-{\theta }_j))}{d{x}^2}$$

$$\frac{\partial^{2}(-\log(p_{2}(\theta_{y_{t}}-\theta_{z_{t}})))}{\partial\theta_{i}\partial\theta_{j}}=\left\{\begin{array}[]{ll}\frac{d^{2}}{dx^{2}}\log(p_{2}(\theta_{y_{t}}-\theta_{z_{t}})),&\hbox{ if }\{i,j\}=\{y_{t},z_{t}\}\\ 0,&\hbox{ otherwise}\end{array}\right.$$

$$\frac{{\partial }^2(-\log (p_2({\theta }_{y_t}-{\theta }_{z_t})))}{\partial {\theta }_i^2}=-\sum _{j\ne i}\frac{{\partial }^2(-\log (p_2({\theta }_{y_t}-{\theta }_{z_t})))}{\partial {\theta }_i\partial {\theta }_j}.$$

$$\frac{{\partial }^2(-\log (p_2({\theta }_{y_t}-{\theta }_{z_t})))}{\partial {\theta }_i\partial {\theta }_j}\le -A{1}_{\{y_t,z_t\}=\{i,j\}}.$$

$${\nabla }^2(-\ell (\theta ))⪰4A\frac{m}{n}{L}_{{\mathbf{M}}_{1/4}},for\; all\theta \in [-b,b{]}^n$$

Full text

We provide a tight characterization of the mean squared error of the maximum likelihood parameter estimator. We also provide similar characterizations for parameter estimators defined by a rank-breaking method, which amounts to deducing one or more pair comparisons from a comparison of two or more items, assuming independence of these pair comparisons, and maximizing a likelihood function derived under these assumptions. We also consider a related binary classification problem where each individual parameter takes value from a set of two possible values and the goal is to correctly classify all items within a prescribed classification error.

The results of this paper shed light on how the parameter estimation accuracy depends on given Thurstone choice model and the structure of comparison sets. In particular, we found that for unbiased input comparison sets of a given cardinality, when in expectation each comparison set of given cardinality occurs the same number of times, for a broad class of Thurstone choice models, the mean squared error decreases with the cardinality of comparison sets, but only marginally according to a diminishing returns relation. On the other hand, we found that there exist Thurstone choice models for which the mean squared error of the maximum likelihood parameter estimator can decrease much faster with the cardinality of comparison sets.

We report empirical evaluation of some claims and key parameters revealed by theory using both synthetic and real-world input data from some popular sport competitions and online labor platforms.

1 Introduction

We consider the statistical inference problem of estimating individual strength or skill parameters of items based on observed choices of items from sets of two or more items. This accommodates the case of pair comparisons as a special case, where each comparison set consists of two items. In our more general case, each observation consists of a comparison set of two or more items and the identity of the chosen item from this set. In other words, each observation is a partial ranking, which is often referred to as a top-1 list. Many applications are accommodated by this framework, e.g., choices indicated by user clicks in various information retrieval systems, outcomes of single-winner contests in crowdsourcing services such as TopCoder or Taskcn, outcomes of hiring decisions where one applicant is hired among those who applied for a job, e.g., in online labour markets such as Fiverr and Upwork, as well as numerous sport competitions and online gaming platforms.

We consider the parameter estimation for the statistical choice model known as the Thurstone choice model; also referred to as the random utility model. According to the Thurstone choice model, items are associated with latent performance random variables that are independent across different items and different comparisons. For any given comparison set, the choice is the item from this set with the largest performance random variable. For any given item, the performance random variable is equal to the sum of a strength parameter, whose value can be specific to this item, and a noise random variable according to a given cumulative distribution function. The values of the strength parameters are unknown and have to be estimated from the observed comparisons, and the distribution of noise random variables is assumed to be known. The Thurstone choice model accommodates many known choice models by admitting different distribution for noise random variables, e.g., Luce choice model (Luce (1959)) for comparison sets of two or more items, and its special case for pair comparisons, often referred to as the Bradley-Terry model (Bradley and Terry (1954)).

In this paper, we study the accuracy of the maximum likelihood estimator of the parameter of the Thurstone choice model. Our goal is to characterize the accuracy of the maximum likelihood estimator and shed light on how it depends on the given Thurstone choice model and properties of the observed input data such as the number of observations and the structure of comparison sets. In particular, we consider the following statistical inference question. Suppose that the input observations are such that all comparison sets are of the same cardinality $k\geq 2$ and are unbiased, meaning that in expectation every comparison set of given cardinality occurs the same number of times in the input data. Then, we would like to understand how does the accuracy of the maximum likelihood estimator of the strength parameters depend on the cardinality of comparison sets. Notice that from any comparison set of cardinality $k$, we can deduce at most $k(k-1)/2$ pair comparisons. Intuitively, we would expect that the parameter estimation accuracy would increase with the cardinality of comparison sets. However, it is not a priori clear how fast the accuracy would improve and whether any significant gains can be achieved by increasing the sizes of comparison sets. Moreover, it is also not a priori clear whether or not there can be any significant difference between different Thurstone choice models, with respect to how the parameter estimation accuracy is related to the cardinality of comparison sets. We also consider these questions for parameter estimators that are derived by rank breaking methods, which amount to deducing one or more pair comparisons from each comparison of two or more items, assuming independence of these pair comparisons, and defining the estimator as the maximum likelihood estimator under these assumptions.

The main contributions of this paper can be summarized as follows.

We provide upper bounds on the mean squared error of the maximum likelihood estimator, and a lower bound that establishes their minimax optimality. We show that the effect of the structure of comparison sets on the mean squared error is captured by one key parameter: algebraic connectivity of a suitably defined weighted-adjacency matrix. The elements of this matrix correspond to distinct pairs of items and are equal to a weighted sum of the number of input comparisons of different cardinalities that contain the corresponding pair of items, where the weights are specific to given Thurstone choice model.

For the statistical inference question of how the estimation accuracy improves with the cardinality of comparison sets, we derive a tight characterization of the mean squared error in terms of the cardinality of comparison sets (Corollary 12). This characterization reveals that for a broad class of Thurstone choice models, which includes the well-known cases such as the Luce choice model, there is a diminishing returns decrease of the mean squared error with the cardinality of comparison sets. For this class of Thurstone choice models, the mean squared error tends to be largely insensitive to the cardinality of comparison sets. On the other hand, we show that there exist Thurstone choice models for which the mean squared error decreases much faster with the cardinality of comparison sets. Perhaps suprisingly, in these cases, the amount of information extracted from a comparison set of cardinality $k$ is in the order of $k^{2}$ independent pair comparisons, which yields a $1/k^{2}$ reduction of the mean squared error of the maximum likelihood estimator. Section 6 provides more discussion.

We consider two natural rank-breaking methods, one that deduces $k-1$ pair comparisons and one that deduces $1$ pair comparison from a comparison set of cardinality $k$ (Section 4). We derive mean squared error upper bounds when choices are according to the Luce choice model for these two rank-breaking methods in Theorem 17 and Theorem 18, respectively. These results show that both estimators are consistent. Interestingly, both mean squared error upper bounds are equal to that of the maximum likelihood estimator up to a constant factor. Hence, they both inherit all the properties that we established to hold for the mean squared error upper bound for the maximum likelihood estimator.

We also consider a binary classification problem where all strength parameters associated with items take one of two possible values (separating items into two classes), and the goal is to correctly classify each item within a prescribed probability of classification error (Section 5). We identify sufficient conditions for correctness of a simple point score classification algorithm (Theorem 19) and establish their tightness (Theorem 20). These conditions are of the same form as those that we imposed for deriving upper bounds on the mean squared error of the maximum likelihood parameter estimator.

We present experimental results using both simulation and real-world data (Section 7). In particular, we validate the claim that the mean squared error can decrease with the cardinality of comparison sets in a qualitatively different way depending on the given Thurstone choice model. We also evaluate algebraic connectivity of weighted-adjacency matrices for several real-world input data, demonstrating that it can cover a wide range of values depending on specific application scenario.

1.1 Related Work

A model of comparative judgement for pair comparisons was introduced by Thurstone (1927), which is a special case of a model that we refer to as a Thurstone choice model, for the case of pair comparisons and Gaussian random noise variables. A statistical model of pair comparisons that postulates that an item is chosen from a set of two items with probability proportional to the strength parameter of this item was introduced by Zermelo (1929), and was then popularized by the work of Bradley and Terry (1952, 1954) and others, and is often referred to as the Bradley-Terry model. The statistical model of choice where for any set of two or more items an item is chosen with probability proportional to its strength parameter was shown to be a unique model satisfying a set of axioms introduced by Luce (1959), and is referred as the Luce choice model. The Bradley-Terry model is the special case of the Luce choice model for pair comparisons. The choice probabilities of the Luce choice model correspond to those of a Thurstone choice model with noise random variables according to a double-exponential distribution. Relationships between the Luce choice model and Thurstone choice model were studied by Yellott (1977). A statistical model for full ranking outcomes (the outcome of a comparison is an ordered list of the compared items) where the ranking is in the order of sampling of items without replacement from the set of compared items and the sampling probabilities are proportional to the strengths of items is referred as the Plackett-Luce model (Luce (1959) and Plackett (1975)).

The Thurstone choice models have been used in the design of several popular skill rating systems, e.g., Elo rating system by Elo (1978) used for skill rating of chess players and in other sports, and TrueSkill by Graepel et al. (2006) used for skill rating of gamers of a popular online gaming platform. All these models are instances of Thurstone models, and are special instances of generalized linear models, see, e.g., Nelder and Wedderburn (1972), McCullagh and Nelder (1989), and Chapter 9 in  Murphy (2012). An exposition to the principles of skill rating systems is available in Chapter 9 in Vojnović (2016).

The parameter estimation problem for the Bradley-Terry model of pair comparisons has been studied by many. The iterative methods for computing a maximum likelihood parameter estimate have been studied in the early work by Hunter (2004) and the recent work by Maystre and Grossglauser (2015). Simons and Yao (1999) shown that the maximum likelihood parameter is consistent and asymptotically normal as the number of items $n$ grows large, under assumption that each pair is compared the same number of times and that the true parameter vector is such that the maximum distance between any of its coordinates is $o(\log(n))$. Maydeau-Olivares (1999) studied Thurstonian model parameter estimation with noise random variables according to a Gaussian distribution.

The accuracy of the parameter estimation for various instances of Thurstone models has been studied in recent work. Negahban et al. (2012) found a sufficient number of input pair comparisons to achieve a given mean squared error of a parameter estimator for the Bradley-Terry model and the input comparisons such that in expectation each distinct pair is compared the same number of times. In particular, they shown that under given assumptions, it suffices to observe $O(n\log(n))$ input pair comparisons. Rajkumar and Agarwal (2014) studied a statistical convergence property of ranking aggregation algorithms for pair comparisons not only for the Bradley-Terry model but also under some more general conditions referred to as low-noise and generalized low-noise. Hajek et al. (2014) provided a characterization of the mean squared error of the maximum likelihood parameter estimate for the Plackett-Luce model of full ranking outcomes. This work found that the algebraic connectivity of a weighted-adjacency matrix captures the effect of the structure of input comparison sets on the mean squared error of the maximum likelihood parameter estimator. Shah et al. (2016) established similar characterization results for the case of pair comparisons according to Thurstone choice models.

Our work differs from previous work in that we consider characterization of the estimation accuracy for a general class of Thurstone choice models for arbitrary sizes of comparisons. Specifically, our work provides a first characterization of the estimation accuracy with respect to the cardinality of comparison sets for unbiased input comparisons, which reveals an insight into the fundamental limits of statistical inference for given cardinality of comparison sets. It is a folklore that different models of pair comparisons yield similar performance with respect to the prediction error, e.g. Stern (1992), which suggests that the precise choice of a Thurstone choice model does not matter much in applications. Our results show that there can be significant difference between two Thurstone choice models with respect to the statistical inference of model parameters.

Parameter estimators derived by various rank-breaking methods have been studied by various authors. For instance, Soufiani et al. (2013) and Soufiani et al. (2014) studied rank-breaking methods for full ranking data and Khetan and Oh (2016b) studied rank-breaking methods for partial rankings. Recently, Khetan and Oh (2016a) characterized a trade-off between the amount of information used per comparison and the mean squared error of a parameter estimate based on rank breaking. Our work is different in that we are interested in top-1 list observations and the effect of the structure of comparison sets for rank-breaking methods. For the top-1 list observations and any given comparison structure, our work provides upper bounds for the mean squared error of two natural ranking breaking methods and shows the optimality of them.

1.2 Organization of the Paper

Section 2 introduces problem formulation, some basic concepts and key technical results used to establish our main results. Section 3 provides a characterization of the mean squared error for the maximum likelihood parameter estimation, including both upper and lower bounds. Section 4 shows the same type of characterizations for two rank-breaking based parameter estimators. Section 5 establishes tight conditions for correct classification of items, when the strength parameters of items are of two possible types. Section 6 discusses how the estimation accuracy depends on the cardinality of comparison sets. Section 7 contains results of our experiments. Finally, Section 8 concludes the paper. Appendix contains some background facts and proofs of our theorems.

2 Problem Formulation and Notation

Let $N=\{1,2,\dots,n\}$ be a set of two or more items. The input data consists of a sequence of one or more observations $(S_{1},y_{1})$, $(S_{2},y_{2})$, $\ldots$, $(S_{m},y_{m})$, where each observation $t$ consists of (a) a comparison set $S_{t}\subseteq N$ and (b) a choice of an item $y_{t}$ from $S_{t}$. The case of pair comparisons is accommodated as a special case when each comparison set consists of two items. Let $w_{i,S}$ denote the number of observations for which the comparison set is $S$ and the chosen item is $i\in S$. With a slight abuse of notation, for pair comparisons, let $w_{i,j}$ be the number of observations for which the comparison set is $\{i,j\}$ and the chosen item is $i$.

A Thurstone choice model, denoted as ${\mathcal{T}}_{F}$, is defined by a cumulative distribution function $F$ of a zero mean random variable and parameter vector $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{n})$ that takes value in $\Theta\subseteq\mathbf{R}^{n}$. We can interpret $\theta_{i}$ as the strength of item $i\in N$. The cumulative distribution function $F$ is assumed to have a density function, denoted by $f$.

According to Thurstone choice model ${\mathcal{T}}_{F}$, for any given sequence of comparison sets, choices are independent random variables according to the following distribution: conditional on that the comparison set is $S$, with $k$ denoting the cardinality of $S$, the distribution of choice is

[TABLE]

where $\mathbf{x}_{i,S}(\theta)$ is a vector in $\mathbf{R}^{k-1}$ with elements $\theta_{i}-\theta_{j}$, for $j\in S\setminus\{i\}$, and

[TABLE]

With a slight abuse of notation, for the case of pair comparisons, we let $p_{i,j}(\theta)$ denote the probability that item $i$ is chosen from the comparison set $\{i,j\}$. In this case, we have

[TABLE]

where

[TABLE]

A Thurstone choice model ${\mathcal{T}}_{F}$ corresponds to the following probabilistic generative model of choice, also referred to as a random utility model. The items of each comparison set are associated with latent performance random variables, which are independent across different items and different comparison sets. For any comparison set $S$ and item $i\in S$, the performance random variable $X_{i}$ is equal to the sum of a strength parameter $\theta_{i}$ and a noise random variable with distribution $F$. The given probabilistic generative model assumes that for any given comparison set $S$, the chosen item is the one with the largest performance. Hence, the distribution of choice is $p_{i,S}(\theta)=\mathbf{P}[X_{i}\geq\max_{j\in S}X_{j}]$ for $i\in S$, which can be expressed as asserted in (1).

It is noteworthy that under a Thurstone choice model, the probability distribution of choice depends only on pairwise differences of the strength parameters. This implies that the probability distribution of choice is shift-invariant with respect to the parameter vector. In order to allow for identifiability of the parameter vector, we assume that $\theta$ is such that $\sum_{i=1}^{n}\theta_{i}=0$.

We refer to several examples of Thurstone choice models as follows: (i) Gaussian distribution : $f(x)=1/(\sqrt{2\pi}\sigma)\,\exp(-x^{2/(2\sigma^{2}}))$ with variance $\sigma^{2}$; (ii) Double-exponential distribution: $F(x)=\exp(-\exp(-(x+\beta\gamma)/\beta))$ with parameter $\beta>0$ and $\gamma$ denoting the Euler-Mascheroni constant, which has variance $\sigma^{2}=\pi^{2}\beta^{2}/6$; (iii) Laplace distribution: $F(x)=1/2\,\exp(x/\beta)$, for $x<0$, and $F(x)=1-1/2\,\exp(-x/\beta)$, for $x\geq 0$, with parameter $\beta$, which has variance $\sigma^{2}=2\beta^{2}$; and (iv) Uniform distribution: $f(x)=1/(2a)$ for $x\in[-a,a]$, which has variance $\sigma^{2}=a^{2}/3$.

For the general case of comparison sets of two or more items, the distribution of choice admits an explicit form only for some special cases. For example, when noise random variables are according to a double-exponential distribution, we have

[TABLE]

This amounts to the choice probabilities $p_{i,S}(\theta)=\exp(\theta_{i}/\beta)/\sum_{v\in S}\exp(\theta_{v}/\beta)$, for $i\in S$ and $S\subseteq N$, which under suitable re-parametrisation corresponds to the well-known Luce choice model. In particular, for pair comparisons, we have $p_{2}(x)=1/(1+\exp(-x/\beta))$, which under suitable re-parametrisation corresponds to the well-known Bradley-Terry model. For pair comparisons, the choice probabilities admit an explicit form also for some other Thurstone choice models; for example, when noise has a Gaussian distribution, we have $p_{2}(x)=\Phi(x/(2\sigma))$ where $\Phi$ is the cumulative distribution function of a standard normal random variable.

Maximum Pseudo Likelihood Estimation

We consider parameter estimators that are defined as maximizers of a pseudo log-likelihood function $\widetilde{\ell}:\Theta\rightarrow\mathbf{R}$. We refer to the parameter estimator $\widehat{\theta}\in\arg\max_{\theta\in\Theta}\widetilde{\ell}(\theta)$ as a maximum pseudo likelihood estimator.

We devote a special attention to maximum likelihood estimator, defined as a maximizer of the log-likelihood function, which for a Thurstone choice model is given by

[TABLE]

The log-likelihood function can be written as follows

[TABLE]

where recall $w_{i,S}$ is the number of observations for which the comparison set is $S$ and $i$ is the choice from $S$. In particular, for pair comparisons, we have

[TABLE]

Evaluating the value of the log-likelihood function in (4) for given parameter vector requires evaluating a sum that in the worst-case consists of exponentially many elements in $n$ (all possible combinations of two or more elements from the ground set of $n$ elements). On the other hand, for pair comparisons, the log-likelihood function in (5) is a sum of at most $n^{2}$ elements; thus, polynomially many elements in $n$. A common approach to reduce computational complexity is to use the so-called rank breaking, which amounts to deducing pair comparisons from any given comparison set of two or more items, and assuming that these pair comparisons are independent (if this is not the case). Using these pair comparisons, one then defines a pseudo log-likelihood function as the log-likelihood function under the assumption that the deduced pair comparisons are independent.

We shall consider two natural rank-breaking methods. The first rank-breaking method deduces $k-1$ pair comparisons from each comparison set of $k$ items, by taking all pairs that consist of the chosen item and each other item in the given comparison set. The pseudo log-likelihood function in this case is given by

[TABLE]

The second rank-breaking method that we consider deduces $1$ pair comparison from each comparison set of $k$ items, by taking a pair that consists of the chosen item and a randomly picked item from the remaining set of items in the given comparison set. The pseudo log-likelihood function in this case is given by

[TABLE]

The first rank-breaking method uses maximum amount of information that is contained in a comparison; by observing choice of one item from a comparison set of $k$ items, we can indeed deduce at most $k-1$ pair comparisons (between the chosen item and each other item in the given comparison set). In general, these pair comparisons are not mutually independent. The second rank-breaking method is conservative in deducing only one pair comparison from each comparison set of two or more items.

Parameter Estimation Accuracy

We study the accuracy of a maximum pseudo log-likelihood estimator $\widehat{\theta}$ of the true parameter vector $\theta^{\star}$ by using the mean squared error defined as follows:

[TABLE]

We also consider the probability of classification error for the case when the strength parameters belong to one of two classes and the goal is to correctly classify each item.

2.1 Eigenvalues, Adjacency, and Laplacian Matrices

Here we review some basic definitions that are used throughout the paper. We denote eigenvalues of a matrix $\mathbf{A}\in\mathbf{R}^{n\times n}$ as $\lambda_{1}(\mathbf{A}),\lambda_{2}(\mathbf{A}),\ldots,\lambda_{n}(\mathbf{A})$. By convention, we assume that $\lambda_{1}(\mathbf{A})\leq\lambda_{2}(\mathbf{A})\leq\cdots\leq\lambda_{n}(\mathbf{A})$. The spectral norm $||\mathbf{A}||_{2}$ of matrix $\mathbf{A}\in\mathbf{R}^{n\times n}$ is defined by $||\mathbf{A}||_{2}=\sqrt{\lambda_{n}(\mathbf{A}^{\top}\mathbf{A})}$. The spectral norm of $\mathbf{A}$ is induced by the Euclidean vector norm as follows $||\mathbf{A}||_{2}=\max\{||\mathbf{A}\mathbf{x}||_{2}:\mathbf{x}\in\mathbf{R}^{n},||\mathbf{x}||=1\}$. If $\mathbf{A}\in\mathbf{R}^{n\times n}$ is a real symmetric matrix, then $||\mathbf{A}||_{2}=\lambda_{n}(\mathbf{A})$.

For any weighted-adjacency matrix $\mathbf{A}\in\mathbf{R}^{n\times n}_{+}$, we consider a Laplacian matrix $L_{\mathbf{A}}$ defined by

[TABLE]

where for any given vector $\mathbf{a}$, $\mathrm{diag}(\mathbf{a})$ denotes the diagonal matrix with diagonal $\mathbf{a}$.

For any weighted-adjacency matrix $\mathbf{A}\in\mathbf{R}^{n\times n}_{+}$, we refer to $\lambda_{2}(L_{\mathbf{A}})$ as the Fiedler value of $\mathbf{A}$ (Fiedler (1973, 1989)). The Fielder eigenvalue of a weighted-adjacency matrix quantifies its algebraic connectivity. A weighted-adjacency matrix $\mathbf{A}$ corresponds to a connected graph if and only if it has strictly positive Fielder value, i.e., $\lambda_{2}(L_{\mathbf{A}})>0$.

For any given observations and given weight function $w:\{1,2,\ldots,n\}\rightarrow\mathbf{R}_{+}$, we define the weighted-adjacency matrix $\mathbf{M}_{w}$ as follows. Let $m_{i,j}(k)$ be the number of comparisons of cardinality $k$ that contain the pair of items $\{i,j\}$. Then, we define $\mathbf{M}_{w}$ to be the matrix in $\mathbf{R}_{+}^{n\times n}$ with zero diagonal elements and other elements given by

[TABLE]

With a slight abuse of notation, let $\mathbf{M}_{a}$ be the weighted-adjacency matrix defined for the weight function that takes constant value $a>0$, and let $\mathbf{M}$ be written in lieu of $\mathbf{M}_{1}$.

If all comparison sets have identical cardinalities, then each element of the weighted-adjacency matrix is equal to the number of comparisons that contain the corresponding pair of items up to a multiplicative factor. The factor $n/m$ can be interpreted as a normalization with the mean number of comparison sets per item. This normalization is admitted so that for the canonical case of pair comparisons when each pair is compared the same number of times, $\lambda_{2}(L_{\mathbf{M}_{a}})$ is a constant independent of the number of observations $m$ and the number of items asymptotically for large $n$. Indeed, in this case, $\lambda_{2}(L_{\mathbf{M}_{a}})=\cdots=\lambda_{n}(L_{\mathbf{M}_{a}})=(1-1/n)a$, which is equal to constant $a$, asymptotically for large $n$.

We say that comparison sets are unbiased if for any given cardinality, each set of the given cardinality occurs the same number of times. In particular, for pair comparisons, this means that each distinct pair is compared the same number of times. For any unbiased comparison sets, the weighted-adjacency matrix can be expressed as follows. Let $\mu(k)$ be the fraction of comparison sets of cardinality $k$. Then, for every integer $k\geq 2$ and pair of items $\{i,j\}$, $m_{i,j}(k)=\left(\binom{n-2}{k-2}/\binom{n}{k}\right)\mu(k)m=\left(k(k-1)/[n(n-1)]\right)\mu(k)m$. Hence, for every pair of items, we have

[TABLE]

It follows that for any unbiased comparison sets, we have

[TABLE]

which is a constant independent of $n$, asymptotically for large $n$.

We shall also consider comparison sets that are assumed to be an independent random sequence according to a given distribution. Specifically, we shall consider the case where all comparison sets are of the same cardinality, and are independent samples according to uniform random sampling without replacement from the set of all items. We denote with $\overline{\mathbf{M}}_{w}$ the expected weighted-adjacency matrix, where the expectation is with respect to the distribution of the sequence of comparison sets. We say that comparison sets are a priori unbiased if all non-diagonal elements of $\overline{\mathbf{M}}_{w}$ are equal.

2.2 A Key Lemma and Probability Tail Bounds

All upper bounds for the mean squared error of a maximum pseudo log-likelihood estimator in this paper are established by using the following key lemma.

Lemma 1

Suppose that $g:\mathbf{R}^{n}\rightarrow\mathbf{R}$ satisfies (i) $\nabla^{2}g(\mathbf{\theta})\mathbf{1}=\mathbf{0}$ and (ii) $\lambda_{2}(\nabla^{2}g(\mathbf{\theta}))>0$ for all $\mathbf{\theta}\in\Theta$, where $\Theta=\{\mathbf{\theta}\in[-b,b]^{n}:\mathbf{\theta}^{\top}\mathbf{1}=0\}$ for $b>0$. Let $\mathbf{\theta}^{\star}$ be an arbitrary vector in $\Theta$ and $\widehat{\mathbf{\theta}}\in\arg\min_{\mathbf{\theta}\in\Theta}g(\mathbf{\theta})$. Then, we have

[TABLE]

We shall apply this lemma to the case where $g$ is a negative pseudo log-likelihood function, $\widehat{\theta}$ is a maximizer of the pseudo log-likelihood function, and $\theta^{\star}$ is the true parameter vector. We upper bound the mean squared error by the following two steps:

• S1

find $\alpha>0$ such that $\left\|\nabla g(\mathbf{\theta}^{\star})\right\|_{2}\leq\alpha$, and

• S2

find $\beta>0$ such that $\min_{\mathbf{\theta}\in\Theta}\lambda_{2}(\nabla^{2}g(\mathbf{\theta}))\geq\beta$

which imply the upper bound $\|\widehat{\mathbf{\theta}}-\mathbf{\theta}^{\star}\|_{2}\leq 2\alpha/\beta$.

All our proofs of the mean squared estimation error upper bounds follow the above two-step procedure, including the proof of Theorem 4 in Section 3.1 and other proofs provided in Appendix.

In step S1, $\nabla g(\mathbf{\theta}^{\star})$ is a sum of random vectors. We will make use of the following vector version of Azuma-Hoeffding bound (Theorem 1.8 in Hayes (2003)) for a sum of random vectors.

Lemma 2** (vector Azuma-Hoeffding bound)**

Suppose that $S_{m}=\sum_{t=1}^{m}X_{t}$ is a martingale where $X_{1},X_{2},\ldots,X_{m}$ are random variables that take values in $\mathbf{R}^{n}$ and are such that $\mathbf{E}[X_{t}]=\mathbf{0}$ and $\|X_{t}\|_{2}\leq\sigma$ for all $t\in\{1,2,\ldots,m\}$, for $\sigma>0$. Then, for every $x>0$,

[TABLE]

In step S2, we need to find a lower bound for the second-smallest eigenvalue of the Hessian matrix $\nabla^{2}g(\mathbf{\theta})$ for all $\theta\in\Theta$. For pair comparisons according to a Thurstone choice model or comparisons of two or more items according to the Luce choice model, $\nabla^{2}g(\mathbf{\theta})$ is determined by comparison sets and does not depend on the choices. We can find $\beta$ from a Laplacian matrix when the comparison sets are given. In other cases, $\nabla^{2}g(\mathbf{\theta})$ is a sum of random matrices. We will make use the following matrix version of Chernoff’s bound along with properties of eigenvalues of a Laplacian matrix (which are given in Appendix A).

Lemma 3** (matrix Chernoff bound)**

Let $S_{m}=\sum_{t=1}^{m}X_{t}$ where $X_{1},X_{2},\ldots,X_{m}$ are random independent real symmetric matrices in $\mathbf{R}^{n\times n}$ such that $\lambda_{1}(X_{t})\geq 0$ and, $\|X_{t}\|_{2}\leq\sigma$ for $t\in\{1,2,\ldots,m\}$, for $\sigma>0$. Then, for $\epsilon\in[0,1)$,

[TABLE]

3 Maximum Likelihood Estimation

In this section, we present upper and lower bounds for the mean squared error of the maximum likelihood parameter estimator for the Thurstone choice model. We will first consider the case of pair comparisons. We then consider the more general case when each comparison set consists of two or more items. For this more general case, we first give an upper bound for the Luce choice model, and then present similar characterization for a class of Thurstone choice models. We end this section with a lower bound on the mean squared error of the maximum likelihood parameter estimator, which establishes minimax optimality of our upper bounds.

3.1 Pair Comparisons

We consider pair comparisons according to a Thurstone choice model ${\mathcal{T}}_{F}$ with parameter vector $\theta^{\star}$ that takes value in $\Theta=[-b,b]^{n}$ and that satisfies the following conditions:

• P1

There exists $A>0$ such that

[TABLE]

• P2

There exists $B>0$ such that

[TABLE]

• P3

The weighted-adjacency matrix $\mathbf{M}$ is irreducible, i.e. $\lambda_{2}(L_{\mathbf{M}})>0$.

Condition P1 means that $p_{2}$ is a strictly log-concave function on $[-2b,2b]$. Condition P2 means that $\log(p_{2}(x))$ has a bounded derivative on $[-2b,2b]$. Notice that this condition is equivalent to $dp_{2}(x)/dx\leq Bp_{2}(x)$ for all $x\in[-2b,2b]$. Constants $A$ and $B$ are specific for given Thurstone choice model and the value of the parameter $b$. In particular, for the Bradley-Terry model, it can be easily checked that P1 and P2 hold with $A=e^{-2b/\beta}/[\beta^{2}(1+e^{-2b/\beta})^{2}]$ and $B=1/[\beta(1+e^{-2b/\beta})]$. Condition P3 means that the observations are such that the graph defined by the weighted-adjacency matrix is connected. Equivalently, there exists no partition of the set of items $N$ into two non-empty sets $S$ and $N\setminus S$ such that some pair of items $i\in S$ and $j\in N\setminus S$ is not compared in the input observations.

Theorem 4

Under conditions P1, P2 and P3, with probability at least $1-2/n$,

[TABLE]

where $D=B/A$.

Before we show a proof of the theorem, we note the following remarks.

First, notice that $D$ is a constant whose value is specific to given Thurstone choice model and the value of parameter $b$. In particular, for the Bradley-Terry model, we have $D=\beta(e^{2b/\beta}+1)$.

Second, from (14), for the mean squared error to be smaller than or equal to $\epsilon^{2}$, for given $\epsilon>0$, it suffices that the number of observations is such that

[TABLE]

Third, and last, if each pair of items is compared the same number of times, then, from (10), we have $m_{i,j}=1/(2(n-1))$ for all $i\neq j$. Hence, in this case $\lambda_{2}(L_{\mathbf{M}_{1/4}})=n/(2(n-1))$, and, from (15), it suffices that the number of observations $m$ is such that

[TABLE]

Proof of Theorem 4

We now go on to provide a proof of Theorem 4. For pair comparisons, the log-likelihood function (3) can be written as

[TABLE]

where $y_{t}$ denotes the choice from the comparison pair $S_{t}=\{y_{t},z_{t}\}$ for each observation $t$. The negative log-likelihood function satisfies the relation in Lemma 1, which combined with the following two lemmas, establishes the statement of the theorem.

Lemma 5

The following relation holds:

[TABLE]

**Proof ** It is easy to verify that for all $\theta\in\mathbf{R}^{n}$ and $i,j\in N$, we have the following identities

[TABLE]

For all $i,j\in N$ such that $i\neq j$, we have

[TABLE]

and

[TABLE]

Combining with condition P1, we have

[TABLE]

It follows that

[TABLE]

where for two matrices $\mathbf{A}$ and $\mathbf{B}$, $\mathbf{A}\succeq\mathbf{B}$ is equivalent to saying that $\mathbf{A}-\mathbf{B}$ is positive definite; see Appendix A. In (16), both $\nabla^{2}(-\ell(\theta))$ and $4A\frac{m}{n}L_{\mathbf{M}_{1/4}}$ are positive-definite matrices. Hence, by the elementary fact stated in Lemma 23 (Appendix), we obtain the assertion of the lemma.

Lemma 6

With probability at least $1-2/n$,

[TABLE]

**Proof ** $\nabla(-\ell(\theta))$ is a sum of independent random vectors in $\mathbf{R}^{n}$ given by

[TABLE]

The elements of $\nabla(-\log(p_{2}(\theta_{y_{t}}-\theta_{z_{t}})))$ can be expressed as follows

[TABLE]

If $i\notin\{y_{t},z_{t}\}$, then clearly

[TABLE]

Otherwise, we have

[TABLE]

where the last equation is by the fact that $dp_{2}(x)/dx$ is an even function.

By condition P2, we have

[TABLE]

Hence, by the vector Azuma-Hoeffding bound in Lemma 2, we have

[TABLE]

3.2 Comparison Sets of Two or More Items

We now consider a more general case were each comparison set consists of two or more items. We first show an upper bound for the mean squared error of the maximum likelihood parameter estimator when the choices are according to the Luce choice model and comparison sets are of identical sizes. We then present a similar characterization under more general assumption that allow for a broader set of Thurstone choice models and non-identical sizes of comparison sets.

Theorem 7

Suppose that choices are according to the Luce choice model, all comparison sets are of cardinality $k\geq 2$, and $\lambda_{2}(L_{\mathbf{M}})>0$. Then, with probability at least $1-2/n$,

[TABLE]

where $D=4k^{2}e^{4b}$.

The proof of Theorem 7 is provided in Appendix C. The mean squared error upper bound in Theorem 7 corresponds that in Theorem 4 up to a constant factor. If the comparisons sets are unbiased, from (11), we have that $\lambda_{2}(L_{\mathbf{M}})=(1-1/n)k(k-1)$. Hence, the mean squared error upper bound in Theorem 7 depends on $k$ only through the factor $1/(1-1/k)^{2}$, which decreases to value $1$ with $k$ in a diminishing returns fashion. This suggests that there is a limited dependence of the mean squared error on the size of comparison sets.

We now go on to establish a mean-squared upper bound for a class of Thurstone choice models. We will allow for comparison sets of different cardinalities taking values in a set $K$. We will admit the following conditions:

• A1

There exists $A>0$ such that for all $S\subseteq N$ with $|S|\in K$, all $y\in S$, all $i,j\in S$ with $i\neq j$, and all $\theta\in[-b,b]^{n}$,

[TABLE]

and, moreover, the following holds

[TABLE]

• A2

There exists $B>0$ such that for all $S\subseteq N$ with $|S|\in K$, all $y\in S$, and all $\theta\in[-b,b]^{n}$,

[TABLE]

• A3

There exists $C>0$ such that for all $S\subseteq N$ with $|S|\in K$, all $y\in S$, and all $\theta\in[-b,b]^{n}$,

[TABLE]

Condition A1 ensures that $\nabla^{2}(-\log(p_{y,S}(\mathbf{0})))$ is a Laplacian matrix with non-negative weights, and that $\nabla^{2}(-\log(p_{y,S}(\mathbf{\theta})))\succeq A\nabla^{2}(-\log(p_{y,S}(\mathbf{0})))$ for all $\theta\in[-b,b]^{n}$. Condition A1 also ensures that the expected value of $\nabla^{2}(-\log(p_{y,S}(\mathbf{\theta})))$ is a positive definite matrix where $\nabla^{2}(-\log(p_{y,S}(\mathbf{0})))$ is a random matrix when $|S|>2$. Condition A2 requires that $\|\nabla p_{y,S}(\theta)\|_{2}$ is bounded for all $\theta\in[-b,b]^{n}$, while condition A3 ensures that the choice probabilities are not too much imbalanced. Conditions A1 and A2 may be seen as generalizations of conditions P1 and P2 for the case of pair comparisons.

Conditions A1, A2 and A3 can be easily shown to hold for the Luce choice model. For the Luce choice model, we have

[TABLE]

hence, A1 holds with $A=e^{-4b/\beta}$ and $\partial^{2}(-\log(p_{y,S}(\mathbf{0})))/\partial\theta_{i}\partial\theta_{j}=-1/(|S|\beta)^{2}$. Conditions A2 and A3 hold with $B=4$ and $C=e^{-2b/\beta}$.

Note that constants $A$, $B$ and $C$ that appear in A1, A2 and A3, respectively, may depend on $F$, the cardinalities of comparison sets, and the parameter $b$, but are independent of any other parameters. In particular, these constants are independent of the number of observations. For any Thurstone choice model, the constants $A$, $B$, and $C$ can be taken to have values arbitrarily near to value $1$ by taking $b$ small enough.

We next show an upper bound for the mean squared error of the maximum likelihood parameter estimator for a class of Thurstone choice models that satisfy the above stated conditions. Before we do that, we need to introduce some new definitions.

Definition 8** (weight function)**

Let $w^{*}$ be a function defined on positive integers greater than or equal to $2$, we refer to as a weight function, which is defined by

[TABLE]

where

[TABLE]

Notice that for the Luce choice model, $\partial p_{k}(\mathbf{0})/\partial x_{1}=1/(\beta k)^{2}$. Hence, in this case, $w^{*}(k)=1/(\beta k)^{2}$. We will see later in Section 6 that for a broad class of Thurstone choice models, which includes well-known cases with noise according to either Gaussian or double-exponential distribution, $\partial p_{k}(\mathbf{0})/\partial x_{1}=\Theta(1/k^{2})$ and, hence, $w^{*}(k)=\Theta(1/k^{2})$.

Definition 9** ($\gamma_{F,k}$ parameter)**

Let $\gamma_{F,k}$ be a parameter defined by

[TABLE]

We note that for any comparison set $S$ of cardinality $k$ and all $y\in S$, we have

[TABLE]

which is discussed in more detail in the proof of Lemma 30. In particular, for the Luce choice model, we have $\gamma_{F,k}=(1-1/k)/\beta^{2}$.

Theorem 10

Assume A1, A2 and A3, let $\sigma_{F,K}$ be such that $1/\gamma_{F,k}\leq\sigma_{F,K}$ for all $k\in K$, and $\lambda_{2}(L_{\overline{\mathbf{M}}_{w^{\star}}})\geq 32(\sigma_{F,K}/C)n\log(n)/m$. Then, with probability at least $1-3/n$,

[TABLE]

where $D=B/(AC)$.

The proof of the theorem is provided in Appendix D. The main technical difference of the proof with respect to that of Theorem 7 is that $\nabla^{2}(-\ell(\theta))$ is a sum of random matrices. Every $\nabla^{2}(-\log(p_{y_{t},S_{t}}(\theta)))$ is a random matrix for the following two reasons: (a) $\nabla^{2}(-\log(p_{y_{t},S_{t}}(\theta)))$ depends on the randomly chosen $y_{t}$ and (b) $S_{t}$ is allowed to be a random set of items. We use the matrix Chernoff bound in the proof of Theorem 10.

We next show two corollaries of Theorem 10, which cover two interesting special cases.

Corollary 11

Suppose that all comparison sets are of identical cardinality of value $k\geq 2$, A1, A2, A3 hold, and $\lambda_{2}(L_{\overline{\mathbf{M}}_{1/k^{2}}})\geq 32(k-1)/(Ck)$. Then, with probability at least $1-3/n$,

[TABLE]

where $\gamma_{F,k}$ is given by (19), and $\mathbf{M}_{1/k^{2}}$ is the weighted-adjacency matrix with the weight function $w(k)=1/k^{2}$.

Corollary 12

Suppose that comparison sets are independent with each comparison set being a sample without replacement from the set of all items, conditions A1, A2, A3 hold, and $m\geq 32(1-1/n)/Cn\log(n)$. Then, with probability at least $1-3/n$,

[TABLE]

3.3 Lower Bound

In this section, we present a lower bound for the mean squared error of the maximum likelihood parameter estimator, which establishes minimax optimality of the established upper bounds. We define the following conditions:

• A1’

There exists $\widetilde{A}>0$ such that for all $S\subseteq N$ with $|S|\in K$, all $y\in S$, all $i,j\in S$ such that $i\neq j$, and all $\theta\in[-b,b]^{n}$, it holds

[TABLE]

• A3’

There exists $\widetilde{C}>0$ such that for all $S\subseteq N$ with $|S|\in K$, all $y\in S$, and all $\theta\in[-b,b]^{n}$, it holds

[TABLE]

Notice that, in particular, for the special case of $F$ being a double-exponential distribution with parameter $\beta$, we have that A1’ and A3’ hold with $\widetilde{A}=e^{4b/\beta}$ and $\widetilde{C}=e^{2b/\beta}$.

Theorem 13

Under conditions A1’ and A3’, for any unbiased estimator $\widehat{\theta}$, we have

[TABLE]

The following two corollaries follow from the last theorem.

Corollary 14

If all comparison sets are of cardinality $k\geq 2$, then any unbiased estimator $\widehat{\theta}$ satisfies

[TABLE]

Corollary 15

If, in addition, each comparison set is drawn independently, uniformly at random from the set of all items, then any unbiased estimator $\widehat{\theta}$ satisfies

[TABLE]

The last corollary implies that under the given assumptions, for the mean squared error to be smaller than a constant, it is necessary that the number of observed comparisons is $m=\Omega(\gamma_{F,k}n)$.

Proof of Theorem 13 We denote by $cov[Y]$ the covariance matrix of a multivariate random variable $Y$ i.e.,

[TABLE]

The proof uses the Cramér-Rao inequality, which is stated as follows.

Lemma 16** (Cramér-Rao bound)**

Let $X$ be a multivariate random variable with distribution $p(x;\theta)$, for parameter $\theta\in\Theta$, and let $\psi:\Theta\rightarrow\mathbf{R}^{r}$ be a differentiable function. Then, for any unbiased estimator $\mathbf{T}(X)=(T_{1}(X),\dots,T_{r}(X))^{\top}$ of $\mathbf{\psi}(\theta)=(\psi_{1}(\theta),\dots,\psi_{r}(\theta))^{\top}$, we have

[TABLE]

where $\frac{\partial\mathbf{\psi}(\theta)}{\partial\theta}$ is the Jacobian matrix of $\psi$ and $F(\theta)$ is the Fisher information matrix given by

[TABLE]

Let us define $\psi_{i}(\theta)=\theta_{i}-\frac{1}{n}\sum_{l=1}^{n}\theta_{l}$ for all $i=1,2,\ldots,n$. Since $\sum_{i=1}^{n}\theta_{i}=0$, we have $\sum_{i=1}^{n}\psi_{i}(\theta)=0$. Note that we can write

[TABLE]

Let $F(\theta)$ be the Fisher information matrix given by

[TABLE]

By conditions A1’ and A3’, and Lemma 36, we have

[TABLE]

where $\mathbf{M}_{S}$ is a matrix in $\mathbf{R}^{n\times n}$ that has each element $(i,j)$ such $\{i,j\}\subseteq S$ equal to $1$ and all other elements equal to [math].

From (21) and (25),

[TABLE]

Note that

[TABLE]

By the Cramér-Rao bound and (26), we have

[TABLE]

4 Rank-Breaking Parameter Estimation Methods

The maximum likelihood parameter estimation requires to find a maximizer of a log-likelihood function which for comparison sets of cardinality $k$ has the form of a sum of $\binom{n}{k}$ elements in the worst case. For pair comparisons, this sum consists of at most $n^{2}$ elements. It is common for parameter estimators to be defined as maximizers of a pseudo log-likelihood function, which is defined as the log-likelihood function of pair comparisons deduced from the input observations under assumption that these pair comparisons are independent (which in general is false under a Thurstone choice model for comparison sets of three or more items). This is commonly referred as rank breaking. In what follows, we consider two different rank-breaking methods: (a) one that deduces $k-1$ pair comparisons from a choice from a comparison set of cardinality $k$, we refer to as rank-breaking method $\mathsf{ALL}$, and (b) one that deduces $1$ pair comparison from a choice from a comparison set of cardinality $k$, we refer to as rank-breaking method $\mathsf{ONE}$.

Rank-breaking method $\mathsf{ALL}$

This rank-breaking method deduces $k-1$ pair comparisons from each comparison set of cardinality $k$. Specifically, for every comparison set $S_{t}$, the method uses pair comparisons between the chosen item $y_{t}$ and each non-chosen item $v\in S_{t}\setminus\{y_{t}\}$. Notice that for any comparison set of three or more items, the pair comparisons selected from this set by the given rank-breaking method are not independent.

For the Luce choice model, the pseudo log-likelihood function is given by (6), which can be written in a more explicit form as follows

[TABLE]

We consider the maximum pseudo log-likelihood estimator $\widehat{\theta}_{k-1}:=\arg\max_{\theta\in\Theta}\ell_{k-1}(\theta)$.

Theorem 17

If $\lambda_{2}\left(L_{\mathbf{M}}\right)\geq 128(k-1)^{2}e^{2b}n\log(n)/m$, then with probability at least $1-3/n$,

[TABLE]

where $D=16\sqrt{2}\sqrt{k(k-1)^{3}}e^{2b}$.

The proof of Theorem 17 is given in Appendix E.

The mean squared error upper bound in Theorem 17 implies that the given parameter estimator is consistent. For any fixed size of a comparison set, the bound in Theorem 17 is equal to that in Theorem 7 up to a constant factor. Both bounds have the same scaling with parameter $k$.

Rank-breaking method $\mathsf{ONE}$

This rank-breaking method deduces $1$ pair comparison from a comparison set of cardinality $k$. From each comparison set $S_{t}$, this rank-breaking methods selects a pair that consists of the chosen item $y_{t}$ and an item $z_{t}$ selected uniformly at random from the set of non-chosen items $S_{t}\setminus\{y_{t}\}$.

For the Luce choice model, the pseudo log-likelihood function is given by (7), which can be written in a more explicit form as follows

[TABLE]

We consider the maximum pseudo log-likelihood estimator $\widehat{\theta}_{1}:=\arg\max_{\theta\in\Theta}\ell_{1}(\theta)$.

Theorem 18

If $\lambda_{2}\left(L_{\mathbf{M}}\right)\geq 8k(k-1)e^{2b}n\log(n)/m$, then with probability at least $1-3/n$,

[TABLE]

where $D=4k(k-1)e^{2b}$.

The proof of Theorem 18 is given in Appendix F.

It is noteworthy that the mean squared error upper bounds in Theorem 17 and Theorem 18 are equal up to a constant factor. Intuitively, one would expect that rank-breaking method $\mathsf{ALL}$ would yield a smaller mean squared error than rank-breaking method $\mathsf{ONE}$ because it uses more information from each observed choice. The reason why the two mean squared error upper bounds are equal up to a constant factor is as follows. When applying Lemma 1 we need to find an upper bound $\alpha$ for the norm of the gradient of the negative pseudo log-likelihood function and a lower bound $\beta$ for the second-smallest eigenvalue of the Hessian matrix of the negative pseudo log-likelihood function. In our proofs, for the case of Theorem 18, we obtained $\alpha$ and $\beta$ that scale with $k$ as $1$ and $1/k^{2}$, respectively. On the other hand, for the case of Theorem 17, we obtained $\alpha$ and $\beta$ that scale with $k$ as $k$ and $1/k$, respectively. For both cases, it follows that the ratio $\alpha/\beta$ scales as $k^{2}$.

5 Binary Classification

We now consider a Thurstone choice model ${\mathcal{T}}_{F}$ where the strength parameter vector $\theta$ takes value in $\Theta=\{-b,b\}^{n}$, for a parameter $b>0$. Here each individual strength parameter takes either a low or a high value. We say that each item is either of a low or a high class. We consider a binary classification problem, where the goal is to correctly classify all items with a prescribed probability of classification error. Let $N_{1}$ be the set of high class items and $N_{2}$ be the set of low class items. We shall consider the case when the total number of items is even and the number of high class items is equal to the number of low class items.

We consider a simple classification algorithm that uses point scores defined as follows: each item is associated with a point score equal to the number of comparison sets in which the given item is the chosen item. The algorithm outputs a classification of items with $\widehat{N}_{1}$ and $\widehat{N}_{2}$ denoting the sets of items classified to be high class or low class, respectively. The algorithm follows the following three steps: (a) for each item compute its point score, (b) sort the items in decreasing order of point scores, and (c) let $\widehat{N}_{1}$ contain $n/2$ items with highest point scores and $\widehat{N}_{2}$ contain remaining items. We refer to this algorithm as a point score ranking method.

Theorem 19

Suppose that $b\leq 4/(k^{2}\partial p_{k}(\mathbf{0})/\partial x_{1})$ and

[TABLE]

Then, for every $\delta\in(0,1]$, if

[TABLE]

then, the point score ranking method correctly classifies all items with probability at least $1-\delta$.

The proof of Theorem 19 is provided in Appendix J.

The sufficient condition in (28) for the point score ranking method to correctly classify all the items with probability at least $1-\delta$ is shown to be necessary up to a constant-factor for any classification algorithm, which is given in the following theorem.

Theorem 20

Suppose that $b\leq 1/(6k^{2}\partial p_{k}(\mathbf{0})/\partial x_{1})$ and that condition (27) holds. Then, for every even $n\geq 16$ and $\delta\in(0,1/4]$, for any algorithm to correctly classify all the items with probability at least $1-\delta$, it is necessary that the following condition holds

[TABLE]

The proof of Theorem 20 is given in Appendix K. In the proof, we use the statistical difference between the case when all the items are correctly classified and the case that an item is incorrectly classified. This proof strategy is motivated by that in Yun and Proutiere (2014) where it was used to analyze the classification error of the stochastic block model.

6 Discussion of Results

In this section, we discuss how the number of observations needed to attain a prescribed parameter estimation error depends on the cardinality of comparison sets.

In Section 3, we found that for a priori unbiased input comparisons, where each comparison set is of cardinality $k$ and is drawn uniformly at random from the set of all items, the number of observations needed for the mean squared error to be within a prescribed tolerance is of the order $\gamma_{F,k}$, defined by (19). In Section 5, we found that this also so to ensure that the probability of classification error is within a prescribed tolerance.

In Table 1, we show the values of the parameter $\gamma_{F,k}$ for several special instances of the Thurstone choice model, along with the values of $\partial p_{k}(\mathbf{0})/\partial x_{1}$. From the expressions in Table 1, we observe that for all the cases but the case of uniform distribution of noise, $\gamma_{F,k}$ decreases with the cardinality of comparison sets, but in a slow manner according to a diminishing returns relation. In particular, observe that for both double-exponential and Laplace distribution of noise, $\gamma_{F,k}=\Theta(1)$ and for Gaussian distribution of noise $\gamma_{F,k}=O(1/k^{\epsilon})$. On the other hand, for uniform distribution of noise, $\gamma_{F,k}=\Theta(1/k^{2})$.

It is noteworthy that $\gamma_{F,k}$ satisfies the following bounds.

Lemma 21

For any cumulative distribution function $F$ with density function $f$:

If $f$ is even and continuously differentiable, then $\gamma_{F,k}=O(1)$. 2. 2.

If $f$ is such that $f(x)\leq C$ for all $x\in\mathbf{R}$, for a constant $C>0$, then $\gamma_{F,k}=\Omega(1/k^{2})$.

We observe that both double-exponential distribution and Laplace distribution of noise are extremal in the sense that they achieve the upper bound $\gamma_{F,k}=O(1)$. On the other hand, uniform distribution of noise is extremal in the sense that it achieves the lower bound $\gamma_{F,k}=\Omega(1/k^{2})$.

7 Experimental Results

In this section, we present our experimental results using both simulations and real-world data. Our first goal is to provide experimental validation of the claim that the mean squared error can depend on the cardinality of comparison sets in different ways for different Thurstone choice models, which is suggested by our theory. Our second goal is to evaluate Fiedler value for different weighted-adjacency matrices observed in real-world data, which demonstrates that it can assume a wide range of values depending on the application scenario.

7.1 MSE versus Cardinality of Comparison Sets

We consider the following simulation experiment. We fix the number of items $n$ and the number of observations $m$. We then run experiments for different values of the cardinality of comparison sets $k$. For each given value of parameter $k$, we generate comparison sets as independent uniform random sets of cardinality $k$ from the set of all items. We then draw choices according to a Thurstone choice model ${\mathcal{T}}_{F}$ for the value of parameter vector $\theta^{\star}=\mathbf{0}$. For every fixed value of $k$, we run $100$ repetitions to estimate the mean squared error. We do this for the distribution of noise according to a double-exponential distribution (Bradley-Terry model) and according to a uniform distribution, both with unit variance.

Figure 1 shows the results for the case of $n=10$ and $m=100$. The results clearly demonstrate that the mean squared error exhibits qualitatively different decay with the cardinality of comparison sets for the two Thurstone choice models under consideration. Our theoretical results in Section 3.2 suggest that the mean squared error should decrease with the cardinality of comparison sets as $1/(1-1/k)$ for the double-exponential distribution, and as $1/k^{2}$ for the uniform distribution of noise. Observe that the latter two terms decrease with $k$ to a strictly positive value and to zero value, respectively. The empirical results in Figure 1 confirm these claims.

7.2 Fiedler Values of Weighted-Adjacency Matrices

We found that Fiedler value of a weighted-adjacency matrix plays a key role in upper bounds on the mean squared error of parameter estimator in Section 3.1 and Section 3.2. Here we evaluate Fiedler value for different weighted-adjacency matrices of different schedules of comparisons. Throughout this section, we use the definition of a weighted-adjacency matrix in (9) with the weight function $w(k)=1/k^{2}$. Our first two examples are representative of schedules in sport competitions, which are typically carefully designed by sport associations and exhibit a large degree of regularity. Our second two examples are representative of comparisons that are induced by user choices in the context of online services, which exhibit much more irregularity.

Sport competitions

We consider the fixtures of games for the season 2014-2015 for (i) football Barclays premier league and (ii) basketball NBA league. In the Barclays premier league, there are 20 teams, each team plays a home and an away game with each other team; thus there are 380 games in total. In the NBA league, there are 30 teams, 1,230 regular games, and 81 playoff games.111The NBA league consists of two conferences, each with three divisions, and the fixture of games has to obey constraints on the number of games played between teams from different divisions. We evaluate Fiedler value of weighted-adjacency matrices defined for first $m$ matches of each season; see Figure 2.

For the Barclays premier league dataset, at the end of the season, the Fiedler value of the weighted-adjacency matrix is of value $n/[2(n-1)]\approx 1/2$. The schedule of matches is such that at the middle of the season, each team played against each other team exactly once, at which point the Fiedler value is $n/[4(n-1)]\approx 1/4$. The Fiedler value is of a strictly positive value after the first round of matches. For most part of the season, its value is near to $1/4$ and it grows to the highest value of approximately $1/2$ in the last round of the matches.

For the NBA league dataset, at the end of the season, the Fiedler value of the weighted-adjacency matrix is approximately $0.375$. It grows more slowly with the number of games played than for the Barclays premier league; this is intuitive as the schedule of games is more irregular, with each team not playing against each other team the same number of times.

Crowdsourcing contests

We consider participation of users in contests of two competition-based online labour platforms: (i) online platform for software development TopCoder and (ii) online platform for various kinds of tasks Tackcn. We refer to coders in TopCoder and workers in Taskcn as users. We consider contests of different categories observed in year 2012. In both these systems, the participation in contests is according to choices made by users.

For each set of tasks of given category, we conduct the following analysis. We consider a conditioned dataset that consists only of a set of top-$n$ users with respect to the number of contests they participated in given year, and of all contests attended by at least two users from this set. We then evaluate Fiedler value of the weighted-adjacency matrix for parameter $n$ ranging from $2$ to the smaller of $100$ or the total number of users. Our analysis reveals that the Fiedler value tends to decrease with $n$. This indicates that the larger the number of users included, the less connected the weighted-adjacency matrix is. See the top plots in Figure 3.

We also evaluated the smallest number of contests from the beginning of the year that is needed for the Fiedler value of the weighted-adjacency matrix to assume a strictly positive value. See the bottom plots in Figure 3. We observe that this threshold number of contests tends to increase with the number of top users considered. There are instances for which this threshold substantially increases for some number of the top users. This, again, indicates that the algebraic connectivity of the weighted-adjacency matrices tends to decrease with the number of top users considered.

8 Conclusion

The results of this paper elucidate how the estimation accuracy of a Thurstone choice model parameter depends on the given model and the structure of comparison sets. They show that a key factor is the algebraic connectivity of a weighted-adjacency matrix, which is specific to given model. It is shown that for a large class of Thurstone choice models, including well-known instances, there is a diminishing returns decrease of the estimation error with the cardinality of comparison sets at a slow rate for comparison sets of three of more items.

The results provide guidelines to the designers of competition schedules such as to ensure that a schedule has a well-connected weighted-adjacency matrix and to expect limited estimation accuracy gains by enlarging the size of comparison sets.

A Background Material

Location of Eigenvalues

We make note of the well-known Geršgorin circles theorem, which we state as at the following lemma:

Lemma 22

Let $\mathbf{A}\in\mathbf{R}^{n\times n}$, then all eigenvalues of $\mathbf{A}$ are located in the union of $n$ discs

[TABLE]

Properties of Positive Definite and Laplacian Matrices

A symmetric matrix $\mathbf{A}\in\mathbf{R}^{n\times n}$ is said to be positive semidefinite if $\mathbf{x}^{\top}\mathbf{A}\mathbf{x}\geq 0$ for all nonzero $\mathbf{x}\in\mathbf{R}^{n}$. If the inequality is replaced with strict inequality, $\mathbf{A}$ is said to be positive definite.

Each eigenvalue of a positive definite matrix is a positive real number. Each eigenvalue of a positive semidefinite matrix is a nonnegative real number.

For two matrices $\mathbf{A},\mathbf{B}\in\mathbf{R}^{n\times n}$, we write $\mathbf{A}\succeq\mathbf{B}$ for the positive semidefinite ordering, which means that $\mathbf{A}-\mathbf{B}$ is a positive semidefinite matrix. Similarly, we write $\mathbf{A}\succ\mathbf{B}$ for the positive definite ordering, which means that $\mathbf{A}-\mathbf{B}$ is a positive definite matrix. Note that $\mathbf{A}\succeq\mathbf{0}$ means that $\mathbf{A}$ is a positive semidefinite matrix, and $\mathbf{A}\succ\mathbf{0}$ means that $\mathbf{A}$ is a positive definite matrix.

We note the following ordering relations for eigenvalues of two positive definite matrices that satisfy the positive semidefinite ordering (e.g., see Corollary 7.7.4 Horn and Johnson (1985)).

Lemma 23

For any two positive definite matrices $\mathbf{A},\mathbf{B}\in\mathbf{R}^{n\times n}$ such that $\mathbf{A}\succeq\mathbf{B}$, $\lambda_{i}(\mathbf{A})\geq\lambda_{i}(\mathbf{B})$ for all $i=1,2,\ldots,n$.

For any $\mathbf{A}\in\mathbf{R}^{n\times n}$ with zero diagonal such that $\mathbf{A}\geq\mathbf{0}$, the Laplacian matrix $L_{\mathbf{A}}$ is positive semidefinite. This follows from the localization of eigenvalues by the Geršgorin circles theorem.

Lemma 24

If $\mathbf{A}\in\mathbf{R}^{n}$ is a symmetric matrix with zero diagonal, then

1. (i)

$\lambda_{1}(L_{\mathbf{A}})=0$, and 2. (ii)

$\lambda_{i}(L_{\mathbf{A}})=\lambda_{i-1}(\mathbf{U}^{\top}L_{\mathbf{A}}\mathbf{U})$, for $i\in\{2,3,\ldots,n\}$

where $\mathbf{U}\in\mathbf{R}^{n\times n-1}$ is any matrix such that $\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}$.

We shall use the following property of real symmetric matrices:

Lemma 25

If $\mathbf{A},\mathbf{B}\in\mathbf{R}^{n\times n}$ are real symmetric matrices with zero diagonals such that $\mathbf{B}\geq\mathbf{A}$ where the inequality holds element-wise, then $L_{\mathbf{B}}\succeq L_{\mathbf{A}}$.

Chernoff Tail Bounds

The following bounds follow from the Chernoff bound:

Lemma 26

Suppose that $X$ is a sum of $m$ independent Bernoulli random variables each with mean $p$, then if $q\leq p\leq 2q$,

[TABLE]

and, if $p\leq q$,

[TABLE]

**Proof ** We prove only prove (30) as (29) follows by similar arguments. By the Chernoff’s bound, for every $s>0$,

[TABLE]

where $h(s)=qs-\log(1-p+pe^{s})$. Since $\log(1-x)\leq-x$ for all $x\in\mathbf{R}$, we have $h(s)\geq qs+p-pe^{s}$. Take $s=s^{*}:=\log(q/p)$ to obtain $h(s^{*})\geq q\log(q/p)+p-q$.

Now, let $\epsilon=q-p$, and note that $q\log(q/p)+p-q:=g(\epsilon)$ where $g(\epsilon)=q\log(q/(q-\epsilon)-\epsilon$. Since $g^{\prime}(\epsilon)=q/(q-\epsilon)-1=\epsilon/(q-\epsilon)\geq\epsilon/(2q)$, we have $g(\epsilon)=\int_{0}^{\epsilon}g^{\prime}(x)dx\geq\epsilon^{2}/(4q)$.

Hence, it follows that $h(s^{*})\geq(p-q)^{2}/(4q)$, and, thus

[TABLE]

B Proof of Lemma 1

Let $\Delta=\widehat{\mathbf{x}}-\mathbf{x}^{\star}$. By the Taylor expansion, we have

[TABLE]

Since $g(\widehat{\mathbf{x}})\leq g(\mathbf{x}^{\star})$, we have

[TABLE]

Hence,

[TABLE]

Fix an arbitrary $\mathbf{x}\in{\mathcal{X}}$. From condition (i) it follows that $\nabla^{2}g(\mathbf{x})$ has eigenvalue [math] with eigenvector $\mathbf{1}$. Combining with condition (ii), we have

[TABLE]

Let $\mathbf{U}=[\mathbf{u}_{1},\mathbf{u}_{2},\ldots,\mathbf{u}_{n}]\in\mathbf{R}^{n\times n}$ where $\mathbf{u}_{1},\mathbf{u}_{2},\ldots,\mathbf{u}_{n}$ are ortonormal eigenvectors of $\nabla^{2}g(\mathbf{x})$, which correspond to eigenvalues $\lambda_{1}(\nabla^{2}g(\mathbf{x})),\lambda_{2}(\nabla^{2}g(\mathbf{x})),\ldots,\lambda_{n}(\nabla^{2}g(\mathbf{x}))$, respectively. Note that

[TABLE]

Let $\Lambda=\mathrm{diag}(\lambda_{1}(\nabla^{2}g(\mathbf{x})),\lambda_{2}(\nabla^{2}g(\mathbf{x})),\ldots,\lambda_{n}(\nabla^{2}g(\mathbf{x})))$.

We have the following relations:

[TABLE]

where we use the properties in (33) and (34).

Hence, it follows that

[TABLE]

By combining (35) and (32), we conclude the proof of the lemma.

C Proof of Theorem 7

The log-likelihood function in (3) can be written in the following more explicit form:

[TABLE]

Since $\nabla^{2}(-\ell(\theta))\mathbf{1}=\mathbf{0}$, for all $\theta\in\mathbf{R}^{n}$, under assumption that $\min_{\theta\in\Theta}\lambda_{2}\left(-\nabla^{2}\ell(\theta)\right)>0$, the upper bound in Lemma 1 holds. This combined with the following two lemmas yields the statement of the theorem.

Lemma 27

The following lower bound holds:

[TABLE]

**Proof ** From (36), we have for $i,j\in\{1,2,\ldots,n\}$ such that $j\neq i$,

[TABLE]

Since for all $i\neq j$,

[TABLE]

we have that $\nabla^{2}(-\ell(\theta))-\frac{1}{k^{2}e^{4b}}\frac{m}{n}L_{\mathbf{M}}$ is a Laplacian matrix of a matrix with nonnegative elements. Every such Laplacian matrix is positive semidefinite, hence

[TABLE]

From (38), we conclude (37).

Lemma 28

With probability at least $1-2/n$,

[TABLE]

**Proof ** From (36), we have

[TABLE]

where $p_{y_{t},S_{t}}(\theta)=e^{\theta^{\star}_{y_{t}}}/\sum_{v\in S_{t}}e^{\theta_{v}}$.

It is straightforward to derive that $\nabla(-\log(p_{y_{t},S_{t}}(\theta)))$ has elements given by

[TABLE]

From (40), $\nabla\ell(\theta)$ is a sum of independent random vectors in $\mathbf{R}^{n}$ that satisfy for all $t\in\{1,2,\ldots,m\}$,

[TABLE]

and

[TABLE]

The last two relations are easy to establish using (41) as follows. Equation (42) holds because for every $i\in N$,

[TABLE]

Equation (42) follows from

[TABLE]

By the vector Azuma-Hoeffding bound in Lemma 2, we have

[TABLE]

which completes the proof of (39).

D Proof of Theorem 10

Since $\nabla^{2}(-\ell(\theta))$ is a Laplacian matrix, by condition A1,

[TABLE]

Hence, in particular,

[TABLE]

We have the following two lemmas.

Lemma 29

If $\lambda_{2}(L_{\overline{\mathbf{M}}_{w^{\star}}})\geq 32(\sigma_{F,K}/C)(n\log(n))/m$, then with probability at least $1-1/n$,

[TABLE]

**Proof ** $\nabla^{2}(-\ell({\bm{0}}))$ is a sum of independent random matrices given by

[TABLE]

We have the following relations:

[TABLE]

where we use Lemma 36.

Hence,

[TABLE]

By Lemma 37, for all $t\in\{1,2,\ldots,m\}$,

[TABLE]

Using the matrix Chernoff bound with $\epsilon=1/2$, (48) and (49, we conclude the statement of the lemma.

Lemma 30

With probability at least $1-2/n$,

[TABLE]

**Proof ** For every $S\subseteq N$ and $i,j\in S$ such that $j\neq i$,

[TABLE]

and

[TABLE]

From (50) and (51), for every $t\in\{1,2,\ldots,m\}$,

[TABLE]

From (50) and (51), for every $S\subseteq N$ such that $|S|=k\geq 2$ and $i,j\in S$ such that $i\neq j$,

[TABLE]

[TABLE]

Hence, for every $t\in\{1,2,\ldots,m\}$ such that $|S_{t}|=k$

[TABLE]

By condition A3 and (53), for every $t\in\{1,2,\ldots,m\}$,

[TABLE]

Using (52) and (54) and the vector Azuma-Hoeffding bound in Lemma 2, with probability at least $1-2/n$,

[TABLE]

The negative log-likelihood function satisfies the bound in Lemma 1. This combined with relation (44), Lemma 29 and Lemma 30 implies that if $\lambda_{2}(L_{\overline{\mathbf{M}}_{w^{\star}}})\geq 32(\sigma_{F,K}/C)n\log(n)/m$, then with probability at least $1-3/n$,

[TABLE]

which proves the theorem.

E Proof of Theorem 17

Under condition that $\lambda_{2}\left(\nabla^{2}(-\ell_{k-1}(\theta)\right)>0$, the negative pseudo log-likelihood function satisfies the bound in Lemma 1. This, combined with the following two lemmas implies the statement of the theorem.

Lemma 31

If $\lambda_{2}(L_{\mathbf{M}})\geq 8k(k-1)e^{2b}n\log(n)/m$, then with probability at least $1-1/n$,

[TABLE]

**Proof ** We will establish the lemma by using the matrix Chernoff bound in Lemma 3 as follows. Note that $\nabla^{2}(-\ell_{k-1}(\theta))$ is a sum independent random matrices given by:

[TABLE]

The nonzero elements of $\nabla^{2}\left(\sum_{v\in S_{t}\setminus\{y_{t}\}}-\log(p_{y_{t},v}(\theta))\right)$ are, for every $i,j\in S_{t}$ such that $i\neq j$,

[TABLE]

and

[TABLE]

From (55), we have

[TABLE]

Hence, we have

[TABLE]

and, in particular,

[TABLE]

To apply the matrix Chernoff bound in Lemma 3, we use the following identities that follow by Lemma 24,

[TABLE]

and

[TABLE]

and the following fact:

[TABLE]

where the first equation follows by $\|\mathbf{U}^{\top}\|_{2}=\|\mathbf{U}\|_{2}=1$, and the second inequality is by the Geršgorin circles theorem (Lemma 22).

[TABLE]

where the first inequality is by (56), the equality is by (57) and (58), the second inequality is by the matrix Chernoff bound in Lemma 3 and (61), the third inequality is by (56), and the last inequality is by the condition $\lambda_{2}(L_{\mathbf{M}})\geq 8k(k-1)e^{2b}n\log(n)/m$.

Lemma 32

With probability at least $1-2/n$,

[TABLE]

**Proof ** $\nabla(-\ell_{k-1}(\theta))$ is a sum of independent random vectors given by:

[TABLE]

It is straightforward to show that $\nabla\left(\sum_{v\in S_{t}\setminus\{y_{t}\}}-\log(p_{y_{t},v}(\theta))\right)$ has the elements given by

[TABLE]

For every $t\in\{1,2,\ldots,m\}$ and $i\in\{1,2,\ldots,n\}$, we have

[TABLE]

The last equation obviously holds for every $i\notin S_{t}$, and it holds for $i\in S_{t}$ by the following derivations

[TABLE]

From (63), we have

[TABLE]

Hence,

[TABLE]

Therefore, by vector Azuma-Hoeffding bound in Lemma 2, (62) holds with probability at least $1-2/n$.

F Proof of Theorem 18

The proof follows by the same steps as that of Theorem 17, and the following two lemmas.

Lemma 33

If $\lambda_{2}\left(L_{\mathbf{M}}\right)\geq 8k(k-1)e^{2b}n\log(n)/m$, then with probability at least $1-1/n$,

[TABLE]

**Proof ** $\nabla^{2}(-\ell_{1}(\theta))$ is a sum of random matrices given by:

[TABLE]

where $p_{y_{t},z_{t}}(\theta)=e^{\theta_{y_{t}}}/(e^{\theta_{y_{t}}}+e^{\theta_{z_{t}}})$.

It is easy to establish that for all $t\in\{1,2,\ldots,m\}$ and $i\neq j$,

[TABLE]

and

[TABLE]

Hence, for every $i,j\in\{1,2,\ldots,n\}$ such that $i\neq j$,

[TABLE]

Hence,

[TABLE]

and, in particular,

[TABLE]

By (66) and (67) and Grešgorin circle theorem,

[TABLE]

The rest of the proof follows by the same arguments as in the proof of Lemma 31.

Lemma 34

With probability at least $1-2/n$,

[TABLE]

**Proof ** $\nabla(-\ell_{1}(\theta))$ is a sum of independent random vectors in $\mathbf{R}^{n}$ with elements given by

[TABLE]

It follows that for all $t\in\{1,2,\ldots,m\}$ and $i\in\{1,2,\ldots,n\}$,

[TABLE]

and

[TABLE]

The statement of the lemma than follows by vector Azuma-Hoeffding bound in Lemma 2.

G Proof of Lemma 3

From Theorem 5.1.1 in Tropp (2015),

[TABLE]

which combined with the fact

[TABLE]

yields the lemma.

G.1 Proof of Lemma 24

Since $L_{\mathbf{A}}=\hbox{diag}(\mathbf{A}\mathbf{1})-\mathbf{A}$, we have

[TABLE]

Hence, [math] is an eigenvalue of $L_{\mathbf{A}}$ for eigenvector $\mathbf{1}$.

Let $\mathbf{y}\in\mathbf{R}^{n}$ be an eigenvector with corresponding eigenvalue $\lambda$ of $L_{\mathbf{A}}$. Since the columns of $\mathbf{U}$ are independent, $\mathbf{U}$ is nonsingular and it has inverse $\mathbf{U}^{-1}$ (e.g., Section 0.5 Horn and Johnson (1985). Let $\mathbf{x}=\mathbf{U}^{-1}\mathbf{y}$).

Note that the following equations hold:

[TABLE]

Hence, it follows that $\lambda$ is an eigenvalue of $\mathbf{U}^{\top}L_{\mathbf{A}}\mathbf{U}$ with corresponding eigenvector $\mathbf{x}$.

H Proof of Lemma 25

Let $\mathbf{C}=\mathbf{B}-\mathbf{A}\geq\mathbf{0}$. Note that

[TABLE]

Since $L_{\mathbf{C}}$ is positive semidefinite, it follows that $L_{\mathbf{B}}-L_{\mathbf{A}}$ is positive semidefinite, i.e. $L_{\mathbf{B}}\succeq L_{\mathbf{A}}$.

Lemma 35

For all $S\subseteq\{1,2,\ldots,n\}$ such that $|S|=k\geq 2$ we have: for all $i,j\in\{1,2,\ldots,n\}$ such that $i\neq j$, for $v\in\{1,2,\ldots,n\}\setminus\{i,j\}$,

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

**Proof ** Let $S\subseteq N$ be such that $|S|=k$, for an integer $2\leq k\leq n$. Without loss of generality, let $S=\{1,2,\ldots,k\}$. Let $\mathbf{x}_{v}(\theta)=(\theta_{v}-\theta_{u},u\in S\setminus\{v\})$, for $v\in S$. We first consider $\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}(-\log(p_{k}(\mathbf{x}_{v}(\theta))))$ for $i\neq j$. It is easy to note that

[TABLE]

We separately consider two different cases.

Consider first the case when $\{i\}\cap\{j\}\cap\{v\}=\emptyset$. By differentiation, we have

[TABLE]

which establishes (73).

Consider now the case when $i\neq v$ and $j=v$. First, note

[TABLE]

For every $u\in S$, $p_{k}(\mathbf{x}_{u}(\theta))$ does not change its value by changing $\theta$ with $\theta+c\mathbf{1}$, for every constant $c\in\mathbf{R}$. Hence, by full differentiation, we have

[TABLE]

From (77),

[TABLE]

Now, note that

[TABLE]

Hence,

[TABLE]

From this it follows

[TABLE]

From (78) and (79),

[TABLE]

From (77),

[TABLE]

Combining (76), (80) and (81), we obtain (74).

Since for every $v\in S$, $\log(p_{k}(x_{v}(\theta)))$ does not change its value by changing $\theta$ to $\theta+c\mathbf{1}$ for all $c\in\mathbf{R}$, by full differentiation

[TABLE]

Taking partial derivative with respect to $\theta_{i}$ on both sides implies (75).

Lemma 36

Let $S\subseteq N$ be such that $|S|=k\geq 2$ and $Y$ be a random variable according to distribution $p_{y,S}(\theta)$, for $\theta\in\mathbf{R}^{n}$. Then,

[TABLE]

where $\mathbf{M}_{S}=[m_{i,j}^{S}]\in\mathbf{R}^{n\times n}$ is such that $m_{i,j}^{S}=1$ if $i,j\in S$ and $i\neq j$, and $m_{i,j}^{S}=0$, otherwise.

**Proof ** By (73) and (74) in Lemma 75, we have for all $i,j\in\{1,2,\ldots,n\}$ such that $i\neq j$,

[TABLE]

By (75) in Lemma 75, we have

[TABLE]

Lemma 37

If for $S\subseteq N$ such that $|S|=k\geq 2$,

[TABLE]

then

[TABLE]

where $1/\gamma_{F,k}=\left(k^{2}\partial p_{k}(\mathbf{0})/\partial x_{1}\right)^{2}$.

**Proof ** Without loss of generality, let $y=1$ and $S=\{1,2,\dots,k\}$. By Lemma 75, we have

[TABLE]

and for $i\neq 1$, $j\neq 1$ and $j\neq i$,

[TABLE]

Since by condition of the lemma the left-hand side in (83) is non positive, we have

[TABLE]

and, since by condition of the lemma the left-hand side in (84) is non positive, we have

[TABLE]

From (83) and (86),

[TABLE]

Consider now the case when $i\neq 1$, $j\neq 1$ and $i\neq j$. If $k=2$, then obviously

[TABLE]

Otherwise, if $k>2$, from (84) and (85), we have

[TABLE]

Hence, it follows

[TABLE]

Therefore, we conclude

[TABLE]

I Remark for Theorem 10

For the special case of noise according to the double-exponential distribution with parameter $\beta$, we have

[TABLE]

For every $\theta\in\theta_{n}$ and every $S\subseteq N$ of cardinality $k$ and $i,j,y\in S$, we can easily check that

[TABLE]

Furthermore, the following two relations hold

[TABLE]

Since

[TABLE]

and

[TABLE]

we have that

[TABLE]

and

[TABLE]

J Proof of Theorem 19

Let $p^{e}$ denote the probability that the point score ranking method incorrectly classifies at least one item:

[TABLE]

Let $R_{i}$ denote the point score of item $i\in N$. If the point scores are such that $R_{v}>m/n$ for every $v\in N_{1}$ and $R_{v}<m/n$ for every $v\in N_{2}$, then this implies a correct classification. Hence, it must be that in the event of a misclassification of an item, $R_{v}\leq m/n$ for some $v\in N_{1}$ or $R_{v}\geq m/n$ for some $v\in N_{2}$. Combining this with the union bound, we have

[TABLE]

Let $i$ and $j$ be arbitrarily fixed items such that $i\in N_{1}$ and $j\in N_{2}$. We will show that for $t\in\{1,2,\ldots,m\}$,

[TABLE]

and

[TABLE]

By the Chernoff bound (29) we have

[TABLE]

Similarly, by the Chernoff bound (30), we have

[TABLE]

Combining with (90), we have $p^{e}\leq\delta$.

In the remainder of the proof we show that (91) and (92) hold.

Let ${\mathcal{A}}$ contain all $A\subseteq N$ such that $|A|=k-1$ and $A\cap\{i,j\}=\emptyset$ and ${\mathcal{B}}$ contain all $B\subseteq N$ such that $|B|=k-2$ and $B\cap\{i,j\}=\emptyset$. Then, we have

[TABLE]

where

[TABLE]

and

[TABLE]

Let $\mathbf{b}$ be a $k-1$-dimensional vector with all elements equal to $b$. Then, note that

[TABLE]

By limited Taylor series development, we have

[TABLE]

where

[TABLE]

Hence, it follows that for every $\theta_{A}\in\{-b,b\}^{k-1}$,

[TABLE]

Under the condition of the theorem, we have

[TABLE]

Hence, combining with (97), for every $\theta_{A}\in\{-b,b\}^{k-1}$,

[TABLE]

By the same arguments, we can show that

[TABLE]

where $\theta^{(b)}_{B}\in\{-b,b\}^{k-1}$ and $\theta^{(-b)}_{B}\in\{-b,b\}^{k-1}$ are $(k-1)$-dimensional vectors with the first two elements equal to $b$ and $-b$, respectively, and other elements equal to the parameters of items $B$.

Since comparison sets are sampled uniformly at random without replacement,

[TABLE]

and

[TABLE]

From (93), (98), (99), (100) and (101), we have

[TABLE]

Using this inequality together with the following facts (i) $\mathbf{P}[y_{t}=v]=\mathbf{P}[y_{t}=i]$ for every $v\in N_{1}$, (ii) $\mathbf{P}[y_{t}=v]=\mathbf{P}[y_{t}=j]$ for every $v\in N_{2}$, (iii) $\sum_{v\in N}\mathbf{P}[y_{t}=v]=1$, and (iv) $|N_{1}|=|N_{2}|=n/2$, it can be readily shown that

[TABLE]

which establishes (91). By the same arguments one can establish (92).

K Proof of Theorem 20

Suppose that $n$ is a positive even integer and $\theta$ is the parameter vector such that $\theta_{i}=b$ for $i\in N_{1}$ and $\theta_{i}=-b$ for $i\in N_{2}$, where $N_{1}=\{1,2,\ldots,n/2\}$ and $N_{2}=\{n/2+1,\ldots,n\}$. Let $\theta^{\prime}$ be the parameter vector that is identical to $\theta$ except for swapping the first and the last item, i.e. $\theta^{\prime}_{i}=b$ for $i\in N_{1}^{\prime}$ and $\theta^{\prime}_{i}=-b$ for $i\in N_{2}^{\prime}$, where $N_{1}^{\prime}=\{n,2,\ldots,n/2\}$ and $N_{2}^{\prime}=\{n/2+1,\ldots,n-1,1\}$.

We denote with $\mathbf{P}_{\theta}[A]$ and $\mathbf{P}_{\theta^{\prime}}[A]$ the probabilities of an event $A$ under hypothesis that the generalized Thurstone model is according to parameter $\theta$ and $\theta^{\prime}$, respectively. We denote with $\mathbf{E}_{\theta}$ and $\mathbf{E}_{\theta^{\prime}}$ the expectations under the two respective distributions.

Given observed data $(\mathbf{S},\mathbf{y})=(S_{1},y_{1}),\ldots,(S_{m},y_{m})$, we denote the log-likelihood ratio statistic $L(\mathbf{S},\mathbf{y})$ as follows

[TABLE]

where $\rho_{t}(S)$ is the probability that $S$ is drawn at time $t$.

The proof follows the following two steps:

Step 1:

We show that for given $\delta\in[0,1]$, for the existence of an algorithm that correctly classifies all the items with probability at least $1-\delta$, it is necessary that the following condition holds

[TABLE]

Step 2:

We show that

[TABLE]

where $\sigma^{2}_{\theta^{\prime}}[L(\mathbf{S},\mathbf{y})]$ denotes the variance of random variable $L(\mathbf{S},\mathbf{y})$ under a generalized Thurstone model with parameter $\theta^{\prime}$.

By Chebyshev’s inequality, for every $g\in\mathbf{R}$,

[TABLE]

Using this for $g=\log(n/\delta)-\mathbf{E}_{\theta^{\prime}}[L(\mathbf{S},\mathbf{y})]$, it follows that (103) implies the following condition:

[TABLE]

Further combining with (104) and (105), we obtain

[TABLE]

which is the condition asserted in the theorem.

Proof of Step 1.

Let us define the following two events

[TABLE]

and

[TABLE]

Let $B^{c}$ denote the complement of event $B$.

Note that

[TABLE]

where the second equation holds because $B\subseteq A$ and every possible partition in $A$ has the same probability under $\theta$.

For every $g\in\mathbf{R}$, we have

[TABLE]

Now, note

[TABLE]

where in the second equation we use the standard change of measure argument.

Since the algorithm correctly classifies all the items with probability at least $1-\delta$, we have

[TABLE]

For $g=\log(n/\delta)$, from (K) and (106), it follows that

[TABLE]

where the last inequality is by the conditions of the theorem.

Proof of Step 2.

If the observed comparison sets $S_{1},S_{2},\ldots,S_{m}$ are such that $S_{t}\cap\{1,n\}=\emptyset$, for every observation $t$, then we obviously have

[TABLE]

We therefore consider the case when $S_{t}\cap\{1,n\}\neq\emptyset$.

Using (94), (95), and (96), we have for every $S$ and $i\in S$,

[TABLE]

where the last inequality is obtained from the condition of this theorem.

From (107), for every comparison set $S$ such that $S\cap\{1,n\}\neq\emptyset$, we have

[TABLE]

which is because for every comparison set $S$ such that $1\in S$,

[TABLE]

and, for every comparison set $S$ such that $n\in S$,

[TABLE]

From (107) and the assumption of the theorem, we have

[TABLE]

For simplicity of notation, let

[TABLE]

Then, for all $S$ such that $S\cap\{1,n\}\neq\emptyset$, we have

[TABLE]

which is obtained from

• (i)

$p_{i,S}(\theta)\geq 1/(2k)$ for all $i\in S$ that holds by (110),

• (ii)

$\sum_{i\in S}(p_{i,S}(\theta^{\prime})-p_{i,S}(\theta))^{2}=2D^{2}$ from (109),

• (iii)

$a\log\frac{a}{b}\leq\frac{(a-b)^{2}}{2b}+a-b$.

Similarly to (112), from (i) and (ii) and $a\left(\log\frac{a}{b}\right)^{2}\leq\frac{(a-b)^{2}}{a\wedge b}\left(1+\frac{|a-b|}{3(a\wedge b)}\right)$, we have

[TABLE]

Since

[TABLE]

and according to the model, the input observations are independent, from (112) and (113), we have

[TABLE]

and

[TABLE]

L Characterizations of $\partial p_{k}(\mathbf{0})/\partial x_{1}$

In this section, we note several different representations of the parameter $\partial p_{k}(\mathbf{0})/\partial x_{1}$.

First, note that

[TABLE]

The integral corresponds to $\mathbf{E}[f(X)]$ where $X$ is a random variable whose distribution is equal to that of a maximum of $k-1$ independent and identically distributed random variables with cumulative distribution $F$.

Second, suppose that $F$ is a cumulative distribution function with its support contained in $[-a,a]$, and that has a differentiable density function $f$. Then, we have

[TABLE]

where

[TABLE]

and

[TABLE]

The identity (122) is shown to hold as follows. Note that

[TABLE]

By integrating over $[-a,a]$, we obtain

[TABLE]

Combining with the fact

[TABLE]

we obtain (122).

Note that $B_{F,k}=\mathbf{E}[-f^{\prime}(X)]/(k(k-1))$ where $X$ is a random variable with distribution that corresponds to that of a maximum of $k$ independent samples from the cumulative distribution function $F$. Note also that if, in addition, $f$ is an even function, then (i) $B_{F,k}\geq 0$ and (ii) $B_{F,k}$ is increasing in $k$.

Third, for any cumulative distribution function $F$ with an even density function $f$, we have $F(-x)=1-F(x)$ for all $x\in\mathbf{R}$. In this case, we have the identity

[TABLE]

M Proof of Lemma 21

The upper bound follows by noting that that $B_{F,k}$ in (122) is such that $B_{F,k}=\Omega(1/k^{2})$. Hence, it follows that

[TABLE]

The lower bound follows by noting that for every cumulative distribution function $F$ such that there exists a constant $C>0$ such that $f(x)\leq C$ for all $x\in\mathbf{R}$,

[TABLE]

Hence, $\gamma_{F,k}\geq(1/C)(k-1)/k^{3}=\Omega(1/k^{2})$.

N Derivations of parameter $\gamma_{F,k}$

We derive explicit expressions for parameter $\gamma_{F,k}$ for our example generalized Thurstone choice models introduced in Section 2

Recall from (19) that we have that

[TABLE]

where

[TABLE]

Gaussian distribution

A cumulative distribution function $F$ is said to have a type-3 domain of maximum attraction if the maximum of $r$ independent and identically distributed random variables with cumulative distribution function $F$ has as a limit a double-exponential cumulative distribution function:

[TABLE]

where

[TABLE]

and

[TABLE]

It is a well known fact that any Gaussian cumulative distribution function has a type-$3$ domain of maximum attraction. Let $\Phi$ denote the cumulative distribution function of a standard normal random variable, and let $\phi$ denotes its density.

Note that

[TABLE]

Now, note that

[TABLE]

It is readily checked that $e^{-a_{r}^{2}/2}\sim 1/r$ and $\Gamma(a_{r}b_{r}+1)=O(r^{\epsilon})$ for every constant $\epsilon>0$. Hence, we have that

[TABLE]

and thus, $\partial p_{k}(\mathbf{0})/\partial x_{1}=O(1/k^{2-\epsilon})$. Hence,

[TABLE]

Double-exponential distribution

Note that $f(x)=\frac{1}{\beta}e^{-\frac{x+\beta\gamma}{\beta}}F(x)$. Hence, we have

[TABLE]

Laplace distribution

Let $\beta=\sigma/\sqrt{2}$. Note that

[TABLE]

[TABLE]

and

[TABLE]

Combining with (123), we obtain

[TABLE]

Uniform distribution

Note that

[TABLE]