Spectral Method and Regularized MLE Are Both
Optimal for Top-K Ranking
Yuxin Chen∗

Jianqing Fan†

arXiv:1707.09971v2 [stat.ML] 28 Sep 2017

August 2017;

Cong Ma†

Kaizheng Wang†

Revised September 2017

Abstract
This paper is concerned with the problem of top-K ranking from pairwise comparisons. Given a
collection of n items and a few pairwise binary comparisons across them, one wishes to identify the set
of K items that receive the highest ranks. To tackle this problem, we adopt the logistic parametric
model—the Bradley-Terry-Luce model, where each item is assigned a latent preference score, and where
the outcome of each pairwise comparison depends solely on the relative scores of the two items involved.
Recent works have made significant progress towards characterizing the performance (e.g. the mean
square error for estimating the scores) of several classical methods, including the spectral method and
the maximum likelihood estimator (MLE). However, where they stand regarding top-K ranking remains
unsettled.
We demonstrate that under a random sampling model, the spectral method alone, or the regularized
MLE alone, is minimax optimal in terms of the sample complexity—the number of paired comparisons
needed to ensure exact top-K identification. This is accomplished via optimal control of the entrywise
error of the score estimates. We complement our theoretical studies by numerical experiments, confirming that both methods yield low entrywise errors for estimating the underlying scores. Our theory is
established based on a novel leave-one-out trick, which proves effective for analyzing both iterative and
non-iterative optimization procedures. Along the way, we derive an elementary eigenvector perturbation
bound for probability transition matrices, which parallels the Davis-Kahan sin Θ theorem for symmetric
matrices. This further allows us to close the gap between the `2 error upper bound for the spectral
method and the minimax lower limit.
Keywords. top-K ranking, pairwise comparisons, spectral method, regularized MLE, perturbation
analysis, Erdős–Rényi random graph, leave-one-out analysis, reversible Markov chain.

1

Introduction

Imagine we have a large collection of n items, and we are given partially revealed comparisons between pairs
of items. These paired comparisons are collected in a non-adaptive fashion, and could be highly noisy and
incomplete. The aim is to aggregate these partial preferences so as to identify the K items that receive the
highest ranks. This problem, which is called top-K rank aggregation, finds applications in numerous contexts,
including web search (Dwork et al., 2001), recommendation systems (Baltrunas et al., 2010; Park et al., 2015),
sports competition (Masse, 1997), and so on. The challenge is both statistical and computational: how can
one achieve reliable top-K ranking from a minimally possible number of pairwise comparisons, while retaining
computational efficiency?

1.1

Popular approaches

To address the aforementioned challenge, many prior approaches have been put forward based on certain
statistical models. Arguably one of the most widely used parametric models is the Bradley-Terry-Luce
Author names are sorted alphabetically.
of Electrical Engineering, Princeton University, Princeton, NJ 08544, USA; Email: yuxin.chen@princeton.edu.
† Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ 08544, USA; Email:
{jqfan, congm, kaizheng}@princeton.edu.
∗ Department

1

(BTL) model (Bradley and Terry, 1952; Luce, 1959; Agresti, 2013), which assigns a latent preference score
wi∗ (1 ≤ i ≤ n) to each of the n items. The BTL model posits that: the chance of each item winning a paired
comparison is determined by the relative scores of the two items involved, or more precisely,
P {item j is preferred over item i} =

wj∗
wi∗ + wj∗

(1)

in each comparison of item i against item j. The items are repeatedly compared in pairs according to this
parametric model. The task then boils down to identifying the K items with the highest preference scores,
given these pairwise comparisons.
Among the ranking algorithms tailored to the BTL model, there are two procedures that have received
particular attention, both of which rank the items in accordance with appropriate estimates of the latent
preference scores:
1. The spectral method. By connecting the winning probability in (1) with the transition probability of
a reversible Markov chain, the spectral method attempts recovery of {wi∗ }1≤i≤n via the leading left
eigenvector of a sample transition matrix. This spectral ranking procedure, which is also known as
Rank Centrality (Negahban et al., 2017a), bears some similarity to the PageRank algorithm (Page
et al., 1999).
2. The maximum likelihood estimator (MLE). This approach proceeds by finding the score assignment
that maximizes the likelihood function (Ford, 1957). When parameterized appropriately, solving the
MLE becomes a convex program, and hence is computationally feasible. There are also important
variants of the MLE that enforce additional regularization; see Negahban et al. (2017a).
Details are postponed to Section 2.2. In addition to their remarkable practical applicability, these two
ranking paradigms are appealing in theory as well. For instance, both of them provably achieve near-optimal
mean square errors when estimating the latent preference scores (Negahban et al., 2017a; Hajek et al., 2014).
Nevertheless, the `2 error for estimating the latent scores merely serves as a “meta-metric” for the ranking
task, which does not necessarily reveal the accuracy of top-K identification. Given that the `2 loss only
reflects the estimation error in some average sense, it is certainly possible that an algorithm obtains minimal
`2 estimation loss but incurs (relatively) large errors when estimating the scores of the highest ranked items.
Interestingly, a recent work Chen and Suh (2015) demonstrates that: a careful combination of the spectral
method and the coordinate-wise MLE is optimal for top-K ranking. This leaves open the following natural
questions: where does the spectral alone, or the MLE alone, stand in top-K ranking? Are they capable of
attaining exact top-K recovery from minimal samples? These questions are poorly understood and form the
primary objectives of our study.
As we will elaborate later, the spectral method part of the preceding question has recently been explored
in an interesting paper Jang et al. (2016), for a regime where a relatively large fraction of item pairs have been
compared. However, it remains unclear how well these two algorithms can perform in a much broader—and
often more challenging—regime.

1.2

Main contributions

The central focal point of the current paper is to assess the accuracy of both the spectral method and the
regularized MLE in top-K identification. Assuming that the pairs of items being compared are randomly
selected, our paper delivers a somewhat surprising message:
Both the spectral method and the regularized MLE achieve perfect identification of top-K ranked
items under optimal sample complexity (up to some constant factor)!
It is worth emphasizing that these two algorithms work even under the sparsest possible regime, a
scenario where only an exceedingly small fraction of pairs of items have been compared. This calls for
precise control of the entrywise error—as opposed to the `2 loss—for estimating the scores. To this end,
our theory is established upon a novel leave-one-out argument, which might shed light on how to analyze
the entrywise error for more general optimization problems. As a byproduct of the analysis, we derive an
2

elementary eigenvector perturbation bound for probability transition matrices, which parallels Davis-Kahan’s
sin Θ theorem for symmetric matrices. This simple perturbation bound leads to an improved `2 error bound
for the spectral method, which matches the minimax lower limit derived in the prior literature.
Finally, we remark that both methods can be performed in time nearly proportional to reading all
samples. In fact, the near-optimal computational cost of the spectral method has already been demonstrated
in Negahban et al. (2017a). When it comes to the regularized MLE, we will show that this method—which
only exhibits a fairly weak level of strong convexity—can be computed by the gradient descent in nearly linear
time. In summary, both methods provably achieve optimal efficiency from both statistical and computational
perspectives.

1.3

Notation

Before continuing, we introduce a few notations that will be useful throughout. To begin with, vectors and
matrices are in boldface whereas scalars are not.
For any strictly positive probability vector π ∈ Rn , we define the inner product space indexed by π as a
real vector space in Rn endowed with the inner product (Brémaud, 1999)
Xn
hx, yiπ =
πi xi yi .
i=1

The corresponding vector norm and the induced matrix norm are defined respectively as
q
and
kAkπ = sup kx> Akπ .
kxkπ = hx, xiπ
kxkπ =1

If 0 < πmin ≤ πi ≤ πmax for all 1 ≤ i ≤ n, then one has the following elementary inequalities
r
r
√
√
πmin
πmax
kAk2 ≤ kAkπ ≤
kAk2 .
πmin kxk2 ≤ kxkπ ≤ πmax kxk2
and
πmax
πmin
Additionally, the standard notation f (n) = O (g(n)) or f (n) . g(n) means that there exists a constant
c > 0 such that |f (n)| ≤ c|g(n)|, f (n) = Ω (g(n)) or f (n) & g(n) means that there exists a constant c > 0
such that |f (n)| ≥ c |g(n)|, f (n) = Θ (g(n)) or f (n)  g(n) means that there exist constants c1 , c2 > 0 such
(n)
that c1 |g(n)| ≤ |f (n)| ≤ c2 |g(n)|, and f (n) = o(g(n)) means that limn→∞ fg(n)
= 0.
P
>
Given a graph G with vertex set {1, 2, . . . , n} and edge set E, we denote by LG = (i,j)∈E,i>j (ei − ej ) (ei − ej )
the (unnormalized) Laplacian matrix (Chung, 1997) associated with it, where {ei }1≤i≤n are the standard
basis vectors in Rn . For a matrix A ∈ Rn×n with n real eigenvalues, we let λ1 (A) ≥ λ2 (A) ≥ · · · ≥ λn (A)
be the eigenvalues sorted in descending order.

2
2.1

Statistical models and main results
Problem setup

We begin with a formal introduction of the Bradley-Terry-Luce parametric model for binary comparisons.
Preference scores. As introduced earlier, we assume the existence of a positive latent score vector
w∗ = [w1∗ , · · · , wn∗ ]>

(2)

that comprises the underlying preference score {wi∗ }1≤i≤n assigned to each of the n items. Alternatively,
it is sometimes more convenient to reparameterize the score vector by
θ ∗ = [θ1∗ , · · · , θn∗ ]> ,

θi∗ = log wi∗ .

(3)

These scores are assumed to fall within a fixed dynamic range, namely,
wi∗ ∈ [wmin , wmax ],

or
3

θi∗ ∈ [θmin , θmax ]

(4)

for all 1 ≤ i ≤ n, where wmin > 0, wmax > 0, θmin = log wmin , and θmax = log wmax are some absolute
constants independent of n. Without loss of generality, it is assumed that
wmax ≥ w1∗ ≥ w2∗ ≥ . . . ≥ wn∗ ≥ wmin ,

(5)

meaning that items 1 through K are the desired top-K ranked items. For notational simplicity, we
define
wmax
,
(6)
κ :=
wmin
which is assumed to be a constant independent of n.
Comparison graph. Let G = (V, E) stand for a comparison graph, where the vertex set V = {1, 2, . . . , n}
represents the n items of interest. The items i and j are compared if and only if (i, j) falls within the
edge set E. Unless otherwise noted, we assume that G is drawn from the Erdős–Rényi random graph
Gn,p (Durrett, 2007), such that an edge between any pair of vertices is present independently with some
probability p. In words, p captures the fraction of item pairs being compared. In order to ensure that the
comparison graph is connected (otherwise there is no basis to rank items belonging to two disconnected
components of G), we assume throughout that
p≥

c0 log n
n

(7)

for some constant c0 > 1.
Pairwise comparisons. For each (i, j) ∈ E, we obtain L independent paired comparisons between items i
(l)
and j. Let yi,j be the outcome of the `-th comparison, which is independently drawn as

∗
1, with probability wj∗ = eθj ∗ ,
(l) ind.
θ
θ∗
wi∗ +wj∗
e i +e j
yi,j =
(8)
0, else.
(l)

(l)

By convention, we set yi,j = 1 − yj,i for all (i, j) ∈ E throughout the paper. This is also known as the
logistic pairwise comparison model, due to its strong resemblance to logistic regression (Agresti, 2013).
It is self-evident that the sufficient statistics under this model are given by
L

y := {yi,j | (i, j) ∈ E} ,

yi,j :=

1 X (l)
yi,j .
L

(9)

l=1

To simplify the notation, we shall also take
∗

∗
yi,j

wj∗
eθ j
:= ∗
=
∗ .
∗
wi + wj∗
eθi + eθj

(10)

Goal. The goal is to identify the set of top-K ranked items—that is, the set of K items that exhibit the
largest preference scores—from the pairwise comparisons y.

2.2
2.2.1

Algorithms
The spectral method: Rank Centrality

The spectral ranking algorithm, or Rank Centrality (Negahban et al., 2017a), is motivated by the connection
between the pairwise comparisons and a random walk over a directed graph. The algorithm starts by
converting the pairwise comparison data y into a transition matrix P = [Pi,j ]1≤i,j≤n in such a way that

1

if (i, j) ∈ E,
 d yi,j , P
1
Pi,j = 1 − d k:(i,k)∈E yi,k , if i = j,
(11)


0,
otherwise,
4

Algorithm 1 Spectral method (Rank Centrality).
Input the comparison graph G, the sufficient statistics y, and the normalization factor d.
Define the probability transition matrix P = [Pi,j ]1≤i,j≤n as in (11).
Compute the leading left eigenvector π of P .
Output the K items that correspond to the K largest entries of π.

for some given normalization factor d, and then proceeds by computing the stationary distribution π ∈ Rn
of the Markov chain induced by P . As we shall see later, the parameter d is taken to be on the same order
of the maximum vertex degree of G while ensuring the non-negativity of P . As asserted by Negahban et al.
(2017a), π is a faithful estimate of w∗ up to some global scaling. The algorithm is summarized in Algorithm
1.
To develop some intuition regarding why this spectral algorithm gives a reasonable estimate of w∗ , it is
∗
perhaps more convenient to look at the population transition matrix P ∗ = [Pi,j
]1≤i,j≤n :

wj∗

if (i, j) ∈ E,
 d1 wi∗ +wj∗ ,

∗
P
∗
w
Pi,j = 1 − 1 k:(i,k)∈E ∗ k ∗ ,
if i = j,
d
wi +wk



0,
otherwise,
which coincides with P by taking L → ∞. It can be seen that the normalized score vector
π ∗ := Pn

1

∗
i=1 wi

>

[w1∗ , w2∗ , . . . , wn∗ ]

(12)

is the stationary distribution of the Markov chain induced by P ∗ , since P ∗ and π ∗ are in detailed balance
(Brémaud, 1999), namely,
∗
∗
πi∗ Pi,j
= πj∗ Pj,i
,
∀ (i, j) .
(13)
As a result, one expects the stationary distribution of the sample version P to form a good estimate of w∗ ,
provided that the sample size is sufficiently large.
2.2.2

The regularized MLE

Under the BTL model, the negative log-likelihood function conditioned on the graph is given by (up to some
global scaling)

X 
eθ j
eθi
+ (1 − yj,i ) log θi
L (θ; y) : = −
yj,i log θi
e + eθ j
e + eθ j
(i,j)∈E,i>j
X 

=
−yj,i (θi − θj ) + log 1 + eθi −θj .
(14)
(i,j)∈E,i>j

The regularized MLE (Negahban et al., 2017a, Section 3.3) then amounts to solving the following convex
program
minimizeθ∈Rn

1
Lλ (θ; y) := L (θ; y) + λkθk22 ,
2

(15)

q
n
for a regularization parameter λ > 0. As will be discussed later, we shall adopt the choice λ  np log
L
throughout this paper. For the sake of brevity, we let θ represent the resulting penalized maximum likelihood
estimate whenever it is clear from the context.

2.3

Main results

The most challenging part of top-K ranking is to distinguish the K-th and the (K + 1)-th items. In
fact, the score difference of these two items captures the distance between the item sets {1, . . . , K} and
5

{K + 1, . . . , n}. Unless their latent scores are sufficiently separated, the finite-sample nature of the model
would make it infeasible to distinguish these two critical items. With this consideration in mind, we define
the following separation measure
∗
w∗ − wK+1
.
(16)
∆K := K
wmax
This metric turns out to play a crucial role in determining the minimal sample complexity for perfect top-K
identification.
The main finding of this paper concerns the optimality of both the spectral method and the regularized
MLE. Recall that under the BTL model introduced before, the total number N of samples we collect
concentrates sharply around its mean, namely,
1
N = (1 + o(1))E [N ] = (1 + o(1)) n2 pL
2

(17)

occurs with high probability. Our main result is stated in terms of the sample complexity required for exact
top-K identification.
Theorem 1. Consider the pairwise comparison model specified in Section 2.1. Suppose that p >
that
c1 n log n
n2 pL
≥
2
∆2K

c0 log n
n

and
(18)

for some sufficiently large constants c0 , c1 > 0. With probability exceeding 1 − O(n−5 ), the set of top-K
ranked items can be recovered exactly by the spectral method given in Algorithmq1, and by the regularized

MLE given in (15). Here, we take d = cd np in the spectral method and λ = cλ
MLE, where cd ≥ 2 and cλ > 0 are some absolute constants.

np log n
L

in the regularized

Theorem 1 asserts that both the spectral method and the regularized MLE achieve a sample complexity
log n
on the order of n ∆
. Encouragingly, this sample complexity coincides with the minimax limit identified in
2
K
(Chen and Suh, 2015, Theorem 2) up to a constant factor:
Theorem 2 (Chen and Suh (2015)). Fix  ∈ (0, 1/2), and suppose that
n2 pL ≤ 2c2

(1 − )n log n − 2
,
∆2K

(19)

4
4
). Then for any ranking procedure ψ, one can find a score vector w∗ with separation
where c2 = wmin
/(4wmax
∆K such that ψ fails to retrieve the top-K items with probability at least .

Moving beyond the sample complexity issue, we remark that the spectral method can be carried out
by means of the power method, while the regularized MLE can be computed via the standard gradient
descent (as we will elaborate in Section 4.2.1). Both procedures are computationally feasible and can be
accomplished within nearly linear time—in time proportional to reading all data. Consequently, these two
algorithms are optimal in terms of both computational complexity and sample complexity. We emphasize
that it is nontrivial to establish the nearly linear-time computational complexity of the regularized MLE,
given that the log-likelihood function (or the logistic link function) is not strongly convex. We will elaborate
on this point in Section 4.2.1.
We are now positioned to compare our results with Jang et al. (2016), which also investigates the accuracy
of the spectral method for top-K ranking. Specifically, (Jang et al., 2016, Theorem 3) establishes the
optimality of the spectral method for the relatively dense regime where
r
log n
p&
.
n
In this regime, however, the total sample size necessarily exceeds
n2 pL
n2 p p 3
≥
& n log n,
2
2
6

(20)

which rules out the possibility of achieving minimal sample complexity if ∆K is sufficiently large. For
instance, consider the case where ∆K  1, then the optimal sample size—as revealed by Theorem 1 or
(Chen and Suh, 2015, Theorem 1)—is on the order of
n log n
 n log n,
∆2K
p
which is a factor of n/ log n lower than the bound in (20). In contrast, our results hold all the way down to
the sparsest possible regime where p  logn n , confirming the optimality of the spectral method even for the
most challenging scenario. Furthermore, we establish that the regularized MLE shares the same optimality
guarantee as the spectral method, which was previously out of reach.

2.4

Optimal control of entrywise estimation errors

In order to establish the ranking accuracy as asserted by Theorem 1, the key is to obtain precise control of
the `∞ loss of the score estimates. Our results are as follows.
n
Theorem 3 (Entrywise error of the spectral method). Suppose that p > c0 log
for some sufficiently
n
large constant c0 > 0. Choose d = cd np for some constant cd ≥ 2 in Algorithm 1. Then the spectral estimate
π satisfies
s
kπ − π ∗ k∞
log n
(21)
.
kπ ∗ k∞
npL

with probability exceeding 1 − O(n−5 ), where π ∗ is the normalized score vector as defined in (12).
n
Theorem 4 (Entrywise error of the regularized MLE). Suppose that p > c0 log
for some sufficiently
n
q
np log n
large constant c0 > 0, and that the regularization parameter is λ = cλ
for some absolute constant
L
cλ > 0. Then the regularized MLE θ satisfies
s
∗
∗
keθ − eθ −θ 1 k∞
log n
.
(22)
∗ −θ ∗ 1
θ
npL
ke
k∞
∗

with probability exceeding 1 − O(n−5 ), where θ :=

1 > ∗
n1 θ

and eθ := [eθ1 , · · · , eθn ]> .

Theorems 3-4 indicate that if the number of comparisons associated with each item—which concentrates
around npL—exceeds log n to some sufficiently large extent, then both methods are able to achieve a small
`∞ error when estimating the scores.
Recall that the `2 estimation error of the spectral method has been characterized by Negahban et al.
(2017a) (or Theorem 6 of this paper with an improved upper bound), which obeys
s
kπ − π ∗ k2
log n
.
(23)
kπ ∗ k2
npL
with high probability. Similar theoretical guarantees have been derived for another variant of the MLE
(the constrained version) under a uniform sampling model as well (Negahban et al., 2017a). In comparison,
our results indicate that the estimation errors for both algorithms are almost evenly spread out across all
coordinates rather than being localized or clustered. Notably, the pointwise errors revealed by Theorems 3-4
immediately lead to exact top-K identification as claimed by Theorem 1.
Proof of Theorem 1. In what follows, we prove the theorem for the spectral method part. The regularized
MLE part follows from an almost identical argument and hence is omitted.
Since the spectral algorithm ranks the items in accordance with the score estimate π, it suffices to
demonstrate that
πi − πj > 0,
∀1 ≤ i ≤ K, K + 1 ≤ j ≤ n.

7

To this end, we first apply the triangle inequality to get
πj − πj∗
πi∗ − πj∗
πi − πj
|πi − πi∗ |
2kπ − π ∗ k∞
≥
−
−
≥ ∆K −
.
∗
∗
∗
∗
kπ k∞
kπ k∞
kπ k∞
kπ k∞
kπ ∗ k∞

(24)

In addition, it follows from Theorem 3 as well as our sample complexity assumption that
s
kπ − π ∗ k∞
n log n
log n
.
and
n2 pL &
.
∗
kπ k∞
npL
∆2K
These conditions taken collectively imply that
kπ − π ∗ k∞
1
< ∆K
∗
kπ k∞
2
as long as

2.5

npL∆2K
log n

is sufficiently large. Substitution into (24) reveals that πi − πj > 0, as claimed.

Heuristic arguments

We pause here to develop some heuristic explanation as to why the score estimation errors are expected to
be spread out across all entries. For simplicity, we focus on the case where p = 1 and L is sufficiently large,
so that y and P sharply concentrate around y ∗ and P ∗ , respectively. For simplicity of presentation, we
.
introduce the notation  to represent asymptotic equivalence up to logarithmic factors.
We begin with the spectral algorithm. Since π and π ∗ are respectively the invariant distributions of the
Markov chains induced by P and P ∗ , we can decompose
>

>

(π − π ∗ ) = π > P − π ∗> P ∗ = (π − π ∗ ) P + π ∗> (P − P ∗ ).
|
{z
}

(25)

:=ξ

∗
∗
∗
max
When p = 1 and w
wmin  1, the entries of π (resp. the off-diagonal entries of P and P − P ) are all of the
same order and, as a result, the energy of the uncertainty term ξ is spread out. In fact, we will demonstrate
in Section 3.2 that
s
kξk∞
log n . kπ − π ∗ k2

.
,
(26)
kπ ∗ k∞
npL
kπ ∗ k2

which coincides with the optimal rate. Further, if we look at each entry of (25), then for all 1 ≤ m ≤ n,


>
∗
πm − πm
= (π − π ∗ ) P m + ξm ,

π1 − π1∗


..
∗
= (πm − πm
) Pm,m + [P1,m , · · · , Pm−1,m , 0, Pm+1,m , · · · , Pn,m ] 
 + ξm .
.
∗
πn − πn
|
{z
} |
{z
}


contraction

(27)

error averaging

By construction of the transition matrix, one can easily verify that Pm,m is bounded away from 1 and
∗
Pj,m  1/n for all j 6= m. As a consequence, the identity π > P = π > allows one to treat each πm − πm
as a mixture of three effects: (i) the first term of (27) behaves as a entrywise contraction of the error; (ii)
the second term of (27) is a (nearly uniformly weighted) average of the errors over all coordinates, which
can essentially be treated as a smoothing operator applied to the error components; and (iii) the uncertainty
term ξm . Rearranging terms in (27), we are left with
∗
(1 − Pm,m ) |πm − πm
|.

8

1 Xn
|πi − πi∗ | + ξm ,
i=1
n

(28)

or equivalently,
kπ − π ∗ k∞ .

1 Xn
|πi − πi∗ | + kξk∞ .
i=1
n

(29)
Pn

There are two possibilities compatible with this bound: (1) kπ − π ∗ k∞ . n1 i=1 |πi − πi∗ |, and (2) kπ − π ∗ k∞ .
∗
k2
∗
kξk∞ . kπ−π
kπ ∗ k2 kπ k∞ by (26). In either case, the errors are fairly delocalized, revealing that
kπ − π ∗ k∞ . kπ − π ∗ k2

.
kπ ∗ k∞
kπ ∗ k2
We now move on to the regularized MLE, which follows a very similar argument. By the optimality
condition that ∇Lλ (θ) = 0, one can derive (for some η to be specified later)
θ − θ ∗ = θ − η∇Lλ (θ) − θ ∗
= θ − η∇Lλ (θ) − (θ ∗ − η∇Lλ (θ ∗ )) − η∇Lλ (θ ∗ )
| {z }
:=ζ

≈ I − η∇2 Lλ θ ∗



(θ − θ ∗ ) − ζ.

Write ∇2 Lλ (θ ∗ ) = D − A, where D and A denote respectively the diagonal and off-diagonal parts. Under
our assumptions, one can check that Dm,m  n for all 1 ≤ m ≤ n and Aj,m  1 for any j 6= m. With these
notations in place, one can write the entrywise error as follows
X

∗
∗
θm − θm
= (1 − ηDm,m ) (θm − θm
)+
ηAj,m θj − θj∗ − ζm .
j:j6=m

By choosing η = c2 /n for some sufficiently small constant c2 > 0 , we get 1 − ηDm,m < 1 and ηAj,m  1/n.
Therefore, the right-hand side of the above relation also comprises a contraction term as well as an error
smoothing term, similar to (27). Carrying out the same argument as for the spectral method, we see that
the estimation errors of the regularized MLE are expected to be spread out.

2.6

Numerical experiments

It is worth noting that extensive numerical experiments have already been conducted in Negahban et al.
(2017a) to justify the practicability of both the spectral method and the regularized MLE. This section
provides some additional simulations to complement their experimental results and our theory. Throughout
the experiments, we set the number of items n to be 200, while the number of repeated comparisons L and
the edge probability p can vary with the experiments. For the tuning parameters,
q we choose d = 2dmax in

n
the spectral method where dmax is the maximum degree of the graph and λ = 2 np log
in the regularized
L
MLE, which are consistent with the configurations considered in the main theorems. All of the results are
averaged over 100 Monte Carlo simulations.
We first investigate the relative `∞ error of the spectral method and the regularized MLE when estimating
the preference scores. To this end, we generate the latent scores wi∗ (1 ≤ i ≤ n) independently and uniformly
at random over the interval [0.5, 1]. Figure 1(a) (resp. Figure 1(b)) displays the entrywise error in the score
estimation as the number of repeated comparisons L (resp. the edge probability p) varies. As can be seen
from the plots, the `∞ error of both methods gets smaller as p and L increase, confirming our results in
Theorems 3-4. Next, we show in Figure 2 the relative `∞ error and top-K ranking accuracy while fixing the
total number of samples (i.e. n2 pL). It can be seen that the performance almost does not change if the total
sample complexity n2 pL remains the same. In addition, Figure 3 illustrates the relative `∞ error and the
relative `2 error in score estimation for both methods. As we can see, the relative `∞ errors are not much
larger than the relative `2 errors (recall that n = 200), thus providing empirical evidence that the errors in
the score estimates are spread out across all entries.
Furthermore, we examine the top-K ranking accuracy of both methods. Here, we fix p = 0.25 and L = 20,
set K = 10, and let wi∗ = 1 for all 1 ≤ i ≤ K and wj∗ = 1 − ∆ for all K + 1 ≤ j ≤ n. By construction,
the score separation satisfies ∆K = ∆. Figure 4 illustrates the accuracy of both methods in identifying the
top-K ranked items. The performance of both methods improves when the score separation becomes larger,
which matches our theory in Theorem 1.

9

0.35

relative `1 estimation error

relative `1 estimation error

0.4
Spectral Method p = 0:25
Spectral Method p = 0:5
Spectral Method p = 0:75
Regularized MLE p = 0:25
Regularized MLE p = 0:5
Regularized MLE p = 0:75

0.3

Spectral Method L = 10
Spectral Method L = 30
Spectral Method L = 50
Regularized MLE L = 10
Regularized MLE L = 30
Regularized MLE L = 50

0.3

0.25

0.2

0.2

0.15

0.1

0.1

0.05
0

0

5

10

15

20

25

30

0

0.2

0.4

0.6

0.8

1

p: edge probability

L: number of repeated comparisons

(a)

(b)

Figure 1: Empirical performance of both the spectral method and the regularized MLE: (a) `∞ error vs. L,
and (b) `∞ error vs. p.

1

relative `1 estimation error

0.25
Spectral Method pL = 20
Spectral Method pL = 40
Spectral Method pL = 60
Regularized MLE pL = 20
Regularized MLE pL = 40
Regularized MLE pL = 60

0.2

0.8

0.15

0.6

0.1

0.4

0.05

0.2

0

0
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

p: edge probability

(a) relative `∞ error

(b) top-K ranking accuracy

Figure 2: Empirical performance of both the spectral method and the regularized MLE when fixing the
number of samples n2 pL: (a) relative `∞ error, and (b) top-K ranking accuracy.

0.25

`1 error, p = 0:25
`1 error, p = 0:5
`1 error, p = 0:75
`2 error, p = 0:25
`2 error, p = 0:5
`2 error, p = 0:75

0.3

relative estimation error

relative estimation error

0.4

0.2

0.1

0

5

10

15

20

25

0.2

0.15

0.1

0.05

30

L: number of repeated comparisons

`1 error, p = 0:25
`1 error, p = 0:5
`1 error, p = 0:75
`2 error, p = 0:25
`2 error, p = 0:5
`2 error, p = 0:75

5

10

15

20

25

30

L: number of repeated comparisons

(a) spectral method

(b) regularized MLE

Figure 3: Comparisons between the relative `∞ error and the relative `2 error for (a) the spectral method,
and (b) the regularized MLE.

10

top-K ranking accuracy

1
Spectral Method
Regularized MLE

0.8
0.6
0.4
0.2
0
0.1

0.2

0.3

0.4

0.5

"K : score separation

Figure 4: The top-K ranking accuracy of both the spectral method and the regularized MLE.

2.7

Other related works

The problem of ranking based on partial preferences has received much attention during the past decade.
Two types of observation models have been considered: (1) the cardinal-based model, where users provide
explicit numerical ratings of the items being compared, (2) the ordinal-based model, where users are asked to
make comparative measurements. See Ammar and Shah (2011); Shah et al. (2014) for detailed experimental
and theoretical comparisons between these two distinct models.
In terms of the ordinal-based model—and in particular, ranking from pairwise comparisons—both parametric and nonparametric models have been extensively studied. For example, Hunter (2004) examined
variants of the parametric BTL model, and established the convergence properties of the minorizationmaximization algorithm for solving the corresponding MLE problem. The BTL model was also investigated
in Borkar et al. (2016), where the authors proposed a randomized Kaczmarz algorithm to infer the latent
scores with provable `2 error guarantees. Moreover, the BTL model falls under the category of low-rank
parametric models, since the preference matrix is generated by passing a rank-2 matrix through the logistic
link function (Rajkumar and Agarwal, 2016). Additionally, the work Jiang et al. (2011); Xu et al. (2012)
proposed a least-squares type method, called HodgeRank, to estimate the full ranking, which generalizes the
simple Borda count algorithm (Borda, 1784; Ammar and Shah, 2011; Rajkumar and Agarwal, 2014). For
many of these algorithms, the sample complexities needed for perfect total ranking have been determined
by Rajkumar and Agarwal (2014), although the top-K ranking accuracy was not considered there.
Going beyond the parametric models, a recent line of works Shah et al. (2017); Shah and Wainwright
(2015); Chen et al. (2017); Pananjady et al. (2017) considered the nonparametric stochastically transitive
model, where the only model assumption is that the comparison probability matrix follows certain transitivity
rules. This type of models subsumes the BTL model as a special case. For instance, Shah and Wainwright
(2015) suggested an exceedingly simple counting-based algorithm which can reliably recover the top-K ranked
items for a wide range of models. However, the sampling paradigm considered therein is quite different from
ours in the sparse regime; for instance, their model does not come close to the setting where p is small but
L is large, which is the most challenging regime of the model adopted in our paper and Negahban et al.
(2017a); Chen and Suh (2015); Jang et al. (2016).
All of the aforementioned papers concentrate on the case where there is a single ground-truth ordering.
It would also be interesting to investigate the scenarios where different users might have different preference
scores (Negahban et al., 2017b; Lu and Negahban, 2014; Wu et al., 2015; Oh and Shah, 2014). To this
end, Negahban et al. (2017b); Lu and Negahban (2014) imposed the low-rank structure on the underlying
preference matrix and adopted the nuclear-norm relaxation approach to recover the users’ preferences. The
paper Wu et al. (2015) proposed a two-step procedure, where a clustering step was first applied to label
different users, followed by an estimation step built upon the previous clustering outcome. Additionally,
several papers explored the ranking problem for the more general Plackett-Luce model (Hajek et al., 2014;
Soufiani et al., 2013, 2014), in the presence of adaptive sampling (Braverman and Mossel, 2008; Jamieson
and Nowak, 2011; Ailon, 2012; Busa-Fekete et al., 2013; Mohajer and Suh, 2016; Heckel et al., 2016; Agarwal
et al., 2017), for the crowdsourcing scenario (Chen et al., 2013; Osting et al., 2016), and in the adversarial
11

setting (Suh et al., 2017). These are beyond the scope of the present paper.
When it comes to the technical tools, it is worth noting that the leave-one-out idea has been invoked to
analyze random designs for other high-dimensional problems as well, e.g. robust M-estimators (El Karoui
et al., 2013; El Karoui, 2017), confidence intervals for Lasso (Javanmard and Montanari, 2015), likelihood
ratio test (Sur et al., 2017). In particular, Zhong and Boumal (2017) and Abbe et al. (2017) use this technique
to precisely characterize entrywise behavior of eigenvectors of a large class of symmetric random matrices,
which improves upon the prior analysis on `∞ eigenvector perturbation (Fan et al., 2016). As a consequence,
they are able to show the sharpness of spectral methods—a type of methods based on the eigen-structure of
a certain data matrix—in many popular models. Our introduction of leave-one-out auxiliary quantities are
motivated by these papers.
Finally, we remark that the family of spectral methods has been successfully applied in numerous applications, e.g. matrix completion (Keshavan et al., 2010; Jain et al., 2013), phase retrieval (Candès et al.,
2015; Chen and Candès, 2017), graph clustering (Rohe et al., 2011; Chaudhuri et al., 2012), synchronization
and joint alignment (Singer, 2011; Chen and Candes, 2016). All of them are designed based on the eigenvectors of some symmetric matrix, or the singular vectors if the matrix of interest is asymmetric. Our paper
contributes to this growing literature by establishing a sharp eigenvector perturbation analysis framework
for an important class of asymmetric matrices—the probability transition matrices.

3

Analysis for the spectral method

This section is devoted to proving Theorem 3, which characterizes the pointwise error of the spectral estimate.

3.1

Preliminaries

Here, we gather some preliminary facts about reversible Markov chains as well as the Erdős–Rényi random
graph.
The first important result concerns the eigenvector perturbation for probability transition matrices, which
can be treated as the analogue of the celebrated Davis-Kahan sin Θ theorem (Davis and Kahan, 1970) for
asymmetric matrices. Due to its potential importance for other problems, we promote it to a theorem as
follows.
Theorem 5 (Eigenvector perturbation). Suppose that P , P̂ , and P ∗ are probability transition matrices
with stationary distributions π, π̂, π ∗ , respectively. Also, assume that P ∗ represents a reversible Markov
chain. Then it holds that
kπ − π̂kπ∗ ≤

π > (P − P̂ )

π∗

1 − max {λ2 (P ∗ ), −λn (P ∗ )} − P − P̂

.
π∗

Proof. See Appendix A.1.
Several remarks regarding Theorem 5 are in order. First, in contrast to standard perturbation results
like Davis-Kahan’s sin Θ theorem, our theorem involves three matrices in total, where P , P̂ , and P ∗ can
all be arbitrary. For example, one may choose P ∗ to be the population transition matrix, and P and P̂ as
two finite-sample versions associated with P ∗ . Second, we only impose reversibility on P ∗ , whereas P and
P̂ need not induce reversible Markov Chains. Third, Theorem 5 allows one to derive the `2 estimation error
in Negahban et al. (2017a) directly without resorting to the power method; in fact, our `2 estimation error
bound improves upon Negahban et al. (2017a) by some logarithmic factor.
Theorem 6. Suppose p ≥ c0 logn n for some sufficiently large constant c0 > 0 and d ≥ cd np for cd ≥ 2 in
Algorithm 1. With probability exceeds 1 − O(n−5 ), one has
kπ − π ∗ k2
1
.√
.
kπ ∗ k2
npL
Proof. See Appendix D.

12

Notably, Theorem 6 matches the minimax lower bound derived in (Negahban et al., 2017a, Theorem 3).
As far as we know, this is the first result to demonstrate the orderwise optimality of the spectral method
when measured by the `2 loss.
The next result is concerned with the concentration of the vertex degrees in an Erdős–Rényi random
graph.
Lemma 1 (Degree concentration). Suppose that G ∼ Gn,p . Let di be the degree of node i, dmin = min1≤i≤n di
n
for some sufficiently large constant c0 > 0, then the following event
and dmax = max1≤i≤n di . If p ≥ c0 log
n


3np
np
≤ dmin ≤ dmax ≤
(30)
A0 =
2
2
obeys
P (A0 ) ≥ 1 − O(n−10 ).
Proof. This lemma follows from the standard Chernoff bound, and hence is omitted.
Since d is chosen to be cd np for some constant cd ≥ 2, we have, by Lemma 1, that the maximum vertex
degree obeys dmax < d with high probability.

3.2

Proof outline of Theorem 3
>

>

In this subsection, we outline the proof of Theorem 3. Recall that π = [π1 , · · · , πn ] and π ∗ = [π1∗ , · · · , πn∗ ]
are the stationary distributions associated with P and P ∗ , respectively. This gives
π> P = π>

and

π ∗> P ∗ = π ∗> .

For each 1 ≤ m ≤ n, one can decompose
>

∗
∗
∗
πm − πm
= π > P·m − π ∗> P·m
= π ∗> (P·m − P·m
) + (π − π ∗ ) P·m
X
X


∗
∗
=
πj∗ Pj,m − Pj,m
+ (πm − πm
) Pm,m +
πj − πj∗ Pj,m ,
j

|

j:j6=m

{z

:=I1m

}

{z

|

:=I2m

}

where P·m denotes the m-th column of P . Then it boils down to controlling I1m , I2m and

P

j:j6=m


πj − πj∗ Pj,m .

1. Since π ∗ is deterministic while P is random, we can easily control I1m using Hoeffding’s inequality.
The bound is the following.
Lemma 2. With probability exceeding 1 − O(n−5 ), one has
r
log n ∗
m
kπ k∞ .
max |I1 | .
m
Ld
Proof. See Appendix A.2.
∗
2. Next, the term I2m behaves as a contraction of |πm − πm
|, as demonstrated below.

Lemma 3. With probability exceeding 1 − O(n−5 ), there exists some constant c > 0 such that for all
1 ≤ m ≤ n,
!
r
np
log
n
∗
|I2m | ≤ 1 −
+c
|πm − πm
|.
2(1 + κ)d
Ld
Here, we recall that κ = wmax /wmin .
Proof. See Appendix A.3.
13

3. The statistical dependency
between
π and P introduces difficulty in obtaining a sharp estimate of

P
the third term j:j6=m πj − πj∗ Pj,m . Nevertheless, the leave-one-out technique helps us decouple
the dependency and obtain effective control of this term. The key component of the analysis is the
introduction of a new probability transition matrix P (m) , which is a leave-one-out version of the original
matrix P . More precisely, P (m) replaces all of the transition probabilities involving the m-th item with
their expected values (unconditional on G); that is, for any i 6= j,
(
Pi,j ,
i 6= m, j 6= m,
(m)
Pi,j := p ∗
y
,
i
= m or j = m,
d i,j
and for any 1 ≤ i ≤ n,
(m)

Pi,i := 1 −

X

(m)

Pi,j

j:j6=i
(m)

in order to ensure that P
is a probability transition matrix. In addition, we let π (m) be the
stationary distribution of the Markov chain induced by P (m) . As to be demonstrated later, the main
advantages of introducing π (m) are two-fold: (1) the original spectral estimate π can be very well
approximated by π (m) , and (2) π (m) is statistically independent of the connectivity
of the m-th node
P
and the comparisons with regards to the m-th item. Now we further decompose j:j6=m πj − πj∗ Pj,m :
X (m)
X
X

(m)
(πj − πj∗ )Pj,m .
πj − πj∗ Pj,m =
(πj − πj )Pj,m +
j:j6=m

j:j6=m

j:j6=m

{z

|

}

:=I3m

|

{z

:=I4m

}

4. For I3m , we apply the Cauchy-Schwarz inequality to obtain that with probability at least 1 − O(n−10 ),
|I3m | ≤ kπ (m) − πk2

X
j:j6=m

2
Pj,m

 12

(i)

1
≤ √ kπ (m) − πk2 ,
d

where (i) follows from the fact that Pj,m ≤ d1 for all j 6= m and dmax ≤ d on the event A0 . Consequently,
it suffices to control the `2 difference between the original spectral estimate π and its leave-one-out
version π (m) . This is accomplished in the following lemma.

Lemma 4. With probability at least 1 − O n−5 ,
√ r
16 κ log n ∗
(m)
kπ
− πk2 ≤
kπ k∞ + kπ − π ∗ k∞ ,
(31)
γ
Ld
where κ = wmax /wmin and γ = 1 − max {λ2 (P ∗ ), −λn (P ∗ )} − kP − P ∗ kπ∗ .
Proof. See Appendix A.4.
As shown in (Negahban et al., 2017a, Section 4.2), γ is lower bounded by some positive constant.
Lemma 5 (Spectral gap Negahban et al. (2017a)). Under the model specified in Section 2.1, if
p ≥ c0 logn n for some sufficiently large constant c0 > 0, then with probability at least 1 − O(n−5 ), one
has
1 dmin
1
γ := 1 − max {λ2 (P ∗ ), −λn (P ∗ )} − kP − P ∗ kπ∗ ≥ 2
≥ 2.
2κ dmax
6κ
With this spectral gap in mind, Lemma 4 reads
r
log n ∗
(m)
kπ
− πk2 .
kπ k∞ + kπ − π ∗ k∞ .
Ld

(32)

5. In order to control I4m , we exploit the statistical independence between π (m) and P·m . Specifically, we
demonstrate that:
14

Lemma 6. Suppose
 that p >
least 1 − O n−10 ,
|I4m |

c0 log n
n

for some sufficiently large constant c0 > 0. With probability at

1
. √ kπ (m) − πk2 +
n

r

log n ∗
kπ k∞ +
Ld

√

np log n + log n (m)
kπ
− π ∗ k∞ .
d

Proof. See Appendix A.5.
The above bound depends on both kπ (m) − πk2 and kπ (m) − π ∗ k∞ . We can invoke Lemma 4 and the
inequality kπ (m) − π ∗ k∞ ≤ kπ (m) − πk2 + kπ − π ∗ k∞ to reach
r
√
√


np log n + log n
log n ∗
np log n + log n
1
(m)
m
√
+
kπ
− πk2 +
kπ k∞ +
kπ − π ∗ k∞
|I4 | .
d
Ld
d
n
√

√
r
1
np log n + log n
κ
log n ∗
√ +
.
+1
kπ k∞
d
γ
Ld
n
√


1
np log n + log n
+ √ +
kπ − π ∗ k∞ .
d
n
6. Finally, putting the preceding bounds together, we conclude that with high probability, for some
absolute constants c1 , c2 , c3 > 0 one has
!
r
√


 √ r
np
logn
1
nplogn+logn 16 κ
logn ∗
c3
∗
−c1
|πm −πm |≤ c2 + √ + √ +c3
kπ k∞
2(1+κ)d
Ld
d
γ
Ld
n
d
√


1
c3
nplogn+logn
+ √ + √ +2c3
kπ−π ∗ k∞ .
d
n
d
simultaneously for all 1 ≤ m ≤ n. By taking the maximum over m on the left-hand side and combining
terms, we get
!
r
√
np
log n
np log n + log n
1
c3
− c1
− √ − √ − 2c3
kπ − π ∗ k∞
2(1 + κ)d
Ld
d
n
d
|
{z
}
:=α1


≤


c2 +

1
c
√ + √3 + c3
n
d

|

√

np log n + log n
d
{z

:=α2



√ r
16 κ
log n ∗
kπ k∞ .
γ
Ld
}

√
Hence, as long as np/ log n is sufficiently large, one has d  np  np log n > log n. This implies that
α1 > 0 and α1  α2  1, which further leads to
s
r
log
n
log n ∗
kπ − π ∗ k∞ .
kπ ∗ k∞ 
kπ k∞ .
Ld
npL
This finishes the proof of Theorem 3.

4

Analysis for the regularized MLE

The goal of this section is to establish the `∞ error of the regularized MLE as claimed in Theorem 4. Recall
∗
that in Theorem 4, we compare the regularized MLE θ with θ ∗ − θ 1, and hence without loss of generality
we can assume that
1> θ ∗ = 0.
(33)

15

4.1

Preliminaries and notation

Before proceeding to the proof, we gather some simple facts. To begin with, the gradient and the Hessian of
L (·; y) in (14) can be computed as

X 
eθ i
(ei − ej )
(34)
∇L (θ; y) =
−yj,i + θi
e + eθ j
(i,j)∈E,i>j

and
eθ i eθ j

X

∇2 L (θ; y) =

(i,j)∈E,i>j

(eθi

+

>

2
eθj )

(ei − ej ) (ei − ej ) ,

(35)

where e1 , · · · , en stand for canonical basis vectors in Rn .
When evaluated at the truth θ ∗ , the size of the gradient can be controlled as follows.
Lemma 7. The following event
(
A2 :=

r
∗

k∇Lλ (θ ; y)k2 .

n2 p log n
L

)
.

(36)

occurs with probability exceeding 1 − O(n−10 ).
Proof. See Appendix B.1.
The following lemmas characterize the smoothness and the strong convexity of the function Lλ (·; y).
P
>
In the sequel, we denote by LG =
the (unnormalized) Laplacian matrix
(i,j)∈E,i>j (ei − ej ) (ei − ej )
(Chung, 1997) associated with G, and for any matrix A we let

λmin,⊥ (A) := µ | z > Az ≥ µkzk2 for all z with 1> z = 0 .
(37)
Lemma 8. Suppose that p >
defined in (30), one has

c0 log n
n

for some sufficiently large constant c0 > 0. Then on the event A0 as


λmax ∇2 Lλ (θ; y) ≤ λ + np,
Proof. Note that

eθi eθj
θi +eθj 2
e
)
(

∀θ ∈ Rn .

≤ 41 . It follows immediately from the Hessian in (35) that

1
1
λmax ∇2 Lλ (θ; y) ≤ λ + kLG k ≤ λ + dmax ,
4
2

where dmax is the maximum vertex degree in the graph G. In addition, on the event A0 we have dmax ≤ 2np,
which completes the proof.
Lemma 9. For all kθk∞ ≤ C where C > 0 is some absolute constant, we have

λmin,⊥ ∇2 Lλ (θ; y) ≥ λ +

e−2C
(1 + e−2C )

2 λmin,⊥

(LG ) .

Proof. The proof is straightforward and hence is omitted.
n
for some sufficiently large constant c0 > 0. Then
Lemma 10. Let G ∼ Gn,p , and suppose that p > c0 log
n
one has



1
P λmin,⊥ (LG ) ≥ np ≥ 1 − O n−10 .
2

Proof. Note that λmin,⊥ (LG ) is exactly the spectral gap of the Laplacian matrix. See (Tropp, 2015, Sec
5.3.3) for the derivation of this lemma.
16

Algorithm 2 Gradient descent for regularized MLE
Initialize θ 0 = 0.
for t = 0, 1, 2, . . . , T − 1 do

θ t+1 = θ t − ηt ∇L (θ t ; y) + λθ t ;
end for
Output θ T .

By combining Lemma 9 with Lemma 10, we reach the following result.

Corollary 1. Under the assumptions of Lemma 10, with probability exceeding 1 − O n−10 one has

λmin,⊥ ∇2 Lλ (θ; y) ≥ λ +

e−2C

2 np

2 (1 + e−2C )

simultaneously for all θ obeying kθk∞ ≤ C.

4.2

Proof outline of Theorem 4

This subsection outlines the main steps for establishing Theorem 4. Before continuing, we single out an
important fact that 1> θ = 0, which is an immediate consequence of the shift-invariance property of L (θ; y).
Fact 1. The regularized MLE θ obeys 1> θ = 0.
Proof. Suppose instead that 1> θ 6= 0. Then setting θ =
kθk22 > kθ − θ1k22

1 >
n1 θ

one has


L (θ; y) = L θ − θ1; y ,

and

where the last
 identity arises since L (θ; y) is invariant under global shift. This implies that Lλ (θ; y) >
Lλ θ − θ1; y , which contradicts the optimality of θ.
The proof relies heavily on the strong convexity of the regularized negative log-likelihood function. Notably, while Lλ is a strongly convex function, its curvature might be arbitrarily close to λ for those points far
away from the optimizer. In contrast, if we restrict our attention to a set of parameters that are reasonably
close to the optimizer, a much stronger level of convexity can be exploited. In light of this observation, we
divide our analysis into two parts: (1) a coarse-level analysis demonstrating that θ is not far from θ ∗ in
an entrywise manner; (2) a fine-scale analysis that reveals the optimal performance of the regularized MLE
using refined strong convexity.
4.2.1

A coarse-level analysis

In this section, our goal is to obtain the crude bound
kθ − θ ∗ k∞ . 1.

(38)

Rather than directly resorting to the optimality condition, we adopt an algorithmic perspective to analyze
the regularized MLE θ. Specifically, we consider the standard gradient descent algorithm that is expected to
converge to the minimizer θ, and analyze the trajectory of this iterative algorithm instead. The algorithm
is stated in Algorithm 2.
Our proof can be separated into three parts. In what follows, we shall take
ηt ≡ η =

1
,
λ + np

∞

t = 0, 1, 2, · · · .

(39)

1. It is seen that the sequence {θ t }t=1 converges to the regularized MLE θ, a property that is standard
in convex optimization literature. This claim is summarized in the following lemma.

17

Lemma 11. On the event A0 as defined in (30), one has
kθ t − θk22 ≤ ct1 kθ 0 − θk22 ,
where c1 = 1 −

λ
λ+np .

Proof. The result directly follows from Lemma 8, Corollary 1, and the property of the gradient descent
algorithm (e.g. (Bubeck, 2015, Theorem 3.10)).
With this convergence property in place, if we can further show that kθ t − θ ∗ k∞ ≤ C for all t ∈ N,
then we obtain the desired result kθ − θ ∗ k∞ ≤ C.
2. Before moving to the analysis for `∞ errors, we develop an upper bound on kθ t − θ ∗ k2 .
Lemma 12. On the event A2 as defined in (36), there exists some constant C2 > 0 such that
√
t = 0, 1, 2, · · · .
kθ t − θ ∗ k2 ≤ C2 n,
Proof. See Appendix B.2.
3. We then turn to the `∞ error of θ t . To this end, we proceed by building a contraction relationship
t+1
∗
t
∗
between θm
− θm
and |θm
− θm
|. More precisely, the aim is to prove that
t+1
∗
t
∗
θm
− θm
≤ ρ θm
− θm
+ c2

for some ρ < 1 and c2 > 0, as detailed in the following lemma.
Lemma 13. On the event A0 ∩ A2 that occurs with probability at least 1 − O(n−10 ), we have
!
√
C2
e−2 3C2
t
∗
∗
t+1
∗
ηnp
θm − θm ≤ 1 − 0.3ηnp
√
2 θm − θm + 8ηnp + ηλkθ k∞ +
−2
3C
4
2
1+e
simultaneously for all 1 ≤ m ≤ n and t ∈ N. Here, C2 > 0 is the same constant as in Lemma 12.
Proof. See Appendix B.3.
Notice that θ 0 = 0 and hence kθ 0 − θ ∗ k∞ = kθ ∗ k∞ . 1. In view of Lemma 13, by picking C > 0 to
be any constant satisfying
√

e−2

1 − 0.3ηnp

1+

3C2

√
2
e−2 3C2

C + 8ηnp + ηλkθ ∗ k∞ +

C2
ηnp ≤ C,
4

(40)

one can conclude that
kθ t − θ ∗ k∞ ≤ C

=⇒

kθ t+1 − θ ∗ k∞ ≤ C.

4. By induction, the above bounds taken together guarantee that on the event A0 ∩ A2 , the iterates
satisfy kθ t − θ ∗ k∞ ≤ C throughout all iterations. In fact, this holds for all C > 0 compatible with the
condition (40), or equivalently,
8 + λ/np +
C≥
0.3

C2
4

This establishes our claim (38).

18

√

1 + e−2
e−2


3C2 2

√
3C2

 1.

Remark 1 (Near-optimal computational cost). Before we move on to refine the analysis, we point out
a byproduct of the preceding analysis: the near-optimal computational complexity of the gradient descent
algorithm. Specifically, the above analysis reveals that
kθ t − θ ∗ k∞ . 1

t = 0, 1, · · · ,

provided that the starting point is θ 0 = 0. This, combined with Corollary 1, suggests a much stronger level
of convexity within a region of interest, namely,



n
o
λmin,⊥ ∇2 Lλ θ̃; y
& λ + np,
∀θ̃ ∈ B := θ̃ | kθ̃ − θ ∗ k∞ . 1 ,
where the entire trajectory of {θ t | t ≥ 0} stays within this region B. Therefore, one can immediately improve
Lemma 11 to
kθ t − θk22 ≤ ct2 kθ 0 − θk22
for some absolute constant 0 < c2 < 1. This reveals an interesting phenomenon: even though the negative
log-likelihood function is not strongly convex (due to the possibly flat curvature of the logistic function) and
the regularization parameter λ is small, the gradient descent algorithm converges geometrically fast to the
optimizer, provided that we start from θ 0 = 0. Given that the gradient update in each iteration takes linear
time—in time proportional to reading the data, the whole algorithm has nearly linear runtime.
4.2.2

A fine-scale analysis

This subsection further refines our analysis using our results in Section 4.2.1. A key ingredient of the analysis
is to introduce a leave-one-out version of the MLE that replaces all log-likelihood components involving the
m-th item with their expected values (unconditional on G); that is,
(m)

θ (m) := arg minn Lλ
θ∈R

(θ; y) ,

where
(m)

Lλ

X

(θ; y) :=



−yj,i (θi − θj ) + log 1 + eθi −θj



(i,j)∈E,i>j,i6=m,j6=m


∗

1
eθ i
θi −θm
+ λkθk22 .
+
p − θ∗
∗ (θi − θm ) + log 1 + e
i + eθ m
2
e
i:i6=m
X



(41)

As before, this auxiliary vector θ (m) serves as a reasonably good proxy for θ, while being statistically
independent of {yi,m | (i, m) ∈ E}.
Our proof proceeds with the following steps, all of which rely on the fact that θ and θ (m) live in a “good”
region near θ ∗ where strong convexity holds, guaranteed by our coarse-level analysis in Section 4.2.1.
1. Show that both θ and θ (m) are close to the truth θ ∗ when measured by the `2 loss. Specifically, we
prove that:
Lemma 14. With probability exceeding 1 − O(n−5 ), one has
s
log n
kθ − θ ∗ k2 .
,
pL
and

s
θ (m) − θ ∗
2

.

log n
,
pL

1 ≤ i ≤ m.

(42)

(43)

q
√
n
When kθ ∗ k2  n, Lemma 14 yields bound on relative `2 error kθ − θ ∗ k2 /kθ ∗ k2 . log
npL of MLE,
the same rate (up to log factors) as the one of spectral method in Theorem 6.
19

Proof. See Appendix B.4.
2. Control the `2 difference between the MLE and its leave-one-out counterpart θ (m) as follows:
n
for some sufficiently large constant c0 > 0. With probability
Lemma 15. Suppose that p > c0 log
n
exceeding 1 − O(n−10 ), there exists some constant c1 > 0 such that
s
log n 1
(m)
kθ
− θk2 ≤ c1
+ kθ − θ ∗ k∞ ,
1 ≤ m ≤ n.
npL
3

Proof. See Appendix B.5.
Notably, while kθ (m) − θk2 depends also on the pointwise error kθ − θ ∗ k∞ , we emphasize the presence
of a contraction factor (i.e. 1/3) in front of kθ − θ ∗ k∞ . This factor plays a crucial role in showing error
contraction.
3. Show that the leave-one-out estimate θ (m) is close to the truth when it comes to the m-th coordinate:
Lemma 16. With probability exceeding 1 − O(n−10 ), we have
s
log n
(m)
∗
,
1 ≤ i ≤ m.
θm − θm .
npL
Proof. See Appendix B.6.
4. Taken collectively, the above bounds imply the existence of some absolute constant c2 > 0 such that
s
log n 1
∗
(m)
∗
(m)
+ kθ − θ ∗ k∞
(44)
θm − θm ≤ θm − θm + kθ − θ k2 ≤ c2
npL
3
with high probability. Since this holds for each 1 ≤ m ≤ n, we can take the maximum over m in the
left-hand side and apply the union bound to derive
s
log n 1
kθ − θ ∗ k∞ ≤ c2
+ kθ − θ ∗ k∞ ,
(45)
npL
3
or equivalently,
∗

kθ − θ k∞

3c2
≤
2

5. It remains to show that
∗

keθ − eθ k∞
.
keθ∗ k∞

s

s

log n
.
npL

(46)

log n
.
npL

Toward this end, we observe that for each 1 ≤ m ≤ n,
∗

∗

∗
∗
|eθm − eθm |
|eθ̃m (θm − θm
)|
eθmax +kθ−θ k∞ · |θm − θm
|
=
≤
,
eθmax
eθmax
eθmax
∗
where θ̃m is between θm and θm
, and θmax is the largest entry of θ∗ . Continuing the derivation and
using (46), we arrive at
s
∗
∗
|eθm − eθm |
eθmax +kθ−θ k∞ kθ − θ ∗ k∞
log n
max
≤
.
1≤m≤n
eθmax
eθmax
npL

as long as npL/ log n is sufficiently large. This completes the proof of Theorem 4.
20

Connection with the spectral method. Before concluding this section, we remark that the proof
strategy for the regularized MLE shares a lot of similarity with that for the spectral method. For instance,
the proofs in both cases proceed by introducing appropriate leave-one-out estimates, which are exceedingly
close to the original estimates in the sense that
s
s
log
n
log n
kπ ∗ k∞
and
kθ (m) − θk2 .
.
kπ (m) − πk2 .
npL
npL
Furthermore, a crucial ingredient of the proofs is some sort of the perturbation bound. This corresponds
to the counterpart of the Davis-Kahan theorem (Theorem 5) for the spectral method, and the optimality
condition for the regularized MLE (see, e.g. (60)) under strong convexity. Both cases take the following
form: for an optimal solution z and its perturbed version z̃, one has
kz − z̃k2 .

kf (z̃)k2
,
γ

(47)

where γ represents some stability constant (spectral gap in spectral method or strong convexity parameter
in MLE). Moreover, f (·) is some vector-valued function obeying f (z) = 0 that captures the optimality of z;
more precisely, one can take
(
z > P − z, for the spectral method;
f (z) : =
∇Lλ (z),
for the regularized MLE.
We expect this form (47) to extend to a much broader class of optimization problems.

5

Discussion

In this paper, we justify the optimality of both the spectral method and the regularized MLE for top-K
rank aggregation. Our theoretical studies are by no means exhaustive, and there are numerous directions
that would be of interest for future investigations. We point out a few possibilities as follows.
• Unregularized MLE.
We have studied the optimality of the regularized MLE with the regularization
q
np log n
. Our analysis relies on the regularization term to obtain the coarse bound on
parameter λ 
L
kθ − θ ∗ k∞ in Section 4.2.1. It is natural to ask whether such a regularization term is necessary or not.
This question remains open.
• More general comparison graphs. So far we have focused on a tractable but somewhat restrictive
comparison graph, namely, the Erdős–Rényi random graph. It would certainly be interesting and
important to understand the performance of both methods under a broader family of comparison
graphs, and to see which algorithms would enable optimal sample complexities under general sampling
patterns.
• Beyond pairwise comparisons. In real world applications, we often encounter the case when we are
given more than two items to compare. For instance, in online advertisement scenarios, the users might
be asked to compare three to five different advertisements. How to extend the methods and analyses
to handle more general listwise comparisons is of significant practical interest.
• Entrywise perturbation analysis for general convex optimization. The current paper provides the `∞
perturbation analysis for the regularized MLE using the leave-one-out trick as well as an inductive
argument along the algorithmic updates. We expect this analysis framework to carry over to a much
broader family of convex optimization problems, which may in turn offer a powerful tool for showing
the stability of optimization procedures in an entrywise fashion.

21

A

Proofs in Section 3

This section collects proofs of the theorems and lemmas that appear in Section 3. By Lemma 1, the event


1
3
A0 =
np ≤ dmin ≤ dmax ≤ np
2
2
happens with probability at least 1 − O(n−10 ). Throughout this section, we shall assume that we are on this
event without explicitly referring to it each time. An immediate consequence is that dmax ≤ d on this event.

A.1

Proof of Theorem 5

To begin with, we write
>

π > − π̂ > = π > P − π̂ > P̂ = π > (P − P̂ ) + (π − π̂) P̂ .

(48)

The last term of the above identity can be further decomposed as
>

>

>

(π − π̂) P̂ = (π − π̂) P ∗ + (π − π̂) (P̂ − P ∗ )

>
>
= (π − π̂) P ∗ − 1π ∗> + (π − π̂) (P̂ − P ∗ ),

(49)

>

where we have used the fact that (π − π̂) 1π ∗> = 0. Combining (48) and (49) we get

>
>
π > − π̂ > = π > (P − P̂ ) + (π − π̂) P ∗ − 1π ∗> + (π − π̂) (P̂ − P ∗ ),
which together with a little algebra gives
kπ − π̂kπ∗ ≤ π > (P − P̂ )

π∗

+ kπ − π̂kπ∗ P ∗ − 1π ∗>

kπ − π̂kπ∗ ≤

=⇒

π∗

π > (P − P̂ )
1−

kP ∗

−

1π ∗> k

π∗

+ kπ − π̂kπ∗ P̂ − P ∗

π∗

− P̂ − P ∗

π∗

.
π∗

The theorem follows by recognizing that
P ∗ − 1π ∗>

A.2

π∗

= max {λ2 (P ∗ ), −λn (P ∗ )} .

Proof of Lemma 2

Observe that
I1m =

X



∗
∗
∗
πj∗ Pj,m − Pj,m
+ πm
Pm,m − Pm,m

j:j6=m
(i)

=

(
X

∗
∗
+ πm
πj∗ Pj,m − Pj,m



!
1−

j:j6=m

=

X

X

Pm,j

j:j6=m

πj∗

+

∗
πm



Pj,m −

∗
Pj,m

!)
−

1−

X

∗
Pm,j

j:j6=m



j:j6=m

=

L



1 X X ∗
(l)
∗
∗
πj + πm
1(j,m)∈E yj,m − yj,m
,
Ld

(50)

j:j6=m l=1

where (i) follows from the fact that P and P ∗ are both probability transition matrices. By Lemma 17 and
the fact that dmax ≤ d, one can derive


L
 X X

n
o



(l)
∗
∗
P |I1m | ≥ t G = P
πj∗ + πm
1(j,m)∈E yj,m − yj,m
≥ Ldt G


j:j6=m l=1


2(Ldt)2
≤ 2 exp −
.
Ldmax (2kπ ∗ k∞ )2
22

This together with the high probability event A0 indicates


n
o
−Ldt2
m
,
P |I1 | ≥ t ≤ 2 exp
2kπ ∗ k2∞
which further implies
(

r

|I1m |

P

≥4

log n ∗
kπ k∞
Ld

)
≤ 2n−8 .

The lemma is established by taking the union bound.

A.3

Proof of Lemma 3

Applying Lemma 17 to the quantity
X

∗
Pm,m − Pm,m
=−

∗
(Pm,j − Pm,j
)=−

L


1 X X
(l)
∗
,
1(j,m)∈E yj,m − yj,m
Ld
j:j6=m l=1

j:j6=m

we get
(

r

∗
P max |Pm,m − Pm,m
|≥2
m

log n
Ld

)
≤ 2n−7 .

On the other hand, we have for all 1 ≤ m ≤ n
∗
Pm,m
=1−

X

∗
Pm,j
≤1−

j:j6=m

dmin
1
np
1
·
≤1−
·
.
d
1+κ
2d 1 + κ

Combining these two pieces completes the proof.

A.4

Proof of Lemma 4

First, by the relationship between k · k2 and k · kπ∗ , we have
kπ (m) − πk2 ≤ p

1
kπ (m) − πkπ∗ ,
∗
πmin

∗
where πmin
:= mini πi∗ . Invoking Theorem 5, we obtain

kπ (m)> (P (m) − P )kπ∗
1 − max {λ2 (P ∗ ), −λn (P ∗ )} − kP − P ∗ kπ∗
1p ∗
πmax kπ (m)> (P (m) − P )k2 ,
≤
γ

kπ (m) − πkπ∗ ≤

∗
where we define γ := 1 − max {λ2 (P ∗ ), −λn (P ∗ )} − kP − P ∗ kπ∗ and πmax
:= maxi πi∗ .
(m)>
(m)
To facilitate the analysis of kπ
(P
− P )k2 , we introduce another Markov chain with transition
probability matrices P (m),G , which is also a leave-one-out version of the transition matrix P . Similar to
P (m) , P (m),G replaces all the transition probabilities involving the m-th item with their expected values
(conditional on G). Concretely, for i 6= j,
(
Pi,j ,
i 6= m, j 6= m,
(m),G
Pi,j
= 1 ∗
i = m or j = m.
d yi,j 1(i,j)∈E ,

And for each 1 ≤ i ≤ n, we define
(m),G

Pi,i

=1−

X
j:j6=i

23

(m),G

Pi,j

to make P (m),G a valid probability transition matrix. Hence by the triangle inequality, we see that
kπ (m)> (P (m) − P )k2 ≤ kπ (m)> (P − P (m),G )k2 + kπ (m)> (P (m) − P (m),G )k2 .
{z
} |
{z
}
|
:=J1m

:=J2m

The next step is then to bound J1m and J2m separately.
For J1m , similar to (50), one has
h
i (i) h
i
π (m)> (P − P (m),G )
= π (m)> (P − P ∗ )
m

m

1 X
Ld

=

L 
X

(m)

πj




(l)
∗
(m)
,
1(j,m)∈E yj,m − yj,m
+ πm

j:j6=m l=1

(m),G

∗
where (i) comes from the fact that P·m
= P·m
. Recognizing that π (m) is statistically independent of
{yj,m }j6=m , by Hoeffding’s inequality (cf. Lemma 17), we get
(
)
r
log n (m)
(m)>
(m),G
(P − P
)]m ≥ 4
P [π
kπ k∞ ≤ 2n−8 .
(51)
Ld

And for j 6= m, we have
h
X (m)
i
(m),G 
π (m)> P − P (m),G
=
πi
Pi,j − Pi,j
j

i
(m)
πj

(m),G 
(m),G 
(m)
Pj,j − Pj,j
+ πm
Pm,j − Pm,j
(m) 1 ∗
(m) 1
∗
(y
− yj,m ) 1(j,m)∈E + πm
(ym,j − ym,j
) 1(j,m)∈E
= πj
d j,m
d


(m)
∗
(m)
∗
= πj
Pj,m
− Pj,m + πm
Pm,j − Pm,j
.

=

In addition by Hoeffding’s inequality in Lemma 17, we have
max Pj,m −
j6=m

∗
Pj,m

2
≤
d

r

log n
L

with probability at least 1 − O(n−5 ). As a consequence,
( q
h
log n
4
(m)
i
k∞ ,
(m)>
(m),G
L kπ
≤ d
π
P −P
j
0,

if (j, m) ∈ E,
else.

Combining (51) and (52) yields
r
r
r
√
(i)
log n (m)
4 dmax − 1 log n (m)
log n (m)
m
kπ k∞ +
kπ k∞ ≤ 8
kπ k∞ ,
J1 ≤ 4
Ld
d
L
Ld
where (i) comes from the fact that dmax ≤ d.
Regarding J2m , we invoke the identity π ∗> (P (m) − P (m),G ) = 0 to get

>
π (m)> (P (m) − P (m),G ) = π (m) − π ∗ (P (m) − P (m),G ).
Therefore, for j 6= m we have


>

(m)
∗
(m)
(m),G
π
−π
P
−P
j

X (m)
(m)
(m),G 
=
(πi − πi∗ ) Pi,j − Pi,j
i
(m)

(m)
(m),G 
(m)
∗
+ (πm
− πm
) Pm,j − Pm,j
(m)
(m)
(m),G 
(m)
(m),G 
(m)
∗
= −(πj − πj∗ ) Pj,m − Pj,m
+ (πm
− πm
) Pm,j − Pm,j .

= (πj

(m)

(m),G 

− πj∗ ) Pj,j − Pj,j

24

(52)

(m)

(m),G

Recognizing that Pj,m − Pj,m

≤

(m)

2
d

(m),G

for (j, m) ∈ E and Pj,m − Pj,m
≤
(

i
h
4
kπ (m) − π ∗ k∞ ,
≤ d2p (m)
π (m) − π ∗ (P (m) − P (m),G )
j
− π ∗ k∞ ,
d kπ

p
d

for (j, m) ∈
/ E, we have
if (j, m) ∈ E,
else.

(53)

And for j = m, it holds that


>

π (m) − π ∗
P (m) − P (m),G
m

X (m)

(m)
(m),G 
(m)
∗
(m)
(m),G
+
(πj − πj∗ ) Pj,m − Pj,m
− πm
) Pm,m
− Pm,m
= (πm

(54)

j:j6=m
(m)
∗
(m)
(m),G
≤ (πm
− πm
) Pm,m
− Pm,m



X

+

(m)

(m)

(m),G 

− πj∗ ) Pj,m − Pj,m

(πj

j:j6=m

=

X

(m)

(m),G 

(m)
∗
(πm
− πm
) Pm,j − Pm,j

+

j:j6=m

(m)

(m),G

(m)

(πj

(m)

(m),G 

− πj∗ ) Pj,m − Pj,m

.

(55)

j:j6=m

|
Given that Pm,j − Pm,j

X

{z

}

:=J3m

=

∗
ym,j
d (p

|

{z

:=J4m

}

− 1(m,j)∈E ), we have

J3m =

∗
ym,j
(m)
∗
(πm
− πm
)
(p − 1(m,j)∈E ).
d}
{z
j:j6=m |

X

(56)

(m)

:=ξj

Since kξ (m) k∞ ≤ d1 kπ (m) − π ∗ k∞ and kξ (m) k2 ≤ d1 kπ (m) − π ∗ k2 , Lemma 18 implies that
√
np log n + log n (m)
|J3m | .
kπ
− π ∗ k∞
d
with high probability. The same bound holds for J4m . Combine (53), (55) and (56) to arrive at
√ !
√
√
np
log
n
+
log
n
p
n
d
J2m .
+
+
kπ (m) − π ∗ k∞ .
d
d
d
Combining all, we deduce that
√
κ m m
(m)
(J1 +J2 )
kπ −πk2 ≤
γ
!
r
√ !
√
√
√
κ
logn (m)
nplogn+logn p n
d
(m)
∗
≤
8
kπ k∞ +C
+
+
kπ −π k∞
γ
Ld
d
d
d
!
r
r
√ !
√
√
√
κ
logn ∗
logn
nplogn+logn p n
d
≤
8
kπ k∞ +C 8
+
+
+
kπ (m) −π ∗ k∞
γ
Ld
Ld
d
d
d
√ r
(i) 8 κ
logn ∗
1
≤
kπ k∞ + kπ (m) −π ∗ k∞ ,
γ
Ld
2
where (i) holds as long as p ≥ c0 n log n for c0 sufficiently large. The triangle inequality
kπ (m) − π ∗ k∞ ≤ kπ (m) − πk2 + kπ − π ∗ k∞
yields
kπ

(m)

√ r
16 κ log n ∗
kπ k∞ + kπ − π ∗ k∞ ,
− πk2 ≤
γ
Ld

which concludes the proof.
25

(57)

A.5

Proof of Lemma 6

For ease of presentation, we define
L

ỹi,j :=

1 X (l)
ỹi,j
L
l=1

for all i 6= j, where
(
(l) ind.
ỹi,j =

1, with probability

wj∗
wi∗ +wj∗ ,

0, else.

This allows us to write y as yi,j = ỹi,j 1(i,j)∈E for all i 6= j. With this notation in place, we can obtain
!
L
X (m)
X
1
(l)
I4m =
(πj − πj∗ )
ỹj,m 1(j,m)∈E .
Ld
j:j6=m

We can further decompose

I4m

l=1

into
I4m = E [I4m | G−m , ỹ] + I4m − E [I4m | G−m , ỹ] ,

where G−m represent the graph without the m-th node, and ỹ = {ỹi,j |i 6= j} represents all the binary
outcomes.
We start with the expectation term
!
L
X (m)
1 X (l)
m
∗
ỹj,m P {(j, m) ∈ E}
E [I4 | G−m , ỹ] =
(πj − πj )
Ld
l=1
j:j6=m
√ (ii)
(i)
1
n
(m)
∗
≤ √ kπ (m) − π ∗ k2
≤ pkπ
− π k2
d
2 n
(iii)

≤

1
√ (kπ (m) − πk2 + kπ − π ∗ k2 ),
2 n

where (i) comes from the Cauchy-Schwarz inequality, (ii) follows from the choice d = cd np ≥ 2np and (iii)
results from the triangle inequality. By Theorem 6, with high probability we have
r
r
log n ∗
log n √
∗
nkπ ∗ k∞ ,
kπ − π k2 ≤ CN
kπ k2 ≤ CN
Ld
Ld
thus indicating that
E [I4m

1
CN
| G−m , ỹ] ≤ √ kπ (m) − πk2 +
2
2 n

r

log n ∗
kπ k∞ .
Ld

(58)

When it comes to the fluctuation term, one can write
I4m − E[I4m | G−m , ỹ] =

X

(m)
(πj

j:j6=m

|

!
L
X

1
(l)
− πj∗ )
ỹj,m 1(j,m)∈E −P {(j, m) ∈ E} .
Ld
l=1
{z
}
(m)

:=βj

Since kβ (m) k2 ≤ d1 kπ (m) − π ∗ k2 and kβ (m) k∞ ≤ d1 kπ (m) − π ∗ k∞ , one can apply Lemma 18 to derive
√


I4m − E I4m | G−m , ỹ .

np log n + log n (m)
kπ
− π ∗ k∞
d

with high probability.
The bounds (58) and (59) taken together complete the proof.
26

(59)

B

Proofs in Section 4

This section gathers the proofs of the lemmas in Section 4.

B.1

Proof of Lemma 7

Observe that
∗



X

∗

∇Lλ (θ ; y) = λθ +

(i,j)∈E,i>j

1
= λθ +
L
∗

∗

eθi
−yj,i + θ∗
∗
e i + eθ j


(ei − ej )


∗
L 
X
eθ i
(l)
(ei − ej ).
−yj,i + θ∗
∗
e i + eθ j
(i,j)∈E,i>j l=1 |
{z
}
X

(l)

:=zi,j
(l)

(l)

It is seen that E[zi,j ] = 0, kzi,j k ≤

√

2,

h
i
h i
(l) (l)>
(l)
>
>
E zi,j zi,j = Var yi,j (ei − ej ) (ei − ej )  (ei − ej ) (ei − ej )

and

h
i
 h
i
(l)> (l)
(l) (l)>
E zi,j zi,j = Tr E zi,j zi,j
≤ 2.

This implies that with high probability (note that the randomness comes from G),
X

L
h
i
X
(l) (l)>
≤L
E zi,j zi,j

(i,j)∈E,i>j l=1

X

>

(ei − ej ) (ei − ej )

= LkLG k . Lnp

(i,j)∈E,i>j

and
X

L
h
i
X
(l)> (l)
E zi,j zi,j ≤ 2L

(i,j)∈E,i>j l=1

n P

X

1 . Ln2 p.

(i,j)∈E,i>j

h

i

h
io
PL
(l) (l)>
(l)> (l)
(l)
Letting V := L12 max
E
z
z
,
E
z
z
and B := maxi,j,l kzi,j k,
i,j i,j
i,j
i,j
(i,j)∈E
l=1
(i,j)∈E
l=1
we can invoke the matrix Bernstein inequality (Tropp, 2012, Theorem 1.6) to reach
r
p
n2 p log n log n
∗
∗
k∇Lλ (θ ; y) − E [∇Lλ (θ ; y) | G]k2 . V log n + B log n .
+
L
L
r
2
n p log n
.
L
PL

P

with probability at least 1 − O(n−10 ). Combining this with the identity E [∇Lλ (θ ∗ ; y) | G] = λθ ∗ yields
k∇Lλ (θ ∗ ; y)k2 ≤ kE [∇Lλ (θ ∗ ; y) | G]k2 + k∇Lλ (θ ∗ ; y) − E [∇Lλ (θ ∗ ; y) | G]k2
r
r
n2 p log n
n2 p log n
∗
. λkθ k +

,
L
L
thus concluding the proof.

B.2

Proof of Lemma 12

First, Lemma 11 reveals that
kθ t − θk2 ≤ kθ 0 − θk2 = kθk2 ≤ kθ − θ ∗ k2 + kθ ∗ k2
27

for all t ≥ 0. Thus, by the triangle inequality we have
kθ t − θ ∗ k2 ≤ kθ t − θk2 + kθ − θ ∗ k2 ≤ 2kθ − θ ∗ k2 + kθ ∗ k2 .
Next, it follows from the optimality of θ as well as the mean value theorem that
Lλ (θ ∗ ; y) ≥ Lλ (θ; y) = Lλ (θ ∗ ; y) + h∇Lλ (θ ∗ ; y) , θ − θ ∗ i +


1
>
(θ − θ ∗ ) ∇2 Lλ θ̃; y (θ − θ ∗ ) ,
2

where θ̃ is between θ and θ ∗ . This together with the Cauchy-Schwarz inequality gives

1
>
(θ − θ ∗ ) ∇2 Lλ θ̃; y (θ − θ ∗ ) ≤ − h∇Lλ (θ ∗ ; y) , θ − θ ∗ i ≤ k∇Lλ (θ ∗ ; y)k2 kθ − θ ∗ k2 .
2
Recall from Fact 1 and the assumption (33) that 1> θ = 1> θ ∗ = 0. Then the above inequality gives
kθ − θ ∗ k2 ≤

2 k∇Lλ (θ ∗ ; y)k2

 ,
λmin,⊥ ∇2 Lλ θ̃; y

(60)

where λmin,⊥ (A) is the largest quantityof µ obeying
x> Ax ≥ µkxk22 for all x obeying 1> x = 0.


2
From the trivial lower bound λmin,⊥ ∇ Lλ θ̃; y
≥ λ, the preceding inequality gives
kθ − θ ∗ k2 ≤

2 k∇Lλ (θ ∗ ; y)k2
λ

(61)

and therefore

4 k∇Lλ (θ ∗ ; y)k2
+ kθ ∗ k2 .
λ


q
q
2
np log n
On the event A2 = k∇Lλ (θ ∗ ; y)k2 . n pLlog n and in the presence of the choice λ 
, we
L
√
obtain kθ t − θ ∗ k2 ≤ C2 n for all t ≥ 0.
kθ t − θ ∗ k2 ≤

B.3

Proof of Lemma 13

The update rule of Algorithm 2 is given by



t+1
∗
t
∗
θm
− θm
= θm
− θm
− η ∇L θ t ; y + λθ t m .

(62)

From the gradient (34), we can rewrite (62) as
(
t+1
θm

−

∗
θm

=

t
θm

=

I1m

−

∗
θm

+

I2m

X

−η

i:(m,i)∈E

+

I3m

t

eθ m
−yi,m + θt
t
e i + eθ m

)
t
− ηλθm

+ I4m ,

where the four terms are defined as
X 
∗
I1m = −η
−yi,m + yi,m
,
i:(m,i)∈E

)
(
)
t
t
θm
θm
X
e
e
∗
∗
,
= −η 
−yi,m
+ θt
−p
−yi,m
+ θt
t
t
i + eθ m
i + eθm
e
e
i6=m
i:(m,i)∈E
(
)
t
θ
X

em
t
∗
∗
t
∗
= θm
− θm
− ηp
−yi,m
+ θt
− ηλ θm
− θm
,
t
θ
e i +e m
i6=m


I2m
I3m

(

X

∗
I4m = −ηλθm
.

28

(63)

Here, we set

∗

eθm
:= θ∗
∗ .
e i + eθm
It then comes down to analyzing the
 preceding four terms.
3np
we have
First of all, on the event A0 := np
2 ≤ dmin ≤ dmax ≤ 2
∗
yi,m

I1m ≤ 2ηdmax ≤ 3ηnp
and
I2m ≤ η (2dmax + 2np) ≤ 5ηnp.
In addition, I4m is simple and bounded by
I4m ≤ ηλkθ ∗ k∞ .
Hence, we are only left with I3m .
To analyze I3m , we define the population loss
X

∗
L∗ (θ) :=
−yj,i
(θi − θj ) + log 1 + eθi −θj ,
i>j

whose gradient and Hessian can be computed as
X
∗
∇L∗ (θ) =
−yj,i
+

eθ i
θ
e i + eθ j

i>j

∇2 L∗ (θ) =

X

eθi eθj

i>j

(eθi + eθj )


(ei − ej ) ;

(64)

>

(ei − ej ) (ei − ej ) .

(65)

With these notation in place, one can simplify I3m as



t
∗
t
∗
I3m = θm
− θm
− ηp ∇L∗ (θ t ) m − ηλ θm
− θm
.

(66)

2

∗

∗

Furthermore, it is straightforward to see that ∇L (θ ) = 0, which together with Taylor’s theorem yields


[∇L∗ (θ t )]m = [∇L∗ (θ t )]m − [∇L∗ (θ ∗ )]m = [∇2 L∗ θ̃ t ]m,· θ t − θ ∗ ,
where θ˜t lies between θ t and θ ∗ . Combine the preceding two identities to get



h
h
i
 X
i
m
2 ∗
t
t
∗
−ηp ∇2 L∗ θ̃ t
I3 = (1 − ηλ) − ηp ∇ L θ̃
θm − θm +
m,m

m,i

i6=m

{z

|

}

:=J1m



{z

|

:=J2m


θit − θi∗ .
}

Therefore, it remains to control J1m and J2m .
For J1m , we have
h
i
J1m ≤ (1 − ηλ) − ηp ∇2 L∗ θ̃ t
h
i
Hence we need to lower bound ∇2 L∗ θ̃ t
m,m

on the event A2 and, as a result,

m,m

t
∗
· θm
− θm
.

√
. Recall from Lemma 12 that we have kθ t − θ ∗ k2 ≤ C2 n

√
kθ˜t − θ ∗ k2 ≤ C2 n.

A straightforward consequence is that1 at least 60% of the entries in θ̃ t − θ ∗ must obey θ̃it − θi∗ ≤
With this observation in mind, one has
t

h
i
∇2 L∗ θ̃ t
m,m

=

t

i:i6=m

t

eθ̃i + eθ̃m

we would have kθ˜t − θ ∗ k22 ≥ 0.4n ×
√
kθ˜t − θ ∗ k2 ≤ C2 n.
1 Otherwise,

√

√

t

eθ̃i eθ̃m

X

3C2

2 ≥ 0.6 (n − 1)

2

29

e−2
1+

√

3C2 .

3C2

√
2 .
e−2 3C2

= 1.2C22 n > C22 n, which is contradictory to the fact that

This gives

√

J1m

≤

1 − 0.3ηnp

e−2

!

3C2

√

1 + e−2

t
∗
θm
.
− θm


3C2 2

When it comes to J2m , we see that
J2m



X h
i
∇2 L∗ θ̃ t
= ηp
m,i

i:i6=m



1/2 


≤ ηp 

θit − θi∗

X h

∇2 L∗ θ̃

i 2
t
m,i

i:i6=m

1/2



h

∇2 L∗ θ̃ t

i



i:i6=m

1/2
X h

∇2 L∗ θ̃

i 2
t
m,i

i:i6=m

Hence we need to upper bound

θit − θi




≤ ηpkθ t − θ ∗ k2 


∗ 2

X

.



for i 6= m. By (65), we have
m,i
t

h

∇2 L∗ θ̃ t

i

=
m,i

t

eθ̃i eθ̃m
t
eθ̃i

+

2
t
eθ̃m

≤

1
,
4

which gives us

√
C2
1 √
ηp n − 1C2 n ≤
ηnp.
4
4
Combing all of the above bounds, we deduce that
!
√
e−2 3C2
C2
t+1
∗
t
∗
∗
θm − θm ≤ 1 − 0.3ηnp
ηnp
√
2 θm − θm + 8ηnp + ηλkθ k∞ +
−2
3C
4
2
1+e
|J2m | ≤

on the event A0 ∩ A2 . This finishes the proof.

B.4

Proof of Lemma 14

We first provide here a proof for the bound on kθ − θ ∗ k2 .
Recall from (60) that
2 k∇Lλ (θ ∗ ; y)k2


 ,
kθ − θ ∗ k2 ≤
λmin,⊥ ∇2 Lλ θ̃; y

(67)

where θ̃ is between θ and θ ∗ . We have also learned from Section 4.2.1 that kθ − θ ∗ k∞ . 1, and hence
θ̃ − θ ∗

∞

≤ kθ − θ ∗ k∞ . 1.

(68)

This together with Corollary 1 reveals that



λmin,⊥ ∇2 Lλ θ̃; y
& np,
and therefore
kθ − θ ∗ k2 .

k∇Lλ (θ ∗ ; y)k2
.
np

Therefore, combining Lemma 7 with (67) leads to
k∇Lλ (θ ∗ ; y)k2
kθ ∗ − θk2 .
.
np
30

s

log n
pL

(69)

as claimed.
The proof concerning θ (m) − θ ∗
that

2

follows an almost identical argument, except that we need to establish
kθ (m) − θ ∗ k∞ . 1.

And this can be easily checked using the same argument in Section 4.2.1.

B.5

Proof of Lemma 15

Repeating the argument in Appendix B.4 gives
kθ

(m)



∇Lλ θ (m) ; y
2 ∇Lλ θ (m) ; y 2


 .
− θk2 ≤
np
λmin ∇2 Lλ θ̃ (m) ; y

It then suffices to upper bound ∇Lλ θ (m) ; y

2

.

(70)



. Recognizing that

(m)
∇Lλ θ (m) ; y = 0,
2

we can deduce that



(m)
∇Lλ θ (m) ; y = ∇Lλ θ (m) ; y − ∇Lλ θ (m) ; y
(
!
!)
(m)
(m)
∗
X
eθ i
eθi
eθ i
=
−ym,i + (m)
1{(i,m)∈E} −p − θ∗
(ei − em )
∗ +
(m)
(m)
(m)
e i + eθ m
eθi + eθm
eθi + eθm
i:i6=m

∗
L 
1 X X
eθi
(l)
=
−ym,i + θ∗
(ei − em )
∗
L
e i + eθ m
i:(i,m)∈E l=1
|
{z
}
:=um
(
!
)
(m)
∗
X

eθ i
eθ i
(71)
+
− θ∗
1{(i,m)∈E} −p (ei − em ).
∗ +
(m)
(m)
e i + eθ m
eθi + eθm
i:i6=m
|
{z
}
:=v m

In the sequel, we control the two terms of (71) separately.
For the first term um in (71), we make the observation that

 P 
∗
(l)
L
1
eθi

,
if (i, m) ∈ E;
∗
 L l=1 −ym,i + eθi∗ +eθm



∗
P
P
θ
m
(l)
e i
ui = 1 i:(i,m)∈E L
, if i = m;
∗
l=1 ym,i − eθi∗ +eθm
L



0,
else.
(l)

θ∗

Since ym,i − θi∗e i θm
≤ 1 and card ({i : (i, m) ∈ E})  np, we can apply Hoeffding’s inequality and union
∗
e +e
bounds to get for all 1 ≤ m ≤ n,
r
r
np log n
log n
m
m
|um | .
and
|ui | .
∀i obeying (i, m) ∈ E,
L
L
which further gives
m

ku k2 ≤

|um
m|

+

s X

r
2
(um
i ) .

i:(i,m)∈E

np log n
, ∀1 ≤ m ≤ n.
L

We then turn to the second term v m of (71). This is a zero-mean random vector that satisfies


− p) ,
if (i, m) ∈ E,
ξi (1

P
m
vi = − i:i6=m ξi 1{(i,m)∈E} −p , if i = m,


−ξi p,
else,
31

where

(m)

∗

1
eθ i
1
eθ i
=−
+
ξi := − θ∗
.
∗ −θ ∗ +
∗
(m)
(m)
(m)
(m)
θm
i
θ
θ
θ
e i + eθ m
1
+
e
e i +e m
1 + e m −θi
The first step is to bound the size of the coefficient ξi . Define g(x) = (1 + ex )−1 for x ∈ R. We have
|g 0 (x)| ≤ 1 and thus
(m)

(m)
|ξi | = |g(θm
− θi

(m)

∗
(m)
) − g(θm
− θi∗ )| ≤ |(θm
− θi

This indicates that
|ξi | ≤ 2kθ (m) − θ ∗ k∞

(m)

∗
) − (θm
− θi∗ )| ≤ θi∗ − θi

n
X

and

∗
(m)
.
− θm
+ θm

ξi2 ≤ 4nkθ (m) − θ ∗ k2∞ .

i=1

Applying the Bernstein inequality in Lemma 18 we obtain
v
!
u
n
p

u X
m
2
t
|vm | .
p
np log n + log n kθ (m) − θ ∗ k∞
ξi log n + max |ξi | log n .
1≤i≤n

i=1

with high probability. As a consequence,
s
s X
m
(vim )2 +
kv m k2 ≤ |vm
|+

X

(vim )2

i:(i,m)∈E
/ and i6=m

i:(i,m)∈E


p
√
√
.
np log n + log n kθ (m) − θ ∗ k∞ + npkθ (m) − θ ∗ k∞ + p nkθ (m) − θ ∗ k∞
p

.
np log n + log n kθ (m) − θ ∗ k∞ .
Putting the above results together, we see that
kθ (m) − θk2 .
.

kum k2 + kv m k2
np
q

√
np log n
+
np log n + log n kθ (m) − θ ∗ k∞
L
np

(72)

.

(73)

Upper bounding kθ (m) − θ ∗ k∞ by kθ − θ ∗ k∞ + kθ − θ (m) k, one has

q

√
 (m)
√

 np log n +
np log n + log n kθ − θ ∗ k∞
np
log
n
+
log
n
kθ
−
θk
L
2
kθ (m) − θk2 ≤ c1
+


np
np
c0 log n
for some
for some sufficiently large constant c0 > 0 so that np >
n
√ constant c1 > 0. When p >
4c1
np log n + log n , the above inequality guarantees that
c1

√ np log n
L

kθ (m) − θk2 ≤
1

B.6

√
+c1 ( np log n+log n)kθ−θ ∗ k∞
np
√
c1 ( np log n+log n)
−
np

≤

4c1

q

np log n
L

3np

1
+ kθ − θ ∗ k∞ .
3

Proof of Lemma 16
(m)

From the definition of θ (m) , we know that θm
entries of θ (m) , that is,

is also the coordinate-wise optimizer when given all other

(m)
θm
= arg min `(m) (τ ; y) ,

32

where
`(m) (τ ; y) = p

X 

−

i:i6=m

∗
 λ



(m)
eθ i
(m)
θi −τ
+ τ2
θ
−
τ
+
log
1
+
e
∗
∗
2
eθi + eθm i

is the coordinate-wise MLE when fixing all other entries of θ (m) . Similar to (60), we obtain from the
optimality condition that
d`(m) (τ ;y)
dτ

2
(m)
∗
≤
θm
− θm

dτ 2

where θ̃m is some quantity between
Observing that
d`(m) (τ ; y)
dτ
and that

∗
θm

and

(74)

τ =θ̃m

(m)
θm .

(
X

=p
∗
τ =θm

∗
τ =θm

d2 `(m) (τ ;y)

i:i6=m

(m)

∗

eθ i
eθ i
∗
∗ −
(m)
θ
θ
∗
e i +e m
eθ i + eθ m

)
∗
+ λθm

(m)

∗

eθ i
1
1
eθi
(m)
−
=
≤ θi∗ − θi ,
∗
∗ −θ ∗ −
∗
(m)
θm
∗ −θ (m)
∗
i
θ
θ
θ
eθi + eθm
1
+
e
e i +e m
1+e m i
we arrive at
d`(m) (τ ; y)
dτ

≤p
∗
τ =θm


X

(m)

θi∗ − θi


i:i6=m
r
(ii)
np log n
,
.
L




(i)

∗
+ λθm
≤ pn



kθ ∗ − θ (m) k2
∗
√
+ λθm
n

where (i) follows from Cauchy-Schwarz and (ii) arises from Lemma 14. Additionally, `(m) (τ ; y) is strongly
convex in τ since
(m)
X
d2 `(m) (τ ; y)
eτ eθ i
λ
=p

2 + & np.
2
(m)
dτ
2
τ =θ̃m
i:i6=m eτ + eθi
These taken collectively with (74) give
s
(m)
∗
θm
− θm
.

C

log n
.
npL

Hoeffding’s and Bernstein’s inequalities

This section collects two standard concentration inequalities used throughout the paper, which can be easily
found in textbooks such as Boucheron et al. (2013). The proofs are omitted.
Lemma 17 (Hoeffding’s inequality). Let P
{Xi }1≤i≤n be a sequence of independent random variables where
n
Xi ∈ [ai , bi ] for each 1 ≤ i ≤ n, and Sn = i=1 Xi . Then
P(|Sn − E [Sn ] | ≥ t) ≤ 2e−2t

2

/

Pn

2
i=1 (bi −ai )

.

The next lemma is about a user-friendly version of the Bernstein inequality.
Lemma 18 (Bernstein’s inequality). Consider n independent random variables zl (1 ≤ l ≤ n), each satisfying |zl | ≤ B. For any a ≥ 2, one has
v
" n #
u
n
n
X
X
X
u
2a
E [zl2 ] + B log n
zl − E
zl ≤ t2a log n
3
l=1

l=1

l=1

with probability at least 1 − 2n−a .
33

D

Proof of Theorem 6

By Theorem 5, we obtain
kπ ∗ − πkπ∗ ≤

π ∗> (P ∗ − P )

(i)

π∗

1 − max {λ2 (P ∗ ), −λn (P ∗ )} − P ∗ − P

(ii)

(iii)

. kπ ∗> (P ∗ − P )k2 . √

. π ∗> (P ∗ − P )

π∗

π∗

1
kπ ∗ k2 ,
npL

where (i) is a consequence of Lemma 5, (ii) follows from the relationship between k · k2 and k · kπ∗ , and (iii)
follows as long as one can justify that
π ∗> (P − P ∗ )

2

.√

1
kπ ∗ k2 .
npL

(75)

Therefore, the rest of the proof is devoted to establishing (75). To simplify the notations hereafter, we denote
∆ := P − P ∗ . In fact, it is easy to check that for any i 6= j,
!
L
1 1 X (l)
∗
∗
yi,j − yi,j
(76)
∆i,j = Pi,j − Pi,j = 1{(i,j)∈E}
d L
l=1

and for 1 ≤ i ≤ n, one has



X

∗
∆i,i = Pi,i − Pi,i
= 1 −



Pi,j  − 1 −

j:j6=i

=−

X
j:j6=i

∗ 
Pi,j

j:j6=i

∗
Pi,j − Pi,j
=−




X

X

∆i,j

j:j6=i

= −∆lower
− ∆upper
.
i,i
i,i

(77)

where
∆lower
:=
i,i

X

∆i,j

and

∆upper
:=
i,i

j:j<i

X

∆i,j .

j:j>i

Towards proving (75), we decompose ∆ into four parts ∆ = ∆lower + ∆upper + ∆diag,l + ∆diag,u , where
∆lower is the lower triangular part (excluding the diagonal) of ∆, ∆upper is the upper triangular part, and




∆diag,l = −diag [∆lower
]1≤i≤n
and
∆diag,u = −diag [∆upper
]1≤i≤n .
i,i
i,i
The triangle inequality then gives
π ∗> ∆

2

≤ π ∗> ∆lower 2 + π ∗> ∆upper 2 + π ∗> ∆diag,l 2 + π ∗> ∆diag,u 2 .
{z
} |
{z
} |
{z
} |
{z
}
|
:=Ilower

:=Iupper

:=Idiag,l

:=Idiag,u

In what follows, we will focus on controlling the first term Ilower . The other three terms can be bounded
using nearly identical arguments.
Note that the j component of π ∗> ∆lower can be expressed as
X
 ∗>

π ∆lower j =
πi∗ ∆i,j .
i:i>j

Recall
∈ E, ∆i,j is a sum of L independent zero-mean random variables, and hence
 ∗> that for any pair (i, j)lower
π ∆lower j is a sum of Ldj
independent zero-mean random variables, where
dlower
:= {(i, j) | (i, j) ∈ E and i > j} .
j
34

In view of Hoeffding’s inequality (Lemma 17), one has, when conditional on G, that
!


2

2t
π ∗> ∆lower j ≥ t ≤ 2 exp − 1 lower
.
P
2
L (2 kπ ∗ k∞ )
(Ld)2 dj


Hence π ∗> ∆lower j can be treated as a sub-Gaussian random variable with variance proxy
σ2 

dlower
1
j
2
2
kπ ∗ k∞ .
kπ ∗ k∞ .
2
d L
dL

Given that the entries of π ∗> ∆lower are independent, we see that
2

(Ilower ) =

n
X
 ∗>
2
π ∆lower j ,
j=1

is a quadratic form of a sub-Gaussian vector. On the one hand, E[Ilower ]2 . nσ 2 . On the other hand, we
invoke (?, Theorem 1.1) to reach


 2
t
t
2
2
,
P{Ilower
− E[Ilower
] ≥ t} ≤ 2 exp −c min
nσ 4 σ 2
√
for some constant c > 0. By choosing t  σ 2 n log n, we see that with probability at least 1 − O(n−10 ),
2
2
Ilower
. E[Ilower
] + σ2

p
n
1
2
2
n log n . nσ 2 .
kπ ∗ k∞ .
kπ ∗ k2 .
dL
dL

The same upper bounds can be derived for other terms using the same arguments. We have thus established
(75) by recognizing that d & np.

References
Abbe, E., Fan, J., Wang, K. and Zhong, Y. (2017). Eigenvectors are sharp: Entry-wise perturbation
analysis. working paper .
Agarwal, A., Agarwal, S., Assadi, S. and Khanna, S. (2017). Learning with limited rounds of adaptivity: Coin tossing, multi-armed bandits, and ranking from pairwise comparisons. In Conference on
Learning Theory.
Agresti, A. (2013). Categorical data analysis. 3rd ed. Wiley Series in Probability and Statistics, WileyInterscience [John Wiley & Sons], Hoboken, NJ.
Ailon, N. (2012). An active learning algorithm for ranking from pairwise preferences with an almost optimal
query complexity. J. Mach. Learn. Res. 13 137–164.
Ammar, A. and Shah, D. (2011). Ranking: Compare, don’t score. In 2011 49th Annual Allerton Conference
on Communication, Control, and Computing (Allerton).
Baltrunas, L., Makcinskas, T. and Ricci, F. (2010). Group recommendations with rank aggregation
and collaborative filtering. In Proceedings of the Fourth ACM Conference on Recommender Systems.
RecSys ’10, ACM, New York, NY, USA.
Borda, J. d. (1784). Mémoire sur les élections au scrutin. Histoire de l’Academie Royale des Sciences pour
1781 (Paris, 1784) .
Borkar, V. S., Karamchandani, N. and Mirani, S. (2016). Randomized kaczmarz for rank aggregation
from pairwise comparisons. In 2016 IEEE Information Theory Workshop (ITW).

35

Boucheron, S., Lugosi, G. and Massart, P. (2013). Concentration inequalities. Oxford University
Press, Oxford. A nonasymptotic theory of independence, With a foreword by Michel Ledoux.
Bradley, R. A. and Terry, M. E. (1952). Rank analysis of incomplete block designs. I. The method of
paired comparisons. Biometrika 39 324–345.
Braverman, M. and Mossel, E. (2008). Noisy sorting without resampling. In Proceedings of the Nineteenth
Annual ACM-SIAM Symposium on Discrete Algorithms. ACM, New York.
Brémaud, P. (1999). Markov chains, vol. 31 of Texts in Applied Mathematics. Springer-Verlag, New York.
Gibbs fields, Monte Carlo simulation, and queues.
Bubeck, S. (2015). Convex optimization: Algorithms and complexity. Found. Trends Mach. Learn. 8
231–357.
Busa-Fekete, R., Szörényi, B., Weng, P., Cheng, W. and Hüllermeier, E. (2013). Top-k selection
based on adaptive sampling of noisy preferences. In Proceedings of the 30th International Conference on
International Conference on Machine Learning - Volume 28. ICML’13, JMLR.org.
Candès, E. J., Li, X. and Soltanolkotabi, M. (2015). Phase retrieval via Wirtinger flow: theory and
algorithms. IEEE Trans. Inform. Theory 61 1985–2007.
Chaudhuri, K., Chung, F. and Tsiatas, A. (2012). Spectral clustering of graphs with general degrees in
the extended planted partition model. In Conference on Learning Theory.
Chen, X., Bennett, P. N., Collins-Thompson, K. and Horvitz, E. (2013). Pairwise ranking aggregation in a crowdsourced setting. In Proceedings of the Sixth ACM International Conference on Web Search
and Data Mining. WSDM ’13, ACM, New York, NY, USA.
Chen, X., Gopi, S., Mao, J. and Schneider, J. (2017). Competitive analysis of the top-k ranking
problem. In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms.
SODA ’17, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.
Chen, Y. and Candes, E. (2016). The projected power method: An efficient algorithm for joint alignment
from pairwise differences. arXiv preprint arXiv:1609.05820 .
Chen, Y. and Candès, E. J. (2017). Solving random quadratic systems of equations is nearly as easy as
solving linear systems. Comm. Pure Appl. Math. 70 822–883.
Chen, Y. and Suh, C. (2015). Spectral mle: Top-k rank aggregation from pairwise comparisons. In
Proceedings of the 32Nd International Conference on International Conference on Machine Learning Volume 37. ICML’15, JMLR.org.
Chung, F. R. K. (1997). Spectral graph theory, vol. 92 of CBMS Regional Conference Series in Mathematics.
Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American
Mathematical Society, Providence, RI.
Davis, C. and Kahan, W. M. (1970). The rotation of eigenvectors by a perturbation. III. SIAM J. Numer.
Anal. 7 1–46.
Durrett, R. (2007). Random graph dynamics, vol. 20 of Cambridge Series in Statistical and Probabilistic
Mathematics. Cambridge University Press, Cambridge.
Dwork, C., Kumar, R., Naor, M. and Sivakumar, D. (2001). Rank aggregation methods for the web.
In Proceedings of the 10th International Conference on World Wide Web. WWW ’01, ACM, New York,
NY, USA.
El Karoui, N. (2017). On the impact of predictor geometry on the performance on high-dimensional
ridge-regularized generalized robust regression estimators. Probability Theory and Related Fields .

36

El Karoui, N., Bean, D., Bickel, P. J., Lim, C. and Yu, B. (2013). On robust regression with
high-dimensional predictors. Proceedings of the National Academy of Sciences 110 14557–14562.
Fan, J., Wang, W. and Zhong, Y. (2016). An `∞ eigenvector perturbation bound and its application to
robust covariance estimation. arXiv preprint arXiv:1603.03516 .
Ford, L. R., Jr. (1957). Solution of a ranking problem from binary comparisons. Amer. Math. Monthly
64 28–33.
Hajek, B., Oh, S. and Xu, J. (2014). Minimax-optimal inference from partial rankings. In Proceedings
of the 27th International Conference on Neural Information Processing Systems. NIPS’14, MIT Press,
Cambridge, MA, USA.
Heckel, R., Shah, N. B., Ramchandran, K. and Wainwright, M. J. (2016). Active ranking from
pairwise comparisons and when parametric assumptions don’t help. arXiv preprint arXiv:1606.08842 .
Hunter, D. R. (2004). MM algorithms for generalized Bradley-Terry models. Ann. Statist. 32 384–406.
Jain, P., Netrapalli, P. and Sanghavi, S. (2013). Low-rank matrix completion using alternating minimization (extended abstract). In STOC’13—Proceedings of the 2013 ACM Symposium on Theory of
Computing. ACM, New York, 665–674.
Jamieson, K. G. and Nowak, R. D. (2011). Active ranking using pairwise comparisons. In Proceedings of
the 24th International Conference on Neural Information Processing Systems. NIPS’11, Curran Associates
Inc., USA.
Jang, M., Kim, S., Suh, C. and Oh, S. (2016). Top-k ranking from pairwise comparisons: When spectral
ranking is optimal. arXiv preprint arXiv:1603.04153 .
Javanmard, A. and Montanari, A. (2015). De-biasing the lasso: Optimal sample size for gaussian designs.
arXiv preprint arXiv:1508.02757 .
Jiang, X., Lim, L.-H., Yao, Y. and Ye, Y. (2011). Statistical ranking and combinatorial Hodge theory.
Math. Program. 127 203–244.
Keshavan, R. H., Montanari, A. and Oh, S. (2010). Matrix completion from noisy entries. J. Mach.
Learn. Res. 11 2057–2078.
Lu, Y. and Negahban, S. N. (2014). Individualized rank aggregation using nuclear norm regularization.
arXiv preprint arXiv:1410.0860 .
Luce, R. D. (1959). Individual choice behavior: A theoretical analysis. John Wiley & Sons, Inc., New York;
Chapman & Hall, Ltd., London.
Masse, K. (1997). Statistical models applied to the rating of sports teams. Technical report Bluefield College
.
Mohajer, S. and Suh, C. (2016). Active top-k ranking from noisy comparisons. In 2016 54th Annual
Allerton Conference on Communication, Control, and Computing (Allerton).
Negahban, S., Oh, S. and Shah, D. (2017a). Rank centrality: ranking from pairwise comparisons. Oper.
Res. 65 266–287.
Negahban, S., Oh, S., Thekumparampil, K. K. and Xu, J. (2017b). Learning from comparisons and
choices. arXiv preprint arXiv:1704.07228 .
Oh, S. and Shah, D. (2014). Learning mixed multinomial logit model from ordinal data. In Proceedings
of the 27th International Conference on Neural Information Processing Systems. NIPS’14, MIT Press,
Cambridge, MA, USA.

37

Osting, B., Xiong, J., Xu, Q. and Yao, Y. (2016). Analysis of crowdsourced sampling strategies for
HodgeRank with sparse random graphs. Appl. Comput. Harmon. Anal. 41 540–560.
Page, L., Brin, S., Motwani, R. and Winograd, T. (1999). The pagerank citation ranking: Bringing
order to the web. Tech. rep., Stanford InfoLab.
Pananjady, A., Mao, C., Muthukumar, V., Wainwright, M. J. and Courtade, T. A. (2017).
Worst-case vs average-case design for estimation from fixed pairwise comparisons. arXiv preprint
arXiv:1707.06217 .
Park, D., Neeman, J., Zhang, J., Sanghavi, S. and Dhillon, I. S. (2015). Preference completion:
Large-scale collaborative ranking from pairwise comparisons. In Proceedings of the 32Nd International
Conference on International Conference on Machine Learning - Volume 37. ICML’15, JMLR.org.
Rajkumar, A. and Agarwal, S. (2014). A statistical convergence perspective of algorithms for rank
aggregation from pairwise data. In Proceedings of the 31st International Conference on International
Conference on Machine Learning - Volume 32. ICML’14, JMLR.org.
Rajkumar, A. and Agarwal, S. (2016). When can we rank well from comparisons of o(n log n) non-actively
chosen pairs? In 29th Annual Conference on Learning Theory.
Rohe, K., Chatterjee, S. and Yu, B. (2011). Spectral clustering and the high-dimensional stochastic
blockmodel. Ann. Statist. 39 1878–1915.
Shah, N. B., Balakrishnan, S., Bradley, J., Parekh, A., Ramchandran, K. and Wainwright, M.
(2014). When is it better to compare than to score? arXiv:1406.6618 .
Shah, N. B., Balakrishnan, S., Guntuboyina, A. and Wainwright, M. J. (2017). Stochastically
transitive models for pairwise comparisons: statistical and computational issues. IEEE Trans. Inform.
Theory 63 934–959.
Shah, N. B. and Wainwright, M. J. (2015). Simple, robust and optimal ranking from pairwise comparisons. arXiv preprint arXiv:1512.08949 .
Singer, A. (2011). Angular synchronization by eigenvectors and semidefinite programming. Appl. Comput.
Harmon. Anal. 30 20–36.
Soufiani, H. A., Chen, W. Z., Parkes, D. C. and Xia, L. (2013). Generalized method-of-moments for
rank aggregation. In Proceedings of the 26th International Conference on Neural Information Processing
Systems. NIPS’13, Curran Associates Inc., USA.
Soufiani, H. A., Parkes, D. C. and Xia, L. (2014). Computing parametric ranking models via rankbreaking. In Proceedings of the 31st International Conference on International Conference on Machine
Learning - Volume 32. ICML’14, JMLR.org.
Suh, C., Tan, V. Y. F. and Zhao, R. (2017). Adversarial top-K ranking. IEEE Trans. Inform. Theory
63 2201–2225.
Sur, P., Chen, Y. and Candès, E. J. (2017). The likelihood ratio test in high-dimensional logistic
regression is asymptotically a rescaled Chi-square. arXiv preprint arXiv:1706.01191 .
Tropp, J. A. (2012). User-friendly tail bounds for sums of random matrices. Found. Comput. Math. 12
389–434.
Tropp, J. A. (2015). An introduction to matrix concentration inequalities. Found. Trends Mach. Learn. 8
1–230.
Wu, R., Xu, J., Srikant, R., Massoulie, L., Lelarge, M. and Hajek, B. (2015). Clustering and
inference from pairwise comparisons. In Proceedings of the 2015 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems. SIGMETRICS ’15, ACM, New York, NY,
USA.
38

Xu, Q., Huang, Q., Jiang, T., Yan, B., Lin, W. and Yao, Y. (2012). Hodgerank on random graphs for
subjective video quality assessment. IEEE Transactions on Multimedia 14 844–857.
Zhong, Y. and Boumal, N. (2017). Near-optimal bounds for phase synchronization. arXiv preprint
arXiv:1703.06605 .

39

