# Optimal Permutation Estimation in Crowd-Sourcing problems

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen

► To cite this version:

Emmanuel Pilliat, Alexandra Carpentier, Nicolas Verzelen. Optimal Permutation Estimation in Crowd-Sourcing problems. 2023. hal-03837935v3

HAL Id: hal-03837935

https://hal.science/hal-03837935v3

Preprint submitted on 29 Mar 2023

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.
# Optimal Permutation Estimation in Crowd-Sourcing problems

Emmanuel Pilliat, Alexandra Carpentier, and Nicolas Verzelen

**Abstract:** Motivated by crowd-sourcing applications, we consider a model where we have partial observations from a bivariate isotonic $n \times d$ matrix with an unknown permutation $\pi^*$ acting on its rows. Focusing on the twin problems of recovering the permutation $\pi^*$ and estimating the unknown matrix, we introduce a polynomial-time procedure achieving the minimax risk for these two problems, this for all possible values of $n$, $d$, and all possible sampling efforts. Along the way, we establish that, in some regimes, recovering the unknown permutation $\pi^*$ is considerably simpler than estimating the matrix.

## 1. Introduction

We consider a crowd-sourcing problem with $n$ experts and $d$ questions. For an unknown matrix $M$, $M_{i,j} \in [0,1]$ stands for the ability of expert $i$ at question $j$. For the purpose of calibrating the model, we receive partial and noisy observations of the matrix $M$ and our goal is to rank the experts according to their ability. Earlier models in crowd-labelling problems or in the related problems of pairwise comparisons typically assumed that the matrix $M$ belongs to a parametric model [4, 16, 28, 11, 5], a prominent example being Bradley-Luce-Terry model. While there has been significant progress in this direction, such models do not tend to fit well real-world data [19, 2].

To address this issue, there has been a recent interest in the class of permutation-based models [6, 24, 25, 17, 8, 12, 21, 27] where it is only assumed that the matrix $M$ satisfies some shape-constrained conditions before one (or two) permutations acts on the rows (and possibly on the columns) of $M$. Quite surprisingly, it has been established in [25] that, at least in some settings, the matrix $M$ can be estimated at the same rate in those non-parametric models as in classical parametric models by relying on the least-square estimator on the class of permuted bi-isotonic matrices. Unfortunately, the corresponding class of matrices is highly non-convex and no polynomial-time algorithm is known for computing this least-square estimator. Furthermore, known computationally efficient procedures such as spectral estimators [6, 7] only achieve sub-optimal convergence rates. This has led several authors to conjecture the existence of computational-statistical trade-offs [12, 26]. Despite recent progress in this direction [17, 14], the fundamental limits of polynomial-time algorithms for this class of problems remain largely unknown.

Arguably, for most applications, the primary objective is to recover the underlying permutation $\pi^*$ acting on the rows or equivalently to rank the experts accordingly. While estimation of the full matrix $M$ is closely related to ranking, it is also of a quite different nature as argued below. In this work, we investigate the estimation of the permutation $\pi^*$ by characterizing the minimax risk for estimating $\pi^*$ in a permuted shape-constrained model, introducing a polynomial-time procedure nearly achieving this risk bound. As a byproduct, we also disprove the existence of a computational-statistical gap for the reconstruction of the matrix $M$.

### 1.1. Problem formulation

A bounded matrix $B \in [0,1]^{n \times d}$ is said to be bi-isotonic if it satisfies $B_{i,j} \le B_{i+1,j}$ and $B_{i,j} \le B_{i,j+1}$ for any $i \in [n-1]$ and $j \in [d-1]$. Henceforth, we write $C_{\text{BISO}}$ for the collection of such $n \times d$ bounded bi-isotonic matrices.
In this work, we assume that the matrix M is a row-permuted bi-isotonic matrix as in [17, 14]. In other words, up to a single permutation $\pi^*$ of $[n]$, the matrix $M_{\pi^{-1}}$ defined by $(M_{\pi^{-1}})_{i,j} = (M_{\pi^{-1}(i),j})$ is bi-isotonic. From a modeling viewpoint, this amounts to assuming that the $d$ questions are ordered from the most difficult question to the most simple question. The permutation $\pi^*$ is not necessarily unique, but the corresponding permuted matrix $M_{\pi^{-1}}$ is unique. Despite that, we refer, with a slight abuse of terminology, to $\pi^*$ as the oracle permutation. With this definition, $\pi^{-1}(i)$ corresponds to any $i$-th smallest row (or equivalently expert to use the crowd-sourcing terminology) in the matrix M. In the following, the $i$th row of M is referred to as expert $i$, whereas the $k$th column is referred to as question $k$.

We consider an observation-scheme where the statistician has partial access to noisy observations Y of M such that

$$Y = M + E, \qquad (1)$$

where the entries of E are centered, independent, subGaussian - see definition 2.2 of [32] - with Orlicz norm at most $\zeta$, but are not necessarily identically distributed. In particular, this model encompasses binary observations $Y_{i,k} \sim Ber(M_{i,k})$ which arise in crowd-labelling problems, in which case we have $\zeta = 1$. In the following, we refer to $\zeta$ as the noise level.

As usual in the literature -e.g. [17], we use the Poissonization trick to model the partial observations. Given some $\lambda > 0$, which is henceforth referred as the sampling effort, we have $N = Poi(\lambda nd)$ observations of the form

$$(x_t, y_t), \quad t=, 1, \dots, N, \qquad (2)$$

where the position $x_t$ is sampled uniformly in $[n] \times [d]$, and $y_t = M_{x_t} + E_{x_t}$ is an independent observation of matrix Y of (1) at position $x_t$. Conditionally to N, this scheme is equivalent to a uniform sampling scheme with replacement [18]. If $\lambda < 1$, then a specific entry of Y is sampled at least once with probability $1 - e^{-\lambda}$ which is close to $\lambda$. More generally, $\lambda$ corresponds to the expected number of times a specific entry of Y is observed, so that $\lambda > 1$ would correspond to the situation where entries are sampled multiple times.

Since our aim is to recover the permutation $\pi^*$ from the partial observations $(x_t, y_t)$, we consider, for some estimator $\hat{\pi}$, the following error metric

$$l(\hat{\pi}; \pi^*) = \|M_{\hat{\pi}^{-1}} - M_{\pi^{-1}}\|_F^2, \qquad (3)$$

where $\|\cdot\|_F$ stands for the Frobenius norm. This loss quantifies the distance between the matrix M ordered according to the oracle permutation $\pi^*$ and the matrix M ordered according to the estimated permutation. When $\pi^*$ is not unique, the error $l(\pi^*, \pi^*)$ between any two oracle permutations is zero. If $\hat{\pi}$ and $\pi^*$ only differ by a transposition or equivalently if the ranking $\hat{\pi}$ and $\pi^*$ only differ on two experts, then $l(\hat{\pi}; \pi^*)$ is twice the square Euclidean distance between the corresponding rows of M. More generally, $l(\hat{\pi}; \pi^*)$ interprets as the sum over all $i = 1, \dots, n$ of the square Euclidean distance between the $i$-th smallest row of M according to $\hat{\pi}$ and according to the oracle ranking $\pi^*$.

The loss (3) is ubiquitous when one aims at estimating the matrix M in Frobenius norm, that is building an estimator $\widehat{M}$ such that $\|\widehat{M} - M\|_F^2$ is as small as possible -see e.g. [25, 17, 14]. Indeed, estimating $\pi^*$ by $\hat{\pi}$ is a first step towards building an estimator of M by doing as if $M_{\hat{\pi}^{-1}}$ was bi-isotonic. It turns out that the error in $\|\widehat{M} - M\|_F^2$ decomposes as the sum of two terms, one of them being $l(\hat{\pi}, \pi^*)$ while the other one does not really depend on $\hat{\pi}$. Conversely, an estimator $\widehat{M}$ can be easily transformed into an estimator $\hat{\pi}$ whose loss $l(\hat{\pi}, \pi^*)$ is controlled by $\|\widehat{M} - M\|_F^2$. See [25, 17] for further discussions. In summary, controlling $l(\hat{\pi}; \pi^*)$ is important in order to evaluate to what extent $\pi^*$ is well estimated, but it is also the key stepping stone towards a good estimation of the matrix M.
In some works, the authors directly consider distances on the symmetric group of permutations. Examples include the Kendall tau distance $d_{KT}(\pi, \pi') = \sum_{(i,j):\pi(i)<\pi(j)} \mathbf{1}\{\pi'(i) > \pi'(j)\}$ or the $l_\infty$ distance $\|\pi - \pi'\|_\infty = \max_{i \in [n]} |\pi(i) - \pi'(i)|$ -see [5, 18] in the noisy sorting model. However, those distances are not well suited to handle the non-parametric class of bi-isotonic matrices, because to control them we would need to make assumptions on the separation between the rows of the matrix $M$ -see Appendix A of [25].

Equipped with this notation, we consider the minimax risk of permutation recovery as a function of the number $n$ of experts, the number $d$ of question, the sampling effort $\lambda$, and the noise level $\zeta$.

$$ \mathcal{R}^*[n, d, \lambda, \zeta] = \inf_{\hat{\pi}} \sup_{\pi^* \in \Pi_n} \sup_{M: M_{\pi^{*-1}} \in \mathbb{C}_{\text{BISO}}} \mathbb{E} \left[ \| M_{\hat{\pi}^{-1}} - M_{\pi^{*-1}} \|_F^2 \right], \quad (4) $$

where $\Pi_n$ stands for the collection of all permutations of $[n]$. In particular, our general aim is to tightly control this minimax risk and, if possible, to provide a computationally efficient procedure achieving this minimax risk. Although our primary interest lies in the estimation of $\pi^*$, we also consider the minimax estimation risk of $M$

$$ \mathcal{R}_{est}^*[n, d, \lambda, \zeta] = \inf_{\hat{M}} \sup_{\pi^* \in \Pi_n} \sup_{M: M_{\pi^{*-1}} \in \mathbb{C}_{\text{BISO}}} \mathbb{E} \left[ \| \hat{M} - M \|_F^2 \right]. \quad (5) $$

as studied in [25, 14, 17, 20] in order to assess the performances of our computationally efficient procedures.

## 1.2. Related work and open questions

The most relevant body of work to the current paper is that on estimating square matrices $M$ satisfying the so-called strong stochastic transitivity class (SST) [6, 25]. A matrix $M$ belongs to the SST class if (i) $M$ is skew-symmetric that is $M + M^T = ee^T$ where $e$ is the constant vector of size $n$ and (ii) there exists a common permutation $\pi^*$ of $[n]$ such that row and column-permuted matrix $M_{\pi^{*-1}\pi^{*-1}}$ is bi-isotonic. This class is suited for considering pairwise comparisons problems. Shah et al. [25] consider the full observation setting, namely a setting where each entry of the matrix $M$ is observed once in noise - which is in some sense akin to $\lambda = 1$ in our Poissonian scheme<fnref n="1" />. They proved that the minimax risk for estimating $M$ in square Frobenius distance is, up to logarithmic terms, of the order of $n$ and is achieved by the corresponding least-square estimator over the SST class. Unfortunately, this estimator cannot be efficiently computed. They also analyzed an efficient spectral estimator achieving the rate $n^{3/2}$. This rate is also achieved [25] by the near-linear time Borda count algorithm CRL that simply ranks the individuals according to the row sums of the observations and then plugs the corresponding permutation to estimate $M$. See also [8] for related results. This led some authors [12, 26] to conjecture the existence of a $\sqrt{n}$ computational gap for SST matrices and for other shape-constrained matrices with unknown permutation.

In crowd-sourcing problems where $M \in [0, 1]^{n \times d}$, non-parametric models [17] assume that the matrix $M$ is bi-isotonic up to a permutation $\pi^*$ of the rows (experts) - and sometimes also up to a permutation $\tau^*$ of the columns (questions)<fnref n="2" />. In this paper as in this literature review, we focus however solely on the case where $M$ is bi-isotonic up to a permutation $\pi^*$ of the rows (experts). Mao et al. [17] have established the minimax risk $\mathcal{R}_{est}^*[n, d, \lambda, 1]$ for estimating $M$ in the specific case where $n \ge d$. In the arguably most interesting regime of partial observations $\lambda \le 1$, they prove that this minimax risk is of the order of $n/\lambda \wedge (nd)$. This rate is achieved by

<footnote n="1">¹In the Poissonian scheme, each entry is observed at least once with probability $1 - e^{-\lambda}$.</footnote>

<footnote n="2">²This would correspond to the situation where the corresponding ordering of the questions is also unknown.</footnote>
the inefficient least-square estimator. Furthermore, Mao et al. [17] were the first to narrow the conjectured computational gap by introducing a new efficient procedure called one-dimensional sorting. In the square case $n = d$ with full observations, these procedures achieve (up to log terms) the rate $n^{5/4}$ for estimating the matrix $M$, thereby improving over the previous $n^{3/2}$ barrier.

Recently, this rate was improved by Liu and Moitra [14] in a specific instance of the problem where $n = d$ and one has access to a sub-polynomial number of noisy independent samples of the complete matrix $M$ from (1) – which is akin to our Poissonian scheme for $\lambda$ being sub-polynomial in $n, d$. They introduce a polynomial-time procedure achieving the rate $n^{1+o(1)}$ for permutation recovery and matrix estimation which, up to the factor $n^{o(1)}$, turns out to be minimax optimal for both problems. As a consequence, in this very specific instance, the computational gap turns out to be nonexistent.

There remain important open problems to characterize the estimation of $\pi^*$ and $M$ in crowdsourcing problems.

* Beyond the case $n \ge d$ handled by Mao et al. [17], the minimax risk $\mathcal{R}^*[n, d, \lambda, 1]$ of estimation of the permutation $\pi^*$ - as well as the minimax risk $\mathcal{R}_{est}^*[n, d, \lambda, 1]$ of estimation of the matrix $M$ - are unknown. In particular, in the rectangular case where $n \ll d$, the number of questions exceeds the number of experts is both relevant from a practical [27] and a conceptual perspective. Indeed, the analysis of the least-square estimator of [17] and related works is based on entropy calculation of the class of permuted bi-isotonic matrices. While the minimax risk $\mathcal{R}_{est}^*[n, d, \lambda, 1]$ turns out to be (up to logarithm terms), characterized by this entropy, this is not always the case for the estimation of $\pi^*$ as many matrices $M$ share the same permutation $\pi^*$. As a consequence, even if we leave aside computational constraints, pinpointing the optimal risk $\mathcal{R}^*[n, d, \lambda, 1]$ for estimating $\pi^*$ requires quite different arguments.

* Beyond the toy "over-complete" observation model in the square case $n = d$ of Liu and Moitra [14], it remains unclear whether there is a computational gap for general rectangular settings with partial observations.

## 1.3. Our Contributions

Echoing with these open problems, we make the following contributions in this work:

* First, we characterize (up to polylogarithmic multiplicative terms) the minimax risk $\mathcal{R}^*[n, d, \lambda, \zeta]$ of permutation recovery, this, for all possible number of experts $n \ge 1$, number of questions $d \ge 1$, noise level $\zeta \ge 0$, and almost all sampling efforts $\lambda > 0$. When $n \ll d$, we prove in particular that $\mathcal{R}^*[n, d, \lambda, 1] \ll \mathcal{R}_{est}^*[n, d, \lambda, 1]$ in all non-trivial regimes, highlighting that when $n \ll d$, the problem of permutation recovery is statistically easier than the problem of matrix estimation.

* Moreover, we introduce a polynomial-time procedure achieving this risk bound, thereby establishing that there does not exist any significant computational-statistical trade-off for the problem of recovering a single permutation $\pi^*$. While our procedure borrows some of the ingredients of Liu and Moitra [14], we need to introduce several new ideas to deal with the significantly more involved case $n \ll d$. Since an estimator $\hat{\pi}$ of $\pi^*$ can be easily combined with a least-square estimator of a bi-isotonic matrix to estimate the matrix $M$ –see e.g. [25, 17] – we also deduce a polynomial time estimator $\widehat{M}$ which nearly achieves the minimax estimation risk $\mathcal{R}_{est}^*[n, d, \lambda, 1]$, thereby proving that this problem
<table><caption>Figure 1: Summary of the minimax risks (up to poly-logarithmic terms) for permutation estimation ($R^*[n, d, 1, 1]$) and matrix estimation ($R_{est}^*[n, d, 1, 1]$) in the specific cases where $\lambda, \zeta = 1$.</caption><thead><tr><th></th><th>n &le; d<sup>1/3</sup></th><th>d<sup>1/3</sup> &le; n &le; d</th><th>n &ge; d</th></tr></thead><tbody><tr><td>Permutation estimation: R<sup>*</sup>[n, d, 1, 1]</td><td>nd<sup>1/6</sup></td><td>n<sup>3/4</sup>d<sup>1/4</sup></td><td>n</td></tr><tr><td>Matrix estimation: R<sub>est</sub><sup>*</sup>[n, d, 1, 1]</td><td>nd<sup>1/3</sup></td><td>&radic;nd</td><td>n</td></tr></tbody></table>

does not either exhibit any computational-statistical trade-off, thereby answering the open problem of [17].

To provide a glimpse of our results, let us describe the minimax risks on the arguably most interesting case where the noise level $\zeta$ is of order 1 as in the Bernoulli observation setting and where $\lambda < 1$ which corresponds to a partially observed matrix. In Section 4, we establish that the minimax risk $R^*[n, d, \lambda, 1]$ of permutation recovery is (up to polylogarithmic multiplicative terms) of the order of

$$ \left[ \frac{nd^{1/6}}{\lambda^{5/6}} \wedge \frac{n^{3/4}d^{1/4}}{\lambda^{3/4}} + \frac{n}{\lambda} \right] \wedge nd, \qquad (6) $$

whereas the minimax reconstruction risk $R_{est}^*[n, d, \lambda, 1]$ is of the order

$$ \left[ \sqrt{\frac{nd}{\lambda}} \wedge \frac{nd}{\lambda^{2/3}(n \vee d)^{2/3}} + \frac{n}{\lambda} \right] \wedge nd. \qquad (7) $$

We display in Figure 1 a summary of our results in the specific case where we also have $\lambda = 1$ on top of $\zeta = 1$, and will discuss this case more in details, as it highlights one of our main findings.

A first comment is that the minimax risk of matrix estimation $R_{est}^*[n, d, 1, 1]$ can be interpreted through the covering numbers of the space of permuted bi-isotonic matrices as in [17]. For $n \ge d$ both minimax risks - $R^*[n, d, 1, 1]$, $R_{est}^*[n, d, 1, 1]$ - are of the order of $n$ so that recovering the permutation $\pi^*$ is as hard as estimating the matrix $M$ (up to logarithmic factors). This is the regime studied in the literature, see [17, 14]. When the number $d$ of questions is large - $n \ll d$ - then the regimes are more tricky. There are two of them, depending on whether $n$ is larger than $d^{1/3}$ or not, and in both regimes $R_{est}^*[n, d, 1, 1]$ is significantly larger than $R^*[n, d, 1, 1]$. More regimes appear when we do not restrict ourselves to $\lambda = 1, \zeta = 1$. This complex picture, as well as the fact that $R^*[n, d, 1, 1] \ll R_{est}^*[n, d, 1, 1]$ for $n \ll d$ - and also in many other configurations of $\lambda, \zeta$ - highlights the fact that the difficulty of estimating $\pi^*$ is not governed by the size of the space of permuted bi-isotonic matrices. As a consequence, even if we leave computational aspects aside, it is not clear that the least-square estimator of [17] achieves optimal risk for estimating the permutation $\pi^*$ and, in any case, entropy-based arguments would lead to suboptimal bounds, at least if we use the same arguments as in [17].

As a byproduct of our results, we also establish the minimax risk - and prove that it is achievable in polynomial time - for another loss function termed $l_\infty(\hat{\pi}, \pi^*)$ (see (30)) put forward in [8, 26, 17] - and we also disprove a conjecture regarding a computational-statistical gap for this loss. See Subsection 4.4.

As our minimax results remain valid in the noiseless case ($\zeta = 0$) where one has access to partial observation of the matrix $M$ itself, we are able to tightly decipher the approximation error which is due to the partial sampling of the matrix $M$ from the stochastic error stemming from noisy observations. In some way, this complements the works of Pananjady et al. [20] on the effect of the design in the specific case where the sampled entries are sampled uniformly.
<figure><img src="image_1.png" /><figcaption>Figure 2: Example of a hierarchical sorting tree.</figcaption></figure>

1.4. *Proof techniques and further comparison with the literature*

In order to build a polynomial-time procedure nearly achieving the minimax permutation risk
in the partial observation setting (2), we first consider the so-called full observation setting
where one has access to poly-logarithmic number $\Upsilon$ of samples $Y^{(0)}, Y^{(1)}, \dots, Y^{(\Upsilon)}$ of the
complete $n \times d$ matrix. This setting is akin to that of Liu and Moitra [14] when they handled
the specific square case $n=d$ with noise level $\zeta = 1$.

For this reason, our estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$ introduced in Section 3 share some features
with the procedure of [14]. From a broad perspective, our procedure and theirs build a hier-
archical sorting tree using a top-down approach as depicted in Figure 2. We start from the
complete set $[n]$ of all experts and build a trisection $(O, P, I)$ of $[n]$, where $O$ (resp. $I$) contains
experts that provably are below (resp. above) the median expert, whereas $P$ contains all the
experts for which we cannot certify with high confidence whether they are above or below the
median. Then, we recursively trisect the sets $O$ and $I$ as depicted in Figure 2. At the end of
the algorithm, we obtain a partial ordering on all the experts which can be used to estimate
the oracle permutation $\pi^*$.

Then, the problem of building a suitable estimator boils down to introducing a suitable
trisection procedure. We could naively do this by comparing the row-sums of the observed
matrix $Y$ which amounts to comparing the mean ability of each expert, but this is well known
to lead to suboptimal performances by a factor $\sqrt{d}$ –see e.g. [17]. To improve over this rate,
we need to compare the experts according to convex combinations of suitable questions. As
in [14], we start by selecting suitable blocks of questions by detecting the high-variation regions
of the mean empirical expert and combine them with spectral algorithms to select suitable
convex combinations of questions. Still, we have to refine significantly their spectral procedure
to handle the rectangular case $n \ll d$. Equipped with these refinements, which are involved
technically, but are built on the ideas developed in [14], we arrive at the estimator $\hat{\pi}_{HT}$ (see
Section 3) that turns out to be minimax optimal in some regimes of $(n, d, \zeta)$.

Unfortunately, this method turns out to be sub-optimal in many regimes, for instance for
mild values of $n \in [d^{1/3}, d]$. Informally, this is due to the fact that our first estimator $\hat{\pi}_{HT}$ as
well as that of Liu and Moitra [14] build an oblivious hierarchical sorting tree. This means that
the trisection method decomposes a group $G^{(0)}$ of experts in the hierarchical sorting tree in
$(O, P, I)$ only using the experts in $G^{(0)}$ of the matrix $Y$. In the related problem of hierarchical
clustering, most top-down procedures also share this feature. It turns out that the observations
of other experts can help improving the trisection of $G^{(0)}$. In particular, sets of experts that are
close in the ordering –such as $G^{(1)}$ and $G^{(-1)}$ in Figure 2– are sometimes valuable to improve
the selection of a suitable convex combination of questions. We emphasize this phenomenon
and provide more intuition on it in Section 3, when we introduce a new estimator $\hat{\pi}_{WM}$ that
builds upon the memory of the sorting tree. This new procedure $\hat{\pi}_{WM}$ turns out to be near
minimax optimal for all values of $(n, d, \zeta)$.

Coming back to the partial observation setting (2), we introduce in Section 4 a reduction
scheme which boils down to reducing the number of questions in order to come back to a full
observation model for a sub-matrix of size $n \times d_-$ where $d_-$ is possibly much smaller than $d$.
Then, relying on the full observation setting described above, we estimate the permutation
$\pi^*$ based on the corresponding reduced matrix. In comparison to the full observation model,
we can suffer from an additional bias terms which arises in the reduction process. To handle
this, we develop a slight variant $\hat{\pi}_{WM-SR}$ of $\hat{\pi}_{WM}$ –see Appendix H for details. The resulting
procedure turns out to nearly achieve minimax permutation recovery risk for all values $(n, d, \zeta)$
and all values of $\lambda$. Plugging our procedure to estimating the matrix $M$, we close all the
computational gaps pointed out in Mao et al. [17] for the problem of matrix estimation with
a single unknown permutation – see Subsection 4.3.

**1.5. Notation and organization of the manuscript**

In the following, $c$, $c_1$, ... stand for numerical positive constants that may change from line to line. Given a vector $u$ and $p \in [1, \infty]$, we write $\|u\|_p$ for its $l_p$ norm. For a matrix $A$, $\|A\|_F$ and $\|A\|_{op}$ stand for its Frobenius and its operator norm. We write $[x]$ (resp. $[x]$) for the largest (resp. smallest) integer smaller than (resp. larger than) or equal to $x$.

Although $M$ stands for an $n \times d$ matrix, we extend it sometimes in an infinite matrix by
setting $M_{i,k} = 0$ when either $i \le 0$ or $k \le 0$ and $M_{i,k} = 1$ when either $i \ge n+1$ and $k > 0$ or
$k \ge d+1$ and $i > 0$. The corresponding infinite matrix $M_{\pi^{*-1}}$ which is obtained by permuting
the $n$ original rows is still bi-isotonic and takes values in $[0,1]$. We shall often work with
sub-matrices of $M$ that are restricted to a subset $P \subset [n]$ and $Q \subset [d]$ of rows and columns, in
which case we write that the corresponding matrix $M'$ belongs to $\mathbb{R}^{P \times Q}$. More precisely, $M'$
is such that, $M'_{i,j} = M_{i,j}$ for any $i \in P$ and any $j \in Q$.

In the following, we write that two sequences or functions $u$ and $v$ satisfy $u \le v$, if there
exists a universal constant such that $u \le cv$.

In Section 2, we first consider the complete observation problem, where one has access
to a poly-logarithmic number of independent samples of the complete noisy matrix $Y$. We
characterize the minimax risk for permutation recovery and prove that it is achieved by a
polynomial-time procedure. In Section 3, we describe the corresponding polynomial-time pro-
cedure. In Section 4, we deal with the problem of partially observed matrix in the model (2).
All the proofs are deferred to the appendix.

**2. Analysis of the full observation problem**

As explained in the introduction, and following [14], we first consider a slightly different prob-
lem where we fully observe a $\Upsilon$-sample $\mathcal{Y} = (Y^{(0)}, \dots, Y^{(\Upsilon-1)})$ of the noisy matrix according to
the model $Y = M + E$ in (1). Here, $\Upsilon$ should be considered as a polylogarithms in $n$ and $d$. This
is of course not very realistic in applications, but it is simpler to first present our algorithmic
procedure in this setting, and it also enables more direct comparison to [14]. We will explain
later in Section 4, how one can transform data in the more realistic partial observation scheme from (2) to this full observation scheme. We will then prove that the algorithm applied to the transformed data is near minimax optimal.

We recall that $M$ is a bi-isotonic matrix, up to an unknown permutation $\pi^*$ of its rows. Besides, the noise matrix $E$ is made of independent mean zero subGaussian entries, with Orlicz norm less than or equal to $\zeta$.

## 2.1. Minimax lower bounds

Before considering ranking procedures, we characterize the minimax risk for the problem of ranking with full observations. For the purpose of the minimax lower bound, we assume that the noise matrix $E$ in (1) is made of independent normal random variables with variance $\zeta^2$. For a permutation $\pi^*$ and a matrix $M$ such that $M_{\pi^*} \in \mathbb{C}_{\text{BISO}}$, we respectively denote $\mathbb{P}_{(\pi^*, M)}$ and $\mathbb{E}_{(\pi^*, M)}$ the corresponding probability and expectations with respect to the $\Upsilon$ independent observations of $Y$. Define

$$ \mathcal{R}_F(n, d, \zeta) = \zeta^2 \left[ \frac{nd^{1/6}}{\zeta^{1/3}} \wedge \frac{n^{3/4}d^{1/4}}{\zeta^{1/2}} \wedge n\sqrt{d} \wedge \frac{n^{2/3}\sqrt{d}}{\zeta^{1/3}} + n \right] . \quad (8) $$

The following minimax lower bound is stated in a setting where one has access to a polylogarithmic number $\Upsilon$ of full samples to be in line with the analysis of the next subsection. Still, we can forget about the dependency in $\Upsilon$ at first reading.

**Theorem 2.1.** There exists a universal constant $c$ such that the following holds for any $n \ge 2$, $d \ge 1$, $\zeta > 0$, and $\kappa > 2$. Provided that the sample size $\Upsilon$ is less than or equal to $\log^\kappa(2nd/\zeta)$, we have

$$ \inf_{\hat{\pi}} \sup_{\pi^* \in \Pi_n} \sup_{M: M_{\pi^{*-1}} \in \mathbb{C}_{\text{BISO}}} \mathbb{E}_{(\pi^*, M)} \|M_{\pi^{-1}} - M_{\pi^{*-1}}\|_F^2 \ge c[\log^{-\kappa}(nd/\zeta)\mathcal{R}_F[n, d, \zeta] \wedge nd] . \quad (9) $$

In fact, this theorem turns out to be a consequence of the minimax lower bound in the partial observation scheme –see Section 4. Together with the risk upper bounds of the next section, (9) characterizes, up to polylogarithmic terms, the minimax risk for estimating $\pi^*$. The term $nd$ in (9) is related to the fact that the loss $\|M_{\pi^{-1}} - M_{\pi^{*-1}}\|_F^2$ cannot be larger than $nd$ because the entries of $M$ are in $[0, 1]$.

**Mild noise level.** The risk bound $\mathcal{R}_F(n, d, \zeta)$ involves five different terms, some of them being significant only when $\zeta$ is small in comparison to $n$ and $d$. As these regimes with very small $\zeta$ are arguably quite specific, and to simplify the discussion, we will now detail the minimax lower bound in the specific case when $\zeta = 1$.

$$ \overline{\mathcal{R}}_F(n, d) = \mathcal{R}_F(n, d, 1) = (nd^{1/6}) \wedge (n^{3/4}d^{1/4}) + n . \quad (10) $$

In particular, we recognize three main regimes in (10) that depend on $n$ and $d$. When the number of experts is relatively small ($n \le d^{1/3}$), the risk is proportional to $nd^{1/6}$. Specifying the result to $n = 2$, one checks that a square distance $d^{1/6}$ is necessary to distinguish two experts. As a consequence, a suitable estimator $\hat{\pi}$ should be able to coherently rank experts that are distant by more than $d^{1/6}$ in squared Frobenius norm, and then to achieve a risk smaller than $nd^{1/6}$. For larger $n \ge d^{1/3}$, it is in fact possible to build upon the large number of experts to improve the comparisons between experts using in particular spectral methods. For this reason, the optimal risk is proportional to $n$ for $n \ge d$. For an intermediary number of experts $n \in [d^{1/3}, d]$, the risk is of the order of $n^{3/4}d^{1/4}$. Our main contribution is the construction of a polynomial-time procedure that achieves these risk bounds, see below.
**Low noise level.** For mild values of $\zeta$, the minimax risk $\mathcal{R}_F(n, d, \zeta)$ has the same form as $\mathcal{R}_F(n, d)$, up to some factors that depend on $\zeta$. However, for very small $\zeta$, the risk becomes qualitatively different. For example, we have $\mathcal{R}_F(n, d, \zeta) \approx \zeta^2 n\sqrt{d}$ when $\zeta \in (0, \frac{1}{n\sqrt{d}}]$. In fact, this rate is quite easy to achieve by a polynomial time algorithm in this extreme case. It is proven in various works – see e.g. [25] that ranking the experts according to the row sum of the matrix correctly compares two experts as long as their square distance is at least $\zeta^2\sqrt{d}$ (up to logarithmic terms). As a consequence, this simple procedure leads to an error $\zeta^2 n\sqrt{d}$. While $\zeta^2 n\sqrt{d}$ is highly suboptimal in most realistic regimes, it turns out to be tight for extremely low level of noise. Finally, the intermediary rate $\mathcal{R}_F(n, d, \zeta) \approx \zeta^{5/3} n^{2/3}\sqrt{d}$ is achieved for slightly larger values of $\zeta$, but it is less clear how to interpret it.

## 2.2. Minimax upper bounds

In the following, we fix a parameter $\delta \in (0, 1)$ that will correspond to a small probability. We write $\zeta_-= \zeta \wedge 1$, where $\zeta$ is the noise level. In this section, we analyze two estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$ of $\pi^*$ that are described in Section 3 and more formally defined in Appendix A. The first estimator $\hat{\pi}_{HT}$ is based on the construction of an *oblivious* hierarchical sorting tree. We will later explain all the ingredients of this procedure. In contrast, the second estimator $\hat{\pi}_{WM}$ relies on the construction of a hierarchical sorting tree with memory. Both procedures have a computational complexity of the order of $\log^c(\frac{nd}{\zeta-\delta})(n^3 + nd^2)$, for some $c > 0$, which makes them polynomial time - unlike the least square procedure e.g. from [17].

**Theorem 2.2.** There exist three numerical constants $c$, $c'$, and $c_0$ such that the following holds. Fix $\delta > 0$ and assume that $\Upsilon \ge c_0 \log^8(\text{nd}/\delta)$. For any permutation $\pi^* \in \Pi_n$ and any matrix $M$ such that $M_{\pi^{*-1}} \in \mathbb{C}_{\text{BISO}}$, the *oblivious* hierarchical sorting tree estimator $\hat{\pi}_{HT}$ defined in the next section satisfies

$$ \| M_{\hat{\pi}_{HT}^{-1}} - M_{\pi^{*-1}} \|_F^2 \le c \zeta^2 \log^{10.5} \left( \frac{2nd}{\delta\zeta_-} \right) \left[ \frac{n^{2/3}d^{1/3}}{\zeta^{2/3}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right], $$

with probability at least $1 - c'n \log^9(\frac{nd}{\delta\zeta})\delta$.

If we take $\delta = \zeta^2(\text{nd})^{-1}$ in the above expression, we easily deduce - reminding that the entries of $M$ are in $[0, 1]$ - the following risk bound

$$ \mathbb{E} \left[ \| M_{\hat{\pi}_{HT}^{-1}} - M_{\pi^{*-1}} \|_F^2 \right] \le c \zeta^2 \log^{10.5} \left( \frac{2nd}{\zeta_-} \right) \left[ \frac{n^{2/3}d^{1/3}}{\zeta^{2/3}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right]. $$

Comparing this bound with (10) in the specific case where $\zeta = 1$, we observe that $\hat{\pi}_{HT}$ achieves the optimal risk $\text{nd}^{1/6}$ for small $n \le d^{1/3}$ and the optimal risk $n$ for large $n \ge d$. Unfortunately, for mild $n \in [d^{1/3}, d]$, the risk bound is of the order of $n^{2/3}d^{1/3}$, which is significantly higher than the minimax lower bound $n^{3/4}d^{1/4}$. To close this gap, we turn to the more refined estimator $\hat{\pi}_{WM}$.

**Theorem 2.3.** There exist three numerical constants $c$, $c'$, and $c_0$ such that the following holds. Fix $\delta > 0$ and assume that $\Upsilon \ge c_0 \log^8(\text{nd}/(\delta\zeta_-))$. For any permutation $\pi^* \in \Pi_n$ and any matrix $M$ such that $M_{\pi^{*-1}} \in \mathbb{C}_{\text{BISO}}$, the hierarchical sorting tree estimator with memory $\hat{\pi}_{WM}$ satisfies

$$ \| M_{\hat{\pi}_{WM}^{-1}} - M_{\pi^{*-1}} \|_F^2 \le \left[ c \log^{11} \left( \frac{2nd}{\delta\zeta_-} \right) \mathcal{R}_F[n, d, \zeta] \right] \wedge nd, \quad (11) $$

with probability at least $1 - c'n \log^9(\frac{nd}{\delta\zeta})\delta$.
As for the previous theorem, this high probability result can be turned into a risk bound by taking $\delta = \zeta^2/(nd)$. In particular, this risk bound matches, up to polylogarithmic terms, the minimax lower bound (9) for all possible values of $n$, $d$, and $\zeta$. As a consequence, the estimator $\hat{\pi}_{WM}$ is nearly minimax and this ranking problem does not exhibit any computational gap.

In [14], the polynomial-time estimator $\hat{\pi}_{LM}$ of Liu and Moitra achieves the minimax risk in the specific square where $n = d$ and $\zeta = 1$. In all the other regimes, no polynomial-time procedure was previously proved to achieve the minimax risk. In fact, even if we do not restrict our attention to polynomial-time procedures, least-square type procedures studied e.g. in [17] provably achieve the minimax risk only in the regime when $n \ge d$. As alluded in the introduction - see Equations (6) and (7), the minimax risks for estimating $\pi^*$ and $M$ differ when $n \le d$, so that achieving the optimal risk for $\pi^*$ is not possible using the classical entropy arguments as in [25, 17]. This highlights the fact that estimating the permutation $\pi^*$ is significantly more challenging in the regime $n \le d$ - both from a statistical and computational perspective - than in the regime $n \ge d$ handled in [14, 17].

**Consequences for the estimation of the matrix $M$.** Provided that we have estimated $\pi^*$ with $\Upsilon-1$ independent samples, we could use the last sample $Y^{(\Upsilon)}$ to estimate the matrix $M$ by minimizing the least-square criterion $\hat{B} = \arg \min_{B \in \mathcal{B}_{\text{BISO}}} \|Y_{\hat{\pi}_{WM}}^{(\Upsilon)} - B\|_F^2$ and setting $\tilde{M} = \hat{B}_{\hat{\pi}_{WM}^{-1}}$. Since the set of bi-isotonic matrices is convex, this estimator is computable efficiently [13]. As argued in Proposition 3.3 of [17] and often used in the ranking literature [26, 8, 21], it turns out that, with high probability, the reconstruction error $\|M - \tilde{M}\|_F^2$ is (up to polylogarithmic terms) the sum of the expected permutation loss $\mathbb{E}[\|M_{\hat{\pi}_{WM}} - M_{\pi^*}\|_F^2]$ and the minimax reconstruction risk of a bi-isotonic matrix $\inf_{\hat{B}} \sup_{B \in \mathcal{B}_{\text{BISO}}} \mathbb{E}[\|\hat{B} - B\|_F^2]$ where $Y = B + E'$ and $E'$ is made of independent subGaussian random variables. Hence, based on $\hat{\pi}_{WM}$ and Theorem 2.3, it is easy to construct a polynomial-time estimator of $M$ that is also near minimax-optimal in the sense of Equation (5). We will further build upon this remark in Section 4 when we come back to the problem of partial observations of the matrix.

### 3. Description of the hierarchical sorting estimators

Let us now describe the construction of the estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$ of $\pi^*$. The construction is quite long and involves several subroutines. For this reason and to ease the understanding of proof details, we also provide a more formal and longer definition in Appendix A. Afterwards, we comment on the different steps of the procedure and on their connection to the literature in Subsection 3.3.

Define $\tau_\infty = [4 \cdot 10^7 \log^7(\frac{nd}{\delta(\zeta)^2})]$ and $t_\infty = [\log(n)/\log(2)]$. We define $\Upsilon^* = 6\tau_\infty t_\infty$ for the total number of independent samples required for the computation of these two estimators.

Hence, we are given independent samples $\mathcal{Y} = (Y^{(0)}, \dots, Y^{(\Upsilon^*-1)})$. From a broad perspective, both procedures are based on the construction of the recursive sorting tree as illustrated in Figure 2. Starting from the root of the tree which corresponds to the set $[n]$ of all experts, we build a partition $O, P, I$, of $[n]$ in such a way that, with high probability, all the experts in $O$ are below the median expert of $[n]$, all the experts in $I$ are above the median expert of $[n]$, while the remaining experts in $P$ are those for which we are not able to decipher whether they are below or above the median expert of $[n]$.

Having trisected $[n]$, we recursively trisect the subsets $O$ and $I$- see Figure 2. Each time, the size of the groups $O$ and $I$ is divided by at least 2. Hence, at depth $t_\infty$, all the groups of $O$ and $I$ have size at most 1. For each depth $t = 0, \dots, t_\infty - 1$, we use $6\tau_\infty$ new samples. The construction of the tree is described in **TreeSort** -see Algorithm 1 and is based on the routine **BlockSort** which performs the trisection of a group into $(O, P, I)$.
Let us now explain how to deduce an estimator $\hat{\pi}$ from the final hierarchical sorting tree $\mathcal{T}$. Indeed, the hierarchical sorting tree $\mathcal{T}$ induces an order on its leaves as follows. For any groups $(O, P, I)$ sharing the same parent, we say that any descendent of $O$ in the tree $\mathcal{T}$ is below $P$, which, in turn, is below any descendent of $I$ in $\mathcal{T}$. This endows a complete ordering on the leaves of the tree $\mathcal{T}$. Denote $\mathcal{G} = (G_1, \dots, G_\alpha)$ the sequence of leaves of the final tree ranked according to this complete order. For any $a \in [\alpha]$, we define the lower bound $\pi_{\bar{\mathcal{G}}}(G_a)$ and the upper bound $\pi_{\bar{\mathcal{G}}}^+(G_a)$ of the ranks of experts in $G_a$ by $\pi_{\bar{\mathcal{G}}}(G_a) := \sum_{a' < a} |G_{a'}|$ and $\pi_{\bar{\mathcal{G}}}(G_a) := \sum_{a' \le a} |G_{a'}|$. Finally, we sample $\hat{\pi}$ arbitrarily in such a way that

$$ \hat{\pi}(G_a) = [\pi_{\bar{\mathcal{G}}}(a) + 1, \pi_{\bar{\mathcal{G}}}^+(a)] . \qquad (12) $$

In other words, the estimator $\hat{\pi}$ ranks the groups $G_a$ according to the ordering of the groups endowed by $\mathcal{T}$ and, given that, ranks the experts $G_a$ uniformly at random. See Appendix A for a more formal definition of the ordering.

### 3.1. Description of the trisection of a leaf G into (O, P, I) with BlockSort

The purpose of **BlockSort** is to build a trisection of a group $G$ of experts into $(O, P, I)$ where $O$ is made of experts that are, with high probability, below the median expert in $G$ and $I$ is made of experts which are, with high probability, above this median expert. It turns out that this construction is based on $\tau_\infty$ iterations of a procedure called **DoubleTrisection** which is the backbone of our procedure. Intuitively, we shall iteratively detect subgroups of experts that are below (resp. above) the median expert of $G$ which, after $\tau_\infty$ iterations, will allow us to obtain $O$ and $I$.

For technical reasons, our definition is slightly more intricate. We shall simultaneously build two collections $(O_\tau, I_\tau)$ and $(\bar{O}_\tau, \bar{I}_\tau)$ of groups, the second one being more conservative. We start with empty sets for $(O_0, I_0, \bar{O}_0, \bar{I}_0) = \emptyset$. Then, at each step $\tau$, we will consider the remaining set of experts $G \setminus (\bar{O}_\tau \cup \bar{I}_\tau)$. Define $\gamma = \lfloor |G|/2 \rfloor - |\bar{O}_\tau|$ for the presumed rank of the median expert of $G$ inside $G \setminus (\bar{O}_\tau \cup \bar{I}_\tau)$. Then, using 6 independent samples, we apply **DoubleTrisection**(Y, T, G, $(\bar{O}_\tau \cup \bar{I}_\tau), \gamma$) to compute four subsets $(L_\tau, U_\tau)$ and $(\bar{L}_\tau, \bar{U}_\tau)$. With high probability, it turns out that $\bar{L}_\tau \subset L_\tau$ is made of experts below the median expert of $G$ and $\bar{U}_\tau \subset U_\tau$ is made of experts above the median expert of $G$. This allows us to update as follows

$$ O_{\tau+1} = O_{\tau} \cup L_{\tau}, \quad I_{\tau+1} = I_{\tau} \cup U_{\tau}, \quad \bar{O}_{\tau+1} = \bar{O}_{\tau} \cup \bar{L}_{\tau}, \quad \bar{I}_{\tau+1} = \bar{I}_{\tau} \cup \bar{U}_{\tau} . \qquad (13) $$

The procedure is summarized in Algorithm 2 below.

**Algorithm 1 TreeSort(Y)**
**Require:** $6\tau_\infty t_\infty$ samples $Y = (Y^{(0)}, \dots, Y^{(6\tau_\infty t_\infty - 1)})$
**Ensure:** A tree $\mathcal{T}$ and an estimator $\hat{\pi}$
1: Initialize $\mathcal{T}$ as the tree with only the root $[n]$
2: **for** $t = 0, \dots, t_\infty - 1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Take $6\tau_\infty$ samples $Y_t = (Y^{(6t\tau_\infty)}, \dots, Y^{(6(t+1)\tau_\infty - 1)})$
4: &nbsp;&nbsp;&nbsp;&nbsp;Initialize $\mathcal{T}' = \mathcal{T}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** All the leaves $G$ at depth $t$ corresponding to $O$ or $I$ as in Figure 2 **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $(O_G, P_G, I_G) = \text{BlockSort}(Y_t, \mathcal{T}, G)$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Add $(O_G, P_G, I_G)$ to the tree $\mathcal{T}'$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update $\mathcal{T} = \mathcal{T}'$
10: **end for**
11: Set $\hat{\pi} := \hat{\pi}(\mathcal{T})$ as in (12)
12: **return** $\mathcal{T}$ and $\hat{\pi}$

**Algorithm 2 BlockSort(Y, T, G)**
**Require:** $6\tau_\infty$ samples $Y = (Y^{(0)}, \dots, Y^{(6\tau_\infty - 1)})$, the tree $\mathcal{T}$, a leaf $G$ in $\mathcal{T}$
**Ensure:** A partition of $G$ into $(O, P, I)$
1: Set $\gamma = \lfloor |G|/2 \rfloor$ and $O_0, I_0, \bar{O}_0, \bar{I}_0 = \emptyset$
2: **for** $\tau = 0, \dots, \tau_\infty - 1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Take 6 samples $Y_\tau = (Y^{(6\tau)}, \dots, Y^{(6\tau+5)})$
4: &nbsp;&nbsp;&nbsp;&nbsp;set $\gamma = \lfloor |G|/2 \rfloor - |\bar{O}_\tau|$ 
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$(L_\tau, U_\tau), (\bar{L}_\tau, \bar{U}_\tau) = \text{DoubleTrisection}(Y_\tau, T, G \setminus (\bar{O}_\tau \cup \bar{I}_\tau), \gamma)$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;as in Algorithm 3
7: &nbsp;&nbsp;&nbsp;&nbsp;Update $O_{\tau+1} = O_\tau \cup L_\tau$, $I_{\tau+1} = I_\tau \cup U_\tau$, $\bar{O}_{\tau+1} = \bar{O}_\tau \cup \bar{L}_\tau$, $\bar{I}_{\tau+1} = \bar{I}_\tau \cup \bar{U}_\tau$
8: **end for**
9: **return** $(O_{\tau_\infty}, G \setminus (O_{\tau_\infty} \cup I_{\tau_\infty}), I_{\tau_\infty})$
## 3.2. Description of the double trisection procedure

We now describe the trisection procedure **DoubleTrisection**. For this purpose, we first provide a few definitions.

### 3.2.1. Definitions

In this subsection, we write $Y$ for one data set sampled according to (1). For the sake of simplicity, we often omit the dependence of $Y$ in the definitions. We write $\mathcal{D}$ for the set of all dyadic numbers: $\mathcal{D} = \{2^k : k \in \mathbb{Z}\}$ and we define the sets $\mathcal{R} = \mathcal{D} \cap [1, d]$ and $\mathcal{H} = \mathcal{D} \cap [\frac{\zeta^2}{nd}, 1]$. The collection $\mathcal{R}$ corresponds to the possible scales, that is the number of questions under consideration, whereas the collection $\mathcal{H}$ corresponds to the possible heights of variations.

For all $r \in \mathcal{R}$, we write $\mathcal{Q}_r = \{1, r+1, 2r+1, \dots, \lfloor \frac{d}{r} \rfloor r+1\}$ for the regular grid of questions with spacing $r$. If $\bar{P} \subset [n]$ is a set of experts, we denote $\bar{y}(\bar{P})$ as the mean of the vectors $Y_{i,}$ for $i \in \bar{P}$, that is, for all $k \in [d]$, we have $\bar{y}_k(\bar{P}) = \frac{1}{|\bar{P}|} \sum_{i \in \bar{P}} Y_{i,k}$. For any $r \in \mathcal{R}$, we define $Z(Y, \bar{P}, r)$ as the aggregation of the matrix $Y$ on blocks of questions of size $r$ and with lines restricted to $\bar{P}$. More formally, for any $i \in \bar{P}$ and $l \in \mathcal{Q}_r$, we have

$$ Z_{i,l}(Y, \bar{P}, r) = \frac{1}{\sqrt{r}} \sum_{k=l}^{l+r-1} Y_{i,k} \quad \text{and} \quad \bar{Z}_{i,l}(Y, \bar{P}, r) = \frac{1}{\sqrt{r}} \sum_{k=l}^{l+r-1} \bar{y}_k(\bar{P}) . \quad (14) $$

Both matrices are of size $|\bar{P}| \times |\mathcal{Q}_r|$. Note that, in the above definition, $Z_{i,l}(Y, \bar{P}, r)$ and $\bar{Z}_{i,l}(Y, \bar{P}, r)$ are rescaled by $\sqrt{r}$ so that the subGaussian norm remains at most $\zeta$. For any subset $Q \subset \mathcal{Q}_r$, we also write $Z(Y, \bar{P}, Q, r)$ for the sub-matrix of $Z(Y, \bar{P}, r)$ restricted to columns in $Q$.

Given a matrix $Z \in \mathbb{R}^{\bar{P} \times Q}$, a vector $w \in \mathbb{R}_+^Q$ with non-negative components and $i, j$ in $\bar{P}$, we say $i$ is $(Z, w)$-above $j$ (or equivalently that $j$ is $(Z, w)$-below $i$) if the projection of $Z_{i,}$ on the direction $w$ is larger than the projection of $Z_{j,}$ on $w$, that is $\langle Z_{i,}, -Z_{j,}, w \rangle > 0$, where $\langle \cdot, \cdot \rangle$ stands for the standard inner product between vectors. Now, for $\gamma \in \{1, \dots, |\bar{P}|\}$, we can consider the $\gamma$-th expert $i_\gamma \in \bar{P}$ such that there are exactly $\gamma-1$ experts which are $(Z, w)$-below $i_\gamma$. Given a tuning parameter $\beta > 0$ to be fixed below, we then define the $(Z, w, \gamma, \beta)$-trisection of $\bar{P}$ on direction $w$ with respect to pivot index $\gamma$ and matrix $Z$ as the sets:

$$ \left\{ \begin{array}{l} U_w := U(Z, w, \gamma, \beta) = \left\{ i \in \bar{P} : \langle Z_{i,}, -Z_{i_\gamma,}, \frac{w}{\|w\|_2} \rangle \ge \beta \sqrt{\log\left(2\frac{|\bar{P}|}{\delta}\right)} \right\} \\ L_w := L(Z, w, \gamma, \beta) = \left\{ i \in \bar{P} : \langle Z_{i,}, -Z_{i_\gamma,}, \frac{w}{\|w\|_2} \rangle \le -\beta \sqrt{\log\left(2\frac{|\bar{P}|}{\delta}\right)} \right\} \end{array} \right. . \quad (15) $$

Hence a $(Z, w, \gamma, \beta)$-trisection on direction $w$ and pivot $\gamma$ consists of two possibly empty disjoint subsets $U$ and $L$ which are respectively taken among the $\gamma-1$ experts (resp. the $|\bar{P}| - \gamma$) which are $(Z, w)$-above (resp. $(Z, w)$-below) the expert $i_\gamma$, with a margin of the order of $\sqrt{\log(|\bar{P}|/\delta)}$. Remark that if $\beta < \bar{\beta}$ then $U(Z, w, \gamma, \bar{\beta}) \subset U(Z, w, \gamma, \beta)$, which means that the trisection of $\bar{P}$ on direction $w$ becomes more conservative as $\beta$ increases.

In fact, (15) turns out to be the cornerstone or our procedure. Since the coordinates of $w$ are non-negative, the corresponding row-wise weighted sums of the aggregation $\mathbb{E}[Z(Y, \bar{P}, r)]$ of the signal matrix $M$ are also ordered according to the oracle permutation. In other words, the $k$-th expert in $\bar{P}$ has the $k$-th highest value of the expectation of this weighted sum.

For $r \in \mathcal{R}$ and $Q \subset \mathcal{Q}_r$, choosing $w = \mathbf{1}_Q$ in (15) amounts to trisecting $\bar{P}$ according to the average of the observations over all questions in $\cup_{l \in Q}[l, l+r)$. In that case, we write for simplicity $(L_Q, U_Q) = (L_{\mathbf{1}_Q}, U_{\mathbf{1}_Q})$. When $Q = \mathcal{Q}_r$ and $w = \mathbf{1}_Q$, then (15) simply amounts to ranking
experts according to their average over all the questions. As explained in the introduction, the global average does not lead to optimal performances. This is why most following steps in the algorithm amount to selecting suitable blocks $Q$ of questions and directions $w$.

In the following, the tuning parameters $\beta$ are set as follows.

$$ \beta_{\text{tris}} = 4\sqrt{2}\zeta, \quad \bar{\beta}_{\text{tris}} = 8\sqrt{2}\zeta. \qquad (16) $$

### 3.2.2. Description of the double trisection procedure

Recall that the purpose of **DoubleTrisection** is to select subsets $(L, U)$ and $(\overline{L}, \overline{U})$ of a group $\overline{P}$ of experts in such a way that $\overline{L} \subset L$, $\overline{U} \subset U$, and experts in $L$ (resp. in $U$) are with high probability below (resp. above) the $\gamma$-th expert of $\overline{P}$.

For that purpose, we have 6 independent samples $(Y^{(s)})_{s=1,...,6}$ sampled from (1) at our disposal. Fix any height $h \in \mathcal{H}$ and any scale $r \in \mathcal{R}$. **DoubleTrisection** relies on the following steps also described in Algorithm 3.

1. **Selection of a suitable subset of questions.** Using the first sample $Y^{(1)}$, we first select a subset $\overline{Q} \subset Q_r$. We postpone the definition of the selection procedure to the next subsection. We will introduce two approaches for this $\overline{Q} := \overline{Q}_{cp}(\overline{P}, h, r)$ as in (20) or $\overline{Q} := \overline{Q}_{WM}(\overline{T}, \overline{P}, h, r)$ as in (26). These two definitions respectively correspond to the *oblivious* estimator $\hat{\pi}_{HT}$ and to the estimator with *memory* $\hat{\pi}_{WM}$.

2. **Average-based trisection.** Using the second sample $Y^{(2)}$, we consider the corresponding aggregated matrix $Z^{(2)} := Z(Y^{(2)}, \overline{P}, \overline{Q}, r)$ as defined in (14) which focuses on the selected blocks of questions $\overline{Q}$. Then, we consider experts whose corresponding row sums on $Z^{(2)}$ is unusually large or small. More formally, we compute the $(Z^{(2)}, \mathbf{1}_{\overline{Q}}, \gamma, \beta_{\text{tris}})$-trisection and the $(Z^{(2)}, \mathbf{1}_{\overline{Q}}, \gamma, \bar{\beta}_{\text{tris}})$-trisection of $\overline{P}$ as defined in (15) and where the tuning parameters $\beta_{\text{tris}}$ and $\bar{\beta}_{\text{tris}}$ are defined in (16). This allows us to obtain $(L_{\overline{Q}}, U_{\overline{Q}})$ and $(\overline{L}_{\overline{Q}}, \overline{U}_{\overline{Q}})$.

3. **PCA-based trisection.** Then, we focus on the conservative subset of remaining experts $\tilde{P} = \overline{P} \setminus \overline{L}_{\overline{Q}} \cup \overline{U}_{\overline{Q}}$. Relying on the samples $Y^{(3)}, Y^{(4)}, Y^{(5)}, Y^{(6)}$, we build the corresponding aggregated matrices $Z^{(s)} := Z(Y^{(s)}, \tilde{P}, \overline{Q}, r)$ restricted to the subset $\tilde{P}$ for $s = 3, 4, 5$. In principle, we would like to aim at the right singular value of $\mathbb{E}[Z^{(3)} - \overline{Z}^{(3)}]$ as this would give us a nice direction $w$ on which we could apply (15). For technical reasons to be explained later, we take a roundabout way, by first computing a vector $\hat{v}$ indexed by $\tilde{P}$ which, in principle, is not too far from the left singular value of $\mathbb{E}[Z^{(3)} - \overline{Z}^{(3)}]$. More precisely, we compute $\hat{v}$ as follows

$$ \hat{v} := \hat{v}(\tilde{P}, \overline{Q}, r) = \arg \max_{\|v\|_2 \le 1} \left[ \|v^T(Z^{(3)} - \overline{Z}^{(3)})\|_2^2 - \frac{1}{2} \|v^T(Z^{(3)} - \overline{Z}^{(3)} - Z^{(4)} + \overline{Z}^{(4)})\|_2^2 \right]. \quad (17) $$

The right-hand side term in (17) allows us to deal with the fact that the entries of the noise matrix $E$ in (1) are possibly heteroskedastic. Although there exist more elegant workarounds for heteroskedastic noise (e.g. PCA [34]), the analysis in those works does not apply in our non-parametric setting. Moreover, $\hat{v}$ in (17) corresponds to the leading eigenvector of a square symmetric matrix and can therefore be computed efficiently. Then, we consider the image $\hat{z} = \hat{v}^T(Z^{(5)} - \overline{Z}^{(5)}) \in \mathbb{R}^{\overline{Q}}$ of $\hat{v}$. After this, we threshold $\hat{z}$ and take the absolute values of the components. Thus, we get $\hat{w}^+ \in \mathbb{R}^Q$ defined by $(\hat{w}^+)_l = |\hat{z}_l| \mathbf{1}_{|\hat{z}_l| \ge 2\zeta\sqrt{2\log(2|\overline{Q}|/\delta)}}$ for any $l \in \overline{Q}$. Finally, we consider the last aggregated sample $Z^{(6)} := Z(Y^{(6)}, \overline{P}, \overline{Q}, r)$ on the set $\overline{P} \supset \tilde{P}$ of experts. We apply these weights $\hat{w}^+$
to compute the row-wise weighted sums of $Z^{(6)}$ and discard experts whose corresponding weighted sums is unusually small or large. More formally, we apply $(Z^{(6)}, \hat{w}^+, \gamma, \beta)$-trisection and $(Z^{(6)}, \hat{w}^+, \gamma, \bar{\beta})$-trisection of $\bar{P}$ as defined in (15). Doing so we obtain $(L_{\hat{w}^+}, U_{\hat{w}^+})$ and $(\bar{L}_{\hat{w}^+}, \bar{U}_{\hat{w}^+})$ respectively.

In the definition of $Z^{(6)}$ we consider the whole set of experts $\bar{P}$ instead of the remaining of experts $\bar{P}$ that have not been discarded because otherwise we should have needed to update the value of $\gamma$ when applying (15).

Finally, we define the trisections $(L, U)$ (resp. $(\bar{L}, \bar{U})$) as the union of the corresponding discarded subsets of experts based on $\mathbb{1}_{\bar{Q}}$ and $\hat{w}^+$, this for all possible height $h \in \mathcal{H}$ and scale $r \in \mathcal{R}$. We recall that the definition of $\bar{Q}$ was depending on $h$ and $r$.

$$ \left\{ \begin{aligned} (L, U) &= \left( \bigcup_{(h,r) \in \mathcal{H} \times \mathcal{R}} L_{\bar{Q}}(h, r) \cup L_{\hat{w}^+}(h, r), \bigcup_{(h,r) \in \mathcal{H} \times \mathcal{R}} U_{\bar{Q}}(h, r) \cup U_{\hat{w}^+}(h, r) \right) \\ (\bar{L}, \bar{U}) &= \left( \bigcup_{(h,r) \in \mathcal{H} \times \mathcal{R}} \bar{L}_{\bar{Q}}(h, r) \cup \bar{L}_{\hat{w}^+}(h, r), \bigcup_{(h,r) \in \mathcal{H} \times \mathcal{R}} \bar{U}_{\bar{Q}}(h, r) \cup \bar{U}_{\hat{w}^+}(h, r) \right). \end{aligned} \right. \quad (18) $$

This whole routine for computing $(L, U)$, $(\bar{L}, \bar{U})$ is referred to as **DoubleTrisection** and is summarized in Algorithm 3. We underline that $\bar{L} \subset L \subset \bar{P}$ and $\bar{U} \subset U \subset \bar{P}$ as we took $\beta_{\text{tris}} < \bar{\beta}_{\text{tris}}$.

<caption>**Algorithm 3 DoubleTrisection((Y(s))_{s=1,...,6}, T, \bar{P}, \gamma)**</caption>

**Require:** 6 samples $(Y^{(s)})_{s=1,...,6}$, a set $\bar{P}$, a tree $T$, a pivot index $\gamma \in [1,...,|\bar{P}|]$ 
**Ensure:** Two trisections $(L, U)$ and $(\bar{L}, \bar{U})$ of $\bar{P}$

1: Start from $L, U, \bar{L}, \bar{U} = \emptyset$

2: **for** $h \in \mathcal{H}$, $r \in \mathcal{R}$ **do**

3: Compute $\bar{Q} := \bar{Q}_{\text{cp}}(\bar{P}, h, r)$ as in (20) or $\bar{Q} := \bar{Q}_{\text{WM}}(T, \bar{P}, h, r)$ as in (26) using sample $Y^{(1)}$

4: Let $Z^{(s)} := Z(Y^{(s)}, \bar{P}, \bar{Q}, r)$, for $s \in \{2,6\}$ be the aggregated matrices of samples defined as in (14)

5: Let $(L_{\bar{Q}}, U_{\bar{Q}})$, $(\bar{L}_{\bar{Q}}, \bar{U}_{\bar{Q}})$ be resp. the $(Z^{(2)}, \mathbb{1}_{\bar{Q}}, \gamma, \beta)$ and the $(Z^{(2)}, \mathbb{1}_{\bar{Q}}, \gamma, \bar{\beta})$-trisections of $\bar{P}$ as in (15)

6: Define $\tilde{P} = \bar{P} \setminus (\bar{L}_{\bar{Q}} \cup \bar{U}_{\bar{Q}})$ and the aggregated samples $Z^{(s)} := Z^{(s)}(Y^{(s)}, \tilde{P}, \bar{Q}, r)$ for $s \in \{3,4,5\}$

7: Compute the PCA-like direction $\hat{v} := \hat{v}(\tilde{P}, \bar{Q}, r)$ as in (17)

8: Compute $\hat{z} = \hat{v}^T(Z^{(5)} - \bar{Z}^{(5)})$ and define $\hat{w}^+$ by $(\hat{w}^+)_l = |\hat{z}_l|\mathbb{1}_{|\hat{z}_l| \ge 2\zeta\sqrt{2\log(2|\bar{Q}|/\delta)}}$ for any $l \in \bar{Q}$

9: Let $(L_{\hat{w}^+}, U_{\hat{w}^+})$, $(\bar{L}_{\hat{w}^+}, \bar{U}_{\hat{w}^+})$ be resp. the $(Z^{(6)}, \hat{w}^+, \gamma, \beta)$ and the $(Z^{(6)}, \hat{w}^+, \gamma, \bar{\beta})$-trisections of $\bar{P}$ as in (15)

10: Update $L = L \cup L_{\hat{w}^+} \cup L_{\bar{Q}}$, $U = U \cup U_{\hat{w}^+} \cup U_{\bar{Q}}$, $\bar{L} = \bar{L} \cup \bar{L}_{\hat{w}^+} \cup \bar{L}_{\bar{Q}}$, $\bar{U} = \bar{U} \cup \bar{U}_{\hat{w}^+} \cup \bar{U}_{\bar{Q}}$

11: **end for**

12: **return** $(L, U)$, $(\bar{L}, \bar{U})$

To finish the definition of the estimator, it remains to describe the selection procedures for the suitable blocks of questions that are used in Line 3 of Algorithm 3. As explained above, we consider two procedures $\bar{Q} := \bar{Q}_{\text{cp}}(\bar{P}, h, r)$ as in (20) or $\bar{Q} := \bar{Q}_{\text{WM}}(T, \bar{P}, h, r)$ as in (26) - which respectively apply to the **oblivious estimator** $\hat{\pi}_{\text{HT}}$ and to the estimator $\hat{\pi}_{\text{WM}}$ with **memory**.

### 3.2.3. Definition of $\bar{Q}_{\text{cp}}$

We start with $\bar{Q}_{\text{cp}}(\bar{P}, h, r)$. The corresponding estimator $\hat{\pi}_{\text{HT}}$ is called an **oblivious hierarchical sorting tree estimator** because $\bar{Q}_{\text{cp}}$ only depends on the restriction of the data to $\bar{P}$. As a consequence, the corresponding **BlockSort** procedure (see Algorithm 2) which builds a trisection of a group $G$ of experts into three subgroups $(O, P, I)$ only depends on the observations on this set $G$ of experts. In other words, the recursive construction of the hierarchical sorting tree estimator is completely oblivious of the rest of the tree. Up to our knowledge, this feature is shared by most hierarchical clustering algorithms.
Fix some height $h \in \mathcal{H}$ and $r \in \mathcal{R}$. Intuitively, $\tilde{Q}_{cp}(\bar{P}, h, r)$ amounts to focusing on the subset of questions around which the empirical mean expert $\bar{y}(\bar{P})$ has a high-variation. We provide some intuition on the rationale of this approach in the next subsection. More precisely, we define the CUSUM statistic:

$$ \tilde{r} = 8 \left[ \left( \frac{32 \zeta^2}{|\bar{P}| h^2} \log\left(\frac{2d}{\delta}\right) \right) \vee r \right] \quad \text{and} \quad \tilde{\mathbf{C}}_{k, \tilde{r}}(\bar{P}) = \frac{1}{\tilde{r}} \left( \sum_{k'=k}^{k+\tilde{r}-1} \bar{y}_{k'}(\bar{P}) - \sum_{k'=k-\tilde{r}}^{k-1} \bar{y}_{k'}(\bar{P}) \right) . \quad (19) $$

In a nustshell, $\tilde{\mathbf{C}}_{k, \tilde{r}}(\bar{P})$ is the empirical variation of $\bar{y}(\bar{P})$ at question $k$ and at scale $\tilde{r} \ge r$. Then, we define $\tilde{D}_{cp}$ as the set of questions where the CUSUM statistic is larger than $h/4$, and $\tilde{Q}_{cp} \subset Q_r$ for the corresponding subset of blocks or questions of size $r$.

$$ \tilde{D}_{cp} = \{k \in [d] : \tilde{\mathbf{C}}_{k, \tilde{r}}(\bar{y}(\bar{P})) \ge \frac{h}{4}\} \quad \text{and} \quad \tilde{Q}_{cp} = \{l \in Q_r : \tilde{D}_{cp} \cap [l, l+r) \ne \emptyset\} . \quad (20) $$

In (19), the choice of $\tilde{r} \ge r$ is due to the fact that we need to compute an empirical mean $\mathbf{C}_{k, \tilde{r}}(\bar{y}(\bar{P}))$ on enough questions so that its standard deviation is small compared to $h$.

### 3.2.4. Definition of $\tilde{Q}_{WM}$

Finally, we describe $\tilde{Q}_{WM}(\mathcal{T}, \bar{P}, h, r)$ which corresponds to the estimator $\hat{\pi}_{WM}$. The set $\bar{P}$ is a subset of a leaf $G$ of the tree $\mathcal{T}$ and we write $t$ for its depth. As illustrated in Figure 2, there is a natural order on the nodes of $\mathcal{T}$ at depth $t$ that have been either obtained as subsets of type $O$ or $I$ in BlockSort (Algorithm 2). We can index these nodes according to the ordering by setting $G^{(0)} = G$ and then $G^{(1)}, G^{(2)}, \dots$ as the following groups. Similarly, $G^{(-1)}, G^{(-2)}, \dots$ stand for the groups preceding $G^{(0)}$. See Figure 2 for an illustration. In fact, with high probability, for any $a$, all the experts in $G^{(a+1)}$ are above the expert in $G^{(a)}$. As a consequence, the observations in $G^{(1)}$ and $G^{(-1)}$ can bring some informations on the behaviour of the experts in $\bar{P} \subset G^{(0)}$.

Fix $r \in \mathcal{R}$ and $h \in \mathcal{H}$. Define $\tilde{r} \in \mathcal{R}$ as $\tilde{r} = 4([2^9 \log(4d|\mathcal{R}|/\delta) \frac{\zeta^2}{|\bar{P}|h^2}]^{dya} \vee r)$, where $[x]^{dya} = 2^{[\log_2(x)]}$. As before, $\tilde{r} \ge r$ stands for the scale which is required if we want to estimate the variation of $\bar{y}(\bar{P})$ with a standard error small compared to $h$.

Now, we consider any scale $r_{cp} \in [4r, 2\tilde{r}] \cap \mathcal{R}$. The rationale is that, if $r_{cp} < \tilde{r}$, we can reduce the standard deviations of the empirical means by considering an average over experts in neighboring groups. Define the upper neighborhood $\mathcal{V}_{r_{cp}}^+$ and lower neighborhood $\mathcal{V}_{r_{cp}}^-$ as the set of groups above $G$ and below $G$ that are necessary to have enough experts at scale $r_{cp}$.

$$ a_{WM}^{+} = \min\{a : |G^{(1)}| + \dots + |G^{(a)}| \ge \frac{2^{11} \zeta^2 \log(4d|\mathcal{R}|/\delta)}{r_{cp}h^2}\} \quad \text{and} \quad \mathcal{V}_{r_{cp}}^{+} = \bigcup_{a=1}^{a_{WM}^{+}} G^{(a)}; \quad (21) $$

$$ a_{WM}^{-} = \min\{a : |G^{(-1)}| + \dots + |G^{(-a)}| \ge \frac{2^{11} \zeta^2 \log(4d|\mathcal{R}|/\delta)}{r_{cp}h^2}\} \quad \text{and} \quad \mathcal{V}_{r_{cp}}^{-} = \bigcup_{a \in a_{WM}^{-}} G^{(a)}. \quad (22) $$

For a given subset $\bar{P} \subset G$, we define the corresponding CUSUM statistic $\tilde{\mathbf{C}}_{k, r_{cp}}^{(\text{ext})}$ computed on the questions $[k - r_{cp}, k + r_{cp})$ and using the empirical mean observations in $\mathcal{V}_{r_{cp}}^+ \cup \mathcal{V}_{r_{cp}}^-$ if $r < 2\tilde{r}$ and in $\bar{P}$ if $r_{cp} = 2\tilde{r}$:

$$ \tilde{\mathbf{C}}_{k, r_{cp}}^{(\text{ext})} = \frac{1}{r_{cp}} \begin{cases} \sum_{k'=k}^{k+r_{cp}-1} \bar{y}_{k'}(\mathcal{V}_{r_{cp}}^{+} \cup \mathcal{V}_{r_{cp}}^{-}) - \sum_{k'=k-r_{cp}}^{k-1} \bar{y}_{k'}(\mathcal{V}_{r_{cp}}^{+} \cup \mathcal{V}_{r_{cp}}^{-}) & \text{if } r_{cp} \in [8r, 2\tilde{r}) \\ \sum_{k'=k}^{k+r_{cp}-1} \bar{y}_{k'}(\bar{P}) - \sum_{k'=k-r_{cp}}^{k-1} \bar{y}_{k'}(\bar{P}) & \text{if } r_{cp} = 2\tilde{r} \end{cases} \quad (23) $$
If $r_{cp} = 2\tilde{r}$, this new definition of the CUSUM with memory matches the definition (19) in the previous paragraph. For $r_{cp} < 2\tilde{r}$, we are not able to average on enough expert in $\bar{P}$. To deal with this issue, we average on a suitable number of neighboring experts.

Beside considering questions around which the variations of $\bar{y}(\bar{P})$ are large enough, we also check whether, on the corresponding regions, the width of $\bar{P}$, that is the difference between the best expert and the worst expert in $\bar{P}$ is high enough. Given a question $k \in [d]$, we define $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})}$ as the difference between the locals average on $[k - r_{cp}, k + r_{cp})$ of the neighbourhoods of $G$ that is

$$ \hat{\Delta}_{k,r_{cp}}^{(\text{ext})} = \frac{1}{2r_{cp}} \sum_{k'=k-r_{cp}}^{k+r_{cp}-1} \bar{y}_{k'}(\mathcal{V}_{r_{cp}}^{+}) - \bar{y}_{k'}(\mathcal{V}_{r_{cp}}^{-}) . \quad (24) $$

Since the groups $G^{(1)}, \dots, G^{(2)}$ are above the best expert in $\bar{P}$ and since the groups $G^{(-1)}, G^{(-2)}, \dots$ are below the worst expert in $\bar{P}$, this statistic $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})}$ overestimates the width of $\bar{P}$. In the next subsection, we will explain why it is relevant to consider the width of $\bar{P}$.

We are now equipped to define the subsets $\hat{D}_{WM}(T, \bar{P}, h, r)$ of suitable questions and the corresponding $\hat{Q}_{WM}(T, \bar{P}, h, r)$ of corresponding blocks.

$$ \hat{D}_{WM} = \{k \in [d] : \exists r_{cp} \in [4r, \tilde{r}] \cap \mathcal{R} \text{ s.t. } \hat{C}_{k, 2r_{cp}}^{(\text{ext})} \ge \frac{h}{16} \text{ and } \hat{\Delta}_{k, r_{cp}}^{(\text{ext})} \ge \frac{h}{16}\} ; \quad (25) $$

$$ \hat{Q}_{WM} = \{l \in Q_r : \hat{D}_{WM} \cap [l, l+r) \neq \emptyset\} . \quad (26) $$

In other words, $\hat{D}_{WM}$ is made of questions for which there exists a scale $r_{cp}$ such that simultaneously the empirical variations $\hat{C}_{k, 2r_{cp}}^{(\text{ext})}$ at scale $2r_{cp}$ is at least of order $h$ and the empirical width at scale $r_{cp}$ is at least of order $h$.

### 3.3. Comments on the procedure and relation to the literature

These twin procedure are quite involved and combine several ingredients, some of them being already used by Liu and Moitra [15]. In particular, they introduced the key ideas of localization of the suitable blocks of questions through change-point detection on the mean expert and of a spectral clustering scheme for dividing blocks of experts. Still, we need to add several key elements in order to deal with the arguably more involved setting where $n \le d$. We describe below how our procedure compares to [15] and highlight also the main differences and new ideas. Also, despite the fact that our procedure is very involved, it remains computationally efficient. Overall, the full procedure requires $O[\log^c(\frac{nd}{\tilde{c}-\delta})(nd^2 + n^3)]$ operations for some $c > 0$. Indeed, each of the main steps of the algorithm correspond to matrix multiplications and computations of the largest eigenvector of a square symmetric matrix.

In this subsection, we discuss three key steps of the algorithm: (i) the selection of blocks of questions corresponding to the high-variation regions of the average expert in the group as in the definition of $\hat{Q}_{cp}$, (ii) construction of the weights vector $w^+$ by a spectral procedure, (iii) the use of neighboring groups in $\hat{Q}_{WM}$.

#### 3.3.1. Detecting high-variation regions of the average expert

Recall that, for a fixed $r \in \mathcal{R}$ and $h \in \mathcal{H}$, $\hat{Q}_{cp}$ selects blocks of questions in which the variations (19) $\hat{C}_{k, \tilde{r}}(\bar{P})$ of the average $\bar{y}(\bar{P})$ at question $k$ and at scale $\tilde{r} \ge r$ is higher than $h/4$.

To explain the rationale behind this choice, let us first consider a toy example depicted in Figure 3. Assume that the group $\bar{P}$ is made of two subgroups of experts $U^*$ and $L^*$ and that all the experts in $U^*$ and all the experts in $L^*$ are identical. Also, assume that the corresponding
rows only differ on *r* consecutive questions by *h* and are otherwise identical. As illustrated in Figure 3, it turns out that the expected average expert $\bar{m}(\bar{P}) = E[\bar{y}(\bar{P})]$ needs to vary by *h* at scale *r* near the block of questions on which the two groups of experts are differing. This is due to the fact that both the rows corresponding to $U^*$ and $L^*$ are isotonic and that the row of $U^*$ is always larger or equal to that of $L^*$. As a consequence, by restricting our attention to the blocks of questions corresponding to high-variation regions of $\bar{m}(\bar{P})$ (or in practice $\bar{y}(\bar{P})$), we are able to much reduce the dimension of the problem and thereby to improve our ability to distinguish different experts.

Beyond this toy example, we show in Lemma D.3 that there exists a suitable scale $r \in \mathbb{R}$ and a suitable height $h \in \mathcal{H}$ such that, by restricting our attention to blocks of questions of size $r$ such that the expected average expert $\bar{m}(\bar{P})$ varies by at least $h/2$, we are able to retain a significant proportion of the differences between experts in $\bar{P}$. In other words, focusing on regions of high-variation of $\bar{y}(\bar{P})$ in the blocks $\tilde{Q}_{cp}$ is, at least for some scale and some height, a suitable dimension reduction technique. This phenomenon was already observed in [14] and their procedure also uses such dimension detection techniques. In our paper, we also build upon this idea, which has also important consequences, in a related yet different manner, in the rectangular case where $n \le d$.

If we do not apply the spectral clustering sorting steps in $\hat{\pi}_{HT}$, that is, if we do not compute $\hat{v}$ and $\hat{w}^+$ in **DoubleTrisection**, then we would get a risk bound for $\hat{\pi}_{HT}$ of the order of $\zeta^{5/3}nd^{1/6}$ instead of that of Theorem 2.2. In other words, the dimension reduction in $\tilde{Q}_{cp}$ is alone sufficient to recover the optimal risk in the case where $d$ is quite large and $\zeta$ is mild - namely $n \le \zeta^{2/3}d^{1/3}$ and $\zeta \in [1/d, \sqrt{d}]$.

<figure><img src="image_2.png" /><figcaption>Figure 3: In this toy example, the group $\bar{P}$ is only made of two types of experts, those in $U^*$ and those in $L^*$. The high-variation region of $\bar{m}(\bar{P})$ corresponds to the questions on which $U^*$ and $L^*$ differ.</figcaption></figure>

### 3.3.2. On the spectral estimation of the weights

In this subsection, we explain how the computation of $\hat{z}$ in (17) and the corresponding weights $\hat{w}^+$ allow to improve over the $\zeta^{5/3}nd^{1/6}$ rate. Again, we start with a motivating toy example depicted in Figure 4. As previously, we consider a situation where $\bar{P}$ can be decomposed into two subgroups $U^*$ and $L^*$ of the same size. The corresponding rows $L^*$ are block-constant with blocks of questions of size $r$ and increased by $h$ at the end of each block of questions. On the other hand, the corresponding lines of $U^*$ are, in each block of questions, either equal to the rows of $L^*$, or are exactly at a distance $h$ above. These last blocks of questions are
the only ones which are informative when it comes to distinguishing the best experts in the group - namely $U^*$ - from the worst experts in the group - namely $L^*$. Some of the blocks corresponding to high-variation regions of the expected average row $\bar{m}(\bar{P})$ do not convey any information on the difference between $U^*$ and $L^*$ - see Figure 4. In this example, at scale $r$, all the blocks of size $r$ are to be detected by the high variation dimension reduction step, that is $\tilde{Q}_{cp} = Q_r$. At the second step, we consider the corresponding aggregated matrix $Z - \bar{Z}$ at scale $r$ as defined in (14). To be more specific, let us assume that $|L^*| = |U^*| = 3$. Then, $Z - \bar{Z}$ is a $6 \times 8$ matrix whose expectation is of the form of the right panel in Figure 4.

<figure><img src="image_3.png" /><figcaption>Figure 4: In this toy example, the group $\bar{P}$ is only made of two types of experts $U^*$ and $L^*$ with $|L^*| = |U^*| = 3$.</figcaption></figure>

In this specific example, the rank of this expected matrix is one and some of its columns are completely useless to decipher experts in $U^*$ from experts in $L^*$. In contrast, taking $w$ as a right singular vector associated to the largest singular value of this matrix would allow us to select the significant blocks of questions while discarding the irrelevant ones. While this example is very specific, this still sheds some light on why spectral clustering procedure can be of interest for this problem and how it can help recover blocks of questions that are the most informative for dividing the experts.

Let us come back to a general matrix $M$ and to the spectral step of **DoubleTrisection** as described in the previous section. Up to a permutation of its rows, the expectation $\Theta^{(3)} - \bar{\Theta}^{(3)}$ of $Z^{(3)} - \bar{Z}^{(3)}$ is isotonic in each column. It turns out that the entries of any left singular vector associated to the largest singular value of $\Theta^{(3)} - \bar{\Theta}^{(3)}$ is, up to the permutation, either non-increasing or non-decreasing. As a consequence, the left-singular value of $\Theta^{(3)} - \bar{\Theta}^{(3)}$ can bring information on the underlying ranking. This property is at the heart of spectral ranking algorithms [31]. Unfortunately, contrary to the analysis of spectral methods in the Bradley-Luce-Terry model [10, 9], we cannot control the entry-wise deviations of the left singular eigenvector of $Z^{(3)} - \bar{Z}^{(3)}$ because the matrix $\Theta^{(3)} - \bar{\Theta}^{(3)}$ is non-parametric and does not necessarily exhibit any spectral gap. To handle this, Liu and Moitra [14] suggest to compute a right singular vector of $Z^{(3)} - \bar{Z}^{(3)}$ and, using another independent sample, to compare the experts based on the corresponding weighted average of the experts. Unfortunately, while their analysis provides near optimal results for $n = d$, this would not work for $n \le d$. In **DoubleTrisection**, we apply a more involved workaround (i) to handle possible heteroskedastic noise and (ii) to improve the convergence rates in comparison to Liu and Moitra [14]. Indeed, we first compute in (17) a debiased version $\hat{v}$ of the left-singular vector of $Z^{(3)} - \bar{Z}^{(3)}$. Then, we compute the image $[Z^{(5)} - \bar{Z}^{(5)}]^T \hat{v}$, threshold it, and take its absolute value to obtain our estimated weights $\hat{w}^+$. In principle, $\hat{w}^+$ aims at being close to the right first singular vector of $\Theta^{(3)} - \bar{\Theta}^{(3)}$. In comparison
to Liu and Moitra [14], $\hat{w}^+$ better handles the situation where the matrix $\Theta^{(3)} - \bar{\Theta}^{(3)}$ is highly rectangular (with many columns) and where its corresponding right singular vector is nearly sparse.

### 3.3.3. On the tree information and the definition of $\hat{Q}_{WM}$

The oblivious estimator $\hat{\pi}_{HT}$ based on $\hat{Q}_{cp}$ is only proved to achieve the suboptimal error of Theorem 2.2. In this section, we explain how $\hat{Q}_{WM}$ improves the performances of the procedure by relying on the neighboring experts to fix one possible weakness of $\hat{Q}_{cp}$ and so, improve the dimension reduction step.

Indeed, $\hat{Q}_{cp}$ selects spurious blocks of questions. In the previous toy example (Figure 4), some of the blocks corresponding to high-variation values of the expected mean expert $\bar{m}(\bar{P})$ do not bring any suitable information for ordering the experts in $\bar{P}$ because, in these blocks, all the experts are close to each other. In other words, the width of $\bar{P}$, that is the difference between the best and worst experts in $\bar{P}$, is small. It is not possible to easily estimate this width from the observations in $\bar{P}$ since this would require to have sorted the experts in $\bar{P}$ in the first place. Still, we can estimate this width by comparing the average of experts that are above $\bar{P}$ with average of experts that are below $\bar{P}$. A first idea would therefore be to consider a large enough number of experts above and below $\bar{P}$ in order to estimate the width with a small variance and to exclude regions such that the estimated width on a window of size $r$ is small compared to $h$. This is exactly the purpose of the statistic $\hat{\Delta}_{k,r}^{(\text{ext})}$. The selected blocks $\hat{Q}_{WM}$ only contain regions such that the estimated width $\hat{\Delta}_{k,r}^{(\text{ext})}$ is large enough compared to $h$ –see the left panel in Figure 5. Unfortunately, the statistic $\hat{\Delta}_{k,r}^{(\text{ext})}$ may suffer from a large positive bias if the experts above or below $\bar{P}$ are away from $\bar{P}$. Moreover, considering only the scale $r$ is not sufficient because we are forced to average over many experts above and below $\bar{P}$ in order to have a small variance at this small scale, leading to a large bias. For this reason, we consider all possible scales $r_{cp}$ (in a dyadic grid) between $r$ and $\tilde{r}$.

Another important idea in the dimension reduction scheme is the following: If there is a region of questions in which, not only the mean experts of the group $\bar{P}$ but also the mean experts in neighboring groups of $\bar{P}$ have a high variation, it is interesting to detect this high-variation region by relying on all these neighboring groups in order to decrease the variance of the CUSUM statistic. With this idea, we are able to consider the CUSUM statistic at a smaller scale $r_{cp} \le \tilde{r}$ –see the right panel in Figure 5. This is exactly the purpose of the statistic $\bar{C}_{k,r_{cp}}^{(\text{ext})}$.

In our procedure, we build a collection $\hat{D}_{WM}$ that selects a question $k$ if there exists a scale $r_{cp}$ in $[4r, \tilde{r}]$ such that both the CUSUM statistic $\bar{C}_{k,2r_{cp}}^{(\text{ext})}$ at scale $2r_{cp}$ is large and the empirical width $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})}$ is large. This combines the two ideas described in the previous paragraphs which, in turn, allows us to further reduce the dimension in comparison to $\hat{D}_{cp}$ while ensuring that the selected questions in $\hat{D}_{WM}$ contains all the relevant regions to trisect $\bar{P}$, namely regions of size $r$, on which $\bar{P}$ has a variation at least of the order of $h$ and the width of $\bar{P}$ is at least of the order of $h$.

Interestingly, in the square case where $n = d$ considered in [14] or more generally when $n \ge d$, this dimension reduction variant is not necessary to achieve the minimax risk as the oblivious estimator $\hat{\pi}_{HT}$ is already optimal. The dimension reduction scheme $\hat{Q}_{WM}$ allows us to improve the risk bound from that Theorem 2.2 to that of Theorem 2.3. In the specific case where the noise level $\zeta$ is equal to one, the term $n^{2/3}d^{1/3}$ in the risk bound is improved to the optimal one $n^{3/4}d^{1/4}$. Hence, building upon the neighboring groups in $\hat{Q}_{WM}$ turns out to be the key ingredient to recover the minimax risk in the large $d$ regime where $n \in [d^{1/3}, d]$.
<figure><img src="image_4.png" /><figcaption>Figure 5: In these two panels, the group $\bar{P}$ is only made of two types of experts, those in $U^*$ and those in $L^*$. The curves $\bar{m}(V_{r_{cp}}^+)$ and $\bar{m}(V_{r_{cp}}^-)$ respectively correspond to the expected average experts in the neighboring groups $V_{r_{cp}}^+$ and $V_{r_{cp}}^-$ defined in (21). In the left panel, the third and fourth blocks are not selected because the corresponding statistic $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})}$ is small. In the right panel, both the statistic $\hat{C}_{k,r_{cp}}^{(\text{ext})}$ and $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})}$ are large compared to $h$. If $\bar{P}$ is small, this block is selected using a scale $r_{cp} \lesssim \tilde{r}$.</figcaption></figure>

# 4. Partial observations

We now come back to the partial observation setting. Given $\lambda > 0$, we are given $Poi(\lambda nd)$ independent observations $(x_t, y_t)$ where $x_t$ is sampled uniformly in $[n] \times [d]$ and, conditionally to $x_t$, $y_t = M_{x_t} + E_{x_t}$ is an observation of the full model (1) at position $x_t$. As noted above, $\lambda$ stands for the sampling effort and the larger $\lambda$, the more samples on average.

## 4.1. Minimax Lower bound

As in Section 2.1, we first state a minimax lower bound in the case where the noise matrix $E$ is made of independent Gaussian random variables with variance $\zeta^2$. Note that the following minimax lower bound also handle the noise case where $\zeta = 0$, i.e. the noiseless case.

**Theorem 4.1.** There exist universal constants $c$ and $c'$ such that the following holds for any $n \ge 2$, any $d \ge 1$, $\lambda > 0$, and $\zeta \ge 0$:

$$ \inf_{\hat{\pi}} \sup_{M: M_{\pi^{*-1}} \in \mathcal{C}_{BISO}} \mathbb{E}_{(\pi^*, M)} [\|M_{\hat{\pi}^{-1}} - M_{\pi^{*-1}}\|_F^2] \ge c \left[ \left( \mathcal{R}_F[n, d, \zeta/\sqrt{\lambda}] + \frac{n}{\lambda} e^{-2\lambda} \right) \wedge nd \right]. \quad (27) $$

As in the previous minimax lower bound, the quantity $nd$ simply appears because the entries of $M$ lie in $[0, 1]$. In (27), we recognize two terms. First, $\mathcal{R}_F[n, d, \zeta/\sqrt{\lambda}]$ corresponds to the minimax risk for recovering $\pi^*$ in a full observation model with noise $\zeta/\sqrt{\lambda}$. The second term $\frac{n}{\lambda}e^{-2\lambda}$ does not depend on $\zeta$ and is also present in the noiseless setting. It simply quantifies the fact that, for $\lambda < 1$, observations are lacking so that it is impossible to correctly rank experts if there are no observations on the questions on which they are distinct. As the minimax lower bound in (27) turns out to be nearly tight in light of Theorem 4.3, we refer to (27) as the minimax risk in the following. For the purpose of the discussion, we will first focus on the case where $\zeta = 1$ and $\lambda < 1$, which corresponds to the case where we really have partial observations on the matrix. We will then turn to $\zeta = 1$ and $\lambda > 1$, which corresponds to the case where we observe several times each entry of the matrix. Finally, we discuss the noiseless case where $\zeta = 0$.
**Low-sample size.** We first focus on the case where $\zeta = 1$ and $\lambda < 1$, which corresponds to the case where we really have partial observations on the matrix. If $\lambda \le 1/d$, then the minimax risk is of the order of $nd$ and it is impossible to perform significantly better than a random guess. This is not surprising as there are, in expectation, less than one observation on each row. For $\lambda \in [1/d, 1]$, the minimax risk is of the order of

$$ \frac{nd^{1/6}}{\lambda^{5/6}} \wedge \frac{n^{3/4}d^{1/4}}{\lambda^{3/4}} + \frac{n}{\lambda} $$

In the rectangular case where $n \ge d$, the minimax risk is then of the order of $n/\lambda$ for $\lambda \in [1/d, 1]$. When $n \in [d^{1/3}, d]$, the minimax risk is of the order of $n/\lambda$ for $\lambda \in [1/d, n/d]$, and of the order of $\frac{n^{3/4}d^{1/4}}{\lambda^{3/4}}$ for $\lambda \in [n/d, 1]$. For even smaller $n \le d^{1/3}$, there is one more regime since

$$ \mathcal{R}_F[n, d, 1/\sqrt{\lambda}] \asymp \begin{cases} \frac{n}{\lambda} & \text{if } \lambda \in [\frac{1}{d}, \frac{n}{d}] ; \\ \frac{n^{3/4}d^{1/4}}{\lambda^{3/4}} & \text{if } \lambda \in [\frac{n}{d}, \frac{n^3}{d}] ; \\ \frac{nd^{1/6}}{\lambda^{5/6}} & \text{if } \lambda \in [\frac{n^3}{d}, 1] . \end{cases} $$

**Large-sample size.** In the setting where $\lambda > 1$ and $\zeta = 1$, there are several observations per entries. In this case, there are many regimes in (27) that depend on $n$, $d$, $\zeta$, and $\lambda$. To simplify the discussion, we focus here on the case $n = d$ and $\zeta = 1$. Then, the minimax risk is of the order of $\frac{n^{3/4}d^{1/4}}{\lambda^{3/4}}$ for $\lambda \in [1, n^2]$ and is of the order of $n\sqrt{d}/\lambda$ for $\lambda \ge n^2$. This 'easy rate' $n\sqrt{d}/\lambda$ is achieved by the simple procedure that ranks the experts according to the row sums [25, 17]. This simple method turns out to be optimal in the regime where there are more than $n^2$ observations per entry.

**Noiseless case.** In the extreme case where $\zeta = 0$ and $\lambda \ge 1/d$, the minimax risk is of the order of $(n/\lambda) \times e^{-2\lambda}$, which, for some small $\lambda$ is of the order of $n/\lambda$. This minimax lower bound is quite simple to prove. Without loss of generality, suppose that $1/\lambda$ is an integer. Consider a matrix $M$ such that all its columns, except its $1/\lambda$ first ones are constant and equal to one, so that it boils down to considering a ranking problem of size $n \times (1/\lambda)$. In this reduced model, there are two types of experts: (a) experts that are constant and equal to zero and (b) experts that are constant and equal to one. Obviously, if one is given at least one noiseless observation on a row, then it is possible to assign it to a group. However, on each row there is a probability $e^{-1}$ of having no observations. Hence, on expectations there are $n/e$ experts that are impossible to classify. For this reason, any estimator must suffer from a risk at least of the order $n/\lambda$.

## 4.2. Reduction to the full observation model

We now describe a scheme to adapt the estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$ that we developed in the full observation setting of Section 2, to this more general Poissonian setting (2), which encompasses the partial observation setting as well as the over-complete observation setting where each entry is sampled several times. Roughly, if $\lambda$ is small, we simply decrease the number of columns of the matrix $M$ in order to obtain a reduced matrix with full observations. Conversely, if $\lambda$ is really large, which corresponds to the case of multiple observations per entry, we simply average the multiple observations per entry to reduce the noise levels.

As in Section 2, we fix $\delta \in (0, 1)$ that will correspond to a small probability. Given this $\delta > 0$, we denote $\Upsilon^* = \Upsilon^*(n, d, \zeta/\sqrt{\lambda \vee 1})$ the number of independent samples required in Section 2 for the estimation through $\hat{\pi}_{HT}$ or $\hat{\pi}_{WM}$ of the $n \times d$ matrix $M$ with a noise level equal to $\zeta/\sqrt{\lambda \vee 1}$. Recall that $\Upsilon^*$ is of the order of $\log^8(nd(\lambda \vee 1)/(\delta\zeta_-))$.
Define $\lambda_-$ = $\lambda/[4\Upsilon^*]$. For any $i \in [n]$ and any $S \subset [d]$, we write $n_{i,S}$ the number of observations in the sample falling in $\{i\} \times S$, that is $n_{i,S} = |\{t: x_t \in \{i\} \times S\}|$. The following lemma is a simple consequence of Chernoff inequality for Poisson random variables.

**Lemma 4.2.** Assume that $\lambda_- \in [2/d, 1]$, we fix $l(\lambda) = [1/\lambda_-]$. With probability higher than $1 - \delta$, we have

$$ \min_{i \in [n]} \min_{j \in [d/l(\lambda)]} n_{i, [(j-1)l(\lambda)+1, jl(\lambda)]} \geq \Upsilon^* . $$

Now assume that $\lambda_- > 1$. With probability higher than $1 - \delta$, we have

$$ \min_{i \in [n]} \min_{j \in [d]} n_{i, \{j\}} \geq \lambda_- \Upsilon^* . $$

Henceforth, we work under the event introduced in the previous lemma. If this event does not hold, we choose $\hat{\pi}_{WMP}$ arbitrarily. To build $\hat{\pi}_{WMP}$, we consider three subcases that depend on the value of $\lambda$:

1. **Very small sample size.** If $\lambda_- \le 2/d$, then we simply choose $\hat{\pi}_{WMP}$ uniformly at random over the set of all possible permutations. While this choice does not depend on the data and could therefore seem sub-optimal, it is not the case, as the minimax lower bound states that it is impossible to perform better than random guess in this setting.

2. **Small sample size.** If $\lambda_- \in [2/d, 1]$, then we build $\Upsilon^*$ matrices $\mathcal{Y}^\downarrow = (Y^{(0)}, Y^{(1)}, Y^{(\Upsilon^*-1)})$ of size $n \times [d/l(\lambda)]$ in the following way. For any $i \in [n]$, $j \in [[d/l(\lambda)]]$ and $s \in [0, \Upsilon^* - 1]$, $Y_{i,j}^{(s)} = y_t$ where $t$ is the $(s+1)$-th observation such that $x_t \in \{i\} \times [l(\lambda)j+1, l(\lambda)(j+1)]$. On the event of Lemma 4.2, this definition is valid as we observe enough samples for any $i, j$. Then, we compute $\hat{\pi}_{WMP}$ as the variant $\hat{\pi}_{WM-SR}$, introduced in Section H, applied to this sample of reduced matrices.

3. **Large sample size.** If $\lambda_- \ge 1$, then we build $\Upsilon^*$ matrices $\mathcal{Y}^\downarrow = (Y^{(0)}, Y^{(1)}, Y^{(\Upsilon^*-1)})$ of size $n \times d$ in the following way. For any $i \in [n]$, $j \in [d]$, $l \in [[\lambda_-]]$, and $s \in [0, \Upsilon^* - 1]$, define $Y_{i,j}^{(s)} = \frac{1}{[\lambda_-]} \sum_t y_t$ where the $y_t$'s are the $z$-th observations such that $x_t = (i, j)$ with $z \in [1 + (s-1)[\lambda_-], s[\lambda_-]]$. In other words, we build the samples $\mathcal{Y}^\downarrow$ by averaging $[\lambda_-]$ observations on each entries. Again, on the event of Lemma 4.2, this definition is valid as we observed enough samples for any $i, j$. Then, we define $\hat{\pi}_{WMP}$ as $\hat{\pi}_{WM}$ applied to this sample of averaged matrices. By averaging the independent observations, we reduce the noise level of each entry from $\zeta$ to $\zeta/\sqrt{\lambda_-}$.

For $\lambda_- \le 2/d$, there are very few observations on each row so that it is very difficult to compare the experts. For $\lambda_- \in [2/d, 1]$, we have access to less than $\Upsilon^*$ noisy observations of the matrix $M$. The rationale of our procedure is to group together $l(\lambda)$ consecutive questions together in such a way that there are enough observations on each of these groups. The resulting matrices of observations $Y^{(s)}$ have around $\lambda d/\Upsilon^*$ columns. We could have applied the procedure $\hat{\pi}_{WM}$ defined in the previous section to $\mathcal{Y}^\downarrow$, but the corresponding subGaussian norm of the noise would be $1 + \zeta$ (instead of $\zeta$) because there is additional variability coming from the fact that any entry in the reduced matrices has been sampled uniformly among $l(\lambda)$ entries in the original matrices. This would lead us to a procedure achieving the minimax rate with respect to $n, d$, and $\lambda$ but with a suboptimal dependency with respect to $\zeta$ since $\zeta$ would be replaced by $\zeta + 1$. This is the reason why, for $\lambda_- \in [2/d, 1]$, we rely on a slight variant $\hat{\pi}_{WM-SR}$ (see Section H) of $\hat{\pi}_{WM}$ that builds upon the fact that the variations that are due to the aggregation of $M$ are very specific.

**Theorem 4.3.** There exist four numerical constants $c_1 - c_4$ such that the following holds. Fix $\delta = \zeta^2[(\lambda \vee 1)nd]^{-2}$. For any permutation $\pi^* \in \Pi_n$ and any matrix $M$ such that $M_{\pi^{*-1}} \in \mathbb{C}_{BISO}$,
the sorting tree estimator $\hat{\pi}_{WMP}$ defined above satisfies

$$ \mathbb{E} \left[ \| M_{\hat{\pi}_{WMP}^{-1}} - M_{\pi^{*-1}} \|_F^2 \right] \le c_1 \log^{c_2} \left( \frac{nd(\lambda v_1)}{\zeta_-} \right) \left[ \mathcal{R}_F(n, d, \zeta \lambda^{-1/2}) + \frac{n}{\lambda} e^{-c_3 \lambda \log^{-c_4} \left( \frac{nd(\lambda v_1)}{\zeta_-} \right)} \right] . \quad (28) $$

Up to logarithmic terms and up to the logarithmic term inside the exponential term (28), both the minimax upper bound (28) and lower bound (27) match for all values of $n, d, \lambda$ and $\zeta$. As a consequence, this problem of estimating a single permutation $\pi^*$ does not exhibit any significant computational gap.

Let us further discuss and compare the exponential term $n\lambda^{-1}e^{-c_3\lambda\log^{-c_4}(\frac{nd(\lambda v_1)}{\zeta_-})}$ in (27) and $n\lambda^{-1}e^{-2\lambda}$ in (28). First, observe that these two terms are larger than $R_F(n, d, \zeta\lambda^{-1/2})$ only when the noise level $\zeta$ is small, so that it is relevant to discuss them only when $\zeta \ll 1$. Second, note that there is a significant mismatch between these exponential terms only when $\lambda$ is close to one, up to a polylogarithmic factor, since otherwise, either the exponential is close to one (for $\lambda \le 1$) or the exponential is so small that it becomes negligible in comparison to $R_F(n, d, \zeta\lambda^{-1/2})$. One may object that the logarithmic term $\log^{-c_4}(\frac{nd(\lambda v_1)}{\zeta_-})$ may be large in case $\zeta$ is really small –think e.g. of $\zeta = e^{-nd}$. Let us consider this extremely low noise setting where, say $\zeta \le 1/(nd)^2$. If one applies the procedure $\hat{\pi}_{WMP}$ with $\zeta_0 = 1/(nd)^2 \ge \zeta$, then the logarithmic terms become bounded inside the exponential. Since $\mathcal{R}_F(n, d, \zeta_0\lambda^{-1/2})$ is always smaller than $nd \wedge \frac{n}{\lambda}e^{-c_3\lambda}$ provided that $\lambda \le 1$, this estimator achieves the risk bound $\frac{n}{\lambda} \wedge nd$, which is optimal for all $\zeta \in [0, 1/(nd)^2]$ and all $\lambda \le 1$. To sum up, there is gap between our minimax lower and upper bounds only either (i) in the low-noise level with large but mild sampling effort, that is $\zeta = o(\log^{-c}(nd))$, $\zeta \ge (nd)^{-2}$, and $\lambda \in [\log^c \log(nd), \log^{c'}(nd)]$ for some $c$ and $c' > 0$ or (ii) in the extremely low noise level with large sampling effort, that is $\zeta \le (nd)^{-2}$ and $\lambda \ge 1$.

In $\hat{\pi}_{WMP}$, we have plugged in the hierarchical sorting tree estimator with memory $\hat{\pi}_{WM}$. If we had plugged in the oblivious hierarchical sorting tree estimator $\hat{\pi}_{HT}$, then the resulting estimator would satisfy a similar rate similar to (28) except that the term $n^{3/4}d^{1/4}/\lambda^{3/4}$ would be replaced by the slower rate $n^{2/3}d^{1/3}/\lambda^{2/3}$.

### 4.3. Reconstruction of the matrix M

In this subsection, we assume again that the noise level $\zeta = 1$ to simplify the exposition. As alluded in Section 2, it is quite straightforward to estimate the matrix $M$ and control the corresponding loss $\|\hat{M}-M\|_F^2$ by a simple subsampling step explained e.g. in [17] that we recall here. First, we split the sample into two part by assigning independently each observation to the first subsample with probability $1/2$ and the second subsample with probability $1/2$. Then, we use the first subsample to estimate the permutation $\hat{\pi}$ of the experts. As for the second subsample $(x_t^{(2)}, y_t^{(2)})$, we define the empirical observed matrix $Y^{(2)}$ by

$$ Y_{i,j}^{(2)} = \frac{1}{\lambda} \sum_t y_t^{(2)} \mathbf{1}_{x_t^{(2)}=(i,j)}. $$

Then, we compute the least-square estimator $\hat{M}_{\hat{\pi}}$ of $M_{\hat{\pi}}$ in the class of bi-isotonic matrix $\hat{M}_{\hat{\pi}} = \arg\min_{B \in \mathcal{B}_{\text{BISO}}} \|B - Y_{\hat{\pi}}^{(2)}\|_F^2$. This estimator can be computed in near linear-time [13]. Then, Proposition 3.3 in [17] states, that with high probability, the loss $\|\hat{M}-M\|_F^2$ is, up to logarithmic terms, smaller than the sum of the minimax risk for estimating a bi-isotonic matrix $B$ and the loss $\|M_{\hat{\pi}^{-1}} - M_{\pi^{*-1}}\|_F^2$. Plugging this proposition with our estimator $\hat{\pi}_{WMP}$ with
$\delta = (\lambda \vee 1)/(np)$, we readily arrive to the following risk bound for the corresponding estimator $\widehat{M}_{WMP}$.

Define $\mathcal{R}_1(n, d, \lambda) = \sqrt{\frac{nd}{\lambda}} \wedge \frac{nd}{\lambda^{2/3}(nvd)^{2/3}} \wedge \frac{nd}{\lambda}$. Mao et al. [17] have proved that, up to polylogarithmic factor and up to a possible additive term $(n \wedge d)/\lambda$, the minimax risk in square Frobenius norm for estimating a bi-isotonic matrix with partial observations is $\mathcal{R}_1(n, d, \lambda)$.

**Corollary 4.4.** There exist two numerical constants $c$ and $c'$ such that the following holds. For any permutation $\pi^* \in \Pi_n$ and any matrix $M$ such that $M_{\pi^*} \in \mathbb{C}_{BISO}$, we have

$$ \begin{aligned} \mathbb{E}[\|\widehat{M}_{WMP} - M\|_F^2] &\le (nd) \wedge \left[c \log^c((\lambda \vee 1)nd) (\mathcal{R}_1(n, d, \lambda) + \mathcal{R}_F(n, d, \lambda^{-1/2}))\right] \\ &\le (nd) \wedge \left[c \log^c((\lambda \vee 1)nd) (\mathcal{R}_1(n, d, \lambda) + \frac{n}{\lambda})\right]. \end{aligned} \quad (29) $$

The proof is a straightforward consequence of Proposition 3.3 in [17] and Theorem 4.3 and is therefore omitted. It turns out that $\mathcal{R}_F(n, d, \lambda^{-1/2})$ is always smaller than $\mathcal{R}_1(n, d, \lambda) + \frac{n}{\lambda}$, so that the cost of reconstruction for not knowing $\pi^*$ is $n/\lambda$.

This risk bound (29) is minimax optimal, up to polylogarithms, and this for all possible values of $n \ge 2$, $d$, and $\lambda > 0$. Indeed, in their Theorem 3.1, Mao et al. [17] provide a matching minimax lower bound in $\mathcal{R}_1(n, d, \lambda)$ in the specific case where $n \ge d$, but their proof easily extends to the case where $n \le d$. Besides, our proof of the minimax lower bound $\frac{n}{\lambda}$ in Theorem 4.1 for the problem of estimating $\pi^*$ straightforwardly extends to the problem of matrix estimation (recall that we consider $\zeta = 1$ here).

The least-square estimator $\hat{\pi}_{LS}$ of Mao et al. has also been proved to achieve the minimax risk for $n \ge d$ –see their theorem 3.1 in [17]. However, no efficient algorithm is known for computing this estimator in $\hat{\pi}_{LS}$, so that our estimator $\widehat{M}_{WMP}$ is, to the best of our knowledge, the first efficient minimax-optimal estimator for estimating $M$ in this context, for any values of $n, d, \lambda$.

### 4.4. Bounds for the max loss of Mao et al. [17]

In [17], Mao et al. control, for an estimator $\hat{\pi}$ of the permutation, a different loss from ours. Up to normalization factors, they indeed focus on the maximum $l_2$ norm of the rows of $(M_{\hat{\pi}^{-1}})_{i, \cdot} - (M_{\pi^{*-1}})_{i, \cdot}$, that is

$$ l_{\infty}(\hat{\pi}, \pi^*) = \sup_{i \in [n]} \|(M_{\hat{\pi}^{-1}})_{i, \cdot} - (M_{\pi^{*-1}})_{i, \cdot}\|_2^2. \quad (30) $$

This loss also considered in [26, 8] corresponds to some maximum error of the estimated permutation so that $l_{\infty}(\hat{\pi}, \pi^*) \ge \|M_{\hat{\pi}^{-1}} - M_{\pi^{*-1}}\|_F^2/n$. Alternatively, we can define the loss $l_{err}$

$$ l_{err}(\hat{\pi}, \pi^*) = \max_{\hat{\pi}(i) < \hat{\pi}(j) \text{ and } \pi^*(i) > \pi^*(j)} \|M_{i, \cdot} - M_{j, \cdot}\|_2^2, $$

which quantifies the maximum distance between two experts that have not been ranked in a consistent manner. The loss $l_{\infty}$ and $l_{err}$ turn out to be equivalent as stated in the following lemma.

**Lemma 4.5.** For any permutation $\hat{\pi}$, we have

$$ l_{\infty}(\hat{\pi}, \pi^*) \le l_{err}(\hat{\pi}, \pi^*) \le 4l_{\infty}(\hat{\pi}, \pi^*) \quad (31) $$

To simplify the discussion in this section, we assume again that the noise level $\zeta$ equals one. Mao et al. [17] provide a simple polynomial time $\hat{\pi}_{ref}$ achieving

$$ \mathbb{E}[l_{\infty}(\hat{\pi}_{\text{ref}}, \pi^*)] \lesssim d \wedge \frac{d^{1/4}}{\lambda^{3/4}} \log^{3/4}(n). \quad (32) $$
Conversely, they prove in their Theorem 3.7 that any estimator $\hat{\pi}$ that only ranks the experts i and j according to the differences of the observations on the rows i and j must incur this risk bound- see [17] for further details. Besides, they conjecture that the risk bound (32) cannot be improved. In [14], Liu and Moitra already pointed out that the max loss $l_{\infty}(\hat{\pi}, \pi^*)$ is less suited than the loss $\|M_{\hat{\pi}} - M_{\pi^*}\|_F^2$ for the purpose of estimating the matrix M -see the discussion in the previous subsection. Still, controlling the max loss $l_{\infty}(\hat{\pi}, \pi^*)$ may be an objective per se, and the study of its minimax value and of the existence of related minimax estimators is relevant. In the following proposition, which is mainly a consequence of our results and proof techniques, we disprove Mao et al.'s conjecture by introducing an estimator $\hat{\pi}_{PC}$ achieving a faster rate than (32). Besides, this rate turns out to be minimax-optimal.

**Proposition 4.6.** There exist numerical constants c, c', and c'' such that the following result holds. There exists a polynomial-time estimator $\hat{\pi}_{PC}$ that performs pair-wise comparisons between the experts and that achieves the risk bound

$$ \mathbb{E}[l_{\infty}(\hat{\pi}_{PC}, \pi^*)] \le c \log^{c'} (nd(\lambda \lor 1)) \left[ \frac{d^{1/6}}{\lambda^{5/6}} \wedge \frac{\sqrt{d}}{\lambda} \right] \wedge d. \quad (33) $$

Conversely, for any $n \ge 2$, any $d \ge 1$, and $\lambda > 0$, we have

$$ \inf_{\hat{\pi}} \sup_{\pi^* \in \Pi_n} \sup_{M: M_{\pi^{*-1}} \in \mathbb{C}_{BISO}} \mathbb{E}_{(\pi^*, M)}[l_{\infty}(\hat{\pi}, \pi^*)] \ge c'' \left[ \frac{d^{1/6}}{\lambda^{5/6}} \wedge \frac{\sqrt{d}}{\lambda} \wedge d \right]. \quad (34) $$

For $\lambda \le 1/d$, it is not possible to perform significantly better than random guess. Then, in the interesting regime $\lambda \in [1/d, d^2]$, the risk is of the order of $d^{1/6}/\lambda^{5/6}$. It turns out that this rate corresponds, up to polylogarithmic terms, to the minimal distance between two experts so that one is able to consistently compare them. For very large sample size $\lambda \ge d^2$, we arrive at the easy regime which is of the order of $\sqrt{d}/\lambda$.

The estimator $\hat{\pi}_{PC}$ is based on pairwise comparisons. For any two experts i and j, we apply the procedure $\hat{\pi}_{WMP}$ to i and j with $\delta = [(\lambda \lor 1)(n^2d)]^{-2}$. If the trisection $(O, P, I)$ is of the form $(\emptyset, \{i\}, \{j\})$, we return $i < j$. If the trisection $(O, P, I)$ is of the form $(\emptyset, \{j\}, \{i\})$, we return $j < i$. Otherwise, we return nothing. Applying this comparison algorithm to all $(i, j)$, we recover a set of pairwise comparisons $\mathcal{P}C = \{(i, j) : i < j\}$. With high probability -see the proof for more details-, it turns that $\mathcal{P}C$ satisfies two properties:

(i) $\mathcal{P}C$ is consistent. For any $(i, j) \in \mathcal{P}C$, we have $\pi^*(i) < \pi^*(j)$.

(ii) $\mathcal{P}C$ contains all 2-tuple of experts that are far apart. More precisely, $\mathcal{P}C$ contains all $(i, j)$ such that $\pi^*(i) < \pi^*(j)$, and

$$ \|M_{i,\cdot} - M_{j,\cdot}\|_2^2 \ge c \log^{c'} (nd(\lambda \lor 1)) \left[ \frac{d^{1/6}}{\lambda^{5/6}} \wedge \frac{\sqrt{d}}{\lambda} \right] \wedge d, \quad (35) $$

for suitable constants c and c'.

Then, define the function $\phi: [n] \mapsto \mathbb{N}$ by $\phi(j) = |\{(i, j) : (i, j) \in \mathcal{P}C\}|$ which simply counts the number of experts i that are detected to be lower than j. Finally, we build $\hat{\pi}_{PC}$ as any permutation that ranks the experts consistently with $\phi$.

In fact, the procedure for computing $\hat{\pi}_{PC}$ could be greatly simplified. Indeed, as we only perform pairwise comparisons, some parts of **BlockSort** turn out to be irrelevant. For instance, the PCA steps are not required. Besides, the sample splits could be avoided and it could even be possible to work with a single observation. As the problem of optimal permutation recovery with respect to the $l_{\infty}$ loss is not the main scope of this paper, we do not provide a simplified and dedicated algorithm. Besides, we conjecture that our original estimator $\hat{\pi}_{WMP}$ also achieves the minimax risk (33) with respect to the $l_{\infty}$ loss.
**Acknowledgements.** The work of A. Carpentier is partially supported by the Deutsche Forschungsgemeinschaft (DFG) Emmy Noether grant MuSyAD (CA 1488/1-1), by the DFG - 314838170, GRK 2297 MathCoRe, by the FG DFG, by the DFG CRC 1294 'Data Assimilation', Project A03, by the Forschungsgruppe FOR 5381 "Mathematical Statistics in the Information Age - Statistical Efficiency and Computational Tractability", Project TP 02, by the Agence Nationale de la Recherche (ANR) and the DFG on the French-German PRCI ANR ASCAI CA 1488/4-1 "Aktive und Batch-Segmentierung, Clustering und Seriation: Grundlagen der KI" and by the UFA-DFH through the French-German Doktorandenkolleg CDFA 01-18 and by the SFI Sachsen-Anhalt for the project RE-BCI. The work of E. Pilliat and N. Verzelen has been partially supported by ANR-21-CE23-0035 (ASCAI).

## Appendix A: Full description of the procedures

In this section, we provide a fuller description of the estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$ as a collection of algorithms. We will rely on this description in the analysis of these estimators. To ease its understanding, we make this section completely self-contained. As a consequence, the material presented here is partly redundant with Section 3.

### A.1. Sorting a group of experts

Some of the notation have already been introduced in Section 3. Still we define them again here for the sake of completeness. We write $\mathcal{D}$ for the set of all dyadic numbers, that is $\mathcal{D} = \{2^k : k \in \mathbb{Z}\}$. Equipped with $\mathcal{D}$, let

$$ \mathcal{R} = \mathcal{D} \cap [1, d] \quad \text{and} \quad \mathcal{H} = \mathcal{D} \cap \left[ \frac{\zeta^2}{nd}, 1 \right], \qquad (36) $$

respectively denote the dyadic collection of numbers between 1 and $d$ and the dyadic collection of numbers between $1/nd$ and 1.

Besides for an integer $r \in \mathcal{R}$, we write $Q_r$ for the regular grid of $[d]$ of width $r$:

$$ Q_r = \{1, r+1, 2r+1, \dots \left\lfloor \frac{d}{r} \right\rfloor r+1\}. $$

In contrast to Section 3, we start by describing the simple comparison routine before moving to the dimension reduction techniques and to the general architecture of the procedures.

Given a collection $\overline{P}$ of experts, some data $Z \in \mathbb{R}^{P \times Q}$ and a direction $w \in (\mathbb{R}^+)^Q$ and a pivot $\gamma \in [1 : |\overline{P}|]$, the following pivoting algorithm sorts the experts in $\overline{P}$ according to the projection of the data onto the vector $w$. More precisely, it returns four subsets $\overline{L} \subset L$ and $\overline{U} \subset U$ of experts such that the $\gamma$-th best expert according to the $(Z, w)$-order - as defined above Equation (15) - is significantly above all experts in $L$ and below all experts in $U$. The subsets $L$ and $\overline{L}$ (resp. $U$ and $\overline{U}$) differ in the level of significance we require. We define the tuning parameters $\beta_{\text{tris}}$ and $\bar{\beta}_{\text{tris}}$ for **Pivot**

$$ \beta_{\text{tris}} = 4\sqrt{2}\zeta, \quad \bar{\beta}_{\text{tris}} = 8\sqrt{2}\zeta. \qquad (37) $$
### Algorithm 4 Pivot(Z, w, γ)

**Require:** A matrix $Z \in \mathbb{R}^{\bar{P} \times Q}$ with $\bar{P} \subset [n]$ a set of experts and $Q \subset [d]$ a set of blocks, a direction $w \in \mathbb{R}_+^Q$, $w \neq 0$ and a pivot index $\gamma$

**Ensure:** Two couples of subsets $(L, U)$ and $(\bar{L}, \bar{U})$ of $\bar{P}$

1: **for** $i \in \bar{P}$ **do**

2: Compute the statistic $\psi(i, w) = (Z_{i, \cdot}, \frac{w}{\|w\|_2})$

3: **end for**

4: Sort the statistics $\psi(i, w)$: $\psi(i_1, w) \le \dots \le \psi(i_{|\bar{P}|}, w)$

5: $U = \{i \in \bar{P} : \psi(i, w) > \psi(i_{\gamma}, w) + \beta_{\text{tris}} \sqrt{\log(2|\bar{P}|/\delta)}\}$

6: $\bar{U} = \{i \in \bar{P} : \psi(i, w) > \psi(i_{\gamma}, w) + \bar{\beta}_{\text{tris}} \sqrt{\log(2|\bar{P}|/\delta)}\}$

7: $L = \{i \in \bar{P} : \psi(i, w) < \psi(i_{\gamma}, w) - \beta_{\text{tris}} \sqrt{\log(2|\bar{P}|/\delta)}\}$

8: $\bar{L} = \{i \in \bar{P} : \psi(i, w) < \psi(i_{\gamma}, w) - \bar{\beta}_{\text{tris}} \sqrt{\log(2|\bar{P}|/\delta)}\}$

9: **return** $(L, U), (\bar{L}, \bar{U})$

When the vector $w$ is equal to $\mathbf{1}_Q$, we simply write $\text{Pivot}(Z, \gamma)$ instead of $\text{Pivot}(Z, \mathbf{1}_Q, \gamma)$ for the sake of simplicity.

$$ \text{Pivot}(Z, \gamma) = \text{Pivot}(Z, w = \mathbf{1}_Q, \gamma). \qquad (38) $$

In fact, $\text{Pivot}(Z, \gamma)$ simply amounts to comparing the row sums of $Z$ for each of the experts in $P$.

In the next two pages, we redefine in more detail the Double Trisection algorithm of Section 3. First, **DoubleTrisection – PCA** relies on a PCA-type argument to find a suitable direction $\hat{w}^+$ and then provides two trisections of the subset $\bar{P}$ of experts using the **Pivot** sub-routine.

### Algorithm 5 DoubleTrisection – PCA($\mathcal{Z}$, $\gamma$)

**Require:** 4 reduced samples $\mathcal{Z} = (Z^{(1)}, Z^{(2)}, Z^{(3)}, Z^{(4)})$ where $Z^{(1)}, Z^{(2)}, Z^{(3)} \in \mathbb{R}^{\bar{P} \times Q}$ and $Z^{(4)} \in \mathbb{R}^{\bar{P} \times Q}$ with some $\bar{P} \subset \bar{P}$, and a pivot index $\gamma$

**Ensure:** Four subsets $(L_{\text{pca}}, U_{\text{pca}})$ and $(\bar{L}_{\text{pca}}, \bar{U}_{\text{pca}})$ of $\bar{P}$

1: Compute the following vector with coefficients in $\bar{P}$:

$$ \hat{v} = \arg \max_{\|v\| \le 1} \left[ \|v^T(Z^{(1)} - \bar{Z}^{(1)})\|_2^2 - \frac{1}{2} \|v^T(Z^{(1)} - \bar{Z}^{(1)} - Z^{(2)} + \bar{Z}^{(2)})\|_2^2 \right] \in \mathbb{R}^{\bar{P}} $$

2: $\hat{z} = v^T Z^{(3)} \in \mathbb{R}^Q$

3: Define $\hat{w}^+$ by $(\hat{w}^+)_{l} = |\hat{z}_l| \mathbf{1}_{|\hat{z}_l| \ge 2\zeta\sqrt{2\log(2|Q|/\delta)}}$

4: $(L_{\text{pca}}, U_{\text{pca}}), (\bar{L}_{\text{pca}}, \bar{U}_{\text{pca}}) = \text{Pivot}(Z^{(4)}, \hat{w}^+, \gamma)$

5: **return** $(L_{\text{pca}}, U_{\text{pca}}), (\bar{L}_{\text{pca}}, \bar{U}_{\text{pca}})$

Next, **DoubleTrisection – Local($\mathcal{Z}$, $\gamma$)** builds two trisections of $\bar{P}$ based on the reduced samples. First, it builds these trisections by simply using the row sums on the data and then it improves them thanks to **DoubleTrisection – PCA**.
**Algorithm 6 DoubleTrisection – Local($\mathcal{Z}$, $\gamma$)**

**Require:** 5 reduced samples $\mathcal{Z} = (Z^{(1)}, Z^{(2)}, Z^{(3)}, Z^{(4)}, Z^{(5)})$ in $\mathbb{R}^{\bar{P} \times Q}$, a pivot index $\gamma$ and a threshold $\beta_{\text{tris}}$

**Ensure:** Two couples of subsets $(L, U)$ and $(\bar{L}, \bar{U})$ of $\bar{P}$

1: $(L_{\text{cp}}, U_{\text{cp}}), (\bar{L}_{\text{cp}}, \bar{U}_{\text{cp}}) = \textbf{Pivot}(Z^{(1)}, 1_Q, \gamma)$

2: Set $\tilde{P} = \bar{P} \setminus (\bar{L}_{\text{cp}} \cup \bar{U}_{\text{cp}})$

3: Set $\mathcal{Z}' = (Z^{(2)}(\tilde{P}), Z^{(3)}(\tilde{P}), Z^{(4)}(\tilde{P}), Z^{(5)}(\bar{P}))$ be the sequence of reduced samples where the three first samples are restricted to $\tilde{P}$.

4: $(L_{\text{pca}}, U_{\text{pca}}), (\bar{L}_{\text{pca}}, \bar{U}_{\text{pca}}) = \textbf{DoubleTrisection} - \textbf{PCA}(\mathcal{Z}', \gamma)$

5: Set $L = L_{\text{cp}} \cup L_{\text{pca}}$ and $\bar{L} = \bar{L}_{\text{cp}} \cup \bar{L}_{\text{pca}}$ and $U = U_{\text{cp}} \cup U_{\text{pca}}$ and $\bar{U} = \bar{U}_{\text{cp}} \cup \bar{U}_{\text{pca}}$

6: **return** $(L, U), (\bar{L}, \bar{U})$

To finish defining **DoubleTrisection**, we simply need to plug a dimension reduction procedure to select a subset of questions $Q \subset [d]$ and then to sum the data on these questions.

The two following algorithms are mainly definitions. For some data $Y \in \mathbb{R}^{[n] \times [d]}$, a set of experts $\bar{P}$ and a set of blocks $Q \subset Q_r$, and a scale $r$, the $|\bar{P}| \times |Q|$ matrix **Encode – Matrix**($Y, \bar{P}, Q, r$) is simply a reduced matrix where we consider the normalized row sums of $Y$ around the questions of $Q$ at scale $r$.

**Algorithm 7 Encode – Matrix($Y, \bar{P}, Q, r$)**

**Require:** A matrix $Y \in \mathbb{R}^{[n] \times [d]}$, a set of experts $\bar{P}$ and a set of blocks $Q \subset Q_r$, a scale $r$

**Ensure:** A reduced matrix $Z \in \mathbb{R}^{\bar{P} \times Q}$

1: **for** $i \in \bar{P}$ and $l \in Q$ **do**

2: Define $Z_{i,l} = \frac{1}{\sqrt{r}} \sum_{k \in [l, l+r)} Y_{i,k}$

▷ $Y_{i,k} = 1$ for $k \ge d+1$

3: **end for**

4: **return** $Z \in \mathbb{R}^{\bar{P} \times Q}$

▷ the restriction of $Z$ to $\bar{P}$ and $Q$

Second, **Encode – Set($D, r$)** transforms a subset $[d]$ of questions into a subset $Q \subset Q_r$ of blocks of questions at scale $r$.

**Algorithm 8 Encode – Set($D, r$)**

**Require:** A set of questions $D \subset [d]$, a scale $r \in \mathcal{R}$

**Ensure:** A set of blocks $Q \subset Q_r$

return $Q = \{l \in Q_r : [l, l+r) \cap D \neq \emptyset\}$

Then, we are in position to redefine this version of Algorithm 3. As in the original definition in Section 3, there are two variations of this procedure depending on whether we are building the estimator $\hat{\pi}_{HT}$ or the estimator $\hat{\pi}_{WM}$ that uses the memory of the tree. Algorithm **DoubleTrisection**($\mathcal{Y}, \mathcal{T}, \bar{P}, \gamma$) takes some original data and then reduces the dimension of the problem to build two trisections of the set $\bar{P}$ of experts.
**Algorithm 9 DoubleTrisection(\mathcal{Y}, \mathcal{T}, \bar{P}, \gamma)**

**Require:** 6 samples $\mathcal{Y} = (Y^{(1)}, \dots, Y^{(6)})$, a tree $\mathcal{T}$, a set of expert $\bar{P}$ included in a leaf $G$ of $\mathcal{T}$ at maximal depth and a pivot index $\gamma$

**Ensure:** Two couples of subsets $(L, U)$, $(\bar{L}, \bar{U})$ of $\bar{P}$

1: Initialize $L, U, \bar{L}, \bar{U} = \emptyset$

2: **for** $h \in \mathcal{H}, r \in \mathcal{R}$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**if** Not using the memory of the tree **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\tilde{Q} := \tilde{Q}_{cp}(h, r) = \text{DimensionReduction}(Y^{(1)}, \bar{P}, h, r)$ - see Algorithm 12 or (20)

5: &nbsp;&nbsp;&nbsp;&nbsp;**else if** Using the memory of the tree **then**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\tilde{Q} := \tilde{Q}_{WM}(h, r) = \text{DimensionReduction} - \text{WM}(Y^{(1)}, \mathcal{T}, \bar{P}, h, r)$ - See Algorithm 13 or (26)

7: **end if**

8: &nbsp;&nbsp;&nbsp;&nbsp;Consider the five samples $\mathcal{Y}' = (Y^{(2)}, Y^{(3)}, Y^{(4)}, Y^{(5)}, Y^{(6)})$

9: &nbsp;&nbsp;&nbsp;&nbsp;Consider the five reduced samples $\mathcal{Z} = \text{Encode} - \text{Matrix}(\mathcal{Y}', \bar{P}, \tilde{Q}_{WM}, r)$

10: Compute $(L_{\text{loc}}, U_{\text{loc}}), (\bar{L}_{\text{loc}}, \bar{U}_{\text{loc}}) = \text{DoubleTrisection} - \text{Local}(\mathcal{Z}, \gamma)$

11: Update $L = L \cup L_{\text{loc}}$ and $\bar{L} = \bar{L} \cup \bar{L}_{\text{loc}}$ and $U = U \cup U_{\text{loc}}$ and $\bar{U} = \bar{U} \cup \bar{U}_{\text{loc}}$

12: **end for**

13: **return** $(L, U), (\bar{L}, \bar{U})$

Finally, we reproduce **BlockSort** here that was originally defined in Algorithm 2. We recall that **BlockSort** iteratively applies a logarithmic number of times the procedure **DoubleTrisection** to build two suitable trisections of a set $G$ of experts. Although implicit in this description, there are two different versions of the corresponding procedure whether we use the memory of the tree - estimator $\hat{\pi}_{WM}$ - or not - estimator $\hat{\pi}_{HT}$ in **DoubleTrisection**. In the following, $\tau_{\infty} = [4 \cdot 10^7 \log^7(\frac{nd}{\delta(\zeta^{-})^2})]$ stands for the number of iterations in **BlockSort**.

**Algorithm 10 BlockSort(\mathcal{Y}, \mathcal{T}, G)**

**Require:** 6$\tau_{\infty}$ samples $\mathcal{Y} = (Y^{(0)}, \dots, Y^{(6\tau_{\infty}-1)})$, the tree $\mathcal{T}$, a leaf $G \subset [n]$ in $\mathcal{T}$ at maximal depth

**Ensure:** A partition of $G$ into three groups $(O, P, I)$

1: Set $\gamma = \lfloor|G|/2\rfloor$ and $O_0, I_0, \bar{O}_0, \bar{I}_0 = \emptyset$

2: **for** $\tau = 0, \dots, \tau_{\infty} - 1$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;Consider 6 fresh samples $\mathcal{Y}_{\tau} = (Y^{(6\tau)}, \dots, Y^{(6\tau+5)})$

4: &nbsp;&nbsp;&nbsp;&nbsp;set $\gamma = \lfloor|G|/2\rfloor - |\bar{O}_{\tau}|$

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$(L_{\tau}, U_{\tau}), (\bar{L}_{\tau}, \bar{U}_{\tau}) = \text{DoubleTrisection}(\mathcal{Y}_{\tau}, \mathcal{T}, G \setminus (\bar{O}_{\tau} \cup \bar{I}_{\tau}), \gamma)$ as in Algorithm 9

6: &nbsp;&nbsp;&nbsp;&nbsp;Update $O_{\tau+1} = O_{\tau} \cup L_{\tau}, I_{\tau+1} = I_{\tau} \cup U_{\tau}, \bar{O}_{\tau+1} = \bar{O}_{\tau} \cup \bar{L}_{\tau}, \bar{I}_{\tau+1} = \bar{I}_{\tau} \cup \bar{U}_{\tau}$

7: **end for**

8: **if** $O_{\tau_{\infty}} \cap I_{\tau_{\infty}} \neq \emptyset$ **then**

9: &nbsp;&nbsp;&nbsp;&nbsp;Set $O_{\tau_{\infty}} := O_{\tau_{\infty}} \setminus I_{\tau_{\infty}}$ and $I_{\tau_{\infty}} := I_{\tau_{\infty}} \setminus O_{\tau_{\infty}}$

10: **end if**

11: **return** $(O_{\tau_{\infty}}, G \setminus (O_{\tau_{\infty}} \cup I_{\tau_{\infty}}), I_{\tau_{\infty}})$

Under an event of high probability (to be later discussed), we have $O_{\tau_{\infty}} \cap I_{\tau_{\infty}} = \emptyset$. The correction at the end of the algorithm simply forces the algorithm to return a partition of $G$.

## A.2. Hierarchical sorting Trees and TreeSort Algorithm

In this subsection, we formally describe how we build and navigate into a hierarchical tree. In the following, a node $G \in \text{Nodes}$ is a labelled subset of $[n]$. Its label belongs to $\{0, p, 1\}$. For a node $G$, we write **Type**(G) for the label (also called type) of G.

**Definition 1.** (Hierarchical sorting Trees) A hierarchical sorting tree $\mathcal{T}$ is a rooted tree that satisfies the three following properties:

* The root $G$ of $\mathcal{T}$ corresponds to the set $[n]$ and its label is 0.
* Any node $G$ of type **p** is a leaf.

* Any node $G$ of type in $\{0, 1\}$ is either a leaf or has three children $(O, P, I)$ with type **0**, **p**, **1** respectively. Besides, $(O, P, I)$ correspond to a partition of $G$.

We write $\mathcal{T}_0$ for the tree of depth 0. The procedure **TreeSort** iteratively builds a hierarchical sorting tree. Hence, we need to define the operation of adding children to a leaf in a tree $\mathcal{T}$. For a specific leaf $G$ of type **0** or **1**, we consider three labelled subsets $O, P, I$ of type **0**, **p**, **1**, respectively. Besides, those subsets satisfy the third condition in Definition 1. Then, $\mathcal{T}' = \text{AddChild}(\mathcal{T}, G, (O, P, I))$ is the supertree of $\mathcal{T}$ where we have added the nodes $(O, P, I)$ as children of $G$. Finally, we observe that for any $t > 0$, all the nodes at depth $t$ of a hierarchical sorting tree $\mathcal{T}$ are disjoint.

In fact, we shall prove in Proposition B.1 and in Corollary B.4 that, with high probability, the final tree $\mathcal{T}_{t_\infty}$ turns out to be a valid hierarchical sorting tree as defined below.

**Definition 2.** (Valid hierarchical sorting Tree) A hierarchical sorting tree $\mathcal{T}$ is valid if non-terminal nodes $G$ of $\mathcal{T}$ satisfy the two following additional properties: if we denote $(O, P, I)$ their children of type **0**, **p**, **1** respectively, then

* All the experts in $O$ are below those of $I$. In other words, for any $i \in O$ and any $j \in I$, we have $\pi^*(i) < \pi^*(j)$.

* $|O| < |G|$ and $|I| < |G|$.

The second property ($|O| < |G|$ and $|I| < |G|$) forces the tree to be finite.

For a node $G$ in a such valid hierarchical sorting tree $\mathcal{T}$, $\text{Depth}(\mathcal{T}, G)$ stands for the depth of $G$ in $\mathcal{T}$. In light of this definition of valid hierarchical sorting trees, a labelled subset $G$ cannot appear twice in a tree $\mathcal{T}$, so that $\text{Depth}(\mathcal{T}, G)$ is well-defined.

We are now equipped to provide a more formal definition of **TreeSort**, although the procedure is in fact the same as the one described in Algorithm 1. Let $t_\infty = \lceil \log(n)/\log(2) \rceil$.

<caption>**Algorithm 11 TreeSort($\mathcal{Y}$)**</caption>

**Require:** $6\tau_\infty t_\infty$ samples $\mathcal{Y} = (Y^{(0)}, \dots, Y^{(6\tau_\infty t_\infty - 1)})$
**Ensure:** A final tree $\mathcal{T}$
1: $\mathcal{T} = \mathcal{T}_0$ \quad $\triangleright$ The root is at depth 0
2: for $t = 0, \dots, t_\infty - 1$ do
3: &nbsp;&nbsp;&nbsp;&nbsp;Consider $6\tau_\infty$ fresh samples $\mathcal{Y} = (Y^{(6t\tau_\infty)}, \dots, Y^{(6(t+1)\tau_\infty - 1)})$
4: &nbsp;&nbsp;&nbsp;&nbsp;for $G \in \mathcal{L}^{(0,1)}(\mathcal{T})$ do
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$(O_G, P_G, I_G) = \text{BlockSort}(\mathcal{Y}, \mathcal{T}, G)$
6: &nbsp;&nbsp;&nbsp;&nbsp;Set $\text{Type}(O_G) = \mathbf{0}$ and $\text{Type}(P_G) = \mathbf{p}$ and $\text{Type}(I_G) = \mathbf{1}$
7: &nbsp;&nbsp;&nbsp;&nbsp;end for
8: for $G \in \mathcal{L}^{(0,1)}(\mathcal{T})$ do
9: &nbsp;&nbsp;&nbsp;&nbsp;\text{AddChild}(\mathcal{T}, G, (O_G, P_G, I_G))
10: end for
11: end for
12: return $\mathcal{T}$

As explained in Section 3, the final estimators $\hat{\pi}_{HT}$ or $\hat{\pi}_{WM}$ are computed from their corresponding hierarchical sorting tree $\mathcal{T}$.

In order to define the **DimensionReduction** – **WM** algorithm in the next subsection, we need to introduce a few more notation. First, we define

$$ \mathcal{L}^{(0,1)}(\mathcal{T}) = \{G \in \text{Leaves}(\mathcal{T}) : \text{Type}(G) \in \{0, 1\}\}, \quad (39) $$

as the collection of leaves of $\mathcal{T}$ that are either of type **0** or of type **1**. In the algorithm **TreeSort**, these are the leaves to be partitionned. In particular at step $t$ of **TreeSort**, $\mathcal{L}^{(0,1)}(\mathcal{T})$ is only made of leaves at depth $t$.
For a subset $P \subset [n]$, **Leaf**($\mathcal{T}$, $P$) is defined as the leaf $G \in \text{Leaves}(\mathcal{T})$ containing $P$ (if it exists). Finally, the groups $G \in \mathcal{L}^{(0,1)}(\mathcal{T})$ inherit from a natural order provided that $\mathcal{T}$ is a valid hierarchical sorting tree. We can enumerate the groups $G_1, G_2, \dots, G_{|\mathcal{L}^{(0,1)}(\mathcal{T})|}$ in such a way that all the experts in $G_s$ are below those of $G_{s'}$ for $s < s'$. To ease the presentation, we also introduce, for any positive integer $s$ the groups $G_{|\mathcal{L}^{(0,1)}(\mathcal{T})|+s} = \{n+s\}$. The corresponding data and signal for the $n+s$-th expert satisfies $Y_{n+s,j} = 1 = M_{n+s,j}$ almost-surely for any $j \in [d]$. Also, for any positive integer $s$ we introduce the groups $G_{1-s} = \{1-s\}$. The corresponding data and signal for this synthetic expert satisfy $Y_{1-s,j} = 0 = M_{1-s,j} = 0$ almost-surely for any $j \in [d]$.

Then, for a specific leaf $G_s \in \mathcal{L}^{(0,1)}(\mathcal{T})$, **Order**($\mathcal{T}$, $G$) stands for the collection $(G^{(a)})$, $a \in \mathbb{Z}$ of leaves where $G^{(a)} = G_{a+s}$. In other words, we have $G^{(0)} = G_s$ and $G^{(1)}$ is the following group, and so on.

### A.3. Dimension Reduction Algorithms

To finish the description of the two procedures, we fully describe the two dimension reduction algorithms both for the oblivious estimator $\hat{\pi}_{HT}$ and for the estimator $\hat{\pi}_{WM}$ with memory. These procedures were already introduced in Section 3. First, **DimensionReduction**($Y, \bar{P}, h, r$) considers the columns-wise mean of the restriction of $Y$ to the group $\bar{P}$ and detects high-variation regions of this vector.

**Algorithm 12 DimensionReduction**($Y, \bar{P}, h, r$)

**Require:** A sample $Y \in \mathbb{R}^{\bar{P} \times [d]}$, a set of experts $\bar{P}, h \in \mathcal{H}$ and $r$ in $\mathcal{R}$

**Ensure:** An encoded set including the high-variation regions $\tilde{Q}_{cp} := \tilde{Q}_{cp}(Y, \bar{P}, h, r) \subset Q_r$

1: $\bar{y}(\bar{P}) = \frac{1}{|\bar{P}|} \sum_{i \in \bar{P}} Y_{i, \cdot}$

2: $\tilde{r} = 8 \left[ \left( \frac{32 \zeta^2}{|\bar{P}| h^2} \log\left(\frac{2d}{\delta}\right) \right) \vee r \right]$

3: Initialize $\tilde{D}_{cp} = \emptyset$

4: **for** $k \in [d]$ **do**

5: Compute

$$ \tilde{C}_k(\bar{y}(\bar{P})) = \frac{1}{\tilde{r}} \left( \sum_{k'=k}^{k+\tilde{r}-1} \bar{y}_{k'}(\bar{P}) - \sum_{k'=k-\tilde{r}}^{k-1} \bar{y}_{k'}(\bar{P}) \right); \quad (40) $$

6: **end for**

7: $\tilde{D}_{cp} = \{k \in [d] : \tilde{C}_k(\bar{y}(\bar{P})) \geq h/4\}$

8: $\tilde{Q}_{cp} = \text{Encode} - \text{Set}(\tilde{D}_{cp}, r)$

9: **return** $\tilde{Q}_{cp}$

For the more involved dimension reduction procedure with memory **DimensionReduction** – **WM**, we compute the CUSUM statistic in larger groups $V \supset \bar{P}$ to reduce its variance and we also require that the estimated "width" of the group of experts is high enough. More precisely, given three sets of expert $V, V^+$ and $V^-$ and a sample $Y$, we consider the two following statistics, for any $k = 1, \dots, d$ and $r' \in \mathcal{R}$:

$$ \tilde{\Delta}_{k,r'}^{(\text{ext})}(V^+, V^-) = \frac{1}{2r'} \sum_{k'=k-r'}^{k+r'-1} \bar{y}_{k'}(V^+) - \bar{y}_{k'}(V^-); \quad \tilde{C}_{k,r'}^{(\text{ext})}(V) = \frac{1}{r'} \left( \sum_{k'=k}^{k+r'-1} \bar{y}_{k'}(V) - \sum_{k'=k-r'}^{k-1} \bar{y}_{k'}(V) \right). \quad (41) $$

Here, $\tilde{\Delta}_{k,r'}^{(\text{ext})}(V^+, V^-)$ computes the width - i.e. the difference - between the mean of experts in $V^+$ and the mean of experts in $V^-$. Since $V^+$ and $V^-$ are built in the algorithm below in such a way that experts in $\bar{P}$ are below those of $V^+$ and above those of $V^-$, $\tilde{\Delta}_{k,r'}^{(\text{ext})}(V^+, V^-)$ provides an upper bound of the width between the best expert in $\bar{P}$ and the worst expert in $\bar{P}$.
The algorithm **DimensionReduction-WM** described below builds a collection of sets $V^+$, $V^-$, and $V$ and detects questions such that both the CUSUM $C_{k,r'}^{(\text{ext})}(V)$ and the width $\hat{\Delta}_{k,r'}^{(\text{ext})}(V^+, V^-)$ are large enough. Further explanations are postponed to the analysis of the algorithm in Section F. Below, we write $[x]^{dya}$ for $2^{\lfloor \log_2(x) \rfloor}$.

**Algorithm 13 DimensionReduction-WM(Y, T, P̄, h, r)**

**Require:** A sample $Y \in \mathbb{R}^{n \times d}$, a tree $T$, a set $P̄$ included in a leaf $G$ of $T$ of type **0** or **1**, $h \in \mathcal{H}$ and $r \in \mathcal{R}$

**Ensure:** A set of blocks $\tilde{Q}_{WM} := \tilde{Q}_{WM}(Y, T, P̄, h, r) \subset Q_r$

1: $r_0 = 2^9 \log(4d|\mathcal{R}|/\delta) \frac{\zeta^2}{|P̄|h^2}$ and $\tilde{r} = 4([r_0]^{dya} \lor r)$

2: $(G^{(a)})_{a \in \mathbb{Z}} = \text{Order}(T, G)$

3: **for** $r_{cp} \in [4r, 2\tilde{r}] \cap \mathcal{R}$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;Set $a_{WM}^+ = \min\{a : |G^{(1)}| + \dots + |G^{(a)}| \ge 2^{11} \log(4d|\mathcal{R}|/\delta) \frac{\zeta^2}{r_{cp}h^2}\}$

5: &nbsp;&nbsp;&nbsp;&nbsp;Set $a_{WM}^- = \min\{a : |G^{(-1)}| + \dots + |G^{(-a)}| \ge 2^{11} \log(4d|\mathcal{R}|/\delta) \frac{\zeta^2}{r_{cp}h^2}\}$

6: Set

$$V_{r_{cp}}^{+} := V_{r_{cp}}^{+}(T, G, h) = \bigcup_{a=1}^{a_{WM}^{+}} G^{(a)} \quad \text{and} \quad V_{r_{cp}}^{-} := V_{r_{cp}}^{-}(T, G, h) = \bigcup_{a \in -a_{WM}^{-}}^{-1} G^{(a)} \quad (42)$$

7: **if** $r_{cp} > \tilde{r}$ **then**

8: &nbsp;&nbsp;&nbsp;&nbsp;Set $V_{r_{cp}} := V_{r_{cp}}(T, P̄, h) = P̄$

9: **else if** $r_{cp} \le \tilde{r}$ **then**

10: &nbsp;&nbsp;&nbsp;&nbsp;Set $V_{r_{cp}} := V_{r_{cp}}(T, P̄, h) = V_{r_{cp}}^{-} \cup V_{r_{cp}}^{+}$

11: **end if**

12: **end for**

13: $\tilde{Q}_{WM} = \emptyset$

14: **for** $r_{cp} \in [4r, \tilde{r}] \cap \mathcal{R}$ **do**

15: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{D}_{WM} = \emptyset$

16: **for** $k = 1, \dots, d$ **do**

17: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\hat{\Delta}_{k,r_{cp}}^{(\text{ext})} := \hat{\Delta}_{k,r_{cp}}^{(\text{ext})}(V_{r_{cp}}^{+}, V_{r_{cp}}^{-})$

18: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\tilde{C}_{k,2r_{cp}}^{(\text{ext})} := \tilde{C}_{k,2r_{cp}}^{(\text{ext})}(V_{2r_{cp}})$

19: **end for**

20: Update $\tilde{D}_{WM} = \{k \in [d] : \hat{\Delta}_{k,r_{cp}}^{(\text{ext})} \ge h/16 \text{ and } \tilde{C}_{k,2r_{cp}}^{(\text{ext})} \ge h/16\}$

21: Update $\tilde{Q}_{WM} = \tilde{Q}_{WM} \cup \text{Encode-Set}(\tilde{D}_{WM}, r)$

22: **end for**

23: **return** $\tilde{Q}_{WM}$

▷ The same set is defined in (26)

# Appendix B: Overview and organization of the proofs of Theorems 2.2 and 2.3

In this section, we divide the analysis of the procedures into several properties that will be proved to hold with high probability in the next sections.

## B.1. Definitions

Since we build our estimator using a hierarchical tree, we need to quantify the error that we suffer at each depth of the tree. For $i \in [n]$, we write $M_i = M_{i,}$ for the expert $i$. By definition of $\pi^*$, we recall that

$$M_{\pi^{*-1}(1)} \le M_{\pi^{*-1}(2)} \le \dots \le M_{\pi^{*-1}(n)}.$$

For a given group of experts $G$, we write $\pi_{\{G\}}^*$ for the oracle ordering in $[1, |G|]$ of the group $G$ according to $\pi^*$, that is for all $i, j \in G$, $\pi_{\{G\}}^*(i)$ and $\pi_{\{G\}}^*(j)$ belong to $[1, |G|]$ and

$$\pi_{\{G\}}^*(i) < \pi_{\{G\}}^*(j) \text{ iff } \pi^*(i) < \pi^*(j).$$
We say that a sequence of sets $\mathcal{G} = (G_1, \dots, G_\alpha)$ is an ordered partition of a set $S$ if $\{G_1, \dots, G_\alpha\}$ is a partition of $S$. For a given ordered partition $\{G_1, \dots, G_\alpha\}$ and $a \in [1, \alpha]$ and any $i \in G_a$ we write

$$ \pi_{\bar{\mathcal{G}}}(G_a) = \pi_{\bar{\mathcal{G}}}(i) := \sum_{a' < a} |G_{a'}| \quad \text{and} \quad \pi_{\mathcal{G}}^{+}(G_a) = \pi_{\mathcal{G}}^{+}(i) := \sum_{a' \le a} |G_{a'}|. \qquad (43) $$

If we are to build a permutation $\pi$ which is consistent with this ordered partition, then one easily checks that $\pi(i) \in [\pi_{\bar{\mathcal{G}}}(i) + 1, \pi_{\mathcal{G}}^{+}(i)]$. For simplicity, we write $G(i)$ for the group $G_a$ such that $i \in G_a$. For a given ordered partition $\mathcal{G} = (G_1, \dots, G_\alpha)$, we define the oracle permutation associated to $\mathcal{G}$ by

$$ \pi_{\bar{\mathcal{G}}}^{*}(i) = \pi_{\bar{\mathcal{G}}}(i) + \pi_{\{G(i)\}}^{*}(i). \qquad (44) $$

For example, $\pi_{\{[n]\}}^{*} = \pi^{*}$ is simply the true permutation. By definition, we have $\pi_{\bar{\mathcal{G}}}^{*}(G(i)) = [\pi_{\bar{\mathcal{G}}}(i) + 1, \pi_{\mathcal{G}}^{+}(i)]$. Given an ordered partition, $\pi_{\bar{\mathcal{G}}}^{*}$ is the best permutation we could hope for after any statistical treatment.

Given an ordered partition $\mathcal{G} = (G_1, \dots, G_\alpha)$, we define the random estimation of $\pi^{*}$ given $\mathcal{G}$ as $\hat{\pi}_{\mathcal{G}}(i)$ which is uniformly distributed in $[\pi_{\bar{\mathcal{G}}}(i) + 1, \pi_{\mathcal{G}}^{+}(i)]$:

$$ \hat{\pi}_{\mathcal{G}}(i) \in [\pi_{\bar{\mathcal{G}}}(i) + 1, \pi_{\mathcal{G}}^{+}(i)]. $$

Note that $\hat{\pi}_{\mathcal{G}}$ is not necessarily bijective.

## B.2. Deterministic Analysis

In this subsection, we analyze **TreeSort** (Algorithm 11) and we characterize the loss of the estimator $\hat{\pi}$ in terms of that of the trisections that are computed inside the subroutine **BlockSort**($\mathcal{Y}$, $\mathcal{T}$, $G$). This algorithm takes a subset $G$ of experts and computes two trisections of $G$. The first one

$$ (O, P, I) = (O_{\tau_{\infty}}, G \setminus (O_{\tau_{\infty}} \cup I_{\tau_{\infty}}), I_{\tau_{\infty}}) $$

is returned by the algorithm. The second one

$$ (\overline{O}, \overline{P}, \overline{I}) = (\overline{O}_{\tau_{\infty}}, G \setminus (\overline{O}_{\tau_{\infty}} \cup \overline{I}_{\tau_{\infty}}), \overline{I}_{\tau_{\infty}}) $$

is important for our analysis. From the definitions of the different procedures, one readily checks that $\overline{O} \subset O$ and $\overline{I} \subset I$. In fact, we shall prove later that, with high probability, the subsets $(O, P, I)$ and $(\overline{O}, \overline{P}, \overline{I})$ satisfy the following stronger property.

### Property 1.

1. $\{O, P, I\}$ and $\{\overline{O}, \overline{P}, \overline{I}\}$ are partitions of the leaf $G$ with $\overline{O} \subset O$, $\overline{I} \subset I$, and $P \subset \overline{P}$,

2. For $\omega = \pi_{\{O,P,I\}}^{*-1}\pi_{\{G\}}^{*}$, we have $\omega(i) = i$ for any $i \in \overline{O} \cup \overline{I}$.

3. For any $i \in O$ and $j \in I$, we have $\pi^{*}(i) < \pi^{*}(j)$.

4. We have $|O| \le |G|/2$ and $|I| \le |G|/2$.

The last claim states that all experts in $O$ are below all experts of $I$. The second claim can be understood as the fact that, if an expert $i$ belongs to $\overline{O}$, then all experts below $i$ belong to $O$.

Let $\mathcal{Y} = (Y^{(0)}, \dots, Y^{(6\tau_{\infty}-1)})$ be a sequence of $6\tau_{\infty}$ matrices in $\mathbb{R}^{n \times d}$. We say that **BlockSort** satisfies Property 1 on $(\mathcal{Y}, \mathcal{T}, G)$ if the two partitions $(O, P, I)$ and $(\overline{O}, \overline{P}, \overline{I})$ worked out in Algorithm 11 satisfy Property 1. We recall that by definition $\overline{O} \subset O$, $\overline{I} \subset I$, $P \subset \overline{P}$ so that $\overline{P}$ corresponds to the collection of experts that are either not sorted by **TreeSort** or are sorted with a small confidence.
For each $t = 0, \dots, t_\infty$, we write $\mathcal{T}_t$ for the hierarchical sorting tree at the beginning of step $t$ of **TreeSort**. Besides, we write $\mathcal{G}_t$ for the corresponding ordered partition obtained by taking the leaves of the tree $\mathcal{T}_t$ in increasing order in the ternary base $\{0, \mathbf{p}, \mathbf{1}\}$. We define the tree $\overline{\mathcal{T}}_{t_\infty}$ as the tree $\mathcal{T}_{t_\infty}$ where we replaced all the leaves $P$ of type $\mathbf{p}$ - at any depth - by $\overline{P}$, where $(\overline{O}, \overline{P}, \overline{I})$ has been worked out by **TreeSort** at the same time as $(O, P, I)$.

We also define

$$ 
\begin{aligned}
\mathcal{L}_t(\mathcal{T}_{t_\infty}) &= \{P \in \mathcal{G}_{t_\infty} : P \text{ is a nonempty leaf at depth } t \text{ of } \mathcal{T}_{t_\infty}\}; \\
\overline{\mathcal{L}}_t(\overline{\mathcal{T}}_{t_\infty}) &= \{\overline{P} \in \mathcal{G}_{t_\infty} : \overline{P} \text{ is a nonempty leaf at depth } t \text{ of } \overline{\mathcal{T}}_{t_\infty}\}.
\end{aligned} 
\quad (45) $$

For simplicity, we sometimes write $\mathcal{L}_t$ for $\mathcal{L}_t(\mathcal{T}_{t_\infty})$ and $\overline{\mathcal{L}}_t$ for $\overline{\mathcal{L}}_t(\overline{\mathcal{T}}_{t_\infty})$. $\mathcal{L}_t$ stands for the collection of experts that have not been sorted at the $t$-th iteration **TreeSort**. The sets in the collection $\overline{\mathcal{L}}_t$ are strictly larger and correspond to the collections of experts in $\overline{P}$ that are either not sorted by **TreeSort** or are sorted with less confidence. Let $M(\overline{P})$ be defined as the restriction of $M$ to the experts in $\overline{P}$, and $\overline{M}(\overline{P})$ the $|\overline{P}| \times d$ matrix with constant columns which correspond to the mean row of $M(\overline{P})$. The following proposition characterizes the loss of the final estimator estimator $\hat{\pi}_{\mathcal{G}_{t_\infty}}$ which is obtained from a hierarchical sorting tree in terms of the variance of the experts $M$ within the groups $\overline{P}$ in $\overline{\mathcal{L}}_t(\overline{\mathcal{T}}_{t_\infty})$.

**Proposition B.1** (Deterministic Analysis of **TreeSort**). Assume that at each step of **TreeSort**, the routine **BlockSort** applied to the data satisfies Property 1. Then, the error of $\hat{\pi} = \hat{\pi}_{\mathcal{G}_{t_\infty}}^{-1}$ is controlled as follows

$$ \| M_{\hat{\pi}^{-1}} - M_{\pi^{*-1}} \|_F^2 \le 10t_\infty \sum_{t=1}^{t_\infty} \sum_{\overline{P} \in \overline{\mathcal{L}}_t} \| M(\overline{P}) - \overline{M}(\overline{P}) \|_F^2. \quad (46) $$

Besides, the hierarchical tree $\mathcal{T}_{t_\infty}$ is valid (as in Definition 2) and all its non-empty leaves are of type $\mathbf{p}$.

Up to a normalization, $\|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2$ corresponds to the variance of $M$ within the group $\overline{P}$. The bound (46) expresses that the loss of a hierarchical sorting tree is controlled by the variance of the set $\overline{P}$ that are not sorted with confidence at each step of the algorithm. Also, we recall that $t_\infty = [\log(n)/\log(2)]$. This proposition only relies on Property 1 and on the construction of the tree. Hence, it applies both to the estimators $\hat{\pi}_{HT}$ and $\hat{\pi}_{WM}$.

The sets $(O, \overline{O}, I, \overline{I})$ built in **BlockSort** arise as unions of set $(L, U)$ and $(\overline{L}, \overline{U})$ that are computed by **DoubleTrisection** for a set $\overline{P}$ and a pivot $\gamma \in [1, |\overline{P}|]$. For this reason, we now state a desired property of the result of the algorithm that will enforce Property 1.

**Property 2** (Property on $(L, U)$ and $(\overline{L}, \overline{U})$). For $\overline{P}' = \overline{P} \setminus (\overline{L} \cup \overline{U})$ and $P' = \overline{P} \setminus (L \cup U)$, we have

1. $\overline{L} \subset L$, and $\overline{U} \subset U$,

2. if $\omega = \pi_{\{L, P', U\}}^{*-1} \pi_{\{\overline{P}\}}^*$ then for any $i \in \overline{L} \cup \overline{U}$ it holds that $\omega(i) = i$,

3. For any $i \in L$ and $j \in U$ we have $\pi_{\{\overline{P}\}}^*(i) < \gamma < \pi_{\{\overline{P}\}}^*(j)$.

We say that **DoubleTrisection** with $(\mathcal{Y}, \mathcal{T}, \overline{P}, \gamma)$ satisfies Property 2 if the corresponding subsets $(L, U)$ and $(\overline{L}, \overline{U})$ satisfy Property 2.

**Proposition B.2** (Deterministic Analysis of **BlockSort**). **BlockSort** satisfies Property 1 on $(\mathcal{Y}, \mathcal{T}, G)$ if, at each step of Algorithm 10, each call of **DoubleTrisection** satisfies Property 2.

In light of Propositions B.1 and B.2, it suffices to show that, with high probability, all applications of **DoubleTrisection** in the construction of the hierarchical sorting tree satisfy Property 2, and then to control the sum of within-group variances in (46).
### B.3. High probability Control of Property 2

We write in this part of the proof (this subsection), for simplicity, $\mathcal{Y} = (Y^{(1)}, \dots, Y^{(6)})$ for 6 independent matrices that are identically distributed as $Y = M + E$ in (1), where we recall that the entries of $E$ are centered, independent and $\zeta$-subgaussian.

Fix a hierarchical sorting tree $\mathcal{T}$ (recall Definition 1), a leaf $G$ of $\mathcal{T}$, a set $\overline{P} \subset G$, a pivot $\gamma \in \{1, \dots, |\overline{P}|\}$. Let $\mathcal{P}_2 := \mathcal{P}_2(\mathcal{T}, \overline{P}, \gamma, \beta_{\text{tris}}, \bar{\beta}_{\text{tris}})$ be the event holding true if **DoubleTrisection** satisfies Property 2 on $\mathcal{Y}$ for $(\mathcal{T}, \overline{P}, \gamma, \beta_{\text{tris}}, \bar{\beta}_{\text{tris}})$. The following proposition states that $\mathcal{P}_2$ holds with uniformly high probability.

**Proposition B.3.** For any $\mathcal{T}$, any leaf $G$, any $\overline{P} \subset G$, any pivot $\gamma \in [|\overline{P}|]$, we have $\mathbb{P}(\mathcal{P}_2) \ge 1 - 3|\mathcal{H}||\mathcal{R}|\delta$.

This result is valid for both versions of **DoubleTrisection** where we use the memory of the tree (estimator $\hat{\pi}_{WM}$) or not (estimator $\hat{\pi}_{HT}$). Recall that in **BlockSort** there are at most $\tau_{\infty}$ calls of **DoubleTrisection**. Since the construction of the hierarchical tree requires at most $2^{t_{\infty}+1}$ applications of **BlockSort**, we arrive at the following straightforward corollary of Propositions B.2, B.1 and B.3.

**Corollary B.4.** There exists an event $\xi$ of probability higher than $1 - 2^{t_{\infty}+1}3\tau_{\infty}|\mathcal{H}||\mathcal{R}|\delta$ such that all results of **BlockSort** within **TreeSort** satisfy Property 1. In particular, the tree $\mathcal{T}_{t_{\infty}}$ is a valid hierarchical sorting tree (as in Definition 2) whose non-empty leaves are all of type **p**. Besides, on this event we also have

$$ \|M_{\hat{\pi}_{G_{t_{\infty}}}^{-1}} - M_{\pi^{*-1}}\|_{F}^{2} \le 10t_{\infty} \sum_{t=1}^{t_{\infty}} \sum_{\overline{P} \in \mathcal{L}_{t}} \|M(\overline{P}) - \overline{M}(\overline{P})\|_{F}^{2}. \quad (47) $$

Again, this results applies to both variants of our procedure - with or without memory.

### B.4. Control of the Loss function

In contrast to the previous subsection, we now need to specify the dimension reduction scheme **DimensionReduction** (which corresponds to $\hat{\pi}_{HT}$) or **DimensionReduction-WM** (which corresponds to $\hat{\pi}_{WM}$) inside **DoubleTrisection** as the convergence rates depend on these quantities.

First we state the results for the method without memory: $\hat{\pi}_{HT}$.

**Proposition B.5.** Consider the oblivious hierarchical sorting tree estimator $\hat{\pi}_{HT}$. On the intersection of event $\xi$ (defined in Corollary B.4) and an event of probability higher than $1 - 5 \cdot 2^t \tau_{\infty} \delta$, it holds that

$$ \sum_{\overline{P} \in \mathcal{L}_t} \|M(\overline{P}) - \overline{M}(\overline{P})\|^2 \lesssim \zeta^2 \log^{8.5} \left( \frac{6nd}{\delta\zeta_-} \right) \left[ \frac{n^{2/3}d^{1/3}}{\zeta^{2/3}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right]. $$

Then we state the results for the method with memory: $\hat{\pi}_{WM}$.

**Proposition B.6.** Consider the hierarchical sorting tree estimator $\hat{\pi}_{WM}$. On the intersection of event $\xi$ (defined in Corollary B.4) and an event of probability higher than $1 - 5 \cdot 2^t \tau_{\infty} \delta$, it holds that

$$ \sum_{\overline{P} \in \mathcal{L}_t} \|M(\overline{P}) - \overline{M}(\overline{P})\|^2 \lesssim \zeta^2 \log^9 \left( \frac{6nd}{\delta\zeta_-} \right) \left[ \left( \frac{n^{3/4}d^{1/4}}{\zeta^{1/2}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} \wedge \frac{n^{2/3}\sqrt{d}}{\zeta^{1/3}} \right) + n \right]. $$

Now, we are in position to easily conclude the proof of Theorems 2.2 and 2.3.
**Proof of Theorem 2.2.** Let $\hat{\pi}_{HT} := \hat{\pi}_{G_{t_{\infty}}}$ denote the oblivious hierarchical sorting tree estimator. Combining Corollary B.4 with Proposition B.5 and a union bound over all $t = 0, \dots, t_{\infty} - 1$, it holds with probability higher than $1 - 8 \cdot 2^{t_{\infty}+1} \tau_{\infty} |\mathcal{H}||\mathcal{R}|\delta$ that

$$ \begin{aligned} \| M_{\hat{\pi}_{HT}^{-1}} - M_{\pi^{*-1}} \|_{F}^{2} &\lesssim t_{\infty}^{2} \zeta^{2} \log^{8.5} \left( \frac{2nd}{\delta\zeta_{-}} \right) \left[ \frac{n^{2/3}d^{1/3}}{\zeta^{2/3}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right] \\ &\lesssim \zeta^{2} \log^{10.5} \left( \frac{2nd}{\delta\zeta_{-}} \right) \left[ \frac{n^{2/3}d^{1/3}}{\zeta^{2/3}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right]. \end{aligned} $$

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

**Proof of Theorem 2.3.** Let $\hat{\pi}_{WM} := \hat{\pi}_{G_{t_{\infty}}}$ denote the hierarchical sorting tree where we use the memory to reduce the dimension (Algorithm **DimensionReduction-WM**). Combining Corollary B.4 with Proposition B.6 and a union bound on $t = 0, \dots, t_{\infty} - 1$, it holds with probability higher than $1 - 8 \cdot 2^{t_{\infty}+1} \tau_{\infty} |\mathcal{H}||\mathcal{R}|\delta$ that

$$ \| M_{\hat{\pi}_{WM}^{-1}} - M_{\pi^{*-1}} \|_{F}^{2} \lesssim \zeta^{2} \log^{11} \left( \frac{6nd}{\delta\zeta_{-}} \right) \left[ \left( \frac{n^{3/4}d^{1/4}}{\zeta^{1/2}} \wedge \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} \wedge \frac{n^{2/3}\sqrt{d}}{\zeta^{1/3}} \right) + n \right]. $$

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

In the next four sections, we prove the intermediary results. Propositions B.1-B.3 are relatively simple. The main difficulty and the key arguments lie in the proofs of Proposition B.5 and B.6 which are respectively in Sections D and F.

## Appendix C: Proofs of Propositions B.1, B.2, and B.3

**Proof of Proposition B.1.** First, we prove by induction that $\mathcal{T}_{t_{\infty}}$ is a valid hierarchical sorting tree. Besides, the last part of Property 1 enforces that the cardinality of any non-terminal node $G$ of $\mathcal{T}_{t_{\infty}}$ of depth $t$ is at most $n/2^t$. As a consequence, the cardinality of any non-terminal node at depth $t_{\infty}-1$ is at most 1 and its children $O$ and $I$ are therefore empty.

We control the error using a telescopic sum. Recall that, by convention, $\pi_{G_0}^* = \pi^*$. We start with the following inequality:

$$ \| M_{\hat{\pi}_{G_{t_{\infty}}^{-1}}} - M_{\pi^{*-1}} \|_{F}^{2} \le 2 \| M_{\hat{\pi}_{G_{t_{\infty}}}} - M_{\pi^{*-1}} \|_{F}^{2} + 2t_{\infty} \sum_{t=1}^{t_{\infty}} \| M_{\pi_{G_t}^{*-1}} - M_{\pi_{G_{t-1}}^{*-1}} \|_{F}^{2}. \quad (48) $$

Since, for any group $P$ in $G_{t_{\infty}}$, $\hat{\pi}_{G_{t_{\infty}}}$ sorts the elements of $P$ uniformly at random and $\pi_{G_{t_{\infty}}}^{*-1}$ acts as another permutation of $P$, we deduce from the triangular inequality that

$$ \begin{aligned} \| M_{\hat{\pi}_{G_{t_{\infty}}^{-1}}} - M_{\pi_{G_{t_{\infty}}^{*-1}}} \|_{F}^{2} &\le \sum_{P \in G_{t_{\infty}}} 2 \| M(P) - \bar{M}(P) \|_{F}^{2} = 2 \sum_{t=1}^{t_{\infty}} \sum_{P \in \mathcal{L}_{t}} \| M(P) - \bar{M}(P) \|_{F}^{2} \\ &\le 2 \sum_{t=1}^{t_{\infty}} \sum_{\bar{P} \in \mathcal{L}_{t}} \| M(\bar{P}) - \bar{M}(\bar{P}) \|_{F}^{2}, \end{aligned} \quad (49) $$

where we used in the last line that $P \subset \bar{P}$. For the second term in (48), remark that $\pi_{G_{t-1}}^*(P) = \pi_{G_t}^*(P)$ for any $P \in \mathcal{L}_{t-1}$ so that the error at step $t$ in the telescopic sum can be restricted to the groups $G$ that are trisected at step $t-1$:

$$ \sum_{t=1}^{t_{\infty}} \| M_{\pi_{G_t}^{*-1}} - M_{\pi_{G_{t-1}}^{*-1}} \|_{F}^{2} = \sum_{t=1}^{t_{\infty}} \sum_{G \in G_{t-1} \setminus (\cup_{s \ge 1} \mathcal{L}_{t-s})} \sum_{i \in G} \| M_{\pi_{G_t}^{*-1}(\pi_{G_{t-1}}^*(i))} - M_i \|_{2}^{2}. $$
Let $(O, P, I)$ be the trisection obtained at the $t$-th iteration when we apply **BlockSort** to a group $G \in \mathcal{G}_{t-1} \setminus (\cup_{s \ge 1} \mathcal{L}_{t-s})$. We also write $(\overline{O}, \overline{P}, \overline{I})$ for the more conservative trisection obtained at the end of **BlockSort**. For short, we write $\omega = \pi_{\mathcal{G}_t}^{*-1} \pi_{\mathcal{G}_{t-1}}^*$. We decompose the sum over $i \in G$:

$$ \sum_{i \in G} \|M_{\omega(i)} - M_i\|^2 = \sum_{i \in \overline{O}} \|M_{\omega(i)} - M_i\|^2 + \sum_{i \in \overline{I}} \|M_{\omega(i)} - M_i\|^2 + \sum_{i \in \overline{P}} \|M_{\omega(i)} - M_i\|^2, $$

By Property 1, all the experts in $\overline{O}$ and in $\overline{I}$ are perfectly sorted within $G$ by $\pi_{\mathcal{G}(t)}^{*-1}$. As a consequence, the two first sums in the right-hand side term of the above equality are equal to zero. To handle the last term, we introduce the row vector $m(\overline{P})$ as the mean of the experts of $M$ over $\overline{P}$:

$$ \begin{aligned} \sum_{i \in G} \|M_{\omega(i)} - M_i\|_2^2 &= \sum_{i \in \overline{P}} \|M_{\omega(i)} - M_i\|_2^2 \le 2 \sum_{i \in \overline{P}} (\|M_i - m(\overline{P})\|_2^2 + \|m(\overline{P}) - M_{\omega(i)}\|_2^2) \\ &= 4 \|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2, \end{aligned} $$

where we used in the last line that $\omega$ acts as a permutation of $\overline{P}$. Since $\overline{P} \in \mathcal{L}_t$, we obtain

$$ \sum_{t=1}^{t_{\infty}} \| M_{\pi_{\mathcal{G}_t}^{*-1}} - M_{\pi_{\mathcal{G}_{t-1}}^{*-1}} \|^2 \le 4 \sum_{t=1}^{t_{\infty}} \sum_{\overline{P} \in \mathcal{L}_t} \| M(\overline{P}) - \overline{M}(\overline{P}) \|_F^2. $$

Together with (48) and (49), this concludes the proof since $t_{\infty} \ge 1$. $\square$

**Proof of Proposition B.2.** Consider any data $\mathcal{Y}$, any tree $\mathcal{T}$ and any leaf $G$ of $\mathcal{T}$. Let $(O, P, I)$ and $(\overline{O}, \overline{P}, \overline{I})$ denote the trisections built in **BlockSort**($\mathcal{Y}$, $\mathcal{T}$, $G$). For any $\tau < \tau_{\infty}$, let $(L_{\tau}, U_{\tau})$, $(\overline{L}_{\tau}, \overline{U}_{\tau})$, $(O_{\tau}, I_{\tau})$ and $(\overline{O}_{\tau}, \overline{I}_{\tau})$ be defined as in Algorithm 10. We also write $P_{\tau} = G \setminus (O_{\tau} \cup I_{\tau})$ and $\overline{P}_{\tau} = G \setminus (\overline{O}_{\tau} \cup \overline{I}_{\tau})$. We only need to prove that, for all $\tau$, $(O_{\tau}, P_{\tau}, I_{\tau})$, and $(\overline{O}_{\tau}, \overline{P}_{\tau}, \overline{I}_{\tau})$ satisfy Property 1. Since

$$ O_{\tau} = \bigcup_{\tau' < \tau} L_{\tau'} \quad \text{and} \quad I_{\tau} = \bigcup_{\tau' < \tau} U_{\tau'} \quad \text{and} \quad \overline{O}_{\tau} = \bigcup_{\tau' < \tau} \overline{L}_{\tau'} \quad \text{and} \quad \overline{I}_{\tau} = \bigcup_{\tau' < \tau} \overline{U}_{\tau'}, $$

we easily deduce from Property 2 for $(L_{\tau}, U_{\tau})$ and $(\overline{L}_{\tau}, \overline{U}_{\tau})$ that the first part of Property 1 is satisfied for $(O_{\tau}, P_{\tau}, I_{\tau})$, and $(\overline{O}_{\tau}, \overline{P}_{\tau}, \overline{I}_{\tau})$.

Let us turn to the third and fourth parts of Property 1. Let us call $i_m$ the expert such that $\pi_{\{G\}}^*(i_m) = \lfloor|G|/2\rfloor$. In fact, we only need to prove that $\max_{i \in O_{\tau}} \pi_{\{G\}}^*(i) \le |G|/2$ and $\min_{i \in I_{\tau}} \pi_{\{G\}}^*(i) \ge |G|/2$. For this purpose, we prove by induction on $\tau$ that the pivot always satisfies $\pi_{\{\overline{P}\}}^{*-1}(\gamma) = i_m$ and that all the experts of $O_{\tau}$ (resp. $I_{\tau}$) are below (resp. above) $i_m$, where $\gamma$ depends on $\tau$ and is defined in Algorithm 10. Assume that this property holds at step $\tau$. Since $\overline{O}_{\tau}$ only contains experts that are below the median expert and since $\gamma = \lfloor|G|/2\rfloor - |\overline{O}_{\tau}|$, it follows that $\pi_{\{\overline{P}\}}^{*-1}(\gamma) = i_m$. Consider any $i \in O_{\tau+1}$. If $i \in O_{\tau}$, then $\pi_{\{G\}}^*(i) \le |G|/2$ by induction. If $i \in L_{\tau}$, then it follows from Property 2 that $i$ is below $i_m$, which in turn implies that $\pi_{\{G\}}^*(i) \le |G|/2$. By symmetry, the property also holds for $I_{\tau}$. We have proved the third and the fourth parts of Property 1.

Finally, we consider the second part of Property 1. Assume that the property holds at step $\tau$. This implies that, for any $i \in \overline{O}_{\tau}$, all experts below $i$ belong to $O_{\tau}$. Consider any expert $i \in \overline{O}_{\tau+1}$. If $i \in \overline{O}_{\tau}$, then, by induction, we have $\pi_{\{G\}}^*(i) = \pi_{\{O_{\tau}\}}^*(i) = \pi_{\{O_{\tau+1}\}}^*(i)$. Then, we turn to the case where $i$ belongs to $\overline{L}_{\tau} \subset \overline{P}_{\tau}$. Consider any $j \in G$ such that $\pi_{\{G\}}^*(j) \le \pi_{\{G\}}^*(i)$. If $j \in \overline{O}_{\tau}$, then we obviously have $j \in \overline{O}_{\tau+1}$. If $j \in \overline{P}_{\tau}$, then the second part of property 2 enforces
that $j \in L_\tau$ and therefore $j \in O_{\tau+1}$. Finally, it is not possible that $j \in \bar{I}_\tau$ since this enforces
that $\pi_{\{G\}}^*(j) > |G|/2 > \pi_{\{G\}}^*(i)$ and contradicts the hypothesis. We prove similarly that, for
any expert $i \in \bar{I}_{\tau+1}$, all experts $j$ above $i$ belong to $I_{\tau+1}$.

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

*Proof of Proposition B.3.* As **DoubleTrisection** is based on multiple applications of the **Pivot** algorithm, we start by considering the latter procedure.

Consider two sets $|\bar{P}| \subset [n]$ and $Q \subset [d]$ and a matrix $\Theta \in \mathbb{R}^{|\bar{P}| \times |Q|}$ which, up to the permutation $\pi_{\{\bar{P}\}}^*$, is bi-isotonic. Let $Z$ be a noisy observation of $\Theta$,

$$Z = \Theta + N, \qquad (50)$$

where the noise matrix $N$ is made of independent, centered, $\zeta$-subGaussian random variables. Let $w \in \mathbb{R}_+^Q$ be a non-zero vector with nonegative coordinates. we write $(L, U)$ and $(\bar{L}, \bar{U})$ for the result of **Pivot**($Z, w, \gamma$). We define the event $\mathcal{P}_3 := \mathcal{P}_3(\bar{P}, Q, w, \gamma)$ as the event on $Z$ such that $(L, U)$ and $(\bar{L}, \bar{U})$ satisfy Property 2.

We remind that $P' = \bar{P} \setminus (L \cup U)$, and $\bar{P}' = \bar{P} \setminus (\bar{L} \cup \bar{U})$.

**Lemma C.1.** For any non-zero vector $w \in \mathbb{R}_+^Q$, any pivot $\gamma \in \{1, \dots, |\bar{P}|\}$, we have $\mathbb{P}[\mathcal{P}_3] \ge 1 - \delta$. Besides, on the same event of probability at least $1 - \delta$, we have

$$\left| \langle \Theta_{i_{\gamma}}, -\Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \right| \le (2\zeta\sqrt{2} + \bar{\beta}_{\text{tris}}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} \quad \text{if } i \in P'. \qquad (51)$$

Before proving the lemma, let us explain why Proposition B.3 is easily deduced from it. The procedure **DoubleTrisection** calls at most $3|\mathcal{H}||\mathcal{R}|$ times **Pivot**. Note that, each time, we rely on an independent sample to choose the direction $w$ and to apply **Pivot**. Then, applying the Lemma, we derive that, with probability higher than $1 - 3|\mathcal{H}||\mathcal{R}|\delta$, each of these $3|\mathcal{H}||\mathcal{R}|$ sets $(L, U)$ and $(\bar{L}, \bar{U})$ satisfy Property 2. Hence, we only need to check that Property 2 is stable by union. If, both $(L^{(1)}, U^{(1)})$ and $(\bar{L}^{(1)}, \bar{U}^{(1)})$ and $(L^{(2)}, U^{(2)})$ and $(\bar{L}^{(2)}, \bar{U}^{(2)})$ satisfy Property 2, then one easily checks that the first and third part of Property 2 are also true for $(L, U) = (L^{(1)} \cup L^{(2)}, U^{(1)} \cup U^{(2)})$ and $(\bar{L}, \bar{U}) = (\bar{L}^{(1)} \cup \bar{L}^{(2)}, \bar{U}^{(1)} \cup \bar{U}^{(2)})$. Consider any expert $i$ in $\bar{L}$. Without loss of generality, we may assume that $i \in \bar{L}^{(1)}$ so that all experts below $i$ in $\bar{P}$ belong to $L^{(1)}$ by the second part of Property 2. As a consequence, all these experts below $i$ belong to $L$ and we deduce that the second part of Property 2 holds. Similarly, we deal with experts $i \in \bar{U}$. This concludes the proof of Proposition B.3.

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

*Proof of Lemma C.1.* Since the noise matrix in (50) is made of independent $\zeta$-subGaussian random variables, it follows from a union bound, that with probability higher than $1 - \delta$, we have

$$\left| \langle Z_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle - \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \right| \le \zeta \sqrt{2 \log \left( \frac{2|\bar{P}|}{\delta} \right)}.$$

simultaneously for all $i$ in $\bar{P}$. Since the entries of $w$ are non-negative, the quantities $\langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle$ are ordered according the permutation $\pi_{\{\bar{P}\}}^*$. Denote $i_{\gamma} = \pi_{\{\bar{P}\}}^{*-1}(\gamma)$ and $\hat{i}_{\gamma}$ the index of $\gamma$-th value of $\langle Z_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle$ for $i \in \bar{P}$. Since at least $\gamma$ experts satisfy $\langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \le \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle$, we deduce from the above uniform deviation inequality that

$$\langle Z_{\hat{i}_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \le \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle + \zeta \sqrt{2 \log \left( \frac{2|\bar{P}|}{\delta} \right)}.$$
By symmetry, we deduce that

$$ \left| \langle Z_{\hat{i}_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle - \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \right| \le \zeta \sqrt{2 \log \left( \frac{2|\bar{P}|}{\delta} \right)}. $$

As a consequence, we have

$$ \langle \Theta_{i, \cdot}, \frac{w}{\|w\|_2} \rangle < \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle - (\beta_{\text{tris}} - 2\zeta\sqrt{2}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} \quad \text{if } i \in L; $$

$$ \langle \Theta_{i, \cdot}, \frac{w}{\|w\|_2} \rangle > \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle + (\beta_{\text{tris}} - 2\zeta\sqrt{2}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} \quad \text{if } i \in U; $$

$$ \left| \langle \Theta_{i, \cdot} - \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle \right| \le (2\zeta\sqrt{2} + \beta_{\text{tris}}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} \quad \text{if } i \in P'. \quad (52) $$

The same inequalities hold for $\bar{L}$, $\bar{U}$, and $\bar{P}'$ provided that we replace $\beta_{\text{tris}}$ by $\bar{\beta}_{\text{tris}}$. It remains to show that $(L, U)$ and $(\bar{L}, \bar{U})$ satisfy Property 2. The first part of the property is obvious. Since $\beta_{\text{tris}} \ge 2\sqrt{2\zeta}$, one observes that $\pi_{\{\bar{P}\}}^*(i) < \pi_{\{\bar{P}\}}^*(i_{\gamma}) = \gamma$ if $i \in L$. Similarly, $\pi_{\{\bar{P}\}}^*(i) > \gamma$ if $i \in U$ and the third part of Property 2 follows. Turning to the second part of the property, we consider without loss of generality some $i \in \bar{L}$ and we need to show that all $j$ satisfying $\pi_{\{\bar{P}\}}^*(j) \le \pi_{\{\bar{P}\}}^*(i)$ belong to $L$. First, such a $j$ does not belong to $U$ since $\pi_{\{\bar{P}\}}^*(i) < \gamma$. Since $i \in \bar{L}$, we deduce that

$$ \langle \Theta_{i, \cdot}, \frac{w}{\|w\|_2} \rangle < \langle \Theta_{i_{\gamma}, \cdot}, \frac{w}{\|w\|_2} \rangle - (\bar{\beta}_{\text{tris}} - 2\sqrt{2\zeta}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)}, $$

which implies that

$$ \left| \langle \Theta_{i_{\gamma}, \cdot}, -\Theta_{j, \cdot}, \frac{w}{\|w\|_2} \rangle \right| > (\bar{\beta}_{\text{tris}} - 2\sqrt{2\zeta}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} \ge (2\zeta\sqrt{2} + \beta_{\text{tris}}) \sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)}. $$

which in light of (52) implies that $j \notin P'$. We have proved that $j$ belongs to $L$. Hence, Property 2 holds, which concludes the proof. $\square$

## Appendix D: Proof of Proposition B.5

In this section, we prove Proposition B.5 and thereby control the loss of the estimator $\hat{\pi}$ with simple dimension reduction. For this purpose, we analyze each step of the algorithm. In Appendix D.2, we first prove that, by detecting the high-variation regions of $M$, we are able to aggregate $M$ at some scale $r$ without decreasing much the variation of $M$. This allows us to drastically reduce the dimension of the problem. Then, in Sections D.3 and D.3, we show that, unless this aggregated matrix $\Theta$ has small variations, **DoubleTrisection** – PCA and **Pivot** will remove some experts so that the corresponding new aggregated matrix $\Theta'$ exhibit significantly smaller variations. As a consequence, after a polylogarithmic number of iterations of the procedure, the variations of the matrix $M$ restricted to the remaining experts of $\bar{P}$ is small enough.
## D.1. Notation

As the arguments rely on considering aggregation of the matrix at different scales, we recall some notation. Let $Y = M + E$ denote a sample of the original matrix. For a set $P \subset [n]$ of experts and a set $Q$ of blocks of questions and a scale $r \in R$, we respectively denote

$$ 
\begin{aligned}
Z(P, Q, r) &= \textbf{Encode} - \textbf{Matrix}(Y, P, Q, r) \in \mathbb{R}^{P \times Q} \\
\Theta(P, Q, r) &= \textbf{Encode} - \textbf{Matrix}(M, P, Q, r) \\
N(P, Q, r) &= \textbf{Encode} - \textbf{Matrix}(E, P, Q, r) ,
\end{aligned}
$$

the aggregations of $Y$, $M$, and $E$ at scale $r$ so that

$$ Z(P, Q, r) = \Theta(P, Q, r) + N(P, Q, r) . $$

By definition of **Encode** – **Matrix**, all the entries of $N$ are independent and $\zeta$-subGaussian. For any $p \times q$ matrix $A$, we define $\bar{a}$ as the row vector corresponding to the column-wise mean of $A$, that is $\bar{a}_j = \frac{1}{q} \sum_{i=1}^q A_{i,j}$. Besides, we write $\bar{A}$ for $p \times q$ matrix whose experts are all equal to $\bar{a}$.

## D.2. Analysis of DimensionReduction

In this subsection, we mainly state, that for any $h \in \mathcal{H}$, $r \in \mathcal{R}$, the set $\tilde{Q}_{cp}$ (which depends on $h$, $r$) detects the high-variation regions of $M$ with high probability. Then, we show in Lemma D.2, that for, for some $(h, r) \in \mathcal{H} \times \mathcal{R}$, the aggregation of $M$ at scale $r$ and at these high-variation regions contains most of the variance of $M$. This motivates us to work with this aggregated matrix henceforth. Consider any set $\bar{P}$ of experts and a sample

$$ Y(\bar{P}) = M(\bar{P}) + E(\bar{P}) . $$

Fix any scale $r \in \mathcal{R}$ and any height $h \in \mathcal{H}$. Recall the two quantities $r_0$ and $\tilde{r}$ defined in **DimensionReduction** by

$$ r_0 = 32\zeta^2 \log\left(\frac{2d}{\delta}\right) \frac{1}{|\bar{P}|h^2} \quad \text{and} \quad \tilde{r} = 8(|r_0| \vee r) . \qquad (53) $$

In a nustshell, $r_0$ stands for the minimal scale at which a variation of order $h$ in the mean $\bar{m}(\bar{P}) = \mathbb{E}[\bar{y}(\bar{P})]$ can be statistically detected. This is why we consider empirical variations of $\bar{y}(\bar{P})$ at scale $\tilde{r} \ge (r_0 \vee r)$ in Algorithm **DimensionReduction** to possibly detect variations at scale $r$.

The purpose of this subsection is to prove, that with high probability, the collections $\tilde{Q}_{cp}(h, r) = \textbf{DimensionReduction}(Y, \bar{P}, h, r)$ of selected blocks of length $r$ is not too large and that there exists at least one $(h, r) \in \mathcal{H} \times \mathcal{R}$ such that the aggregation of $M(\bar{P})$ at scale $r$ restricted to the blocks $\tilde{Q}_{cp}(h, r)$ captures most of the variance of $M(\bar{P})$.

For this purpose, we recall the CUSUM statistics introduced in **DimensionReduction** and we introduce its population counterpart. Given positive integers $k \in [d]$ and $r > 0$, consider

$$ \tilde{\mathbf{C}}_{k,r} = \frac{1}{r} \left( \sum_{k'=k}^{k+r-1} \bar{y}_{k'}(\bar{P}) - \sum_{k'=k-r}^{k-1} \bar{y}_{k'}(\bar{P}) \right) \quad \text{and} \quad \mathbf{C}_{k,r}^* = \frac{1}{r} \left( \sum_{k'=k}^{k+r-1} \bar{m}_{k'}(\bar{P}) - \sum_{k'=k-r}^{k-1} \bar{m}_{k'}(\bar{P}) \right) . $$

Equipped with this notation, we define $\tilde{D}_{cp}$ as in the algorithm as the set of positions $d$ such that the association CUSUM statistic is above the threshold, and $D_{cp}^*$ and $\tilde{D}_{cp}^*$ as some
population versions of $\hat{D}_{cp}$, but with different tuning parameters:

$$ \hat{D}_{cp}(h, r) = \left\{ k \in [d] : \hat{C}_{k, \tilde{r}} \geq \frac{1}{4}h \right\}, \qquad (54) $$

$$ D_{cp}^{*}(h, r) = \left\{ k \in [d] : \mathbf{C}_{k, 8r}^{*} \geq \frac{1}{2}h \right\}; \quad \overline{D}_{cp}^{*}(h, r) = \left\{ k \in [d] : \mathbf{C}_{k, \tilde{r}}^{*} \geq \frac{1}{8}h \right\}. \qquad (55) $$

Then, we consider the collection of blocks $\hat{Q}_{cp}(h, r)$, $Q_{cp}^{*}(h, r)$, and $\overline{Q}_{cp}^{*}(h, r)$ of size $r$ that are associated with these positions. In terms of our algorithms, this means that $Q_{cp}^{*}(h, r) = \text{Encode} - \text{Set}(D_{cp}^{*}, r)$, $\overline{Q}_{cp}^{*}(h, r) = \text{Encode} - \text{Set}(\overline{D}_{cp}^{*}, r)$, and $\hat{Q}_{cp}(h, r) = \text{Encode} - \text{Set}(\hat{D}_{cp}, r)$. The first proposition states that, with high probability, $\hat{Q}_{cp}(h, r)$ is sandwiched between $Q_{cp}^{*}(h, r)$ and $\overline{Q}_{cp}^{*}(h, r)$, so that, on the corresponding event, it is sufficient to study these two quantities.

**Lemma D.1.** For all $h, r$, the event $\xi_{cp} := \xi_{cp}(\overline{P}, h, r)$ defined by

$$ Q_{cp}^{*} \subset \hat{Q}_{cp} \subset \overline{Q}_{cp}^{*}, \qquad (56) $$

holds true with probability at least $1 - \delta$.

Then, we show that there are not too many significant blocks in $\overline{Q}_{cp}^{*}$. The proof is based on the fact that the row vector $\overline{m}(\overline{P})$ is isotonic and lies in $[0, 1]$. As a consequence, there cannot exist two many regions where the variations of $\overline{m}(\overline{P})$ is large.

**Lemma D.2.** For all $h \in \mathcal{H}$ and all $r \in \mathcal{R}$, we have

$$ |\overline{Q}_{cp}^{*}| \leq \frac{64\tilde{r}}{rh}. \qquad (57) $$

The next lemma states that, at least for a height $h \in \mathcal{H}$ and a scale $r \in \mathcal{R}$, the aggregation of $M$ at scale $r$ and restricted to the regions $Q_{cp}^{*}(h, r)$ of significant variations contains almost all the variance of the signal.

For any number $\theta$ and any $\eta > 0$, we define $[\theta]_{\eta} = (-1)^{\text{sgn}(\theta)}\eta\mathbf{1}_{|\theta|\geq\eta}$. For any matrix $\Theta$, we write $[\Theta]_{\eta}$ for the thresholded matrix with coefficients $[\Theta_{i,j}]_{\eta}$.

**Lemma D.3.** For any set $\overline{P} \subset [n]$ and any bi-isotonic matrix $M \in [0, 1]^{n \times d}$, there exist $r \in \mathcal{R}$ and $h \in \mathcal{H}$ such that

$$ \|M(\overline{P}) - \overline{M}(\overline{P})\|_{F}^{2} \leq 16\zeta^{2} + 96|\mathcal{R}\|\mathcal{H}| \left\| [\Theta(\overline{P}, Q_{cp}^{*}) - \overline{\Theta}(\overline{P}, Q_{cp}^{*})]_{\sqrt{rh}} \right\|_{F}^{2}, \qquad (58) $$

The proof of the above lemmas is postponed to Appendix E.1.

### D.3. Analysis of Pivot based on the row sums

We consider a specific subset $\overline{P}$ of experts, a subset $Q$ of blocks of questions, the corresponding aggregated model

$$ Z(\overline{P}, Q) = \Theta(\overline{P}, Q) + N(\overline{P}, Q) \in \mathbb{R}^{\overline{P} \times Q}, \qquad (59) $$

and a pivot $\gamma \in [1, |\overline{P}|]$. Let $(\overline{L}, \overline{U}) = \text{Pivot}(Z, \mathbf{1}_{Q}, \gamma)$ be the conservative result of **Pivot** based on the row sums of $Z(\overline{P}, Q)$ and let $\overline{P}' = \overline{P} \setminus (\overline{L} \cup \overline{U})$ be the subgroup of experts which have not been classified by **Pivot**. The following proposition states that, provided that for some $\eta$ the norm $\| [\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)]_{\eta} \|_{F}^{2}$ is large enough compared to $|\overline{P}|\sqrt{|\overline{Q}|}$, the resulting matrix $\Theta(\overline{P}', Q) - \overline{\Theta}(\overline{P}', Q)$ after **Pivot** has a significantly smaller norm. We shall often use the following quantity.

$$ \phi_{l1} = 2(2\zeta\sqrt{2} + \beta_{\text{tris}}) \leq 29\zeta. \qquad (60) $$
**Proposition D.4.** Consider any $\overline{P} \subset [n]$, any $r \in \mathbb{R}$, and any subset $Q \subset Q_r$. Also, fix any $\eta > 0$ and any $\phi > 0$. If

$$ \| [\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)]_{\eta} \|_{F}^{2} \geq \frac{1}{\phi} \| \Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q) \|_{F}^{2} \geq 8 \phi_{l_1} \eta \sqrt{\log(\frac{2|\overline{P}|}{\delta})} |\overline{P}| \sqrt{|Q|} , $$

then, with probability higher than $1 - \delta$, we have

$$ \|\Theta(\overline{P}', Q) - \overline{\Theta}(\overline{P}', Q)\|_{F}^{2} \leq \left(1 - \frac{1}{16\phi}\right) \|\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)\|_{F}^{2}. $$

### D.4. Analysis of DoubleTrisection - PCA

In this subsection, we state the main result regarding the trisection of a set $\overline{P}$ based on the first singular vector of a suitable matrix. We start from a subset $\overline{P}$ of experts. In **DoubleTrisection - Local**, we start applying **Pivot** and define $\tilde{P} = \overline{P} \setminus (\overline{L}_{cp} \cup \overline{U}_{cp})$ as the set of experts that have not been classified by **Pivot**. We are given four independent samples $\mathcal{Z} = (Z^{(1)}, Z^{(2)}, Z^{(3)}, Z^{(4)})$ according to the aggregated model (59). The first three samples are restricted to $\tilde{P}$, whereas the last one concerns $\overline{P}$. Fix $\gamma \in [1, |\overline{P}|]$. We consider $(\overline{L}, \overline{U}) = \text{DoubleTrisection - PCA}(\mathcal{Z}, \gamma)$ and $\overline{P}' = \tilde{P} \setminus (\overline{L}_{pca} \cup \overline{U}_{pca})$ the set of experts that have not been classified by **DoubleTrisection - PCA**.

Recall the definition (60) of $\phi_{l_1}$. Henceforth, the matrix $\Theta(\tilde{P}, Q)$ is said to be undistinguishable in $l_1$-norm if it satisfies

$$ \max_{i,j \in \overline{P}} \| \Theta_{i, \cdot}(\tilde{P}, Q) - \Theta_{j, \cdot}(\tilde{P}, Q) \|_{1} \leq \phi_{l_1} \sqrt{|Q| \log\left(\frac{2|\overline{P}|}{\delta}\right)}. \quad (61) $$

Since, up to permutation of its experts, the matrix $\Theta(\tilde{P}, Q)$ is bi-isotonic, the $l_1$ norm $\|\Theta_{i, \cdot}(\tilde{P}, Q) - \Theta_{j, \cdot}(\tilde{P}, Q)\|_{1}$ is simply the difference of the row sums of $\Theta(\tilde{P}, Q)$. Since $\tilde{P}$ has been deduced from $\overline{P}$ by applying **Pivot**($\mathcal{Z}, \gamma$), we can safely assume that $\Theta(\tilde{P}, Q)$ is undistinguishable in $l_1$-norm with high probability – see the next subsection for a proper justification.

The next result states that, if $\Theta(\tilde{P}, Q)$ is undistinguishable in $l_1$-norm and if the Frobenius norm of $\Theta(\tilde{P}, Q) - \overline{\Theta}(\tilde{P}, Q)$ is large enough, then the corresponding matrix $\Theta(\overline{P}', Q)$ obtained after trisection has a significantly smaller Frobenius norm.

**Proposition D.5.** Let $\overline{P} \subset [n]$ and $Q \subset [d]$. If $\Theta(\tilde{P}, Q)$ is undistinguishable in $l_1$-norm and if

$$ \|\Theta(\tilde{P}, Q) - \overline{\Theta}(\tilde{P}, Q)\|_{F}^{2} \geq 2 \cdot 10^{5} \zeta^{2} \log^{3}\left(\frac{6nd}{\delta\zeta_{-}}\right) \left(\sqrt{|\tilde{P}||Q|} + |\tilde{P}|\right), \quad (62) $$

then, with probability higher than $1 - 3\delta$, we have

$$ \|\Theta(\overline{P}', Q) - \overline{\Theta}(\overline{P}', Q)\|_{F}^{2} \leq \left(1 - \frac{1}{200 \log^{2}(nd/\zeta_{-})}\right) \|\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)\|_{F}^{2}. $$

Then, we gather the two previous results to analyze the routine **DoubleTrisection - Local**. Fix any $\overline{P} \subset [n]$ and $Q \subset [p]$. Let $\mathcal{Z} = (Z^{(1)}(\overline{P}, Q), Z^{(2)}(\overline{P}, Q), Z^{(3)}(\overline{P}, Q), Z^{(4)}(\overline{P}, Q), Z^{(5)}(\overline{P}, Q))$ be five independent samples of the model (50). Fix any $\gamma \in [1, |\overline{P}|]$. Let $(\overline{L}, \overline{U})$ be the conservative result of **DoubleTrisection - Local**($\mathcal{Z}, \gamma$) and $\overline{P}' = \tilde{P} \setminus (\overline{L} \cup \overline{U})$. In the following, we write $\Theta(\overline{P}, Q_r)$ for the aggregation of $M(\overline{P})$ at all blocks of size $r$.
**Corollary D.6.** Fix any $r \in \mathbb{R}$. If, for some $\overline{P} \subset [n]$, $Q \subset Q_r$, and $\eta > 0$, $\Theta(\overline{P}, Q)$ satisfies

$$ \left\{ \begin{array}{l} \| [\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)]_{\eta} \|_{F}^{2} \ge \frac{1}{120 \log^{2}\left(\frac{nd}{\delta\zeta_{-}}\right)} \|\Theta(\overline{P}, Q_{r}) - \overline{\Theta}(\overline{P}, Q_{r})\|_{F}^{2} \\ \|\Theta(\overline{P}, Q) - \overline{\Theta}(\overline{P}, Q)\|_{F}^{2} \ge 4 \cdot 10^{5} \zeta^{2} \log^{3}\left(\frac{6nd}{\delta\zeta_{-}}\right) \left[ \frac{\eta}{\zeta} |\overline{P}| \sqrt{|Q|} \wedge \left(\sqrt{|\overline{P}|\|Q|} + |\overline{P}|\right) \right], \end{array} \right. \quad (63) $$

then, with probability higher than $1 - 4\delta$, we have

$$ \|\Theta(\overline{P}', Q_r) - \overline{\Theta}(\overline{P}', Q_r)\|_{F}^{2} \le \left( 1 - \frac{1}{3 \cdot 10^{5} \log^{4}\left(\frac{nd}{\delta\zeta_{-}}\right)} \right) \|\Theta(\overline{P}, Q_r) - \overline{\Theta}(\overline{P}, Q_r)\|_{F}^{2}. $$

## D.5. Analysis of BlockSort

Next, we combine the results of the previous sections to control the error of **BlockSort**. We are given a collection $\mathcal{Y}$ of $6\tau_\infty$ samples of the model (1) and a valid hierarchical sorting tree $\mathcal{T}$ of depth $t$ that we consider as fixed. Then, we take a leaf $G$ of $\mathcal{T}$ with maximal depth and we consider $(O, P, I) = \textbf{BlockSort}(\mathcal{Y}, \mathcal{T}, G)$ the trisection of $G$, as well as $\overline{P} \supset P$ the more conservative intermediary set. In this section, we provide a high-probability control of $\overline{P}$.

For any height $h \in \mathcal{H}$ and scale $r \in \mathbb{R}$, recall that $Q_{cp}^*$ is the subset(defined in Appendix D.2) of block of questions at scale $r$ such that the mean $\overline{m}(\overline{P})$ increases by at least $h/2$. Also recall the superset $\overline{Q}_{cp}^* \supset Q_{cp}^*$.

At a high level, the next proposition states that, after $\tau_\infty$ iterations of the **DoubleTrisection** – **Local** routines at all scales $r \in \mathbb{R}$ and all heights $h \in \mathcal{H}$, the size $\|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2$ is quite small. This is mainly due to the fact that, by Lemma D.3, at each step $\tau$, there exists some $(r, h) \in \mathbb{R} \times \mathcal{H}$ such that the norm of the thresholded aggregated matrix $[\Theta(\overline{P}_\tau, Q_{cp}^*) - \overline{\Theta}(\overline{P}_\tau, Q_{cp}^*)]_{\sqrt{rh}}$ is of the same order as $\|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2$. By Lemma D.1, the estimated blocks $\overline{Q}_{cp}^*$ contain $Q_{cp}^*$ with high probability. Hence, unless the norm of the thresholded aggregated matrix is small, we derive from corollary D.6 that the norm of this aggregated matrix has contracted at step $\tau + 1$. Hence, after $\tau_\infty$ steps, one could expect that the norm of $\|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2$ is small. In fact, both the statement and the proof of this proposition are slightly more involved because we need to keep track of the scales and heights of interest. Define the function $\Psi(p, r, h, q)$ by

$$ \Psi(p, r, h, q) = \frac{hp\sqrt{rq}}{\zeta} \wedge (\sqrt{pq}) + p. \quad (64) $$

**Proposition D.7.** With probability higher than $1 - 5\tau_\infty\delta$, there exists a subset $\overline{P}^\dagger$ such that $\overline{P} \subseteq \overline{P}^\dagger \subseteq G$ and the following property holds. For some $r^\dagger \in \mathbb{R}$ and some $h^\dagger \in \mathcal{H}$, upon writing $Q_{cp}^\dagger = Q_{cp}^*(\overline{P}^\dagger, h^\dagger, r^\dagger)$ and $\overline{Q}_{cp}^\dagger = \overline{Q}_{cp}^*(\overline{P}^\dagger, h^\dagger, r^\dagger)$, we have simultaneously

$$ \| [\Theta(\overline{P}^\dagger, Q_{cp}^\dagger) - \overline{\Theta}(\overline{P}^\dagger, Q_{cp}^\dagger)]_{\sqrt{r^\dagger h^\dagger}} \|_F^2 \le 4 \cdot 10^5 \zeta^2 \log^3 \left( \frac{6nd}{\delta\zeta_{-}} \right) \Psi(|\overline{P}^\dagger|, r^\dagger, h^\dagger, |\overline{Q}_{cp}^\dagger|) ; \quad (65) $$

$$ \|M(\overline{P}^\dagger) - \overline{M}(\overline{P}^\dagger)\|_F^2 \le 16\zeta^2 + 96|\mathcal{R}|\|\mathcal{H}\|\|[\Theta(\overline{P}^\dagger, Q_{cp}^\dagger) - \overline{\Theta}(\overline{P}^\dagger, Q_{cp}^\dagger)]_{\sqrt{r^\dagger h^\dagger}}\|_F^2. \quad (66) $$

In other words, there exists a superset $\overline{P}^\dagger$ of $\overline{P}$ such that, for a suitable height and scale, at the high-variation regions, both the original matrix $M(\overline{P}^\dagger)$ and the thresholded aggregated matrix are controlled at the level $\Psi(|\overline{P}^\dagger|, r^\dagger, h^\dagger, |\overline{Q}_{cp}^\dagger|)$. The virtue of the above result is that it easily adapts to the block sorting variant with memory. Unfortunately, the rate $\Psi(|\overline{P}^\dagger|, r^\dagger, h^\dagger, |\overline{Q}_{cp}^\dagger|)$ is a bit difficult to handle. In the next corollary, we replace it by a simpler but cruder bound that only depends on $|G|$, $h^\dagger$ and $d$.
**Corollary D.8.** Under the same event of probability higher than $1 - 5\tau_\infty\delta$ as in the previous proposition, the set $\overline{P}^\dagger$, the scale $r^\dagger$, and the height $h^\dagger$ also satisfy

$$ \| [\Theta(\overline{P}^\dagger, Q_{cp}^\dagger) - \overline{\Theta}(\overline{P}^\dagger, Q_{cp}^\dagger)]_{\sqrt{r^\dagger}h^\dagger} \|_F^2 \lesssim \zeta^2 \log^{3.5} \left( \frac{6nd}{\delta\zeta_-} \right) \left[ \frac{h^\dagger|G|\sqrt{d}}{\zeta} \wedge \sqrt{|G|d} \wedge \sqrt{\frac{|G|}{h^\dagger}} + |G| \right], \quad (67) $$

where we recall that $G$ is the initial group.

**D.6. Analysis of the complete procedure TreeSort**

We are now equipped to prove Proposition B.5.

*Proof of Proposition B.5.* Let us fix an integer $t \le t_\infty$ and let us consider the collection $\overline{\mathcal{L}}_t$ of the $2^t$ groups $\overline{P}$ that are not sorted with confidence. Let us apply Proposition D.7 to each of these sets $\overline{P}$. In view of this proposition, we define $(\overline{P}^\dagger, h^\dagger, r^\dagger)$ as well as $Q_{cp}^\dagger$. We also define

$$ s^\dagger = \left| \left\{ l : \exists i, j \in \overline{P}^\dagger \text{ s.t. } [\Theta_{i,l}(\overline{P}^\dagger, Q_{cp}^\dagger) - \Theta_{j,l}(\overline{P}^\dagger, Q_{cp}^\dagger)]_{\sqrt{r^\dagger}h^\dagger} \neq 0 \right\} \right|. $$

In a nustshell, $s^\dagger$ is the number of columns of the thresholded aggregated matrix which are not equal to zero.

**Definition 3.** Define the dyadic collection $\mathcal{S} = \{1, 2, 4, \dots, 2^{\lfloor \log_2(d) \rfloor}\}$. For any $s \in \mathcal{S}$, $r \in \mathcal{R}$, and $h \in \mathcal{H}$, we consider the collection $\mathcal{P}^*(h, r, s) \subset \overline{\mathcal{L}}_t$ satisfying $r^\dagger = r, h^\dagger = h$, and $s^\dagger \in [s, 2s)$.

The following lemma controls the cardinality of $\mathcal{P}^*(h, r, s)$. This bound mainly relies on the facts that the matrix $M$ is, up to a row permutation, bi-isotonic and that its entries lie in $[0, 1]$.

**Lemma D.9.** Assume that there exists an ordering $\sigma$ of $\overline{\mathcal{L}}_t$ that orders all groups $\overline{P}$'s. In other words, for any $r \le s$, any expert $i \in \overline{P}_{\sigma(r)}$ is below any expert $j \in \overline{P}_{\sigma(s)}$. Then, upon this assumption,

$$ |\mathcal{P}^*(h, r, s)| \le \frac{2d}{hrs} \wedge 2^t \le \sqrt{\frac{d2^{t+1}}{hrs}}, $$

for any $h \in \mathcal{H}$, $r \in \mathcal{R}$, and $s \in \mathcal{S}$.

In fact, all the collections $\overline{\mathcal{L}}_t$ with $t = 0, \dots, t_\infty$ satisfy the assumption in the above under the event $\xi$ defined in Corollary B.4 –see the proof of Proposition B.3.

Putting everything together and summing over the groups $\mathcal{P}^*(h, r, s)$, we derive from Proposition D.7 and Corollary D.8 that, with probability higher than $1 - 5 \cdot 2^t \tau_\infty \delta$, we have

$$ \begin{align*} \sum_{\overline{P} \in \overline{\mathcal{L}}_t} \|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2 &\le 16\zeta^2 |\overline{\mathcal{L}}_t| + 96|\mathcal{R}|\|\mathcal{H}\| \sum_{h,r,s} \sum_{\overline{P} \in \mathcal{P}^*(h,r,s)} \| [\Theta(\overline{P}, Q_{cp}^\dagger) - \overline{\Theta}(\overline{P}, Q_{cp}^\dagger)]_{\sqrt{rh}} \|_F^2 \\ &\lesssim \zeta^2 \log^{5.5} \left( \frac{6nd}{\delta\zeta_-} \right) \sum_{h,r,s} \sum_{\overline{P} \in \mathcal{P}^*(h,r,s)} \left[ \left( \frac{n}{2^t \zeta^2} s^\dagger rh^2 \right) \wedge \left( \sqrt{\frac{n}{2^t h}} + \frac{n}{2^t} \right) \right] \\ &\lesssim \zeta^2 \log^{5.5} \left( \frac{6nd}{\delta\zeta_-} \right) \sum_{h,r,s} \left[ \frac{nsrh^2}{\zeta^2} \wedge \sqrt{\frac{dn}{srh^2}} + n \right] \\ &\lesssim \zeta^2 \log^{8.5} \left( \frac{6nd}{\delta\zeta_-} \right) \left[ \frac{n^{2/3} d^{1/3}}{\zeta^{2/3}} + n \right], \end{align*} $$
where in (a), we combined Corollary D.8 with the fact that the size of each group is at most $n/2^t$, the crude bound $\| [A]_\eta \|_F^2 \le \eta^2 d_1 d_2$ for any $d_1 \times d_2$ matrix and that $n/2^t \ge 1$. In (b), we relied on Lemma D.9, whereas in (c) we used that $x \wedge y \le x^{1/3} y^{2/3}$.

We have proved the desired $n^{2/3} d^{1/3} / \zeta^{2/3} + n$ upper bound. The rate $nd^{1/6} / \zeta^{1/3}$ is proved using the same scheme except that we apply Corollary D.8 differently in (a). More precisely, we have

$$
\begin{align*}
\sum_{\overline{P} \in \overline{\mathcal{L}}_t} \|M(\overline{P}) - \overline{M}(\overline{P})\|_F^2 &\le 16\zeta^2 |\overline{\mathcal{L}}_t| + 96 |\mathcal{R}| \|\mathcal{H}\| \sum_{h,r,s} \sum_{\overline{P} \in \mathcal{P}^*(h,r,s)} \| [\Theta(\overline{P}, Q_{\mathrm{cp}}^\dagger) - \overline{\Theta}(\overline{P}, Q_{\mathrm{cp}}^\dagger)]_{\sqrt{rh}} \|_F^2 \\
&\le \zeta^2 \log^{5.5} \left(\frac{6nd}{\delta\zeta_-}\right) \sum_{h,r,s} \sum_{\overline{P} \in \mathcal{P}^*(h,r,s)} \left[ \left(\frac{n}{2t\zeta}h\sqrt{d}\right) \wedge \left(\sqrt{\frac{n}{2t^h}}h\sqrt{d}\right) \wedge \sqrt{\frac{n}{2^t}d} + \frac{n}{2^t} \right] \\
&\lesssim \zeta^2 \log^{5.5} \left(\frac{6nd}{\delta\zeta_-}\right) \sum_h \sum_{r,s} |\mathcal{P}^*(h,r,s)| \left[ \left(\frac{n}{2t\zeta}h\sqrt{d}\right) \wedge \left(\sqrt{\frac{n}{2t^h}}h\sqrt{d}\right) \wedge \sqrt{\frac{n}{2^t}d} + \frac{n}{2^t} \right] \\
&\stackrel{(a')}{\lesssim} \zeta^2 \log^{5.5} \left(\frac{6nd}{\delta\zeta_-}\right) \sum_h \left[ \frac{nh\sqrt{d}}{\zeta} \wedge \sqrt{\frac{n^2}{h}} \wedge \sqrt{n2^td} + n \right] \\
&\stackrel{(b')}{\lesssim} \zeta^2 \log^{6.5} \left(\frac{6nd}{\delta\zeta_-}\right) \left[ \frac{nd^{1/6}}{\zeta^{1/3}} \wedge n\sqrt{d} + n \right],
\end{align*}
$$

where in (a'), we used that $\sum_{r,s} |\mathcal{P}^*(h,r,s)| \le 2^t \le 2n$ and in (b') that $xy \le x^{1/3} y^{2/3}$. $\square$

*Proof of Lemma D.9.* To ease the notation, we write $\mathcal{P}^* = \mathcal{P}^*(h,r,s)$ in this proof. Since $\mathcal{P}^* \subset \overline{\mathcal{L}}_t$, we straightforwardly derive that $|\mathcal{P}^*| \le |\overline{\mathcal{L}}_t| \le 2^t$. Let us introduce the width of a matrix $\Theta \in \mathbb{R}^{[n] \times Q_r}$ on the set $P \subset [n]$ and $Q \subset Q_r$:

$$W_{\infty,1}(\Theta, P, Q) = \max_{i,j \in P} \sum_{l \in Q} |\Theta_{i,l} - \Theta_{j,l}|.$$

Consider any set $\overline{P}$ and the corresponding quantities $\overline{P}^\dagger$, $s^\dagger$, $r^\dagger$, and $h^\dagger$. By definition of $s^\dagger$, we have $W_{\infty,1}(\Theta, \overline{P}^\dagger, Q_{\mathrm{cp}}^\dagger) \ge s^\dagger \sqrt{r^\dagger h^\dagger}$. Recall that the matrix $\Theta$ is, up to a permutation of its rows, bi-isotonic. Besides, all the groups $\overline{P}$ in $\overline{\mathcal{L}}_t$ are perfectly ordered by assumption. As a consequence, the width of $\Theta$ on $[n]$ is larger or equal to the sum of the width on each set $\overline{P}$. Since $\mathcal{P}^*$ is an ordered sub-partition, it holds that

$$W_{\infty,1}(\Theta, [n], Q_r) \ge \sum_{\overline{P} \in \mathcal{P}^*} W_{\infty,1}(\Theta, \overline{P}^\dagger, Q_r) \ge \sum_{\overline{P} \in \mathcal{P}^*} W_{\infty,1}(\Theta, \overline{P}^\dagger, Q_{\mathrm{cp}}^\dagger) \ge |\mathcal{P}^*|s^\dagger\sqrt{rh} \ge |\mathcal{P}^*|s\sqrt{rh}. \quad (68)$$

By definition of $Q_r$, we have $|Q_r| \le 2d/r$. Since the values of $\Theta$ lie in $[0, \sqrt{r}]$, we deduce that $W_{\infty,1}(\Theta, [n], Q_r) \le 2d/\sqrt{r}$. Together with (68), this yields

$$|\mathcal{P}^*| \le \frac{2d}{rsh},$$

which concludes the proof. $\square$

# Appendix E: Remaining proofs for Proposition B.5

## E.1. Proofs of the results on DimensionReduction (Appendix D.2)

*Proof of lemma D.1.* It is sufficient to prove that $D_{\mathrm{cp}}^*(h,r) \subset \tilde{D}_{\mathrm{cp}}(h,r) \subset \overline{D}_{\mathrm{cp}}^*(h,r)$. Recall that we use the convention that $\bar{y}_i = \bar{m}_i = 0$ if $i \le 0$ and $\bar{y}_i = \bar{m}_i = 1$ if $i > d$. Since the CUSUM
statistic is linear, we have the decomposition

$$ \hat{\mathbf{C}}_{k, \tilde{r}} = \mathbf{C}_{k, \tilde{r}}^* + \frac{1}{\tilde{r}} \left( \sum_{k'=k}^{k+\tilde{r}-1} \bar{e}_{k'}(\bar{P}) - \sum_{k'=k-\tilde{r}}^{k-1} \bar{e}_{k'}(\bar{P}) \right), $$

where the latter random variable is centered and $\zeta(|\bar{P}|\tilde{r}/2)^{-1/2}$-subGaussian. By a union bound, we derive that, with probability higher than $1 - \delta$, we have

$$ \max_{k \in [d]} |\hat{\mathbf{C}}_{k, \tilde{r}} - \mathbf{C}_{k, \tilde{r}}^*| \le \zeta \sqrt{\frac{2 \cdot 2}{|\bar{P}|} \log \left( \frac{2d}{\delta} \right)}. $$

Since $\tilde{r}$ is defined in such a way that

$$ \zeta \sqrt{\frac{4}{|\bar{P}|} \log \left( \frac{2d}{\delta} \right)} \le \frac{1}{8} \sqrt{\tilde{r}h}, $$

we deduce that $\hat{D}_{cp}(h, r) \subset \bar{D}_{cp}^*(h, r)$. Conversely, if $k$ belongs to $D_{cp}^*(h, r)$, we have $\mathbf{C}_{k, 8r}^* \ge h/2$. Since $\bar{m}(\bar{P})$ is an isotonic vector and $\tilde{r} \ge 8r$, it follows that $\mathbf{C}_{k, \tilde{r}}^* \ge \mathbf{C}_{k, 8r}^* \ge h/2$. We deduce that

$$ \hat{\mathbf{C}}_{k, \tilde{r}} \ge h \left[ \frac{1}{2} - \frac{1}{8} \right] \ge \frac{h}{4}, $$

which implies that $D_{cp}^*(h, r) \subset \hat{D}_{cp}(h, r)$. $\square$

*Proof of Lemma D.2.* If an index $k$ belongs to $\bar{D}^*(h, r)$, this implies that $\bar{m}_{k+\tilde{r}} - \bar{m}_{k-\tilde{r}} \ge h/8$, since the vector $\bar{m}$ is isotonic. Define $\kappa = 1 + [\tilde{r}/r]$. Since $\bar{m}$ is an isotonic vector, for $l \in \bar{Q}_{cp}^*(h, r)$, we deduce that $\bar{m}_{l+\kappa r} - \bar{m}_{l-\kappa r} \ge h/8$. Consider the regular grid $Q_{\kappa r}$ of width $\kappa r$ and define $\bar{Q}^*(h, r, \kappa) = \{l \in Q_{\kappa r} : \bar{Q}_{cp}^*(h, r) \cap [l, l + \kappa r) \ne \emptyset\}$. Since, for $l \in \bar{Q}^*(h, r, \kappa)$, we have $\bar{m}_{l+2\kappa r} - \bar{m}_{l-2\kappa r} \ge h/8$ and since the total variation of $\bar{m}$ is at most one, this implies

$$ \frac{h}{8} |\bar{Q}^*(h, r, \kappa)| \le \sum_{l \in D(\kappa, r, h)} \bar{m}_{l+2\kappa r} - \bar{m}_{l-2\kappa r} \le \sum_{l \in Q_{\kappa r}} \bar{m}_{l+2\kappa r} - \bar{m}_{l-2\kappa r} \le 4. $$

Since $|\bar{Q}^*(h, r)| \le \kappa |\bar{Q}^*(h, r, \kappa)|$ and since $\kappa \le 2\tilde{r}/r$, we obtain the desired result. $\square$

*Proof of Lemma D.3.* For any height $h \in \mathcal{H}$-recall the definition of the dyadic class $\mathcal{H}$ in (36)– and any expert $i \in \bar{P}$, we consider the $h$-level set $M_{i, \cdot} - \bar{m}$, that is

$$ F(i, h) = \{k \in [d] : M_{i,k} - \bar{m}_k \ge h\}; \quad F(i, -h) = \{k \in [d] : M_{i,k} - \bar{m}_k \le -h\}. \quad (69) $$

Since $F(i, h)$ and $F(i, -h)$ are subsets of $[d]$, we can decompose them into unions of disjoint intervals. For any positive integer $r \in \mathcal{R}$, we write $F(i, h, r)$ as the union of intervals of $F(i, h)$ whose size belongs $[2r - 1, 4r - 1)$. Finally, we consider the subset $F(i, h, r; 2h) \subset F(i, h, r)$ of all intervals of $F(i, h, r)$ that intersect $F(i, 2h)$. In other words, any maximal interval $I$ in $F(i, h, r; 2h)$ is a $h$-level set whose size belongs to $[2r - 1, 4r - 1)$ and such that $M_{i, \cdot} - \bar{m}$ crosses the level $2h$ in $I$. We define similarly $F(i, h, r)$ and $F(i, h, r; 2h)$ when $h$ is negative and $-h \in \mathcal{H}$. It follows from these definitions that, for any $h$ such that either $h \in \mathcal{H}$ or $-h \in \mathcal{H}$, we have

$$ F(i, 2h) \subset \bigcup_{r \in \mathcal{R}} F(i, h, r; 2h). \quad (70) $$
We define $F^*(h, r, 2h)$ as the union of those intervals for $i \in \bar{P}$.

$$F^*(h, r, 2h) = \bigcup_{i \in \bar{P}} F(i, h, r; 2h) .$$

First, we claim that this collection of intervals $F^*(h, r, 2h)$ is contained in the significant regions of variation of $\bar{m}$. This result heavily relies on the monotonicity assumptions.

**Lemma E.1.** For any $h \in \mathcal{H}$ and any $r \in \mathcal{R}$.

$$[F^*(h, r, 2h) \cup F^*(-h, r, -2h)] \subset D_{cp}^*(h, r) .$$

Next, we quantify $\|M(\bar{P}) - \bar{M}(\bar{P})\|_F^2$ using regions of large variation of $M_{i, \cdot} - \bar{m}$.

**Lemma E.2.** For any $\bar{P}$, it holds that

$$\|M(\bar{P}) - \bar{M}(\bar{P})\|_F^2 \le 16 \left[ \zeta^2 + \sum_{i \in \bar{P}} \sum_{r \in \mathcal{R}, h \in \mathcal{H}} h^2 (|F(i, h, r; 2h)| + |F(i, -h, r; -2h)|) \right] .$$

The last lemma connects these sets $|F(i, h, r; 2h)|$ to the norm of the thresholded aggregated matrix.

**Lemma E.3.** For any $r \in \mathcal{R}$ and $h \in \mathcal{H}$, we consider $\Theta(\bar{P}, Q_{cp}^*(h, r))$ the aggregation of $M$ at scale $r$ and at $Q_{cp}^*(h, r)$. We have

$$h^2 \sum_{i \in \bar{P}} [|F(i, h, r; 2h)| + |F(i, -h, r; -2h)|] \le 3 \| [\Theta(\bar{P}, Q_{cp}^*(h, r)) - \bar{\Theta}(\bar{P}, Q_{cp}^*(h, r))]_{\sqrt{rh}} \|_F^2 .$$

Combining Lemmas E.2 and E.3, we conclude that

$$\begin{aligned} \|M(\bar{P}) - \bar{M}(\bar{P})\|_F^2 &\le 16 \left[ \zeta^2 + \sum_{i \in \bar{P}} \sum_{r \in \mathcal{R}, h \in \mathcal{H}} h^2 (|F(i, h, r; 2h)| + |F(i, -h, r; -2h)|) \right] \\ &\le 16 \left[ \zeta^2 + 3 \sum_{r \in \mathcal{R}, h \in \mathcal{H}} \| [\Theta(\bar{P}, Q_{cp}^*(h, r)) - \bar{\Theta}(\bar{P}, Q_{cp}^*(h, r))]_{\sqrt{rh}} \|_F^2 \right] \\ &\le 16 \left[ \zeta^2 + 6 |\mathcal{R}| \max_{r \in \mathcal{R}, h \in \mathcal{H}} \| [\Theta(\bar{P}, Q_{cp}^*(h, r)) - \bar{\Theta}(\bar{P}, Q_{cp}^*(h, r))]_{\sqrt{rh}} \|_F^2 \right], \end{aligned}$$

which concludes the proof of Lemma D.3.

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

*Proof of Lemma E.1.* Consider any $i \in \bar{P}$, any height $h \in \mathcal{H}$, and any scale $r \in \mathcal{R}$. Without loss of generality, we only focus on $F(i, h, r; 2h)$; the case of $F(i, -h, r; -2h)$ being analogous. Let $I$ be an interval of $F(i, h, r; 2h)$. Fix any question $k \in I$ such that $|M_{i,k} - \bar{m}_k| \ge 2h$. Since $k \in F(i, h, r)$, it follows that there exists $l < 4r$ such that $M_{i,k+l} - \bar{m}_{k+l} \le h$. Since both the vectors $M_{i, \cdot}$ and $\bar{m}$ are isotonic, it follows that $\bar{m}_{k+l} - \bar{m}_k \ge h$. Now consider any $k_0 \in I$. Using again the monotonicity of $\bar{m}$, we deduce that,

$$C_{k_0, 8r}^* \ge \frac{4rh}{8r} \ge \frac{1}{2}h ,$$

and $k_0$ therefore belongs to $D_{cp}^*(h, r)$. We have proved the desired result.

<figure><img src="image_10.png" /></figure>
Proof of Lemma E.2. Consider any expert $i \in \bar{P}$. We decompose the norm of $[M_{i,} - \bar{m}(\bar{P})]$ using the level sets of this vector. We recall that $\mathcal{H}$ is of the form $\{h_{\min}, 2h_{\min}, 4h_{\min}, \dots\}$ where $h_{\min} \in [\zeta^2/nd, 2\zeta^2/nd]$.

$$
\begin{aligned}
\|M_{i,} - \bar{m}(\bar{P})\|_2^2 &\le \sum_{h \in \mathcal{H}} \sum_{k=1}^d (M_{i,k} - \bar{m}_k(\bar{P}))^2 \mathbf{1}\{2h \le |M_{i,k} - \bar{m}_k(\bar{P})| < 4h\} + 4dh_{\min}^2 \\
&\le 16 \sum_{h \in \mathcal{H}} \sum_{k=1}^d h^2 \mathbf{1}\{2h \le |M_{i,k} - \bar{m}_k(\bar{P})| < 4h\} + \frac{16\zeta^2}{n^2d} \\
&\le 16 \sum_{h \in \mathcal{H}} h^2 |F(i, 2h)| + \frac{16\zeta^2}{n^2d} \\
&\le 16 \sum_{h \in \mathcal{H}} \sum_{r \in \mathcal{R}} h^2 |F(i, h, r; 2h)| + \frac{16\zeta^2}{n^2d},
\end{aligned}
$$

where in the last line, we used (70). Then, we sum over $i \in \bar{P}$ to conclude. $\square$

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

Proof of Lemma E.3. Consider any $i \in \bar{P}$, any height $h \in \mathcal{H}$, and any scale $r \in \mathcal{R}$. Without loss of generality, we only consider $F(i, h, r; 2h)$ the case of $F(i, -h, r; -2h)$ being analogous. Let $I$ be a maximal interval of $F(i, h, r; 2h)$. We deduce from Lemma E.1 that $I$ is included in $D_{\mathrm{cp}}^*(h, r)$. Let $I_0$ be the largest sub-interval of $I$ of the form $[qr, q'r)$ where $q$ and $q' \in \mathcal{Q}_r$. Since $|I| \ge 2r-1$, it follows that $|I| \le 3|I_0|$. We write $L_0$ the subset of columns of the aggregated matrix $\Theta(\bar{P}, Q_{\mathrm{cp}}^*(h, r))$ corresponding to $I_0$ so that $|L_0| = |I_0|/r$. On each column $l$ of $L_0$, we have $\Theta_{i,l}(\bar{P}, Q_{\mathrm{cp}}^*(h, r)) - \bar{\theta}_l(\bar{P}, Q_{\mathrm{cp}}^*(h, r)) \ge \sqrt{rh}$. Putting everything together, we get

$$h^2|I| \le 3h^2|I_0| = 3h^2 r |L_0| \le 3 \sum_{l \in L_0} \left( \left[ \Theta(\bar{P}, Q_{\mathrm{cp}}^*(h, r))_{i,l} - \bar{\theta}(\bar{P}, Q_{\mathrm{cp}}^*(h, r))_l \right]_{\sqrt{rh}} \right)^2 .$$

Summing over all intervals $I$ and over all experts $i \in \bar{P}$ and also accounting for the $F[i, -h, r; -2h]$ concludes the proof. $\square$

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

E.2. Proof of Proposition D.4

To simplify the notation, we define $\Phi_{l_1} = 2(2\zeta\sqrt{2} + \bar{\beta}_{\mathrm{tris}})\sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)} = \phi_{l_1}\sqrt{\log\left(\frac{2|\bar{P}|}{\delta}\right)}$.

For simplicity, we respectively write $\Theta(\bar{P}) = \Theta(\bar{P}, Q)$ and $\Theta(\bar{P}') = \Theta(\bar{P}', Q)$ in this proof. Recall that $\bar{\theta}(\bar{P})$ stands the mean row of $\Theta(\bar{P})$ whereas $\bar{\theta}(\bar{P}')$ stands for the mean row of $\Theta(\bar{P}')$.

Invoking lemma C.1 with $w = \mathbf{1}_Q$ and since the matrix $\Theta$ is isotonic, we deduce that outside an event of probability smaller than $\delta$, we have

$$\max_{i,j \in \bar{P}'} \|\Theta(\bar{P}')_{i,\cdot} - \Theta(\bar{P}')_{j,\cdot}\|_1 \le \Phi_{l_1} \sqrt{|Q|}. \quad (71)$$

since the matrix $\Theta$ is isotonic. We shall deduce from this inequality the desired bound. We consider two cases depending on the difference between $\bar{\theta}(\bar{P})$ and $\bar{\theta}(\bar{P}')$ the mean rows in $\bar{P}$ and $\bar{P}'$.

**Case 1:** $|\bar{P}'| \cdot \|\bar{\theta}(\bar{P}) - \bar{\theta}(\bar{P}')\|_2^2 > \frac{1}{16\phi} \|\Theta(\bar{P}) - \bar{\Theta}(\bar{P})\|_F^2$. Since $\bar{P}' \subset \bar{P}$, we deduce that

$$
\begin{aligned}
\|\Theta(\bar{P}) - \bar{\Theta}(\bar{P})\|_F^2 - \|\Theta(\bar{P}') - \bar{\Theta}(\bar{P}')\|_F^2 &\ge \sum_{i \in \bar{P}'} \|\Theta(\bar{P})_{i,\cdot} - \bar{\theta}(\bar{P})\|_2^2 - \|\Theta(\bar{P})_{i,\cdot} - \bar{\theta}(\bar{P}')\|_2^2 \\
&= |\bar{P}'| \cdot \|\bar{\theta}(\bar{P}) - \bar{\theta}(\bar{P}')\|_2^2 \\
&\ge \frac{1}{16\phi} \|\Theta(\bar{P}) - \bar{\Theta}(\bar{P})\|_F^2,
\end{aligned}
$$
where we used the condition in the last line. We have proved the desired result.

**Case 2:** $|\overline{P}'| \cdot \|\overline{\theta}(\overline{P}) - \overline{\theta}(\overline{P}')\|_2^2 \le \frac{1}{16\phi} \|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2$. We start with the decomposition

$$ \|\Theta(\overline{P}') - \overline{\Theta}(\overline{P}')\|_F^2 \le \|\Theta(\overline{P}') - \overline{\Theta}(\overline{P})\|_F^2 = \|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2 - \|\Theta(\overline{P} \setminus \overline{P}') - \overline{\Theta}(\overline{P})\|_F^2, \quad (72) $$

so that we only have to control $\|\Theta(\overline{P} \setminus \overline{P}') - \overline{\Theta}(\overline{P})\|_F^2$ from below. By definition of the operator $[\cdot]_\eta$, we have

$$ \|\Theta(\overline{P} \setminus \overline{P}') - \overline{\Theta}(\overline{P})\|_F^2 \ge \|[\Theta(\overline{P} \setminus \overline{P}') - \overline{\Theta}(\overline{P})]_\eta\|_F^2 = \|[\Theta(\overline{P}) - \overline{\Theta}(\overline{P})]_\eta\|_F^2 - \|[\Theta(\overline{P}') - \overline{\Theta}(\overline{P})]_\eta\|_F^2. $$

By assumption, we have $\|[\Theta(\overline{P}) - \overline{\Theta}(\overline{P})]_\eta\|_F^2 \ge \phi^{-1} \|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2$. Hence, as long as we prove that

$$ \|[\Theta(\overline{P}') - \overline{\Theta}(\overline{P})]_\eta\|_F^2 \le (2\phi)^{-1} \|[\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2, \quad (73) $$

we can safely conclude from (72) that

$$ \|\Theta(\overline{P}') - \overline{\Theta}(\overline{P}')\|_F^2 \le (1 - (2\phi)^{-1}) \|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2. $$

Thus, we only have to prove (73). Again, by definition of the thresholding operator, we have

$$ \begin{aligned} \|[\Theta(\overline{P}') - \overline{\Theta}(\overline{P})]_\eta\|_F^2 &= \eta^2 \sum_{i \in \overline{P}', l \in Q} \mathbf{1}_{|\Theta_{i,l} - \overline{\theta}(\overline{P})_l| \ge \eta} \\ &\le \eta^2 \sum_{i \in \overline{P}', l \in Q} \mathbf{1}_{|\Theta_{i,l} - \overline{\theta}(\overline{P})_l| \ge \eta/2} + \mathbf{1}_{|\theta(\overline{P}')_l - \overline{\theta}(\overline{P})_l| \ge \eta/2}. \end{aligned} \quad (74) $$

By Markov inequality, the condition that defines Case 2 above implies that

$$ \sum_{l \in Q} \mathbf{1}\{|\overline{\theta}(\overline{P}')_l - \overline{\theta}(\overline{P})_l| \ge \eta/2\} \le \frac{1}{4\phi} \frac{\|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2}{|\overline{P}'|\eta^2}. \quad (75) $$

From (71) and a convexity argument, we deduce that, for any $i \in \overline{P}'$, $\|\Theta(\overline{P}')_{i,\cdot} - \overline{\theta}(\overline{P}')\|_1 \le \Phi_{l_1} \sqrt{|\overline{Q}|}$. Then, applying again Markov inequality, we deduce that, for any expert $i$ in $\overline{P}'$ and any $\eta > 0$, we have

$$ \sum_{l \in Q} \mathbf{1}\{|\Theta(\overline{P}')_{i,l} - \overline{\theta}(\overline{P})_l| \ge \eta/2\} \le 2\Phi_{l_1} \frac{\sqrt{|\overline{Q}|}}{\eta}. $$

Since we assume that $\|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2 \ge 8\phi\Phi_{l_1}\eta|\overline{P}|\sqrt{|\overline{Q}|} \ge 8\phi\Phi_{l_1}\eta|\overline{P}'|\sqrt{|\overline{Q}|}$, we deduce that

$$ \sum_{l \in Q} \mathbf{1}\{|\Theta(\overline{P})_{i,l} - \overline{\theta}(\overline{P})_l| \ge \eta/2\} \le \frac{1}{4\phi} \frac{\|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2}{|\overline{P}'|\eta^2}. \quad (76) $$

So that, combining (74), (75) and (76), we arrive at

$$ \|[\Theta(\overline{P}') - \overline{\Theta}(\overline{P})]_\eta\|_F^2 \le \frac{1}{2\phi} \|\Theta(\overline{P}) - \overline{\Theta}(\overline{P})\|_F^2. $$

We have proved (73).