<figure><img src="image_1.png" /></figure>

HAL
open science

# Ranking and synchronization from pairwise measurements via SVD

Alexandre d'Aspremont, Mihai Cucuringu, Hemant Tyagi

► To cite this version:

Alexandre d'Aspremont, Mihai Cucuringu, Hemant Tyagi. Ranking and synchronization from pairwise measurements via SVD. Journal of Machine Learning Research, 2021, 22 (19), pp.1-63. 10.48550/arXiv.1906.02746 . hal-02340372

HAL Id: hal-02340372

https://hal.science/hal-02340372

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Ranking and synchronization from pairwise measurements via SVD

Alexandre d'Aspremont, Mihai Cucuringu, Hemant Tyagi

► To cite this version:

Alexandre d'Aspremont, Mihai Cucuringu, Hemant Tyagi. Ranking and synchronization from pairwise measurements via SVD. Journal of Machine Learning Research, 2021, 22 (19), pp.1-63. 10.48550/arXiv.1906.02746 . hal-02340372

HAL Id: hal-02340372

https://hal.science/hal-02340372v1

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ranking and synchronization from pairwise measurements via SVD

## Alexandre d'Aspremont

CNRS & Ecole Normale Supérieure, Paris, France

ASPREMON@ENS.FR

## Mihai Cucuringu

Department of Statistics and Mathematical Institute
University of Oxford
The Alan Turing Institute, London, UK

MIHAI.CUCURINGU@STATS.OX.AC.UK

## Hemant Tyagi

Inria, Univ. Lille, CNRS
UMR 8524 - Laboratoire Paul Painlevé, F-59000 *

HEMANT.TYAGI@INRIA.FR

**Editor:** Michael Mahoney

## Abstract

Given a measurement graph $G = (V, E)$ and an unknown signal $r \in \mathbb{R}^n$, we investigate algorithms for recovering $r$ from pairwise measurements of the form $r_i - r_j; \{i, j\} \in E$. This problem arises in a variety of applications, such as ranking teams in sports data and time synchronization of distributed networks. Framed in the context of ranking, the task is to recover the ranking of $n$ teams (induced by $r$) given a small subset of noisy pairwise rank offsets. We propose a simple SVD-based algorithmic pipeline for both the problem of time synchronization and ranking. We provide a detailed theoretical analysis in terms of robustness against both sampling sparsity and noise perturbations with outliers, using results from matrix perturbation and random matrix theory. Our theoretical findings are complemented by a detailed set of numerical experiments on both synthetic and real data, showcasing the competitiveness of our proposed algorithms with other state-of-the-art methods.

**Keywords:** ranking, angular synchronization, spectral algorithms, matrix perturbation theory, singular value decomposition, random matrix theory, low-rank matrix completion.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Related work</td><td>4</td></tr><tr><td>2.1</td><td>Synchronization literature</td><td>5</td></tr><tr><td>2.2</td><td>Ranking literature</td><td>6</td></tr><tr><td>3</td><td>Problem setup and main results</td><td>9</td></tr><tr><td>3.1</td><td>Main idea: SVD-based spectral algorithm</td><td>10</td></tr></table>

<footnote>* Authors are listed in alphabetical order.</footnote>
# Ranking and synchronization from pairwise measurements via SVD

## Alexandre d'Aspremont

CNRS & Ecole Normale Supérieure, Paris, France

ASPREMON@ENS.FR

## Mihai Cucuringu

Department of Statistics and Mathematical Institute
University of Oxford
The Alan Turing Institute, London, UK

MIHAI.CUCURINGU@STATS.OX.AC.UK

## Hemant Tyagi

Inria, Univ. Lille, CNRS
UMR 8524 - Laboratoire Paul Painlevé, F-59000 *

HEMANT.TYAGI@INRIA.FR

**Editor:** Michael Mahoney

## Abstract

Given a measurement graph $G = (V, E)$ and an unknown signal $r \in \mathbb{R}^n$, we investigate algorithms for recovering $r$ from pairwise measurements of the form $r_i - r_j; \{i, j\} \in E$. This problem arises in a variety of applications, such as ranking teams in sports data and time synchronization of distributed networks. Framed in the context of ranking, the task is to recover the ranking of $n$ teams (induced by $r$) given a small subset of noisy pairwise rank offsets. We propose a simple SVD-based algorithmic pipeline for both the problem of time synchronization and ranking. We provide a detailed theoretical analysis in terms of robustness against both sampling sparsity and noise perturbations with outliers, using results from matrix perturbation and random matrix theory. Our theoretical findings are complemented by a detailed set of numerical experiments on both synthetic and real data, showcasing the competitiveness of our proposed algorithms with other state-of-the-art methods.

**Keywords:** ranking, angular synchronization, spectral algorithms, matrix perturbation theory, singular value decomposition, random matrix theory, low-rank matrix completion.

## Contents

1 Introduction 3

2 Related work 4

2.1 Synchronization literature 5

2.2 Ranking literature 6

3 Problem setup and main results 9

3.1 Main idea: SVD-based spectral algorithm 10

<footnote>* Authors are listed in alphabetical order.</footnote>
<table><tbody><tr><td>3.2</td><td>Recovering the global scale: data-driven procedures</td><td>12</td></tr><tr><td>3.3</td><td>Summary of main theoretical results</td><td>14</td></tr><tr><td>4</td><td>Theoretical results for SVD-RS and SVD-NRS</td><td>16</td></tr><tr><td>4.1</td><td>Analysis of Algorithm 1 (SVD-RS)</td><td>16</td></tr><tr><td>4.2</td><td>Analysis of Algorithm 2 (SVD-NRS)</td><td>20</td></tr><tr><td>5</td><td>Proofs</td><td>22</td></tr><tr><td>5.1</td><td>Proof of Theorem 5</td><td>22</td></tr><tr><td>5.2</td><td>Proof of Theorem 6</td><td>24</td></tr><tr><td>5.3</td><td>Proof of Theorem 7</td><td>26</td></tr><tr><td>5.4</td><td>Proof of Theorem 8</td><td>27</td></tr><tr><td>5.5</td><td>Proof of Theorem 9</td><td>28</td></tr><tr><td>5.6</td><td>Proof of Theorem 10</td><td>30</td></tr><tr><td>6</td><td>Matrix completion as a preprocessing step for ranking</td><td>31</td></tr><tr><td>7</td><td>Numerical Experiments</td><td>33</td></tr><tr><td>7.1</td><td>Synthetic data</td><td>34</td></tr><tr><td>7.2</td><td>Real data</td><td>35</td></tr><tr><td>8</td><td>Conclusion and future directions</td><td>46</td></tr><tr><td>A</td><td>Matrix perturbation analysis</td><td>51</td></tr><tr><td>B</td><td>Useful concentration inequalities</td><td>52</td></tr><tr><td>B.1</td><td>Spectral norm of random symmetric matrices</td><td>52</td></tr><tr><td>B.2</td><td>Bernstein inequality</td><td>52</td></tr><tr><td>B.3</td><td>Product of a random matrix (raised to a power) and a fixed vector</td><td>52</td></tr><tr><td>C</td><td>Proof of Theorem 5</td><td>53</td></tr><tr><td>C.1</td><td>Proof of Lemma 11</td><td>53</td></tr><tr><td>C.2</td><td>Proof of Lemma 12</td><td>53</td></tr><tr><td>C.3</td><td>Proof of Lemma 13</td><td>54</td></tr><tr><td>C.4</td><td>Proof of Lemma 14</td><td>54</td></tr><tr><td>D</td><td>Proof of Lemmas from Section 5.2</td><td>56</td></tr><tr><td>D.1</td><td>Proof of Lemma 15</td><td>56</td></tr><tr><td>D.2</td><td>Proof of Lemma 16</td><td>58</td></tr><tr><td>D.3</td><td>Proof of Lemma 17</td><td>60</td></tr><tr><td>D.4</td><td>Proof of Lemma 18</td><td>61</td></tr><tr><td>E</td><td>Proofs of Lemmas from Section 5.5</td><td>62</td></tr><tr><td>E.1</td><td>Proof of Lemma 19</td><td>62</td></tr></tbody></table>
<table><tbody><tr><td>3.2</td><td>Recovering the global scale: data-driven procedures</td><td>12</td></tr><tr><td>3.3</td><td>Summary of main theoretical results</td><td>14</td></tr><tr><td>4</td><td>Theoretical results for SVD-RS and SVD-NRS</td><td>16</td></tr><tr><td>4.1</td><td>Analysis of Algorithm 1 (SVD-RS)</td><td>16</td></tr><tr><td>4.2</td><td>Analysis of Algorithm 2 (SVD-NRS)</td><td>20</td></tr><tr><td>5</td><td>Proofs</td><td>22</td></tr><tr><td>5.1</td><td>Proof of Theorem 5</td><td>22</td></tr><tr><td>5.2</td><td>Proof of Theorem 6</td><td>24</td></tr><tr><td>5.3</td><td>Proof of Theorem 7</td><td>26</td></tr><tr><td>5.4</td><td>Proof of Theorem 8</td><td>27</td></tr><tr><td>5.5</td><td>Proof of Theorem 9</td><td>28</td></tr><tr><td>5.6</td><td>Proof of Theorem 10</td><td>30</td></tr><tr><td>6</td><td>Matrix completion as a preprocessing step for ranking</td><td>31</td></tr><tr><td>7</td><td>Numerical Experiments</td><td>33</td></tr><tr><td>7.1</td><td>Synthetic data</td><td>34</td></tr><tr><td>7.2</td><td>Real data</td><td>35</td></tr><tr><td>8</td><td>Conclusion and future directions</td><td>46</td></tr><tr><td>A</td><td>Matrix perturbation analysis</td><td>51</td></tr><tr><td>B</td><td>Useful concentration inequalities</td><td>52</td></tr><tr><td>B.1</td><td>Spectral norm of random symmetric matrices</td><td>52</td></tr><tr><td>B.2</td><td>Bernstein inequality</td><td>52</td></tr><tr><td>B.3</td><td>Product of a random matrix (raised to a power) and a fixed vector</td><td>52</td></tr><tr><td>C</td><td>Proof of Theorem 5</td><td>53</td></tr><tr><td>C.1</td><td>Proof of Lemma 11</td><td>53</td></tr><tr><td>C.2</td><td>Proof of Lemma 12</td><td>53</td></tr><tr><td>C.3</td><td>Proof of Lemma 13</td><td>54</td></tr><tr><td>C.4</td><td>Proof of Lemma 14</td><td>54</td></tr><tr><td>D</td><td>Proof of Lemmas from Section 5.2</td><td>56</td></tr><tr><td>D.1</td><td>Proof of Lemma 15</td><td>56</td></tr><tr><td>D.2</td><td>Proof of Lemma 16</td><td>58</td></tr><tr><td>D.3</td><td>Proof of Lemma 17</td><td>60</td></tr><tr><td>D.4</td><td>Proof of Lemma 18</td><td>61</td></tr><tr><td>E</td><td>Proofs of Lemmas from Section 5.5</td><td>62</td></tr><tr><td>E.1</td><td>Proof of Lemma 19</td><td>62</td></tr></tbody></table>
# 1. Introduction

Let $r = (r_1, \dots, r_n)^T \in \mathbb{R}^n$ be an unknown signal and $G = ([n], E)$ be an undirected measurement graph. Given a subset of noisy pairwise measurements of the form $r_i - r_j$ for each $\{i, j\} \in E$, the goal is to estimate the original vector $r$. Clearly, this is only possible only up to a global shift. Moreover, when measurements are exact without any measurement noise, one can recover the strength vector $r$ if and only if the graph $G$ is connected, by simply considering a spanning tree, fixing the value of the root node, and traversing the tree while propagating the information by summing the given offsets. For simplicity, we assume the graph is connected, otherwise it is not possible to estimate the offset values between nodes belonging to different connected components of the graph.

Instantiations of the above problem are ubiquitous in engineering, machine learning and computer vision, and have received a great deal of attention in the recent literature. There are two main classes of applications where this problem arises.

* *Time synchronization of wireless networks.* A popular application arises in engineering, and is known as time synchronization of distributed networks (Giridhar and Kumar, 2006; Karp et al., 2003), where clocks measure noisy time offsets $r_i - r_j$, and the goal is to recover $r_1, \dots, r_n \in \mathbb{R}$.

* *Ranking.* A fundamental problem in information retrieval is that of recovering the ordering induced by the latent strengths or scores $r_1, \dots, r_n \in \mathbb{R}$ of a set of $n$ players, that is best reflected by the given set of pairwise comparisons $r_i - r_j$. We refer the reader to Cucuringu (2016) and references therein for a detailed overview.

A naive approach by sequentially propagating the information along spanning trees is doomed to fail in the presence of noise, due to accumulation of the errors. To this end, in order to increase the robustness to noise, one typically aims to simultaneously integrate all the pairwise measurements in a globally consistent framework. This line of thought appears in the literature in the context of the *group synchronization* problem (Singer, 2011), for recovering group elements from noisy pairwise measurements of their ratios. We briefly review the current literature in group synchronization in Section 2.1, which almost exclusively pertains to synchronization over compact groups. The problem we study in this paper can be construed as synchronization over the real line, hence a non-compact group. There exists a very rich literature on ranking, and it is beyond the scope of our work to provide an extensive review of it. Instead, in Section 2.2, we give a brief overview of the literature on ranking and time synchronization, and point out some of the methods that relate to our work. In particular, we focus on spectral methods and approaches that leverage low-rank matrix completion as a pre-processing step.

**Contributions.** We propose SVD-RS, a simple spectral algorithm for ranking and synchronization from pairwise comparisons, along with a normalized version denoted SVD-NRS<fnref n="1" />. We provide a detailed theoretical consistency analysis for both algorithms for a random measurement model (see Section 3) in terms of robustness against sampling sparsity of

<footnote n="1">¹By normalized, we mean that the measurement matrix is normalized by the degree matrix of the graph, see Section 4.</footnote>
# 1. Introduction

Let $r = (r_1, \dots, r_n)^T \in \mathbb{R}^n$ be an unknown signal and $G = ([n], E)$ be an undirected measurement graph. Given a subset of noisy pairwise measurements of the form $r_i - r_j$ for each $\{i, j\} \in E$, the goal is to estimate the original vector $r$. Clearly, this is only possible only up to a global shift. Moreover, when measurements are exact without any measurement noise, one can recover the strength vector $r$ if and only if the graph $G$ is connected, by simply considering a spanning tree, fixing the value of the root node, and traversing the tree while propagating the information by summing the given offsets. For simplicity, we assume the graph is connected, otherwise it is not possible to estimate the offset values between nodes belonging to different connected components of the graph.

Instantiations of the above problem are ubiquitous in engineering, machine learning and computer vision, and have received a great deal of attention in the recent literature. There are two main classes of applications where this problem arises.

* *Time synchronization of wireless networks.* A popular application arises in engineering, and is known as time synchronization of distributed networks (Giridhar and Kumar, 2006; Karp et al., 2003), where clocks measure noisy time offsets $r_i - r_j$, and the goal is to recover $r_1, \dots, r_n \in \mathbb{R}$.

* *Ranking.* A fundamental problem in information retrieval is that of recovering the ordering induced by the latent strengths or scores $r_1, \dots, r_n \in \mathbb{R}$ of a set of $n$ players, that is best reflected by the given set of pairwise comparisons $r_i - r_j$. We refer the reader to Cucuringu (2016) and references therein for a detailed overview.

A naive approach by sequentially propagating the information along spanning trees is doomed to fail in the presence of noise, due to accumulation of the errors. To this end, in order to increase the robustness to noise, one typically aims to simultaneously integrate all the pairwise measurements in a globally consistent framework. This line of thought appears in the literature in the context of the *group synchronization* problem (Singer, 2011), for recovering group elements from noisy pairwise measurements of their ratios. We briefly review the current literature in group synchronization in Section 2.1, which almost exclusively pertains to synchronization over compact groups. The problem we study in this paper can be construed as synchronization over the real line, hence a non-compact group. There exists a very rich literature on ranking, and it is beyond the scope of our work to provide an extensive review of it. Instead, in Section 2.2, we give a brief overview of the literature on ranking and time synchronization, and point out some of the methods that relate to our work. In particular, we focus on spectral methods and approaches that leverage low-rank matrix completion as a pre-processing step.

**Contributions.** We propose SVD-RS, a simple spectral algorithm for ranking and synchronization from pairwise comparisons, along with a normalized version denoted SVD-NRS<fnref n="1" />. We provide a detailed theoretical consistency analysis for both algorithms for a random measurement model (see Section 3) in terms of robustness against sampling sparsity of

<footnote n="1">¹By normalized, we mean that the measurement matrix is normalized by the degree matrix of the graph, see Section 4.</footnote>
the measurement graph and noise level. Additionally, we provide extensive numerical experiments on both synthetic and real data, showing that in certain noise and sparsity regimes, our proposed algorithms perform comparable or better than state-of-the-art methods.

On the theoretical side, our specific contributions can be summarized as follows (see also Section 3.3).

* For SVD-RS, we provide $\ell_2$ and $\ell_\infty$ recovery guarantees for the score vector $r$ (see Theorem 8). For instance, in the setting when $r_i = i$, $\Omega(n \log n)$ measurements suffice for $\ell_2$ recovery, while $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice for $\ell_\infty$ recovery.

* The $\ell_\infty$ analysis of SVD-RS leads to guarantees for rank recovery in terms of the maximum displacement error between the recovered ranking and the ground truth, as in Theorem 7. For e.g., when $r_i = i$, $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice.

* For SVD-NRS, we provide in Theorem 10, $\ell_2$ recovery guarantees for the score vector $r$, and leave the $\ell_\infty$ guarantees for future work, essentially by following a similar, though more intricate, pipeline. Similar to SVD-RS, for $r_i = i$, $\Omega(n \log n)$ measurements suffice for $\ell_2$ recovery.

We remark that all the above recovery results hold with high probability.

**Outline.** The remainder of this paper is organized as follows. Section 2 is a brief survey of the relevant literature, with a focus on group synchronization and ranking. Section 3 starts with the formal setup of the problem, presents the gist of our SVD-based approach along with the two algorithms SVD-RS, SVD-NRS, and summarizes our main theoretical results. Section 4 details and interprets our theoretical results, with the main steps of the proofs outlined in Section 5. Section 6 discusses the low-rank matrix completion problem and its applicability in the setting of this paper. Section 7 contains numerical experiments on various synthetic and real data sets. Finally, Section 8 summarizes our findings along with future research directions. The Appendix contains additional technical details for the theoretical results, as well as further numerical experiments.

**Notation.** Vectors and matrices are denoted by lower case and upper case letters respectively. For a matrix $A \in \mathbb{R}^{m \times n}$ where $m \ge n$, we denote its singular values by $\sigma_1 \ge \cdots \ge \sigma_n$ and the corresponding left (resp. right) singular vectors by $u_i$ (resp. $v_i$). $\|A\|_2$ denotes the spectral norm (largest singular value), $\|A\|_*$ denotes the nuclear norm (sum of the singular values), and $\|A\|_{\max} := \max_{i,j} |A_{ij}|$ denotes the max-norm of $A$. $\mathcal{R}(A)$ denotes the range space of its columns. We denote $e$ to be the all ones column vector. The symbol $\mathcal{U}$ is used to denote the uniform distribution. For positive numbers $a, b$, we denote $a \lesssim b$ (resp. $a \gtrsim b$) to mean that there exists an absolute constant $C > 0$ such that $a \le Cb$ (resp. $a \ge Cb$).

## 2. Related work

This section is a brief survey of the relevant literature for the tasks we consider out in this paper. The first part focuses on the group synchronization problem, a relevant instance of it being that of synchronization of clocks arising often in engineering, while the second part surveys the very rich ranking literature, with an emphasis on spectral methods and approaches that leverage low-rank matrix completion.
the measurement graph and noise level. Additionally, we provide extensive numerical experiments on both synthetic and real data, showing that in certain noise and sparsity regimes, our proposed algorithms perform comparable or better than state-of-the-art methods.

On the theoretical side, our specific contributions can be summarized as follows (see also Section 3.3).

* For SVD-RS, we provide $l_2$ and $l_\infty$ recovery guarantees for the score vector $r$ (see Theorem 8). For instance, in the setting when $r_i = i$, $\Omega(n \log n)$ measurements suffice for $l_2$ recovery, while $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice for $l_\infty$ recovery.

* The $l_\infty$ analysis of SVD-RS leads to guarantees for rank recovery in terms of the maximum displacement error between the recovered ranking and the ground truth, as in Theorem 7. For e.g., when $r_i = i$, $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice.

* For SVD-NRS, we provide in Theorem 10, $l_2$ recovery guarantees for the score vector $r$, and leave the $l_\infty$ guarantees for future work, essentially by following a similar, though more intricate, pipeline. Similar to SVD-RS, for $r_i = i$, $\Omega(n \log n)$ measurements suffice for $l_2$ recovery.

We remark that all the above recovery results hold with high probability.

**Outline.** The remainder of this paper is organized as follows. Section 2 is a brief survey of the relevant literature, with a focus on group synchronization and ranking. Section 3 starts with the formal setup of the problem, presents the gist of our SVD-based approach along with the two algorithms SVD-RS, SVD-NRS, and summarizes our main theoretical results. Section 4 details and interprets our theoretical results, with the main steps of the proofs outlined in Section 5. Section 6 discusses the low-rank matrix completion problem and its applicability in the setting of this paper. Section 7 contains numerical experiments on various synthetic and real data sets. Finally, Section 8 summarizes our findings along with future research directions. The Appendix contains additional technical details for the theoretical results, as well as further numerical experiments.

**Notation.** Vectors and matrices are denoted by lower case and upper case letters respectively. For a matrix $A \in \mathbb{R}^{m \times n}$ where $m \ge n$, we denote its singular values by $\sigma_1 \ge \cdots \ge \sigma_n$ and the corresponding left (resp. right) singular vectors by $u_i$ (resp. $v_i$). $\|A\|_2$ denotes the spectral norm (largest singular value), $\|A\|_*$ denotes the nuclear norm (sum of the singular values), and $\|A\|_{\max} := \max_{i,j} |A_{ij}|$ denotes the max-norm of $A$. $\mathcal{R}(A)$ denotes the range space of its columns. We denote $e$ to be the all ones column vector. The symbol $\mathcal{U}$ is used to denote the uniform distribution. For positive numbers $a, b$, we denote $a \lesssim b$ (resp. $a \gtrsim b$) to mean that there exists an absolute constant $C > 0$ such that $a \le Cb$ (resp. $a \ge Cb$).

## 2. Related work

This section is a brief survey of the relevant literature for the tasks we consider out in this paper. The first part focuses on the group synchronization problem, a relevant instance of it being that of synchronization of clocks arising often in engineering, while the second part surveys the very rich ranking literature, with an emphasis on spectral methods and approaches that leverage low-rank matrix completion.
## 2.1 Synchronization literature

Much of the engineering literature has focused on least-squares approaches for solving the time synchronization problem. In its simplest terms, the approach can be formulated as follows. Let $m = |E|$ denote the number of edges in $G$, and $B$ denote the edge-vertex incidence matrix of size $m \times n$. For $l = \{i, j\} \in E$, with $i < j$, the entries $B_{li}$ and $B_{lj}$ are given by $B_{li} = 1$, $B_{lj} = -1$. If $i' \notin \{i, j\}$ then $B_{li'} = 0$. Let $w \in \mathbb{R}^m$ encode the pairwise rank measurements $C_{ij}$, for all edges $\{i, j\} \in E$. The least-squares solution to the ranking problem can be obtained by solving

$$ \underset{x \in \mathbb{R}^n}{\text{minimize}} \quad \|Bx - w\|_2^2, \qquad (2.1) $$

where $B$ denotes the design matrix with two non-zero entries per row corresponding to the edge $\{i, j\}$ indexed by $l$. In a related line of work, Hirani et al. (2010) show that the problem of least-squares ranking on graphs has far-reaching rich connections with seemingly unrelated areas, such as spectral graph theory and multilevel methods for graph Laplacian systems, Hodge decomposition theory and random clique complexes in topology.

The problem we consider in this paper can also be framed in the context of the group synchronization problem, of finding group elements from noisy measurements of their ratios. For example, synchronization over the special orthogonal group $SO(d)$ consists of estimating a set of $n$ unknown $d \times d$ matrices $R_1, \dots, R_n \in SO(d)$ from noisy measurements of a subset of the pairwise ratios $R_i R_j^{-1}$ via

$$ \underset{R_1, \dots, R_n \in SO(d)}{\text{minimize}} \quad \sum_{\{i,j\} \in E} w_{ij} \|R_i^{-1} R_j - R_{ij}\|_F^2, \qquad (2.2) $$

where $\|\cdot\|_F$ denotes the Frobenius norm, and $w_{ij}$ are non-negative weights representing the confidence in the noisy pairwise measurements $R_{ij}$. Spectral and semidefinite programming (SDP) relaxations for solving an instance of the above synchronization problem were originally introduced and analyzed by Singer (2011), in the context of angular synchronization over the group $SO(2)$ of planar rotations. There, one is asked to estimate $n$ unknown angles $\theta_1, \dots, \theta_n \in [0, 2\pi)$ given $m$ noisy measurements of their offsets $\theta_i - \theta_j \mod 2\pi$. The difficulty of the problem is amplified, on one hand, by the amount of **noise** in the offset measurements, and on the other hand by **sparsity** - the fact that $m \ll \binom{n}{2}$, i.e., only a very small subset of all possible pairwise offsets are measured. In general, one may consider other groups $G$ (such as $SO(d)$, $O(d)$) for which there are available noisy measurements $g_{ij}$ of ratios between the group elements

$$ g_{ij} = g_i g_j^{-1}, g_i, g_j \in G. \qquad (2.3) $$

The set $E$ of pairs $\{i, j\}$ for which a ratio of group elements is available can be realized as the edge set of a graph $G = (V, E)$, $|V| = n$, $|E| = m$, with vertices corresponding to the group elements $g_1, \dots, g_n$, and edges to the available pairwise measurements $g_{ij} = g_i g_j^{-1}$. As long as the group $G$ is compact and has a real or complex representation, one may construct a real or Hermitian matrix (which may also be a matrix of matrices) where the element in the position $\{i, j\}$ is the matrix representation of the measurement $g_{ij}$ (possibly a matrix of size $1 \times 1$, as is the case for $\mathbb{Z}_2$), or the zero matrix if there is no direct measurement for the ratio
## 2.1 Synchronization literature

Much of the engineering literature has focused on least-squares approaches for solving the time synchronization problem. In its simplest terms, the approach can be formulated as follows. Let $m = |E|$ denote the number of edges in $G$, and $B$ denote the edge-vertex incidence matrix of size $m \times n$. For $l = \{i, j\} \in E$, with $i < j$, the entries $B_{li}$ and $B_{lj}$ are given by $B_{li} = 1$, $B_{lj} = -1$. If $i' \notin \{i, j\}$ then $B_{li'} = 0$. Let $w \in \mathbb{R}^m$ encode the pairwise rank measurements $C_{ij}$, for all edges $\{i, j\} \in E$. The least-squares solution to the ranking problem can be obtained by solving

$$ \underset{x \in \mathbb{R}^n}{\text{minimize}} \quad \|Bx - w\|_2^2, \qquad (2.1) $$

where $B$ denotes the design matrix with two non-zero entries per row corresponding to the edge $\{i, j\}$ indexed by $l$. In a related line of work, Hirani et al. (2010) show that the problem of least-squares ranking on graphs has far-reaching rich connections with seemingly unrelated areas, such as spectral graph theory and multilevel methods for graph Laplacian systems, Hodge decomposition theory and random clique complexes in topology.

The problem we consider in this paper can also be framed in the context of the group synchronization problem, of finding group elements from noisy measurements of their ratios. For example, synchronization over the special orthogonal group $SO(d)$ consists of estimating a set of $n$ unknown $d \times d$ matrices $R_1, \dots, R_n \in SO(d)$ from noisy measurements of a subset of the pairwise ratios $R_i R_j^{-1}$ via

$$ \underset{R_1, \dots, R_n \in SO(d)}{\text{minimize}} \quad \sum_{\{i,j\} \in E} w_{ij} \|R_i^{-1} R_j - R_{ij}\|_F^2, \qquad (2.2) $$

where $\|\cdot\|_F$ denotes the Frobenius norm, and $w_{ij}$ are non-negative weights representing the confidence in the noisy pairwise measurements $R_{ij}$. Spectral and semidefinite programming (SDP) relaxations for solving an instance of the above synchronization problem were originally introduced and analyzed by Singer (2011), in the context of angular synchronization over the group $SO(2)$ of planar rotations. There, one is asked to estimate $n$ unknown angles $\theta_1, \dots, \theta_n \in [0, 2\pi)$ given $m$ noisy measurements of their offsets $\theta_i - \theta_j \mod 2\pi$. The difficulty of the problem is amplified, on one hand, by the amount of **noise** in the offset measurements, and on the other hand by **sparsity** - the fact that $m \ll \binom{n}{2}$, i.e., only a very small subset of all possible pairwise offsets are measured. In general, one may consider other groups $\mathcal{G}$ (such as $SO(d)$, $O(d)$) for which there are available noisy measurements $g_{ij}$ of ratios between the group elements

$$ g_{ij} = g_i g_j^{-1}, g_i, g_j \in \mathcal{G}. \qquad (2.3) $$

The set $E$ of pairs $\{i, j\}$ for which a ratio of group elements is available can be realized as the edge set of a graph $G = (V, E)$, $|V| = n$, $|E| = m$, with vertices corresponding to the group elements $g_1, \dots, g_n$, and edges to the available pairwise measurements $g_{ij} = g_i g_j^{-1}$. As long as the group $\mathcal{G}$ is compact and has a real or complex representation, one may construct a real or Hermitian matrix (which may also be a matrix of matrices) where the element in the position $\{i, j\}$ is the matrix representation of the measurement $g_{ij}$ (possibly a matrix of size $1 \times 1$, as is the case for $\mathbb{Z}_2$), or the zero matrix if there is no direct measurement for the ratio
of $g_i$ and $g_j$. For example, the rotation group SO(3) has a real representation using 3 × 3 rotation matrices, and the group SO(2) of planar rotations has a complex representation as points on the unit circle.

The setting we consider in this paper, namely that of recovering points on the real line from a subset of noisy pairwise differences, is essentially synchronization over the non-compact group $\mathbb{R}$. In recent work in the context of ranking from pairwise cardinal and ordinal measurements (Cucuringu, 2016), the real line was compactified by wrapping it over the upper half of the unit circle, making the problem amenable to standard synchronization over the compact group SO(2). The estimated solution allowed for the recovery of the player rankings, after a post-processing step of modding out the best circular permutation. Note that the proposed approach only focused on recovering the individual rankings, and not the magnitude (i.e., strength) of each player, as we propose to do in this paper.

Very recently, Ozyesil et al. (2018) proposed an approach that allows one to consider the synchronization problem in the setting of non-compact groups. The authors leverage a compactification process, that relies on a mapping from a non-compact domain into a compact one, for solving the synchronization problem. The contraction mapping enables one to transform measurements from a Cartan motion group to an associated compact group, whose synchronization solution provides a solution for the initial synchronization problem over the original domain.

## 2.2 Ranking literature

Fogel et al. introduced *SERIAL-RANK* (Fogel et al., 2014, 2016), a ranking algorithm that explicitly leverages the connection with seriation, a classical ordering problem that considers the setting where the user has available a similarity matrix (with ±1 entries) between a set of *n* items, and assumes that there exists an underlying one-dimensional ordering such that the similarity between items decreases with their distance. The authors propose a spectral algorithm *SERIAL-RANK*, and show that under an Erdös-Renyi random graph model and a given noise model, it recovers the underlying ranking with $\Omega(n^{3/2} \log^4 n)$ comparisons where the estimation error is bounded in the $\ell_\infty$ norm.

A popular theme in the literature is the rank aggregation setting (where players meet in multiple matches), where there exists a latent probability matrix $P \in [0, 1]^{n \times n}$, wherein $P_{ij}$ denotes the probability that player *i* defeats player *j*, with $P_{ij} + P_{ji} = 1$. For each pair $(i, j)$, one observes $Y_{ij} \in \{0, 1\} \sim \text{Bern}(P_{ij})$. A number of so-called Random Utility models have been considered in the literature, starting with the seminal work of Bradley-Terry-Luce (BTL) (Bradley and Terry, 1952), which is by far the most popular model considered. In the most basic version of the BTL model, the probability that player *i* beats player *j* is given by $P_{ij} = \frac{w_i}{w_i+w_j}$, where the vector $w \in \mathbb{R}_+^n$ is the parameter vector to be inferred from the data, with $w_i$ being a proxy for the score or strength associated to player *i*. In the Thurstone model (Thurstone, 1927), $P_{ij} = \Phi(s_i - s_j)$, where $\Phi$ denotes the standard normal cumulative distribution function (CDF), and $s \in \mathbb{R}^n$ is the score vector. Negahban et al. proposed *RANK-CENTRALITY* in the context of the rank aggregation problem from multiple ranking systems, under a Bradley-Terry-Luce (BTL) model (Negahban et al., 2012, 2017). *Rank-Centrality* is an iterative algorithm that estimates the ranking scores from the stationary distribution of a certain random walk on the graph of players. Edges encode the outcome of
of $g_i$ and $g_j$. For example, the rotation group SO(3) has a real representation using 3 × 3 rotation matrices, and the group SO(2) of planar rotations has a complex representation as points on the unit circle.

The setting we consider in this paper, namely that of recovering points on the real line from a subset of noisy pairwise differences, is essentially synchronization over the non-compact group $\mathbb{R}$. In recent work in the context of ranking from pairwise cardinal and ordinal measurements (Cucuringu, 2016), the real line was compactified by wrapping it over the upper half of the unit circle, making the problem amenable to standard synchronization over the compact group SO(2). The estimated solution allowed for the recovery of the player rankings, after a post-processing step of modding out the best circular permutation. Note that the proposed approach only focused on recovering the individual rankings, and not the magnitude (i.e., strength) of each player, as we propose to do in this paper.

Very recently, Ozyesil et al. (2018) proposed an approach that allows one to consider the synchronization problem in the setting of non-compact groups. The authors leverage a compactification process, that relies on a mapping from a non-compact domain into a compact one, for solving the synchronization problem. The contraction mapping enables one to transform measurements from a Cartan motion group to an associated compact group, whose synchronization solution provides a solution for the initial synchronization problem over the original domain.

## 2.2 Ranking literature

Fogel et al. introduced *SERIAL-RANK* (Fogel et al., 2014, 2016), a ranking algorithm that explicitly leverages the connection with seriation, a classical ordering problem that considers the setting where the user has available a similarity matrix (with ±1 entries) between a set of *n* items, and assumes that there exists an underlying one-dimensional ordering such that the similarity between items decreases with their distance. The authors propose a spectral algorithm *SERIAL-RANK*, and show that under an Erdös-Renyi random graph model and a given noise model, it recovers the underlying ranking with $\Omega(n^{3/2} \log^4 n)$ comparisons where the estimation error is bounded in the $\ell_\infty$ norm.

A popular theme in the literature is the rank aggregation setting (where players meet in multiple matches), where there exists a latent probability matrix $P \in [0, 1]^{n \times n}$, wherein $P_{ij}$ denotes the probability that player *i* defeats player *j*, with $P_{ij} + P_{ji} = 1$. For each pair $(i, j)$, one observes $Y_{ij} \in \{0, 1\} \sim \text{Bern}(P_{ij})$. A number of so-called Random Utility models have been considered in the literature, starting with the seminal work of Bradley-Terry-Luce (BTL) (Bradley and Terry, 1952), which is by far the most popular model considered. In the most basic version of the BTL model, the probability that player *i* beats player *j* is given by $P_{ij} = \frac{w_i}{w_i+w_j}$, where the vector $w \in \mathbb{R}_+^n$ is the parameter vector to be inferred from the data, with $w_i$ being a proxy for the score or strength associated to player *i*. In the Thurstone model (Thurstone, 1927), $P_{ij} = \Phi(s_i - s_j)$, where $\Phi$ denotes the standard normal cumulative distribution function (CDF), and $s \in \mathbb{R}^n$ is the score vector. Negahban et al. proposed *RANK-CENTRALITY* in the context of the rank aggregation problem from multiple ranking systems, under a Bradley-Terry-Luce (BTL) model (Negahban et al., 2012, 2017). *Rank-Centrality* is an iterative algorithm that estimates the ranking scores from the stationary distribution of a certain random walk on the graph of players. Edges encode the outcome of
pairwise comparisons, and are directed $i \to j$ towards the winner, where the weight captures
the proportion of times $j$ defeated $i$ in case of multiple direct matches. The authors show
that under some assumptions on the connectivity of the underlying graph, their algorithm
estimates the underlying score vector of the BTL model with $\Omega(npoly(\log n))$ comparisons.
Very recently, Agarwal et al. (2018) propose Accelerated Spectral Ranking, a provably faster
spectral ranking algorithm in the setting of the multinomial logit (MNL) and BTL models.
The authors considered a random walk that has a faster mixing time than the random walks
associated with previous algorithms (including Rank-Centrality (Negahban et al., 2017)),
along with improved sample complexity bounds (of the order $npoly(\log n)$) for recovery of
the MNL and BTL parameters.

Cucuringu (2016) introduced SYNC-RANK, formulating the problem of ranking with
incomplete noisy information as an instance of the group synchronization problem over
the group SO(2) of planar rotations. SYNC-RANK starts by making the ansatz that items
are embedded on the upper-half of the unit circle in the complex plane, and transforms
the pairwise input measurements into pairwise angle offsets. It then solves the spectral or
SDP relaxation of synchronization, and finally chooses amongst the circular permutations
of the ranking induced by the synchronized solution by minimizing the number of upsets.
The method was shown to compare favorably with a number of state-of-the-art methods,
including RANK-CENTRALITY (Negahban et al., 2012), SERIAL-RANK (Fogel et al., 2016),
and a SVD-based ranking algorithm similar to the one we consider in this paper, but without
the projection step and its normalized extension, and without any theoretical guarantees
(see Remark 3 for details). Along the same lines, and inspired by the group synchronization
framework, Fanuel and Suykens (2017) propose a certain deformation of the combinatorial
Laplacian, in particular, the so-called dilation Laplacian whose spectrum is leveraged for
ranking in directed networks of pairwise comparisons. The method is reported to perform
well on both synthetic and real data, and enjoys the property that it can place an emphasis
of the top-k items in the ranking.

More recently, De Bacco et al. (2018) proposed SPRING-RANK, an efficient physically-
inspired algorithm, based on solving a linear system of equations, for the task of inferring
hierarchical rankings in directed networks, which also comes with a statistical significance
test for the inferred hierarchy. The model considered also incorporates the assumption
that interactions are more likely to occur between individuals with similar ranks, which is
sometimes the case in real applications, such as sport competitions. The authors compare
SPRING-RANK to a wealth of algorithms, including the above-mentioned SYNC-RANK and
SERIAL-RANK, in terms of an accuracy defined as the fraction of edges whose direction is
consistent with the inferred ranking, and conclude that, on average across various synthetic
and real data sets, SPRING-RANK and SYNC-RANK have the highest accuracy.

Volkovs and Zemel (2014) consider two instances (in the unsupervised and supervised
settings) of the preference aggregation problem where the task is to combine multiple pref-
erences (either in the form of binary comparisons or as score differences) over objects into
a single consensus ranking. Both settings rely on a newly introduced Multinomial Pref-
erence model (MPM) which uses a multinomial generative process to model the observed
preferences. Dalal et al. (2012) consider a multi-objective rank aggregation problem (where
pairwise information is in the form of score differences), and rely on the framework of the
combinatorial Hodge decomposition (Jiang et al., 2011). In particular, they first formulate
pairwise comparisons, and are directed $i \rightarrow j$ towards the winner, where the weight captures the proportion of times $j$ defeated $i$ in case of multiple direct matches. The authors show that under some assumptions on the connectivity of the underlying graph, their algorithm estimates the underlying score vector of the BTL model with $\Omega(npoly(\log n))$ comparisons. Very recently, Agarwal et al. (2018) propose Accelerated Spectral Ranking, a provably faster spectral ranking algorithm in the setting of the multinomial logit (MNL) and BTL models. The authors considered a random walk that has a faster mixing time than the random walks associated with previous algorithms (including Rank-Centrality (Negahban et al., 2017)), along with improved sample complexity bounds (of the order $npoly(\log n)$) for recovery of the MNL and BTL parameters.

Cucuringu (2016) introduced SYNC-RANK, formulating the problem of ranking with incomplete noisy information as an instance of the group synchronization problem over the group $SO(2)$ of planar rotations. SYNC-RANK starts by making the ansatz that items are embedded on the upper-half of the unit circle in the complex plane, and transforms the pairwise input measurements into pairwise angle offsets. It then solves the spectral or SDP relaxation of synchronization, and finally chooses amongst the circular permutations of the ranking induced by the synchronized solution by minimizing the number of upsets. The method was shown to compare favorably with a number of state-of-the-art methods, including RANK-CENTRALITY (Negahban et al., 2012), SERIAL-RANK (Fogel et al., 2016), and a SVD-based ranking algorithm similar to the one we consider in this paper, but without the projection step and its normalized extension, and without any theoretical guarantees (see Remark 3 for details). Along the same lines, and inspired by the group synchronization framework, Fanuel and Suykens (2017) propose a certain deformation of the combinatorial Laplacian, in particular, the so-called dilation Laplacian whose spectrum is leveraged for ranking in directed networks of pairwise comparisons. The method is reported to perform well on both synthetic and real data, and enjoys the property that it can place an emphasis of the top-$k$ items in the ranking.

More recently, De Bacco et al. (2018) proposed SPRING-RANK, an efficient physically-inspired algorithm, based on solving a linear system of equations, for the task of inferring hierarchical rankings in directed networks, which also comes with a statistical significance test for the inferred hierarchy. The model considered also incorporates the assumption that interactions are more likely to occur between individuals with similar ranks, which is sometimes the case in real applications, such as sport competitions. The authors compare SPRING-RANK to a wealth of algorithms, including the above-mentioned SYNC-RANK and SERIAL-RANK, in terms of an accuracy defined as the fraction of edges whose direction is consistent with the inferred ranking, and conclude that, on average across various synthetic and real data sets, SPRING-RANK and SYNC-RANK have the highest accuracy.

Volkovs and Zemel (2014) consider two instances (in the unsupervised and supervised settings) of the preference aggregation problem where the task is to combine multiple preferences (either in the form of binary comparisons or as score differences) over objects into a single consensus ranking. Both settings rely on a newly introduced Multinomial Preference model (MPM) which uses a multinomial generative process to model the observed preferences. Dalal et al. (2012) consider a multi-objective rank aggregation problem (where pairwise information is in the form of score differences), and rely on the framework of the combinatorial Hodge decomposition (Jiang et al., 2011). In particular, they first formulate
the problem of reducing global inconsistencies and then propose techniques for identifying
local observations which can reduce such global inconsistencies.

**Low-rank matrix completion approaches for ranking.** Gleich and Lim (2011) proposed a new method for ranking a set of items from pairwise observations in the form $r_i - r_j$. They first use matrix completion to fill in a partial skew-symmetric matrix, followed by a simple row-sum to recover the rankings. Using standard recovery results from matrix completion literature, it follows that in the absence of noise, and under some coherence assumptions on the score vector $r$, the true ranking can be recovered exactly from a random subset of pairs of size $\Omega(npoly(\log n))$. In case the observations are given as a rating matrix (users and items), they show how to transform it into a pairwise comparison matrix involving score differences so that their method can be applied. Ye et al. (2012) proposed a rank-minimization approach to aggregate the predicted confidence scores of multiple models. The authors cast the score fusion problem as that of finding a shared rank-2 pairwise relationship matrix that renders each of the original score matrices from the different models to be decomposed into the common rank-2 matrix and the sparse error component.

The recent work of Levy et al. (2018) relies on matrix completion in the context of
the rank aggregation problem from a small subset of noisy pairwise comparisons, wherein
the user observes *L* repeated comparisons between the same set of players. The authors
work in the setting of the Bradley-Terry-Luce model (BTL) that assumes that a set of
latent scores (or strengths) underlies all players, and each individual pairwise comparison
is a probabilistic outcome, as a function of the underlying scores. The resulting pipeline is
reported to improve over state-of-the-art in both simulated scenarios and real data. Prior
to that, Kang et al. (2016) also relied on matrix completion and proposed an algorithm for
top-*N* recommender systems. The authors filled in the user-item matrix based on a low-
rank assumption, by considering a nonconvex relaxation, as opposed to the usual nuclear
norm, and argued that it provides a better rank approximation and empirically leads to
an accuracy superior to that of any state-of-the-art algorithm for the top-*N* recommender
problem, on a comprehensive set of experiments on real data sets.

Another relevant line of work is that of Massimino and Davenport (2013), who consider
an adaptation of the one-bit matrix completion framework (Davenport et al., 2014) to
the setting of pairwise comparison matrices. One observes measurements $Y_{ij} = \pm 1$ with
probability $f(M_{ij})$, where $M_{ij} = r_i - r_j$, and $f(M_{ij}) = \mathbb{P}(M_{ij} > 0)$; for example, the authors
consider $f(x) = (1 + e^{-x})^{-1}$. On a related note, Yang and Wakin (2015) considered the
rank aggregation setting, extending the work of Gleich and Lim (2011) to the setup of non-
transitive matrices, for which it does not necessarily hold true that $Y_{i,j} = Y_{i,k} + Y_{k,j}$, $\forall i, j, k$,
as is the case in the special setting when $Y_{i,j} = r_i - r_j$ (and thus $Y = re^T - er^T$), for a
score vector $r$. In particular, Yang and Wakin were interested in modeling and recovering $Y$
itself as opposed to the one-dimensional ranking, and introduced a model for non-transitive
pairwise comparisons $Y_{ij} = r_i a_j - r_j a_i$, where, for example, $r_i$ could denote offensive strength
of player *i*, and $a_j$ defensive strength of player *j*, thus giving $Y_{ij}$ the interpretation of the
anticipated margin of victory for player *i* over player *j*. The authors then propose a low-
rank matrix completion approach based on alternating minimization, along with recovery
guarantees for the estimate $\hat{Y}$ of the form $||\hat{Y} - Y||_F \le \epsilon$, after $\log 1/\epsilon$ iterations.
the problem of reducing global inconsistencies and then propose techniques for identifying
local observations which can reduce such global inconsistencies.

**Low-rank matrix completion approaches for ranking.** Gleich and Lim (2011) proposed a new method for ranking a set of items from pairwise observations in the form $r_i - r_j$. They first use matrix completion to fill in a partial skew-symmetric matrix, followed by a simple row-sum to recover the rankings. Using standard recovery results from matrix completion literature, it follows that in the absence of noise, and under some coherence assumptions on the score vector $r$, the true ranking can be recovered exactly from a random subset of pairs of size $\Omega(npoly(\log n))$. In case the observations are given as a rating matrix (users and items), they show how to transform it into a pairwise comparison matrix involving score differences so that their method can be applied. Ye et al. (2012) proposed a rank-minimization approach to aggregate the predicted confidence scores of multiple models. The authors cast the score fusion problem as that of finding a shared rank-2 pairwise relationship matrix that renders each of the original score matrices from the different models to be decomposed into the common rank-2 matrix and the sparse error component.

The recent work of Levy et al. (2018) relies on matrix completion in the context of
the rank aggregation problem from a small subset of noisy pairwise comparisons, wherein
the user observes *L* repeated comparisons between the same set of players. The authors
work in the setting of the Bradley-Terry-Luce model (BTL) that assumes that a set of
latent scores (or strengths) underlies all players, and each individual pairwise comparison
is a probabilistic outcome, as a function of the underlying scores. The resulting pipeline is
reported to improve over state-of-the-art in both simulated scenarios and real data. Prior
to that, Kang et al. (2016) also relied on matrix completion and proposed an algorithm for
top-*N* recommender systems. The authors filled in the user-item matrix based on a low-
rank assumption, by considering a nonconvex relaxation, as opposed to the usual nuclear
norm, and argued that it provides a better rank approximation and empirically leads to
an accuracy superior to that of any state-of-the-art algorithm for the top-*N* recommender
problem, on a comprehensive set of experiments on real data sets.

Another relevant line of work is that of Massimino and Davenport (2013), who consider
an adaptation of the one-bit matrix completion framework (Davenport et al., 2014) to
the setting of pairwise comparison matrices. One observes measurements $Y_{ij} = \pm 1$ with
probability $f(M_{ij})$, where $M_{ij} = r_i - r_j$, and $f(M_{ij}) = \mathbb{P}(M_{ij} > 0)$; for example, the authors
consider $f(x) = (1 + e^{-x})^{-1}$. On a related note, Yang and Wakin (2015) considered the
rank aggregation setting, extending the work of Gleich and Lim (2011) to the setup of non-
transitive matrices, for which it does not necessarily hold true that $Y_{i,j} = Y_{i,k} + Y_{k,j}$, $\forall i, j, k$,
as is the case in the special setting when $Y_{i,j} = r_i - r_j$ (and thus $Y = re^T - er^T$), for a
score vector $r$. In particular, Yang and Wakin were interested in modeling and recovering $Y$
itself as opposed to the one-dimensional ranking, and introduced a model for non-transitive
pairwise comparisons $Y_{ij} = r_i a_j - r_j a_i$, where, for example, $r_i$ could denote offensive strength
of player $i$, and $a_j$ defensive strength of player $j$, thus giving $Y_{ij}$ the interpretation of the
anticipated margin of victory for player $i$ over player $j$. The authors then propose a low-
rank matrix completion approach based on alternating minimization, along with recovery
guarantees for the estimate $\hat{Y}$ of the form $||\hat{Y} - Y||_F \le \epsilon$, after $\log 1/\epsilon$ iterations.
The recent seminal work of Rajkumar and Agarwal (2016) considered the question of recovering rankings from $O(n \log n)$ comparisons between the items. Their proposed pipeline starts by mapping the input binary entries of the comparison matrix to a low-rank matrix. More specifically, a link function is applied entry-wise to the input comparison matrix, which renders the resulting incomplete matrix to be of low rank. Standard low-rank matrix completion is applied to complete the matrix, followed by the inverse link function applied entry wise, which effectively maps the data back to binary pairwise measurements. Finally, the Copeland ranking procedure (Copeland, 1951) is applied, that ranks items by their Copeland scores, which is effectively given by the number of wins in the comparison matrix after thresholding the entries with respect to 0.5. This approach can be seen as a generalization of the pipeline proposed by Gleich and Lim (2011), in the sense that both lines of work rely on matrix completion followed by row-sum, the main difference being that Rajkumar and Agarwal (2016) uses the link function that maps the binary input matrix to a low-rank one. Borrowing tools from matrix completion, restricted strong convexity and prior work on Random Utility Models, Negahban et al. (2018) propose a nuclear norm regularized optimization problem for learning the parameters of the MultiNomial Logit (MNL) model that best explains the data, typically given as user preference in the form of choices and comparisons. The authors show that the convex relaxation for learning the MNL model is minimax optimal up to a logarithmic factor, by comparing its performance to a fundamental lower bound.

### 3. Problem setup and main results

Our formal setup is as follows. Consider an undirected graph $G = ([n], E)$ and an unknown vector $r \in \mathbb{R}^n$, where $r_i$ is the score associated with node $i$. In particular, $G$ is assumed to be a $G(n, p)$, i.e., the popular Erdős-Rényi random graph model, where edges between vertices are present independently with probability $p$. Moreover, we assume $r_i$ to be bounded uniformly, i.e., $r_i \in [0, M]$ for each $i$. Hence $r_i - r_j \in [-M, M]$ for all $i, j$. Requiring $r_i \ge 0$ is only for convenience and w.l.o.g. Importantly, $M$ is not assumed to be known to the algorithm.

For each $\{i, j\} \in E$, we are given noisy, independent measurements $R_{ij}$ where

$$R_{ij} = \begin{cases} r_i - r_j; & \text{w.p } \eta \\ \sim \mathcal{U}[-M, M]; & \text{w.p } (1 - \eta). \end{cases} \quad (3.1)$$

The parameter $\eta \in [0, 1]$ controls the level of noise; we will denote the noise level explicitly by $\gamma = 1 - \eta$. It is important to note that the parameters $\eta, p$ are not assumed to be known to the algorithm. This model will be referred to as the **Erdős-Rényi Outliers** model, or in short, ERO($n, p, \eta$). It was also considered in previous works in ranking (Cucuringu, 2016) and angular synchronization (Singer, 2011), and is amenable to a theoretical analysis. Alternatives to this noise model include the multiplicative uniform noise model, as considered in Cucuringu (2016). Our goal is two fold - we would like to recover<fnref n="2" /> the score vector $r$, and also the ranking $\pi$ induced by $r$.

<footnote n="2">²Clearly, this is possible only up to a global shift.</footnote>
The recent seminal work of Rajkumar and Agarwal (2016) considered the question of recovering rankings from $O(n \log n)$ comparisons between the items. Their proposed pipeline starts by mapping the input binary entries of the comparison matrix to a low-rank matrix. More specifically, a link function is applied entry-wise to the input comparison matrix, which renders the resulting incomplete matrix to be of low rank. Standard low-rank matrix completion is applied to complete the matrix, followed by the inverse link function applied entry wise, which effectively maps the data back to binary pairwise measurements. Finally, the Copeland ranking procedure (Copeland, 1951) is applied, that ranks items by their Copeland scores, which is effectively given by the number of wins in the comparison matrix after thresholding the entries with respect to 0.5. This approach can be seen as a generalization of the pipeline proposed by Gleich and Lim (2011), in the sense that both lines of work rely on matrix completion followed by row-sum, the main difference being that Rajkumar and Agarwal (2016) uses the link function that maps the binary input matrix to a low-rank one. Borrowing tools from matrix completion, restricted strong convexity and prior work on Random Utility Models, Negahban et al. (2018) propose a nuclear norm regularized optimization problem for learning the parameters of the MultiNomial Logit (MNL) model that best explains the data, typically given as user preference in the form of choices and comparisons. The authors show that the convex relaxation for learning the MNL model is minimax optimal up to a logarithmic factor, by comparing its performance to a fundamental lower bound.

### 3. Problem setup and main results

Our formal setup is as follows. Consider an undirected graph $G = ([n], E)$ and an unknown vector $r \in \mathbb{R}^n$, where $r_i$ is the score associated with node $i$. In particular, $G$ is assumed to be a $G(n, p)$, i.e., the popular Erdős-Rényi random graph model, where edges between vertices are present independently with probability $p$. Moreover, we assume $r_i$ to be bounded uniformly, i.e., $r_i \in [0, M]$ for each $i$. Hence $r_i - r_j \in [-M, M]$ for all $i, j$. Requiring $r_i \ge 0$ is only for convenience and w.l.o.g. Importantly, $M$ is not assumed to be known to the algorithm.

For each $\{i, j\} \in E$, we are given noisy, independent measurements $R_{ij}$ where

$$R_{ij} = \begin{cases} r_i - r_j; & \text{w.p } \eta \\ \sim \mathcal{U}[-M, M]; & \text{w.p } (1 - \eta). \end{cases} \quad (3.1)$$

The parameter $\eta \in [0, 1]$ controls the level of noise; we will denote the noise level explicitly by $\gamma = 1 - \eta$. It is important to note that the parameters $\eta, p$ are not assumed to be known to the algorithm. This model will be referred to as the **Erdős-Rényi Outliers** model, or in short, ERO($n, p, \eta$). It was also considered in previous works in ranking (Cucuringu, 2016) and angular synchronization (Singer, 2011), and is amenable to a theoretical analysis. Alternatives to this noise model include the multiplicative uniform noise model, as considered in Cucuringu (2016). Our goal is two fold - we would like to recover<fnref n="2" /> the score vector $r$, and also the ranking $\pi$ induced by $r$.

<footnote n="2">²Clearly, this is possible only up to a global shift.</footnote>
**Remark 1** The above model is only for the purpose of theoretically analyzing the statistical performance of our methods. Other statistical models could also be considered of course, and we will see that our methods are completely model-independent.

## 3.1 Main idea: SVD-based spectral algorithm

We start by forming the measurement matrix $H \in \mathbb{R}^{n \times n}$, where

* $H_{ii} = 0, \forall i = 1, \dots, n,$

* $H_{ij} = R_{ij}$ and $H_{ji} = -R_{ij}$, if $\{i, j\} \in E$, and

* $H_{ij} = 0$, if $\{i, j\} \notin E$.

If $G$ is the complete graph and the measurements are noise free, it holds true that $H = re^T - er^T$, which is a rank 2 skew-symmetric matrix. Denoting $\alpha = \frac{r^T e}{n}$, one can verify that the two non-zero left singular vectors are $u_1 = e/\sqrt{n}$, $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ with equal non-zero singular values $\sigma_1 = \sigma_2 = \|r - \alpha e\|_2\sqrt{n}$ (see Lemma 11). Therefore, given any orthonormal basis for $\{u_1, u_2\}$, we can simply find a vector orthonormal to $e/\sqrt{n}$ and which lies in $\text{span}\{u_1, u_2\}$; this will give us candidate solutions $\pm \frac{r-\alpha e}{\|r-\alpha e\|_2}$. Multiplying these candidates by $\sigma_1/\sqrt{n}$ recovers the scale information of $r$, while the sign ambiguity is resolved by selecting the candidate which is most consistent with the measurements. If one is interested only in ranking the items, then there is of course no need to estimate the scale above.

If $G$ is not complete (thus there are missing edges) and the measurements are noisy, then $H$ will typically not be rank 2, but as will see shortly, it can be thought of as a perturbation of the rank-2 matrix $\eta p(re^T - er^T)$. Proceeding similarly to the noiseless case, we can find the top two singular vectors of $H$ (denoted by $\hat{u}_1, \hat{u}_2$), project $u_1 = e/\sqrt{n}$ on to $\text{span}\{\hat{u}_1, \hat{u}_2\}$ to obtain $\bar{u}_1$, and then find a unit vector in $\text{span}\{\hat{u}_1, \hat{u}_2\}$ which is orthonormal to $\bar{u}_1$ (call this $\tilde{u}_2$). If the noise level is not too large and we have sufficiently many edges in $G$, then one can imagine that $\tilde{u}_2$ will be close (up to a sign flip) to $u_2$, and $\hat{\sigma}_i \approx \eta p \sigma_i$. Hence we can recover the ranking from $\tilde{u}_2$ after resolving the sign ambiguity (as mentioned earlier). This also implies that the centered version of $\frac{\sigma_1 \tilde{u}_2}{\eta p \sqrt{n}}$ will be close (up to a sign flip) to $r - \alpha e$. But since we do not know $\eta, p$, we will resort to other data-driven approaches for recovering the scale parameter, see Section 3.2.

The above approach is outlined formally as Algorithm 1, namely SVD-RS (SVD-Ranking and Synchronization) which is a spectral method for recovering the ranks and scores (up to a global shift) of a collection of $n$ items. Additionally, we also consider a “normalized” version of SVD-RS, wherein $H$ is replaced by $H_{ss} = \bar{D}^{-1/2} H \bar{D}^{-1/2}$ with $\bar{D}$ being a diagonal matrix and $\bar{D}_{ii} = \sum_{j=1}^n |H_{ij}|$. Clearly, $H_{ss}$ is also skew symmetric. Such a normalization step is particularly useful when the degree distribution is skewed, and is commonly employed in other problem domains involving spectral methods, such as clustering (see for eg. Kunegis et al. (2010)). The ensuing algorithm, namely SVD-NRS (SVD-Normalized Ranking and Synchronization), is outlined as Algorithm 2.

**Remark 2** Algorithms 1 and 2 are model-independent, and only require as input (a) the measurement graph $G$, and (b) the corresponding pairwise measurements $R_{ij}$.
**Remark 1** The above model is only for the purpose of theoretically analyzing the statistical performance of our methods. Other statistical models could also be considered of course, and we will see that our methods are completely model-independent.

## 3.1 Main idea: SVD-based spectral algorithm

We start by forming the measurement matrix $H \in \mathbb{R}^{n \times n}$, where

* $H_{ii} = 0, \forall i = 1, \dots, n,$

* $H_{ij} = R_{ij}$ and $H_{ji} = -R_{ij}$, if $\{i, j\} \in E$, and

* $H_{ij} = 0$, if $\{i, j\} \notin E$.

If $G$ is the complete graph and the measurements are noise free, it holds true that $H = re^T - er^T$, which is a rank 2 skew-symmetric matrix. Denoting $\alpha = \frac{r^T e}{n}$, one can verify that the two non-zero left singular vectors are $u_1 = e/\sqrt{n}$, $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ with equal non-zero singular values $\sigma_1 = \sigma_2 = \|r - \alpha e\|_2\sqrt{n}$ (see Lemma 11). Therefore, given any orthonormal basis for $\{u_1, u_2\}$, we can simply find a vector orthonormal to $e/\sqrt{n}$ and which lies in $\text{span}\{u_1, u_2\}$; this will give us candidate solutions $\pm \frac{r-\alpha e}{\|r-\alpha e\|_2}$. Multiplying these candidates by $\sigma_1/\sqrt{n}$ recovers the scale information of $r$, while the sign ambiguity is resolved by selecting the candidate which is most consistent with the measurements. If one is interested only in ranking the items, then there is of course no need to estimate the scale above.

If $G$ is not complete (thus there are missing edges) and the measurements are noisy, then $H$ will typically not be rank 2, but as will see shortly, it can be thought of as a perturbation of the rank-2 matrix $\eta p(re^T - er^T)$. Proceeding similarly to the noiseless case, we can find the top two singular vectors of $H$ (denoted by $\hat{u}_1, \hat{u}_2$), project $u_1 = e/\sqrt{n}$ on to $\text{span}\{\hat{u}_1, \hat{u}_2\}$ to obtain $\bar{u}_1$, and then find a unit vector in $\text{span}\{\hat{u}_1, \hat{u}_2\}$ which is orthonormal to $\bar{u}_1$ (call this $\tilde{u}_2$). If the noise level is not too large and we have sufficiently many edges in $G$, then one can imagine that $\tilde{u}_2$ will be close (up to a sign flip) to $u_2$, and $\hat{\sigma}_i \approx \eta p \sigma_i$. Hence we can recover the ranking from $\tilde{u}_2$ after resolving the sign ambiguity (as mentioned earlier). This also implies that the centered version of $\frac{\sigma_1 \tilde{u}_2}{\eta p \sqrt{n}}$ will be close (up to a sign flip) to $r - \alpha e$. But since we do not know $\eta, p$, we will resort to other data-driven approaches for recovering the scale parameter, see Section 3.2.

The above approach is outlined formally as Algorithm 1, namely SVD-RS (SVD-Ranking and Synchronization) which is a spectral method for recovering the ranks and scores (up to a global shift) of a collection of $n$ items. Additionally, we also consider a “normalized” version of SVD-RS, wherein $H$ is replaced by $H_{ss} = \bar{D}^{-1/2} H \bar{D}^{-1/2}$ with $\bar{D}$ being a diagonal matrix and $\bar{D}_{ii} = \sum_{j=1}^n |H_{ij}|$. Clearly, $H_{ss}$ is also skew symmetric. Such a normalization step is particularly useful when the degree distribution is skewed, and is commonly employed in other problem domains involving spectral methods, such as clustering (see for eg. Kunegis et al. (2010)). The ensuing algorithm, namely SVD-NRS (SVD-Normalized Ranking and Synchronization), is outlined as Algorithm 2.

**Remark 2** Algorithms 1 and 2 are model-independent, and only require as input (a) the measurement graph $G$, and (b) the corresponding pairwise measurements $R_{ij}$.
**Algorithm 1 SVD-RS**

1: **Input:** Measurement graph $G = ([n], E)$ and pairwise measurements $R_{ij}$ for $\{i, j\} \in E$.

2: **Output:** Rank estimates: $\hat{\pi}$ and score estimates $\hat{r} \in \mathbb{R}^n$.

3: Form measurement matrix $H \in \mathbb{R}^{n \times n}$ using $R_{ij}$ as outlined in Section 3.1.

4: Find the top two left singular vectors (resp. singular values) of $H$, namely $\hat{u}_1, \hat{u}_2$ (resp. $\hat{\sigma}_1, \hat{\sigma}_2$).

5: Obtain vector $\bar{u}_1$ as the orthogonal projection of $u_1 = e/\sqrt{n}$ on to span $\{\hat{u}_1, \hat{u}_2\}$.

6: Obtain a unit vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ such that $\tilde{u}_2 \perp \bar{u}_1$.

7: **Rank recovery:** Obtain ranking $\tilde{\pi}$ induced by $\tilde{u}_2$ (up to global sign ambiguity) and reconcile its global sign by minimizing the number of upsets. Output ranking estimate $\tilde{\pi}$.

8: **Score recovery:** Use $\tilde{u}_2$, $H$ to recover the scale $\tau \in \mathbb{R}$ as in Section 3.2. Output $\hat{r} = \tau \tilde{u}_2 - \frac{e^T(\tau \tilde{u}_2)}{n}e$.

**Algorithm 2 SVD-NRS**

1: **Input:** Measurement graph $G = ([n], E)$ and pairwise measurements $R_{ij}$ for $\{i, j\} \in E$.

2: **Output:** Rank estimates: $\hat{\pi}$ and score estimates $\hat{r} \in \mathbb{R}^n$.

3: Form $H_{ss} = \bar{D}^{-1/2}H\bar{D}^{-1/2}$ with $H$ formed using $R_{ij}$ as outlined in Section 3.1. $\bar{D}$ is a diagonal matrix with $\bar{D}_{ii} = \sum_{j=1}^n |H_{ij}|$.

4: Find the top two left singular vectors (resp. singular values) of $H_{ss}$, namely $\hat{u}_1, \hat{u}_2$ (resp. $\hat{\sigma}_1, \hat{\sigma}_2$).

5: Obtain vector $\bar{u}_1$ as the orthogonal projection of $u_1 = \frac{\bar{D}^{-1/2}e}{\|\bar{D}^{-1/2}e\|_2}$ on to span $\{\hat{u}_1, \hat{u}_2\}$.

6: Obtain a unit vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ such that $\tilde{u}_2 \perp \bar{u}_1$.

7: **Rank recovery:** Obtain ranking $\tilde{\pi}$ induced by $\tilde{D}^{1/2}\tilde{u}_2$ (up to global sign ambiguity) and reconcile its global sign by minimizing the number of upsets. Output ranking estimate $\tilde{\pi}$.

8: **Score recovery:** Use $\tilde{u}_2$, $H$ to recover the scale $\tau \in \mathbb{R}$ as in Section 3.2. Output $\hat{r} = \tau \bar{D}^{1/2}\tilde{u}_2 - \frac{e^T(\tau \bar{D}^{1/2}\tilde{u}_2)}{n}e$.

**Remark 3** Note that the SVD-based algorithm introduced in Cucuringu (2016) considers the four possible rankings induced by the singular vectors $\{\pm\hat{u}_1, \pm\hat{u}_2\}$, and chooses the one that minimizes the number of upsets. In contrast, SVD-RS first computes $\bar{u}_1$ as the orthogonal projection of the all ones vector $e$ on to span $\{\hat{u}_1, \hat{u}_2\}$, finds a vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ perpendicular to $\bar{u}_1$, and finally extracts the ranking induced by $\tilde{u}_2$. Furthermore, unlike Cucuringu (2016), SVD-NRS introduces an additional normalization step of the measurement matrix $H$ to alleviate potential issues arising from skewed degree distributions. Such a normalization step is fairly common for spectral methods in general. For example, this operator was considered in Cucuringu et al. (2012a) and Cucuringu et al. (2012b) in the context of angular synchronization (Singer, 2011) and the graph realization problem, and Singer and Wu (2012) who introduced Vector Diffusion Maps for nonlinear dimensionality reduction and explored the interplay with the Connection-Laplacian operator for vector fields over manifolds.
**Algorithm 1 SVD-RS**

1: **Input:** Measurement graph $G = ([n], E)$ and pairwise measurements $R_{ij}$ for $\{i, j\} \in E$.

2: **Output:** Rank estimates: $\hat{\pi}$ and score estimates $\hat{r} \in \mathbb{R}^n$.

3: Form measurement matrix $H \in \mathbb{R}^{n \times n}$ using $R_{ij}$ as outlined in Section 3.1.

4: Find the top two left singular vectors (resp. singular values) of $H$, namely $\hat{u}_1, \hat{u}_2$ (resp. $\hat{\sigma}_1, \hat{\sigma}_2$).

5: Obtain vector $\bar{u}_1$ as the orthogonal projection of $u_1 = e/\sqrt{n}$ on to span $\{\hat{u}_1, \hat{u}_2\}$.

6: Obtain a unit vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ such that $\tilde{u}_2 \perp \bar{u}_1$.

7: **Rank recovery:** Obtain ranking $\tilde{\pi}$ induced by $\tilde{u}_2$ (up to global sign ambiguity) and reconcile its global sign by minimizing the number of upsets. Output ranking estimate $\tilde{\pi}$.

8: **Score recovery:** Use $\tilde{u}_2$, $H$ to recover the scale $\tau \in \mathbb{R}$ as in Section 3.2. Output $\hat{r} = \tau \tilde{u}_2 - \frac{e^T(\tau \tilde{u}_2)}{n}e$.

**Algorithm 2 SVD-NRS**

1: **Input:** Measurement graph $G = ([n], E)$ and pairwise measurements $R_{ij}$ for $\{i, j\} \in E$.

2: **Output:** Rank estimates: $\hat{\pi}$ and score estimates $\hat{r} \in \mathbb{R}^n$.

3: Form $H_{ss} = \bar{D}^{-1/2}H\bar{D}^{-1/2}$ with $H$ formed using $R_{ij}$ as outlined in Section 3.1. $\bar{D}$ is a diagonal matrix with $\bar{D}_{ii} = \sum_{j=1}^n |H_{ij}|$.

4: Find the top two left singular vectors (resp. singular values) of $H_{ss}$, namely $\hat{u}_1, \hat{u}_2$ (resp. $\hat{\sigma}_1, \hat{\sigma}_2$).

5: Obtain vector $\bar{u}_1$ as the orthogonal projection of $u_1 = \frac{\bar{D}^{-1/2}e}{\|\bar{D}^{-1/2}e\|_2}$ on to span $\{\hat{u}_1, \hat{u}_2\}$.

6: Obtain a unit vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ such that $\tilde{u}_2 \perp \bar{u}_1$.

7: **Rank recovery:** Obtain ranking $\tilde{\pi}$ induced by $\tilde{D}^{1/2}\tilde{u}_2$ (up to global sign ambiguity) and reconcile its global sign by minimizing the number of upsets. Output ranking estimate $\tilde{\pi}$.

8: **Score recovery:** Use $\tilde{u}_2$, $H$ to recover the scale $\tau \in \mathbb{R}$ as in Section 3.2. Output $\hat{r} = \tau \bar{D}^{1/2}\tilde{u}_2 - \frac{e^T(\tau \bar{D}^{1/2}\tilde{u}_2)}{n}e$.

**Remark 3** Note that the SVD-based algorithm introduced in Cucuringu (2016) considers the four possible rankings induced by the singular vectors $\{\pm\hat{u}_1, \pm\hat{u}_2\}$, and chooses the one that minimizes the number of upsets. In contrast, SVD-RS first computes $\bar{u}_1$ as the orthogonal projection of the all ones vector $e$ on to span $\{\hat{u}_1, \hat{u}_2\}$, finds a vector $\tilde{u}_2 \in \text{span}\{\hat{u}_1, \hat{u}_2\}$ perpendicular to $\bar{u}_1$, and finally extracts the ranking induced by $\tilde{u}_2$. Furthermore, unlike Cucuringu (2016), SVD-NRS introduces an additional normalization step of the measurement matrix $H$ to alleviate potential issues arising from skewed degree distributions. Such a normalization step is fairly common for spectral methods in general. For example, this operator was considered in Cucuringu et al. (2012a) and Cucuringu et al. (2012b) in the context of angular synchronization (Singer, 2011) and the graph realization problem, and Singer and Wu (2012) who introduced Vector Diffusion Maps for nonlinear dimensionality reduction and explored the interplay with the Connection-Laplacian operator for vector fields over manifolds.
**Computational complexity.** Our proposed SVD-based methods rely on computing the top two singular vectors of the input matrix which in practice can be done very efficiently. In general, the leading singular values and singular vectors can be computed using iterative techniques at a typical cost of $O(pn^2)$, using for example, a simple power method (since all iterations only require a matrix-vector product at a cost of $pn^2$, where $p$ is the sampling probability). In the sparse setting, the computational cost is essentially linear in the number of nonzero entries in the matrix, thus making the approach scalable to large measurement graphs. The data-driven approach for scale recovery will also require $O(pn^2)$ cost, as will be seen in the next Section 3.2.

**Remark 4** The computational complexity of our proposed SVD-based algorithms is on par with that of other competitive methods, detailed and compared against in the numerical experiments Section 7. For example, the Serial-Rank approach (Fogel et al., 2014) first computes a shifted version of the similarity matrix $HH^T$, then considers the graph Laplacian, and computes its Fiedler eigenvector, much like our SVD-based approach, which can be done at a typical cost of $O(pn^2)$. The Row-Sum method computes the sum of the entries in each row of $H$, which has complexity $O(pn^2)$, though Gleich and Lim (2011) first rely on a matrix completion step to fill in a partial skew-symmetric matrix, before computing the row sums. However, performing a low-rank matrix completion step can be computationally expensive, depending on the approach used; in general, alternating minimization has been proven empirically to be one of the most accurate and efficient methods for matrix completion (Jain et al., 2013). In the numerical experiments reported in Section 7, we ran into computational issues when running matrix completion for values of $n$ significantly larger than 1000, even for very sparse graphs (eg, $p = 0.01$).

For the Least-Squared-based ranking approach, also considered in Cucuringu (2016), the setup is the following. Each (potentially noisy) pairwise comparison $r_i - r_j$ contributes with a row to a design matrix $T$ of size $m \times n$, where $m \approx pn^2$ denotes the number of edges in the graph $G$. Each row of $T$ has only two non-zero elements, in particular, $a+1$ in column $i$, and $a-1$ in column $j$. To estimate the score vector $s$ of length $n$, we solve the linear system $Ts = b$ in the least-squares sense, where $b$ is a column vector with the outputs of the pairwise comparisons. One potential approach for solving such linear least-squares problems of the form $Ts = b$ is to use conjugate gradient iterations applied to the normal equations $T^T Ts = T^T b$ (which can be achieved without explicitly performing the expensive calculation of the matrix $T^T T$). Note that the rate of convergence of the gradient iterations is dictated by the condition number $\kappa$ of the matrix $T^T T$, and the number of iterations required for convergence is $O(\sqrt{\kappa})$ (Trefethen and Bau, 1997). For matrices that are sparse or have exploitable structure (like the sparse matrix $T$ with only two nonzero entries per row), each conjugate gradient iteration has complexity as low as $O(m)$.

## 3.2 Recovering the global scale: data-driven procedures

To recover the scaling factor $\tau \in \mathbb{R}$, we consider two possible approaches: a first one based on a median estimator, and a second one that relies on a least-squares estimator. We briefly review both methods, and remark that we only consider the median-based estimator throughout the rest of the paper, due to its additional robustness. We compare results on
**Computational complexity.** Our proposed SVD-based methods rely on computing the top two singular vectors of the input matrix which in practice can be done very efficiently. In general, the leading singular values and singular vectors can be computed using iterative techniques at a typical cost of $O(pn^2)$, using for example, a simple power method (since all iterations only require a matrix-vector product at a cost of $pn^2$, where $p$ is the sampling probability). In the sparse setting, the computational cost is essentially linear in the number of nonzero entries in the matrix, thus making the approach scalable to large measurement graphs. The data-driven approach for scale recovery will also require $O(pn^2)$ cost, as will be seen in the next Section 3.2.

**Remark 4** The computational complexity of our proposed SVD-based algorithms is on par with that of other competitive methods, detailed and compared against in the numerical experiments Section 7. For example, the Serial-Rank approach (Fogel et al., 2014) first computes a shifted version of the similarity matrix $HH^T$, then considers the graph Laplacian, and computes its Fiedler eigenvector, much like our SVD-based approach, which can be done at a typical cost of $O(pn^2)$. The Row-Sum method computes the sum of the entries in each row of $H$, which has complexity $O(pn^2)$, though Gleich and Lim (2011) first rely on a matrix completion step to fill in a partial skew-symmetric matrix, before computing the row sums. However, performing a low-rank matrix completion step can be computationally expensive, depending on the approach used; in general, alternating minimization has been proven empirically to be one of the most accurate and efficient methods for matrix completion (Jain et al., 2013). In the numerical experiments reported in Section 7, we ran into computational issues when running matrix completion for values of $n$ significantly larger than 1000, even for very sparse graphs (eg, $p = 0.01$).

For the Least-Squared-based ranking approach, also considered in Cucuringu (2016), the setup is the following. Each (potentially noisy) pairwise comparison $r_i - r_j$ contributes with a row to a design matrix $T$ of size $m \times n$, where $m \approx pn^2$ denotes the number of edges in the graph $G$. Each row of $T$ has only two non-zero elements, in particular, $a+1$ in column $i$, and $a-1$ in column $j$. To estimate the score vector $s$ of length $n$, we solve the linear system $Ts = b$ in the least-squares sense, where $b$ is a column vector with the outputs of the pairwise comparisons. One potential approach for solving such linear least-squares problems of the form $Ts = b$ is to use conjugate gradient iterations applied to the normal equations $T^T Ts = T^T b$ (which can be achieved without explicitly performing the expensive calculation of the matrix $T^T T$). Note that the rate of convergence of the gradient iterations is dictated by the condition number $\kappa$ of the matrix $T^T T$, and the number of iterations required for convergence is $O(\sqrt{\kappa})$ (Trefethen and Bau, 1997). For matrices that are sparse or have exploitable structure (like the sparse matrix $T$ with only two nonzero entries per row), each conjugate gradient iteration has complexity as low as $O(m)$.

## 3.2 Recovering the global scale: data-driven procedures

To recover the scaling factor $\tau \in \mathbb{R}$, we consider two possible approaches: a first one based on a median estimator, and a second one that relies on a least-squares estimator. We briefly review both methods, and remark that we only consider the median-based estimator throughout the rest of the paper, due to its additional robustness. We compare results on
two problem instances of varying levels of noise, and also compare with the ground truth, as detailed below.

<figure><img src="image_3.png" /><figcaption>(a) Histogram of pairwise ratios (b) Histogram of pairwise ratios (c) Regression-based estimator in absolute value (on a log scale), in absolute value, for $\gamma = 0.02$. of the scaling factor, for $\gamma = 0.02$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(d) Histogram of pairwise ratios (e) Histogram of pairwise ratios (f) Regression-based estimator in absolute value (on a log scale), (zoom in [-15,15]), with the median estimator, for $\gamma = 0.30$. of the scaling factor, for $\gamma = 0.30$.</figcaption></figure>

Figure 1: Two instances (top: low-noise regime $\gamma = 0.02$, bottom: high-noise regime $\gamma = 0.30$) showcasing the median and least-squares based estimators for scale recovery (shown in blue), along with their respective ground truth values (shown in red). Note that the bottom instance is also reconciliating the global sign. Scores are gamma distributed, with $n = 500$, and edge density $p = 0.25$.

To begin with, we define the entrywise ratio of offsets/comparisons as

$$ \Pi_{i,j} = \begin{cases} \frac{H_{i,j}}{S_{i,j}}, & \{i, j\} \in E \\ 0, & \text{otherwise,} \end{cases} \qquad (3.2) $$

where $H_{i,j}$ denotes the pairwise comparison measurement initially available, and $S_{i,j}$ denotes the recovered/estimated pairwise offset given by $S_{i,j} = s_i - s_j$, where $s_i$ denotes the recovered score by either SVD-RS and SVD-NRS. In particular, $s_i = (\tilde{u}_2)_i$ for SVD-RS, and $s_i = (\bar{D}^{1/2}\tilde{u}_2)_i$ for SVD-NRS. This scale recovery method is broadly applicable to other ranking methods we compare against, and we further detail this in the experiments section. If $|E| = m$, we next consider the vector $\rho$ of length $m$ with the pairwise ratios, i.e. nonzero entries in the matrix $\Pi_{i,j}$ in (3.2) corresponding to the edges in the measurement graph $G$.
two problem instances of varying levels of noise, and also compare with the ground truth, as detailed below.

<figure><img src="image_3.png" /><figcaption>(a) Histogram of pairwise ratios (b) Histogram of pairwise ratios (c) Regression-based estimator in absolute value (on a log scale), in absolute value, for $\gamma = 0.02$. of the scaling factor, for $\gamma = 0.02$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(d) Histogram of pairwise ratios (e) Histogram of pairwise ratios (f) Regression-based estimator in absolute value (on a log scale), (zoom in [-15,15]), with the median estimator, for $\gamma = 0.30$. of the scaling factor, for $\gamma = 0.30$.</figcaption></figure>

Figure 1: Two instances (top: low-noise regime $\gamma = 0.02$, bottom: high-noise regime $\gamma = 0.30$) showcasing the median and least-squares based estimators for scale recovery (shown in blue), along with their respective ground truth values (shown in red). Note that the bottom instance is also reconciliating the global sign. Scores are gamma distributed, with $n = 500$, and edge density $p = 0.25$.

To begin with, we define the entrywise ratio of offsets/comparisons as

$$ \Pi_{i,j} = \begin{cases} \frac{H_{i,j}}{S_{i,j}}, & \{i,j\} \in E \\ 0, & \text{otherwise,} \end{cases} \qquad (3.2) $$

where $H_{i,j}$ denotes the pairwise comparison measurement initially available, and $S_{i,j}$ denotes the recovered/estimated pairwise offset given by $S_{i,j} = s_i - s_j$, where $s_i$ denotes the recovered score by either SVD-RS and SVD-NRS. In particular, $s_i = (\tilde{u}_2)_i$ for SVD-RS, and $s_i = (\bar{D}^{1/2}\tilde{u}_2)_i$ for SVD-NRS. This scale recovery method is broadly applicable to other ranking methods we compare against, and we further detail this in the experiments section. If $|E| = m$, we next consider the vector $\rho$ of length $m$ with the pairwise ratios, i.e. nonzero entries in the matrix $\Pi_{i,j}$ in (3.2) corresponding to the edges in the measurement graph $G$.
The median estimator is given by

$$ \hat{\tau}^{Median} = \text{median}(\rho_1, \dots, \rho_m), \qquad (3.3) $$

as illustrated in the middle plot in Figure 1. From a computational perspective, the median-based estimator only requires forming the pairwise ratios for the $m$ edges in the graph, and computing the median element can be performed in linear time.

An alternative approach is to consider the least-squares formulation

$$ \hat{\tau}^{LS} = \arg \min_{\tau \in \mathbb{R}} \sum_{\{i,j\} \in E} (H_{i,j} - \tau S_{i,j})^2 = \frac{\sum_{\{i,j\} \in E} H_{i,j}}{\sum_{\{i,j\} \in E} S_{i,j}}, \qquad (3.4) $$

as shown in the rightmost column of Figure 1. In the low-noise regime, the least-squares estimator is close to the median-based estimator, and the discrepancy increases as the noise levels becomes larger. As shown in the left plot of Figure 1, which contains a histogram of the nonzero entries $\Pi_{i,j} \neq 0$, $\{i, j\} \in E$ on a log-scale, gross outliers can affect the regression formulation (3.4). As expected, we have observed in our simulations that the median estimator (3.3) leads to lower recovery errors compared to the regression-based estimator. In Figure 1, we show in red the recovered scaling when employing the ground truth data, in particular

$$ \Pi_{i,j}^{(\text{ground truth})} = \begin{cases} \frac{r_i - r_j}{S_{i,j}}, & \{i, j\} \in E \\ 0; & \text{otherwise,} \end{cases} $$

Similarly to (3.3) and (3.4), we compute $\tau^{Median}$, respectively, $\tau^{LS}$, construed as the ground truth, and shown in red in the middle and right plots of Figure 1. Note that, in both instances, the median-based estimator is closer to the ground truth than the least-squares based estimator. In the low-noise regime $\gamma = 0.02$ (top plot of Figure 1), $\hat{\tau}^{Median}$ is within 0.23% from its ground true counterpart, while $\hat{\tau}^{LS}$ is within 1.7%. However, the performance gap widens as we increase the noise level. In the high noise regime $\gamma = 0.30$ (bottom plot of Figure 1), $\hat{\tau}^{Median}$ is off by 4.5% with respect to the ground truth, while $\hat{\tau}^{LS}$ is off by 27.1% from its ground truth counterpart. Note that in the latter problem instance, the recovered scaling is negative, emphasizing that the global sign has been reconciled. Due to its increased robustness, for the remainder of the numerical experiments in this paper, we henceforth rely only on the median-based estimator to perform the scale recovery step.

## 3.3 Summary of main theoretical results

Broadly speaking, the bulk of our analysis revolves around bounding the distance between the unit norm vectors $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ and $\tilde{u}_2$ (up to a global sign), the latter obtained in SVD-RS (Step 6) and SVD-NRS (Step 6). For SVD-RS, these error estimates are stated in Theorem 5 (resp. Theorem 6) for the $\ell_2$ (resp. $\ell_\infty$) norm, while Theorem 9 states the corresponding $\ell_2$ error bound for SVD-NRS. Below, we outline what this implies for rank and score recovery.

1. Theorem 7 gives guarantees for rank recovery (up to a global order reversal) for SVD-RS in terms of the maximum displacement error (defined in (4.6)) between
The median estimator is given by

$$ \hat{\tau}^{Median} = \text{median}(\rho_1, \dots, \rho_m), \qquad (3.3) $$

as illustrated in the middle plot in Figure 1. From a computational perspective, the median-based estimator only requires forming the pairwise ratios for the $m$ edges in the graph, and computing the median element can be performed in linear time.

An alternative approach is to consider the least-squares formulation

$$ \hat{\tau}^{LS} = \arg\min_{\tau \in \mathbb{R}} \sum_{\{i,j\} \in E} (H_{i,j} - \tau S_{i,j})^2 = \frac{\sum_{\{i,j\} \in E} H_{i,j}}{\sum_{\{i,j\} \in E} S_{i,j}}, \qquad (3.4) $$

as shown in the rightmost column of Figure 1. In the low-noise regime, the least-squares estimator is close to the median-based estimator, and the discrepancy increases as the noise levels becomes larger. As shown in the left plot of Figure 1, which contains a histogram of the nonzero entries $\Pi_{i,j} \neq 0$, $\{i, j\} \in E$ on a log-scale, gross outliers can affect the regression formulation (3.4). As expected, we have observed in our simulations that the median estimator (3.3) leads to lower recovery errors compared to the regression-based estimator. In Figure 1, we show in red the recovered scaling when employing the ground truth data, in particular

$$ \Pi_{i,j}^{(\text{ground truth})} = \begin{cases} \frac{r_i - r_j}{S_{i,j}}, & \{i, j\} \in E \\ 0; & \text{otherwise,} \end{cases} $$

Similarly to (3.3) and (3.4), we compute $\tau^{Median}$, respectively, $\tau^{LS}$, construed as the ground truth, and shown in red in the middle and right plots of Figure 1. Note that, in both instances, the median-based estimator is closer to the ground truth than the least-squares based estimator. In the low-noise regime $\gamma = 0.02$ (top plot of Figure 1), $\hat{\tau}^{Median}$ is within 0.23% from its ground true counterpart, while $\hat{\tau}^{LS}$ is within 1.7%. However, the performance gap widens as we increase the noise level. In the high noise regime $\gamma = 0.30$ (bottom plot of Figure 1), $\hat{\tau}^{Median}$ is off by 4.5% with respect to the ground truth, while $\hat{\tau}^{LS}$ is off by 27.1% from its ground truth counterpart. Note that in the latter problem instance, the recovered scaling is negative, emphasizing that the global sign has been reconciled. Due to its increased robustness, for the remainder of the numerical experiments in this paper, we henceforth rely only on the median-based estimator to perform the scale recovery step.

## 3.3 Summary of main theoretical results

Broadly speaking, the bulk of our analysis revolves around bounding the distance between the unit norm vectors $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ and $\tilde{u}_2$ (up to a global sign), the latter obtained in SVD-RS (Step 6) and SVD-NRS (Step 6). For SVD-RS, these error estimates are stated in Theorem 5 (resp. Theorem 6) for the $\ell_2$ (resp. $\ell_\infty$) norm, while Theorem 9 states the corresponding $\ell_2$ error bound for SVD-NRS. Below, we outline what this implies for rank and score recovery.

1. Theorem 7 gives guarantees for rank recovery (up to a global order reversal) for SVD-RS in terms of the maximum displacement error (defined in (4.6)) between
the ranking $\tilde{\pi}$ obtained in Step 7 of Algorithm 1, and the ground truth ranking $\pi$. Denoting $\alpha = \frac{r^T e}{n}$ to be the average score, it states that if

$$ \|r - \alpha e\|_2 \gtrsim \frac{M}{\eta \sqrt{p}}, \quad p \gtrsim \frac{\log n}{n}, \qquad (3.5) $$

and if $n$ is large enough, then $\|\tilde{\pi} - \pi\|_{\infty} \lesssim \frac{\|r - \alpha e\|_2}{\rho} \Upsilon(n, M, \eta, p, \varepsilon, r)$ w.h.p. Here $\Upsilon(n, M, \eta, p, \varepsilon, r)$ is the bound on the $\ell_{\infty}$ error between $\tilde{u}_2$ and $u_2$ (up to a global sign) in Theorem 6. For concreteness, if $r_i = i$, we show in Example 3 that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{(\log n)^{2/3}}{n^{2/3}})$ and $n$ is large enough, then $\|\tilde{\pi} - \pi\|_{\infty} \lesssim \delta \frac{n}{\log n}$ w.h.p. Note that this implies that $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice.

2. Theorem 8 provides $\ell_2$ and $\ell_{\infty}$ bounds (up to a global sign and shift) for SVD-RS for the score vector $r$ and a vector $\hat{r}$ (derived from $\tilde{u}_2$) where

$$ \tilde{w} = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}} \tilde{u}_2; \quad \tilde{r} := \tilde{w} - \frac{e^T \tilde{w}}{n} e. \qquad (3.6) $$

If one knew the product $\eta p$, then (3.6) corresponds to choosing the scale $\tau = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}}$ in Step 8.

(a) If the conditions in (3.5) hold, then for $\tilde{r}$ as in (3.6), $\exists \beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \le \frac{M}{\eta \sqrt{p}} + \sqrt{\frac{M \|r - \alpha e\|_2}{\eta p^{1/2}}} $$

holds w.h.p. If $r_i = i$, we show in Example 4 that if $p \gtrsim \frac{\log n}{\eta^2 \delta^{2/2}}$ for a fixed $\delta \in (0, 1)$, then $\|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n^{3/2}}{(\log n)^{1/4}} \sqrt{\delta}$ w.h.p. Hence $\Omega(n \log n)$ measurements suffice for $\ell_2$ recovery.

(b) If the conditions in (3.5) hold and $n$ is large enough, then for $\hat{r}$ as in (3.6), there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \|r - \alpha e\|_2 [C(n, M, \eta, p, r) + \sqrt{n} C^2(n, M, \eta, p, r)] + \frac{M(M - \alpha)}{\eta \sqrt{p} \|r - \alpha e\|_2}, $$

holds with $C(n, M, \eta, p, r)$ as defined in (4.3). When $r_i = i$, we show in Example 4 that if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{(\log n)^{2/3}}{n^{2/3}})$ for a fixed $\delta \in (0, 1)$, and $n$ is large enough, then $\|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \frac{n}{\log n} \delta$ holds w.h.p.

3. In the same spirit as Theorem 8, Theorem 10 provides a $\ell_2$ error bound (up to a global sign and shift) for SVD-NRS for the score vector $r$ and a vector $\tilde{r}$ (derived from $\tilde{u}_2$ in Step 6) where

$$ \tilde{w} = \frac{\hat{\sigma}_1}{\eta p \|\bar{D}^{-1/2} e\|_2} \tilde{u}_2; \quad \tilde{r} := \tilde{w} - \frac{e^T \tilde{w}}{n} e. \qquad (3.7) $$
the ranking $\tilde{\pi}$ obtained in Step 7 of Algorithm 1, and the ground truth ranking $\pi$. Denoting $\alpha = \frac{r^T e}{n}$ to be the average score, it states that if

$$ \|r - \alpha e\|_2 \gtrsim \frac{M}{\eta\sqrt{p}}, \quad p \gtrsim \frac{\log n}{n}, \qquad (3.5) $$

and if $n$ is large enough, then $\|\tilde{\pi} - \pi\|_{\infty} \lesssim \frac{\|r - \alpha e\|_2}{\rho} \Upsilon(n, M, \eta, p, \varepsilon, r)$ w.h.p. Here $\Upsilon(n, M, \eta, p, \varepsilon, r)$ is the bound on the $\ell_{\infty}$ error between $\tilde{u}_2$ and $u_2$ (up to a global sign) in Theorem 6. For concreteness, if $r_i = i$, we show in Example 3 that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{\log n)^{2/3}}{n^{2/3}}$ and $n$ is large enough, then $\|\tilde{\pi} - \pi\|_{\infty} \lesssim \delta \frac{n}{\log n}$ w.h.p. Note that this implies that $\Omega(n^{4/3}(\log n)^{2/3})$ measurements suffice.

2. Theorem 8 provides $\ell_2$ and $\ell_{\infty}$ bounds (up to a global sign and shift) for SVD-RS for the score vector $r$ and a vector $\hat{r}$ (derived from $\tilde{u}_2$) where

$$ \tilde{w} = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}} \tilde{u}_2; \quad \tilde{r} := \tilde{w} - \frac{e^T \tilde{w}}{n} e. \qquad (3.6) $$

If one knew the product $\eta p$, then (3.6) corresponds to choosing the scale $\tau = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}}$ in Step 8.

(a) If the conditions in (3.5) hold, then for $\tilde{r}$ as in (3.6), $\exists \beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \le \frac{M}{\eta\sqrt{p}} + \sqrt{\frac{M\|r - \alpha e\|_2}{\eta p^{1/2}}} $$

holds w.h.p. If $r_i = i$, we show in Example 4 that if $p \gtrsim \frac{\log n}{\eta^2 \delta^{2/3}}$ for a fixed $\delta \in (0, 1)$, then $\|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n^{3/2}}{(\log n)^{1/4}} \sqrt{\delta}$ w.h.p. Hence $\Omega(n \log n)$ measurements suffice for $\ell_2$ recovery.

(b) If the conditions in (3.5) hold and $n$ is large enough, then for $\hat{r}$ as in (3.6), there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \|r - \alpha e\|_2 [C(n, M, \eta, p, r) + \sqrt{n} C^2(n, M, \eta, p, r)] + \frac{M(M - \alpha)}{\eta \sqrt{p} \|r - \alpha e\|_2}, $$

holds with $C(n, M, \eta, p, r)$ as defined in (4.3). When $r_i = i$, we show in Example 4 that if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{\log n)^{2/3}}{n^{2/3}}$ for a fixed $\delta \in (0, 1)$, and $n$ is large enough, then $\|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \frac{n}{\log n} \delta$ holds w.h.p.

3. In the same spirit as Theorem 8, Theorem 10 provides a $\ell_2$ error bound (up to a global sign and shift) for SVD-NRS for the score vector $r$ and a vector $\tilde{r}$ (derived from $\tilde{u}_2$ in Step 6) where

$$ \tilde{w} = \frac{\hat{\sigma}_1}{\eta p \|\bar{D}^{-1/2} e\|_2} \tilde{u}_2; \quad \tilde{r} := \tilde{w} - \frac{e^T \tilde{w}}{n} e. \qquad (3.7) $$
Again, if one knew the product $ηp$, then (3.7) corresponds to choosing the scale $τ = \frac{\sigma_1}{ηp||D^{-1/2}e||_2}$ in Step 8. While the conditions on $p$ are more convoluted than those in Theorem 8 – in large part due to the normalization – the essence is the same. In the setting where $r_i = i$, we can derive a simplified version of the theorem as shown in Example 6. Let us now denote $α = \frac{r^T(E[D])^{-1}e}{e^T(E[D])^{-1}e}$ and $α' = \frac{e^T(r-αe)}{n}$. Theorem 10 states that if $p \gtrsim \frac{\log n}{\eta^2\delta^2n}$ for fixed $δ \in (0, 1)$, then for $\tilde{r}$ as in (3.7), there exists $β \in \{-1, 1\}$ such that w.h.p,

$$ \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \lesssim n^{3/2}\eta^{1/4}\delta^{1/2} + \delta^{1/2}\frac{n^{3/2}}{(\log n)^{1/4}}. $$

The closest works to our setting are those of Gleich and Lim (2011) and Yang and Wakin (2015) which rely on matrix completion, followed by a row sum to compute rankings. In our experiments, we also consider matrix completion as a preprocessing step (see Section 6) before applying our SVD-based algorithms, and find that in general it does improve the performance of our methods. However, this comes at the cost of an additional computational overhead; our SVD-based methods rely on computing the top two singular vectors of the input matrix, which in practice can be done very efficiently. The theoretical guarantees we provide for our problem setting, are new to the best of our knowledge, and have been lacking in the ranking and synchronization literatures currently.

# 4. Theoretical results for SVD-RS and SVD-NRS

We begin by detailing the theoretical results for SVD-RS in Section 4.1, followed by those for SVD-NRS in Section 4.2. Throughout, we instantiate the main theorems for the special setting where $r_i = i$, for ease of interpretability.

## 4.1 Analysis of Algorithm 1 (SVD-RS)

We will assume throughout this section that for $\{i, j\} \in E$, the measurement $R_{ij}$ corresponds to $i < j$. This is clearly without loss of generality. From (3.1), we can write the entry $H_{ij}$ (for $i < j$) as the following mixture model

$$ H_{ij} = \begin{cases} r_i - r_j & \text{w.p. } ηp \\ N_{ij} \sim U[-M, M] & \text{w.p. } (1-η)p \\ 0 & \text{w.p. } 1-p, \end{cases} $$

with $(H_{ij})_{i<j}$ being independent random variables. Since $H_{ii} = 0$ and $H_{ij} = -H_{ji}$ (by construction of $H$), therefore for all $i, j \in [n]$, we obtain $E[H_{ij}] = ηp(r_i - r_j)$. In particular, we have that

$$ E[H] = ηpC, \qquad (4.1) $$

where $C = re^T - er^T$ is a skew-symmetric matrix of rank 2. We can decompose $H$ as

$$ H = E[H] + Z = ηpC + Z, $$
Again, if one knew the product $ηp$, then (3.7) corresponds to choosing the scale $τ = \frac{\sigma_1}{ηp||D^{-1/2}e||_2}$ in Step 8. While the conditions on $p$ are more convoluted than those in Theorem 8 – in large part due to the normalization – the essence is the same. In the setting where $r_i = i$, we can derive a simplified version of the theorem as shown in Example 6. Let us now denote $α = \frac{r^T(E[D])^{-1}e}{e^T(E[D])^{-1}e}$ and $α' = \frac{e^T(r-αe)}{n}$. Theorem 10 states that if $p \gtrsim \frac{\log n}{\eta^2\delta^2n}$ for fixed $δ \in (0,1)$, then for $\tilde{r}$ as in (3.7), there exists $β \in \{-1, 1\}$ such that w.h.p,

$$ \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \lesssim n^{3/2}\eta^{1/4}\delta^{1/2} + \delta^{1/2}\frac{n^{3/2}}{(\log n)^{1/4}}. $$

The closest works to our setting are those of Gleich and Lim (2011) and Yang and Wakin (2015) which rely on matrix completion, followed by a row sum to compute rankings. In our experiments, we also consider matrix completion as a preprocessing step (see Section 6) before applying our SVD-based algorithms, and find that in general it does improve the performance of our methods. However, this comes at the cost of an additional computational overhead; our SVD-based methods rely on computing the top two singular vectors of the input matrix, which in practice can be done very efficiently. The theoretical guarantees we provide for our problem setting, are new to the best of our knowledge, and have been lacking in the ranking and synchronization literatures currently.

# 4. Theoretical results for SVD-RS and SVD-NRS

We begin by detailing the theoretical results for SVD-RS in Section 4.1, followed by those for SVD-NRS in Section 4.2. Throughout, we instantiate the main theorems for the special setting where $r_i = i$, for ease of interpretability.

## 4.1 Analysis of Algorithm 1 (SVD-RS)

We will assume throughout this section that for $\{i, j\} \in E$, the measurement $R_{ij}$ corresponds to $i < j$. This is clearly without loss of generality. From (3.1), we can write the entry $H_{ij}$ (for $i < j$) as the following mixture model

$$ H_{ij} = \begin{cases} r_i - r_j & \text{w.p. } ηp \\ N_{ij} \sim U[-M, M] & \text{w.p. } (1-η)p \\ 0 & \text{w.p. } 1-p, \end{cases} $$

with $(H_{ij})_{i<j}$ being independent random variables. Since $H_{ii} = 0$ and $H_{ij} = -H_{ji}$ (by construction of $H$), therefore for all $i, j \in [n]$, we obtain $E[H_{ij}] = ηp(r_i - r_j)$. In particular, we have that

$$ E[H] = ηpC, \qquad (4.1) $$

where $C = re^T - er^T$ is a skew-symmetric matrix of rank 2. We can decompose $H$ as

$$ H = E[H] + Z = ηpC + Z, $$
where Z is a random noise matrix with $Z_{ii} = 0, \forall i = 1, \dots, n$. For $1 \le i < j \le n$ the entries of Z are defined by the following mixture model

$$ Z_{ij} = \begin{cases} (r_i - r_j) - \eta p(r_i - r_j) & ; \text{w.p } \eta p \\ N_{ij} - \eta p(r_i - r_j) & ; \text{w.p } (1 - \eta)p \\ -\eta p(r_i - r_j) & ; \text{w.p } (1 - p). \end{cases} \quad (4.2) $$

Note that $(Z_{ij})_{i<j}$ are independent<fnref n="3" /> random variables and $Z_{ij} = -Z_{ji}$. Hence the matrices H, C, Z are all skew-symmetric. We now proceed to show in Theorem 5 that, provided p is large enough, it holds true that $\tilde{u}_2$ (obtained in Step 6 of SVD-RS) is close to $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ in the $l_2$ norm (up to a global sign) with high probability. The proof is outlined in Section 5.1.

**Theorem 5 (l₂ recovery, ERO model)** Denoting $\alpha = \frac{r^T e}{n}$, for a given $0 < \varepsilon \le 1/2$, let $\|r - \alpha e\|_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}}(2 + \varepsilon)$. Let $\tilde{u}_2 \in \mathbb{R}^n$ be the vector obtained in Step 6 of SVD-RS (with $\|\tilde{u}_2\|_2 = 1$) and $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$. Then, there exists $\beta \in \{-1, 1\}$ and a constant $c_\varepsilon > 0$ depending only on $\varepsilon$ such that

$$ \|\tilde{u}_2 - \beta u_2\|_2^2 \le \frac{120M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2 + \varepsilon)}{\|r - \alpha e\|_2}, $$

with probability at least $1 - 2n \exp(-\frac{20pn}{3c_\varepsilon})$.

Theorem 5 says that for $\delta \in (0, 1)$, if $p \ge \max\{\frac{M^2 \log n}{\eta^2 \delta^2 \|r-\alpha e\|_2^2}, \frac{\log n}{n}\}$, then with high probability, we have

$$ \|\tilde{u}_2 - \beta u_2\|_2^2 \lesssim \frac{\delta}{\sqrt{\log n}}. $$

Note that the bounds involved are invariant to the scaling of r, and essentially depend on the variance of the normalized entries of r, wherein each $\frac{r_i}{M} \in [0, 1]$. We can see that as $\frac{\|r-\alpha e\|_2^2}{M^2}$ becomes small, then the corresponding condition on p becomes more stringent.

**Example 1** Consider the case where $r_i = i$ for $i = 1, \dots, n$. Then $M = n$, $\alpha = \frac{n+1}{2}$ and $\|r - \alpha e\|_2 = \Theta(n^{3/2})$. Hence Theorem 5 now says that if $p \ge \frac{\log n}{\eta^2 \delta^2 n}$, then with high probability, we have $\|\tilde{u}_2 - \beta u_2\|_2^2 \lesssim \frac{\delta}{\sqrt{\log n}}$.

If in addition, n is also large enough, then we show in Theorem 6 that $\tilde{u}_2$ is close to $u_2$ (up to a global sign) in the $l_\infty$ norm as well. The proof is outlined in Section 5.2.

**Theorem 6 (l_\infty recovery, ERO model)** With the same notation as in Theorem 5, for a given $0 < \varepsilon \le 1/2$, let $\|r - \alpha e\|_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}}(2 + \varepsilon)$. Assume $p \ge \max\{\frac{1}{2n}, \frac{2 \log n}{15n}\}$. Choose $\xi > 1, 0 < \kappa < 1$ and define $\mu = \frac{2}{\kappa+1}$. Let n satisfy $\frac{16}{\kappa} \le (\log n)^\xi$. Then, there exists $\beta \in \{-1, 1\}$ and constants $C_\varepsilon, c_\varepsilon > 0$ depending only on $\varepsilon$ such that with probability at least

$$ 1 - 2n \exp\left(-\frac{20pn}{3c_{\varepsilon}}\right) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}} $$

<footnote n="3">³Since $Z_{ii} = 0$, clearly $(Z_{ij})_{i \le j}$ are independent as well.</footnote>
where Z is a random noise matrix with $Z_{ii} = 0, \forall i = 1, \dots, n$. For $1 \le i < j \le n$ the entries of Z are defined by the following mixture model

$$ Z_{ij} = \begin{cases} (r_i - r_j) - \eta p(r_i - r_j) & ; \text{w.p } \eta p \\ N_{ij} - \eta p(r_i - r_j) & ; \text{w.p } (1 - \eta)p \\ -\eta p(r_i - r_j) & ; \text{w.p } (1 - p). \end{cases} \quad (4.2) $$

Note that $(Z_{ij})_{i<j}$ are independent<fnref n="3" /> random variables and $Z_{ij} = -Z_{ji}$. Hence the matrices H, C, Z are all skew-symmetric. We now proceed to show in Theorem 5 that, provided p is large enough, it holds true that $\tilde{u}_2$ (obtained in Step 6 of SVD-RS) is close to $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$ in the $l_2$ norm (up to a global sign) with high probability. The proof is outlined in Section 5.1.

**Theorem 5 (l₂ recovery, ERO model)** Denoting $\alpha = \frac{r^T e}{n}$, for a given $0 < \varepsilon \le 1/2$, let $\|r - \alpha e\|_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}}(2 + \varepsilon)$. Let $\tilde{u}_2 \in \mathbb{R}^n$ be the vector obtained in Step 6 of SVD-RS (with $\|\tilde{u}_2\|_2 = 1$) and $u_2 = \frac{r-\alpha e}{\|r-\alpha e\|_2}$. Then, there exists $\beta \in \{-1, 1\}$ and a constant $c_\varepsilon > 0$ depending only on $\varepsilon$ such that

$$ \|\tilde{u}_2 - \beta u_2\|_2^2 \le \frac{120M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2 + \varepsilon)}{\|r - \alpha e\|_2}, $$

with probability at least $1 - 2n \exp(-\frac{20pn}{3c_\varepsilon})$.

Theorem 5 says that for $\delta \in (0, 1)$, if $p \ge \max\{\frac{M^2 \log n}{\eta^2 \delta^2 \|r-\alpha e\|_2^2}, \frac{\log n}{n}\}$, then with high probability, we have

$$ \|\tilde{u}_2 - \beta u_2\|_2^2 \lesssim \frac{\delta}{\sqrt{\log n}}. $$

Note that the bounds involved are invariant to the scaling of r, and essentially depend on the variance of the normalized entries of r, wherein each $\frac{r_i}{M} \in [0, 1]$. We can see that as $\frac{\|r-\alpha e\|_2^2}{M^2}$ becomes small, then the corresponding condition on p becomes more stringent.

**Example 1** Consider the case where $r_i = i$ for $i = 1, \dots, n$. Then $M = n$, $\alpha = \frac{n+1}{2}$ and $\|r - \alpha e\|_2 = \Theta(n^{3/2})$. Hence Theorem 5 now says that if $p \ge \max\{\frac{1}{2n}, \frac{2 \log n}{15n}\}$, then with high probability, we have $\|\tilde{u}_2 - \beta u_2\|_2^2 \lesssim \frac{\delta}{\sqrt{\log n}}$.

If in addition, n is also large enough, then we show in Theorem 6 that $\tilde{u}_2$ is close to $u_2$ (up to a global sign) in the $l_\infty$ norm as well. The proof is outlined in Section 5.2.

**Theorem 6 (l_\infty recovery, ERO model)** With the same notation as in Theorem 5, for a given $0 < \varepsilon \le 1/2$, let $\|r - \alpha e\|_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}}(2 + \varepsilon)$. Assume $p \ge \max\{\frac{1}{2n}, \frac{2 \log n}{15n}\}$. Choose $\xi > 1, 0 < \kappa < 1$ and define $\mu = \frac{2}{\kappa+1}$. Let n satisfy $\frac{16}{\kappa} \le (\log n)^\xi$. Then, there exists $\beta \in \{-1, 1\}$ and constants $C_\varepsilon, c_\varepsilon > 0$ depending only on $\varepsilon$ such that with probability at least

$$ 1 - 2n \exp\left(-\frac{20pn}{3c_{\varepsilon}}\right) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}} $$

<footnote n="3">³Since $Z_{ii} = 0$, clearly $(Z_{ij})_{i \le j}$ are independent as well.</footnote>
we have that

$$ ||\tilde{u}_2 - \beta u_2||_\infty \le 4(2 + \sqrt{2})C(n, M, \eta, p, \varepsilon, r) + 4\sqrt{n}C^2(n, M, \eta, p, \varepsilon, r), $$

where

$$ C(n, M, \eta, p, \varepsilon, r) = C_{\varepsilon} \left[ \left( \frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r - \alpha e\|_2} + \frac{M^2(\log n)^{2\xi}}{\eta^2 p\|r - \alpha e\|_2^2} \right) \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right) + \frac{M^3}{\eta^3 p^{3/2}\|r - \alpha e\|_2^3} \right]. \quad (4.3) $$

Let us look at bounding the admittedly complicated looking term $C(n, M, \eta, p, \varepsilon, r)$. For $\delta \in (0, 1)$, say $p$ satisfies

$$ p \gtrsim \frac{M^2}{\|r - \alpha e\|_2^2} \max \left\{ \frac{n \log^3 n}{\eta^2 \delta^2} \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right)^2, \frac{\sqrt{n}(\log n)^{2\xi+1}}{\eta^2 \delta} \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right), \frac{(\sqrt{n} \log n)^{2/3}}{\eta^2 \delta^{2/3}} \right\}. \quad (4.4) $$

Then we have that $C(\cdot) \lesssim \frac{\delta}{\sqrt{n} \log n}$. Hence if additionally $p \gtrsim \max\left\{\frac{M^2}{\eta^2 \|r - \alpha e\|_2^2}, \frac{\log n}{n}\right\}$ and $n = \Omega(1)$ hold, then Theorem 6 says that

$$ ||\tilde{u}_2 - \beta u_2||_\infty \lesssim \frac{\delta}{\sqrt{n} \log n} $$

holds w.h.p. Note that the condition on $p$ is stricter than that in Theorem 5.

**Example 2** Let us revisit the case where $r_i = i$ for all $i$. The condition (4.4) then simplifies to

$$ p \gtrsim \max \left\{ \frac{1}{\eta^2 \delta^2} \left( \frac{(\log n)^3}{n} \right), \frac{1}{\eta^2 \delta} \left( \frac{(\log n)^{2\xi+1}}{n} \right), \frac{1}{\eta^2 \delta^{2/3}} \left( \frac{(\log n)^{2/3}}{n^{2/3}} \right) \right\}. \quad (4.5) $$

Since $\frac{M^2}{\|r - \alpha e\|_2^2} = \frac{1}{n}$, we conclude that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} \left( \frac{(\log n)^{2/3}}{n^{2/3}} \right)$ and $n = \Omega(1)$, then $||\tilde{u}_2 - \beta u_2||_\infty \lesssim \frac{\delta}{\sqrt{n} \log n}$ holds with high probability.

**Recovering ranks in the $\ell_\infty$ norm.** Let us assume that $\tilde{u}_2$ is aligned with $u_2$, i.e., $||\tilde{u}_2 - u_2||_\infty \le \Upsilon(n, M, \eta, p, \varepsilon, r)$ where $\Upsilon(\cdot)$ is the bound in Theorem 6. For any $s \in \mathbb{R}^n$, we say that the permutation $\pi: [n] \to [n]$ is consistent with $s$ if for all pairs $(i, j)$, $\pi(i) < \pi(j)$ implies $s_i \ge s_j$. For the purpose of recovering the rankings, we will consider the entries of $r$ to be pairwise distinct, i.e. $r_i \ne r_j$ for all $i \ne j$. Hence there is a unique permutation $\pi$ which is consistent with $r$. Let $\tilde{\pi}$ be a permutation that is consistent with $\tilde{u}_2$, note that this is not unique since some entries of $\tilde{u}_2$ could have the same value. Our goal is to bound the maximum displacement error between $\pi$ and $\tilde{\pi}$ defined as

$$ ||\tilde{\pi} - \pi||_\infty := \max_i \left( \sum_{\pi(j) > \pi(i)} \mathbf{1}_{\tilde{\pi}(j) < \tilde{\pi}(i)} + \sum_{\pi(j) < \pi(i)} \mathbf{1}_{\tilde{\pi}(j) > \tilde{\pi}(i)} \right). \quad (4.6) $$

To this end, we have the following theorem the proof of which is provided in Section 5.3. Notice that bounding $||\pi - \tilde{\pi}||_\infty$ requires a bound on $||\tilde{u}_2 - u_2||_\infty$ which is obtained from Theorem 6. The proof technique is essentially the same as that of (Fogel et al., 2016, Theorem 24).
we have that

$$ ||\tilde{u}_2 - \beta u_2||_{\infty} \le 4(2 + \sqrt{2})C(n, M, \eta, p, \varepsilon, r) + 4\sqrt{n}C^2(n, M, \eta, p, \varepsilon, r), $$

where

$$ C(n, M, \eta, p, \varepsilon, r) = C_{\varepsilon} \left[ \left( \frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r - \alpha e\|_2} + \frac{M^2(\log n)^{2\xi}}{\eta^2 p\|r - \alpha e\|_2^2} \right) \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right) + \frac{M^3}{\eta^3 p^{3/2}\|r - \alpha e\|_2^3} \right]. \quad (4.3) $$

Let us look at bounding the admittedly complicated looking term $C(n, M, \eta, p, \varepsilon, r)$. For $\delta \in (0, 1)$, say $p$ satisfies

$$ p \gtrsim \frac{M^2}{\|r - \alpha e\|_2^2} \max \left\{ \frac{n \log^3 n}{\eta^2 \delta^2} \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right)^2, \frac{\sqrt{n}(\log n)^{2\xi+1}}{\eta^2 \delta} \left( \frac{1}{\sqrt{n}} + \frac{M - \alpha}{\|r - \alpha e\|_2} \right), \frac{(\sqrt{n} \log n)^{2/3}}{\eta^2 \delta^{2/3}} \right\}. \quad (4.4) $$

Then we have that $C(\cdot) \lesssim \frac{\delta}{\sqrt{n} \log n}$. Hence if additionally $p \gtrsim \max\left\{\frac{M^2}{\eta^2 \|r - \alpha e\|_2^2}, \frac{\log n}{n}\right\}$ and $n = \Omega(1)$ hold, then Theorem 6 says that

$$ ||\tilde{u}_2 - \beta u_2||_{\infty} \lesssim \frac{\delta}{\sqrt{n} \log n} $$

holds w.h.p. Note that the condition on $p$ is stricter than that in Theorem 5.

**Example 2** Let us revisit the case where $r_i = i$ for all $i$. The condition (4.4) then simplifies to

$$ p \gtrsim \max \left\{ \frac{1}{\eta^2 \delta^2} \left( \frac{(\log n)^3}{n} \right), \frac{1}{\eta^2 \delta} \left( \frac{(\log n)^{2\xi+1}}{n} \right), \frac{1}{\eta^2 \delta^{2/3}} \left( \frac{(\log n)^{2/3}}{n^{2/3}} \right) \right\}. \quad (4.5) $$

Since $\frac{M^2}{\|r - \alpha e\|_2^2} = \frac{1}{n}$, we conclude that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} \left( \frac{(\log n)^{2/3}}{n^{2/3}} \right)$ and $n = \Omega(1)$, then $||\tilde{u}_2 - \beta u_2||_{\infty} \lesssim \frac{\delta}{\sqrt{n} \log n}$ holds with high probability.

**Recovering ranks in the $\ell_{\infty}$ norm.** Let us assume that $\tilde{u}_2$ is aligned with $u_2$, i.e., $||\tilde{u}_2 - u_2||_{\infty} \le \Upsilon(n, M, \eta, p, \varepsilon, r)$ where $\Upsilon(\cdot)$ is the bound in Theorem 6. For any $s \in \mathbb{R}^n$, we say that the permutation $\pi: [n] \to [n]$ is consistent with $s$ if for all pairs $(i, j)$, $\pi(i) < \pi(j)$ implies $s_i \ge s_j$. For the purpose of recovering the rankings, we will consider the entries of $r$ to be pairwise distinct, i.e. $r_i \ne r_j$ for all $i \ne j$. Hence there is a unique permutation $\pi$ which is consistent with $r$. Let $\tilde{\pi}$ be a permutation that is consistent with $\tilde{u}_2$, note that this is not unique since some entries of $\tilde{u}_2$ could have the same value. Our goal is to bound the maximum displacement error between $\pi$ and $\tilde{\pi}$ defined as

$$ ||\tilde{\pi} - \pi||_{\infty} := \max_i \left( \sum_{\pi(j) > \pi(i)} \mathbf{1}_{\tilde{\pi}(j) < \tilde{\pi}(i)} + \sum_{\pi(j) < \pi(i)} \mathbf{1}_{\tilde{\pi}(j) > \tilde{\pi}(i)} \right). \quad (4.6) $$

To this end, we have the following theorem the proof of which is provided in Section 5.3. Notice that bounding $||\pi - \tilde{\pi}||_{\infty}$ requires a bound on $||\tilde{u}_2 - u_2||_{\infty}$ which is obtained from Theorem 6. The proof technique is essentially the same as that of (Fogel et al., 2016, Theorem 24).
**Theorem 7** Assuming $r_i \neq r_j$ for all $i \neq j$, let $\pi$ denote the (unique) ranking consistent with $r$ and also define $\rho := \min_{i \neq j} |r_i - r_j|$. Assuming $\beta = 1$, denote $\Upsilon(n, M, \eta, p, \varepsilon, r)$ to be the bound on $\|\tilde{u}_2 - u_2\|_{\infty}$ in Theorem 6. Then under the notation and assumptions in Theorem 6, we have with probability at least $1 - 2n \exp(-\frac{20pn}{3c_{\varepsilon}}) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}}$ that

$$ \|\tilde{\pi} - \pi\|_{\infty} \le \frac{4\|r - \alpha e\|_2}{\rho} \Upsilon(n, M, \eta, p, \varepsilon, r) \quad (4.7) $$

holds for all rankings $\tilde{\pi}$ which are consistent with $\tilde{u}_2$.

**Example 3** Let us examine the bound in (4.7) for the case where $r_i = i$ for $i = 1, \dots, n$. We have seen that $\|r - \alpha e\|_2 = \Theta(n^{3/2})$. Theorem 6 tells us that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{(\log n)^{2/3}}{n^{2/3}})$ and $n = \Omega(1)$, then $\|\tilde{u}_2 - u_2\|_{\infty} \le \Upsilon(\cdot) \lesssim \frac{\delta}{\sqrt{n \log n}} \text{ w.h.p.}$. Since $\rho = 1$, we obtain from (4.7) that

$$ \|\tilde{\pi} - \pi\|_{\infty} \lesssim \delta \frac{n}{\log n}. $$

**Recovering scores in the $l_2$ and $l_{\infty}$ norms.** We now bound the error between the score vector estimate $\tilde{r}$ as in (3.6), and $r$ up to a global shift and sign. This is shown in the following theorem for the $l_2$ and $l_{\infty}$ norms. The proof is outlined in Section 5.4. Note that both $\tilde{r}$ and $r - \alpha e$ are centered vectors.

**Theorem 8** Recall $\tilde{r} \in \mathbb{R}^n$ as defined in (3.6).

1. Under the notation and assumptions of Theorem 5, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \le \frac{8M}{\eta} \sqrt{\frac{5}{3p}}(2 + \varepsilon) + \sqrt{\frac{120M(2 + \varepsilon)\|r - \alpha e\|_2}{\eta p^{1/2}}} \left(\frac{5}{3}\right)^{1/4} \quad (4.8) $$

with probability at least $1 - 2n \exp(-\frac{20pn}{3c_{\varepsilon}})$.

2. Under the notation and assumptions of Theorem 6, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \le \frac{16}{3}\|r - \alpha e\|_2[(2 + \sqrt{2})C(n, M, \eta, p, \varepsilon, r) + \sqrt{n}C^2(n, M, \eta, p, \varepsilon, r)] \\ + 8\sqrt{\frac{5}{3}}(2 + \varepsilon)\frac{M(M - \alpha)}{\eta\sqrt{p}\|r - \alpha e\|_2} \quad (4.9) $$

with probability at least

$$ 1 - 2n \exp \left( -\frac{20pn}{3c_{\varepsilon}} \right) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}}. $$

**Example 4** When $r_i = i$ for all $i$, the bound in (4.8) is of the form

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n}{\eta\sqrt{p}} + \frac{n^{5/4}}{\sqrt{\eta p}^{1/4}}. \quad (4.10) $$
**Theorem 7** Assuming $r_i \neq r_j$ for all $i \neq j$, let $\pi$ denote the (unique) ranking consistent with $r$ and also define $\rho := \min_{i \neq j} |r_i - r_j|$. Assuming $\beta = 1$, denote $\Upsilon(n, M, \eta, p, \varepsilon, r)$ to be the bound on $\|\tilde{u}_2 - u_2\|_{\infty}$ in Theorem 6. Then under the notation and assumptions in Theorem 6, we have with probability at least $1-2n \exp(-\frac{20pn}{3c_{\varepsilon}}) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}}$ that

$$ \|\tilde{\pi} - \pi\|_{\infty} \le \frac{4\|r - \alpha e\|_2}{\rho} \Upsilon(n, M, \eta, p, \varepsilon, r) \quad (4.7) $$

holds for all rankings $\tilde{\pi}$ which are consistent with $\tilde{u}_2$.

**Example 3** Let us examine the bound in (4.7) for the case where $r_i = i$ for $i = 1, \dots, n$. We have seen that $\|r - \alpha e\|_2 = \Theta(n^{3/2})$. Theorem 6 tells us that for a fixed $\delta \in (0, 1)$, if $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}} (\frac{(\log n)^{2/3}}{n^{2/3}})$ and $n = \Omega(1)$, then $\|\tilde{u}_2 - u_2\|_{\infty} \le \Upsilon(\cdot) \lesssim \frac{\delta}{\sqrt{n \log n}} \text{ w.h.p.}$. Since $\rho = 1$, we obtain from (4.7) that

$$ \|\tilde{\pi} - \pi\|_{\infty} \lesssim \delta \frac{n}{\log n}. $$

**Recovering scores in the $l_2$ and $l_{\infty}$ norms.** We now bound the error between the score vector estimate $\tilde{r}$ as in (3.6), and $r$ up to a global shift and sign. This is shown in the following theorem for the $l_2$ and $l_{\infty}$ norms. The proof is outlined in Section 5.4. Note that both $\tilde{r}$ and $r - \alpha e$ are centered vectors.

**Theorem 8** Recall $\tilde{r} \in \mathbb{R}^n$ as defined in (3.6).

1. Under the notation and assumptions of Theorem 5, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \le \frac{8M}{\eta} \sqrt{\frac{5}{3p}} (2 + \varepsilon) + \sqrt{\frac{120M(2 + \varepsilon)\|r - \alpha e\|_2}{\eta p^{1/2}}} \left(\frac{5}{3}\right)^{1/4} \quad (4.8) $$

with probability at least $1 - 2n \exp(-\frac{20pn}{3c_{\varepsilon}})$.

2. Under the notation and assumptions of Theorem 6, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \le \frac{16}{3}\|r - \alpha e\|_2[(2 + \sqrt{2})C(n, M, \eta, p, \varepsilon, r) + \sqrt{n}C^2(n, M, \eta, p, \varepsilon, r)] \\ + 8\sqrt{\frac{5}{3}}(2 + \varepsilon)\frac{M(M - \alpha)}{\eta\sqrt{p}\|r - \alpha e\|_2} \quad (4.9) $$

with probability at least

$$ 1 - 2n \exp \left( -\frac{20pn}{3c_{\varepsilon}} \right) - \frac{4}{n} - 2n^{1-\frac{1}{4}(\log_{\mu} n)^{\xi-1}(\log_{\mu} e)^{-\xi}}. $$

**Example 4** When $r_i = i$ for all $i$, the bound in (4.8) is of the form

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n}{\eta\sqrt{p}} + \frac{n^{5/4}}{\sqrt{\eta p}^{1/4}}. \quad (4.10) $$
Recall from Theorem 5 that $p$ is required to satisfy $p \gtrsim \frac{\log n}{\eta^2 \delta^{2/2} n}$ for $\delta \in (0, 1)$. Then, (4.10) simplifies to

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n^{3/2}}{(\log n)^{1/4}} \sqrt{\delta}. $$

Turning our attention to (4.9), recall from Theorem 6 that we require $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}}$ $(\frac{(\log n)^{2/3}}{n^{2/3}})$ and $n = \Omega(1)$. This leads to

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \frac{n}{\log n} \delta + \frac{\sqrt{n}}{\eta\sqrt{p}} \lesssim \frac{n}{\log n} \delta + \frac{n^{5/6}}{(\log n)^{1/3}} \delta^{1/3} \lesssim \frac{n}{\log n} \delta. $$

## 4.2 Analysis of Algorithm 2 (SVD-NRS)

We now analyze Algorithm 2, namely SVD-NRS, under the ERO model. The random matrix $\bar{D}^{-1/2} H \bar{D}^{-1/2}$ will concentrate around $(\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2}$. Recall from (4.1) that $\mathbb{E}[H] = \eta p(re^T - er^T)$, hence $(\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2}$ is a rank-2 matrix. It is easy to verify that

$$ \mathbb{E}[\bar{D}_{ii}] = \eta p \sum_{j=1}^{n} |r_i - r_j| + (1 - \eta) p \frac{M}{2}. $$

It will be useful to denote

$$ \max_i \mathbb{E}[\bar{D}_{ii}] = p \left[ \eta \max_i \left( \sum_{j=1}^n |r_i - r_j| \right) + (1 - \eta) \frac{M}{2} \right] = p \lambda_{\max}, \quad (4.11) $$

$$ \min_i \mathbb{E}[\bar{D}_{ii}] = p \left[ \eta \min_i \left( \sum_{j=1}^n |r_i - r_j| \right) + (1 - \eta) \frac{M}{2} \right] = p \lambda_{\min}. \quad (4.12) $$

Recall $\tilde{u}_2 \in \mathbb{R}^n$ obtained in Step 6 of SVD-NRS, and denote $u_2 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e)}{\| (\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e) \|_2}$. The following theorem bounds the $\ell_2$ error between $\tilde{u}_2$ and $u_2$ up to a global sign, its proof is given in Section 5.5.

**Theorem 9 ($\ell_2$ recovery, ERO model) Denote**

$$ \alpha = \frac{r^T (\mathbb{E}[\bar{D}])^{-1} e}{e^T (\mathbb{E}[\bar{D}])^{-1} e}, \quad A(\eta, M) = \eta M^2 + (1 - \eta) \frac{M}{2} \quad \text{and} \quad C_1(\eta, M) = 4A^{1/4}(\eta, M). $$

For $\varepsilon > 0$, denote

$$ \sigma_{\min} = \frac{\eta \|r - \alpha e\|_2 \sqrt{n}}{\lambda_{\max}}, \quad \sigma_{\max} = \frac{\eta \|r - \alpha e\|_2 \sqrt{n}}{\lambda_{\min}}, $$

$$ \tilde{\Delta} = 16M \sqrt{\frac{5}{3}pn} \frac{(2+\varepsilon)}{p\lambda_{\min}} + \frac{C_1(\eta, M)(np \log n)^{1/4}\sigma_{\max}}{p^{3/2}\lambda_{\min}^{3/2}} \left( \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right). $$

Let $\tilde{u}_2 \in \mathbb{R}^n$ be the vector obtained in Step 6 of SVD-NRS (with $\|\tilde{u}_2\|_2 = 1$) and denote $u_2 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e)}{\| (\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e) \|_2}$. If

$$ p \ge \frac{M^2}{9A(\eta, M)} \frac{\log n}{n}, \quad p \ge \frac{16(\sqrt{2} + 1)^2 A(\eta, M) n \log n}{\lambda_{\min}^2} \quad \text{and} \quad \tilde{\Delta} \le \frac{\sigma_{\min}}{3} \quad (4.13) $$
Recall from Theorem 5 that $p$ is required to satisfy $p \gtrsim \frac{\log n}{\eta^2 \delta^{2/2} n}$ for $\delta \in (0, 1)$. Then, (4.10) simplifies to

$$ \|\tilde{r} - \beta(r - \alpha e)\|_2 \lesssim \frac{n^{3/2}}{(\log n)^{1/4}} \sqrt{\delta}. $$

Turning our attention to (4.9), recall from Theorem 6 that we require $p \gtrsim \frac{1}{\eta^2 \delta^{2/3}}$ $(\frac{(\log n)^{2/3}}{n^{2/3}})$ and $n = \Omega(1)$. This leads to

$$ \|\tilde{r} - \beta(r - \alpha e)\|_{\infty} \lesssim \frac{n}{\log n} \delta + \frac{\sqrt{n}}{\eta\sqrt{p}} \lesssim \frac{n}{\log n} \delta + \frac{n^{5/6}}{(\log n)^{1/3}} \delta^{1/3} \lesssim \frac{n}{\log n} \delta. $$

## 4.2 Analysis of Algorithm 2 (SVD-NRS)

We now analyze Algorithm 2, namely SVD-NRS, under the ERO model. The random matrix $\bar{D}^{-1/2} H \bar{D}^{-1/2}$ will concentrate around $(\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2}$. Recall from (4.1) that $\mathbb{E}[H] = \eta p (re^T - er^T)$, hence $(\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2}$ is a rank-2 matrix. It is easy to verify that

$$ \mathbb{E}[\bar{D}_{ii}] = \eta p \sum_{j=1}^{n} |r_i - r_j| + (1 - \eta) p \frac{M}{2}. $$

It will be useful to denote

$$ \max_i \mathbb{E}[\bar{D}_{ii}] = p \left[ \eta \max_i \left( \sum_{j=1}^n |r_i - r_j| \right) + (1 - \eta) \frac{M}{2} \right] = p \lambda_{\max}, \quad (4.11) $$

$$ \min_i \mathbb{E}[\bar{D}_{ii}] = p \left[ \eta \min_i \left( \sum_{j=1}^n |r_i - r_j| \right) + (1 - \eta) \frac{M}{2} \right] = p \lambda_{\min}. \quad (4.12) $$

Recall $\tilde{u}_2 \in \mathbb{R}^n$ obtained in Step 6 of SVD-NRS, and denote $u_2 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e)}{\| (\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e) \|_2}$. The following theorem bounds the $\ell_2$ error between $\tilde{u}_2$ and $u_2$ up to a global sign, its proof is given in Section 5.5.

**Theorem 9 ($\ell_2$ recovery, ERO model) Denote**

$$ \alpha = \frac{r^T (\mathbb{E}[\bar{D}])^{-1} e}{e^T (\mathbb{E}[\bar{D}])^{-1} e}, \quad A(\eta, M) = \eta M^2 + (1 - \eta) \frac{M}{2} \quad \text{and} \quad C_1(\eta, M) = 4A^{1/4}(\eta, M). $$

For $\varepsilon > 0$, denote

$$ \sigma_{\min} = \frac{\eta \|r - \alpha e\|_2 \sqrt{n}}{\lambda_{\max}}, \quad \sigma_{\max} = \frac{\eta \|r - \alpha e\|_2 \sqrt{n}}{\lambda_{\min}}, $$

$$ \tilde{\Delta} = 16M \sqrt{\frac{5}{3}pn} \frac{(2+\varepsilon)}{p\lambda_{\min}} + \frac{C_1(\eta, M)(np \log n)^{1/4}\sigma_{\max}}{p^{3/2}\lambda_{\min}^{3/2}} \left( \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right). $$

Let $\tilde{u}_2 \in \mathbb{R}^n$ be the vector obtained in Step 6 of SVD-NRS (with $\|\tilde{u}_2\|_2 = 1$) and denote $u_2 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e)}{\| (\mathbb{E}[\bar{D}])^{-1/2}(r-\alpha e) \|_2}$. If

$$ p \ge \frac{M^2}{9A(\eta, M)} \frac{\log n}{n}, \quad p \ge \frac{16(\sqrt{2}+1)^2 A(\eta, M) n \log n}{\lambda_{\min}^2} \quad \text{and} \quad \tilde{\Delta} \le \frac{\sigma_{\min}}{3} \quad (4.13) $$
hold, then there exists $\beta \in \{-1, 1\}$ and a constant $c_{\varepsilon} > 0$ depending only on $\varepsilon$ such that

$$||\tilde{u}_2 - \beta u_2||_2^2 \le \frac{15\tilde{\Delta}}{\sigma_{\min}}$$

holds with probability at least $1 - \frac{2}{n} - 2n \exp(-\frac{20pn}{3c_{\varepsilon}})$.

The quantities $\sigma_{\min}, \sigma_{\max}$ are lower and upper bounds on the non-zero singular values of the rank 2 matrix $(\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$, while $\tilde{\Delta}$ is an upper bound on the quantity $||\bar{D}^{-1/2}H\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}||_2$. The theorem essentially states that if $p$ is sufficiently large, then $\tilde{u}_2$ is close to $u_2$ up to a sign flip. In order to get a sense of the precise scaling of the quantities involved, we consider the special case $r_i = i$.

**Example 5** Consider $r_i = i$ for all $i$ with $\eta \in (0, 1)$ fixed, hence $M = n$. Denoting $S_i = \sum_{j=1}^n |r_i - r_j|$, one can verify that $S_i = i^2 - i(n + 1) + \frac{n^2+n}{2}$. This implies that $\lambda_{\max}, \lambda_{\min} = \Theta(\eta n^2) = A(\eta, M)$, and so, $C(\eta, M) = \Theta(\eta^{1/4}n^{1/2})$. Now let us find the scaling of $\alpha$. Since $\mathbb{E}[\bar{D}_{ii}] = \Theta(p\eta n^2)$, we obtain

$$\alpha = \frac{\sum_{i=1}^{n} i (\mathbb{E}[\bar{D}_{ii}])^{-1}}{\sum_{i=1}^{n} (\mathbb{E}[\bar{D}_{ii}])^{-1}} = \Theta\left(\frac{\sum_{i=1}^{n} i}{n}\right) = \Theta(n).$$

This in particular implies that $||r - \alpha e||_2^2 = \sum_{j=1}^n (j - \alpha)^2 = \Theta(n^3)$. We then obtain $\sigma_{\min}, \sigma_{\max} = \Theta(1)$; plugging the aforementioned scalings in the expression for $\tilde{\Delta}$ leads to

$$\tilde{\Delta} = O\left(\frac{1}{\eta\sqrt{np}} + \frac{(\log n)^{1/4}}{(\eta p)^{5/4}n^{9/4}}\left(\frac{(\log n)^{1/4}}{(\eta np)^{1/4}} + 2\sqrt{2}\right)\right).$$

The first two conditions on $p$ in (4.13) translate to $p \gtrsim \frac{1}{\eta} \frac{\log n}{n}$. Since $\sigma_{\min} = \Theta(1)$, hence the third condition in (4.13) holds provided $\tilde{\Delta} \lesssim 1$. In fact, for any $\delta \in (0, 1)$, we have $\tilde{\Delta} \lesssim \frac{\delta}{\sqrt{\log n}}$ if $p$ additionally satisfies

$$\frac{1}{\eta\sqrt{np}} \lesssim \frac{\delta}{\sqrt{\log n}} \Leftrightarrow p \gtrsim \frac{\log n}{\eta^2\delta^2n} \quad \text{and} \quad \frac{(\log n)^{1/4}}{(\eta p)^{5/4}n^{9/4}} \lesssim \frac{\delta}{\sqrt{\log n}} \Leftrightarrow p \gtrsim \frac{(\log n)^{3/5}}{\eta n^{9/5}\delta^{4/5}}.$$

Clearly, the first condition on $p$ dominates since $\frac{\log n}{\eta^2\delta^2n} \gtrsim \frac{(\log n)^{3/5}}{\eta n^{9/5}\delta^{4/5}}$. To summarize, Theorem 9 states that if $p$ satisfies

$$p \gtrsim \frac{\log n}{\eta^2\delta^2n}, \qquad (4.14)$$

then $\exists \beta \in \{-1, 1\}$ such that $||\tilde{u}_2 - \beta u_2||_2^2 \lesssim \delta$.

**Recovering scores in the $\ell_2$ norm.** The following theorem bounds the $\ell_2$ error between the score vector estimate $\tilde{r}$ defined in (3.7), and $r$ up to a global sign and shift. Its proof is outlined in Section 5.6. Note that both $\tilde{r}$ and $r - (\alpha + \alpha')e$ are centered vectors.
hold, then there exists $\beta \in \{-1, 1\}$ and a constant $c_{\varepsilon} > 0$ depending only on $\varepsilon$ such that

$$||\tilde{u}_2 - \beta u_2||_2^2 \le \frac{15\tilde{\Delta}}{\sigma_{\min}}$$

holds with probability at least $1 - \frac{2}{n} - 2n \exp(-\frac{20pn}{3c_{\varepsilon}})$.

The quantities $\sigma_{\min}, \sigma_{\max}$ are lower and upper bounds on the non-zero singular values of the rank 2 matrix $(\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$, while $\tilde{\Delta}$ is an upper bound on the quantity $||\bar{D}^{-1/2}H\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}||_2$. The theorem essentially states that if $p$ is sufficiently large, then $\tilde{u}_2$ is close to $u_2$ up to a sign flip. In order to get a sense of the precise scaling of the quantities involved, we consider the special case $r_i = i$.

**Example 5** Consider $r_i = i$ for all $i$ with $\eta \in (0, 1)$ fixed, hence $M = n$. Denoting $S_i = \sum_{j=1}^n |r_i - r_j|$, one can verify that $S_i = i^2 - i(n + 1) + \frac{n^2+n}{2}$. This implies that $\lambda_{\max}, \lambda_{\min} = \Theta(\eta n^2) = A(\eta, M)$, and so, $C(\eta, M) = \Theta(\eta^{1/4}n^{1/2})$. Now let us find the scaling of $\alpha$. Since $\mathbb{E}[\bar{D}_{ii}] = \Theta(p\eta n^2)$, we obtain

$$\alpha = \frac{\sum_{i=1}^{n} i (\mathbb{E}[\bar{D}_{ii}])^{-1}}{\sum_{i=1}^{n} (\mathbb{E}[\bar{D}_{ii}])^{-1}} = \Theta\left(\frac{\sum_{i=1}^{n} i}{n}\right) = \Theta(n).$$

This in particular implies that $||r - \alpha e||_2^2 = \sum_{j=1}^n (j - \alpha)^2 = \Theta(n^3)$. We then obtain $\sigma_{\min}, \sigma_{\max} = \Theta(1)$; plugging the aforementioned scalings in the expression for $\tilde{\Delta}$ leads to

$$\tilde{\Delta} = O\left(\frac{1}{\eta\sqrt{np}} + \frac{(\log n)^{1/4}}{(\eta p)^{5/4}n^{9/4}}\left(\frac{(\log n)^{1/4}}{(\eta np)^{1/4}} + 2\sqrt{2}\right)\right).$$

The first two conditions on $p$ in (4.13) translate to $p \gtrsim \frac{1}{\eta} \frac{\log n}{n}$. Since $\sigma_{\min} = \Theta(1)$, hence the third condition in (4.13) holds provided $\tilde{\Delta} \lesssim 1$. In fact, for any $\delta \in (0, 1)$, we have $\tilde{\Delta} \lesssim \frac{\delta}{\sqrt{\log n}}$ if $p$ additionally satisfies

$$\frac{1}{\eta\sqrt{np}} \lesssim \frac{\delta}{\sqrt{\log n}} \Leftrightarrow p \gtrsim \frac{\log n}{\eta^2\delta^2n} \quad \text{and} \quad \frac{(\log n)^{1/4}}{(\eta p)^{5/4}n^{9/4}} \lesssim \frac{\delta}{\sqrt{\log n}} \Leftrightarrow p \gtrsim \frac{(\log n)^{3/5}}{\eta n^{9/5}\delta^{4/5}}.$$

Clearly, the first condition on $p$ dominates since $\frac{\log n}{\eta^2\delta^2n} \gtrsim \frac{(\log n)^{3/5}}{\eta n^{9/5}\delta^{4/5}}$. To summarize, Theorem 9 states that if $p$ satisfies

$$p \gtrsim \frac{\log n}{\eta^2\delta^2n}, \qquad (4.14)$$

then $\exists \beta \in \{-1, 1\}$ such that $||\tilde{u}_2 - \beta u_2||_2^2 \lesssim \delta$.

**Recovering scores in the $\ell_2$ norm.** The following theorem bounds the $\ell_2$ error between the score vector estimate $\tilde{r}$ defined in (3.7), and $r$ up to a global sign and shift. Its proof is outlined in Section 5.6. Note that both $\tilde{r}$ and $r - (\alpha + \alpha')e$ are centered vectors.
**Theorem 10** Recall $\tilde{r} \in \mathbb{R}^n$ as defined in (3.7). Under the notation and assumptions of Theorem 9, and denoting $\alpha' = \frac{e^T(r-\alpha e)}{n}$, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \le \frac{2}{\eta p} \left[ \sqrt{\frac{3}{n}} (\sigma_{\max} + \tilde{\Delta}) (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4} \left( \sqrt{p\lambda_{\max}} + \frac{\lambda_{\max}}{\lambda_{\min}} \right) \right. \\ \left. + \frac{\tilde{\Delta}p\lambda_{\max}}{\sqrt{n}} + \sqrt{\frac{15\tilde{\Delta}}{\sigma_{\min}n}} \sigma_{\max}p\lambda_{\max} \right] $$

holds with probability at least $1 - \frac{2}{n} - 2n \exp(-\frac{20pn}{3c_\epsilon})$.

**Example 6** Consider the scenario where $r_i = i$ for each $i$. Then as discussed above for Theorem 9, if $p$ satisfies (4.14), we obtain the bound

$$ \begin{aligned} \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 &\lesssim \frac{1}{\eta p} \left[ \frac{(\eta n^3 p \log n)^{1/4}}{\sqrt{n}} (\sqrt{p\eta n^2} + 1) + \frac{p\eta n^2 \delta}{\sqrt{n} \log n} + \sqrt{\frac{\delta}{n\sqrt{\log n}}} p\eta n^2 \right] \\ &\lesssim \frac{1}{\eta p} \left[ \eta^{3/4} n^{5/4} (\log n)^{1/4} p^{3/4} + \eta p \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \right] \\ &= \frac{n^{5/4} (\log n)^{1/4}}{(\eta p)^{1/4}} + \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \\ &\lesssim n^{3/2} \eta^{1/4} \delta^{1/2} + \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \quad (\text{using (4.14)}). \end{aligned} $$

# 5. Proofs

## 5.1 Proof of Theorem 5

We now outline the proof of Theorem 5 by highlighting the steps involved, along with the intermediate Lemmas (the proofs of which are in Appendix C). The proof is broken into the following steps.

**Step 1: Singular values and singular vectors of C.** We begin by finding the singular values and singular vectors of C.

**Lemma 11** For $C = re^T - er^T$, it holds true that $C = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T$, with $\sigma_1 = \sigma_2 = \|r - \alpha e\|_2 \sqrt{n}$, and

$$ u_1 = v_2 = \frac{e}{\sqrt{n}}; \quad v_1 = -\frac{r - \alpha e}{\|r - \alpha e\|_2}; \quad u_2 = \frac{r - \alpha e}{\|r - \alpha e\|_2}, $$

where $\alpha = \frac{r^T e}{n}$.

**Step 2: Bounding the spectral norm of Z.** Recall that $(Z_{ij})_{i \le j}$ are independent zero-mean random variables with $Z_{ii} = 0$ and moreover, $Z_{ji} = -Z_{ij}$. Denote $Z = S - S^T$, where
**Theorem 10** Recall $\tilde{r} \in \mathbb{R}^n$ as defined in (3.7). Under the notation and assumptions of Theorem 9, and denoting $\alpha' = \frac{e^T(r-\alpha e)}{n}$, there exists $\beta \in \{-1, 1\}$ such that

$$ \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \le \frac{2}{\eta p} \left[ \sqrt{\frac{3}{n}} (\sigma_{\max} + \tilde{\Delta}) (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4} \left( \sqrt{p\lambda_{\max}} + \frac{\lambda_{\max}}{\lambda_{\min}} \right) \right. \\ \left. + \frac{\tilde{\Delta}p\lambda_{\max}}{\sqrt{n}} + \sqrt{\frac{15\tilde{\Delta}}{\sigma_{\min}n}} \sigma_{\max}p\lambda_{\max} \right] $$

holds with probability at least $1 - \frac{2}{n} - 2n \exp(-\frac{20pn}{3c_\epsilon})$.

**Example 6** Consider the scenario where $r_i = i$ for each $i$. Then as discussed above for Theorem 9, if $p$ satisfies (4.14), we obtain the bound

$$ \begin{aligned} \|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 &\lesssim \frac{1}{\eta p} \left[ \frac{(\eta n^3 p \log n)^{1/4}}{\sqrt{n}} (\sqrt{p\eta n^2} + 1) + \frac{p\eta n^2 \delta}{\sqrt{n} \log n} + \sqrt{\frac{\delta}{n\sqrt{\log n}}} p\eta n^2 \right] \\ &\lesssim \frac{1}{\eta p} \left[ \eta^{3/4} n^{5/4} (\log n)^{1/4} p^{3/4} + \eta p \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \right] \\ &= \frac{n^{5/4} (\log n)^{1/4}}{(\eta p)^{1/4}} + \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \\ &\lesssim n^{3/2} \eta^{1/4} \delta^{1/2} + \delta^{1/2} \frac{n^{3/2}}{(\log n)^{1/4}} \quad (\text{using (4.14)}). \end{aligned} $$

# 5. Proofs

## 5.1 Proof of Theorem 5

We now outline the proof of Theorem 5 by highlighting the steps involved, along with the intermediate Lemmas (the proofs of which are in Appendix C). The proof is broken into the following steps.

**Step 1: Singular values and singular vectors of C.** We begin by finding the singular values and singular vectors of C.

**Lemma 11** For $C = re^T - er^T$, it holds true that $C = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T$, with $\sigma_1 = \sigma_2 = \|r - \alpha e\|_2 \sqrt{n}$, and

$$ u_1 = v_2 = \frac{e}{\sqrt{n}}; \quad v_1 = -\frac{r - \alpha e}{\|r - \alpha e\|_2}; \quad u_2 = \frac{r - \alpha e}{\|r - \alpha e\|_2}, $$

where $\alpha = \frac{r^T e}{n}$.

**Step 2: Bounding the spectral norm of Z.** Recall that $(Z_{ij})_{i \le j}$ are independent zero-mean random variables with $Z_{ii} = 0$ and moreover, $Z_{ji} = -Z_{ij}$. Denote $Z = S - S^T$, where
$S$ is an upper-triangular matrix with $S_{ij} = Z_{ij}$ for $i \le j$. Note that $||Z||_2 \le ||S||_2 + ||S^T||_2 \le 2||S||_2$. We define a symmetric matrix

$$ \tilde{S} := \begin{bmatrix} 0 & S \\ S^T & 0 \end{bmatrix} \in \mathbb{R}^{2n \times 2n}. $$

It is easy to check that $||\tilde{S}||_2 = ||S||_2$, and thus $||Z||_2 \le 2||\tilde{S}||_2$. Since $\tilde{S}$ is a random symmetric matrix with the entries $(\tilde{S}_{ij})_{i \le j}$ being centered, independent and bounded random variables, we can bound $||\tilde{S}||_2$ using a recent result of (Bandeira and van Handel, 2016, Corollary 3.12, Remark 3.13) that gives non-asymptotic bounds on the spectral norm of such random symmetric matrices. This altogether leads to the following Lemma.

**Lemma 12** For $Z \in \mathbb{R}^{n \times n}$ as defined in (4.2), we have that

$$ ||Z||_2 \le 8M \sqrt{\frac{5}{3}} pn(2 + \varepsilon) \qquad (5.1) $$

holds with probability at least $1 - 2n \exp(-\frac{20pn}{3c_{\varepsilon}})$.

**Step 3: Using Wedin's bound.** Let us denote $U = [u_1 u_2] \in \mathbb{R}^{n \times 2}$ (recall $u_1, u_2$ from Lemma 11). Denote $\hat{U} = [\hat{u}_1 \hat{u}_2] \in \mathbb{R}^{n \times 2}$ where $\hat{u}_1, \hat{u}_2$ are the left singular vectors corresponding to the top two singular values of $H$. The following Lemma states that $\mathcal{R}(U)$ is close to $\mathcal{R}(\hat{U})$ if $||Z||_2$ is small. The proof uses Wedin's bound (Wedin, 1972) for perturbation of singular subspaces.

**Lemma 13** Given $H = \eta p C + Z$, where $||Z||_2 \le \Delta$, if $\Delta < \eta p ||r - \alpha e||_2 \sqrt{n}$, then

$$ ||(I - \hat{U}\hat{U}^T)U||_2 \le \frac{\Delta}{\eta p ||r - \alpha e||_2 \sqrt{n} - \Delta} (=: \delta). \qquad (5.2) $$

**Step 4: Analyzing the projection step.** Next, we project $u_1$ on $\text{span}(\hat{U})$ to obtain $\tilde{u}_1$, find a unit vector orthogonal to $\tilde{u}_1$ lying in $\mathcal{R}(\hat{U})$ (denote this by $\tilde{u}_2$), and show that $\tilde{u}_2$ is close to $u_2$ up to a sign flip.

**Lemma 14** With $\delta$ as defined in Lemma 13 let $\delta \le 1/2$ hold. Then, there exists $\beta \in \{-1, 1\}$ such that

$$ ||\tilde{u}_2 - \beta u_2||_2^2 \le 10\delta. $$

**Step 5: Putting it together.** We now use the above Lemma's to readily obtain Theorem 5 (see appendix for details). Using $||Z||_2 \le \Delta = 8M \sqrt{\frac{5}{3}} pn(2 + \varepsilon)$ in the expression for $\delta$, one can verify that

$$ \delta \le 1/2 \Leftrightarrow 3\Delta \le \eta p ||r - \alpha e||_2 \sqrt{n} \qquad (5.3) $$

holds if $||r - \alpha e||_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}} (2 + \varepsilon)$. By using the bound on $\Delta$ in (5.3) in the expression for $\delta$, we finally obtain

$$ \delta \le \frac{12M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2 + \varepsilon)}{||r - \alpha e||_2}. $$

Plugging this in Lemma 14 yields the stated $\ell_2$ bound and completes the proof.
$S$ is an upper-triangular matrix with $S_{ij} = Z_{ij}$ for $i \le j$. Note that $||Z||_2 \le ||S||_2 + ||S^T||_2 \le 2||S||_2$. We define a symmetric matrix

$$ \tilde{S} := \begin{bmatrix} 0 & S \\ S^T & 0 \end{bmatrix} \in \mathbb{R}^{2n \times 2n}. $$

It is easy to check that $||\tilde{S}||_2 = ||S||_2$, and thus $||Z||_2 \le 2||\tilde{S}||_2$. Since $\tilde{S}$ is a random symmetric matrix with the entries $(\tilde{S}_{ij})_{i \le j}$ being centered, independent and bounded random variables, we can bound $||\tilde{S}||_2$ using a recent result of (Bandeira and van Handel, 2016, Corollary 3.12, Remark 3.13) that gives non-asymptotic bounds on the spectral norm of such random symmetric matrices. This altogether leads to the following Lemma.

**Lemma 12** For $Z \in \mathbb{R}^{n \times n}$ as defined in (4.2), we have that

$$ ||Z||_2 \le 8M \sqrt{\frac{5}{3}} pn(2 + \varepsilon) \quad (5.1) $$

holds with probability at least $1 - 2n \exp(-\frac{20pn}{3c_\varepsilon})$.

**Step 3: Using Wedin's bound.** Let us denote $U = [u_1 u_2] \in \mathbb{R}^{n \times 2}$ (recall $u_1, u_2$ from Lemma 11). Denote $\hat{U} = [\hat{u}_1 \hat{u}_2] \in \mathbb{R}^{n \times 2}$ where $\hat{u}_1, \hat{u}_2$ are the left singular vectors corresponding to the top two singular values of $H$. The following Lemma states that $\mathcal{R}(U)$ is close to $\mathcal{R}(\hat{U})$ if $||Z||_2$ is small. The proof uses Wedin's bound (Wedin, 1972) for perturbation of singular subspaces.

**Lemma 13** Given $H = \eta p C + Z$, where $||Z||_2 \le \Delta$, if $\Delta < \eta p ||r - \alpha e||_2 \sqrt{n}$, then

$$ ||(I - \hat{U}\hat{U}^T)U||_2 \le \frac{\Delta}{\eta p ||r - \alpha e||_2 \sqrt{n} - \Delta} (=: \delta). \quad (5.2) $$

**Step 4: Analyzing the projection step.** Next, we project $u_1$ on span($\hat{U}$) to obtain $\tilde{u}_1$, find a unit vector orthogonal to $\tilde{u}_1$ lying in $\mathcal{R}(\hat{U})$ (denote this by $\tilde{u}_2$), and show that $\tilde{u}_2$ is close to $u_2$ up to a sign flip.

**Lemma 14** With $\delta$ as defined in Lemma 13 let $\delta \le 1/2$ hold. Then, there exists $\beta \in \{-1, 1\}$ such that

$$ ||\tilde{u}_2 - \beta u_2||_2^2 \le 10\delta. $$

**Step 5: Putting it together.** We now use the above Lemma's to readily obtain Theorem 5 (see appendix for details). Using $||Z||_2 \le \Delta = 8M \sqrt{\frac{5}{3}} pn(2 + \varepsilon)$ in the expression for $\delta$, one can verify that

$$ \delta \le 1/2 \Leftrightarrow 3\Delta \le \eta p ||r - \alpha e||_2 \sqrt{n} \quad (5.3) $$

holds if $||r - \alpha e||_2 \ge \frac{24M}{\eta} \sqrt{\frac{5}{3p}} (2 + \varepsilon)$. By using the bound on $\Delta$ in (5.3) in the expression for $\delta$, we finally obtain

$$ \delta \le \frac{12M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2 + \varepsilon)}{||r - \alpha e||_2}. $$

Plugging this in Lemma 14 yields the stated $\ell_2$ bound and completes the proof.
## 5.2 Proof of Theorem 6

Recall that $U = [u_1 \ u_2]$, $\hat{U} = [\hat{u}_1 \ \hat{u}_2] \in \mathbb{R}^{n \times 2}$ correspond to the two largest singular vectors of $\mathbb{E}[H]$, and $H$ respectively. The proof is broken down into three main steps. The first step involves bounding $\|\tilde{u}_2 - \beta u_2\|_{\infty}$ in terms of $\|\hat{U} - UO^*\|_{\max}$ – where $O^*$ is a $2 \times 2$ orthogonal matrix that “aligns” $U$ with $\hat{U}$ – provided $\|Z\|_2$ is small. The second step involves bounding $\|\hat{U} - UO^*\|_{\max}$, and in the third step we combine the results from the previous steps. The proofs of all lemmas are provided in Appendix D.

**Step 1: Bounding** $\|\tilde{u}_2 - \beta u_2\|_{\infty}$ in terms of $\|\hat{U} - UO^*\|_{\max}$. The following Lemma states this formally.

**Lemma 15** For $\Delta$ as defined in Lemma 13, let $\Delta \le \frac{\eta p \|r - \alpha e\|_2 \sqrt{n}}{3}$ hold. Then for any orthogonal matrix $O \in \mathbb{R}^{2 \times 2}$, there exists $\beta \in \{-1, 1\}$ depending on $O$ such that

$$ \|\tilde{u}_2 - \beta u_2\|_{\infty} \le 4\|\hat{U} - UO\|_{\max} \left(2 + \frac{\sqrt{n}(M-\alpha)}{\|r-\alpha e\|_2}\right) + 4\sqrt{n}\|\hat{U} - UO\|_{\max}^2. \quad (5.4) $$

Moreover, there exists an orthogonal matrix $O^* \in \mathbb{R}^{2 \times 2}$ such that

$$ \|\hat{U} - UO^*\|_2 \le \frac{3\Delta}{\eta p \sqrt{n} \|r - \alpha e\|_2}. \quad (5.5) $$

**Step 2: Bounding** $\|\hat{U} - UO^*\|_{\max}$. Since $H^T = -H$, therefore $-H^2 = HH^T$. Using the SVD $H = \hat{U}\hat{\Sigma}\hat{V}$, we arrive at

$$ \begin{aligned} HH^T \hat{U} &= -H^2 \hat{U} = \hat{U}\hat{\Sigma}^2 \\ \Rightarrow \hat{U} &= -H^2 \hat{U}\hat{\Sigma}^{-2} = -(\mathbb{E}[H]^2 + \mathbb{E}[H]Z + Z\mathbb{E}[H] + Z^2)\hat{U}\hat{\Sigma}^{-2}. \end{aligned} \quad (5.6) $$

For convenience of notation, denote $\sigma = \sigma_i(\mathbb{E}[H]) = \eta p \sqrt{n} \|r - \alpha e\|_2$; $i = 1, 2$. Since $\mathbb{E}[H]^T = -\mathbb{E}[H]$, we get

$$ \begin{aligned} -\mathbb{E}[H]^2 &= \mathbb{E}[H]\mathbb{E}[H]^T = U\Sigma^2 U^T = \sigma^2(UO^*)(UO^*)^T \\ UO^* &= -\frac{1}{\sigma^2}(\mathbb{E}[H])^2(UO^*). \end{aligned} \quad (5.7) $$

Subtracting (5.7) from (5.6) leads to

$$ \begin{aligned} \hat{U} - UO^* &= \underbrace{-\mathbb{E}[H]^2(\hat{U} - UO^*)\hat{\Sigma}^{-2}}_{E_1} + \underbrace{\mathbb{E}[H]^2UO^*(\sigma^{-2}I - \hat{\Sigma}^{-2})}_{E_2} - \underbrace{\mathbb{E}[H]Z\hat{U}\hat{\Sigma}^{-2}}_{E_3} \\ &\quad - \underbrace{Z\mathbb{E}[H]\hat{U}\hat{\Sigma}^{-2}}_{E_4} - \underbrace{Z^2\hat{U}\hat{\Sigma}^{-2}}_{E_5} \\ &= E_1 + E_2 - E_3 - E_4 - E_5 \end{aligned} $$

$$ \Rightarrow \|\hat{U} - UO^*\|_{\max} \le \sum_{i=1}^{5} \|E_i\|_{\max}. \quad (5.8) $$

We now proceed to bound each term in the RHS of (5.8). This is stated precisely in the following Lemma.
## 5.2 Proof of Theorem 6

Recall that $U = [u_1 \ u_2]$, $\hat{U} = [\hat{u}_1 \ \hat{u}_2] \in \mathbb{R}^{n \times 2}$ correspond to the two largest singular vectors of $\mathbb{E}[H]$, and $H$ respectively. The proof is broken down into three main steps. The first step involves bounding $\|\tilde{u}_2 - \beta u_2\|_{\infty}$ in terms of $\|\hat{U} - UO^*\|_{\max}$ – where $O^*$ is a $2 \times 2$ orthogonal matrix that “aligns” $U$ with $\hat{U}$ – provided $\|Z\|_2$ is small. The second step involves bounding $\|\hat{U} - UO^*\|_{\max}$, and in the third step we combine the results from the previous steps. The proofs of all lemmas are provided in Appendix D.

**Step 1: Bounding** $\|\tilde{u}_2 - \beta u_2\|_{\infty}$ in terms of $\|\hat{U} - UO^*\|_{\max}$. The following Lemma states this formally.

**Lemma 15** For $\Delta$ as defined in Lemma 13, let $\Delta \le \frac{\eta p \|r - \alpha e\|_2 \sqrt{n}}{3}$ hold. Then for any orthogonal matrix $O \in \mathbb{R}^{2 \times 2}$, there exists $\beta \in \{-1, 1\}$ depending on $O$ such that

$$ \|\tilde{u}_2 - \beta u_2\|_{\infty} \le 4\|\hat{U} - UO\|_{\max} \left(2 + \frac{\sqrt{n}(M-\alpha)}{\|r-\alpha e\|_2}\right) + 4\sqrt{n}\|\hat{U} - UO\|_{\max}^2. \quad (5.4) $$

Moreover, there exists an orthogonal matrix $O^* \in \mathbb{R}^{2 \times 2}$ such that

$$ \|\hat{U} - UO^*\|_2 \le \frac{3\Delta}{\eta p \sqrt{n} \|r - \alpha e\|_2}. \quad (5.5) $$

**Step 2: Bounding** $\|\hat{U} - UO^*\|_{\max}$. Since $H^T = -H$, therefore $-H^2 = HH^T$. Using the SVD $H = \hat{U}\hat{\Sigma}\hat{V}$, we arrive at

$$ \begin{aligned} HH^T \hat{U} &= -H^2 \hat{U} = \hat{U}\hat{\Sigma}^2 \\ \Rightarrow \hat{U} &= -H^2 \hat{U}\hat{\Sigma}^{-2} = -(\mathbb{E}[H]^2 + \mathbb{E}[H]Z + Z\mathbb{E}[H] + Z^2)\hat{U}\hat{\Sigma}^{-2}. \end{aligned} \quad (5.6) $$

For convenience of notation, denote $\sigma = \sigma_i(\mathbb{E}[H]) = \eta p \sqrt{n} \|r - \alpha e\|_2$; $i = 1, 2$. Since $\mathbb{E}[H]^T = -\mathbb{E}[H]$, we get

$$ \begin{aligned} -\mathbb{E}[H]^2 &= \mathbb{E}[H]\mathbb{E}[H]^T = U\Sigma^2 U^T = \sigma^2(UO^*)(UO^*)^T \\ UO^* &= -\frac{1}{\sigma^2}(\mathbb{E}[H])^2(UO^*). \end{aligned} \quad (5.7) $$

Subtracting (5.7) from (5.6) leads to

$$ \begin{aligned} \hat{U} - UO^* &= \underbrace{-\mathbb{E}[H]^2(\hat{U} - UO^*)\hat{\Sigma}^{-2}}_{E_1} + \underbrace{\mathbb{E}[H]^2UO^*(\sigma^{-2}I - \hat{\Sigma}^{-2})}_{E_2} - \underbrace{\mathbb{E}[H]Z\hat{U}\hat{\Sigma}^{-2}}_{E_3} \\ &\quad - \underbrace{Z\mathbb{E}[H]\hat{U}\hat{\Sigma}^{-2}}_{E_4} - \underbrace{Z^2\hat{U}\hat{\Sigma}^{-2}}_{E_5} \\ &= E_1 + E_2 - E_3 - E_4 - E_5 \end{aligned} $$

$$ \Rightarrow \|\hat{U} - UO^*\|_{\max} \le \sum_{i=1}^{5} \|E_i\|_{\max}. \quad (5.8) $$

We now proceed to bound each term in the RHS of (5.8). This is stated precisely in the following Lemma.
**Lemma 16** If $\Delta \le \frac{\sigma}{3} = \frac{\eta p \sqrt{n} \|r-\alpha e\|_2}{3}$, then the following holds true.

$$1. \quad \|E_1\|_{\max} \le \frac{27\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$2. \quad \|E_2\|_{\max} \le \frac{21\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$3. \quad \|E_3\|_{\max} \le \frac{9\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$4. \quad \|E_4\|_{\max} \le \frac{9}{4\sigma} (\|Zu_1\|_{\infty} + \|Zu_2\|_{\infty}).$$

$$5. \quad \|E_5\|_{\max} \le \frac{27\Delta^3}{4\sigma^3} + \frac{9}{4\sigma^2} (\|Z^2u_1\|_{\infty} + \|Z^2u_2\|_{\infty}).$$

Next, we will bound $\|Zu_i\|_{\infty}$ and $\|Z^2u_i\|_{\infty}$ for $i = 1, 2$. Bounding $\|Zu_i\|_{\infty}$ is a relatively straightforward consequence of Bernstein's inequality (Boucheron et al., 2013, Corollary 2.11), we state this in the following Lemma.

**Lemma 17** If $p \ge \frac{2\log n}{15n}$ then the following holds.

$$1. \quad \mathbb{P}(\|Zu_1\|_{\infty} \ge \frac{2\sqrt{2}+4}{3} M \sqrt{15p \log n}) \le \frac{2}{n}.$$

$$2. \quad \mathbb{P}(\|Zu_2\|_{\infty} \ge \frac{2\sqrt{2}+4}{3} M \sqrt{15pBn \log n}) \le \frac{2}{n} \quad \text{where } B = \frac{(M-\alpha)^2}{\|r-\alpha e\|_2^2}.$$

Bounding $\|Z^2u_i\|_{\infty}$ is a considerably more challenging task since the entries within a row of $Z^2$ will not be independent now. However, leveraging a recent entry-wise concentration result for the product of a random matrix (raised to a power) and a fixed vector (Eldridge et al., 2018, Theorem 15) (see Theorem 25), we are able to do so and arrive at the following Lemma.

**Lemma 18** Let $p \ge \frac{1}{2n}$ hold. Choose $\xi > 1, 0 < \kappa < 1$ and define $\mu = \frac{2}{\kappa+1}$. Then if $\frac{16}{\kappa} \le (\log n)^\xi$, the following holds true.

$$1. \quad \mathbb{P}(\|Z^2u_1\|_{\infty} \ge 8\sqrt{np}M^2(\log n)^{2\xi}) \le n^{1-\frac{1}{4}(\log_\mu n)^{\xi-1}(\log_\mu e)^{-\xi}}.$$

$$2. \quad \mathbb{P}(\|Z^2u_2\|_{\infty} \ge 8npM^2 \frac{M-\alpha}{\|r-\alpha e\|_2} (\log n)^{2\xi}) \le n^{1-\frac{1}{4}(\log_\mu n)^{\xi-1}(\log_\mu e)^{-\xi}}.$$

**Step 3: Putting it together.** We now combine the results of the preceding Lemmas to derive the final approximation bound. Recall from Lemma 12 that $\Delta = 8M\sqrt{\frac{5}{3}pn}(2+\varepsilon)$. Also, $\sigma = \eta p \sqrt{n} \|r-\alpha e\|_2$ and so

$$\frac{\Delta}{\sigma} = \frac{8M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2+\varepsilon)}{\|r-\alpha e\|_2} \quad (5.9)$$
**Lemma 16** If $\Delta \le \frac{\sigma}{3} = \frac{\eta p \sqrt{n} \|r-\alpha e\|_2}{3}$, then the following holds true.

$$1. \quad \|E_1\|_{\max} \le \frac{27\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$2. \quad \|E_2\|_{\max} \le \frac{21\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$3. \quad \|E_3\|_{\max} \le \frac{9\Delta}{4\sigma} \left( \frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2} \right).$$

$$4. \quad \|E_4\|_{\max} \le \frac{9}{4\sigma} (\|Zu_1\|_{\infty} + \|Zu_2\|_{\infty}).$$

$$5. \quad \|E_5\|_{\max} \le \frac{27\Delta^3}{4\sigma^3} + \frac{9}{4\sigma^2} (\|Z^2u_1\|_{\infty} + \|Z^2u_2\|_{\infty}).$$

Next, we will bound $\|Zu_i\|_{\infty}$ and $\|Z^2u_i\|_{\infty}$ for $i = 1, 2$. Bounding $\|Zu_i\|_{\infty}$ is a relatively straightforward consequence of Bernstein's inequality (Boucheron et al., 2013, Corollary 2.11), we state this in the following Lemma.

**Lemma 17** If $p \ge \frac{2\log n}{15n}$ then the following holds.

$$1. \quad \mathbb{P}(\|Zu_1\|_{\infty} \ge \frac{2\sqrt{2}+4}{3} M \sqrt{15p \log n}) \le \frac{2}{n}.$$

$$2. \quad \mathbb{P}(\|Zu_2\|_{\infty} \ge \frac{2\sqrt{2}+4}{3} M \sqrt{15pBn \log n}) \le \frac{2}{n} \quad \text{where } B = \frac{(M-\alpha)^2}{\|r-\alpha e\|_2^2}.$$

Bounding $\|Z^2u_i\|_{\infty}$ is a considerably more challenging task since the entries within a row of $Z^2$ will not be independent now. However, leveraging a recent entry-wise concentration result for the product of a random matrix (raised to a power) and a fixed vector (Eldridge et al., 2018, Theorem 15) (see Theorem 25), we are able to do so and arrive at the following Lemma.

**Lemma 18** Let $p \ge \frac{1}{2n}$ hold. Choose $\xi > 1, 0 < \kappa < 1$ and define $\mu = \frac{2}{\kappa+1}$. Then if $\frac{16}{\kappa} \le (\log n)^\xi$, the following holds true.

$$1. \quad \mathbb{P}(\|Z^2u_1\|_{\infty} \ge 8\sqrt{n}pM^2(\log n)^{2\xi}) \le n^{1-\frac{1}{4}(\log_\mu n)^{\xi-1}(\log_\mu e)^{-\xi}}.$$

$$2. \quad \mathbb{P}(\|Z^2u_2\|_{\infty} \ge 8npM^2 \frac{M-\alpha}{\|r-\alpha e\|_2} (\log n)^{2\xi}) \le n^{1-\frac{1}{4}(\log_\mu n)^{\xi-1}(\log_\mu e)^{-\xi}}.$$

**Step 3: Putting it together.** We now combine the results of the preceding Lemmas to derive the final approximation bound. Recall from Lemma 12 that $\Delta = 8M\sqrt{\frac{5}{3}pn}(2+\varepsilon)$. Also, $\sigma = \eta p \sqrt{n} \|r-\alpha e\|_2$ and so

$$\frac{\Delta}{\sigma} = \frac{8M}{\eta} \sqrt{\frac{5}{3p}} \frac{(2+\varepsilon)}{\|r-\alpha e\|_2} \quad (5.9)$$
Using Lemmas 16,17,18 and (5.9) in (5.8), we get

$$
\begin{align*}
\|\hat{U} - UO^*\|_{\max} &\le \left(\frac{114M}{\eta}\sqrt{\frac{5}{3p}}\frac{(2+\varepsilon)}{\|r-\alpha e\|_2}\right) \left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) \\
&\quad + \left(\frac{3\sqrt{15}(2\sqrt{2}+4)}{4}\right) \left(\frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r-\alpha e\|_2}\right) \left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) \\
&\quad + \frac{18M^2(\log n)^{2\xi}}{\eta^2 p \|r-\alpha e\|_2^2} \left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) + \frac{27}{4} \left(\frac{8\sqrt{5}}{3}\right)^3 \frac{M^3}{\eta^3 p^{3/2}} \frac{(2+\varepsilon)^3}{\|r-\alpha e\|_2^3} \\
&\le C_{\varepsilon} \left[ \left(\frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r-\alpha e\|_2} + \frac{M^2(\log n)^{2\xi}}{\eta^2 p \|r-\alpha e\|_2^2}\right) \left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) + \frac{M^3}{\eta^3 p^{3/2}} \frac{(2+\varepsilon)^3}{\|r-\alpha e\|_2^3} \right] \\
&= C(n, M, \eta, p, \varepsilon, r)
\end{align*}
$$

where $C_\varepsilon > 0$ is a universal constant depending only on $\varepsilon$. Plugging this in (5.4), we arrive
at the stated $\ell_\infty$ bound in the Theorem. The lower bound on the success probability follows
readily via a union bound on the events stated in Lemmas 12, 17, 18.

**5.3 Proof of Theorem 7**

We assume w.l.o.g that $\pi$ is identity so that $\pi(i) = i$. For any $i \in [n]$, we have for all rankings $\tilde{\pi}$ which are consistent with $\tilde{u}_2$ that

$$
|i - \tilde{\pi}(i)| = \sum_{j>i} \mathbf{1}_{\{\tilde{\pi}(j) < \tilde{\pi}(i)\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{\pi}(j) > \tilde{\pi}(i)\}} \le \sum_{j>i} \mathbf{1}_{\{\tilde{u}_{2,j} \ge \tilde{u}_{2,i}\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{u}_{2,j} \le \tilde{u}_{2,i}\}},
$$

leading to the bound

$$
\|\pi - \tilde{\pi}\|_{\infty} = \max_{i} |i - \tilde{\pi}(i)| \le \sum_{j>i} \mathbf{1}_{\{\tilde{u}_{2,j} \ge \tilde{u}_{2,i}\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{u}_{2,j} \le \tilde{u}_{2,i}\}}.
$$

Now for a given $i' \in [n]$, we can decompose $\tilde{u}_{2,i'} - \tilde{u}_{2,j}$ as

$$
\tilde{u}_{2,i'} - \tilde{u}_{2,j} = (\tilde{u}_{2,i'} - u_{2,i'}) + (u_{2,i'} - u_{2,j}) + (u_{2,j} - \tilde{u}_{2,j}),
$$

which in turn implies $\tilde{u}_{2,i'} - \tilde{u}_{2,j} \le (u_{2,i'} - u_{2,j}) + 2\|\tilde{u}_2 - u_2\|_{\infty}$. Hence if $j$ is such that
$i' > j$, then $\|\tilde{u}_2 - u_2\|_{\infty} \le \frac{|u_{2,i'} - u_{2,j}|}{2}$ implies

$$
\tilde{u}_{2,i'} - \tilde{u}_{2,j} \le \underbrace{(u_{2,i'} - u_{2,j})}_{\le 0} + |u_{2,i'} - u_{2,j}| = 0.
$$
Using Lemmas 16,17,18 and (5.9) in (5.8), we get

$$
\begin{align*}
\|\hat{U} - UO^*\|_{\max} &\le \left(\frac{114M}{\eta}\sqrt{\frac{5}{3p}}\frac{(2+\varepsilon)}{\|r-\alpha e\|_2}\right)\left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) \\
&\quad + \left(\frac{3\sqrt{15}(2\sqrt{2}+4)}{4}\right)\left(\frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r-\alpha e\|_2}\right)\left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) \\
&\quad + \frac{18M^2(\log n)^{2\xi}}{\eta^2 p\|r-\alpha e\|_2^2}\left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) + \frac{27}{4}\left(\frac{8\sqrt{5}}{3}\right)^3\frac{M^3}{\eta^3 p^{3/2}}\frac{(2+\varepsilon)^3}{\|r-\alpha e\|_2^3} \\
&\le C_{\varepsilon}\left[\left(\frac{M\sqrt{\log n}}{\eta\sqrt{p}\|r-\alpha e\|_2} + \frac{M^2(\log n)^{2\xi}}{\eta^2 p\|r-\alpha e\|_2^2}\right)\left(\frac{1}{\sqrt{n}} + \frac{M-\alpha}{\|r-\alpha e\|_2}\right) + \frac{M^3}{\eta^3 p^{3/2}}\frac{(2+\varepsilon)^3}{\|r-\alpha e\|_2^3}\right] \\
&= C(n, M, \eta, p, \varepsilon, r)
\end{align*}
$$

where $C_\varepsilon > 0$ is a universal constant depending only on $\varepsilon$. Plugging this in (5.4), we arrive
at the stated $\ell_\infty$ bound in the Theorem. The lower bound on the success probability follows
readily via a union bound on the events stated in Lemmas 12, 17, 18.

5.3 Proof of Theorem 7

We assume w.l.o.g that $\pi$ is identity so that $\pi(i) = i$. For any $i \in [n]$, we have for all rankings $\tilde{\pi}$ which are consistent with $\tilde{u}_2$ that

$$
|i - \tilde{\pi}(i)| = \sum_{j>i} \mathbf{1}_{\{\tilde{\pi}(j) < \tilde{\pi}(i)\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{\pi}(j) > \tilde{\pi}(i)\}} \le \sum_{j>i} \mathbf{1}_{\{\tilde{u}_{2,j} \ge \tilde{u}_{2,i}\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{u}_{2,j} \le \tilde{u}_{2,i}\}},
$$

leading to the bound

$$
\|\pi - \tilde{\pi}\|_{\infty} = \max_{i} |i - \tilde{\pi}(i)| \le \sum_{j>i} \mathbf{1}_{\{\tilde{u}_{2,j} \ge \tilde{u}_{2,i}\}} + \sum_{j<i} \mathbf{1}_{\{\tilde{u}_{2,j} \le \tilde{u}_{2,i}\}}.
$$

Now for a given $i' \in [n]$, we can decompose $\tilde{u}_{2,i'} - \tilde{u}_{2,j}$ as

$$
\tilde{u}_{2,i'} - \tilde{u}_{2,j} = (\tilde{u}_{2,i'} - u_{2,i'}) + (u_{2,i'} - u_{2,j}) + (u_{2,j} - \tilde{u}_{2,j}),
$$

which in turn implies $\tilde{u}_{2,i'} - \tilde{u}_{2,j} \le (u_{2,i'} - u_{2,j}) + 2\|\tilde{u}_2 - u_2\|_{\infty}$. Hence if $j$ is such that
$i' > j$, then $\|\tilde{u}_2 - u_2\|_{\infty} \le \frac{|u_{2,i'} - u_{2,j}|}{2}$ implies

$$
\tilde{u}_{2,i'} - \tilde{u}_{2,j} \le \underbrace{(u_{2,i'} - u_{2,j})}_{\le 0} + |u_{2,i'} - u_{2,j}| = 0.
$$
Using this, we obtain

$$
\begin{align*}
\sum_{j<i'} \mathbf{1}_{\{\tilde{u}_{2,j} < \tilde{u}_{2,i'}\}} &\le \sum_{j<i'} \mathbf{1}_{\{\|\tilde{u}_2 - u_2\|_{\infty} > \frac{|u_{2,i'} - u_{2,j}|}{2}\}} \\
&\le \sum_{j<i'} \mathbf{1}_{\{\Upsilon(n, M, \eta, p, \varepsilon, r) > \frac{|u_{2,i'} - u_{2,j}|}{2}\}} \quad \text{(since } \|\tilde{u}_2 - u_2\|_{\infty} \le \Upsilon(n, M, \eta, p, \varepsilon, r)\text{)} \\
&= \sum_{j<i'} \mathbf{1}_{\{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2 > |r_{i'} - r_j|\}} \\
&= \left\lfloor \frac{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2}{\rho} \right\rfloor \le \frac{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2}{\rho}.
\end{align*}
$$

The same bound holds for $\sum_{j>i'} \mathbf{1}_{\{\tilde{u}_{2,j} > \tilde{u}_{2,i'}\}}$ and hence the statement of the Theorem follows.

5.4 Proof of Theorem 8

Recall that $\tilde{w} = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}} \tilde{u}_2$. Let us denote

$$
w = \frac{\sigma}{\eta p \sqrt{n}} \beta u_2 = \|r - \alpha e\|_2 \beta u_2 = (r - \alpha e) \beta.
$$

For any norm $\|\cdot\|$ on $\mathbb{R}^n$, we have by triangle inequality that

$$
\|\tilde{r} - w\| \le \|\tilde{w} - w\| + \left| \frac{e^T \tilde{w}}{n} \right| \|e\|.
$$

Since $e^T w = 0$, therefore we obtain via Cauchy-Schwarz that

$$
\left| \frac{e^T \tilde{w}}{n} \right| = \left| \frac{e^T (\tilde{w} - w)}{n} \right| \le \frac{\| \tilde{w} - w \|_2}{\sqrt{n}}.
$$

Hence we have that

$$
\|\tilde{r} - w\|_2 \le 2\|\tilde{w} - w\|_2 \quad \text{and} \quad \|\tilde{r} - w\|_{\infty} \le \|\tilde{w} - w\|_{\infty} + \frac{\|\tilde{w} - w\|_2}{\sqrt{n}} \le 2\|\tilde{w} - w\|_{\infty}.
$$

It remains to bound $\|\tilde{w} - w\|_2$ and $\|\tilde{w} - w\|_{\infty}$. Our starting point will be

$$
\begin{align*}
\|\tilde{w} - w\| &= \frac{1}{\eta p \sqrt{n}} \|\hat{\sigma}_1 \tilde{u}_2 - \sigma_1 \beta u_2\| \\
&= \frac{1}{\eta p \sqrt{n}} \|(\hat{\sigma}_1 - \sigma_1) \tilde{u}_2 + \sigma_1 (\tilde{u}_2 - \beta u_2)\| \\
\Rightarrow \|\tilde{w} - w\| &\le \frac{1}{\eta p \sqrt{n}} (|\hat{\sigma}_1 - \sigma| \|\tilde{u}_2\| + \sigma_1 \|\tilde{u}_2 - \beta u_2\|) \\
&\le \frac{1}{\eta p \sqrt{n}} (\Delta \|\tilde{u}_2\| + \eta p \|r - \alpha e\|_2 \sqrt{n} \|\tilde{u}_2 - \beta u_2\|). \tag{5.10}
\end{align*}
$$
Using this, we obtain

$$
\begin{align*}
\sum_{j<i'} \mathbf{1}_{\{\tilde{u}_{2,j} < \tilde{u}_{2,i'}\}} &\le \sum_{j<i'} \mathbf{1}_{\{\|\tilde{u}_2 - u_2\|_{\infty} > \frac{|u_{2,i'} - u_{2,j}|}{2}\}} \\
&\le \sum_{j<i'} \mathbf{1}_{\{\Upsilon(n, M, \eta, p, \varepsilon, r) > \frac{|u_{2,i'} - u_{2,j}|}{2}\}} \quad \text{(since } \|\tilde{u}_2 - u_2\|_{\infty} \le \Upsilon(n, M, \eta, p, \varepsilon, r)\text{)} \\
&= \sum_{j<i'} \mathbf{1}_{\{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2 > |r_{i'} - r_j|\}} \\
&= \left[ \frac{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2}{\rho} \right] \le \frac{2\Upsilon(n, M, \eta, p, \varepsilon, r) \|r - \alpha e\|_2}{\rho}.
\end{align*}
$$

The same bound holds for $\sum_{j>i'} \mathbf{1}_{\{\tilde{u}_{2,j} > \tilde{u}_{2,i'}\}}$ and hence the statement of the Theorem follows.

5.4 Proof of Theorem 8

Recall that $\tilde{w} = \frac{\hat{\sigma}_1}{\eta p \sqrt{n}} \tilde{u}_2$. Let us denote

$$w = \frac{\sigma}{\eta p \sqrt{n}} \beta u_2 = \|r - \alpha e\|_2 \beta u_2 = (r - \alpha e) \beta.$$

For any norm $\|\cdot\|$ on $\mathbb{R}^n$, we have by triangle inequality that

$$\|\tilde{r} - w\| \le \|\tilde{w} - w\| + \left| \frac{e^T \tilde{w}}{n} \right| \|e\|.$$

Since $e^T w = 0$, therefore we obtain via Cauchy-Schwarz that

$$\left| \frac{e^T \tilde{w}}{n} \right| = \left| \frac{e^T (\tilde{w} - w)}{n} \right| \le \frac{\|\tilde{w} - w\|_2}{\sqrt{n}}.$$

Hence we have that

$$\|\tilde{r} - w\|_2 \le 2\|\tilde{w} - w\|_2 \quad \text{and} \quad \|\tilde{r} - w\|_{\infty} \le \|\tilde{w} - w\|_{\infty} + \frac{\|\tilde{w} - w\|_2}{\sqrt{n}} \le 2\|\tilde{w} - w\|_{\infty}.$$

It remains to bound $\|\tilde{w} - w\|_2$ and $\|\tilde{w} - w\|_{\infty}$. Our starting point will be

$$
\begin{align*}
\|\tilde{w} - w\| &= \frac{1}{\eta p \sqrt{n}} \|\hat{\sigma}_1 \tilde{u}_2 - \sigma_1 \beta u_2\| \\
&= \frac{1}{\eta p \sqrt{n}} \|(\hat{\sigma}_1 - \sigma_1) \tilde{u}_2 + \sigma_1 (\tilde{u}_2 - \beta u_2)\| \\
\Rightarrow \|\tilde{w} - w\| &\le \frac{1}{\eta p \sqrt{n}} (|\hat{\sigma}_1 - \sigma| \|\tilde{u}_2\| + \sigma_1 \|\tilde{u}_2 - \beta u_2\|) \\
&\le \frac{1}{\eta p \sqrt{n}} (\Delta \|\tilde{u}_2\| + \eta p \|r - \alpha e\|_2 \sqrt{n} \|\tilde{u}_2 - \beta u_2\|). \tag{5.10}
\end{align*}
$$
1. **Bounding** $\|\tilde{w} - w\|_2$. From (5.10), we have that

$$ \|\tilde{w} - w\|_2 \le \frac{1}{\eta p \sqrt{n}} (\Delta + \eta p \|r - \alpha e\|_2 \sqrt{n} \|\tilde{u}_2 - \beta u_2\|_2) $$

Plugging the expression for $\Delta$ from Lemma 12, and the bound on $\|\tilde{u}_2 - \beta u_2\|_2$ from Theorem 5, we readily obtain the bound in part (i) of the Theorem.

2. **Bounding** $\|\tilde{w} - w\|_{\infty}$. From (5.10), we have that

$$ \begin{aligned} \|\tilde{w} - w\|_{\infty} &\le \frac{1}{\eta p \sqrt{n}} [\Delta \|\tilde{u}_2\|_{\infty} + \sigma_1 \|\tilde{u}_2 - \beta u_2\|_{\infty}] \\ &\le \frac{1}{\eta p \sqrt{n}} [\Delta (\|\tilde{u}_2 - \beta u_2\|_{\infty} + \|u_2\|_{\infty}) + \sigma_1 \|\tilde{u}_2 - \beta u_2\|_{\infty}] \\ &= \frac{1}{\eta p \sqrt{n}} \left[ (\sigma_1 + \Delta) \|\tilde{u}_2 - \beta u_2\|_{\infty} + \Delta \frac{(M - \alpha)}{\|r - \alpha e\|_2} \right]. \end{aligned} \quad (5.11) $$

Plugging the expression for $\sigma_1$ (from Lemma 11) and $\Delta$ (from Lemma 12), followed by some simplification, we arrive at the bound in part (ii) of the Theorem.

## 5.5 Proof of Theorem 9

The outline is similar to the proof of Theorem 5 with some technical changes. We first note that the SVD of $(\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$ is given by $\sigma(u_1v_1^T + u_2v_2^T)$ where

$$ \sigma = \eta p \|(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2, \quad v_2 = u_1 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}e}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2}, \quad (5.12) $$

$$ v_1 = -u_2 = - \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)}{\|(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)\|_2}, $$

and $\alpha = \frac{r^T(\mathbb{E}[\bar{D}])^{-1}e}{e^T(\mathbb{E}[\bar{D}])^{-1}e}$. This is verified easily by proceeding as in the proof of Lemma 11. One can also readily see that $\sigma_{\min} \le \sigma \le \sigma_{\max}$. Now let us write

$$ \bar{D}^{-1/2} H \bar{D}^{-1/2} = (\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2} + \tilde{Z} $$

where $\tilde{Z} = \bar{D}^{-1/2}(H - \mathbb{E}[H])\bar{D}^{-1/2} + \bar{D}^{-1/2}\mathbb{E}[H]\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$. In order to bound $\|\tilde{Z}\|_2$, we will first need to establish the concentration of $\bar{D}$ around $\mathbb{E}[\bar{D}]$.

**Lemma 19** Denote $A(\eta, M) = \eta M^2 + (1 - \eta) \frac{M}{2}$. If $p \ge \frac{M^2}{9A(\eta, M)} \frac{\log n}{n}$ then,

$$ \mathbb{P}(\|\bar{D} - \mathbb{E}[\bar{D}]\|_2 \ge 2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n}) \le 2/n. $$

The proof is deferred to Appendix E. Conditioned on the event in Lemma 19, we have that

$$ \bar{D}_{ii} \in [\mathbb{E}[\bar{D}_{ii}] \pm 2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n}], \quad \forall i = 1, \dots, n. $$

In particular, if $2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n} \le \frac{1}{2} \min_i \mathbb{E}[\bar{D}_{ii}] = \frac{p\lambda_{\min}}{2}$, then

$$ \min_i \bar{D}_{ii} \ge p\lambda_{\min}/2. \quad (5.13) $$
1. **Bounding** $\|\tilde{w} - w\|_2$. From (5.10), we have that

$$ \|\tilde{w} - w\|_2 \le \frac{1}{\eta p \sqrt{n}} (\Delta + \eta p \|r - \alpha e\|_2 \sqrt{n} \|\tilde{u}_2 - \beta u_2\|_2) $$

Plugging the expression for $\Delta$ from Lemma 12, and the bound on $\|\tilde{u}_2 - \beta u_2\|_2$ from Theorem 5, we readily obtain the bound in part (i) of the Theorem.

2. **Bounding** $\|\tilde{w} - w\|_{\infty}$. From (5.10), we have that

$$ \begin{aligned} \|\tilde{w} - w\|_{\infty} &\le \frac{1}{\eta p \sqrt{n}} [\Delta \|\tilde{u}_2\|_{\infty} + \sigma_1 \|\tilde{u}_2 - \beta u_2\|_{\infty}] \\ &\le \frac{1}{\eta p \sqrt{n}} [\Delta (\|\tilde{u}_2 - \beta u_2\|_{\infty} + \|u_2\|_{\infty}) + \sigma_1 \|\tilde{u}_2 - \beta u_2\|_{\infty}] \\ &= \frac{1}{\eta p \sqrt{n}} \left[ (\sigma_1 + \Delta) \|\tilde{u}_2 - \beta u_2\|_{\infty} + \Delta \frac{(M - \alpha)}{\|r - \alpha e\|_2} \right]. \end{aligned} \quad (5.11) $$

Plugging the expression for $\sigma_1$ (from Lemma 11) and $\Delta$ (from Lemma 12), followed by some simplification, we arrive at the bound in part (ii) of the Theorem.

## 5.5 Proof of Theorem 9

The outline is similar to the proof of Theorem 5 with some technical changes. We first note that the SVD of $(\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$ is given by $\sigma(u_1v_1^T + u_2v_2^T)$ where

$$ \sigma = \eta p \|(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2, \quad v_2 = u_1 = \frac{(\mathbb{E}[\bar{D}])^{-1/2}e}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2}, \quad (5.12) $$

$$ v_1 = -u_2 = - \frac{(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)}{\|(\mathbb{E}[\bar{D}])^{-1/2}(r - \alpha e)\|_2}, $$

and $\alpha = \frac{r^T(\mathbb{E}[\bar{D}])^{-1}e}{e^T(\mathbb{E}[\bar{D}])^{-1}e}$. This is verified easily by proceeding as in the proof of Lemma 11. One can also readily see that $\sigma_{\min} \le \sigma \le \sigma_{\max}$. Now let us write

$$ \bar{D}^{-1/2} H \bar{D}^{-1/2} = (\mathbb{E}[\bar{D}])^{-1/2} \mathbb{E}[H] (\mathbb{E}[\bar{D}])^{-1/2} + \tilde{Z} $$

where $\tilde{Z} = \bar{D}^{-1/2}(H - \mathbb{E}[H])\bar{D}^{-1/2} + \bar{D}^{-1/2}\mathbb{E}[H]\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}$. In order to bound $\|\tilde{Z}\|_2$, we will first need to establish the concentration of $\bar{D}$ around $\mathbb{E}[\bar{D}]$.

**Lemma 19** Denote $A(\eta, M) = \eta M^2 + (1 - \eta) \frac{M}{2}$. If $p \ge \frac{M^2}{9A(\eta, M)} \frac{\log n}{n}$ then,

$$ \mathbb{P}(\|\bar{D} - \mathbb{E}[\bar{D}]\|_2 \ge 2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n}) \le 2/n. $$

The proof is deferred to Appendix E. Conditioned on the event in Lemma 19, we have that

$$ \bar{D}_{ii} \in [\mathbb{E}[\bar{D}_{ii}] \pm 2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n}], \quad \forall i = 1, \dots, n. $$

In particular, if $2(\sqrt{2} + 1)\sqrt{A(\eta, M)np \log n} \le \frac{1}{2} \min_i \mathbb{E}[\bar{D}_{ii}] = \frac{p\lambda_{\min}}{2}$, then

$$ \min_i \bar{D}_{ii} \ge p\lambda_{\min}/2. \quad (5.13) $$
Using triangle inequality, we obtain the bound

$$ \|\tilde{Z}\|_2 \le \underbrace{\|\bar{D}^{-1/2}(H - \mathbb{E}[H])\bar{D}^{-1/2}\|_2}_{\tilde{Z}_1} + \underbrace{\|\bar{D}^{-1/2}\mathbb{E}[H]\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}\|_2}_{\tilde{Z}_2}. $$

From sub-multiplicativity of the spectral norm and also (5.13), we have that

$$ \|\tilde{Z}_1\|_2 \le \frac{\|H - \mathbb{E}[H]\|_2}{\min_i \bar{D}_{ii}} \le \frac{2\Delta}{p\lambda_{\min}}. \quad (5.14) $$

Here, we used the bound $\|H - \mathbb{E}[H]\|_2 \le \Delta$ from Lemma 12. In order to bound $\|\tilde{Z}_2\|_2$, we add and subtract $\mathbb{E}[\bar{D}]$ from $\bar{D}$ and apply triangle inequality. This yields

$$ \begin{aligned} \|\tilde{Z}_2\|_2 &\le \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\mathbb{E}[H](\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 \\ &\quad + \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\mathbb{E}[H]\bar{D}^{-1/2}\|_2 + \|\bar{D}^{-1/2}\mathbb{E}[H](\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 \\ &\le \|\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\|_2^2 \|\mathbb{E}[H]\|_2 + 2\|\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\|_2 \|\mathbb{E}[H]\|_2 \|\bar{D}^{-1/2}\|_2. \end{aligned} \quad (5.15) $$

Note that

$$ \begin{aligned} \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 &= \|\bar{D}^{-1/2}(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \\ &\le \|\bar{D}^{-1/2}\|_2 \|(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}\|_2. \end{aligned} $$

Moreover, we can bound $\|(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})\|_2 \le \|\bar{D} - \mathbb{E}[\bar{D}]\|_2^{1/2}$ since $\bar{D}, \mathbb{E}[\bar{D}] > 0$ and $(\cdot)^{1/2}$ is operator monotone (see (Bhatia, 1996, Theorem X.1.1)). Using $\|\bar{D}^{-1/2}\|_2 \le \sqrt{\frac{2}{p\lambda_{\min}}}$, $\|(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \le \sqrt{\frac{1}{p\lambda_{\min}}}$ and the bound from Lemma 19 in (5.15), we get

$$ \|\tilde{Z}_2\|_2 \le C_1^2(\eta, M) \frac{(np \log n)^{1/4}}{p\lambda_{\min}} \|\mathbb{E}[H]\|_2 + \frac{2\sqrt{2}C_1(\eta, M)\|\mathbb{E}[H]\|_2(np \log n)^{1/4}}{(p\lambda_{\min})^{3/2}}, \quad (5.16) $$

where $C_1(\eta, M) = 4A^{1/4}(\eta, M)$. From (5.14), (5.16), we get

$$ \begin{aligned} \|\tilde{Z}\|_2 &\le \frac{2\Delta}{p\lambda_{\min}} + \frac{C_1(\eta, M)\|\mathbb{E}[H]\|_2(np \log n)^{1/4}}{(p\lambda_{\min})^{3/2}} \left[ \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right] \\ &\le 16M \sqrt{\frac{5}{3}}pn \frac{(2+\varepsilon)}{p\lambda_{\min}} + \frac{C_1(\eta, M)(np \log n)^{1/4}\sigma}{(p\lambda_{\min})^{3/2}} \left( \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right) \\ &= \tilde{\Delta} \end{aligned} $$

where the second inequality follows using Lemma 12 and (5.12). Let $\hat{U} = [\hat{u}_1 \hat{u}_2]$ denote the top two left singular vectors of $\bar{D}^{-1/2}H\bar{D}^{-1/2}$. Then from Wedin's bound, we know that

$$ \tilde{\Delta} < \sigma \Rightarrow \|(I - \hat{U}\hat{U}^T)U\|_2 \le \frac{\tilde{\Delta}}{\sigma - \tilde{\Delta}} = \delta. \quad (5.17) $$

Finally, Lemma 14 can be used here unchanged. Hence if $\delta \le 1/2 \Leftrightarrow \tilde{\Delta} \le \sigma/3$, then there exists $\beta \in \{-1, 1\}$ such that $\|\tilde{u}_2 - \beta u_2\|_2^2 \le 10\delta$. Since $\tilde{\Delta} \le \sigma_{\min}/3 \Leftrightarrow \delta \le \frac{3\tilde{\Delta}}{2\sigma_{\min}} \Rightarrow \tilde{\Delta} \le \sigma/3$, we obtain the stated bound on $\|\tilde{u}_2 - \beta u_2\|_2^2$. The lower bound on the success probability follows by applying the union bound to the events in Lemmas 12, 19.
Using triangle inequality, we obtain the bound

$$ \|\tilde{Z}\|_2 \le \underbrace{\|\bar{D}^{-1/2}(H - \mathbb{E}[H])\bar{D}^{-1/2}\|_2}_{\tilde{Z}_1} + \underbrace{\|\bar{D}^{-1/2}\mathbb{E}[H]\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\mathbb{E}[H](\mathbb{E}[\bar{D}])^{-1/2}\|_2}_{\tilde{Z}_2}. $$

From sub-multiplicativity of the spectral norm and also (5.13), we have that

$$ \|\tilde{Z}_1\|_2 \le \frac{\|H - \mathbb{E}[H]\|_2}{\min_i \bar{D}_{ii}} \le \frac{2\Delta}{p\lambda_{\min}}. \quad (5.14) $$

Here, we used the bound $\|H - \mathbb{E}[H]\|_2 \le \Delta$ from Lemma 12. In order to bound $\|\tilde{Z}_2\|_2$, we add and subtract $\mathbb{E}[\bar{D}]$ from $\bar{D}$ and apply triangle inequality. This yields

$$ \begin{aligned} \|\tilde{Z}_2\|_2 &\le \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\mathbb{E}[H](\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 \\ &\quad + \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\mathbb{E}[H]\bar{D}^{-1/2}\|_2 + \|\bar{D}^{-1/2}\mathbb{E}[H](\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 \\ &\le \|\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\|_2^2 \|\mathbb{E}[H]\|_2 + 2\|\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2}\|_2 \|\mathbb{E}[H]\|_2 \|\bar{D}^{-1/2}\|_2. \end{aligned} \quad (5.15) $$

Note that

$$ \begin{aligned} \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})\|_2 &= \|\bar{D}^{-1/2}(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \\ &\le \|\bar{D}^{-1/2}\|_2 \|(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}\|_2. \end{aligned} $$

Moreover, we can bound $\|(\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2})\|_2 \le \|\bar{D} - \mathbb{E}[\bar{D}]\|_2^{1/2}$ since $\bar{D}, \mathbb{E}[\bar{D}] > 0$ and $(\cdot)^{1/2}$ is operator monotone (see (Bhatia, 1996, Theorem X.1.1)). Using $\|\bar{D}^{-1/2}\|_2 \le \sqrt{\frac{2}{p\lambda_{\min}}}$, $\|(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \le \sqrt{\frac{1}{p\lambda_{\min}}}$ and the bound from Lemma 19 in (5.15), we get

$$ \|\tilde{Z}_2\|_2 \le C_1^2(\eta, M) \frac{(np \log n)^{1/4}}{p\lambda_{\min}} \|\mathbb{E}[H]\|_2 + \frac{2\sqrt{2}C_1(\eta, M)\|\mathbb{E}[H]\|_2(np \log n)^{1/4}}{(p\lambda_{\min})^{3/2}}, \quad (5.16) $$

where $C_1(\eta, M) = 4A^{1/4}(\eta, M)$. From (5.14), (5.16), we get

$$ \begin{aligned} \|\tilde{Z}\|_2 &\le \frac{2\Delta}{p\lambda_{\min}} + \frac{C_1(\eta, M)\|\mathbb{E}[H]\|_2(np \log n)^{1/4}}{(p\lambda_{\min})^{3/2}} \left[ \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right] \\ &\le 16M \sqrt{\frac{5}{3}}pn \frac{(2+\varepsilon)}{p\lambda_{\min}} + \frac{C_1(\eta, M)(np \log n)^{1/4}\sigma}{(p\lambda_{\min})^{3/2}} \left( \frac{C_1(\eta, M)(np \log n)^{1/4}}{\sqrt{p\lambda_{\min}}} + 2\sqrt{2} \right) \\ &= \tilde{\Delta} \end{aligned} $$

where the second inequality follows using Lemma 12 and (5.12). Let $\hat{U} = [\hat{u}_1 \hat{u}_2]$ denote the top two left singular vectors of $\bar{D}^{-1/2}H\bar{D}^{-1/2}$. Then from Wedin's bound, we know that

$$ \tilde{\Delta} < \sigma \Rightarrow \|(I - \hat{U}\hat{U}^T)U\|_2 \le \frac{\tilde{\Delta}}{\sigma - \tilde{\Delta}} = \delta. \quad (5.17) $$

Finally, Lemma 14 can be used here unchanged. Hence if $\delta \le 1/2 \Leftrightarrow \tilde{\Delta} \le \sigma/3$, then there exists $\beta \in \{-1, 1\}$ such that $\|\tilde{u}_2 - \beta u_2\|_2^2 \le 10\delta$. Since $\tilde{\Delta} \le \sigma_{\min}/3 \Leftrightarrow \delta \le \frac{3\tilde{\Delta}}{2\sigma_{\min}} \Rightarrow \tilde{\Delta} \le \sigma/3$, we obtain the stated bound on $\|\tilde{u}_2 - \beta u_2\|_2^2$. The lower bound on the success probability follows by applying the union bound to the events in Lemmas 12, 19.
## 5.6 Proof of Theorem 10

Recall that $\tilde{w} = \frac{\hat{\sigma}_1}{\eta p \|\bar{D}^{-1/2}e\|_2} \bar{D}^{1/2} \tilde{u}_2$ and $\tilde{r} = \tilde{w} - \frac{e^T \tilde{w}}{n} e$. Let us denote

$$w = \frac{\sigma}{\eta p \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \beta u_2 = (r - \alpha e)\beta, \\ w' = w - \frac{e^T w}{n} e = \beta(r - (\alpha + \alpha')e),$$

where $\alpha' = \frac{e^T(r - \alpha e)}{n}$. Then using triangle inequality, we can bound $\|\tilde{r} - w'\|_2$ as

$$\|\tilde{r} - w'\|_2 \le \|w - \tilde{w}\|_2 + \left| \frac{e^T(\tilde{w} - w)}{n} \right| \|e\|_2 \le 2\|w - \tilde{w}\|_2.$$

In order to bound $\|w - \tilde{w}\|_2$, note that we can write $\tilde{w} - w = \frac{1}{\eta p}(a_1 + a_2 + a_3 + a_4)$ where

$$\begin{align*}
a_1 &= \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} \bar{D}^{1/2} \tilde{u}_2 - \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_2 &= \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\hat{\sigma}_1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_3 &= \frac{\hat{\sigma}_1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_4 &= \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \beta u_2.
\end{align*}$$

Since $\|\tilde{w} - w\|_2 \le \frac{1}{\eta p} \sum_{i=1}^4 \|a_i\|_2$ we will now bound $\|a_i\|_2$ for each $i$. Before proceeding, recall from the proof of Theorem 5.5 that $\min_i \bar{D}_{ii} \ge \frac{p\lambda_{\min}}{2}$. We also have

$$\max_i \bar{D}_{ii} \le p\lambda_{\max} + \frac{p\lambda_{\min}}{2} \le \frac{3p\lambda_{\max}}{2},$$

and $\hat{\sigma}_1 \le \sigma + \tilde{\Delta}$ where the latter is due to Weyl's inequality.

**Bounding $\|a_1\|_2$.** Since $\sigma \le \sigma_{\max}$, we have that

$$\|a_1\|_2 \le \frac{\tilde{\Delta} + \sigma_{\max}}{\|\bar{D}^{-1/2}e\|_2} \|\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2}\|_2 \le \frac{\tilde{\Delta} + \sigma_{\max}}{\|\bar{D}^{-1/2}e\|_2} \|\bar{D} - \mathbb{E}[\bar{D}]\|_2^{1/2} \quad (5.18)$$

where the second inequality is due to $\bar{D}, \mathbb{E}[\bar{D}] > 0$ and since $(\cdot)^{1/2}$ is operator monotone (see (Bhatia, 1996, Theorem X.1.1)). Using the bound $\|\bar{D}^{-1/2}e\|_2 \ge \sqrt{\frac{2n}{3p\lambda_{\max}}}$, along with the bound for $\|\bar{D} - \mathbb{E}[\bar{D}]\|_2$ (from Lemma 19) in (5.18), we arrive at

$$\|a_1\|_2 \le \sqrt{\frac{3}{2n} p \lambda_{\max} (\sigma_{\max} + \tilde{\Delta}) (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}}.$$
## 5.6 Proof of Theorem 10

Recall that $\tilde{w} = \frac{\hat{\sigma}_1}{\eta p \|\bar{D}^{-1/2}e\|_2} \bar{D}^{1/2} \tilde{u}_2$ and $\tilde{r} = \tilde{w} - \frac{e^T \tilde{w}}{n} e$. Let us denote

$$w = \frac{\sigma}{\eta p \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \beta u_2 = (r - \alpha e)\beta, \\ w' = w - \frac{e^T w}{n} e = \beta(r - (\alpha + \alpha')e),$$

where $\alpha' = \frac{e^T(r-\alpha e)}{n}$. Then using triangle inequality, we can bound $\|\tilde{r} - w'\|_2$ as

$$\|\tilde{r} - w'\|_2 \le \|w - \tilde{w}\|_2 + \left| \frac{e^T(\tilde{w} - w)}{n} \right| \|e\|_2 \le 2\|w - \tilde{w}\|_2.$$

In order to bound $\|w - \tilde{w}\|_2$, note that we can write $\tilde{w} - w = \frac{1}{\eta p}(a_1 + a_2 + a_3 + a_4)$ where

$$\begin{align*}
a_1 &= \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} \bar{D}^{1/2} \tilde{u}_2 - \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_2 &= \frac{\hat{\sigma}_1}{\|\bar{D}^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\hat{\sigma}_1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_3 &= \frac{\hat{\sigma}_1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2, \\
a_4 &= \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \tilde{u}_2 - \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} (\mathbb{E}[\bar{D}])^{1/2} \beta u_2.
\end{align*}$$

Since $\|\tilde{w} - w\|_2 \le \frac{1}{\eta p} \sum_{i=1}^4 \|a_i\|_2$ we will now bound $\|a_i\|_2$ for each $i$. Before proceeding, recall from the proof of Theorem 5.5 that $\min_i \bar{D}_{ii} \ge \frac{p\lambda_{\min}}{2}$. We also have

$$\max_i \bar{D}_{ii} \le p\lambda_{\max} + \frac{p\lambda_{\min}}{2} \le \frac{3p\lambda_{\max}}{2},$$

and $\hat{\sigma}_1 \le \sigma + \tilde{\Delta}$ where the latter is due to Weyl's inequality.

**Bounding $\|a_1\|_2$.** Since $\sigma \le \sigma_{\max}$, we have that

$$\|a_1\|_2 \le \frac{\tilde{\Delta} + \sigma_{\max}}{\|\bar{D}^{-1/2}e\|_2} \|\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2}\|_2 \le \frac{\tilde{\Delta} + \sigma_{\max}}{\|\bar{D}^{-1/2}e\|_2} \|\bar{D} - \mathbb{E}[\bar{D}]\|_2^{1/2} \quad (5.18)$$

where the second inequality is due to $\bar{D}, \mathbb{E}[\bar{D}] > 0$ and since $(\cdot)^{1/2}$ is operator monotone (see (Bhatia, 1996, Theorem X.1.1)). Using the bound $\|\bar{D}^{-1/2}e\|_2 \ge \sqrt{\frac{2n}{3p\lambda_{\max}}}$, along with the bound for $\|\bar{D} - \mathbb{E}[\bar{D}]\|_2$ (from Lemma 19) in (5.18), we arrive at

$$\|a_1\|_2 \le \sqrt{\frac{3}{2n} p \lambda_{\max} (\sigma_{\max} + \tilde{\Delta}) (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}}.$$
**Bounding $\|a_2\|_2$.** We begin by noting that

$$ \|a_2\|_2 = \widehat{\sigma_1} \left| \frac{1}{\|\bar{D}^{-1/2}e\|_2} - \frac{1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \right| \|(\mathbb{E}[\bar{D}])^{1/2}\tilde{u}_2\|_2 \le (\sigma_{\max} + \tilde{\Delta}) \left| \frac{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 - \|\bar{D}^{-1/2}e\|_2}{\|\bar{D}^{-1/2}e\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \right| $$

Since $\|\bar{D}^{-1/2}e\|_2 \ge \sqrt{\frac{2n}{3p\lambda_{\max}}}$, $\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 \ge \frac{\sqrt{n}}{\sqrt{p\lambda_{\max}}}$, and

$$ \begin{aligned} \left| \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 - \|\bar{D}^{-1/2}e\|_2 \right| & \le \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})e\|_2 \\ & \le \|\bar{D}^{-1/2}\|_2 \|\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2}\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \sqrt{n} \\ & \le \frac{\sqrt{2}}{p\lambda_{\min}} (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}, \end{aligned} $$

we can bound $\|a_2\|_2$ as

$$ \|a_2\|_2 \le \sqrt{\frac{3}{n}} (\sigma_{\max} + \tilde{\Delta}) \frac{\lambda_{\max}}{\lambda_{\min}} (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}. $$

**Bounding $\|a_3\|_2$.** This is easily achieved by noting that

$$ \|a_3\|_2 = \frac{|\widehat{\sigma}_1 - \sigma|}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \|(\mathbb{E}[\bar{D}])^{1/2}\tilde{u}_2\|_2 \le \frac{\tilde{\Delta}\sqrt{p\lambda_{\max}}}{\sqrt{n}/\sqrt{p\lambda_{\max}}} = \frac{\tilde{\Delta}p\lambda_{\max}}{\sqrt{n}}. $$

**Bounding $\|a_4\|_2$.** This is also easily achieved by noting that

$$ \begin{aligned} \|a_4\|_2 & \le \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \|(\mathbb{E}[\bar{D}])^{1/2}\|_2 \|\tilde{u}_2 - \beta u_2\|_2 \\ & \le \frac{\sigma_{\max}\sqrt{p\lambda_{\max}}}{\sqrt{n}/\sqrt{p\lambda_{\max}}} \left( \sqrt{15}\sqrt{\frac{\tilde{\Delta}}{\sigma_{\min}}} \right) \\ & \le \sqrt{\frac{15\tilde{\Delta}}{\sigma_{\min}n}} \sigma_{\max} p \lambda_{\max}. \end{aligned} $$

The stated bound now follows from $\|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \le \frac{2}{\eta p} \sum_{i=1}^4 \|a_i\|_2$.

## 6. Matrix completion as a preprocessing step for ranking

Low-rank matrix completion is the problem of recovering the missing entries of a low-rank matrix given a subset of its entries. This line of research started with the results in Candès and Recht (2009) and Candès and Tao (2010) which showed that given a rank-r matrix C of size $n_1 \times n_2$ (with $n = \max\{n_1, n_2\}$), one can recover it by observing only $O(nr)$ polylog(n)) randomly selected entries (under some assumptions on C) via a simple convex optimization algorithm. This was partly inspired by similar approaches used previously in the compressed sensing literature (Candes et al., 2006; Candès and Tao, 2006). This problem has received tremendous attention in the last decade, due to a number of applications such as in phase
**Bounding $\|a_2\|_2$.** We begin by noting that

$$ \|a_2\|_2 = \widehat{\sigma_1} \left| \frac{1}{\|\bar{D}^{-1/2}e\|_2} - \frac{1}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \right| \|(\mathbb{E}[\bar{D}])^{1/2}\tilde{u}_2\|_2 \le (\sigma_{\max} + \tilde{\Delta}) \left| \frac{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 - \|\bar{D}^{-1/2}e\|_2}{\|\bar{D}^{-1/2}e\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \right| $$

Since $\|\bar{D}^{-1/2}e\|_2 \ge \sqrt{\frac{2n}{3p\lambda_{\max}}}$, $\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 \ge \frac{\sqrt{n}}{\sqrt{p\lambda_{\max}}}$, and

$$ \begin{aligned} \left| \|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2 - \|\bar{D}^{-1/2}e\|_2 \right| &\le \|(\bar{D}^{-1/2} - (\mathbb{E}[\bar{D}])^{-1/2})e\|_2 \\ &\le \|\bar{D}^{-1/2}\|_2 \|\bar{D}^{1/2} - (\mathbb{E}[\bar{D}])^{1/2}\|_2 \|(\mathbb{E}[\bar{D}])^{-1/2}\|_2 \sqrt{n} \\ &\le \frac{\sqrt{2}}{p\lambda_{\min}} (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}, \end{aligned} $$

we can bound $\|a_2\|_2$ as

$$ \|a_2\|_2 \le \sqrt{\frac{3}{n}} (\sigma_{\max} + \tilde{\Delta}) \frac{\lambda_{\max}}{\lambda_{\min}} (2\sqrt{2} + 1)^{1/2} (A(\eta, M)np \log n)^{1/4}. $$

**Bounding $\|a_3\|_2$.** This is easily achieved by noting that

$$ \|a_3\|_2 = \frac{|\widehat{\sigma}_1 - \sigma|}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \|(\mathbb{E}[\bar{D}])^{1/2}\tilde{u}_2\|_2 \le \frac{\tilde{\Delta}\sqrt{p\lambda_{\max}}}{\sqrt{n}/\sqrt{p\lambda_{\max}}} = \frac{\tilde{\Delta}p\lambda_{\max}}{\sqrt{n}}. $$

**Bounding $\|a_4\|_2$.** This is also easily achieved by noting that

$$ \begin{aligned} \|a_4\|_2 &\le \frac{\sigma}{\|(\mathbb{E}[\bar{D}])^{-1/2}e\|_2} \|(\mathbb{E}[\bar{D}])^{1/2}\|_2 \|\tilde{u}_2 - \beta u_2\|_2 \\ &\le \frac{\sigma_{\max}\sqrt{p\lambda_{\max}}}{\sqrt{n}/\sqrt{p\lambda_{\max}}} \left( \sqrt{15}\sqrt{\frac{\tilde{\Delta}}{\sigma_{\min}}} \right) \\ &\le \sqrt{\frac{15\tilde{\Delta}}{\sigma_{\min}n}} \sigma_{\max} p \lambda_{\max}. \end{aligned} $$

The stated bound now follows from $\|\tilde{r} - \beta(r - (\alpha + \alpha')e)\|_2 \le \frac{2}{\eta p} \sum_{i=1}^4 \|a_i\|_2$.

# 6. Matrix completion as a preprocessing step for ranking

Low-rank matrix completion is the problem of recovering the missing entries of a low-rank matrix given a subset of its entries. This line of research started with the results in Candès and Recht (2009) and Candès and Tao (2010) which showed that given a rank-r matrix C of size $n_1 \times n_2$ (with $n = \max\{n_1, n_2\}$), one can recover it by observing only $O(nr)$ polylog(n)) randomly selected entries (under some assumptions on C) via a simple convex optimization algorithm. This was partly inspired by similar approaches used previously in the compressed sensing literature (Candes et al., 2006; Candès and Tao, 2006). This problem has received tremendous attention in the last decade, due to a number of applications such as in phase
retrieval (Candès et al., 2013), computer vision (Tomasi and Kanade, 1992; Özyeşil et al., 2017) and sensor network localization (Cucuringu et al., 2012a) to name a few.

More formally, let us assume for simplicity that $r = O(1)$ and the SVD of $C = \sum_{i \in [r]} \sigma_i u_i v_i^T$ satisfies

$$ \|u_i\|_{\infty} \le \sqrt{\frac{\mu}{n_1}}, \quad \|v_i\|_{\infty} \le \sqrt{\frac{\mu}{n_2}}; \quad \forall i, j, \qquad (6.1) $$

for some $\mu \ge 1$. Here, $\mu$ is a coherence parameter which measures how spread out the entries of $u_i, v_i$ are - the smaller the value of $\mu$, the better. Say we observe $m$ entries of $C$ on a subset $\Omega \subset [n_1] \times [n_2]$, sampled uniformly at random. Denoting $P_{\Omega}: \mathbb{R}^{n_1 \times n_2} \to \mathbb{R}^{n_1 \times n_2}$ to be the projection operator on $\Omega$, it was shown in Candès and Tao (2010) that the solution $\hat{C}$ of

$$ \min \|X\|_{*} \quad \text{s.t} \quad P_{\Omega}(X) = P_{\Omega}(C) \qquad (6.2) $$

equals $C$ with high probability, provided $m = \Omega(\mu^4 n \log^2 n)$. The nuclear norm minimization problem (6.2) is a SDP and hence can be solved in polynomial time using, for eg., interior point methods. In fact, the result holds under a different Bernoulli sampling model too, wherein each entry of $M$ is observed independently with a certain probability (Candès and Recht, 2009, Section 4.1). Moreover, say that the observations are noisy, i.e, we observe

$$ Y_{i,j} = C_{i,j} + Z_{i,j}; \quad (i,j) \in \Omega \Leftrightarrow P_{\Omega}(Y) = P_{\Omega}(C) + P_{\Omega}(Z), $$

where $Z$ is the noise matrix. Say $\|P_{\Omega}(Z)\|_F \le \delta$. Then, it was shown in Candès and Plan (2010) (under certain additional conditions) that the solution of

$$ \min \|X\|_{*} \quad \text{s.t} \quad \|P_{\Omega}(X) - P_{\Omega}(Y)\|_{F} \le \delta, \qquad (6.3) $$

is stable, i.e., the estimation error $\|\hat{C} - C\|_F$ is bounded by a term proportional to $\delta$.

In our setting, we observe (noisy versions of) a subset of the entries of the $n \times n$ matrix $C = re^T - er^T$ where each off-diagonal entry of $C$ is revealed with probability $p$. Since $C$ has rank 2, it is natural to consider estimating $C$ via matrix completion as a preprocessing step, and then subsequently applying SVD-RS or SVD-NRS on the obtained estimate of $C$ for recovering the underlying ranks and scores. In order to understand the sample complexity for successful matrix completion, we need to express $\mu$ defined in (6.1) in terms of the score vector $r$. To this end, we see from Lemma 11 that

$$ \|u_2\|_{\infty}, \|v_1\|_{\infty} = \frac{\|r - \alpha e\|_{\infty}}{\|r - \alpha e\|_2} \le \frac{M - \alpha}{\|r - \alpha e\|_2}, $$

and $\|u_1\|_{\infty}, \|v_2\|_{\infty} = 1/\sqrt{n}$. Hence it follows that

$$ \mu = \max \left\{ \frac{(M - \alpha)\sqrt{n}}{\|r - \alpha e\|_2}, 1 \right\}. $$

In order to get a sense of the scaling of $\mu$, consider the setting $r_i = i$. As seen before, we then have $\|r - \alpha e\|_2 = \Theta(n^{3/2})$, $M = n$, and $\alpha = \Theta(n)$. This in turn implies that $\mu = \Theta(1)$ which is the ideal scenario for matrix completion. In our numerical experiments,