# Random pairing MLE for estimation of item parameters
in Rasch model

Yuepeng Yang  Department of Statistics, University of Chicago  Cong Ma
Department of Statistics, University of Chicago

###### Abstract

The Rasch model, a classical model in the item response theory, is widely used
in psychometrics to model the relationship between individuals’ latent traits
and their binary responses on assessments or questionnaires. In this paper, we
introduce a new likelihood-based estimator—random pairing maximum likelihood
estimator ($\mathsf{RP\text{-}MLE}$) and its bootstrapped variant multiple
random pairing MLE ($\mathsf{MRP\text{-}MLE}$) that faithfully estimate the
item parameters in the Rasch model. The new estimators have several appealing
features compared to existing ones. First, both work for sparse observations,
an increasingly important scenario in the big data era. Second, both
estimators are provably minimax optimal in terms of finite sample
$\ell_{\infty}$ estimation error. Lastly, $\mathsf{RP\text{-}MLE}$ admits
precise distributional characterization that allows uncertainty quantification
on the item parameters, e.g., construction of confidence intervals of the item
parameters. The main idea underlying $\mathsf{RP\text{-}MLE}$ and
$\mathsf{MRP\text{-}MLE}$ is to randomly pair user-item responses to form
item-item comparisons. This is carefully designed to reduce the problem size
while retaining statistical independence. We also provide empirical evidence
of the efficacy of the two new estimators using both simulated and real data.

##  1 Introduction

The item response theory (IRT) [ER13] is a framework widely used in
psychometrics to model the relationship between individuals’ latent traits
(such as ability or personality) and their responses on assessments or
questionnaires. It is particularly useful in the development, analysis, and
scoring of tests and assessments; see a recent survey [CLLY21] for a
statistical account of IRT.

Among statistical models in IRT, the Rasch model [Ras60] is a simple but
fundamental one for modeling binary responses. Specifically, for a user $t$
(e.g., test-taker) and an item $i$ (e.g., test problem), Rasch model assumes
that the response user $t$ has to item $i$ is binary and obeys

$\mathbb{P}\left[\text{user $t$ ``beats'' item
$i$}\right]=\frac{e^{\zeta_{t}^{\star}}}{e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}}},$

where $\zeta_{t}^{\star},\theta_{i}^{\star}\in\mathbb{R}$ are latent traits of
user $t$ and item $i$, respectively. The term “beats” here refers to positive
responses, such as answering an exam question correctly, writing a favorable
review of a product, etc. We also call this response a comparison between user
$t$ and item $i$.

In this paper, we focus on estimating the item parameters
$\bm{\theta}^{\star}$, which is one of the four main statistical tasks
surrounding the Rasch model (or IRT more generally) listed in [CLLY21].
Estimating the item parameters is also quite useful in practice. For instance,
in education testing, $\bm{\theta}^{\star}$ could reveal the difficulty of the
exam questions, while in product reviews, $\bm{\theta}^{\star}$ could reveal
the popularity of the products. Various methods have been proposed for
estimating the item parameters $\bm{\theta}^{\star}$, including the joint
maximum likelihood estimator (JMLE), the marginal maximum likelihood
estimator, the conditional maximum likelihood estimator (CMLE), as well as the
spectral estimator recently proposed in [NZ22, NZ23]. We refer readers to a
recent article [Rob21] for comparisons between different item parameter
estimation methods. However, three main gaps remain in tackling item
estimation in the Rasch model:

  * •

Non-asymptotic guarantee. Apart from the recently proposed spectral estimator
[NZ22, NZ23], most theoretical guarantees for the likelihood-based estimators
are asymptotic. Since all the estimation procedures are necessarily applied
with finite samples, the asymptotic guarantee alone fails to inform
practitioners about the performance of different estimators when working with
a limited number of samples.

  * •

Sparse observations. It is not uncommon to encounter situations where each
user only responds to a handful of questions (or items). This brings the
challenge of incomplete and sparse observations. Among the likelihood based
estimators, CMLE is known to have numerical issues with incomplete
observations [Mol95]. While the spectral estimator [NZ22, NZ23] is capable of
handling incomplete observations, it still requires the observations to be
relatively dense. We will elaborate on this point later.

  * •

Uncertainty quantification. Beyond estimation, uncertainty quantification on
the item parameters is central to realizing the full potential of the Rasch
model. However, existing results do not address this problem under sparse
observations. An exception is the recent work [CLOX23], which is based on
joint estimation and inference on the item parameters $\bm{\theta}^{\star}$
and the user parameters $\bm{\zeta}^{\star}$. Their sampling scheme is
restrictive and requires a relatively dense sampling rate.

In light of the above gaps, we raise the following natural question:

_Can we develop an estimator for the item parameters_
$\bm{\theta}^{\star}$_that (1) enjoys optimal estimation guarantee in finite
sample, and (2) is amenable to tight uncertainty quantification, when the
observations are sparse?_

###  1.1 Main contributions

The main contribution of our work is the proposal of a novel estimator named
random pairing maximum likelihood estimator ($\mathsf{RP\text{-}MLE}$ in
short) that achieves the two desiderata listed above.

In essence, $\mathsf{RP\text{-}MLE}$ compiles user-item comparisons to item-
item comparisons by randomly pairing responses of the same user to different
items. This pairing procedure is carefully designed to extract information of
the item parameters while retaining statistical independence. After this
compilation step, item parameters $\bm{\theta}^{\star}$ are estimated by the
MLE $\widehat{\bm{\theta}}$ given the item-item comparisons.

Even when the observations are extremely sparse, $\mathsf{RP\text{-}MLE}$
achieves the following:

  * •

Regarding estimation, we show that both $\mathsf{RP\text{-}MLE}$ and its
bootstrapped version enjoy optimal finite sample $\ell_{\infty}$ error
guarantee. Compared to the conventional $\ell_{2}$ error guarantee, the
$\ell_{\infty}$ guarantee, as an entrywise guarantee is more fine-grained.
Consequently, we also show that $\mathsf{RP\text{-}MLE}$ can recover the
top-$K$ items with minimal sample complexity.

  * •

While the optimal $\ell_{\infty}$ error guarantee directly yields optimal
$\ell_{2}$ guarantee, such guarantee is only correct in an order-wise sense.
We provide a refined finite-sample $\ell_{2}$ error guarantee of
$\mathsf{RP\text{-}MLE}$ that is precise even in the leading constant.

  * •

Supplementing the estimation guarantee, we also build an inferential framework
based on $\mathsf{RP\text{-}MLE}$ $\widehat{\bm{\theta}}$. More specifically,
we precisely characterize the asymptotic distribution of
$\widehat{\bm{\theta}}$. This result facilitates several inferential tasks
such as hypothesis testing and construction of confidence regions of
$\bm{\theta}^{\star}$.

We test our methods on both synthetic and real data, which clearly show
competitive empirical estimation performance. The inferential result on
synthetic data also closely matches our theoretical predictions.

###  1.2 Prior art

##### Item response theory.

The item response theory is a popular statistical framework for modeling
response data. It often involves a probabilistic model that links categorical
responses to latent traits of both users and items. Early endeavors include
[Ras60] that introduces the Rasch model studied herein and [LNB68] that
describes a more general framework using parametric models. Popular IRT models
include the Rasch model, the two-parameter model (2PL), and the three-
parameter logistic model (3PL). As response data widely appears in real life,
IRT finds application in numerous fields including educational assessment
[DC10], psychometrics [LNB68], political science [VHSA20], and medical
assessment [FBC05]. See also [CLLY21] for an overview of IRT from
statisticians’ perspective.

##### Latent score estimation for Rasch model.

An important statistical question in the Rasch model is to estimate the latent
parameters of the items. As the Rasch model is an explicit probabilistic
model, many methods are based on the principle of maximum likelihood
estimation. For instance, marginal MLE assumes a prior on the user parameters
that is either given or optimized within a parametric distribution family. The
item parameter is then estimated by maximizing the marginal likelihood. A
drawback is that MMLE relies on a good prior. On the other hand, joint MLE
(JMLE) makes no distributional assumption and maximizes the joint likelihood
w.r.t. both the item and user parameters. However, it is not consistent for
estimating the item parameters when the number $m$ of items is fixed [Gho95].
Interested readers may also consult [Lin99] for an overview of other classical
estimators.

Several methods are more relevant to our proposed estimator
$\mathsf{RP\text{-}MLE}$ as they follow a similar philosophy to form item-item
pairs from user-item responses. Pseudo MLE (PMLE) [Zwi95] maximizes the sum of
the log-likelihood of all pairs of responses from the same users to different
items. However, due to the dependency issue, no satisfying finite sample
performance guarantee has been established. Another related approach is the
spectral method, in which a Markov chain on the space of items is formed and
the item parameters are estimated via the stationary distribution of the
Markov chain. The most recent works in this category are [NZ22, NZ23], which
essentially use the same idea as pseudo MLE in forming item-item comparisons.

##### The Bradley-Terry-Luce model with sparse comparisons.

An informed reader may realize that the Rasch model resembles the Bradley-
Terry-Luce (BTL) model [Luc59, BT52] in the ranking literature. Indeed, one
can view the Rasch model as a special case of the BTL model that distinguishes
the two groups of users and items, and only allows inter-group comparisons.
There has been a recent surge in interest in studying top-$K$ ranking in the
BTL model [CS15, JKSO16, CFMW19, CGZ22, GSZ23] and its extensions [FHY22,
FLWY22, FLWY23], especially under sparse observations of the pairwise
responses. Most notably, under a uniform sampling scheme, [CFMW19] shows that
(regularized) MLE and spectral methods are both optimal in top-$K$ ranking and
[GSZ23] provides inference results for both methods. Going beyond uniform
sampling, [Che23, LSR22] investigate the performance of MLE in the BTL model
with a general comparison graph and later [YCOM24] improves the analysis to
show the optimality of MLE for the BTL model in both uniform and semi-random
sampling.

##### Notation.

For a positive integer $n$, we denote $[n]=\\{1,2,\ldots,n\\}$. For any
$a,b\in\mathbb{R}$, $a\wedge b$ means the minimum of $a,b$ and $a\vee b$ means
the maximum of $a,b$. We use $\overset{\mathrm{d}}{\rightarrow}$ to denote
convergence in distribution. For a symmetric matrix
$\bm{A}\in\mathbb{R}^{n\times n}$, we use
$\lambda_{1}(\bm{A})\geq\lambda_{2}(\bm{A})\geq\cdots\geq\lambda_{n}(\bm{A})$
to denote its eigenvalues and $\bm{A}^{\dagger}$ to denote its Moore-Penrose
pseudo-inverse. For symmetric matrices $\bm{A},\bm{B}\in\mathbb{R}^{n\times
n}$, $\bm{A}\preceq\bm{B}$ means $\bm{B}-\bm{A}$ is positive semidefinite,
i.e., $\bm{v}^{\top}(\bm{B}-\bm{A})\bm{v}\geq 0$ for any
$\bm{v}\in\mathbb{R}^{n}$. We use $\bm{e}_{i}$ to denote the standard unit
vector with $1$ at $i$-th coordinate and 0 elsewhere. Unless specified
otherwise, $\log(\cdot)$ denotes the natural log.

##  2 Problem setup and new estimator $\mathsf{RP\text{-}MLE}$

In this section, we first introduce the formal setup of the item parameter
estimation problem in the Rasch model. Then we present the newly proposed
estimator $\mathsf{RP\text{-}MLE}$ along with the rationale behind its
development.

###  2.1 Problem setup

The Rasch model considers pairwise comparisons between elements of two groups:
users and items. Let $n$ (resp. $m$) be the number of users (resp. items).
Rasch assumes a user parameter $\bm{\zeta}^{\star}\in\mathbb{R}^{n}$ and an
item parameter $\bm{\theta}^{\star}\in\mathbb{R}^{m}$ that measures the latent
traits (e.g., difficulty of a problem) of users and items, respectively. For a
subset of possible user-item pairs $\mathcal{E}_{X}\subset[n]\times[m]$, we
observe binary responses $\\{X_{ti}\\}_{(t,i)\in\mathcal{E}_{X}}$ obeying

$\mathbb{P}[X_{ti}=1]=\frac{e^{\theta_{i}^{\star}}}{e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}}}.$
(1)

Here $X_{ti}=1$ means user $t$ has negative response against item $i$ (e.g.,
unable to solve a problem). The goal is to estimate $\bm{\theta}^{\star}$, the
item parameters.

To model sparse observations/comparisons, we assume that
$\mathbb{P}[(t,i)\text{ is compared}]=p$ independently for every
$(t,i)\in[n]\times[m]$. To put it in the language of graph theory, we denote
the associated bipartite comparison graph to be
$\mathcal{G}_{X}=(\mathcal{V}_{X},\mathcal{E}_{X})$, where $\mathcal{V}_{X}$
consists of $n$ users and $m$ items. Then essentially, we are assuming that
the bipartite graph follows an Erdős-Rényi random model.111Alternatively, we
can assume each user responds to $mp$ items uniformly at random. Our estimator
and performance guarantee continue to work in this sampling scheme.

Before moving on, we introduce the condition numbers that characterize the
range of the latent traits. Let $\kappa_{1}$, $\kappa_{2}$, and $\kappa$ be
defined by
$\log\;(\kappa_{1})=\max_{ij}\\{|\theta_{i}^{\star}-\theta_{j}^{\star}|\\}$,
$\log(\kappa_{2})=\max_{ti}\\{|\zeta_{t}^{\star}-\theta_{i}^{\star}|\\}$, and
$\kappa=\max\\{\kappa_{1},\kappa_{2}\\}$, respectively.

###  2.2 Random pairing maximum likelihood estimator

Algorithm 1 Random Pairing Maximum Likelihood Estimator
($\mathsf{RP\text{-}MLE}$)

  1. 1.

For each tester $t$,

    1. (a)

Randomly split the $n_{t}$ problems taken by tester $t$ into $\lfloor
n_{t}/2\rfloor$ pairs of problems.

    2. (b)

For each $(i,j)\in[m]\times[m]$, do the following:

      1. i.

If $(i,j)$ is selected as a pair in Step 1(a), $R_{ij}^{t}=1$. Furthermore, if
$X_{ti}\neq X_{tj}$, let $Y_{ij}^{t}=\mathds{1}\\{X_{ti}<X_{tj}\\}$ and
$L_{ij}^{t}=1$; if $X_{ti}=X_{tj}$, let $L_{ij}^{t}=0$.

      2. ii.

If $(i,j)$ is not selected as a pair in Step 1(a), let $L_{ij}^{t}=0$ and
$R_{ij}^{t}=0$.

  2. 2.

Let $\mathcal{E}_{Y}$ be a set of edges defined by
$\mathcal{E}_{Y}\coloneqq\\{(i,j):\sum_{t=1}L_{ij}^{t}\geq 1\\}$ and let
$\mathcal{G}_{Y}=([m],\mathcal{E}_{Y})$. For each $(i,j)\in\mathcal{E}_{Y}$,
let $L_{ij}\coloneqq\sum_{t=1}^{n}L_{ij}^{t}$ and
$Y_{ij}\coloneqq(1/L_{ij})\sum_{\\{t:L_{ij}^{t}=1\\}}Y_{ij}^{t}$.

  3. 3.

Compute MLE on $Y_{ij}$, i.e.,

$\widehat{\bm{\theta}}\coloneqq\arg\min_{\bm{\theta}:\bm{1}_{m}^{\top}\bm{\theta}=0}\mathcal{L}(\bm{\theta})$

where

$\mathcal{L}(\bm{\theta})=\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}\left(-Y_{ji}(\theta_{i}-\theta_{j})+\log(1+e^{\theta_{i}-\theta_{j}})\right).$
(2)

  4. 4.

Return the top-$K$ items by selecting the top-$K$ entries of
$\widehat{\bm{\theta}}$.

In this section, we present our main method $\mathsf{RP\text{-}MLE}$; see
Algorithm 1. The algorithm can be divided into two parts. The first part—Steps
1 and 2—uses random pairing to compile the observed user-item responses
$\bm{X}\in\mathbb{R}^{n\times m}$ to item-item comparisons
$\bm{Y}\in\mathbb{R}^{m\times m}$. The second part—Steps 3 and 4—computes a
standard MLE on the item-item comparisons. Some intuitions regarding the
development of $\mathsf{RP\text{-}MLE}$ are in order.

##### Random pairing to construct item-item comparisons.

The idea of pairing is that by matching the responses $X_{ti}$ with $X_{tj}$,
we form a comparison between items $i$ and $j$ to directly extract information
of item parameters $\theta_{i}^{\star}$ and $\theta_{j}^{\star}$. More
specifically, the item-item comparisons $\bm{Y}$ follow the Bradley-Terry-Luce
model [BT52, Luc59],i.e.,
$\mathbb{P}[Y_{ij}^{t}]=e^{\theta_{j}^{\star}}/(e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}})$;
see Section 4.1 for a formal argument.

By compiling user-item responses to item-item comparisons, we reduce the size
of the data matrix from $n\times m$ to $m\times m$, and also the number of
intrinsic parameters from $n+m$ to $m$, since the likelihood function (2) of
$Y_{ij}^{t}$ is completely independent with the user parameter
$\zeta_{t}^{\star}$.

More importantly, the pairing is performed in a disjoint fashion. This ensures
that all constructed item-item comparisons $Y_{ij}^{t}$ are independent with
each other; see Section 4.1 for a formal statement. This is the key ingredient
that enables us to improve over previous implementation of item-item
comparisons, such as pseudo-likelihood [Cho82, Zwi95] and spectral methods
[NZ22, NZ23].

##### A variant via bootstrapping.

A drawback of this random pairing is that it potentially induces a loss of
information since not every possible pairing is considered. Once $X_{ti}$ is
paired with $X_{tj}$, we cannot pair $X_{ti}$ with another response $X_{tl}$.
That being said, we will later show that the $\ell_{\infty}$ error of
$\mathsf{RP\text{-}MLE}$ is still rate-optimal up to logarithmic factors.
Hence the loss of information can at most incur a small constant factor in
terms of estimation error. Nevertheless, we provide a remedy to this
phenomenon in $\mathsf{MRP\text{-}MLE}$ (Algorithm 2) by running (in other
words, bootstrapping) the $\mathsf{RP\text{-}MLE}$ multiple times with
different random data splitting and averaging the resulting estimates.
$\mathsf{MRP\text{-}MLE}$ trivially enjoys the same theoretical guarantee (cf.
Theorem 1) while improving the estimation error in practice over
$\mathsf{RP\text{-}MLE}$. See Figure 4 in Section 5.1 for the empirical
evidence.

Algorithm 2 Multiple Random Pairing Maximum Likelihood Estimator
($\mathsf{MRP\text{-}MLE}$)

  1. 1.

Let $n_{\mathrm{split}}$ be the number of runs. For
$i=1,\ldots,n_{\mathrm{split}}$, run $\mathsf{RP\text{-}MLE}$ (Algorithm 1),
each time with an independent random splitting in Step 1. Let the estimated
latent scores be $\widehat{\bm{\theta}}^{(i)}$.

  2. 2.

Estimate the latent score with

$\widehat{\bm{\theta}}=\frac{1}{n_{\mathrm{split}}}\sum_{i=1}^{n_{\mathrm{split}}}\widehat{\bm{\theta}}^{(i)}.$

  3. 3.

Estimate the top-$K$ items by selecting the top-$K$ entries of
$\widehat{\bm{\theta}}$.

##  3 Main results

In this section, we collect the main theoretical guarantees for
$\mathsf{RP\text{-}MLE}$ and its variant $\mathsf{MRP\text{-}MLE}$. Section
3.1 focuses on the finite sample $\ell_{\infty}$ error bound. While one can
translate $\ell_{\infty}$ error into an $\ell_{2}$ bound, in Section 3.2, we
present a much sharper characterization of the $\ell_{2}$ error of
$\mathsf{RP\text{-}MLE}$. Lastly in Section 3.3, we provide a distributional
characterization of $\mathsf{RP\text{-}MLE}$.

###  3.1 $\ell_{\infty}$ estimation error and top-$K$ recovery

Without loss of generality, we assume that the scores of the items are
ordered, i.e.,
$\theta_{1}^{\star}\geq\theta_{2}^{\star}\geq\cdots\geq\theta_{m}^{\star}$,
and denote $\Delta_{K}\coloneqq\theta_{K}^{\star}-\theta_{K+1}^{\star}$. In
words, $\Delta_{K}$ measures the difference between the difficulty levels of
items $K$ and $K+1$. The following theorem provides $\ell_{\infty}$ error
bounds and top-$K$ recovery guarantee for both $\mathsf{RP\text{-}MLE}$ and
$\mathsf{MRP\text{-}MLE}$. We defer its proof to Section 4.2.

######  Theorem 1.

Suppose that $mp\geq 2$ and $np\geq
C_{1}\kappa_{1}^{4}\kappa_{2}^{5}\log^{3}(n)$ for some sufficiently large
constant $C_{1}>0$. Suppose that there exists some constant $\alpha>0$ such
that $m\leq n^{\alpha}$. Let $\widehat{\bm{\theta}}$ be the
$\mathsf{RP\text{-}MLE}$ estimator. With probability at least $1-O(n^{-10})$,
$\widehat{\bm{\theta}}$ satisfies

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}\leq
C_{2}\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{\log(n)}{np}}.$

Consequently, the estimator is able to exactly recover the top-$K$ items as
soon as

$np\geq\frac{C_{3}\kappa_{1}^{2}\kappa_{2}^{1}\log(n)}{\Delta_{K}^{2}}.$

Here $C_{2},C_{3}>0$ are some universal constants. All the claims continue to
hold for $\mathsf{MRP\text{-}MLE}$ as long as there exists some constant
$\beta>0$ such that $n_{\mathrm{split}}\leq n^{\beta}$.

Some remarks are in order.

##### Finite sample minimax optimality.

Intuitively, as we have $\Theta(mnp)$ sampled comparisons and $m$ parameters
to estimate, the optimal $\ell_{\infty}$ error should be
$O(\sqrt{m/(mnp)})=O(1/\sqrt{np})$. This is formalized in the following result
from [NZ23].

######  Proposition 1 (Minimax lower bound, Theorems 3.3 and 3.4 in [NZ23]).

Assume that $np\geq C_{1}$ for some suffiently large constant $C_{1}>0$. For
any $n$ and $m$, there exists a class of user and item parameters $\Theta$
such that

$\inf_{\hat{\bm{\theta}}}\sup_{(\bm{\zeta}^{\star},\bm{\theta}^{\star})\in\Theta}\mathbb{E}\left\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\right\|_{2}^{2}\geq\frac{C_{2}m}{np},\qquad\text{and}\qquad\inf_{\hat{\bm{\theta}}}\sup_{(\bm{\zeta}^{\star},\bm{\theta}^{\star})\in\Theta}\mathbb{E}\left\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\right\|_{\infty}^{2}\geq\frac{C_{2}}{np},$

where $C_{2}>0$ is some constant. Moreover if $np\leq
C_{K}\log(m)/\Delta_{K}^{2}$ for some constant $C_{K}>0$, we have

$\inf_{\hat{\bm{\theta}}}\sup_{(\bm{\zeta}^{\star},\bm{\theta}^{\star})\in\Theta}\mathbb{P}\left[\widehat{\bm{\theta}}\text{
fails to identify all top-}K\text{ items}\right]\geq\frac{1}{2}.$

Comparing our upper bounds with the lower bound in the proposition, we can see
that both $\mathsf{RP\text{-}MLE}$ and $\mathsf{MRP\text{-}MLE}$ are rate-
optimal in $\ell_{\infty}$ estimation error and top-$K$ recovery sample
complexity, up to logarithmic and $\kappa$ factors.

##### Sample size requirement.

While the rates are optimal, it is worth noting that in Theorem 1 we have made
several sample size requirements. We now elaborate on them.

First, the assumption $mp\geq 2$ is a mild requirement on the expected number
of items compared by each user. This is required as we need user $t$ to
compare at least two items to form a comparison between items. In fact, if a
user only responds to one item, it is clear that this data point is not useful
at all for item parameter estimation.

Second, it is a standard and necessary requirement to have $np\gtrsim\log(n)$
to make sure that each item is compared to at least one user with high
probability. In Theorem 1 we require an extra $\log^{2}(n)$ factor to suppress
a quadratic error term that comes up in the analysis. This cubic log factor
can possibly be loose, but it is a minor issue and we leave it as a problem
for future research.

Lastly, $m\leq n^{\alpha}$ and $n_{\mathrm{split}}\leq n^{\beta}$ are both
minor as we only need these to allow union bounds over $m$ and
$n_{\mathrm{split}}$.

##### Comparison with [NZ23].

The closest result to our paper in terms of $\ell_{\infty}$ guarantee for the
Rasch model appears in the recent work [NZ23]. Their spectral method uses a
similar construction of the item-item comparisons but without disjoint
pairing. To provide detailed comparisons, we restate their results below.

######  Proposition 2 (Informal, Theorem 3.1 in [NZ23]).

Assume that $p\gtrsim\log(m)/\sqrt{n}$ and $mp\gtrsim\log(n)$, with
probability at least $1-O(m^{-10}+n^{-10})$, spectral estimator
$\widehat{\bm{\theta}}_{\mathrm{spectral}}$ satisfies

$\|\widehat{\bm{\theta}}_{\mathrm{spectral}}-\widehat{\bm{\theta}}\|_{\infty}\lesssim\kappa^{9}\sqrt{\frac{\log(m)}{np}}.$

This error rate is similar to ours. However, the required sample size is much
larger as they require $p\gtrsim\log(m)/\sqrt{n}$. Our result makes a
significant improvement by allowing a much smaller sampling rate $p$, cf.
$mp\geq 2$ and $np\gtrsim\log^{3}(n)$. In fact, as we have argued earlier, it
is nearly the sparsest possible. In addition, our methods enjoy a
significantly better error rate dependency on $\kappa$. In Section 5.1, we
provide empirical evidence for this improvement: when $\kappa$ is large,
$\mathsf{RP\text{-}MLE}$ and $\mathsf{MRP\text{-}MLE}$ outperform the spectral
methods in [NZ23].

###  3.2 Refined $\ell_{2}$ error characterization

The $\ell_{\infty}$ error guarantee in Theorem 1 immediately implies an
$\ell_{2}$ error bound

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|\leq
C\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{m\log(n)}{np}},$ (3)

which is rate-optimal compared to the minimax lower bound in Proposition 1.
However, this guarantee is only correct in an order-wise sense. In this
section, we present a refined characterization of
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ that is precise in the leading
constant.

We start with necessary notation. Let $z_{ij}\coloneqq
e^{\theta_{i}^{\star}}e^{\theta_{j}^{\star}}/(e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}})^{2}$
and $\widehat{z}_{ij}\coloneqq
e^{\hat{\theta}_{i}}e^{\hat{\theta}_{j}}/(e^{\hat{\theta}_{i}}+e^{\hat{\theta}_{j}})^{2}$.
Define a matrix

$\bm{L}_{L\tilde{z}}=\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}\widetilde{z}_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}$
(4)

for $\tilde{z}\in\\{z,\widehat{z}\\}$ and let $\bm{L}_{L\tilde{z}}^{\dagger}$
be its pseudo-inverse. One can view $\bm{L}_{L\tilde{z}}$ as a weighted graph
Laplacian. In the following theorem, we show that
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ concentrates tightly around
$[\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})]^{1/2}$. We defer the analysis
to Section 4.3 and the complete proof to Section D.

######  Theorem 2.

Instate the assumptions of Theorem 1. Then for some constants $C_{1},C_{2}>0$,
with probability at least $1-O(n^{-10})$, we have

$\displaystyle\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\right|$
$\displaystyle\leq
C_{1}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{\log(n)}{np}}+\frac{C_{2}\kappa_{1}^{6}\sqrt{m}\log^{2}(n)}{np},$
(5)
$\displaystyle\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{L\widehat{z}}^{\dagger})}\right|$
$\displaystyle\leq
C_{1}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{\log(n)}{np}}+\frac{C_{2}(\kappa_{1}^{6}+\kappa_{1}^{7/2}\kappa_{2}^{2})\sqrt{m}\log^{2}(n)}{np}.$
(6)

Theorem 2 is more refined compared to (3). First, it provides both upper and
lower bounds for $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$. Second,
there is no hidden constant in front of the leading term
$[\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})]^{1/2}$. In addition, inspecting the
proof of Theorem 2, we see that

$\sqrt{\frac{m-1}{np}}\leq\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\leq
4\kappa_{1}^{1/2}\kappa_{2}^{1/2}\sqrt{\frac{m}{np}},$

and the same holds for
$[\mathrm{Trace}(\bm{L}_{L\widehat{z}}^{\dagger})]^{1/2}$. Consequently, the
right hand sides of both (5) and (6) are lower order terms compared to
$[\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})]^{1/2}$ when $n,m\rightarrow\infty$.
Indeed this recovers the naive $\ell_{2}$ bound (3) under appropriate sample
size assumptions.

###  3.3 Uncertainty quantification

Now we move on to uncertainty quantification about the
$\mathsf{RP\text{-}MLE}$ estimator $\widehat{\bm{\theta}}$ in order to
perform, say, hypothesis testing on the item parameters $\bm{\theta}^{\star}$.
Our starting point is the following characterization of the limiting
distribution of $\widehat{\bm{\theta}}-\bm{\theta}^{\star}$. The analysis is
deferred to Section 4.3 and the complete proof is deferred to Section E.

######  Theorem 3.

Instate the assumptions of Theorem 1 and fix $m$. As
$np/\log^{4}(n)\rightarrow\infty$, for $\tilde{z}\in\\{z,\widehat{z}\\}$, we
have

$\bm{L}_{L\tilde{z}}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}\right).$

Though simple, we record the following immediate consequence on the asymptotic
normality of $\widehat{\theta}_{i}-\theta_{i}^{\star}$ for any single
coordinate $i$.

######  Corollary 1.

Instate the assumptions of Theorem 1. As $np/\log^{4}(n)\rightarrow\infty$,
for $\tilde{z}\in\\{z,\widehat{z}\\}$, we have

$\frac{\widehat{\theta}_{i}-\theta_{i}^{\star}}{\sqrt{[\bm{L}_{L\tilde{z}}^{\dagger}]_{ii}}}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,1\right).$
(7)

We compare our result with that in [CLOX23]. Theorem 8 therein studies the
Rasch model and provides an inference result for a joint estimator of
$\bm{\theta}^{\star}$ and $\bm{\zeta}^{\star}$. However, it requires a dense
sampling scheme with $p\gtrsim 1/\sqrt{\log(n)/n}$ and $p\gtrsim
n\log^{2}(n)/m^{2}$. It also requires that both $m$ and $n$ tend to infinity.
These two assumptions are more restrictive than ours.

Theorem 3 is quite general as it allows various applications to inference. As
an illustrative example, we demonstrate how to construct confidence intervals
for a particular item parameter $\theta_{i}^{\star}$.

Fix any $\alpha>0$, the level of significance. Define

$\displaystyle\mathcal{C}_{i}^{+}(\alpha)$
$\displaystyle\coloneqq\widehat{\theta}_{i}+\Phi^{-1}(1-\alpha)\left[\bm{L}_{L\hat{z}}^{\dagger}\right]_{ii}^{1/2},$
$\displaystyle\mathcal{C}_{i}^{-}(\alpha)$
$\displaystyle\coloneqq\widehat{\theta}_{i}-\Phi^{-1}(1-\alpha)\left[\bm{L}_{L\hat{z}}^{\dagger}\right]_{ii}^{1/2},$

where $\Phi$ is the cdf for the standard normal distribution. Corollary 1
implies the following proposition.

######  Proposition 3.

Suppose the statement in Theorem 3 holds. Then for any $(i,j)$, we have

$\displaystyle\lim_{np\rightarrow\infty}\mathbb{P}\left[\theta_{i}^{\star}\geq\mathcal{C}_{i}^{+}(\alpha)\right]$
$\displaystyle=\alpha,$ (8a)
$\displaystyle\lim_{np\rightarrow\infty}\mathbb{P}\left[\theta_{i}^{\star}\leq\mathcal{C}_{i}^{-}(\alpha)\right]$
$\displaystyle=\alpha,$ (8b)

and

$\lim_{np\rightarrow\infty}\mathbb{P}\left[\mathcal{C}_{i}^{-}(\alpha/2)\leq\theta_{i}^{\star}\leq\mathcal{C}_{i}^{+}(\alpha/2)\right]=\alpha.$
(9)

This proposition can be extended to confidence intervals for any linear
combination of $\bm{\theta}^{\star}$, such as
$\theta_{i}^{\star}-\theta_{j}^{\star}$. These confidence intervals naturally
enable hypothesis testings such as

$\mathcal{H}_{0}:\theta_{j}^{\star}\geq\theta_{i}^{\star},\qquad\mathcal{H}_{1}:\theta_{j}^{\star}<\theta_{i}^{\star}$

and

$\mathcal{H}_{0}:\text{item $i$ is a top-$K$
item},\qquad\mathcal{H}_{1}:\text{item $i$ is not a top-$K$ item}.$

Various other inference tasks have been considered in the ranking literature
[LFL23][FLWY22, FLWY23]. We expect our general distributional characterization
in Theorem 3 to be readily applied for these tasks.

##  4 Analysis

In this section, we present the main steps in the analysis to obtain
theoretical results in the previous section. Section 4.1 gives a complete
argument on the reduction to the BTL model we alluded in Section 2.2, Section
4.2 provides the analysis of the $\ell_{\infty}$ error, and Section 4.3
provides the analysis of both the $\ell_{2}$ error and the asymptotic
distribution as they share similar components in the analysis.

###  4.1 Reduction to Bradley-Terry-Luce model

A key idea in $\mathsf{RP\text{-}MLE}$ is the randomized pairing in Steps 1
and 2 of Algorithm 1. It compiles the user-item responses $\bm{X}$ to item-
item comparisons $\bm{Y}$. In this section, we make a detailed argument that
$\bm{Y}$ follows the Bradley-Terry-Luce model with a non-uniform sampling
scheme.

Recall that $L_{ij}^{t}\coloneqq\mathds{1}\\{X_{ti}\neq X_{tj}\\}$ and
$Y_{ij}^{t}\coloneqq\mathds{1}\\{X_{ti}<X_{tj}\\}$. The following fact
provides the distribution of $Y_{ij}^{t}$ conditional on $L_{ij}^{t}=1$. We
defer its proof to Section A.

######  Fact 1.

Let $i,j$ be two items and $t$ be a user. Suppose that user $t$ has responded
to both items $i$ and $j$. Let $X_{ti}$ and $X_{tj}$ be the responses sampled
from the probability model (1). Then we have

$\mathbb{P}[X_{ti}<X_{tj}\mid
L_{ij}^{t}=1]=\frac{e^{\theta_{j}^{\star}}}{e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}}},$

and

$\mathbb{P}[L_{ij}^{t}=1]\geq\frac{2\kappa_{2}}{(1+\kappa_{2})^{2}}.$ (10)

Fact 1 shows that conditional on $L_{ij}^{t}=1$, $Y_{ij}^{t}$ follows the BTL
model with parameters $\bm{\theta}^{\star}$. More importantly, as we deploy
random pairing (cf. Step 1a), each response $X_{ti}$ is used at most once. As
a result, conditional on $\\{L_{ij}^{t}\\}_{ijt}$, $Y_{ij}^{t}$’s are jointly
independent across users and items. In light of these, we can equivalently
describe the data generating process of $\bm{Y}$ as follows:

  1. 1.

For each user-item pair $(t,i)$, there is a comparison between them with
probability $p$ independently.

  2. 2.

Randomly split the $n_{t}$ problems taken by tester $t$ into $\lfloor
n_{t}/2\rfloor$ pairs of problems. (Step 1(a) of Algorithm 1)

  3. 3.

For all $(t,i,j)$, items $i$ and $j$ are compared by user $t$ if
$L_{ij}^{t}\coloneqq\mathds{1}\\{X_{ti}\neq X_{tj}\\}=1$.

  4. 4.

Conditioned on $L_{ij}^{t}=1$, one observes the outcome
$Y_{ij}^{t}\coloneqq\mathds{1}\\{X_{ti}<X_{tj}\\}$.

Steps 1–3 generates a non-uniform comparison graph $\mathcal{E}_{Y}$ between
items. Step 4 reveals the independent outcomes of these comparisons following
the BTL model, conditional on the graph $\mathcal{E}_{Y}$. This justifies that
(2) is truly the likelihood function of the BTL model conditional on the
comparison graph $\mathcal{E}_{Y}$.

###  4.2 Analysis of $\ell_{\infty}$ error

We have seen that analyzing $\mathsf{RP\text{-}MLE}$ under the Rasch model can
be reduced to analyzing the MLE under the BTL model. This reduction allows us
to invoke the general theory of MLE in the BTL model established in the recent
work [YCOM24].

To facilitate the presentation, we introduce the necessary notation. For any
$i\in[m]$, let $d_{i}\coloneqq\sum_{j:j\neq i}L_{ij}$ be the weighted degree
of item $i$ in $\mathcal{G}_{Y}$ and $d_{\max}=\max_{i\in[m]}d_{i}$. Let the
weighted graph Laplacian $\bm{L}_{L}$ be

$\bm{L}_{L}\coloneqq\sum_{i,j:i>j}L_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}.$

The following lemma adapts Theorem 3 of the recent work [YCOM24] to our
setting.

######  Lemma 1 (Theorem 3 in [YCOM24]).

Assume that $\mathcal{G}_{Y}$ is connected, and that

$[\lambda_{m-1}(\bm{L}_{L})]^{5}\geq
C_{1}\kappa_{1}^{4}(d_{\max})^{4}\log^{2}(n)$ (11)

for some large enough constant $C_{1}>0$. Then with probability at least
$1-n^{-10}$, we have

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}\leq
C_{2}\kappa_{1}\sqrt{\frac{\log(n)}{\lambda_{m-1}(\bm{L}_{L})}}$

for some constant $C_{2}>0$.

To leverage this general result, we need to characterize the spectral and
degree properties of the comparison graph $\mathcal{G}_{Y}$, which is achieved
in the following two lemmas. The proofs are deferred to Section B.

######  Lemma 2 (Degree bound in $\mathcal{G}_{Y}$).

Suppose that $np\geq C\kappa_{2}^{2}\log(n)$ for some large enough constant
$C>0$ and $m\leq n^{\alpha}$ for some sufficiently large constant $\alpha>0$.
With probability at least $1-2n^{-10}$, for all $i\in[m]$,

$\frac{1}{24\kappa_{2}}np\leq d_{i}\leq\frac{3}{2}np.$ (12)

######  Lemma 3.

Suppose $mp\geq 2$, $np\geq C\kappa_{2}^{2}\log(n)$ for some large enough
constant $C$, and $m\leq n^{\alpha}$ for some constant $\alpha>0$. With
probability at least $1-10n^{-10}$, we have

$\frac{np}{4\kappa_{2}}\leq\lambda_{m-1}(\bm{L}_{L})\leq\lambda_{1}(\bm{L}_{L})\leq
3np$ (13)

and

$\frac{np}{16\kappa_{1}\kappa_{2}}\leq\lambda_{m-1}(\bm{L}_{Lz})\leq\lambda_{1}(\bm{L}_{Lz})\leq
np.$ (14)

####  4.2.1 Proof of Theorem 1

Now we are ready to prove Theorem 1. We focus on analyzing
$\mathsf{RP\text{-}MLE}$, as the analysis of $\mathsf{MRP\text{-}MLE}$ follows
immediately from the union bound of the different data splitting and the
triangular inequality:

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}\leq\frac{1}{n_{\mathrm{split}}}\sum_{i=1}^{n_{\mathrm{split}}}\|\widehat{\bm{\theta}}^{(i)}-\bm{\theta}^{\star}\|_{\infty}.$

By assumption we have $mp\geq 2$ and $np\geq
C_{1}\kappa_{1}^{4}\kappa_{2}^{5}\log^{3}(n)$ for some constant $C_{1}>0$.
Then we can apply Lemmas 3 and 2 to see that

$\frac{np}{4\kappa_{2}}\leq\lambda_{m-1}(\bm{L}_{L}),\qquad\text{and}\qquad
d_{\max}\leq\frac{3}{2}np.$

We observe that (11) is satisfied as long as $np\geq
C_{1}\kappa_{1}^{4}\kappa_{2}^{5}\log^{3}(n)$ for some constant $C_{1}$ that
is large enough. Invoking Lemma 1, we conclude that

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}\leq
C_{2}\kappa_{1}\sqrt{\frac{\log(n)}{\lambda_{m-1}(\bm{L}_{L})}}\leq
2C_{2}\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{\log(n)}{np}}.$

It remains to show the top-$K$ recovery sample complexity. As
$\theta_{1}^{\star}\geq\cdots\geq\theta_{K}^{\star}>\theta_{K+1}^{\star}\geq\cdots\geq\theta_{m}^{\star}$
by assumption, it suffices to show
$\widehat{\theta}_{i}-\widehat{\theta}_{j}>0$ for any $i\leq K$ and $j>K$.
Using the $\ell_{\infty}$ error bound, we have that

$\displaystyle\widehat{\theta}_{i}-\widehat{\theta}_{j}$
$\displaystyle\geq\left(\theta_{i}^{\star}-\theta_{j}^{\star}\right)-\left|\widehat{\theta}_{i}-\theta_{i}^{\star}\right|-\left|\widehat{\theta}_{j}-\theta_{j}^{\star}\right|\geq\Delta_{K}-4C_{2}\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{\log(n)}{np}}.$

Then $\widehat{\theta}_{i}-\widehat{\theta}_{j}>0$ as long as

$np\geq\frac{16C_{2}^{2}\kappa_{1}^{2}\kappa_{2}\log(n)}{\Delta_{K}^{2}}.$

###  4.3 Analysis for refined $\ell_{2}$ error and uncertainty quantification

To make the main text concise, we sketch the key ideas that lead to both
Theorems 2 and 3.

We start with some necessary notation. For convenience, for any $(i,j)$,
renumber the set of comparisons $\\{Y_{ji}^{t}\\}_{t:L_{ij}^{t}=1}$ as
$\\{Y_{ji}^{(l)}\\}_{l=1,\ldots,L_{ij}^{t}}$ in an arbitrary order. Let
$\epsilon_{ij}^{(l)}\coloneqq
Y_{ji}^{(l)}-\sigma(\theta_{i}^{\star}-\theta_{j}^{\star})$ where
$\sigma(x)=e^{x}/(1+e^{x})$ is the sigmoid function. Let
$\bm{B}\coloneqq\left[\cdots,\sqrt{z_{ij}}(\bm{e}_{i}-\bm{e}_{j}),\cdots\right]_{i>j:(i,j)\in\mathcal{E}_{Y}}\in\mathbb{R}^{m\times
L_{\mathrm{total}}}$ (repeat $L_{ij}$ times for edge $(i,j)$), where
$L_{\mathrm{total}}\coloneqq\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}L_{ij}$ is the
total number of observed effective comparisons in $\mathcal{G}_{Y}$. Define

$\widehat{\bm{\epsilon}}\coloneqq\left[\cdots,\epsilon_{ij}^{(l)}/\sqrt{z_{ij}},\cdots\right]_{i>j:(i,j)\in\mathcal{E}_{Y},l=1,\ldots,L_{ij}^{t}}\in\mathbb{R}^{L_{\mathrm{total}}}.$

The key step in establishing Theorems 2 and 3 is the following refined
characterization of the estimation error
$\widehat{\bm{\theta}}-\bm{\theta}^{\star}$.

######  Lemma 4.

Instate the assumptions of Theorem 1. With probability at least
$1-O(n^{-10})$, the estimator $\widehat{\bm{\theta}}$ given by the Algorithm
(1) can be written as

$\widehat{\bm{\theta}}-\bm{\theta}^{\star}=-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\bm{r},$
(15)

where $\bm{r}\in\mathbb{R}^{m}$ is a random vector obeying
$\|\bm{r}\|_{\infty}\leq C\kappa_{1}^{6}\log^{2}(n)/(np)$ for some constant
$C>0$.

As the residual term $\bm{r}$ has small magnitude, we may analyze the
properties of $\widehat{\bm{\theta}}-\bm{\theta}^{\star}$ by focusing on the
leading term $-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}$.

##### Proof outline of Theorem 2.

Using the approximation
$\widehat{\bm{\theta}}-\bm{\theta}^{\star}\approx-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}$,
we can turn our attention to the concentration of

$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|^{2}\approx\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}=\widehat{\bm{\epsilon}}^{\top}\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}.$

This is a quadratic form of sub-Gaussian random variables. The Hanson-Wright
inequality tells us that
$\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}$ concentrates
around its expectation—$\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})$. The complete
proof is deferred to Section D for Theorem 2.

##### Proof outline of Theorem 3.

With Lemma 4 in hand, we can expand
$\bm{L}_{Lz}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})$ to reach

$\displaystyle\bm{L}_{Lz}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})$
$\displaystyle\approx-(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{B}\widehat{\bm{\epsilon}}$
$\displaystyle=-\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\underbrace{\left[Y_{ji}^{(l)}-\sigma(\theta_{i}^{\star}-\theta_{j}^{\star})\right](\bm{L}_{Lz}^{\dagger})^{1/2}(\bm{e}_{i}-\bm{e}_{j})}_{\eqqcolon\bm{u}_{ij}^{(l)}}.$

The main component of the last expression is a sum of bounded random variables
$\bm{u}_{ij}^{(l)}$ that are independent conditional on $\mathcal{G}_{Y}$.
This sum has zero mean and variance

$\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}z_{ij}(\bm{L}_{Lz}^{\dagger})^{1/2}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}(\bm{L}_{Lz}^{\dagger})^{1/2}=\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}.$

Therefore applying multivariate CLT, we obtain

$-(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{B}\widehat{\bm{\epsilon}}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}\right).$

See the complete proof in Section E.

###  4.4 Proof outline of Lemma 4

Due to the importance and usefulness of Lemma 4, we provide a sketch of its
proof and leave the full one to Section C.

The proof is inspired by the proof of Theorem 1 in [Che23], which analyzes MLE
via the trajectory of the preconditioned gradient descent (PGD) dynamic
starting from ground truth. More precisely, letting
$\bm{\theta}^{0}=\bm{\theta}^{\star}$, we consider the PGD iterates defined by

$\bm{\theta}^{t+1}=\bm{\theta}^{t}-\eta\bm{L}_{Lz}^{\dagger}\nabla\mathcal{L}(\bm{\theta}^{t}),$

where $\eta>0$ is the step size of PGD. [Che23] shows that this dynamic
converges to $\widehat{\bm{\theta}}$. We proceed one step further by
establishing precise distributional characterization of $\hat{\bm{\theta}}$
via analyzing PGD. With Taylor expansion, the gradient can be decomposed into

$\nabla\mathcal{L}(\bm{\theta}^{t})=\bm{L}_{Lz}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\bm{B}\widehat{\bm{\epsilon}}+\bm{r}^{t},$

where $\bm{r}^{t}$ is a residual vector with small magnitude. Then the PGD
update becomes

$\bm{\theta}^{t+1}-\bm{\theta}^{\star}=\left(1-\eta\right)(\bm{\theta}^{t}-\bm{\theta}^{\star})-\eta\left(\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}-\bm{L}_{Lz}^{\dagger}\bm{r}^{t}\right).$

We establish Lemma 4 by solving this recursive relation. More specifically, as
$\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}$ does not depend on $t$
and $\|\bm{L}_{Lz}^{\dagger}\bm{r}^{t}\|_{\infty}$ can be controlled for each
step $t$, taking $t\rightarrow\infty$, we see that

$\widehat{\bm{\theta}}-\bm{\theta}^{\star}=\lim_{t\rightarrow\infty}\bm{\theta}^{t}-\bm{\theta}^{\star}=-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\bm{r}$

for some residual term $\bm{r}$ that is well controlled in $\ell_{\infty}$
norm.

##  5 Experiment

In this section, we demonstrate the empirical performance of
$\mathsf{RP\text{-}MLE}$ and its variant $\mathsf{MRP\text{-}MLE}$ using both
simulated and real data.

###  5.1 Simulations

We begin by using simulated data to validate our theoretical results and
compare our estimators with existing ones for the Rasch model.

The data generating process follows the model specified in Section 2.1. Unless
specified otherwise, in each trial, the user and item parameters are randomly
drawn from

$\widetilde{\bm{\zeta}}^{\star}\sim\mathcal{N}(0,\bm{I}_{n}),\qquad\text{and}\qquad\widetilde{\bm{\theta}}^{\star}\sim\mathcal{N}(0,\bm{I}_{m}).$

Afterwards, $\bm{\zeta}^{\star}$ and $\bm{\theta}^{\star}$ is computed by
shifting $\widetilde{\bm{\zeta}}^{\star}$ and
$\widetilde{\bm{\theta}}^{\star}$ to zero mean.

####  5.1.1 $\ell_{\infty}$ estimation error

Here, we investigate the $\ell_{\infty}$ estimation error. First, we validate
the theoretical result in $\ell_{\infty}$ estimation error of
$\mathsf{RP\text{-}MLE}$. Second, we show how much advantage the
$\mathsf{MRP\text{-}MLE}$ brings through multiple runs of data splitting.
Lastly, we compare our methods with existing comparison-based algorithms,
including the case where $\kappa_{1},\kappa_{2}$ are large.

##### Validating the theoretical result.

Theorem 1 tells us that the $\ell_{\infty}$ error scales as $1/\sqrt{np}$. As
shown in Figure 1, $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ is
very close to being linear w.r.t. both $1/\sqrt{n}$ and $1/\sqrt{p}$. This is
consistent with our theoretical predictions.

(a) $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ with v.s.
$1/\sqrt{n}$. The parameter is chosen to be $m=50,p=0.1$ and $n$ varies from
$10000$ to $40000$. Each point represents the average of 1000 trials. (b)
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ v.s. $1/\sqrt{p}$.
The parameter is chosen to be $m=50,n=10000$ and $p$ varies from $1/9$ to $1$.
Each point represents the average of 1000 trials.

Figure 1: Estimation error
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ of
$\mathsf{RP\text{-}MLE}$ with varying $n$ and $p$.

##### Multiple runs in $\mathsf{MRP\text{-}MLE}$.

As we have discussed after Theorem 1, the random data splitting could incur a
small loss of information. We have provided a remedy $\mathsf{MRP\text{-}MLE}$
(Algorithm 2) by averaging over multiple runs with independent data splitting.
In Figure 2, we show the empirical evidence of its effectiveness. For each
trial, we run 50 independent data splitting and obtain the 50 estimates
$\widehat{\bm{\theta}}^{(i)}$. We then report the $\ell_{\infty}$ error of the
running average, i.e.,
$\|\frac{1}{k}\sum_{i=1}^{k}\widehat{\bm{\theta}}^{(i)}-\bm{\theta}^{\star}\|$
for $k=1,\ldots,50$. As a baseline for comparison, we also run pseudo MLE on
the same datasets, which uses all possible pairs of comparison. We can see
that by averaging over more runs of data splittings, $\mathsf{MRP\text{-}MLE}$
gives a better estimation. In particular, one can see that with sufficient
number of independent data splittings, the average estimation error matches
and even beats pseudo MLE.

Figure 2: The $\ell_{\infty}$ estimation error of $\mathsf{MRP\text{-}MLE}$
with different numbers of data splitting. The dashed line is the performance
of Pseudo MLE. The parameters are chosen to be $m=50,p=0.1,n=10000$. Each
point is averaged over 1000 trials.

##### Comparison with existing estimator.

We then compare our algorithms with two other comparison-based algorithms:
pseudo MLE (PMLE) and the spectral method in [NZ23]. We plot the
$\ell_{\infty}$ error $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$
against $n$ in Figure 4. We can see that all four algorithms follow a similar
trend. The performance of $\mathsf{MRP\text{-}MLE}$ is on par with PMLE and
better than the spectral method. Our proposed algorithm not only enjoys an
improved theoretical guarantee but is also competitive in practice.

##### Performance with large $\kappa_{1},\kappa_{2}$.

While we assume $\kappa=\max\\{\kappa_{1},\kappa_{2}\\}=O(1)$ is most of this
article, the scenario with large $\kappa$ can be practically relevant. Here we
compare the $\ell_{\infty}$ error of different methods when the condition
number is large. Fixing $\kappa$, we draw the user and item parameters with

$\widetilde{\bm{\zeta}}^{\star}\sim\mathrm{Unif}(0,\log(\kappa))\qquad\text{and}\qquad\widetilde{\bm{\theta}}^{\star}\sim\mathrm{Unif}(0,\log(\kappa))$

and compute $\bm{\zeta}^{\star}$ and $\bm{\theta}^{\star}$ by shifting
$\widetilde{\bm{\zeta}}^{\star}$ and $\widetilde{\bm{\theta}}^{\star}$ to zero
mean. Figure 4 shows the performance of different estimators as $\kappa$
varies. The MLE-based approaches including $\mathsf{RP\text{-}MLE}$ and
$\mathsf{MRP\text{-}MLE}$ achieve better $\ell_{\infty}$ error than the
spectral method when $\kappa$ is large.

Figure 3:  $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ v.s. $n$
using Spectral method, PMLE, $\mathsf{RP\text{-}MLE}$, and
$\mathsf{MRP\text{-}MLE}$ using 20 data splittings. The parameter is chosen to
be $m=50,p=0.1$ and $n$ varies from $5000$ to $20000$. The result is averaged
over 1000 trials.

Figure 4:  $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|_{\infty}$ v.s.
$\log(\kappa)$ using Spectral method, PMLE, $\mathsf{RP\text{-}MLE}$, and
$\mathsf{MRP\text{-}MLE}$ using 20 data splittings. The parameter is chosen to
be $m=50,p=0.1,n=20000$ and $\kappa$ varies from $1$ to $e^{10}$. The result
is averaged over 1000 trials.

####  5.1.2 Top-$K$ recovery

In this part, we investigate the performance of different algorithms in
top-$K$ recovery. We set $\theta_{i}^{\star}=(1-K/m)\Delta_{K}$ for $i\leq K$
and $\theta_{i}^{\star}=(-K/m)\Delta_{K}$ otherwise. For any estimator
$\widehat{\bm{\theta}}$, we define top-$K$ recovery rate to be

$\frac{1}{K}\left|\left\\{i\leq K:i\in A_{K}\right\\}\right|,$

where $A_{K}$ is an arbitrary $K$-element set such that
$\widehat{\theta}_{i}\geq\widehat{\theta}_{j}$ for any $i\in A_{K},j\notin
A_{K}$. We compare the top-$K$ recovery rate of PMLE and the spectral method
in [NZ23] with $\mathsf{RP\text{-}MLE}$ and $\mathsf{MRP\text{-}MLE}$ in
Figure 6. The recovery rate of PMLE, spectral method and
$\mathsf{MRP\text{-}MLE}$ is similar, indicating again that our algorithm
performs well in practice.

Figure 5: Top-$K$ recovery rate using spectral method, PMLE,
$\mathsf{RP\text{-}MLE}$, and $\mathsf{MRP\text{-}MLE}$ using 20 data
splittings. The parameter is chosen to be $m=10000,m=50,p=0.1,K=5$ and
$\Delta_{K}$ varies from 0.1 to 0.7. The result is averaged over 1000 trials.

Figure 6:  The relative deviation of
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ from
$\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}$ v.s. $n$ for both
$\widetilde{z}=\widehat{z}$ and $\widetilde{z}=z$. The parameter is chosen to
be $p=0.1$, $n$ varies from $10000$ to $50000$ and $m=n/500$. The result is
averaged over 1000 trials.

####  5.1.3 Refined $\ell_{2}$ estimation error

In this part of the section, we show the empirical evidence of Theorem 2. In
Theorem 2 we have shown that the $\ell_{2}$ error concentrate around
$\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}$ for
$\widetilde{z}\in\left\\{\widehat{z},z\right\\}$. In each trial we compute the
following quantity

$\frac{\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}\right|}{\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}}$

for both $\widetilde{z}=\widehat{z}$ and $\widetilde{z}=z$. This measures the
relative deviation of $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ from
$\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}$. In Figure 6 we
consider the regime where $p$ and $n/m$ is fixed. In this case, Theorem 2
implies that

$\frac{\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}\right|}{\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}}\lesssim\frac{\sqrt{\frac{1}{np}}+\frac{\sqrt{m}}{np}}{\sqrt{\frac{m}{np}}}\lesssim\frac{1}{\sqrt{n}}.$

In other words, $\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ concentrate
tightly around $\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})$. We can see
that the deviation is very small between $\widetilde{z}=\widehat{z}$ and
$\widetilde{z}=z$. In both cases, the relative deviation of
$\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|$ from
$\sqrt{\mathrm{Trace}(\bm{L}_{L\tilde{z}}^{\dagger})}$ decreases as $n$ and
$m$ increase as expected.

####  5.1.4 Inference

Lastly, we validate our inferential results in Section 3.3.

##### Asymptotic normality.

We first check the asymptotic normality presented in Theorem 3 and Corollary
1. For simplicity, we focus on the first coordinate. Theorem 3 and Corollary 1
have claimed that

$\left[\bm{L}_{Lz}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})\right]_{1}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,1-\frac{1}{m}\right)$

and

$\frac{\widehat{\theta}_{1}-\theta_{1}^{\star}}{\sqrt{[\bm{L}_{Lz}^{\dagger}]_{11}}}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,1\right).$

For simplicity we set both $\bm{\zeta}^{\star}$ and $\bm{\theta}^{\star}$ to
be all zero vectors. In Figure 7, we present the QQ-plot comparing the
generated instances and the theoretical distributions. In both cases, the
spread of the term of interest matches the normality claimed in Theorem 3 and
Corollary 1.

Figure 7: Quantile-quantile plots of
$\left[\bm{L}_{Lz}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})\right]_{1}$
and
$(\widehat{\theta}_{1}-\theta_{1}^{\star})/\sqrt{[\bm{L}_{Lz}^{\dagger}]_{11}}$
v.s. zero-mean normal distributions. The estimate $\widehat{\bm{\theta}}$ is
computed with $\mathsf{RP\text{-}MLE}$. The parameters are set to be
$n=10000,m=10,p=0.5$. Data is collected over 10000 trials.

##### Coverage of the confidence intervals.

We compare the empirical coverage rate of the two-sided confidence intervals
in Proposition 3 with the theoretical ones. Without loss of generality we
consider the confidence interval for the first coordinate. In each trial, we
compute the empirical coverage rate

$1-\widehat{\alpha}=\frac{1}{m}\sum_{i=1}^{m}\mathds{1}\\{\mathcal{C}_{i}^{-}(\alpha/2)\leq\theta_{i}^{\star}\leq\mathcal{C}_{i}^{+}(\alpha/2)\\}$

and report the average $1-\widehat{\alpha}$ over the trials. We do this for
multiple $\alpha$’s and compare $1-\alpha$ with $1-\widehat{\alpha}$ in Table
1. The empirical coverage rate $1-\widehat{\alpha}$ is very close to the
theoretical coverage rate $1-\alpha$.

$1-\alpha$ | 0.8 | 0.9 | 0.95 | 0.98 | 0.99 | 0.995 | 0.999
---|---|---|---|---|---|---|---
$1-\widehat{\alpha}$ | 0.801145 | 0.900635 | 0.95108 | 0.980395 | 0.99025 | 0.99523 | 0.99904
Table 1:  Theoretical coverage rate $1-\alpha$ and empirical coverage rate
$1-\widehat{\alpha}$ of the two-sided confidence intervals
$[\mathcal{C}_{i}^{-}(\alpha/2),\mathcal{C}_{i}^{+}(\alpha/2)]$. The
parameters are set to be $n=10000,m=20,p=0.1$. The result is averaged over
$10000$ trials.

###  5.2 LSAT real data

We study a real dataset (LSAT) on the Law School Admissions Test from [DBL70].
LSAT has full observation of 1000 people answering 5 problems. Foro each
person-item pair it records whether the answer is correct or not. The second
row in Table 2 lists how many people answer each problem correctly. From the
first look, Problem 3 seems to be the hardest question. We proceed to quantify
the hardness of these problems under the Rasch model and infer how confident
we are in claiming it is the hardest. We use $\mathsf{RP\text{-}MLE}$ to get
an estimate of the latent scores and calculate two-sided confidence intervals
at level $\alpha=0.01$ for each coordinate with the framework in Section 3.3.
See Table 2 for the result. Higher latent score here means higher difficulty.
The estimated parameters agree with the total number of correct answers in
reverse order. We can also see that the confidence interval lower bound of
$\theta_{3}^{\star}$ is larger than the largest confidence interval upper
bound of
$\theta_{1}^{\star},\theta_{2}^{\star},\theta_{4}^{\star},\theta_{5}^{\star}$.
With Bonferroni correction, we have a 0.95 confidence level to say Problem 3
is the most difficult problem.

Problem | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
Total correct | 924 | 709 | 553 | 763 | 870
$\theta$ estimate | -1.2139 | 0.4093 | 1.1123 | 0.2977 | -0.6053
CI lower bound | -1.5957 | 0.1449 | 0.8321 | 0.0234 | -0.9050
CI upper bound | -0.8322 | 0.6736 | 1.3926 | 0.5720 | -0.3057
Table 2:  Latent score estimate calculated using $\mathsf{RP\text{-}MLE}$ and
confidence interval calculated with the framework in Section 3.3. The
significance level is chosen to be $\alpha=0.01$ for each coordinate.

Now we assume Problem 3 is the top-$1$ item and investigate the top-$1$
recovery rate of different algorithms on LSAT under incomplete observation. In
each trial we randomly select $\widetilde{n}$ people, and for each of them we
randomly select their outcome on $\widetilde{m}$ problems. We then estimate
$\bm{\theta}^{\star}$ using the subsampled data with different methods and
compare their top-1 recovery rate, i.e., the proportion of trials where the
top-$1$ item is correctly identified. In Figure 8 we see results similar to
the simulation. Our algorithm $\mathsf{MRP\text{-}MLE}$ has a similar recovery
rate compared to PMLE and spectral method.

Figure 8: Top-1 recovery rate using Spectral method, PMLE,
$\mathsf{RP\text{-}MLE}$, and $\mathsf{MRP\text{-}MLE}$ using 20 data
splittings. The parameters are chosen to be $\widetilde{m}=4$ and
$\widetilde{n}$ varies from 100 to 200. The result is averaged over 10000
trials.

##  6 Discussion

In this paper, we propose two new likelihood-based estimators
$\mathsf{RP\text{-}MLE}$ and $\mathsf{MRP\text{-}MLE}$ for item parameter
estimation in the Rasch model. Both enjoy optimal finite sample estimation
guarantee and $\mathsf{RP\text{-}MLE}$ is further amenable to tight
uncertainty quantification. All this is achieved even when the user-item
response data are extremely sparse (cf. [NZ23]). Below, we identify several
questions that are interesting for further investigation:

  * •

Does PMLE achieve optimal theoretical guarantee?__ In our experiments, pseudo
MLE has shown a similar performance to $\mathsf{MRP\text{-}MLE}$. This
naturally leads to the question of whether PMLE can enjoy the same theoretical
guarantee. This is relevant to our work because our methods can be viewed as a
modification of pseudo MLE by incorporating random disjoint pairing to
decouple statistical dependency among paired $Y_{ij}$’s. It remains unclear
whether such dependency is a fundamental bottleneck.

  * •

Inference for multiple data splitting. In $\mathsf{MRP\text{-}MLE}$ (Algorithm
2), we average the estimates over multiple runs of data splitting to improve
the empirical result over $\mathsf{RP\text{-}MLE}$. While this is effective
for estimation, we do not have an immediate analog for inference as the
estimations $\\{\widehat{\bm{\theta}}^{(i)}\\}_{i\in[n_{\mathrm{split}}]}$ are
dependent on each other. It is an interesting question whether we can
characterize the distribution of $\widehat{\bm{\theta}}$ in
$\mathsf{MRP\text{-}MLE}$ and whether it leads to better inference, e.g.,
tighter confidence intervals.

  * •

Extending random pairing to other models in IRT. Some IRT models parameterize
the latent score of users and items differently from the Rasch model. For
instance, the two-parameter logistic model (2PL) assumes that $X_{ti}$, the
response of user $t$ to item $i$, follows the law

$\mathbb{P}[X_{ti}=1]=\frac{\exp(a_{t}^{\star}+b_{t}^{\star}\theta_{i}^{\star})}{1+\exp(a_{t}^{\star}+b_{t}^{\star}\theta_{i}^{\star})},$

where $\bm{\theta}^{\star}$ is the latent score of the items while
$\bm{a}^{\star}$ and $\bm{b}^{\star}$ are two parameter vectors for the users.
Unlike the Rasch model, in the 2PL model,

$\mathbb{P}[X_{ti}>X_{tj}\mid X_{ti}\neq
X_{tj}]=\frac{\exp(b_{t}^{\star}\theta_{i}^{\star})}{\exp(b_{t}^{\star}\theta_{i}^{\star})+\exp(b_{t}^{\star}\theta_{j}^{\star})}$

is not independent of the user parameter $\bm{b}^{\star}$. Therefore reduction
to the BTL model is no longer true in this case. It is interesting and non-
trivial to extend the idea of random pairing to the 2PL model.

### Acknowledgement

YY and CM were partially supported by the National Science Foundation via
grant DMS-2311127.

## References

  * [BT52] Ralph Allan Bradley and Milton E Terry.  Rank analysis of incomplete block designs: I. the method of paired comparisons.  Biometrika, 39(3/4):324–345, 1952.
  * [CFMW19] Yuxin Chen, Jianqing Fan, Cong Ma, and Kaizheng Wang.  Spectral method and regularized mle are both optimal for top-k ranking.  Annals of statistics, 47(4):2204, 2019.
  * [CGZ22] Pinhan Chen, Chao Gao, and Anderson Y Zhang.  Partial recovery for top-k ranking: Optimality of mle and suboptimality of the spectral method.  The Annals of Statistics, 50(3):1618–1652, 2022.
  * [Che23] Yanxi Chen.  Ranking from pairwise comparisons in general graphs and graphs with locality.  arXiv preprint arXiv:2304.06821, 2023.
  * [Cho82] Bruce Choppin.  A fully conditional estimation procedure for rasch model parameters.  1982\.
  * [CLLY21] Yunxiao Chen, Xiaoou Li, Jingchen Liu, and Zhiliang Ying.  Item response theory–a statistical framework for educational and psychological measurement.  arXiv preprint arXiv:2108.08604, 2021.
  * [CLOX23] Yunxiao Chen, Chengcheng Li, Jing Ouyang, and Gongjun Xu.  Statistical inference for noisy incomplete binary matrix.  Journal of Machine Learning Research, 24(95):1–66, 2023.
  * [CS15] Yuxin Chen and Changho Suh.  Spectral mle: Top-k rank aggregation from pairwise comparisons.  In International Conference on Machine Learning, pages 371–380. PMLR, 2015.
  * [DBL70] R Darrell Bock and Marcus Lieberman.  Fitting a response model for $n$ dichotomously scored items.  Psychometrika, 35(2):179–197, 1970.
  * [DC10] André F De Champlain.  A primer on classical test theory and item response theory for assessments in medical education.  Medical education, 44(1):109–117, 2010.
  * [ER13] Susan E Embretson and Steven P Reise.  Item response theory.  Psychology Press, 2013.
  * [FBC05] James F Fries, Bonnie Bruce, and David Cella.  The promise of promis: using item response theory to improve assessment of patient-reported outcomes.  Clinical and experimental rheumatology, 23(5):S53, 2005.
  * [FHY22] Jianqing Fan, Jikai Hou, and Mengxin Yu.  Uncertainty quantification of mle for entity ranking with covariates.  arXiv preprint arXiv:2212.09961, 2022.
  * [FLWY22] Jianqing Fan, Zhipeng Lou, Weichen Wang, and Mengxin Yu.  Ranking inferences based on the top choice of multiway comparisons.  arXiv preprint arXiv:2211.11957, 2022.
  * [FLWY23] Jianqing Fan, Zhipeng Lou, Weichen Wang, and Mengxin Yu.  Spectral ranking inferences based on general multiway comparisons.  arXiv preprint arXiv:2308.02918, 2023.
  * [Gho95] Malay Ghosh.  Inconsistent maximum likelihood estimators for the rasch model.  Statistics & Probability Letters, 23(2):165–170, 1995.
  * [GSZ23] Chao Gao, Yandi Shen, and Anderson Y Zhang.  Uncertainty quantification in the bradley–terry–luce model.  Information and Inference: A Journal of the IMA, 12(2):1073–1140, 2023.
  * [JKSO16] Minje Jang, Sunghyun Kim, Changho Suh, and Sewoong Oh.  Top-$k$ ranking from pairwise comparisons: When spectral ranking is optimal.  arXiv preprint arXiv:1603.04153, 2016.
  * [LFL23] Yue Liu, Ethan X Fang, and Junwei Lu.  Lagrangian inference for ranking problems.  Operations research, 71(1):202–223, 2023.
  * [Lin99] John M Linacre.  Understanding rasch measurement: estimation methods for rasch measures.  Journal of outcome measurement, 3:382–405, 1999.
  * [LNB68] FM Lord, MR Novick, and Allan Birnbaum.  Statistical theories of mental test scores.  1968\.
  * [LSR22] Wanshan Li, Shamindra Shrotriya, and Alessandro Rinaldo.  $\ell_{\infty}$ bounds of the MLE in the BTL model under general comparison graphs.  In Uncertainty in Artificial Intelligence, pages 1178–1187. PMLR, 2022.
  * [Luc59] R Duncan Luce.  Individual choice behavior.  1959\.
  * [Mol95] Ivo W. Molenaar.  Estimation of Item Parameters, pages 39–51.  Springer New York, New York, NY, 1995.
  * [NZ22] Duc Nguyen and Anderson Ye Zhang.  A spectral approach to item response theory.  Advances in neural information processing systems, 35:38818–38830, 2022.
  * [NZ23] Duc Nguyen and Anderson Ye Zhang.  Optimal and private learning from human response data.  In International Conference on Artificial Intelligence and Statistics, pages 922–958. PMLR, 2023.
  * [Ost59] Alexander M Ostrowski.  A quantitative formulation of sylvester’s law of inertia.  Proceedings of the National Academy of Sciences, 45(5):740–744, 1959.
  * [Rai19] Martin Raič.  A multivariate berry–esseen theorem with explicit constants.  2019\.
  * [Ras60] Georg Rasch.  Studies in mathematical psychology: I. probabilistic models for some intelligence and attainment tests.  1960\.
  * [Rob21] Alexander Robitzsch.  A comparison of estimation methods for the rasch model.  2021\.
  * [RV13] Mark Rudelson and Roman Vershynin.  Hanson-wright inequality and sub-gaussian concentration.  2013\.
  * [Sch92] Bernhard A Schmitt.  Perturbation bounds for matrix square roots and pythagorean sums.  Linear algebra and its applications, 174:215–227, 1992.
  * [Spi07] Daniel A Spielman.  Spectral graph theory and its applications.  In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS’07), pages 29–38. IEEE, 2007.
  * [Tro15] Joel A. Tropp.  An introduction to matrix concentration inequalities.  Foundations and Trends® in Machine Learning, 8(1-2):1–230, 2015.
  * [Ver18] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [VHSA20] Steven M Van Hauwaert, Christian H Schimpf, and Flavio Azevedo.  The measurement of populist attitudes: Testing cross-national scales using item response theory.  Politics, 40(1):3–21, 2020.
  * [Wed73] Per-Åke Wedin.  Perturbation theory for pseudo-inverses.  BIT Numerical Mathematics, 13:217–232, 1973.
  * [YCOM24] Yuepeng Yang, Antares Chen, Lorenzo Orecchia, and Cong Ma.  Top-${K}$ ranking with a monotone adversary.  Conference on learning theory, 2024.
  * [Zwi95] Aeilko H Zwinderman.  Pairwise parameter estimation in rasch models.  Applied Psychological Measurement, 19(4):369–375, 1995.

##  Appendix A Proof of Fact 1

Expanding the probability of the random events, we have that

$\displaystyle\mathbb{P}[X_{ti}<X_{tj}\mid X_{ti}\neq X_{tj}]$
$\displaystyle=\frac{\mathbb{P}[X_{ti}=0,X_{tj}=1]}{\mathbb{P}[X_{ti}=0,X_{tj}=1\text{\text{
or }}X_{ti}=1,X_{tj}=0]}$
$\displaystyle=\frac{e^{\zeta_{t}^{\star}}e^{\theta_{j}^{\star}}}{(e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}})(e^{\zeta_{t}^{\star}}+e^{\theta_{j}^{\star}})}\Bigg{/}\left(\frac{e^{\zeta_{t}^{\star}}e^{\theta_{j}^{\star}}}{(e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}})(e^{\zeta_{t}^{\star}}+e^{\theta_{j}^{\star}})}+\frac{e^{\theta_{i}^{\star}}e^{\zeta_{t}^{\star}}}{(e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}})(e^{\zeta_{t}^{\star}}+e^{\theta_{j}^{\star}})}\right)$
$\displaystyle=\frac{e^{\zeta_{t}^{\star}}e^{\theta_{j}^{\star}}}{e^{\zeta_{t}^{\star}}(e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}})}$
$\displaystyle=\frac{e^{\theta_{j}^{\star}}}{e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}}}.$

Now consider $\mathbb{P}[X_{ti}\neq X_{tj}]$, we have

$\displaystyle\mathbb{P}[X_{ti}\neq X_{tj}]$
$\displaystyle=\frac{e^{\zeta_{t}^{\star}}e^{\theta_{j}^{\star}}}{(e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}})(e^{\zeta_{t}^{\star}}+e^{\theta_{j}^{\star}})}+\frac{e^{\theta_{i}^{\star}}e^{\zeta_{t}^{\star}}}{(e^{\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}})(e^{\zeta_{t}^{\star}}+e^{\theta_{j}^{\star}})}$
$\displaystyle=\frac{e^{\theta_{j}^{\star}-\zeta_{t}^{\star}}+e^{\theta_{i}^{\star}-\zeta_{t}^{\star}}}{(1+e^{\theta_{i}^{\star}-\zeta_{t}^{\star}})(1+e^{\theta_{j}^{\star}-\zeta_{t}^{\star}})}.$
(16)

Let $f:[1/\kappa_{2},\kappa_{2}]^{2}\rightarrow\mathbb{R}$ defined by

$f(a,b)\coloneqq\frac{a+b}{(1+a)(1+b)}.$

Its partial derivatives are

$\frac{\partial}{\partial
a}f(a,b)=\frac{b^{2}-1}{(1+a)^{2}(1+b)^{2}}\quad\text{and}\quad\frac{\partial}{\partial
b}f(a,b)=\frac{a^{2}-1}{(1+a)^{2}(1+b)^{2}}.$

It is now easy to see that the minimum or maximum of $f$ can only happen if
$(a,b)=(1,1)$ or $(a,b)\in\\{1/\kappa_{2},\kappa_{2}\\}^{2}$. After comparing
the value of $f$ at these points, we conclude that $f$ achieves minimum at

$f(1/\kappa_{2},1/\kappa_{2})=f(\kappa_{2},\kappa_{2})=\frac{2\kappa_{2}}{(1+\kappa_{2})^{2}}.$

By the definition of $\kappa_{2}$,
$|\theta_{l}^{\star}-\zeta_{t}^{\star}|\leq\log(\kappa_{2})$ for any
$l\in[m]$. Then (16) fits the definition of $f$ and the proof is completed.

##  Appendix B Degree and spectral properties of the comparison graphs

In this section, we present the analysis for lemmas that characterize the
degree and spectral properties of the comparison graphs. We start with a lemma
that controls the degrees in $\mathcal{G}_{X}$, and then prove Lemmas 2 and 3.

###  B.1 Degree range of $\mathcal{G}_{X}$

Recall that $n_{t}$ is the number of neighbors of user $t$ in
$\mathcal{G}_{X}$. Furthermore, we denote $m_{i}$ as the number of users that
is compared with problem $i$ and at least another item, i.e.,

$m_{i}\coloneqq\left|\\{t:(t,i)\in\mathcal{E}_{X},n_{t}\geq 2\\}\right|.$ (17)

The following lemma controls the size of $n_{t}$ and $m_{i}$.

######  Lemma 5 (Degree bounds in $\mathcal{G}_{X}$).

Suppose that $np\geq C\log(n)$ for some large enough constant $C>0$ and that
$m\leq n^{\alpha}$ for some constant $\alpha>0$. Then with probability at
least $1-2n^{-10}$, for all $i\in[m]$, we have

$\frac{1}{4}np\leq m_{i}\leq\frac{3}{2}np.$ (18)

Moreover, with probability at least $1-n^{-10}$, for all $t\in[n]$, we have

$n_{t}\leq\left(\frac{3}{2}mp\right)\vee 165\log(n).$ (19)

###### Proof.

We prove the two claims in the lemma sequentially.

Fix any $t,i$. One has

$\displaystyle\mathbb{P}\left[t:(t,i)\in\mathcal{E}_{X},n_{t}\geq 2\right]$
$\displaystyle=\mathbb{P}\left[(t,i)\in\mathcal{E}_{X}\right]-\mathbb{P}\left[(t,i)\in\mathcal{E}_{X}\text{
and }n_{t}=1\right]$ $\displaystyle=p-p(1-p)^{m-1}$ $\displaystyle\geq
p(1-e^{-(m-1)p})$ $\displaystyle\geq\frac{1}{2}p,$ (20)

as long as $mp\geq 2$. Let $\mu_{i}\coloneqq\mathbb{E}[m_{i}]$. By the
linearity of expectation, we have

$np/2\leq\mu_{i}\leq\sum_{t}\mathbb{P}\left[(t,i)\in\mathcal{E}_{X}\right]=np.$
(21)

Fix $i\in[m]$. Since the sampling is independent with different $t$, by the
Chernoff bound,

$\displaystyle\mathbb{P}[|m_{i}-\mu_{i}|\leq(1/2)\mu_{i}]$ $\displaystyle\leq
2e^{-\frac{1}{12}\mu_{i}}\leq 2e^{-\frac{1}{24}np}\leq m^{-1}n^{-10}$

as long as $np\geq C\log(n)$ for large enough constant $C$. Applying (21) and
union bound on $i\in[m]$ yields (18).

Moving on to (19), we first consider the case where $mp\geq 110\log(n)$. By
Chernoff bound,

$\mathbb{P}[n_{t}\geq(3/2)mp]\leq e^{-\frac{1}{10}mp}\leq n^{-11}.$

In the case of $mp<110\log(n)$, the quantity $\mathbb{P}\left[n_{t}\geq
165\log(n)\right]$ clear decreases as $p$ decreases. So we may use the case
$mp=110\log(n)$ to bound this quantity and conclude that

$\mathbb{P}\left[n_{t}\geq 165\log(n)\right]\leq n^{-11}.$

Finally we apply union bound on $t\in[n]$ to reach (19).

∎

###  B.2 Proof of Lemma 2

The assumption of this lemma allows us to invoke Lemma 5. For the upper bound
of $d_{i}$, since (18) is true,

$\displaystyle d_{i}$ $\displaystyle=\sum_{j:j\neq i}L_{ij}$
$\displaystyle=\sum_{t:(t,i)\in\mathcal{E}_{X}}\sum_{j:j\neq i}L_{ij}^{t}$
$\displaystyle\overset{(\text{i})}{=}\sum_{t:(t,i)\in\mathcal{E}_{X},n_{t}\geq
2}\sum_{j:j\neq i}L_{ij}^{t}$
$\displaystyle\overset{(\text{ii})}{\leq}\sum_{t:(t,i)\in\mathcal{E}_{X},n_{t}\geq
2}\sum_{j:j\neq i}R_{ij}^{t}\leq m_{i}\leq\frac{3}{2}np.$

Here (i) holds since $L_{ij}^{t}$ can only be 0 when $n_{t}\leq 1$, and (ii)
holds since $L_{ij}^{t}\leq R_{ij}^{t}$ by definition. For the lower bound of
$d_{i}$, notice that for any $(t,i)$, $\sum_{j}L_{ij}^{t}$ is either $0$ or 1.
Fix $\mathcal{E}_{X}$ and only consider randomness on $L_{ij}^{t}$. By
Hoeffding’s inequality,

$\displaystyle\mathbb{P}\left\\{d_{i}-\sum_{t:(t,i)\in\mathcal{E}_{X},n_{t}\geq
2}\mathbb{E}\left[\sum_{j:j\neq
i}L_{ij}^{t}\mid(t,i)\in\mathcal{E}_{X}\right]\leq-\frac{1}{12\kappa_{2}}np\right\\}$
$\displaystyle\leq\exp\left(-\frac{(1/72)\kappa_{2}^{-2}n^{2}p^{2}}{m_{i}}\right)$
$\displaystyle\leq\exp\left(-\frac{\kappa_{2}^{-2}np}{108}\right)$
$\displaystyle\leq m^{-1}n^{-10}$ (22)

as long as $np\geq 1200\kappa_{2}^{2}\log(n)$. The second to last inequality
uses (18). For each $(t,i)\in\mathcal{E}_{X}$,

$\displaystyle\mathbb{E}\left[\sum_{j}L_{ij}^{t}\mid(t,i)\in\mathcal{E}_{X}\right]$
$\displaystyle=\sum_{j}\mathbb{P}\left[R_{ij}^{t}=1\mid(t,i)\in\mathcal{E}_{X}\right]\mathbb{P}\left[\sum_{j}L_{ij}^{t}=1\mid
R_{ij}^{t}=1\right]$
$\displaystyle\geq\sum_{j}\mathbb{P}\left[R_{ij}^{t}=1\mid(t,i)\in\mathcal{E}_{X}\right]\frac{2\kappa_{2}}{(1+\kappa_{2})^{2}}$
$\displaystyle\geq\frac{2\lfloor
n_{t}/2\rfloor}{n_{t}}\cdot\frac{2\kappa_{2}}{(1+\kappa_{2})^{2}}\geq\frac{1}{3\kappa_{2}}$

as long as $n_{t}\geq 2$. The first inequality here uses Fact 1. Then by
definition of $m_{i}$ in (17),

$\sum_{t:(t,i)\in\mathcal{E}_{X},n_{t}\geq
2}\mathbb{E}\left[\sum_{j}L_{ij}^{t}\mid(t,i)\in\mathcal{E}_{X}\right]\geq\frac{1}{3\kappa_{2}}m_{i}.$
(23)

Combining (22), (23) and (18),

$d_{i}\geq\frac{1}{3\kappa_{2}}m_{i}-\frac{1}{12\kappa_{2}}np\geq\frac{1}{24\kappa_{2}}np.$

Applying union bound over $i\in[m]$ yields the desired result.

###  B.3 Proof of Lemma 3

We first consider the spectrum of $\bm{L}_{L}$. Recall that

$\displaystyle\bm{L}_{L}$
$\displaystyle=\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}$
$\displaystyle=\sum_{t=1}^{n}\underbrace{\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}^{t}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}}_{\bm{L}_{L}^{t}}.$

For the upper bound, it is clear from Lemmas 10 and 2 that
$\lambda_{1}(\bm{L}_{L})\leq 2\max_{i}d_{i}\leq 3np$. For the lower bound, we
use the matrix Chernoff inequality (see Section 5 of [Tro15]). Let
$\bm{R}\in\mathbb{R}^{(n-1)\times n}$ be a partial isometry such that
$\bm{R}\bm{R}^{\top}=\bm{I}_{m-1}$ and $\bm{R}\bm{1}=\bm{0}$. Then
$\lambda_{m-1}(\bm{L}_{L})=\lambda_{m-1}(\bm{R}\bm{L}_{L}\bm{R}^{\top})$. For
any $t\in[n]$, by (19),

$0\leq\lambda_{m-1}(\bm{R}\bm{L}_{L}^{t}\bm{R}^{\top})\leq\lambda_{1}(\bm{R}\bm{L}_{L}^{t}\bm{R}^{\top})=\lambda_{1}(\bm{L}_{L}^{t})\leq
2.$

The last inequality follows from Lemma 10 since $\sum_{j}\bm{L}_{ij}^{t}\leq
1$ for any $t\in[n]$ and $i\in[m]$. By Fact 1,

$\mathbb{P}[L_{ij}^{t}=1\mid
R_{ij}^{t}=1]\geq\frac{2\kappa_{2}}{(1+\kappa_{2})^{2}}\geq 1/(2\kappa_{2}).$

Then

$\displaystyle\lambda_{m-1}(\mathbb{E}\bm{R}\bm{L}_{L}\bm{R}^{\top})$
$\displaystyle=\lambda_{m-1}\left(\bm{R}\sum_{t=1}^{n}\sum_{i>j}\mathbb{E}L_{ij}^{t}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{R}^{\top}\right)$
$\displaystyle\geq\frac{1}{2\kappa_{2}}\lambda_{m-1}\left(\bm{R}\sum_{t=1}^{n}\sum_{i>j}\mathbb{E}R_{ij}^{t}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{R}^{\top}\right)$

Moreover
$\sum_{i>j}\mathbb{E}R_{ij}^{t}\geq\frac{1}{2}(\mathbb{E}[n_{t}]-1)=(mp-1)/2$,
where the $-1$ accounts for possible unpaired $X_{ti}$. By symmetry
$\mathbb{E}R_{ij}^{t}$ is the same for any $(i,j)$. Then for any $(i,j)$,

$\mathbb{E}R_{ij}^{t}\geq\frac{mp-1}{2}/\binom{m}{2}\geq\frac{p}{2m}$

as long as $mp\geq 2$. Thus

$\displaystyle\lambda_{m-1}(\mathbb{E}\bm{R}\bm{L}_{L}\bm{R}^{\top})$
$\displaystyle\geq\frac{1}{2\kappa_{2}}\lambda_{m-1}\left(\bm{R}\sum_{t=1}^{n}\sum_{i>j}\frac{p}{2m}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{R}^{\top}\right)$
$\displaystyle=\frac{1}{2\kappa_{2}}\cdot n\cdot\frac{p}{2m}\cdot m$
$\displaystyle=\frac{np}{4\kappa_{2}}.$

Now invoke the matrix Chernoff inequality, we have

$\displaystyle\mathbb{P}\left\\{[\lambda_{m-1}(\bm{R}\bm{L}_{L}\bm{R}^{\top})]\leq\frac{np}{8\kappa_{2}}\right\\}$
$\displaystyle\leq
m\cdot\left[\frac{e^{-1/2}}{(1/2)^{1/2}}\right]^{\frac{np}{4\kappa_{2}}/2}\leq
n^{-10}$

as long as $np\geq C\kappa_{2}\log(n)$ for some large enough constant $C$.

The spectrum of $\bm{L}_{Lz}$ comes directly from the spectrum of
$\bm{L}_{L}$. Recall

$\bm{L}_{Lz}=\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}z_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}.$

By Lemma 9,

$\displaystyle\lambda_{1}(\bm{L}_{Lz})$
$\displaystyle=\max_{\bm{v}:\|\bm{v}\|=1}\bm{v}^{\top}\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}z_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle\leq\frac{1}{4}\max_{\bm{v}:\|\bm{v}\|=1}\bm{v}^{\top}\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle=\frac{1}{4}\lambda_{1}(\bm{L}_{L})\leq np,$

and

$\displaystyle\lambda_{m-1}(\bm{L}_{Lz})$
$\displaystyle=\min_{\bm{v}:\|\bm{v}\|=1,\bm{v}^{\top}\bm{1}_{m}=0}\bm{v}^{\top}\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}z_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle\geq\frac{1}{4\kappa_{1}}\min_{\bm{v}:\|\bm{v}\|=1,\bm{v}^{\top}\bm{1}_{m}=0}\bm{v}^{\top}\sum_{(i,j)\in\mathcal{E}_{Y},i>j}L_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle=\frac{1}{4\kappa_{1}}\lambda_{m-1}(\bm{L}_{L})\geq\frac{np}{16\kappa_{1}\kappa_{2}}.$

##  Appendix C Proof of Lemma 4

We study the MLE $\widehat{\bm{\theta}}$ by analyzing the iterates of
preconditioned gradient descent starting from the ground truth. Let
$\bm{\theta}_{0}=\bm{\theta}^{\star}$ be the starting point and $\eta>0$ be
the stepsize that is small enough. The preconditioned gradient descent
iterates are given by

$\bm{\theta}^{t+1}=\bm{\theta}^{t}-\eta\bm{L}_{Lz}^{\dagger}\nabla\mathcal{L}(\bm{\theta}^{t}).$
(24)

Consider the Taylor expansion of $\nabla\mathcal{L}(\bm{\theta}^{t})$, we have

$\displaystyle\nabla\mathcal{L}(\bm{\theta}^{t})$
$\displaystyle=\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\left(\left(-Y_{ji}^{(l)}+\sigma(\theta_{i}^{t}-\theta_{j}^{t})\right)(\bm{e}_{i}-\bm{e}_{j})\right)$
$\displaystyle=\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\left(\left(-\epsilon_{ji}^{(l)}-\sigma(\theta_{i}^{\star}-\theta_{j}^{\star})+\sigma(\theta_{i}^{t}-\theta_{j}^{t})\right)(\bm{e}_{i}-\bm{e}_{j})\right)$
$\displaystyle=\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\left(\left(-\epsilon_{ji}^{(l)}+\sigma^{\prime}(\theta_{i}^{\star}-\theta_{j}^{\star})(\delta_{i}^{t}-\delta_{j}^{t})+\frac{1}{2}\sigma^{\prime\prime}(\xi_{ij}^{t})(\delta_{i}^{t}-\delta_{j}^{t})^{2}\right)(\bm{e}_{i}-\bm{e}_{j})\right).$

Here $\bm{\delta}^{t}\coloneqq\bm{\theta}^{t}-\bm{\theta}^{\star}$ and for all
$(i,j)$, $\xi_{ij}^{t}\in\mathbb{R}$ is some number that lies between
$\theta_{i}^{\star}-\theta_{j}^{\star}$ and $\theta_{i}^{t}-\theta_{j}^{t}$.
As $\sigma^{\prime}(\theta_{i}^{\star}-\theta_{j}^{\star})=z_{ij}$ and
$\delta_{i}^{t}-\delta_{j}^{t}=(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{\delta}^{t}$,
we have

$\nabla\mathcal{L}(\bm{\theta}^{t})=\bm{L}_{Lz}\bm{\delta}^{t}-\bm{B}\widehat{\bm{\epsilon}}+\bm{r}^{t},$

where
$\bm{r}^{t}=\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}[\frac{1}{2}\sigma^{\prime\prime}(\xi_{ij}^{t})(\delta_{i}^{t}-\delta_{j}^{t})^{2}(\bm{e}_{i}-\bm{e}_{j})]$.
Feeding this into (24), we have

$\bm{\delta}^{t+1}=\left(1-\eta\right)\bm{\delta}^{t}-\eta\left(\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}-\bm{L}_{Lz}^{\dagger}\bm{r}^{t}\right)$
(25)

By definition $\bm{\delta}^{0}=\bm{0}$. Applying this recursive relation $t-1$
times we obtain

$\displaystyle\bm{\delta}^{t}$
$\displaystyle=-\eta\sum_{i=0}^{t-1}(1-\eta)^{i}\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\sum_{i=0}^{t-1}(1-\eta)^{t-1-i}\bm{L}_{Lz}^{\dagger}\bm{r}^{i}$
$\displaystyle=-\left[1-(1-\eta)^{t}\right]\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\sum_{i=0}^{t-1}(1-\eta)^{t-1-i}\bm{L}_{Lz}^{\dagger}\bm{r}^{i}.$

At this point, we invoke an existing result on these terms that have been
studied in [YCOM24] for a more general setting. The proof of Lemma 2 combined
with Lemmas 3 and 4 in [YCOM24] reveal the following properties of (24).

######  Lemma 6.

Instate the assumptions of Theorem 1. Suppose that

$\kappa_{1}^{3}\frac{(d_{\max})^{2}\log^{2}(n)}{(\lambda_{m-1}(\bm{L}_{Lz}))^{3}}\leq
C_{1}\kappa_{1}\sqrt{\frac{\log(n)}{\lambda_{m-1}(\bm{L}_{Lz})}},$ (26)

for some constant $C_{1}>0$. Then with probability at least $1-n^{-10}$, the
precondition gradient descent dynamic satisfies the following properties:

  1. 1.

There exists a unique minimizer $\widehat{\bm{\theta}}$ of (2).

  2. 2.

There exist some $\alpha_{1},\alpha_{2}$ obeying $0<\alpha_{1}\leq\alpha_{2}$
such that any $t\in\mathbb{N}$,

$\|\bm{\theta}^{t}-\widehat{\bm{\theta}}\|_{\bm{L}_{Lz}}\leq(1-\eta\alpha_{1})^{t}\|\bm{\theta}^{0}-\widehat{\bm{\theta}}\|_{\bm{L}_{Lz}},$

provided that $0<\eta\leq 1/\alpha_{2}$.

  3. 3.

For any $k,l$ and iteration $t\geq 0$,

$\left|\left(\theta_{k}^{t}-\theta_{l}^{t}\right)-\left(\theta_{k}^{\star}-\theta_{l}^{\star}\right)\right|\leq
C_{2}\kappa_{1}\sqrt{\frac{\log(n)}{\lambda_{m-1}(\bm{L}_{L})}}$ (27)

for some constant $C_{2}>0$.

  4. 4.

For any $k,l$ and iteration $t\geq 0$,

$\left|(\bm{e}_{k}-\bm{e}_{l})^{\top}\bm{L}_{Lz}^{\dagger}\bm{r}^{t}\right|\leq
C_{3}\kappa_{1}^{3}\frac{(d_{\max})^{2}\log^{2}(n)}{(\lambda_{m-1}(\bm{L}_{L}))^{3}}$

for some constant $C_{3}>0$.

Lemmas 3 and 2 imply that

$\frac{np}{4\kappa_{1}\kappa_{2}}\leq\lambda_{m-1}(\bm{L}_{Lz})\qquad\text{and}\qquad
d_{\max}\leq\frac{3}{2}np.$

Then the condition (26) holds as long as $np\geq
C_{4}\kappa_{1}^{2}\kappa_{2}\log^{3}(n)$ for some large enough constant
$C_{4}$. Invoke Lemma 6 to see that for any $(k,l,t)$,

$\left|(\bm{e}_{k}-\bm{e}_{l})^{\top}\bm{L}_{Lz}^{\dagger}\bm{r}^{t}\right|\leq
C_{5}\kappa_{1}^{6}\frac{\log^{2}(n)}{np},$ (28)

where $C_{5}>0$ is some constant. Furthermore, the convergence given by Lemma
6 implies that

$\displaystyle\widehat{\bm{\delta}}$
$\displaystyle=\lim_{t\rightarrow\infty}\bm{\delta}^{t}=-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\underbrace{\eta\lim_{t\rightarrow\infty}\sum_{i=0}^{t-1}(1-\eta)^{t-1-i}\bm{L}_{Lz}^{\dagger}\bm{r}^{i}}_{\eqqcolon\bm{r}}.$
(29)

It remains to control the $\ell_{\infty}$ norm of $\bm{r}$. For any $(k,l$),
(28) shows that

$|r_{k}-r_{j}|\leq\eta\lim_{t\rightarrow
0}\sum_{i=0}^{t-1}(1-\eta)^{t-1-i}C_{5}\kappa_{1}^{6}\frac{\log^{2}(n)}{np}=C_{5}\kappa_{1}^{6}\frac{\log^{2}(n)}{np}.$

As $\bm{1}^{\top}\bm{L}_{Lz}^{\dagger}=0$, the above inequality implies that

$\displaystyle|r_{k}|$ $\displaystyle=\left|\frac{1}{m}\cdot
m\bm{e}_{k}^{\top}\bm{r}\right|=\left|\frac{1}{m}\sum_{l=1}^{m}(\bm{e}_{k}-\bm{e}_{l})^{\top}\bm{r}\right|$
$\displaystyle=\left|\frac{1}{m}\sum_{l=1}^{m}(r_{k}-r_{l})\right|\leq
C_{5}\kappa_{1}^{6}\frac{\log^{2}(n)}{np}.$

The proof is now completed.

##  Appendix D Proof of Theorem 2

We start with the proof of (5). By Lemma 4 we can express the MLE estimation
error $\widehat{\bm{\theta}}-\bm{\theta}^{\star}$ as

$\widehat{\bm{\theta}}-\bm{\theta}^{\star}=-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\bm{r}$
(30)

for $\bm{B},\widehat{\bm{\epsilon}}$ defined in Section 4.3 and
$\bm{r}\in\mathbb{R}^{m}$ is a residual term obeying $\|\bm{r}\|_{\infty}\leq
C_{1}\kappa_{1}^{6}\log^{2}n/(np)$ for some constant $C_{1}$.

We first focus on the main term
$\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}$. Expanding $\bm{B}$ and
$\widehat{\bm{\epsilon}}$, we rewrite it as

$\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}=\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\underbrace{\epsilon_{ij}^{(l)}\bm{L}_{Lz}^{\dagger}(\bm{e}_{i}-\bm{e}_{j})}_{\eqqcolon\bm{u}_{ij}^{(l)}}.$

It is easy to see that conditional on $\mathcal{G}_{Y}$,
$\\{\bm{u}_{ij}^{(l)}\\}_{i,j,l}$ is a set of independent zero-mean random
variables. Thus we can expand
$\mathbb{E}[\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}]$ to be

$\displaystyle\mathbb{E}\left[\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}\right]$
$\displaystyle=\mathbb{E}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\bm{u}_{ij}^{(l)\top}\bm{u}_{ij}^{(l)}$
$\displaystyle=\mathbb{E}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\mathrm{Trace}\left(\bm{u}_{ij}^{(l)}\bm{u}_{ij}^{(l)\top}\right)$
$\displaystyle\overset{\text{(i)}}{=}\mathrm{Trace}\left(\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\bm{L}_{Lz}^{\dagger}z_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{L}_{Lz}^{\dagger}\right)$
$\displaystyle\overset{\text{(ii)}}{=}\mathrm{Trace}(\bm{L}_{Lz}^{\dagger}).$
(31)

Here (i) follows from the equality

$z_{ij}=\mathrm{Var}(\epsilon_{ji}^{(l)})=\sigma^{\prime}(\theta_{i}^{\star}-\theta_{j}^{\star})=\frac{e^{\theta_{i}^{\star}}e^{\theta_{j}^{\star}}}{(e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}})^{2}},$

and (ii) follows from the definition of $\bm{L}_{Lz}$. By Lemma 3,

$\frac{m}{2np}\leq(m-1)\lambda_{m-1}(\bm{L}_{Lz}^{\dagger})\leq\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})\leq
m\|\bm{L}_{Lz}^{\dagger}\|\leq\frac{16\kappa_{1}\kappa_{2}m}{np}.$ (32)

Moreover, $\\{\epsilon_{ij}^{(l)}\\}_{i,j,l}$ is a set of sub-Gaussian random
variable with variance proxy $1/z_{ij}$ (see, e.g., Section 2.5 in [Ver18])
and $1/z_{ij}\leq 4\kappa_{1}$ by Lemma 9. Applying Hanson-Wright inequality
(see [RV13]), for any scalar $a>0$ we have

$\displaystyle\mathbb{P}\left[\left|\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}-\mathbb{E}\left[\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}\right]\right|>a\right]$
(33) $\displaystyle\quad\leq
2\exp\left[-C_{2}\left(\frac{a^{2}}{\left(4\kappa_{1}\right)^{4}\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|_{\mathrm{F}}^{2}}\wedge\frac{a}{\left(4\kappa_{1}\right)^{2}\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|}\right)\right]$

for some constant $C_{2}>0$. For
$\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|_{\mathrm{F}}$
and $\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|$ we
have

$\displaystyle\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|$
$\displaystyle\leq\|\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\|_{\mathrm{F}}$
$\displaystyle=\sqrt{\mathrm{Trace}\left(\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\right)}$
$\displaystyle=\sqrt{\mathrm{Trace}\left(\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\bm{B}^{\top}\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\bm{B}\bm{B}^{\top}\right)}$
$\displaystyle\overset{\text{(i)}}{=}\sqrt{\mathrm{Trace}\left(\bm{L}_{Lz}^{\dagger}\bm{L}_{Lz}^{\dagger}\right)}$
$\displaystyle\overset{\text{(ii)}}{\leq}\frac{\sqrt{m}}{\lambda_{m-1}(\bm{L}_{Lz})}\overset{\text{(iii)}}{\leq}\frac{16\kappa_{1}\kappa_{2}\sqrt{m}}{np}.$

Here (i) follows from the fact that $\bm{B}\bm{B}^{\top}=\bm{L}_{Lz}$, (ii)
follows from Lemma 3 and the fact that $\mathrm{Trace}(\bm{M})\leq
m\|\bm{M}\|$ for any $m\times m$ matrix $\bm{M}$, and (iii) follows from Lemma
3. Now substitute $a=C_{3}\kappa_{1}^{3}\kappa_{2}\sqrt{m}\log(n)/(np)$ in
(33) for some large enough constant $C_{3}$. We have that with probability at
least $1-2n^{-10}$,

$\left|\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}-\mathbb{E}\left[\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}\right]\right|\leq\frac{C_{3}\kappa_{1}^{3}\kappa_{2}\sqrt{m}\log(n)}{np}.$
(34)

Combining this with (31) and (32),

$\displaystyle\left|\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|-\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\right|$
$\displaystyle=\frac{\left|\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|^{2}-\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})\right|}{\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|+\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}}$
$\displaystyle\leq\frac{C_{3}\kappa_{1}^{3}\kappa_{2}\sqrt{m}\log(n)/(np)}{\sqrt{m/(2np)}}$
$\displaystyle\leq C_{4}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{\log(n)}{np}}$

for some constant $C_{4}>0$.

Substituting (34) and (31) into (30), we have that for some constant
$C_{1},C_{2}$,

$\displaystyle\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\right|$
$\displaystyle\leq\left|\|\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}\|-\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\right|+\|\bm{r}\|$
$\displaystyle\leq
C_{4}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{\log(n)}{np}}+\frac{C_{1}\kappa_{1}^{6}\sqrt{m}\log^{2}(n)}{np}.$
(35)

The proof of (5) is now completed. For (6), the following lemma connects
$\widehat{z}$ and $z$. The proof is deferred to the end of this section.

######  Lemma 7.

Instate the assumptions of Theorem 1, then with probability at least
$1-C_{5}n^{-10}$ for some constant $C_{5}>0$,

$\|\bm{L}_{Lz}^{\dagger}-\bm{L}_{L\hat{z}}^{\dagger}\|\leq\frac{C_{6}\kappa_{1}^{7/2}\kappa_{2}^{2}}{(np)^{3/2}}$

for some large enough constant $C_{6}$.

Combining this lemma with Weryl’s inequality, we have that

$\left|\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})-\mathrm{Trace}(\bm{L}_{L\hat{z}}^{\dagger})\right|\leq\frac{C_{6}m\kappa_{1}^{7/2}\kappa_{2}^{2}}{(np)^{3/2}}.$

Then by (32),

$\displaystyle\left|\sqrt{\mathrm{Trace}(\bm{L}_{L\hat{z}}^{\dagger})}-\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}\right|$
$\displaystyle=\frac{\left|\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})-\mathrm{Trace}(\bm{L}_{L\hat{z}}^{\dagger})\right|}{\sqrt{\mathrm{Trace}(\bm{L}_{Lz}^{\dagger})}+\sqrt{\mathrm{Trace}(\bm{L}_{L\hat{z}}^{\dagger})}}$
$\displaystyle=\frac{C_{6}m\kappa_{1}^{7/2}\kappa_{2}^{2}/(np)^{3/2}}{\sqrt{m/(2np)}}\leq\frac{\sqrt{2}C_{6}\kappa_{1}^{7/2}\kappa_{2}^{2}\sqrt{m}}{np}.$

Using triangular inequality, we conclude that

$\left|\|\widehat{\bm{\theta}}-\bm{\theta}^{\star}\|-\sqrt{\mathrm{Trace}(\bm{L}_{L\hat{z}}^{\dagger})}\right|\leq
C_{4}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{\log(n)}{np}}+\frac{(C_{1}\kappa_{1}^{6}+\sqrt{2}C_{6}\kappa_{1}^{7/2}\kappa_{2}^{2})\sqrt{m}\log^{2}(n)}{np}.$

##### Proof of Lemma 7.

Recall $\sigma$ is the sigmoid function and its derivative $\sigma^{\prime}$
is $1$-Lipschitz. By Theorem 1, for all $(i,j)$,

$\displaystyle\left|z_{ij}-\widehat{z}_{ij}\right|$
$\displaystyle=\left|\sigma^{\prime}(\theta_{i}^{\star}-\theta_{j}^{\star})-\sigma^{\prime}(\theta_{i}^{\star}-\theta_{j}^{\star})\right|$
$\displaystyle\leq\left|\left(\widehat{\theta}_{i}-\widehat{\theta}_{j}\right)-\left(\theta_{i}^{\star}-\theta_{j}^{\star}\right)\right|\leq
C_{7}\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{\log(n)}{np}}$

for some constant $C_{7}>0$. Then

$\displaystyle\|\bm{L}_{Lz}-\bm{L}_{L\hat{z}}\|$
$\displaystyle=\max_{\bm{v}\in\mathbb{R}^{m}:\|\bm{v}\|=1}\left|\bm{v}^{\top}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}L_{ij}\left(z_{ij}-\widehat{z}_{ij}\right)(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}\right|$
$\displaystyle\leq\max_{\bm{v}\in\mathbb{R}^{m}:\|\bm{v}\|=1}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\left|z_{ij}-\widehat{z}_{ij}\right|\bm{v}^{\top}L_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle\leq
C_{7}\kappa_{1}\kappa_{2}^{1/2}\sqrt{\frac{\log(n)}{np}}\cdot\|\bm{L}_{L}\|$
$\displaystyle\leq 3C_{7}\kappa_{1}\kappa_{2}^{1/2}\sqrt{np\log(n)},$

where the last line follows from Lemma 3. As $np\geq
C_{8}\kappa_{1}^{4}\kappa_{2}^{3}\log^{2}(n)$ for some large enough constant
$C_{8}$, $\|\bm{L}_{Lz}-\bm{L}_{L\hat{z}}\|\leq np/(32\kappa_{1}\kappa_{2})$.
By Weryl’s inequality and Lemma 3,

$\lambda_{m-1}(\bm{L}_{L\hat{z}})\geq\lambda_{m-1}(\bm{L}_{Lz})-\|\bm{L}_{Lz}-\bm{L}_{L\hat{z}}\|\geq\frac{np}{16\kappa_{1}\kappa_{2}}-\frac{np}{32\kappa_{1}\kappa_{2}}=\frac{np}{32\kappa_{1}\kappa_{2}}.$

This implies that $\|\bm{L}_{Lz}^{\dagger}\|\leq 16\kappa_{1}\kappa_{2}/(np)$
and $\|\bm{L}_{L\hat{z}}^{\dagger}\|\leq 32\kappa_{1}\kappa_{2}/(np)$. Using
the perturbation bound of pseudo-inverse (see Theorem 4.1 in [Wed73]), we have

$\displaystyle\|\bm{L}_{Lz}^{\dagger}-\bm{L}_{L\hat{z}}^{\dagger}\|$
$\displaystyle\leq
3\cdot\|\bm{L}_{Lz}^{\dagger}\|\cdot\|\bm{L}_{L\hat{z}}^{\dagger}\|\cdot\|\bm{L}_{Lz}-\bm{L}_{L\hat{z}}\|$
$\displaystyle\leq\frac{C_{9}\kappa_{1}^{3}\kappa_{2}^{5/2}}{(np)^{3/2}}$

for some constant $C_{9}>0$.

##  Appendix E Proofs for inference

In this section, we present the analysis for the inference results we
presented in Section 3.3.

###  E.1 Proof of Theorem 3

We condition the whole analysis on the high probability event when Lemmas 4
and 3 holds. Let
$\widehat{\bm{\delta}}\coloneqq\widehat{\bm{\theta}}-\bm{\theta}^{\star}$.
Expanding (15), we have

$\displaystyle\widehat{\bm{\delta}}$
$\displaystyle=-\bm{L}_{Lz}^{\dagger}\bm{B}\widehat{\bm{\epsilon}}+\bm{r}$
$\displaystyle=-\bm{L}_{Lz}^{\dagger}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\left[Y_{ji}^{(l)}-\sigma(\theta_{i}^{\star}-\theta_{j}^{\star})\right](\bm{e}_{i}-\bm{e}_{j})+\bm{r},$

where $\|\bm{r}\|_{\infty}\leq C_{1}\kappa_{1}^{6}\log^{2}(n)/(np)$ for some
constant $C_{1}$.

We start with the case with $\widetilde{\bm{z}}=\bm{z}$. Consider
$\bm{L}_{Lz}\widehat{\bm{\delta}}$. As $\bm{L}_{Lz}\bm{1}_{m}=\bm{0}$,
$\lambda_{m-1}(\bm{L}_{Lz})>0$, and
$\widehat{\bm{\delta}}^{\top}\bm{1}_{m}=0$,

$\displaystyle\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}$
$\displaystyle=-(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{B}\widehat{\bm{\epsilon}}+\bm{L}_{Lz}^{1/2}\bm{r}$
$\displaystyle=-\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\underbrace{\left[Y_{ji}^{(l)}-\sigma(\theta_{i}^{\star}-\theta_{j}^{\star})\right](\bm{L}_{Lz}^{\dagger})^{1/2}(\bm{e}_{i}-\bm{e}_{j})}_{\eqqcolon\bm{u}_{ij}^{(l)}}+\bm{L}_{Lz}^{1/2}\bm{r}.$

Conditional on $(\mathcal{E}_{Y},\\{L_{ij}\\})$, $\bm{u}_{ij}^{(l)}$ are
independent random variables. It is easy to see that $\bm{u}_{ij}^{(l)}$ is
zero-mean and has covariance

$\displaystyle\mathbb{E}\left[\bm{u}_{ij}^{(l)}\bm{u}_{ij}^{(l)\top}\right]$
$\displaystyle=z_{ij}(\bm{L}_{Lz}^{\dagger})^{1/2}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}(\bm{L}_{Lz}^{\dagger})^{1/2}.$

It is also bounded in third moment of its spectral norm by

$\mathbb{E}\left[\|\bm{u}_{ij}^{(l)}\|^{3}\right]\leq
2^{3/2}\|\bm{L}_{Lz}^{\dagger}\|^{3/2}\leq\frac{2^{15/2}\kappa_{1}^{3/2}\kappa_{2}^{3/2}}{(np)^{3/2}},$

where the last inequality uses Lemma 3. Summing up across $i,j$ and $l$, we
have

$\displaystyle\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}\sum_{l=1}^{L_{ij}}\mathbb{E}\left[\bm{u}_{ij}^{(l)}\bm{u}_{ij}^{(l)\top}\right]$
$\displaystyle=(\bm{L}_{Lz}^{\dagger})^{1/2}\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}L_{ij}z_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}(\bm{L}_{Lz}^{\dagger})^{1/2}$
$\displaystyle=(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{L}_{Lz}(\bm{L}_{Lz}^{\dagger})^{1/2}$
$\displaystyle=\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}.$

The last line holds since $\bm{L}_{Lz}\bm{1}_{m}=\bm{0}$ and
$\lambda_{m-1}(\bm{L}_{Lz})>0$. Now using multivariate CLT (see, e.g.,
[Rai19]), we have that as
$\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}L_{ij}\rightarrow\infty$,

$\displaystyle-(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{B}\widehat{\bm{\epsilon}}$
$\displaystyle\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}\right).$

As $m$ is fixed, Lemma 2 implies that
$\sum_{i>j:(i,j)\in\mathcal{E}_{Y}}L_{ij}\rightarrow\infty$ is equivalent to
$np\rightarrow\infty$. For the residual term, Lemma 4 and Lemma 3 imply that
with probability at least $1-O(n^{-10})$,

$\|\bm{L}_{Lz}^{1/2}\bm{r}\|\leq\|\bm{L}_{Lz}^{1/2}\|\cdot\|\bm{r}\|\leq\|\bm{L}_{Lz}^{1/2}\|\cdot\sqrt{m}\|\bm{r}\|_{\infty}\leq\frac{C_{1}\sqrt{m}\kappa_{1}^{6}\log^{2}(n)}{\sqrt{np}}.$

Then as $np/\log^{4}(n)\rightarrow\infty$,

$\bm{L}_{Lz}^{1/2}\bm{r}\rightarrow\bm{0}$

in probability. Combining the limiting distribution of both two terms, we have

$\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}=-(\bm{L}_{Lz}^{\dagger})^{1/2}\bm{B}\widehat{\bm{\epsilon}}+\bm{L}_{Lz}^{1/2}\bm{r}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}\right).$

In the case of $\widetilde{\bm{z}}=\widehat{\bm{z}}$, we introduce the
following lemma. The proof is deferred to the end of this section.

######  Lemma 8.

Instate the assumptions of Theorem 1. Then with probability at least
$1-C_{1}n^{-10}$ for some constant $C_{1}$, one has

$\|\bm{L}_{Lz}^{1/2}-\bm{L}_{L\hat{z}}^{1/2}\|\leq
C_{2}\kappa_{1}^{2}\kappa_{2}^{1/2}\sqrt{\log(n)}$

for some constant $C_{2}$.

Now the above lemma and Theorem 1 imply

$\left\|\bm{L}_{L\widehat{z}}^{1/2}\widehat{\bm{\delta}}-\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}\right\|\leq\|\bm{L}_{L\widetilde{z}}^{1/2}-\bm{L}_{Lz}^{1/2}\|\cdot\|\widehat{\bm{\delta}}\|\leq
C_{3}\kappa_{1}^{3}\kappa_{2}\sqrt{\frac{m\log^{2}(n)}{np}}$

for some constant $C_{3}>0$. Therefore, as $np/\log^{4}(n)\rightarrow\infty$,

$\bm{L}_{L\widehat{z}}^{1/2}\widehat{\bm{\delta}}-\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}\rightarrow
0$

in probability and hence

$\bm{L}_{L\widehat{z}}^{1/2}\widehat{\bm{\delta}}=\left(\bm{L}_{L\widehat{z}}^{1/2}\widehat{\bm{\delta}}-\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}\right)+\bm{L}_{Lz}^{1/2}\widehat{\bm{\delta}}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top}\right).$

##### Proof of Lemma 8.

Let $\bm{R}\in\mathbb{R}^{(n-1)\times n}$ be a partial isometry such that
$\bm{R}\bm{R}^{\top}=\bm{I}_{m-1}$ and $\bm{R}\bm{1}=\bm{0}$. Then

$\|\bm{L}_{L\hat{z}}^{1/2}-\bm{L}_{Lz}^{1/2}\|=\|\bm{R}(\bm{L}_{L\hat{z}}^{1/2}-\bm{L}_{Lz}^{1/2})\bm{R}^{\top}\|.$

We also observe that
$(\bm{R}\bm{L}_{L\widetilde{z}}\bm{R}^{\top})^{1/2}=\bm{R}\bm{L}_{L\widetilde{z}}^{1/2}\bm{R}^{\top}$
for $\widetilde{z}\in\\{z,\widehat{z}\\}$. As shown in the proof of Lemma 7,

$\|\bm{L}_{L\hat{z}}-\bm{L}_{Lz}\|\leq C_{1}\kappa_{1}^{3/2}\sqrt{np\log(n)}$

for some constant $C_{1}>0$.

Lemma 3 implies that

$\lambda_{m-1}(\bm{R}\bm{L}_{Lz}^{1/2}\bm{R}^{\top})=\left[\lambda_{m-1}(\bm{L}_{Lz})\right]{}^{1/2}\geq\frac{\sqrt{np}}{4\kappa_{1}^{1/2}\kappa_{2}^{1/2}}.$

Then by Lemma 2.1 in [Sch92],

$\|\bm{R}(\bm{L}_{L\widehat{z}}^{1/2}-\bm{L}_{Lz}^{1/2})\bm{R}^{\top}\|\leq\frac{\|\bm{R}(\bm{L}_{L\widehat{z}}-\bm{L}_{Lz})\bm{R}^{\top}\|}{\lambda_{m-1}(\bm{R}\bm{L}_{Lz}^{1/2}\bm{R}^{\top})}\leq
4C_{1}\kappa_{1}^{2}\kappa_{2}^{1/2}\sqrt{\log(n)}.$

###  E.2 Proof of Corollary 1

Observe that

$\frac{\widehat{\theta}_{i}-\theta_{i}^{\star}}{\sqrt{[\bm{L}_{L\tilde{z}}^{\dagger}]_{ii}}}=\frac{\bm{e}_{i}^{\top}(\bm{L}_{L\tilde{z}}^{\dagger})^{1/2}\bm{L}_{L\tilde{z}}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})}{\sqrt{[\bm{L}_{L\tilde{z}}^{\dagger}]_{ii}}}.$

Then Theorem 3 implies

$\frac{\bm{e}_{i}^{\top}(\bm{L}_{L\tilde{z}}^{\dagger})^{1/2}(\bm{I}_{m}-\frac{1}{m}\bm{1}_{m}\bm{1}_{m}^{\top})(\bm{L}_{L\tilde{z}}^{\dagger})^{1/2}\bm{e}_{i}}{[\bm{L}_{L\tilde{z}}^{\dagger}]_{ii}}=1.$

The last inequality follows from the fact that
$\bm{1}_{m}^{\top}(\bm{L}_{L\tilde{z}}^{\dagger})^{1/2}=0$. As a side note, we
compare this with a similar result for the BTL model. Proposition 4.1 in
[GSZ23] studies the uncertainty quantification of MLE in the BTL model with
uniform observation. It says that

$\frac{\widehat{\theta}_{i}-\theta_{i}^{\star}}{\sqrt{\left[\bm{L}_{Lz}\right]_{ii}^{-1}}}\overset{\mathrm{d}}{\rightarrow}\mathcal{N}\left(0,1\right).$

This is close to our result as $[\bm{L}_{Lz}]_{ii}^{-1}$ can be viewed as an
approximation of $[\bm{L}_{Lz}^{\dagger}]_{ii}$ in the context of the BTL
model. To see why, consider the simplified case where
$\bm{\theta}^{\star}=\bm{0}_{n}$ and $L=1$. As $np\rightarrow\infty$, the
Laplacian $\bm{L}_{Lz}$ concentrate to its expectation
$\mathbb{E}[\bm{L}_{Lz}]=np\bm{I}_{n}-p\bm{1}_{n}\bm{1}_{n}^{\top}$. Then
$\bm{L}^{\dagger}\approx(1/np)(\bm{I}_{n}-\frac{1}{n}\bm{1}_{n}\bm{1}_{n}^{\top})$
and $[\bm{L}_{Lz}^{\dagger}]_{ii}\approx[\bm{L}_{Lz}]_{ii}^{-1}$ .

##  Appendix F Auxiliary lemmas

In this section, we gather some auxiliary results that are useful throughout
this paper.

######  Lemma 9 (Range of $z_{ij}$).

Recall

$z_{ij}=\frac{e^{\theta_{i}^{\star}}e^{\theta_{j}^{\star}}}{(e^{\theta_{i}^{\star}}+e^{\theta_{j}^{\star}})^{2}}=\frac{e^{\theta_{i}^{\star}-\theta_{j}^{\star}}}{(1+e^{\theta_{i}^{\star}-\theta_{j}^{\star}})^{2}}.$

For any $(i,j)$,

$\frac{1}{4\kappa_{1}}\leq z_{ij}\leq\frac{1}{4}.$

###### Proof.

Consider the function $f:[0,\infty)\rightarrow\mathbb{R}$ defined by
$f(x)=x/(1+x)^{2}$. It has derivative $(1-x^{2})/(1+x)^{4}$, so it is
increasing at $x\in[0,1)$ and decreasing at $x\in(1,\infty)$. By the
definition of $\kappa_{1}$,
$|\theta_{i}^{\star}-\theta_{j}^{\star}|\leq\log(\kappa_{1})$. Then

$\frac{1}{4\kappa_{1}}\leq f(e^{-\log(\kappa_{1})})\wedge
f(e^{\log(\kappa_{1})})\leq z_{ij}\leq f(1)=\frac{1}{4}.$

∎

######  Lemma 10 (Maximum eigenvalue of Laplacian).

Let
$\bm{L}=\sum_{(i,j):i>j}w_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}$
be a weighted graph Laplacian. Then $\lambda_{1}(\bm{L})\leq
2\max_{i}\sum_{j}w_{ij}$.

###### Proof.

Let $\bm{v}\in\mathbb{R}^{m}$, then

$\displaystyle\bm{v}^{\top}\bm{L}\bm{v}$
$\displaystyle=\bm{v}^{\top}\sum_{(i,j):i>j}w_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}\bm{v}$
$\displaystyle=\sum_{(i,j):i>j}w_{ij}(v_{i}-v_{j})^{2}$ $\displaystyle\leq
2\sum_{(i,j):i>j}w_{ij}(v_{i}^{2}+v_{j}^{2})$ $\displaystyle\leq
2\sum_{i}\sum_{j\neq i}w_{ij}v_{j}^{2}$ $\displaystyle\leq
2\sum_{i}\max_{i}\sum_{j}w_{ij}\|\bm{v}\|^{2}$

So
$\lambda_{1}(\bm{L})=\max_{\bm{v}\in\mathbb{R}^{m},\|\bm{v}\|=1}\bm{v}^{\top}\bm{L}\bm{v}\leq
2\max_{i}\sum_{j}w_{ij}$.∎

######  Lemma 11 (A quantitative version of Sylvester’s law of inertia,
[Ost59]).

For any real symmetric matrix $\bm{A}\in\mathbb{R}^{n\times n}$ and
$\bm{S}\in\mathbb{R}^{n\times n}$ be a non-singular matrix. Then for any
$i\in[n]$, $\lambda_{i}(\bm{S}\bm{A}\bm{S}^{\top})$ lies between
$\lambda_{i}(\bm{A})\lambda_{1}(\bm{S}^{\top}\bm{S})$ and
$\lambda_{i}(\bm{A})\lambda_{n}(\bm{S}^{\top}\bm{S})$.

######  Fact 2.

Let $\mathcal{G}$ be an arbitrary graph with $m$ vertices and let $\bm{L}_{w}$
be a weighted graph Laplacian defined by

$\bm{L}_{w}\coloneqq\sum_{i>j:(i,j)\in\mathcal{G}}w_{ij}(\bm{e}_{i}-\bm{e}_{j})(\bm{e}_{i}-\bm{e}_{j})^{\top}.$

If $w_{ij}>0$ for all $(i,j)\in\mathcal{G}$ and $\mathcal{G}$ is a connected
graph, then $\bm{L}_{w}$ is rank $m-1$, $\bm{L}_{w}\bm{1}_{m}=\bm{0}_{m}$ and
$\bm{L}_{w}^{\dagger}\bm{1}_{m}=0$. Moreover for any $i\in[n-1]$,
$\lambda_{i}(\bm{L}_{w}^{\dagger})=\lambda_{n-i}(\bm{L}_{w})$.

###### Proof.

The fact that $\bm{L}_{w}$ is rank $m-1$ when $\mathcal{G}$ is connected is
well-known. See e.g. [Spi07] for reference. Since $\bm{L}_{w}$ is a real
symmetric matrix, it has an eigendecomposition
$\bm{L}_{w}=\bm{U}\bm{\Sigma}\bm{U}^{\top}$ and then
$\bm{L}_{w}^{\dagger}=\bm{U}\bm{\Sigma}^{\dagger}\bm{U}^{\top}$. The rest
follows from this decomposition and the form of $\bm{L}_{w}$.∎