Accelerated MM Algorithms for Ranking Scores Inference from Comparison Data

TL;DR

This paper analyzes the convergence of MM algorithms for ranking inference models, establishes their linear convergence rates, and introduces an accelerated version that significantly improves efficiency based on theoretical insights and experimental validation.

Contribution

The paper provides a tight characterization of MM algorithm convergence rates for ranking models and proposes an accelerated algorithm with proven efficiency gains.

Findings

MM algorithms converge linearly with rates linked to data connectivity and prior parameters.

The accelerated MM algorithm outperforms classical methods in synthetic and real data.

Theoretical analysis guides the parameter rescaling for improved convergence speed.

Abstract

In this paper, we study a popular method for inference of the Bradley-Terry model parameters, namely the MM algorithm, for maximum likelihood estimation and maximum a posteriori probability estimation. This class of models includes the Bradley-Terry model of paired comparisons, the Rao-Kupper model of paired comparisons allowing for tie outcomes, the Luce choice model, and the Plackett-Luce ranking model. We establish tight characterizations of the convergence rate for the MM algorithm, and show that it is essentially equivalent to that of a gradient descent algorithm. For the maximum likelihood estimation, the convergence is shown to be linear with the rate crucially determined by the algebraic connectivity of the matrix of item pair co-occurrences in observed comparison data. For the Bayesian inference, the convergence rate is also shown to be linear, with the rate determined by a parameter of the prior distribution in a way that can make the convergence arbitrarily slow for small values of this parameter. We propose a simple modification of the classical MM algorithm that avoids the observed slow convergence issue and accelerates the convergence. The key component of the accelerated MM algorithm is a parameter rescaling performed at each iteration step that is carefully chosen based on theoretical analysis and characterisation of the convergence rate. Our experimental results, performed on both synthetic and real-world data, demonstrate the identified slow convergence issue of the classic MM algorithm, and show that significant efficiency gains can be obtained by our new proposed method.

Tables

Table 1. Table 1: Dataset properties.

Dataset	$m$	$n$	$d\left(𝐌\right)$	$a\left(𝐌\right)$
GIFGIF: A (full)	161,584	6,123	83	0
GIFGIF: C (full)	108,126	6,122	56	0
GIFGIF: H (full)	225,695	6,124	153	0
GIFGIF: A (sample)	702	252	15	0.671
GIFGIF: C (sample)	734	256	28	0.569
GIFGIF: H (sample)	1040	251	23	1.357
Chess (full)	65,030	8,631	155	0
Chess (sample)	13,181	985	135	1.773
NASCAR	64,596	83	1,507	39.338

Table 2. Table 2: Number of iterations for the MM algorithm and accelerated MM algorithm (AccMM).

Dataset	Algorithm	$\beta =0$	$0.01$	$0.1$	$1$	$10$
GIFGIF: A (full)	MM	MLE	572	125	70	16
	AccMM	non-existant	509	123	42	13
GIFGIF: C (full)	MM	MLE	733	150	49	13
	AccMM	non-existant	551	93	37	13
GIFGIF: H (full)	MM	MLE	1,127	149	98	21
	AccMM	non-existant	1,044	159	51	18
GIFGIF: A (sample)	MM	145	854	177	26	7
	AccMM		125	81	22	7
GIFGIF: C (sample)	MM	130	694	151	39	9
	AccMM		111	78	36	9
GIFGIF: H (sample)	MM	216	1,234	237	38	8
	AccMM		146	72	26	8
Chess (full)	MM	MLE	2,217	581	113	33
	AccMM	non-existant	2,291	302	49	25
Chess (sample)	MM	121	122	91	74	19
	AccMM		117	93	48	16
NASCAR	MM	11	695	971	58	10
	AccMM		11	11	10	6
NASCAR ($\xi ={10}^{-5}$)	MM	14	1,528	2,069	105	16
	AccMM		14	14	12	7
NASCAR ($\xi ={10}^{-6}$)	MM	17	2,362	3,223	157	23
	AccMM		17	16	14	8
NASCAR ($\xi ={10}^{-8}$)	MM	22	4,029	5,544	261	36
	AccMM		22	21	18	11

Table 3. Table 3: Number of iterations for the MM algorithm with normalization such that ∑ i = 1 n e w i = 1 superscript subscript 𝑖 1 𝑛 superscript 𝑒 subscript 𝑤 𝑖 1 \sum_{i=1}^{n}e^{w_{i}}=1 in each iteration.

Dataset	Algorithm	$\beta =0$	$0.01$	$0.1$	$1$	$10$
GIFGIF: A (sample)	MM norm	141	136	125	154	No convergence
GIFGIF: C (sample)	MM norm	128	130	144	350	No convergence
GIFGIF: H (sample)	MM norm	205	203	184	124	No convergence

Table 4. Table 4: Parameters d ​ ( 𝐌 ) 𝑑 𝐌 d({\mathbf{M}}) and a ​ ( 𝐌 ) 𝑎 𝐌 a({\mathbf{M}}) for matrix 𝐌 𝐌 {\mathbf{M}} being the adjacency matrix of graph G 𝐺 G with n 𝑛 n vertices

$G$	$d\left(𝐌\right)$	$a\left(𝐌\right)$	$\frac{d\left(𝐌\right)}{a\left(𝐌\right)}$
complete	$n-1$	$n$	$\Theta \left(1\right)$
star	$n-1$	$1$	$\Theta \left(n\right)$
circuit	$2$	$2\left(1-\cos \left(\frac{\pi }{n}\right)\right)\sim \frac{4{\pi }^2}{n^2}$	$\Theta \left(n^2\right)$
path	$2$	$2\left(1-\cos \left(\frac{2\pi }{n}\right)\right)\sim \frac{{\pi }^2}{n^2}$	$\Theta \left(n^2\right)$

Equations

$$T^{\mathrm{ML}}=O\left(\frac{d(\mathbf{M})}{a(\mathbf{M})}\log \left(\frac{1}{ϵ}\right)\right)$$

$$T^{\mathrm{MAP}}=O\left(\left(\frac{d(\mathbf{M})}{\beta }+1\right)\log \left(\frac{1}{ϵ}\right)\right)$$

$$T_{\mathrm{Acc}}^{\mathrm{MAP}}=O\left(\min \left\{\frac{d(\mathbf{M})}{a(\mathbf{M})},\frac{d(\mathbf{M})}{\beta }\right\}\log \left(\frac{1}{ϵ}\right)\right).$$

$$\mathrm{Pr}[i\succ j]=\frac{{\theta }_i}{{\theta }_i+{\theta }_j}$$

$$\mathrm{Pr}[i\succ j]=\frac{e^{w_i}}{e^{w_i}+e^{w_j}}.$$

$$\underset{\mathbf{w}\in {\mathrm{I}\mathrm{R}}^n}{\max }\ell (\mathbf{w})$$

$$\ell (\mathbf{w})=\sum _{i=1}^n\sum _{j\ne i}{d}_{i,j}\left(w_i-\log \left(e^{w_i}+e^{w_j}\right)\right)$$

$$\rho (\mathbf{w})=\ell (\mathbf{w})+{\ell }_0(\mathbf{w})$$

$${\ell }_0(\mathbf{w})=\sum _{i=1}^n\left((\alpha -1)w_i-\beta e^{w_i}\right).$$

$${\mathbf{x}}^{(t+1)}=\arg \underset{\mathbf{x}}{\min }g(\mathbf{x};{\mathbf{x}}^{(t)}).$$

$$\underline{\ell }(\mathbf{x};\mathbf{y})=$$

$${\underline{\ell }}_{ij}(\mathbf{x};\mathbf{y})=d_{i,j}\left(x_i-\frac{e^{x_i}+e^{x_j}}{e^{y_i}+e^{y_j}}-\log \left(e^{y_i}+e^{y_j}\right)+1\right).$$

$${\theta }_i^{(t+1)}=\frac{{\sum }_{j=1}^n{d}_{i,j}}{{\sum }_{j=1}^n\frac{m_{i,j}}{{\theta }_i^{(t)}+{\theta }_j^{(t)}}}.$$

$$\underline{\rho }(\mathbf{x};\mathbf{y})=\underline{\ell }(\mathbf{x};\mathbf{y})+{\ell }_0(\mathbf{x})$$

$${\theta }^{(t+1)}=\frac{\alpha -1+{\sum }_{j\ne i}{d}_{i,j}}{\beta +{\sum }_{j\ne i}\frac{m_{i,j}}{{\theta }_i^{(t)}+{\theta }_j^{(t)}}}.$$

$${\mathbf{x}}^{(t+1)}={\mathbf{x}}^{(t)}-\eta \nabla f({\mathbf{x}}^{(t)}).$$

$$f({\mathbf{x}}^{(t+1)})-f({\mathbf{x}}^{\star })\le (1-\alpha )(f({\mathbf{x}}^{(t)})-f({\mathbf{x}}^{\star }))$$

$$f(\mathbf{x})-f(\mathbf{y})\le \nabla f(\mathbf{x}{)}^{\top }(\mathbf{x}-\mathbf{y})-\frac{\gamma }{2}||\mathbf{x}-\mathbf{y}|{|}^2.$$

$$||\nabla f(\mathbf{x})-\nabla f(\mathbf{y})||\le \mu ||\mathbf{x}-\mathbf{y}||.$$

$$|f(\mathbf{x})-f(\mathbf{y})-\nabla f(\mathbf{y}{)}^{\top }(\mathbf{x}-\mathbf{y})|\le \frac{\mu }{2}||\mathbf{x}-\mathbf{y}|{|}^2.$$

$$f(\mathbf{x})-f({\mathbf{x}}^{\ast })\le \frac{1}{2\gamma }||\nabla f(\mathbf{x})|{|}^2$$

$$g(\mathbf{x};\mathbf{y})-f(\mathbf{x})\le \frac{\delta }{2}||\mathbf{x}-\mathbf{y}|{|}^{2}for\; all\mathbf{x},\mathbf{y}\in {\mathcal{X}}_{\mu }.$$

$$\gamma =c_{\omega }{\lambda }_2({\mathbf{L}}_{\mathbf{M}})and\mu =\frac{1}{4}{\lambda }_n({\mathbf{L}}_{\mathbf{M}})$$

$${\Pi }_c(\mathbf{x})=\mathbf{x}+c1.$$

$${\alpha }_{\mathbf{M},\omega }=2c_{\omega }\frac{a(\mathbf{M})}{d(\mathbf{M})}.$$

$$T=O\left(\frac{d(\mathbf{M})}{a(\mathbf{M})}\log \left(\frac{1}{ϵ}\right)\right).$$

$$\delta =\frac{1}{2}{e}^{2\omega }d(\mathbf{M}).$$

$${\alpha }_{\mathbf{M},\omega }=2c_{\omega }^{\prime }\frac{a(\mathbf{M})}{d(\mathbf{M})}$$

$$\gamma =e^{-\omega }\beta and\mu =\frac{1}{4}{\lambda }_n({\mathbf{L}}_{\mathbf{M}})+e^{\omega }\beta .$$

$${\alpha }_{\mathbf{M},\omega }=\frac{2e^{-\omega }\beta }{d(\mathbf{M})+2e^{\omega }\beta }.$$

Full text

\SingleSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough

\RUNAUTHOR

Vojnovic, Yun and Zhou

\RUNTITLE

Accelerated MM Algorithms for Ranking Scores Inference from Comparison Data

\TITLE

Accelerated MM Algorithms for Ranking Scores Inference from Comparison Data

\ARTICLEAUTHORS\AUTHOR

Milan Vojnovic \AFFDepartment of Statistics, London School of Economics, London, United Kingdom, WC2A 2AE, \EMAILm.vojnovic@lse.ac.uk \AUTHORSe-Young Yun \AFFGraduate School of AI, KAIST, Daejeon, South Korea, 34141, \EMAILyunseyoung@kaist.ac.kr \AUTHORKaifang Zhou \AFFDepartment of Statistics, London School of Economics, London, United Kingdom, WC2A 2AE, \EMAILk.zhou3@lse.ac.uk

\ABSTRACT

The problem of assigning ranking scores to items based on observed comparison data, e.g., paired comparisons, choice, and full ranking outcomes, has been of continued interest in a wide range of applications, including information search, aggregation of social opinions, electronic commerce, online gaming platforms, and more recently, evaluation of machine learning algorithms. The key problem is to compute ranking scores, which are of interest for quantifying the strength of skills, relevancies or preferences, and prediction of ranking outcomes when ranking scores are estimates of parameters of a statistical model of ranking outcomes. One of the most popular statistical models of ranking outcomes is the Bradley-Terry model for paired comparisons (equivalent to multinomial logit model), and its extensions to choice and full ranking outcomes. The problem of computing ranking scores under the Bradley-Terry models amounts to estimation of model parameters.

In this paper, we study a popular method for inference of the Bradley-Terry model parameters, namely the MM algorithm, for maximum likelihood estimation and maximum a posteriori probability estimation. This class of models includes the Bradley-Terry model of paired comparisons, the Rao-Kupper model of paired comparisons allowing for tie outcomes, the Luce choice model, and the Plackett-Luce ranking model. We establish tight characterizations of the convergence rate for the MM algorithm, and show that it is essentially equivalent to that of a gradient descent algorithm. For the maximum likelihood estimation, the convergence is shown to be linear with the rate crucially determined by the algebraic connectivity of the matrix of item pair co-occurrences in observed comparison data. For the Bayesian inference, the convergence rate is also shown to be linear, with the rate determined by a parameter of the prior distribution in a way that can make the convergence arbitrarily slow for small values of this parameter. We propose a simple modification of the classical MM algorithm that avoids the observed slow convergence issue and accelerates the convergence. The key component of the accelerated MM algorithm is a parameter rescaling performed at each iteration step that is carefully chosen based on theoretical analysis and characterisation of the convergence rate.

Our experimental results, performed on both synthetic and real-world data, demonstrate the identified slow convergence issue of the classic MM algorithm, and show that significant efficiency gains can be obtained by our new proposed method.

\KEYWORDS

Rank aggregation, Generalized Bradley-Terry model, Maximum likelihood estimation, Bayesian inference, Convex optimization, Gradient descent, MM algorithm, Rate of convergence

1 Introduction

Rank aggregation is an important task that arises in a wide-range of applications, including recommender systems, information retrieval, online gaming, sport competitions, and evaluation of machine learning algorithms. Given a set of items, the goal is to infer ranking scores of items or an ordering of items based on observed data that contains partial orderings of items. A common scenario is that of paired comparisons, e.g., player A defeats player B in a game, product A is preferred over product B by a user, and machine learning algorithm A outperforms machine learning algorithm B in an evaluation. In such scenarios, a common goal is not only to compute an aggregate ranking of items, but also to compute ranking scores, which represent strengths of individual items. Such ranking scores are used for predicting outcomes of future ranking outcomes, such as predicting outcomes of matches in online games and sport contests, and predicting preference of users in product shopping or movie watching scenarios, among others. Note that, importantly, observations are not restricted to paired comparisons, but may also include other types of comparison data, such as choice (e.g., product A chosen from a set of two or more products) or full ranking (e.g., a ranking list of players or teams participating in a competition).

In this paper, our goals are (a) to shed light on the efficiency of one of the most popular iterative optimization methods for inference of ranking scores, namely the MM algorithm, where ranking scores correspond to parameter estimates of popular Bradley-Terry family of models, and (b) to propose an accelerated MM algorithm that resolves a slow convergence issue found to hold for a classic MM algorithm.

Related work

Statistical models of ranking data play an important role in a wide range of applications, including learning to rank in information retrieval (Burges et al. (2006), Li (2011)), skill rating in sport games (Elo (1978)), online gaming platforms (Herbrich et al. (2006)), and evaluation of machine learning algorithms by comparing them with each other (Balduzzi et al. (2018)).

A common class of statistical models of ranking data are generalized Bradley-Terry models, which accommodate paired comparisons with win-lose outcomes (Zermelo (1929), Bradley and Terry (1952), Bradley (1954)), paired comparisons with win-lose-draw outcomes (Rao and Kupper (1967)), choices from comparison sets of two or more items, e.g., Luce choice model (Luce (1959)), full ranking outcomes for comparison sets of two or more items, e.g., Plackett-Luce ranking model (Plackett (1975)), as well as group comparisons (Huang et al. (2006, 2008)). These models can be derived from suitably defined latent variable models, where items are associated with independent latent performance random variables, which is in the spirit of the well-known Thurstone model of comparative judgment (Thurstone (1927)).

An iterative optimization algorithm for the maximum likelihood (ML) parameter estimation (MLE) of the Bradley-Terry model has been known since the original work of Zermelo (1929). Lange et al. (2000) showed that this algorithm belongs to the class of MM optimization algorithms. Here MM refers to either minorize-maximization or majorize-minimization, depending on whether the optimization problem is maximization or minimization of an objective function. Lange (2016) provided a book on MM algorithms and Hunter and Lange (2004) provided a tutorial. Mairal (2015) established some convergence results for incremental MM algorithms.

In a seminal paper, Hunter (2004) derived MM algorithms for generalized Bradley-Terry models as well as sufficient conditions for their convergence to ML estimators using the framework of MM optimization algorithms. For the Bradley-Terry model of paired comparisons, a necessary and sufficient condition for the existence of a ML estimator is that the directed graph whose vertices correspond to items and edges represent outcomes of paired comparisons is connected. In other words, the set of items cannot be partitioned in two sets such that none of the items in one partition won against an item in other partition.

A Bayesian inference method for generalized Bradley-Terry models was proposed by Caron and Doucet (2012), showing that classical MM algorithms can be reinterpreted as special instances of Expectation-Maximization (EM) algorithms associated with suitably defined latent variables and proposed some original extensions. This amounts to MM algorithms for maximum a posteriori probability (MAP) parameter estimation, for a specific family of prior distributions. This prior distribution is a product-form distribution with Gamma$(\alpha,\beta)$ marginal distributions, where $\alpha\geq 1$ is the shape parameter and $\beta>0$ is the rate parameter. Importantly, unlike to the ML estimation, the MAP estimate is always guaranteed to exist, for any given observation data.

Algorithms for fitting Bradley-Terry model parameters are implemented in open source software packages, including BradleyTerry2 (Turner and Firth (2012)), BradleyTerryScalable (Kaye and Firth (2020)), and Choix (Maystre (2018)). The first package uses a Fisher scoring algorithm (a second-order optimization method), while the latter two use MM algorithms (a first-order optimization method). First-order methods are generally preferred over second-order methods for fitting high-dimensional models using large training datasets.

While the conditions for convergence of MM algorithms for generalized Bradley-Terry models are well understood, to the best of our knowledge, not much is known about their convergence rates for either ML or MAP estimation. In this paper, we close this gap by providing tight characterizations of convergence rates. Our results identify key properties of input data that determine the convergence rate, and in the case of MAP estimation, how the convergence rate depends on prior distribution parameters. Our results show that MM algorithms, commonly used for MAP estimation for generalized Bradley-Terry models, can have a slow convergence for some prior distributions.

Recent research on statistical models of paired comparisons focused on characterization of the accuracy of parameter estimators and development of new, scalable parameter estimation methods, e.g., Guiver and Snelson (2009), Wauthier et al. (2013), Hajek et al. (2014), Rajkumar and Agarwal (2014), Chen and Suh (2015), Shah et al. (2016), Vojnovic and Yun (2016), Khetan and Oh (2016a), Negahban et al. (2017), Borkar et al. (2016), Chen et al. (2019). Note that the question about statistical estimation accuracy and computation complexity tradeoff is out of the scope of our paper, and this was studied in the above cited papers. The focus of our work is on convergence properties of first-order iterative optimization methods for parameter estimation of Bradley-Terry models. Here ”first-order” refers to optimization methods that are restricted to value oracle access to gradients of the optimization objective function, thus not allowing access to second-order properties such as values of the Hessian matrix. Specifically, we are interested in convergence properties of first-order methods for ML and MAP estimation objectives. It is noteworthy that some recently proposed algorithms show empirically faster convergence rate than MM, e.g., Negahban et al. (2017), Maystre and Grossglauser (2015), Agarwal et al. (2018), but it is hard to apply them for the MAP estimation objective. We thus restrict our attention to MM and gradient descent algorithms which are able to solve both MLE and MAP optimization problems.

A preliminary version of our paper was published in Vojnovic et al. (2020), which contains results on convergence rates of gradient descent and MM algorithms (Section 3). In the present paper, we extend our prior work by proposing new accelerated algorithms and establishing their theoretical guarantees (Section 4) as well as demonstrating their efficiency through numerical evaluations (Section 5).

Summary of our contributions

We present tight characterizations of the rate of convergence of gradient descent and MM algorithms for ML and MAP estimation for generalized Bradley-Terry models. Our results show that both gradient descent and MM algorithms have linear convergence with convergence rates differing only in constant factors. An iterative optimization algorithm that has linear convergence is generally considered to be fast in the space of first-order optimization algorithms, and many first-order algorithms cannot guarantee a linear convergence. For example, standard stochastic gradient descent algorithm is known to have sub-linear convergence, see, e.g., Bubeck (2015). We provide explicit bounds on convergence rates that provide insights into which properties of observed comparison data play a key role for the rate of convergence.

Specifically, we show that the rate of convergence critically depends on certain properties of the matrix of item pair co-occurrences, ${\mathbf{M}}$, of input comparison data. We found that two key properties are: (a) maximum number of paired comparisons per item (denoted as $d({\mathbf{M}})$) and (b) the algebraic connectivity of matrix ${\mathbf{M}}$ (denoted as $a({\mathbf{M}})$). Intuitively, $a({\mathbf{M}})$ quantifies how well is the graph of paired comparisons connected. Here $a({\mathbf{M}})$ is the Fiedler value (eigenvalue), cf. Fiedler (1973), defined as the second smallest eigenvalue of the Laplacian matrix $\mathbf{L}_{\mathbf{M}}={\mathbf{D}}_{{\mathbf{M}}}-{\mathbf{M}}$, where ${\mathbf{D}}_{{\mathbf{M}}}$ is the diagonal matrix with diagonal elements equal to the row sums of ${\mathbf{M}}$. The Fiedler value of a matrix of paired comparison counts is known to play a key role in determining the MLE accuracy, e.g., Hajek et al. (2014), Shah et al. (2016), Khetan and Oh (2016b), Vojnovic and Yun (2016), Negahban et al. (2017). These works characterized the number of samples needed to estimate the true parameter vector within a statistical estimation error tolerance. This is different from the problem of characterizing the number of iterations needed for an iterative optimization algorithm to compute a ML or a MAP parameter estimate satisfying an error tolerance condition, which is studied in this paper.

Our results reveal the following facts about convergence time, defined as the number of iterations that an iterative optimization algorithm takes to reach the value of the underlying objective function within a given error tolerance parameter $\epsilon>0$ of the optimum value.

For the ML objective, we show that the convergence time satisfies

[TABLE]

which reveals that the rate of convergence critically depends on the connectivity of the graph of paired comparisons in observed data.

For the MAP estimation, we show that the convergence time satisfies

[TABLE]

where, recall, $\beta>0$ is the rate parameter of the Gamma prior distribution. This bound is shown to be tight for some input data instances. We observe that the convergence time for the MAP estimation problem can be arbitrarily large for small enough parameter $\beta$, where small values of parameter $\beta$ correspond to less informative prior distributions.

Our results identify a slow rate of convergence issue for gradient descent and MM algorithms for the MAP estimation problem. While the MAP estimation alleviates the issue of the non-existence of a ML estimator when the graph of paired comparisons is disconnected, it can have a much slower convergence than ML when the graph of paired comparisons is connected. Perhaps surprisingly, the rate of convergence has a discontinuity at $\beta=0$, in the sense that for $\alpha=1$ and $\beta=0$, the MM algorithm for the MAP estimation corresponds to the classic MM algorithm for ML estimation, and in this case, the convergence bound (1) holds, while for the MAP estimation, the convergence time grows arbitrarily large as $\beta$ approaches [math] from above.

We propose an acceleration method for the MAP estimation objective, with convergence time bounded as follows

[TABLE]

This acceleration method resolves the slow convergence issue of classic MM algorithm for the MAP estimation for generalized Bradley-Terry models. This accelerated method does not have a discontinuity at $\beta=0$ with respect to the rate of convergence: as $\beta$ approaches [math] from above, the convergence time bound corresponds to that of the MM algorithm for ML estimation. The acceleration method normalizes the parameter vector estimate in each iteration of the gradient descent or MM algorithm using a particular normalization that ensures (a) that the value of the objective function is non decreasing along the sequence of estimated parameter vectors and (b) that the objective function satisfies certain smoothness and strong convexity properties that ensure high convergence rate. This amounts to a slight modification of the classical MM algorithm that resolves the identified convergence issue. This acceleration method is derived by using a theoretical framework that may be of general interest. This framework can be applied to different statistical models of ranking data and prior distributions for Bayesian inference of parameters of these models.

We present numerical evaluation of the convergence time of different iterative optimization algorithms using input data comparisons from a collection of real-world datasets. These results demonstrate the extent of the slow convergence issue of the existing MM algorithm for MAP estimation and show significant speed ups achieved by our accelerated MM algorithm.

Our theoretical results are established by using the framework of convex optimization analysis and spectral theory of matrices. In particular, the convergence rate bounds are obtained by using concepts of smooth and strongly convex functions. We derive accelerated iterative optimization algorithms based on a general approach that may be of independent interest. This approach transforms the parameter estimator in each iteration so that certain conditions are preserved for the gradient vector and the Hessian matrix of the objective function. For generalized Bradley-Terry models, this transformation turns out to be simple, yielding a practical algorithm.

Organization of the paper

In Section 2, we present problem formulation and some background material. Section 3 contains our main results on characterization of convergence rates of gradient descent and MM algorithms; for simplicity of exposition, we focus only on the Bradley-Terry model of paired comparisons. Section 4 presents our accelerated algorithms for MAP estimation. Section 5 contains our numerical results. We conclude in Section 6. Section 7 contains all our proofs, additional discussions, and extensions to generalized Bradley-Terry models.

2 Problem formulation

According to the Bradley-Terry model of paired comparisons with win-lose outcomes, each comparison of items $i$ and $j$ has an independent outcome: either $i$ wins against $j$ ($i\succ j$) or $j$ wins against $i$ ($j\succ i$). The distribution of outcomes is given by

[TABLE]

where $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{n})^{\top}\in\mathrm{I\!R}_{+}^{n}$ is the parameter vector. The Bradley-Terry model of paired comparisons was studied by many, e.g., Ford (1957), Dykstra, Jr. (1956, 1960), Simons and Yao (1999) and is covered by classic books on categorical data analysis, e.g., Agresti (2002).

We will sometimes use the parametrization $\theta_{i}=e^{w_{i}}$ when it is simpler to express an equation or when we want to make a connection with the literature using this parametrization. Using parameterization ${\mathbf{w}}=(w_{1},w_{2},\ldots,w_{n})^{\top}\in\mathrm{I\!R}^{n}$, we have

[TABLE]

All our convergence results are for the model with parameter ${\mathbf{w}}$. The Bradley-Terry type models for paired comparisons with ties, choice, and full ranking outcomes, we refer to as generalized Bradley-Terry models, are defined in Section 7.15. Our results apply to all these different models. In the main body of the paper, we focus only on the Bradley-Terry model for paired comparisons in order to keep the presentation simple.

Maximum likelihood estimation

The maximum likelihood parameter estimation problem corresponds to finding ${\mathbf{w}}^{\star}$ that solves the following optimisation problem:

[TABLE]

where $\ell({\mathbf{w}})$ is the log-likelihood function,

[TABLE]

with $d_{i,j}$ denoting the number of observed paired comparisons such that $i\succ j$.

The maximum likelihood optimisation problem (2.2) is a convex optimization problem. Note, however, that the objective function is not a strictly concave function as adding a common constant to each element of the parameter vector keeps the value of the objective function unchanged.

MAP estimation problem

An alternative objective is obtained by using a Bayesian inference framework, which amounts to finding a maximum a posteriori estimate of the parameter vector under a given prior distribution. We consider the Bayesian method introduced by Caron and Doucet (2012), which assumes the prior distribution to be of product-form with marginal distributions such that $\theta_{i}(=e^{w_{i}})$ has a Gamma($\alpha,\beta$) distribution where $\alpha\geq 1$ is the shape parameter and $\beta>0$ is the rate parameter. Note that $\alpha$ and $\beta$ affect the scale of the parameter vector as with respect to the Gamma($\alpha,\beta$) prior distribution, $\theta_{i}$ has the expected value and the mode equal to $\alpha/\beta$ and $(\alpha-1)/\beta$, respectively. For any fixed $\alpha\geq 1$, the density of Gamma($\alpha,\beta$) distribution becomes more flat as $\beta$ approaches zero which corresponds to a less informative prior. According to the assumed prior distribution, $\sum_{i=1}^{n}\theta_{i}\sim\mathrm{Gamma}(n\alpha,\beta)$ and, hence, the mode of $\sum_{i=1}^{n}\theta_{i}$ is $(n\alpha-1)/\beta$. We can interpret the mode of $\sum_{i=1}^{n}\theta_{i}$ as the scale of the parameter vector.

The log-a posteriori probability function can be written as

[TABLE]

where $\ell$ is the log-likelihood function in (2.3) and $\ell_{0}$ is the log-likelihood of the prior distribution given by

[TABLE]

Note that for $\alpha=1$ and $\beta=0$, the log-a posteriori probability function corresponds to the log-likelihood function. For these values of parameters $\alpha$ and $\beta$, MAP and ML estimation problems are equivalent.

MM algorithms

The MM algorithm for minimizing a function $f$ is defined by minimizing a surrogate function that majorizes $f$.

A surrogate function $g({\mathbf{x}};{\mathbf{y}})$ is said to be a majorant function of $f$ if $f({\mathbf{x}})\leq g({\mathbf{x}};{\mathbf{y}})$ and $f({\mathbf{x}})=g({\mathbf{x}};{\mathbf{x}})$ for all ${\mathbf{x}}$ and ${\mathbf{y}}$. The MM algorithm is defined by the iterative updates:

[TABLE]

For maximizing a function $f$, we can analogously define the MM algorithm as minimization of a surrogate function $g$ that minorizes function $f$. Majorization surrogate functions are used for minimization of convex functions, and minorization surrogate functions are used for maximization of concave functions.

The classic MM algorithm for the Bradley-Terry model of paired comparisons uses the following minorization function of $\ell({\mathbf{x}})$:

[TABLE]

where

[TABLE]

It is easy to observe that $\underline{\ell}({\mathbf{x}};{\mathbf{y}})$ is a minorization surrogate function of $\ell({\mathbf{x}})$ by noting that $\log(x)\leq x-1$ and that equality holds if, and only if, $x=1$, which is used to break $\log(e^{x_{i}}+e^{x_{j}})$ terms in the log-likelihood function.

The classic MM algorithm for the ML parameter estimation of the Bradley-Terry model of paired comparisons (Ford (1957), Hunter (2004)), is defined by the following iterative updates, for $i=1,2,\ldots,n$,

[TABLE]

Following Caron and Doucet (2012), the MM algorithm for the MAP parameter estimation of the Bradley-Terry model of paired comparisons is derived for the minorant surrogate function $\underline{\rho}$ of function $\rho$ in (2.4), defined as

[TABLE]

where $\underline{\ell}({\mathbf{x}};{\mathbf{y}})$ is the minorant surrogate function of the log-likelihood function (2.7) and $\ell_{0}$ is the prior log-likelihood function (2.5).

The iterative updates of the MM algorithm are defined by, for $i=1,2,\ldots,n$,

[TABLE]

Note that this iterative optimization algorithm corresponds to the classic MM algorithm for ML estimation (2.8) when $\alpha=1$ and $\beta=0$.

We also consider gradient descent algorithm with constant step size $\eta>0$, which has iterative updates as

[TABLE]

Our goal in this paper is to characterize the rate of convergence of MM algorithms for generalized Bradley-Terry models. We also derive convergence rates for gradient descent algorithms. It is natural to consider gradient descent algorithms as they belong to the class of first-order optimization methods (not requiring second-order quantities such as Hessian of the objective function or its approximations). Intuitively, the rate of convergence of an iterative algorithm quantifies how fast the value of the objective function converges to the optimum value with the number of iterations.

For an iterative optimization method for minimizing function $f$, which outputs a sequence of points ${\mathbf{x}}^{(0)},{\mathbf{x}}^{(1)},\ldots$, we say that there is an $\alpha$-improvement with respect to $f$ at time step $t$ if

[TABLE]

where ${\mathbf{x}}^{\star}$ is a minimizer of $f$. An iterative optimization method is said to have linear convergence if there exist positive constants $\alpha$ and $t_{0}$ such that the method yields an $\alpha$-improvement at each time step $t\geq t_{0}$.

Background on convex analysis

We define some basic concepts from convex analysis that we will use throughout the paper.

Function $f:\mathrm{I\!R}^{n}\rightarrow\mathrm{I\!R}$ is $\gamma$-strongly convex on $\mathcal{X}$ if it satisfies the following subgradient inequality, for all ${\mathbf{x}},{\mathbf{y}}\in\mathcal{X}$:

[TABLE]

$f$ is $\gamma$-strongly convex on $\mathcal{X}$ if, and only if, $f({\mathbf{x}})-\frac{\gamma}{2}||{\mathbf{x}}||^{2}$ is convex on $\mathcal{X}$.

Function $f$ is $\mu$-smooth on $\mathcal{X}$ if its gradient $\nabla f$ is $\mu$-Lipschitz on $\mathcal{X}$, i.e., for all ${\mathbf{x}},{\mathbf{y}}\in\mathcal{X}$,

[TABLE]

For any $\mu$-smooth function $f$ on $\mathcal{X}$, we also have that, for all ${\mathbf{x}},{\mathbf{y}}\in\mathcal{X}$,

[TABLE]

A proof of the last claim can be found in Lemma 3.4 in Bubeck (2015).

Function $f$ satisfies the Polyak-Lojasiewicz inequality on $\mathcal{X}$ (Polyak (1963)) if there exists $\gamma>0$ such that for all ${\mathbf{x}}\in\mathcal{X}$,

[TABLE]

where ${\mathbf{x}}^{*}$ is a minimizer of $f$. When the PL inequality holds on $\mathcal{X}$ for a specific value of $\gamma$, we say that $\gamma$-PL inequality holds on $\mathcal{X}$. If $f$ is $\gamma$-strongly convex on $\mathcal{X}$, then $f$ satisfies the $\gamma$-PL inequality on $\mathcal{X}$.

3 Convergence rates

In this section, we present results on the rate of convergence for gradient descent and MM algorithms for ML and MAP estimation for the Bradley-Terry model of paired comparisons. We first show some general convergence theorems that hold for any strongly convex and smooth function $f$, which characterize the rate of convergence in terms of the strong-convexity and smoothness parameters of $f$, and a parameter of the surrogate function used to define the MM algorithm. These results are then used to derive convergence rate bounds for the Bradley-Terry model.

The results of this section can be extended to other instances of generalized Bradley-Terry models, including the Rao-Kupper model of paired comparisons with tie outcomes, the Luce choice model, and the Plackett-Luce ranking model. These extensions are established by following the same main steps as for the Bradley-Terry model of paired comparisons. The differences lie in the characterization of the strong-convexity and smoothness parameters. The results provide characterizations of the convergence rates that are equivalent to those for the Bradley-Terry model of paired comparisons up to constant factors. We provide details in Section 7.15.

3.1 General convergence theorems

We first present a known result that for any $\mu$-smooth function satisfying the $\gamma$-PL inequality, gradient descent algorithm (2.10) with suitable choice of the step size $\eta$ has a linear convergence with the rate of convergence $1-\gamma/\mu<1$. This result is due to Nesterov (2012) and a simple proof can be found in Boyd and Vandenberghe (2004), Chapter 9.3. This result follows from the following theorem.

Theorem 3.1** **(gradient descent)

Assume $f$ is a convex $\mu$-smooth function on $\mathcal{X}_{\mu}$ satisfying the $\gamma$-PL inequality on ${\mathcal{X}}_{\gamma}\subseteq{\mathcal{X}}_{\mu}$, ${\mathbf{x}}^{*}\in{\mathcal{X}}_{\gamma}$ is a minimizer of $f$, and ${\mathbf{x}}^{(t)}\mapsto{\mathbf{x}}^{(t+1)}$ is according to the gradient descent algorithm (2.10) with step size $\eta=1/\mu$.

Then, if ${\mathbf{x}}^{(t)}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}^{(t+1)}\in{\mathcal{X}}_{\mu}$, there is an $\gamma/\mu$-improvement with respect to $f$ at time step $t$.

Note that if there exists $t_{0}\geq 0$ such that ${\mathbf{x}}^{(t)}\in{\mathcal{X}}_{\gamma}$ for all $t\geq t_{0}$, then Theorem 3.1 implies a linear convergence rate with rate $\gamma/\mu$. Such a $t_{0}$ indeed exists as it can be shown that $||{\mathbf{x}}^{(t)}-{\mathbf{x}}^{*}||$ is non-increasing in $t$ and is decreasing for every $t$ such that $||\nabla f({\mathbf{x}}^{(t)})||\neq 0$.

We next present a new result which shows that the MM algorithm also has linear convergence, for any smooth and strongly convex function $f$ that has a surrogate function $g$ satisfying a certain condition.

Theorem 3.2** **(MM)

Assume $f$ is a convex $\mu$-smooth function on $\mathcal{X}_{\mu}$ satisfying the $\gamma$-PL inequality on $\mathcal{X}_{\gamma}\subseteq\mathcal{X}_{\mu}$, ${\mathbf{x}}^{*}\in{\mathcal{X}}_{\gamma}$ is a minimizer of $f$ and ${\mathbf{x}}^{(t)}\mapsto{\mathbf{x}}^{(t+1)}$ is according to the MM algorithm (2.6). Let $g$ be a majorant surrogate function of $f$ such that for some $\delta>0$,

[TABLE]

Then, if ${\mathbf{x}}^{(t)}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}^{(t)}-\frac{1}{\mu+\delta}\nabla f({\mathbf{x}}^{(t)})\in\mathcal{X}_{\mu}$, there is a $\gamma/(\mu+\delta)$-improvement with respect to $f$ at time step $t$.

The proof of Theorem 3.2 is based on separating the gap between the objective function value at an iteration and the optimum function value in two components: (a) one due to using a surrogate function and (b) other due to a virtual gradient descent update. The remaining arguments are similar to those of the proof of Theorem 3.1. It is worth noting that a different set of sufficient conditions for linear convergence of MM algorithms were found in Proposition 2.7 in Mairal (2015). These conditions require that $g$ is a first-order surrogate function of $f$, a notion we discuss further in Section 7.3.

From Theorems 3.1 and 3.2, we observe that the MM algorithm has the same rate of convergence bound as the gradient descent algorithm except for the smoothness parameter $\mu$ being enlarged for value $\delta$. If $\delta\leq c\mu$, for a constant $c>0$, then the MM algorithm has the same rate of convergence bound as the gradient descent algorithm up to a constant factor.

3.2 Maximum likelihood estimation

We consider the rate of convergence for the ML parameter estimation for the Bradley-Terry model of paired comparisons. This estimation problem amounts to finding a parameter vector that minimizes the negative log-likelihood function, with the log-likelihood function given in (2.3). Recall that $\mathbf{M}$ denotes the matrix of the item-pair co-occurrence counts and $\mathbf{L}_{\mathbf{M}}$ denotes the associated Laplacian matrix. For any positive semidefinite matrix $\mathbf{A}$, we let $\lambda_{i}(\mathbf{A})$ denote the $i$-th smallest eigenvalue of $\mathbf{A}$.

Lemma 3.1

For any $\omega\geq 0$, the negative log-likelihood function for the Bradley-Terry model of paired comparisons is $\gamma$-strongly convex on $\mathcal{W}_{\omega,0}=\mathcal{W}_{\omega}\cap\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:{\mathbf{w}}^{\top}{\mathbf{1}}=0\}$, where $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$, and $\mu$-smooth on $\mathrm{I\!R}^{n}$ with

[TABLE]

where $c_{\omega}=1/(e^{-\omega}+e^{\omega})^{2}$.

By Lemma 3.1, the smoothness parameter $\mu$ is proportional to the largest eigenvalue of the Laplacian matrix ${\mathbf{L}}_{{\mathbf{M}}}$. By the Gershgorin circle theorem, e.g., Theorem 7.2.1 in Golub and Loan (2013), we have $\lambda_{n}({\mathbf{L}}_{{\mathbf{M}}})\leq 2d({\mathbf{M}})$. Thus, we can take $\mu=d({\mathbf{M}})/2$. We will express all our convergence time results in terms of $d({\mathbf{M}})$ instead of $\lambda_{n}({\mathbf{L}}_{{\mathbf{M}}})$. This is a tight characterization up to constant factors. When ${\mathbf{M}}$ is a graph adjacency matrix, then $\lambda_{n}({\mathbf{L}}_{{\mathbf{M}}})\geq d({\mathbf{M}})+1$ by Grone et al. (1990). In the context of paired comparisons, $d({\mathbf{M}})$ has an intuitive interpretation as the maximum number of observed paired comparisons involving an item.

The following lemma will be useful for showing that a function $f$ satisfies the $\gamma$-PL inequality if it satisfies a $\gamma$-strong convexity condition.

Lemma 3.2

Assume that $\mathcal{X}$ is a convex set such that $f$ is $\gamma$-strongly convex on $\mathcal{X}_{0}=\mathcal{X}\cap\{{\mathbf{x}}\in\mathrm{I\!R}^{n}:{\mathbf{x}}^{\top}{\mathbf{1}}=0\}$, and that for all $c\in\mathrm{I\!R}$ and ${\mathbf{x}}\in\mathcal{X}$,

(C1)

$f(\Pi_{c}({\mathbf{x}}))=f({\mathbf{x}})$* and*

(C2)

$\nabla f(\Pi_{c}({\mathbf{x}}))=\nabla f({\mathbf{x}})$**

where

[TABLE]

Then, $f$ satisfies the $\gamma$-PL inequality on $\mathcal{X}$.

The proof of Lemma 3.2 follows by noting that if $f$ is $\gamma$-strongly convex on $\mathcal{X}_{0}$, then it satisfies the $\gamma$-PL inequality on $\mathcal{X}_{0}$. Since for every ${\mathbf{x}}\in\mathcal{X}$, ${\mathbf{x}}={\mathbf{x}}^{\prime}+c{\mathbf{1}}$ for some ${\mathbf{x}}^{\prime}\in\mathcal{X}_{0}$ and $c\in\mathrm{I\!R}$, by conditions (C1) and (C2) and definition of $\gamma$-PL inequality in (2.12), it follows that if $\gamma$-PL inequality holds on $\mathcal{X}_{0}$, it holds as well on $\mathcal{X}$. Conditions (C1) and (C2) are satisfied by negative log-likelihood functions for generalized Bradley-Terry models.

Since the negative log-likelihood function of the Bradley-Terry model of paired comparisons satisfies conditions (C1) and (C2) of Lemma 3.2, combining with Lemma 3.1, we observe that it satisfies the $\gamma$-PL inequality on $\mathcal{W}_{\omega}$ with $\gamma=c_{\omega}a({\mathbf{M}})$. Furthermore, by Lemma 3.1, the negative log-likelihood function is $\mu$-smooth on $\mathrm{I\!R}^{n}$ with $\mu=d({\mathbf{M}})/2$. Combining these facts with Theorem 3.1, we have the following corollary:

Corollary 3.1** **(gradient descent)

Assume that ${\mathbf{w}}^{*}$ is the maximum likelihood parameter estimate in $\mathcal{W}_{\omega}=\{{\mathbf{w}}:\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$, for some $\omega\geq 0$, and ${\mathbf{w}}^{(t)}\mapsto{\mathbf{w}}^{(t+1)}$ is according to gradient descent algorithm with step size $\eta=2/d({\mathbf{M}})$.

Then, if ${\mathbf{w}}^{(t)}\in\mathcal{W}_{\omega}$, there is an $\alpha_{{\mathbf{M}},\omega}$-improvement at time step $t$ where

[TABLE]

The result in Corollary 3.1 implies a linear convergence with the rate of convergence bound $1-2c_{\omega}a({\mathbf{M}})/d({\mathbf{M}})$. Hence, we have the following convergence time bound:

[TABLE]

We next consider the classic MM algorithm for the ML estimation problem, which uses the surrogate function in (2.7). This surrogate function satisfies the following property:

Lemma 3.3

For any $\omega\geq 0$, for all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$, $\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})\geq-\frac{\delta}{2}\|{\mathbf{x}}-{\mathbf{y}}\|^{2}$ where

[TABLE]

By Theorem 3.2 and Lemmas 3.1, 3.2, and 3.3, we have the following corollary:

Corollary 3.2** **(MM)

Assume that ${\mathbf{w}}^{*}$ is the maximum likelihood parameter estimate in $\mathcal{W}_{\omega}=\{{\mathbf{w}}:\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$, for some $\omega\geq 0$, and that ${\mathbf{w}}^{(t)}\mapsto{\mathbf{w}}^{(t+1)}$ is according to the MM algorithm.

Then, if ${\mathbf{w}}^{(t)}\in\mathcal{W}_{\omega}$, there is an $\alpha_{{\mathbf{M}},\omega}$-improvement with respect at time step $t$ where

[TABLE]

and $c^{\prime}_{\omega}=1/[(e^{-\omega}+e^{\omega})^{2}(1+e^{2\omega})]$.

From Corollaries 3.1 and 3.2, we observe that both gradient descent and MM algorithms have the rate of convergence bound of the form $1-ca({\mathbf{M}})/d({\mathbf{M}})$ for some constant $c>0$. The only difference is the value of constant $c$. Hence, both gradient descent and MM algorithm have a linear convergence, and the convergence time bound (3.1).

3.3 Maximum a posteriori probability estimation

We next consider the maximum a posteriori probability estimation problem. We first note that the negative log-a posteriori probability function has the following properties.

Lemma 3.4

The negative log-a posteriori probability function for the Bradley-Terry model of paired comparisons and the prior distribution Gamma$(\alpha,\beta)$ is $\gamma$-strongly convex and $\mu$-smooth on $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$ with

[TABLE]

Note that the strong convexity parameter $\gamma$ is proportional to $\beta$ while as shown in Lemma  3.1, for the ML objective $\gamma$ is proportional to $\lambda_{2}(\mathbf{M})$. This has important implications on the rate of convergence which we discuss next.

By Theorem 3.1 and Lemma 3.4, we have the following corollary:

Corollary 3.3** **(gradient descent)

Assume ${\mathbf{w}}^{*}$ is the maximum a posteriori parameter estimate in $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$, for some $\omega\geq 0$, and ${\mathbf{w}}^{(t)}\mapsto{\mathbf{w}}^{(t+1)}$ is according to gradient descent algorithm (2.10) with step size $\eta=2/(d({\mathbf{M}})+2\beta e^{\omega})$.

Then, if ${\mathbf{w}}^{(t)}\in\mathcal{W}_{\omega}$, there is a $\alpha_{{\mathbf{M}},\omega}$-improvement where

[TABLE]

The result in Corollary 3.3 implies a linear convergence with the convergence time bound

[TABLE]

This bound can be arbitrarily large by taking parameter $\beta$ to be small enough. In Section 3.4, we will show a simple instance for which this bound is tight. Hence, the convergence time for MAP estimation can be arbitrarily slow, and much slower than for the ML case.

We next consider the MM algorithm. First, note that since $\underline{\rho}({\mathbf{x}};{\mathbf{y}})-\rho({\mathbf{x}})=\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})$, by Lemma 3.3, we have

Lemma 3.5

For all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$, $\underline{\rho}({\mathbf{x}};{\mathbf{y}})-\rho({\mathbf{x}})\geq-\frac{\delta}{2}\|{\mathbf{x}}-{\mathbf{y}}\|^{2}$ where $\delta=\frac{1}{2}e^{2\omega}d({\mathbf{M}}).$

By Theorem 3.2 and Lemmas 3.4 and 3.5, we have the following corollary:

Corollary 3.4** **(MM)

Assume ${\mathbf{w}}^{*}$ is the maximum a posteriori parameter estimate in $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$, for some $\omega\geq 0$, and ${\mathbf{w}}^{(t)}\mapsto{\mathbf{w}}^{(t+1)}$ is according to the MM algorithm.

Then, if ${\mathbf{w}}^{(t)}\in\mathcal{W}_{\omega}$, there is an $\alpha_{{\mathbf{M}},\omega}$-improvement at time step $t$ where

[TABLE]

From Corollaries 3.2 and 3.4, we observe that both gradient descent algorithm and MM algorithm have the rate of convergence bound $1-\Omega(\beta/(\beta+d({\mathbf{M}})))$, and hence both have linear convergence and both have the convergence time bound (3.2). We next establish tightness of this rate of convergence by showing that it is achieved for a simple instance.

3.4 Tightness of the rate of convergence

Consider an instance with two items, which are compared $m$ times. Note that $d({\mathbf{M}})=m$ and $a({\mathbf{M}})=2m$. This simple instance allows us to express the parameter vector estimate ${\mathbf{w}}^{(t)}$ and the log-a posteriori probability function $\rho({\mathbf{w}}^{(t)})$ in a closed form, with detailed derivations provided in Section 7.8. From this analysis, we have that $\rho({\mathbf{w}}^{(t)})$ satisfies

[TABLE]

Note that the rate of convergence is approximately $1-4(\alpha-1)/m$ for small $\alpha-1$. By taking $\alpha$ such that $\alpha-1=c\beta$, for a constant $c>0$ such that $||{\mathbf{w}}^{*}||_{\infty}\leq\omega$, we have the rate of convergence $1-\Theta(\beta/d({\mathbf{M}}))$. This establishes the tightness of the rate of convergence bound in Corollary 3.4.

3.5 A simple illustrative numerical example

We illustrate the rate of convergence for a simple example, using randomly generated observations of paired comparisons. This allows us to demonstrate how the number of iterations grows as the value of parameter $\beta$ becomes smaller, and how the number of iterations is affected by the value of parameter $\omega$. Later, in Section 5, we provide further validation by using real-world datasets.

Our example is for an instance with $10$ items with each distinct pair of items compared $10$ times and the input data generated according to the Bradley-Terry model of paired comparisons with the parameter vector such that a half of items have parameter value $-\omega$ and the other half of items have parameter value $\omega$, for a parameter $\omega>0$. We define the convergence time $T$ to be the smallest integer $t$ such that $||{\mathbf{w}}^{(t)}-{\mathbf{w}}^{(t-1)}||_{\infty}\leq\xi$, for a fixed parameter $\xi>0$. In our experiments, we set $\xi=0.0001$.

The results in Figure 1, obtained for $\omega=1/2$, demonstrate that the MM algorithm for the MAP estimation problem with $\beta>0$ can be much slower than the MM algorithm for the ML estimation problem.

We further evaluate the convergence time of gradient descent and MM algorithms for different values of parameter $\omega$, for each distinct pair of items compared $100$ times. The numerical results in Figure 2 show the number of iterations versus the value of parameter $\beta$ for gradient descent and MM algorithms, for different values of parameter $\omega$. We observe that for small enough value of $\omega$, convergence times of gradient descent and MM algorithms are nearly identical. Both algorithms have the convergence time increasing by decreasing the value of $\beta$ for strictly positive values of this parameter. We also observe a discontinuity in convergence time, for $\beta=0$ (MLE case) being smaller than for some small positive value of $\beta$ (MAP case).

The discontinuity at $\beta=0$ originates from the fact that the log-likelihood function has infinitely many solutions for $\beta=0$, but has a unique solution whenever $\beta>0$. Consider a simple illustrative example: $f(x_{1},x_{2})=(x_{1}-1)^{2}+\theta(x_{2}-1)^{2}$, for a parameter $\theta\geq 0$. Then, the gradient descent converges to the unique solution $(1,1)$ slowly when $\theta$ is close to [math]. When $\theta=0$, however, we just need to find the minimum point of $(x_{1}-1)^{2}$ which can be solved in a few iterations.

4 Accelerated MAP inference

In this section, we present a new accelerated algorithm for gradient descent and MM algorithms for MAP estimation. The key element is a transformation of the parameter vector estimate in each iteration of an iterative optimization algorithm that (a) ensures monotonic improvement of the optimization objective along the sequence of parameter vector estimates and (b) ensures certain second-order properties of the objective function hold along the sequence of parameter vector estimates.

We first introduce transformed versions of gradient descent and MM algorithms. Given a mapping $\Pi:\mathrm{I\!R}^{n}\rightarrow\mathrm{I\!R}^{n}$, we define the $\Pi$-transformed gradient descent algorithm by

[TABLE]

Similarly, we define the $\Pi$-transformed MM algorithm by the iteration:

[TABLE]

Importantly, function $f$ and mapping $\Pi$ have to satisfy certain conditions in order to provide a convergence rate guarantee, which we discuss in the following section.

4.1 General convergence theorems

Assume that $f$ and $\Pi$ satisfy the following conditions, for a convex set $\mathcal{X}$ that contains optimum point ${\mathbf{x}}^{*}$, and a vector ${\mathbf{d}}\in\mathrm{I\!R}^{n}$:

(F1)

$f$ is $\mu$-smooth on $\mathcal{X}$;

(F2)

$f$ satisfies the $\gamma$-PL inequality on

[TABLE]

and

(P1)

$f(\Pi({\mathbf{x}}))\leq f({\mathbf{x}})$ for all ${\mathbf{x}}\in\mathrm{I\!R}^{n}$;

(P2)

$\Pi({\mathbf{x}})\in\mathcal{X}_{0}$ when $\Pi({\mathbf{x}})\in\mathcal{X}$.

Condition (F1) is a standard smoothness condition imposed on $\mathcal{X}$. Condition (F2) is a standard $\gamma$-PL condition imposed on the subset of points in $\mathcal{X}$ at which the gradient of the function $f$ is orthogonal to vector ${\mathbf{d}}$. Condition (P1) means that applying $\Pi$ to a point cannot increase the value of function $f$. This condition is crucial to ensure a monotonic improvement of the objective function value when transformation $\Pi$ is applied to an iterative optimization method. Condition (P2) is satisfied when at any $\Pi$-transformed point, the gradient of $f$ is orthogonal to vector ${\mathbf{d}}$. This condition is crucial to ensure certain second-order properties hold when $\Pi$ is applied to an iterative optimization method.

We have the following two theorems.

Theorem 4.1** **(Gradient descent)

Assume that $f$ satisfies (F1) and (F2), $\Pi$ satisfies (P1), and $\eta=1/\mu$. Let ${\mathbf{x}}^{(t)}\mapsto{\mathbf{x}}^{(t+1)}$ be according to the $\Pi$-transformed gradient descent algorithm (4.1).

Then, if ${\mathbf{x}}^{(t)},{\mathbf{x}}^{(t+1)}\in\mathcal{X}_{0}$, there is an $\gamma/\mu$-improvement with respect to $f$ at time step $t$.

The theorem establishes the same rate of convergence as for the gradient descent algorithm in Theorem 3.1 but with the strong convexity condition restricted to points at which the gradient of $f$ is orthogonal to vector ${\mathbf{d}}$.

Theorem 4.2** **(MM)

Assume that $f$ satisfies (F1) and (F2), $\Pi$ satisfies (P1), and $g$ is a majorant surrogate function of $f$ such that $g({\mathbf{x}};{\mathbf{y}})-f({\mathbf{x}})\leq\frac{\delta}{2}||{\mathbf{x}}-{\mathbf{y}}||^{2}$. Let ${\mathbf{x}}^{(t)}\mapsto{\mathbf{x}}^{(t+1)}$ be according to the $\Pi$-transformed MM algorithm (4.2).

Then, if ${\mathbf{x}}^{(t)},{\mathbf{x}}^{(t+1)}\in\mathcal{X}_{0}$, there is an $\gamma/(\gamma+\delta)$-improvement with respect to $f$ at time step $t$.

The last theorem establishes the same rate of convergence as for classic MM algorithm in Theorem 3.2, but with a strong convexity condition imposed only at the points at which the gradient of $f$ is orthogonal to vector ${\mathbf{d}}$.

We next present a lemma which will be instrumental in showing that the PL condition in (F1) holds for the MAP estimation problem.

Lemma 4.1

Assume $f$ is a convex, twice-differentiable function. Let $\mathcal{X}$ be a convex set, $\mathcal{X}_{0}$ be defined by (4.3) for a given vector ${\mathbf{d}}$, and ${\mathbf{x}}^{*}\in\mathcal{X}_{0}$ be a minimizer of $f$.

If for some positive semidefinite ${\mathbf{A}}_{\mathcal{X}}$ such that

(A1)

$\nabla^{2}f({\mathbf{x}})\succeq{\mathbf{A}}_{\mathcal{X}}$* for all ${\mathbf{x}}\in\mathcal{X}$, and*

(A2)

${\mathbf{u}}^{\top}{\mathbf{A}}_{\mathcal{X}}{\mathbf{v}}=0$* for all ${\mathbf{u}},{\mathbf{v}}$ such that*

[TABLE]

where

[TABLE]

then, $f$ satisfies the $\gamma$-PL inequality on $\mathcal{X}_{0}$ for all $\gamma\leq\gamma_{0}$ with

[TABLE]

Note that ${\mathbf{P}}_{\mathbf{d}}$ is the projection matrix onto the space orthogonal to vector ${\mathbf{d}}$. The value of $\gamma_{0}$ is maximized when vector ${\mathbf{d}}$ is the eigenvector corresponding to the smallest eigenvalue of ${\mathbf{A}}_{\mathcal{X}}$. In this case, $\gamma_{0}$ is the second smallest eigenvalue of ${\mathbf{A}}_{\mathcal{X}}$.

4.2 Convergence rate for the BT model

In this section, we apply the framework developed in the previous section to characterize the convergence rate for the MAP parameter estimation of the Bradley-Terry model of paired comparisons. The MAP parameter estimation problem amounts to finding a parameter vector that maximizes the log-a posteriori probability function $\rho$ defined in (2.4). Let the transformation $\Pi$ be defined as

[TABLE]

where

[TABLE]

We next show that $f$ and $\Pi$ satisfy conditions (F1), (F2), (P1), and (P2) for the direction vector ${\mathbf{d}}={\mathbf{1}}$. This will allow us to apply Theorems 4.1 and 4.2 to characterize the rate of convergence for $\Pi$-transformed gradient descent and MM algorithms.

We first show that $f$ satisfies conditions (F1) and (F2) for the set $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\}$. Condition (F1) holds because, in Lemma 3.4, we have already shown that $f$ is $\mu$-smooth on $\mathcal{W}_{\omega}$ with $\mu=d({\mathbf{M}})/2+e^{\omega}\beta$. Condition (F2) can be shown to hold by Lemma 4.1 as follows. Note that we have $\nabla^{2}(f({\mathbf{w}}))\succeq{\mathbf{A}}_{\mathcal{W}_{\omega}}$, for all ${\mathbf{w}}\in\mathcal{W}_{\omega}$, where ${\mathbf{A}}_{\mathcal{W}_{\omega}}=c_{\omega}\mathbf{L}_{\mathbf{M}}+e^{-\omega}\beta{\mathbf{I}}$. The assumptions of Lemma 4.1 hold: (A1) holds because ${\mathbf{A}}_{\mathcal{W}_{\omega}}$ is a positive semidefinite matrix, and (A2) holds because ${\mathbf{u}}^{\top}{\mathbf{L}}_{{\mathbf{M}}}{\mathbf{v}}=0$ (which follows from ${\mathbf{L}}_{{\mathbf{M}}}{\mathbf{1}}={\mathbf{0}}$) and ${\mathbf{u}}^{\top}{\mathbf{I}}{\mathbf{v}}={\mathbf{u}}^{\top}{\mathbf{v}}=0$ (${\mathbf{u}}$ and ${\mathbf{v}}$ are orthogonal). Since $\gamma_{0}$ is the smallest eigenvalue of ${\mathbf{A}}_{\mathcal{X}}$ on the subspace orthogonal to vector ${\mathbf{1}}$, we have

[TABLE]

Hence, by Lemma 4.1, it follows that $f$ satisfies condition (F2) with $\gamma=c_{\omega}\lambda_{2}({\mathbf{L}}_{{\mathbf{M}}})+e^{-\omega}\beta$.

We next show that $\Pi$, defined in (4.4), satisfies conditions (P1) and (P2). These two conditions are shown to hold in the following lemma.

Lemma 4.2

For all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$,

[TABLE]

and

[TABLE]

From Theorem 4.1, we have the following corollary:

Corollary 4.1** **(Gradient descent)

Assume that ${\mathbf{w}}^{(t)}\mapsto{\mathbf{w}}^{(t+1)}$ is according to the $\Pi$-transformed gradient descent (4.1) for the negative log-a posteriori probability function of the Bradley-Terry model of paired comparisons, with product-form prior distribution such that $e^{w_{i}}\sim$ Gamma ($\alpha,\beta$), $\alpha\geq 1$ and $\beta>0$, and $\eta=2/(||\mathbf{M}||_{\infty}+2e^{\omega}\beta)$.

Then, there is an $\alpha_{{\mathbf{M}},\omega,\beta}$-improvement with respect to $\rho$ at time step $t$ where

[TABLE]

and $c_{\omega}=1/(e^{-\omega}+e^{\omega})^{2}$.

From Theorem 4.2, we have the following corollary:

Corollary 4.2** **(MM)

Assume that iterates are according to the $\Pi$-transformed MM (4.1) for the negative log-a posteriori probability function of the Bradley-Terry model of paired comparisons, with product-form prior distribution such that $e^{w_{i}}\sim$ Gamma ($\alpha,\beta$), $\alpha\geq 1$ and $\beta>0$.

Then, there is an $\alpha_{{\mathbf{M}},\omega,\beta}$-improvement with respect to $\rho$ at time step $t$ where

[TABLE]

and $c_{\omega}=1/(e^{-\omega}+e^{\omega})^{2}$.

Proof.

Proof Condition (F1) holds for $-\rho$ because we have already shown that $-\rho$ is $\mu$-smooth with $\mu=d({\mathbf{M}})/2+e^{\omega}\beta$ on $\mathcal{W}_{\omega}$ and $\delta=e^{2\omega}d({\mathbf{M}})/2$. Condition (F2) holds by Lemma 4.1 with $\gamma=c_{\omega}a({\mathbf{M}})+e^{-\omega}\beta$. Conditions (P1) and (P2) hold by (4.6) and (4.7) in Lemma 4.7, respectively. \Halmos∎

Note that in the limit of small $\beta$, the convergence rate bounds in Corollaries 4.1 and 4.2 correspond to the bounds for the ML estimation in Corollaries 3.1 and 3.2, respectively. From Corollaries 4.1 and 4.2, it follows that for accelerated gradient descent and accelerated MM algorithm, the convergence time satisfies

[TABLE]

For the Bradley-Terry model of paired comparisons with parametrization $\theta=(\theta_{1},\ldots,\theta_{n})^{\top}$, where $\theta_{i}=e^{w_{i}}$ for $i=1,2,\ldots,n$, the transformation $\Pi$ given by (4.4) is equivalent to a rescaling as shown in a procedural form in Algorithm 1. This algorithm first performs the standard MM update in Eq. (2.9), which is followed by rescaling the resulting intermediate parameter vector such that the parameter vector $\theta$ at every iteration satisfies $\sum_{i=1}^{n}\theta_{i}=c$ where $c=n(\alpha-1)/\beta$. This can be interpreted as fixing the scale of parameters to a carefully chosen scale that is dependent on the choice of the prior distribution. Note that the scaling factor $c$ cannot be arbitrarily fixed while still preserving good convergence properties. In particular, selecting the scale $c=1$ can result in undesired convergence properties. We demonstrate this in Section 5 through numerical examples.

It turns out that the rescaling in Algorithm 1 is roughly of the same order as the random rescaling suggested in Caron and Doucet (2012). Therein, the authors suggested using independent identically distributed rescaling factors across different iteration steps with distribution Gamma($n\alpha,\beta$). This ensures the distribution of $\sum_{i=1}^{n}\theta_{i}$ to remain invariant across different iterations, equal to Gamma($n\alpha,\beta$). The mode of this rescaling factor is $(n\alpha-1)/\beta$. This bears a similarity with the rescaling in Algorithm 1, in particular, with respect to the dependence on parameter $\beta$. Our results show that it suffices to use a simple deterministic rescaling factor to ensure linear convergence. Moreover, using a different rescaling than the one used in Algorithm 1 can result in a lack or slow convergence, which is shown by numerical experiments in Section 5.

Our numerical example revisited

We ran the accelerated MM algorithm for our numerical example and obtained the results shown in Figure 3. By comparing with the corresponding results obtained by the MM algorithm with no acceleration, shown in Figure 1, we observe that the acceleration resolves the slow convergence issue and that it can yield a significant reduction of the convergence time.

5 Numerical results

In this section we present evaluation of convergence times of MM algorithms for different generalized Bradley-Terry models for a collection of real-world datasets. Our goal is to provide empirical validation of some of the hypotheses derived from our theoretical analysis. Overall, our numerical results validate that (a) the convergence of the MM algorithm for MAP estimation can be much slower than for ML estimation, (b) MM algorithm for MAP estimation has convergence time that increases as parameter $\beta$ of the prior distribution decreases, and (c) a significant reduction of the convergence time can be achieved by the accelerated MM algorithm defined in Section 4. The code and datasets for reproducing our experiments are available online at:

https://github.com/GDMMBT/AcceleratedBradleyTerry.

5.1 Datasets

We consider three datasets, which vary in the type of data, size and sparsity. The three datasets are described as follows.

GIFGIF

This dataset contains user evaluations of digital images by paired comparisons with respect to different metrics, such as amusement, content, and happiness. The dataset was collected through an online web service by the MIT Media Lab as part of the PlacePulse project Rich et al. (2018). This service presents the user with a pair of images and asks to select one that better expresses a given metric, or select none. The dataset contains 1,048,576 observations and covers 17 metrics. We used this dataset to evaluate convergence of MM algorithms for the Bradley-Terry model of paired comparisons. We did this for each of the three aforementioned metrics.

Chess

This dataset contains game-by-game results for 65,030 matches among 8,631 chess players. The dataset was used in a Kaggle chess ratings competition Sonas (2010). Each observation contains information for a match between two players including unique identifiers of the two players, information about which one of the two players played with white figures, and the result of the match, which is either win, loss, or draw. This dataset has a large degree of sparsity. We used this dataset to evaluate convergence of the Rao-Kupper model of paired comparisons with ties.

NASCAR

This dataset contains auto racing competition results. Each observation is for an auto race, consisting of the ranking of drivers in increasing order of their race finish times. The dataset is available from a web page maintained by Hunter (2003). This dataset was previously used for evaluation of MM algorithms for the Plackett-Luce ranking model by Hunter (2004) and more recently by Caron and Doucet (2012). We used this dataset to evaluate convergence times of MM algorithms for the Plackett-Luce ranking model.

We summarise some key statistics for each dataset in Table 1. We use the shorthand notation GIFGIF: A, GIFGIF: C, and GIFGIF: H to denote datasets for metrics amusement, contempt, and happiness, respectively. For full GIFGIF and Chess datasets, we can split the items into two groups such that at least one item in one group is not compared with any item in the other group, i.e., the algebraic connectivity $a({\mathbf{M}})$ of matrix ${\mathbf{M}}$ is zero. In this case, there exists no ML estimate, while an MAP estimate always exists. In order to consider cases when an MLE exists, we consider sampled datasets by restricting to the set of items such that the algebraic connectivity for this subset of items is strictly positive. This subsampling was done by selecting the largest connected component of items.

5.2 Experimental results

We evaluated the convergence time defined as the number of iterations that an algorithm takes until a convergence criteria is satisfied. We use the standard convergence criteria based on the difference of successive parameter vector estimates. Specifically, the convergence time $T$ is defined as the smallest integer $t>0$ such that $||{\mathbf{w}}^{(t)}-{\mathbf{w}}^{(t-1)}||_{\infty}\leq\xi$, for fixed value of parameter $\xi>0$, with initial value ${\mathbf{w}}^{(0)}={\mathbf{0}}$. In our experiments, we used $10^{-4}$ as the default value for parameter $\xi$. For NASCAR dataset, we also present results for several other values of $\xi$ to demonstrate how the convergence time changes. In our experiments, we also evaluated the convergence time measured in real processor time units. We noted that they validate all the observations derived from the convergence times measured in the number of iterations, and hence we do not further discuss them.

In our experiments we varied the value of parameter $\beta$ and, unless specified otherwise, we set the value of parameter $\alpha$ such that $\alpha-1=\beta$. This corresponds to fixing the mode of the Gamma prior marginal distributions to value $1$. Note that the case $\beta=0$ corresponds to ML estimation.

Before discussing numerical convergence time results, we first show results validating that the MM algorithm converges and that this convergence is linear. This is shown in Figure 4 for GIFGIF A dataset for three different values of parameter $\beta$. For space reasons, we only include results for this dataset. We observe that in all cases the log-a posteriori probability monotonically increases with the number of iterations, thus validating convergence. We also observe that the gap between the maximum log-a posteriori probability and the log-a posteriori probability decreases with the number of iterations in a linear fashion for sufficiently large number of iterations, when plotted using the logarithmic scale for the $y$ axis, thus validating linear convergence.

We next discuss our numerical results for convergence time evaluated for the MM algorithm and accelerated MM algorithm for different datasets and choice of parameters. Our numerical results are shown in Table 2.

For GIFGIF datasets, we observe that the convergence time increases as the value of parameter $\beta$ decreases for $\beta>0$. For the values of $\beta$ considered, this increase can be for as much as two orders of magnitude. When the ML estimate exists (for sampled data), we observe that the MM algorithm for ML estimation converges much faster than the MM algorithm for MAP estimation for sufficiently small values of parameter $\beta$. We also observe that a significant reduction of the convergence time can be achieved by the accelerated MM algorithm. This reduction can be for as much as order $10\%$ of the convergence time of the MM algorithm without acceleration. These empirical results validate our theoretical results.

For Chess datasets, all the observations derived by using the GIFGIF datasets remain to hold.

For NASCAR dataset, we show results for different values of parameter $\xi$, including the default value of $10^{-4}$. Again, all the observations made for GIFGIF and Chess datasets remain to hold. It is noteworthy that the MM algorithm for ML estimation converges much faster than for MAP estimation for sufficiently small values of parameter $\beta$. This is especially pronounced for smaller values of $\xi$. For the cases considered, this can be for as much as three orders of magnitude. Similarly, the accelerated MM algorithm converges much faster than the classical MM algorithm.

In the reminder of this section we highlight the importance of carefully changing the scale of the parameter vector in each iteration, as done in our accelerated MM algorithm, Algorithm 1, as otherwise the monotonic convergence may not be guaranteed or the convergence may be slow. To demonstrate this, we examine the alternative change of scale such that the parameter vector ${\mathbf{w}}$ in each iteration satisfies $\sum_{i=1}^{n}e^{w_{i}}=1$. We present the results for the GIFGIF (sample) datasets. From Figure 5, we observe that the algorithm does not guarantee a monotonic increase of the a posteriori probability with the number of iterations, which is unlike to our accelerated MM algorithm for which this always holds. In Table 3, we show the same quantities as in Table 2 but for the MM algorithm with the alternative change of scale under consideration. We observe that our acceleration method can converge much faster, and that there are cases for which the alternative change of scale results in no convergence within a bound on the maximum number of iterations.

6 Further discussion

We have shown that for generalized Bradley-Terry models, gradient descent and MM algorithms for the ML estimation problem have a linear convergence with the convergence time bound $O(d({\mathbf{M}})/a({\mathbf{M}}))$, where $d({\mathbf{M}})$ is the maximum number of observed comparisons per item and $a({\mathbf{M}})$ is the algebraic connectivity of the matrix $\mathbf{M}$ of the observed counts of item-pair co-occurrences. We have also shown that for generalized Bradley-Terry models, gradient descent and MM algorithms for the MAP estimation problem, with the prior product-form distribution with Gamma$(\alpha,\beta)$ marginal distributions, the convergence time is also linear but with the convergence time bound $O(d({\mathbf{M}})/\beta)$. This bound is shown to be tight. Our results identify a slow convergence issue for gradient descent and MM algorithms for the MAP estimation problem, which occurs for small values of parameter $\beta$. The small values of parameter $\beta$ correspond to more vague prior distributions. Our results identify a discontinuity of the convergence time at $(\alpha,\beta)=(1,0)$, which corresponds to ML estimation. The proposed acceleration method for the MAP estimation problem resolves the slow convergence issue, and yields a convergence time that is bounded by the best of what can be achieved for the ML and MAP estimation problems.

Our results provide insights into how the observed comparison data affect the rate of convergence of gradient descent and MM algorithms. The two key parameters affecting the rate of convergence are $d({\mathbf{M}})$ and $a({\mathbf{M}})$. For illustration purposes, in Table 4 we show values of $d({\mathbf{M}})$ and $a({\mathbf{M}})$ for examples of matrix ${\mathbf{M}}$ with [math]-$1$ valued entries, which correspond to graph adjacency matrices. We observe that when each distinct pair is compared the same number of times, i.e. for the complete graph case, the convergence time is $T=O(\log(1/\epsilon))$. For other cases, the convergence time is $T=O(n^{c}\log(1/\epsilon))$, for some $c\geq 1$.

We further consider the case of random design matrices where each distinct pair of items is either compared once or not compared at all, and this is according to independent Bernoulli random variables with parameter $p$ across all distinct pairs of items. In other words, the item pair co-occurrence is according to the Erdös-Rényi random graph $G_{n,p}$ and ${\mathbf{M}}$ is its adjacency matrix. $d({\mathbf{M}})$ corresponds to the maximum degree of $G_{n,p}$ which has been extensively studied, with precise results obtained for different scalings of $p$ with $n$. In particular, by Bollobás (2001) (Corollary 3.4), $d({\mathbf{M}})=pn+O(\sqrt{pn\log(n)})$ with probability $1-1/n$ provided that $p=\omega(\log(n)^{3}/n)$. The algebraic connectivity for Erdös-Rényi graphs has been studied as well. By Coja-Oghlan (2007) (Theorem 1.3), $a({\mathbf{M}})=pn+O(\sqrt{pn\log(n)})$ with probability $1-o(1)$, provided that $p=\omega(\log(n)^{2}/n)$. Intuitively, if the expected degree $np$ is large enough, $d({\mathbf{M}})/a({\mathbf{M}})=\Theta(1)$.

We can derive an upper bound for the convergence time, which depends only on some simple properties of the graph associated with matrix ${\mathbf{M}}$. Let ${\mathbf{A}}$ be the adjacency matrix of a graph $G$ which has edge $(i,j)$ if, and only if, $m_{i,j}>0$. Let $r=\overline{m}/\underline{m}$ where $\overline{m}=\max_{i,j}m_{i,j}$ and $\underline{m}=\min\{m_{i,j}:m_{i,j}>0\}$. Let $d(n)$ be the maximum degree and $D(n)$ be the diameter of $G$. Then, for both gradient descent and MM algorithms for the ML estimation, we have the convergence time bound (shown in Section 7.14):

[TABLE]

This implies that $T=O(rn^{3}\log(1/\epsilon))$ for every connected graph $G$, which follows by trivial facts $d(n)\leq n$ and $D(n)\leq n$. The bound in (6.1) follows from the lower bound on the algebraic connectivity of a Laplacian matrix $\lambda_{2}({\mathbf{L}}_{{\mathbf{A}}})\geq 4/(nD(n))$, see Theorem 3.4 in Merris (1994).

7 Proofs and Additional Results

7.1 Proof of Theorem 3.1

Let ${\mathbf{x}}^{\prime}$ be the output of the gradient descent iteration update for input ${\mathbf{x}}$ with step size $\eta$.

If ${\mathbf{x}}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}^{\prime}\in{\mathcal{X}}_{\mu}$, then

[TABLE]

where the first inequality is by the assumption that $f$ is $\mu$-smooth on $\mathcal{X}_{\mu}$ and the second inequality is by the assumption that $f$ satisfies the $\gamma$-PL inequality on $\mathcal{X}_{\gamma}$. Taking $\eta=1/\mu$, which minimizes the above bound, establishes the claim of the theorem.

7.2 Proof of Theorem 3.2

Let ${\mathbf{x}}^{\prime}$ be the output of the MM algorithm iteration update for input ${\mathbf{x}}$.

By the facts $f({\mathbf{x}}^{\prime})\leq g({\mathbf{x}}^{\prime};{\mathbf{x}})$ and $g({\mathbf{x}}^{\prime};{\mathbf{x}})\leq g({\mathbf{z}};{\mathbf{x}})$ for all ${\mathbf{z}}$, for any $\eta\geq 0$,

[TABLE]

Now, by the same arguments as in the proof of Theorem 3.1, if ${\mathbf{x}}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}-\eta\nabla f({\mathbf{x}})\in\mathcal{X}_{\mu}$, we have

[TABLE]

Next, if ${\mathbf{x}}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}-\eta\nabla f({\mathbf{x}})\in\mathcal{X}_{\mu}$,

[TABLE]

where the first inequality is by the smoothness condition on the majorant surrogate function and the second inequality is by the assumption that $f$ satisfies the PL inequality with parameter $\gamma$ on $\mathcal{X}_{\gamma}$.

Putting the pieces together, we have

[TABLE]

Taking $\eta=1/(\mu+\delta)$ (which minimizes the factor involving $\eta$ in the last inequality) yields the asserted result.

7.3 Comparison of Theorem 3.2 with Proposition 2.7 in Mairal (2015)

Theorem 7.1* *(Proposition 2.7 in Mairal (2015))

Suppose that $f$ is a strongly convex function on $\mathcal{X}_{\gamma}$ and ${\mathbf{x}}^{*}$ is a minimizer of $f$ and that it holds ${\mathbf{x}}^{*}\in\mathcal{X}_{\gamma}$. Assume that $g$ is a first-order surrogate function of $f$ on $\mathcal{X}_{\mu}$ with parameter $\mu_{0}>0$. Let ${\mathbf{x}}^{(t+1)}$ be the output of the MM algorithm for input ${\mathbf{x}}^{(t)}$. Then, if ${\mathbf{x}}^{(t)}\in\mathcal{X}_{\gamma}$ and ${\mathbf{x}}^{(t+1)}\in\mathcal{X}_{\mu}$, then we have

[TABLE]

where

[TABLE]

Proof.

Proof If $g$ is a first-order surrogate function on $\mathcal{X}_{\mu}$ with parameter $\mu_{0}$, then

[TABLE]

where ${\mathbf{x}}^{\prime}=\arg\min_{{\mathbf{z}}^{\prime}}g({\mathbf{z}}^{\prime};{\mathbf{y}})$.

From this, it follows that

[TABLE]

where the last inequality is by convexity of $f$.

We have established the following inequality

[TABLE]

By assumption that $f$ is $\gamma$-strongly convex on $\mathcal{X}_{\gamma}$ and ${\mathbf{x}}\in\mathcal{X}_{\gamma}$, we have

[TABLE]

It follows that

[TABLE]

It remains only to note that

[TABLE]

\Halmos

∎

The rate of convergence bound derived from Theorem 3.2 can be tighter than the rate of convergence bound derived from Theorem 7.1.

To show this consider the Bradley-Terry model for which we have shown in Lemma 3.3 that the surrogate function $\underline{\ell}$ of the log-likelihood function $\ell$ satisfies condition of Theorem 3.2 on $[-\omega,\omega]^{n}$ with $\delta=\frac{1}{2}e^{2\omega}d({\mathbf{M}})$. It also holds that surrogate function $\underline{\ell}$ is also a first-order surrogate function of $\ell$ on $[-\omega,\omega]^{n}$ with $\mu_{0}=\frac{1}{2}e^{2\omega}d({\mathbf{M}})$. Hence in this case, we have $\delta=\mu_{0}$.

The convergence rate bound of Theorem 3.2 is tighter than the convergence rate bound of Theorem 7.1 if and only if $\mu+\delta<4\mu_{0}$. Since $\delta=\mu_{0}$, this is equivalent to $\mu<3\delta$. Since by Lemma 3.1 we can take $\mu=\frac{1}{2}d({\mathbf{M}})$, the latter condition reads as

[TABLE]

which indeed holds true.

7.4 Proof of Lemma 3.1

The Hessian of the negative log-likelihood function has the following elements:

[TABLE]

We will show that for all $i\neq j$,

[TABLE]

and

[TABLE]

From (7.2), we have $\nabla^{2}(-\ell({\mathbf{w}}))\succeq c_{\omega}\mathbf{L}_{\mathbf{M}}$ for all ${\mathbf{w}}\in[-\omega,\omega]^{n}$. Hence, for all ${\mathbf{w}}\in[-\omega,\omega]^{n}$ and ${\mathbf{x}}\in\mathcal{X}$,

[TABLE]

where $\mathcal{X}=\{{\mathbf{x}}\in\mathrm{I\!R}^{n}:{\mathbf{x}}^{\top}{\mathbf{1}}=0\}$. This shows that $-\ell$ is $c_{\omega}\lambda_{2}(\mathbf{L}_{\mathbf{M}})$-strongly convex on $\mathcal{X}$.

From (7.3), we have $\frac{1}{4}\mathbf{L}_{\mathbf{M}}\succeq\nabla^{2}(-\ell({\mathbf{w}}))$ for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$. Hence, for all ${\mathbf{x}}\in\mathrm{I\!R}^{n}$,

[TABLE]

This shows that $-\ell$ is $\frac{1}{4}\lambda_{n}(\mathbf{L}_{\mathbf{M}})$-smooth on $\mathrm{I\!R}^{n}$.

It remains to show that (7.2) and (7.3) hold. For (7.2), we need to show that $c_{\omega}\leq x_{i}x_{j}/(x_{i}+x_{j})^{2}$ for all ${\mathbf{x}}\in[-\omega,\omega]^{n}$. Note that $x_{i}x_{j}/(x_{i}+x_{j})^{2}=z(1-z)$ where $z:=x_{i}/(x_{i}+x_{j})$. Note that $z\in\Omega:=[e^{-\omega}/(e^{-\omega}+e^{\omega}),1-e^{-\omega}/(e^{-\omega}+e^{\omega})]$ for all ${\mathbf{x}}\in[-\omega,\omega]^{n}$. The function $z(1-z)$ achieves its minimum over the interval $\Omega$ at a boundary of $\Omega$. Thus, it holds $\min_{z\in\Omega}z(1-z)=c_{\omega}$. For (7.3), we can immediately note that for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$,

[TABLE]

7.5 Proof of Lemma 3.3

Let ${\mathbf{y}}$ be an arbitrary vector in $[-\omega,\omega]^{n}$. Let $r({\mathbf{x}};{\mathbf{y}})=\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})$ for ${\mathbf{x}}\in[-\omega,\omega]^{n}$. Then, we have

[TABLE]

where $A$ is a $n\times n$ diagonal matrix with diagonal elements

[TABLE]

Since $\nabla^{2}(-\ell({\mathbf{x}}))$ is a positive semi-definite matrix and $A$ is a diagonal matrix, for all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$ and ${\mathbf{w}}\in[-\omega,\omega]^{n}$, we have

[TABLE]

By limited Taylor expansion, for all ${\mathbf{x}}\in[-\omega,\omega]^{n}$,

[TABLE]

By the definition of $r({\mathbf{x}};{\mathbf{y}})$, we have $\bar{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})\geq-\frac{\delta}{2}||{\mathbf{x}}-{\mathbf{y}}||^{2}$.

7.6 Surrogate function (2.7) for the Bradley-Terry model is a first-order surrogate function

We show that the surrogate function $\underline{\ell}$ of the log-likelihood function $\ell$ of the Bradley-Terry model, given by (2.7), is a first-order surrogate function on $\mathcal{X}_{\omega}=[-\omega,\omega]^{n}$ with $\mu_{0}=\frac{1}{2}e^{2\omega}d({\mathbf{M}})$.

We need to show that the error function $h({\mathbf{x}};{\mathbf{y}})=\ell({\mathbf{x}})-\underline{\ell}({\mathbf{x}};{\mathbf{y}})$ is a $\mu_{0}$-smooth function on $\mathcal{X}_{\omega}$.

By a straightforward calculus, we note

[TABLE]

where $D({\mathbf{x}},{\mathbf{y}})$ is a diagonal matrix with diagonal elements

[TABLE]

We can take

[TABLE]

For any $A=B+D$ where $B$ is a $n\times n$ matrix and $D$ is a $n\times n$ diagonal matrix with diagonal elements $d_{1},d_{2},\ldots,d_{n}$, we have

[TABLE]

It thus follows that

[TABLE]

Now note that for all ${\mathbf{x}},{\mathbf{y}}\in\mathcal{X}_{\omega}$,

[TABLE]

and

[TABLE]

We have

[TABLE]

and

[TABLE]

7.7 Proof of Lemma 3.4

We consider the log-a posteriori probability function $\rho({\mathbf{w}})=\ell({\mathbf{w}})+\ell_{0}({\mathbf{w}})+\mathrm{const}$ where $\ell$ is the log-likelihood function given by (2.3) and $\ell_{0}$ is the prior log-likelihood function given by (2.5). Note that $\nabla^{2}(-\ell_{0}({\mathbf{w}}))$ is a diagonal matrix with diagonal elements equal to $\beta e^{w_{i}}$, for $i=1,2,\ldots,n$. It can be readily shown that for ${\mathbf{w}}\in{\mathcal{W}}_{\omega}$,

[TABLE]

From (7.5), for all ${\mathbf{w}}\in{\mathcal{W}}_{\omega}$ and ${\mathbf{x}}\in\mathrm{I\!R}^{n}$,

[TABLE]

Hence, $-\rho$ is $e^{-\omega}\beta$-strongly convex on ${\mathcal{W}}_{\omega}$.

Similarly, from (7.5), for all ${\mathbf{w}}\in{\mathcal{W}}_{\omega}$, and ${\mathbf{x}}\in\mathrm{I\!R}^{n}$,

[TABLE]

Hence, $-\rho$ is $\mu$-smooth on ${\mathcal{W}}_{\omega}$ with $\mu=\frac{1}{4}\lambda_{n}(\mathbf{L}_{\mathbf{M}})+e^{\omega}\beta$.

7.8 The asymptote in Section 3.4

We consider the case of two items, compared $m$ times. Suppose that the observed data is such that the number of comparisons won by items $1$ and $2$ are $d_{1}$ and $d_{2}$, respectively.

The MM algorithm iterates ${\mathbf{w}}^{(t)}$ are such that

[TABLE]

where $s^{(t)}=e^{w_{1}^{(t)}}+e^{w_{2}^{(t)}}$. Observe that $s^{(t)}$ evolves according to the following autonomous nonlinear dynamical system:

[TABLE]

The limit point of $s^{(t)}$ as $t$ goes to infinity is $2(\alpha-1)/\beta$. Note that $(\alpha-1)/\beta$ is the mode of Gamma$(\alpha,\beta)$.

Let us define $a^{(t)}$ by $s^{(t)}=[2(\alpha-1)/\beta](1+a^{(t)})$. Note that $a^{(t)}$ goes to [math] as $t$ goes to infinity. By a tedious but straightforward calculus, we can show that

[TABLE]

From (7.6), note that $1/a^{(t)}$ evolves according to a linear dynamical system, which allows us to derive the solution for $a^{(t)}$ in the explicit form given as follows:

[TABLE]

From (7.7), $\rho({\mathbf{w}}^{*})-\rho({\mathbf{w}}^{(t)})=(\alpha-1)(a^{(t)})^{2}(1+o(1))$ for large $t$, and thus

[TABLE]

It follows that the rate of convergence of the log-a posteriori probability function is given as follows:

[TABLE]

7.9 Proof of Theorem 4.1

Since $f(\Pi({\mathbf{x}}))\leq f({\mathbf{x}})$ for all ${\mathbf{x}}\in\mathrm{I\!R}^{n}$,

[TABLE]

By the same steps as those in the proof of Theorem 3.1, we can show that

[TABLE]

Hence, it follows that

[TABLE]

7.10 Proof of Lemma 4.1

By a limited Taylor expansion, for any ${\mathbf{x}},{\mathbf{y}}\in\mathrm{I\!R}^{n}$, we have

[TABLE]

Let

[TABLE]

where

[TABLE]

Notice that

1. (i)

${\mathbf{u}}+{\mathbf{v}}={\mathbf{y}}-{\mathbf{x}}$, and 2. (ii)

${\mathbf{u}}$ and ${\mathbf{v}}$ are orthogonal, i.e., ${\mathbf{u}}^{\top}{\mathbf{v}}=0$.

From now on, assume that ${\mathbf{x}}$ and ${\mathbf{y}}$ are such that ${\mathbf{x}},{\mathbf{y}}\in\mathcal{X}_{0}$ and ${\mathbf{y}}={\mathbf{x}}^{*}$.

By definition of $\mathcal{X}_{0}$, we have ${\mathbf{d}}^{\top}\nabla f({\mathbf{x}})=0$, which together with ${\mathbf{u}}+{\mathbf{v}}={\mathbf{x}}^{*}-{\mathbf{x}}$, implies

[TABLE]

Now, note that for any $a\in[0,1]$, we have the following relations:

[TABLE]

where $(a)$ is by assumption (A1) and $(b)$ is by assumption that ${\mathbf{A}}_{\mathcal{X}}$ is a positive semidefinite matrix and (A2). Hence, we have shown that, for all $a\in[0,1]$,

[TABLE]

Next, note that

[TABLE]

Combining with (7.8)-(7.10), we obtain

[TABLE]

7.11 Proof of Lemma 4.7

Proof of (4.6)

Since $\ell({\mathbf{w}})=\ell(\Pi({\mathbf{w}}))$ for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$, we have that $\rho(\Pi({\mathbf{w}}))\geq\rho({\mathbf{w}})$ is equivalent to $\ell_{0}(\Pi({\mathbf{w}}))\geq\ell_{0}({\mathbf{w}})$.

Now, note

[TABLE]

where the last inequality holds by the fact $x-1+e^{-x}\geq 0$ for all $x\in\mathrm{I\!R}$.

Proof of (4.7)

Indeed, $\nabla\rho({\mathbf{w}})=\nabla\ell({\mathbf{w}})+\nabla\ell_{0}({\mathbf{w}})$. It is readily checked that $\nabla\ell({\mathbf{w}})^{\top}{\mathbf{1}}=0$ for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$. We next show that $\nabla\ell_{0}(\Pi({\mathbf{w}}))^{\top}{\mathbf{1}}=0$ for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$.

Note that

[TABLE]

Hence,

[TABLE]

Now, by definition of the mapping $\Pi$ given by (4.4) and (4.5), for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$,

[TABLE]

7.12 Proof of Lemma 7.1

Let $t_{i,j}$ be the number of paired comparisons in the input data with tie outcome for items $i$ and $j$. Note that $t_{i,j}=t_{j,i}$. The log-likelihood function can be written as follows:

[TABLE]

Let $\bar{d}_{i,j}$ be the number of paired comparisons of items $i$ and $j$ such that $i\succeq j$, i.e., $\bar{d}_{i,j}=d_{i,j}+t_{i,j}$. By a straightforward calculus, we can write

[TABLE]

Now, we note when $i\neq j$,

[TABLE]

and

[TABLE]

For any $i\neq j$, it indeed holds

[TABLE]

Hence, when $i\neq j$,

[TABLE]

It follows that $\frac{1}{2}{\mathbf{L}}_{\mathbf{M}}\succeq\nabla^{2}(-\ell({\mathbf{w}}))$ for all ${\mathbf{w}}\in\mathrm{I\!R}^{n}$. Hence,

[TABLE]

This implies that $-\ell$ is a $\frac{1}{2}\lambda_{n}({\mathbf{L}}_{\mathbf{M}})$-smooth function on $\mathrm{I\!R}^{n}$.

On the other hand, we can show that for all ${\mathbf{w}}\in[-\omega,\omega]^{n}$,

[TABLE]

This can be noted as follows. Let $z=\theta e^{w_{j}}/(e^{w_{i}}+\theta e^{w_{j}})$. Note that

[TABLE]

and that $z\in\Omega:=[1/(1+\theta e^{2\omega}),1/(1+\theta e^{-2\omega})]$. The function $z(1-z)$ is convex and thus achieves its minimum value over the interval $\Omega$ at one of its boundary points. It can be readily checked that the minimum is achieved at $z^{*}=1/(1+\theta e^{2\omega})$, which yields $z^{*}(1-z^{*})=c_{\theta,\omega}$.

Hence, when $i\neq j$,

[TABLE]

It follows that $\nabla^{2}(-\ell({\mathbf{w}}))\succeq c_{\theta,\omega}{\mathbf{L}}_{\mathbf{M}}$. From this, we have that for all ${\mathbf{w}}\in[-\omega,\omega]^{n}$ and ${\mathbf{x}}\in\mathcal{X}$,

[TABLE]

where $\mathcal{X}=\{{\mathbf{x}}\in\mathrm{I\!R}^{n}:||{\mathbf{x}}||_{\infty}\leq\omega\hbox{ and }{\mathbf{x}}^{\top}{\mathbf{1}}=0\}$. This implies that $-\ell$ is $c_{\theta,\omega}\lambda_{2}({\mathbf{L}}_{\mathbf{M}})$-strongly convex on $\mathcal{X}$.

7.13 Proof of Lemma 7.5

It can be easily shown that for all ${\mathbf{w}}\in[-\omega,\omega]^{n}$, $S\subseteq N$ such that $|S|\geq 2$, and $u,v\in S$ such that $u\neq v$, we have

[TABLE]

Combining with (7.1), we have

[TABLE]

From this it follows that for all ${\mathbf{x}}\in\mathrm{I\!R}^{n}$ such that ${\mathbf{x}}^{\top}{\mathbf{1}}=0$,

[TABLE]

Similarly, we have

[TABLE]

From this it follows that for all ${\mathbf{x}}$,

[TABLE]

7.14 Derivation of the convergence time bound (6.1)

First note that

[TABLE]

where the inequalities hold elementwise. From this, it follows that ${\mathbf{L}}_{\mathbf{M}}\succeq\underline{m}{\mathbf{L}}_{\mathbf{A}}$ and $\bar{m}{\mathbf{L}}_{\mathbf{A}}\succeq{\mathbf{L}}_{\mathbf{M}}$, where recall ${\mathbf{A}}$ is the adjacency matrix induced by matrix ${\mathbf{M}}$. Now, note

[TABLE]

and

[TABLE]

where $d(n)$ is the maximum degree of a node in graph $G$.

Hence, we have

[TABLE]

By Theorem 3.4 in Merris (1994), for any graph $G$ with adjacency matrix ${\mathbf{A}}$ and diameter $D(n)$, $\lambda_{2}({\mathbf{L}}_{\mathbf{A}})\geq 4/(nD(n))$.

It thus follows that

[TABLE]

which implies the convergence time bound $T=O(rd(n)D(n)n\log(1/\epsilon))$.

7.15 Generalized Bradley-Terry models

In this section, we discuss how the results for Bradley-Terry model of paired comparisons can be extended to other instances of generalized Bradley-Terry models. In particular, we show this for the Rao-Kupper model of paired comparisons with tie outcomes, the Luce choice model and the Plackett-Luce ranking model.

We discuss only the characterization of the strong-convexity and smoothness parameters as the convergence rate bounds for gradient descent and MM algorithms follow similarly as in Section 3.1, from Theorems 3.1 and 3.2, respectively. Similarly, the rate of convergence bounds for accelerated gradient descent and MM algorithms follow readily, similarly to as in Section 4.1, from Theorems, 4.1 and 4.2, respectively.

7.15.1 Model definitions

Bradley-Terry model of paired comparisons

According to the Bradley-Terry model, each paired comparison of items $i$ and $j$ has two possible outcomes: either $i$ wins against $j$ ($i\succ j$) or $j$ wins against $i$ ($j\succ i$). The distribution of the outcomes is given by

[TABLE]

where ${\mathbf{w}}=(w_{1},w_{2},\ldots,w_{n})^{\top}\in\mathrm{I\!R}^{n}$ are model parameters.

Rao-Kupper model of paired comparisons with ties

The Rao-Kupper model is such that each paired comparison of items $i$ and $j$ has three possible outcomes: either $i\succ j$ or $j\succ i$ or $i\equiv j$ (tie). The model is defined by the probability distribution of outcomes that is given by

[TABLE]

and

[TABLE]

where ${\mathbf{w}}=(w_{1},w_{2},\ldots,w_{n})^{\top}\in\mathrm{I\!R}^{n}$ and $\theta\geq 1$ are model parameters.

The larger the value of parameter $\theta$, the more mass is put on the tie outcome. For the value of parameter $\theta=1$, the model corresponds to the Bradley-Terry model for paired comparisons.

Luce choice model

The Luce choice model is a natural generalization of the Bradley-Terry model of paired comparisons to comparison sets of two or more items. For any given comparison set $S\subseteq N=\{1,2,\dots,n\}$ of two or more items, the outcome is a choice of one item $i\in S$ (an event we denote as $i\succeq S$) which occurs with probability

[TABLE]

where ${\mathbf{w}}=(w_{1},w_{2},\ldots,w_{n})^{\top}\in\mathrm{I\!R}_{n}$ are model parameters.

We will use the following definitions and notation. Let $T$ be the set of ordered sequences of two or more items from $N$ such that for each $y=(y_{1},y_{2},\ldots,y_{k})\in T$, $y_{1}$ is an arbitrary item and $y_{2},\ldots,y_{k}$ are sorted in lexicographical order. We can interpret each $y=(y_{1},y_{2},\ldots,y_{k})\in T$ as a choice of item $y_{1}$ from the set of items $\{y_{1},y_{2},\ldots,y_{k}\}$. According to the Luce’s choice model, the probability of outcome $y$ is given by

[TABLE]

We denote with $d_{y}$ the number of observed outcomes $y$ in the input data. For each $y\in T$, let $|y|$ denote the number of items in $y$.

Plackett-Luce ranking model

The Plackett-Luce ranking model is a model of full rankings: for each comparison set of items $S\subseteq N=\{1,2,\ldots,n\}$, the set of possible outcomes contains all possible permutations of items in $S$. The distribution over possible outcomes is defined as follows. Let $T$ be the set of all possible permutations of subsets of two or more items from $N$. Each $y=(y_{1},y_{2},\ldots,y_{k})\in T$ corresponds to a permutation of the set of items $S=\{y_{1},y_{2},\ldots,y_{k}\}$. The probability of outcome $y$ is given by

[TABLE]

where ${\mathbf{w}}=(w_{1},w_{2},\ldots,w_{n})^{\top}\in\mathrm{I\!R}^{n}$ are model parameters.

The model has an intuitive explanation as a sampling of items without replacement proportional to the item weights $e^{w_{i}}$. The Plackett-Luce ranking model corresponds to the Bradley-Terry model of paired comparisons when the comparison sets consist of two items. We denote with $d_{y}$ the number of observed outcomes $y$ in the input data.

In this section, we discuss how the results for Bradley-Terry model of paired comparisons can be extended to other instances of generalized Bradley-Terry models. In particular, we show this for the Rao-Kupper model of paired comparisons with tie outcomes, the Luce choice model and the Plackett-Luce ranking model.

7.15.2 Rao-Kupper model

The probability distribution of outcomes according to the Rao-Kupper model is defined in Section 7.15.1. The log-likelihood function can be written as

[TABLE]

where $\bar{d}_{i,j}$ is the number of observed paired comparisons of items $i$ and $j$ such that either $i$ wins against $j$ or there is a tie outcome, and $t_{i,j}$ is the number of observed paired comparisons of items $i$ and $j$ with tie outcomes.

Lemma 7.1

The negative log-likelihood function for the Rao-Kupper model of paired comparisons with parameter $\theta>1$ is $\gamma$-strongly convex on $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\hbox{ and }{\mathbf{w}}^{\top}{\mathbf{1}}=0\}$ and $\mu$-smooth on $\mathrm{I\!R}^{n}$ with

[TABLE]

where $c_{\theta,\omega}=\theta/(\theta e^{-\omega}+e^{\omega})^{2}$.

Proof of Lemma 7.1 is provided in Section 7.12.

A surrogate minorant function for the log-likelihood function of the Rao-Kupper model is given as follows:

[TABLE]

The MM algorithm is defined by, for $i=1,2,\ldots,n$,

[TABLE]

Lemma 7.2

For all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$, $\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})\geq-\frac{\delta}{2}||{\mathbf{x}}-{\mathbf{y}}||^{2}$ where

[TABLE]

7.15.3 Luce choice model

The probability distribution of outcomes according to the Luce choice model is defined in Section 7.15.1. The log-likelihood function can be written as:

[TABLE]

Lemma 7.3

The negative log-likelihood function for the Luce choice model with comparison sets of size $k\geq 2$ is $\gamma$-strongly convex and $\mu$-smooth on $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\hbox{ and }{\mathbf{w}}^{\top}{\mathbf{1}}=0\}$ with

[TABLE]

where

[TABLE]

and

[TABLE]

Note that for every fixed $\omega>0$, (a) $c_{\omega,k}/d_{\omega,k}$ is decreasing in $k$, (b) $1/e^{8\omega}\leq c_{\omega,k}/d_{\omega,k}\leq 1/e^{2\omega}$, and (c) $1/e^{8\omega}$ is the limit value of $c_{\omega,k}/d_{\omega,k}$ as $k$ goes to infinity.

A minorant surrogate function for the log-likelihood function of the Luce choice model is given by

[TABLE]

The MM algorithm iteration can be written as: for $i=1,2,\ldots,n$,

[TABLE]

where $\sum_{y\in T}d_{y}\mathrm{1\!\!I}_{i=y_{1}}$ is the number of observed comparisons in which item $i$ is the chosen item.

Lemma 7.4

For all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$, $\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})\geq-\frac{\delta}{2}||{\mathbf{x}}-{\mathbf{y}}||^{2}$ where

[TABLE]

7.15.4 Plackett-Luce ranking model

The probability distribution of outcomes according to the Plackett-Luce ranking model is defined in Section 7.15.1. The log-likelihood function can be written as follows:

[TABLE]

Lemma 7.5

The negative log-likelihood function for the Plackett-Luce ranking model with comparison sets of size $k\geq 2$ is $\gamma$-strongly convex and $\mu$-smooth on $\mathcal{W}_{\omega}=\{{\mathbf{w}}\in\mathrm{I\!R}^{n}:||{\mathbf{w}}||_{\infty}\leq\omega\hbox{ and }{\mathbf{w}}^{\top}{\mathbf{1}}=0\}$ with

[TABLE]

where

[TABLE]

Proof of Lemma 7.5 is provided in Section 7.13.

Note that for fixed values of $\omega$ and $k$, Lemma 7.5 implies the convergence time $\log(d({\mathbf{M}})/a({\mathbf{M}}))$. Note, however, that for fixed $\omega>0$, $\tilde{c}_{\omega,k}/\tilde{d}_{\omega,k}$ decreases to [math] with $k$ and is of the order $1/k^{2}$. This is because in the derivation of parameters $\tilde{c}_{\omega,k}$ and $\tilde{d}_{\omega,k}$ we use (conservative) deterministic bounds. Following Hajek et al. (2014), one can derive bounds for $\gamma$ and $\mu$ that hold with high probability, which are such that $\tilde{c}_{\omega,k}$ and $\tilde{d}_{\omega,k}$ scale with $k$ in the same way.

The log-likelihood function of the Plackett-Luce ranking model admits the following minorization function:

[TABLE]

The MM algorithm is given by: for $i=1,2,\ldots,n$,

[TABLE]

where $S_{a,b}(y)=\{y_{a},y_{a+1},\ldots,y_{b}\}$.

Lemma 7.6

For all ${\mathbf{x}},{\mathbf{y}}\in[-\omega,\omega]^{n}$, $\underline{\ell}({\mathbf{x}};{\mathbf{y}})-\ell({\mathbf{x}})\geq-\frac{\delta}{2}||{\mathbf{x}}-{\mathbf{y}}||^{2}$ where

[TABLE]