Breaking the $1/\sqrt{n}$ Barrier: Faster Rates for Permutation-based Models in Polynomial Time

TL;DR

This paper introduces a polynomial-time algorithm for estimating permutation-based isotonic matrices that achieves faster convergence rates, narrowing the gap between statistical and computational efficiency in matrix estimation.

Contribution

The authors develop a new algorithm that improves the estimation rate for permutation-based models, achieving a rate of O(n^{-3/4}) in polynomial time, surpassing previous methods.

Findings

Achieves O(n^{-3/4}) rate in Frobenius norm

Efficient polynomial-time algorithm for permutation-based matrix estimation

Narrowing the gap between statistical and computational rates

Abstract

Many applications, including rank aggregation and crowd-labeling, can be modeled in terms of a bivariate isotonic matrix with unknown permutations acting on its rows and columns. We consider the problem of estimating such a matrix based on noisy observations of a subset of its entries, and design and analyze a polynomial-time algorithm that improves upon the state of the art. In particular, our results imply that any such $n \times n$ matrix can be estimated efficiently in the normalized Frobenius norm at rate $\widetilde{\mathcal O}(n^{-3/4})$, thus narrowing the gap between $\widetilde{\mathcal O}(n^{-1})$ and $\widetilde{\mathcal O}(n^{-1/2})$, which were hitherto the rates of the most statistically and computationally efficient methods, respectively.

Figures

Equations

$$\mathbb{E}\left[\frac{1}{n^2}\Vert \hat{M}-M^{\ast }{\Vert }_F^2\right]\le C{\left(\frac{\log n}{n}\right)}^{3/4}.$$

$${\left[M(\pi ,\sigma )\right]}_{i,j}=M_{\pi (i),\sigma (j)}\text{ for }i\in [n_1],j\in [n_2].$$

$${\mathbb{C}}_{\mathsf{Perm}}:=\bigcup _{\begin{array}{c}\pi \in {\mathfrak{S}}_{n_1}\\ \sigma \in {\mathfrak{S}}_{n_2}\end{array}}{\mathbb{C}}_{\mathsf{BISO}}(\pi ,\sigma ).$$

$$y_{\ell }=\mathrm{tr}(X_{\ell }^{\top }{M}^{\ast })+z_{\ell },$$

$$\displaystyle y_{\ell}\sim\mbox{\sf Ber}\big{(}\operatorname{\mathsf{tr}}(X_{\ell}^{\top}M)\big{)}$$

$$Y_{i,j}=\frac{1}{p_{\mathsf{obs}}}\frac{1}{1\vee {\sum }_{\ell =1}^{N^{\prime }}1\{X_{\ell }=E^{(i,j)}\}}\sum _{\ell =1}^{N^{\prime }}{y}_{\ell }1\{X_{\ell }=E^{(i,j)}\}.$$

$${\hat{M}}_{\mathsf{LS}}(Y):=\arg \underset{M\in {\mathbb{C}}_{\mathsf{Perm}}}{\min }\Vert Y-M{\Vert }_F^2.$$

$$\frac{1}{n_1{n}_2}\Vert {\hat{M}}_{\mathsf{LS}}(Y)-M^{\ast }{\Vert }_F^2\lesssim ({\zeta }^2\vee 1)\frac{n_1{\log }^2{n}_1}{N}$$

$$\underset{M^{\ast }\in {\mathbb{C}}_{\mathsf{Perm}}}{\sup }\mathbb{E}\left[\frac{1}{n_1{n}_2}\Vert \hat{M}-M^{\ast }{\Vert }_F^2\right]\gtrsim \frac{n_1}{N}.$$

$$\displaystyle\frac{1}{n_{1}n_{2}}\big{\|}\widehat{M}(\widehat{\pi},\widehat{\sigma})-M^{*}\big{\|}_{F}^{2}\lesssim(\zeta^{2}\lor 1)\frac{n_{1}\log^{2}n_{1}}{N}$$

$$\displaystyle+\frac{1}{n_{1}n_{2}}\big{\|}M^{*}(\mathsf{id},\sigma^{-1}\circ\widehat{\sigma})-M^{*}\big{\|}_{F}^{2}.$$

$$C(j)=\sum _{i=1}^{n_1}{Y}_{i,j}.$$

$${\mathsf{bl}}_1$$

$${\mathsf{bl}}_k$$

$${\mathsf{bl}}_K$$

$$S(i)=\sum _{j\in [n_2]}{Y}_{i,j}^{(2)}\text{ for each }i\in [n_1],$$

$$S_{{\mathsf{BL}}_k}(i)=\sum _{j\in {\mathsf{BL}}_k}{Y}_{i,j}^{(2)}\text{ for each }i\in [n_1],k\in [K].$$

$$S(v)-S(u)$$

$$S_{{\mathsf{BL}}_k}(v)-S_{{\mathsf{BL}}_k}(u)$$

$$\displaystyle\frac{1}{n_{1}n_{2}}\big{\|}\widehat{M}(\widehat{\pi}_{\sf tds},\widehat{\sigma}_{\sf tds})-M^{*}\big{\|}_{F}^{2}\lesssim(\zeta^{2}\lor 1)\left[\Big{(}\frac{n_{1}\log n_{1}}{N}\Big{)}^{3/4}+\frac{n_{1}\log^{2}n_{1}}{N}\right]$$

$$\displaystyle\frac{1}{n_{1}n_{2}}\big{\|}\widehat{M}(\widehat{\pi}_{\sf tds},\widehat{\sigma}_{\sf tds})-M^{*}\big{\|}_{F}^{2}=\widetilde{O}\left(n_{2}^{-3/4}\right),$$

$${\mathbb{C}}_{\mathsf{SST}}(n):=\left(\bigcup _{\pi \in {\mathfrak{S}}_n}{\mathbb{C}}_{\mathsf{BISO}}^{\prime }(\pi ,\pi )\right)\bigcap {\mathbb{C}}_{\mathsf{skew}}(n).$$

$$\mathrm{Pr}(i\succ k)\ge \max \{\mathrm{Pr}(i\succ j),\mathrm{Pr}(j\succ k)\}.$$

$$\displaystyle\Pr\left\{\Big{\|}\widehat{M}_{{\sf LS}}(Y)-M^{*}\Big{\|}_{F}^{2}\geq\frac{c_{3}}{p_{\mathsf{obs}}}(\zeta^{2}\lor 1)n_{1}\log^{2}n_{1}\right\}\leq(n_{1}n_{2})^{-3}.$$

$$\displaystyle\Pr\left\{\bigg{|}\sum_{(i,j)\in\mathcal{S}}(Y_{i,j}-M^{*}_{i,j})\bigg{|}\geq 8(\zeta+1)\bigg{(}\sqrt{\frac{|\mathcal{S}|n_{1}n_{2}}{N}\log(n_{1}n_{2})}+2\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})\bigg{)}\right\}\leq 2(n_{1}n_{2})^{-4}.$$

$$W_{i,j}=Y_{i,j}-M_{i,j}^{\ast }=\{\begin{array}{l}-M_{i,j}^{\ast }\text{ if entry (i,j) is not observed, and}\\ M_{i,j}^{\ast }/p_{\mathsf{obs}}-M_{i,j}^{\ast }+\frac{W_{i,j}^{\prime }}{p_{\mathsf{obs}}},\text{ otherwise,}\end{array}$$

$$Z_{i,j}^{(1)}=\{\begin{array}{ll}{M}_{i,j}^{\ast }(p_{\mathsf{obs}}^{-1}-1)& \text{ with probability }{p}_{\mathsf{obs}},\\ -M_{i,j}^{\ast }& \text{ with probability }1-p_{\mathsf{obs}},\end{array}$$

$$Z_{i,j}^{(2)}\text{ is }\{\begin{array}{ll}\frac{\zeta }{p_{\mathsf{obs}}}\text{-sub-Gaussian}& \text{ with probability }{p}_{\mathsf{obs}},\\ 0& \text{ with probability }1-p_{\mathsf{obs}}.\end{array}$$

$$\displaystyle\Pr\left\{\bigg{|}\sum_{(i,j)\in\mathcal{S}}Z^{(1)}_{i,j}\bigg{|}\geq t\right\}\leq 2\exp\Big{(}-\frac{t^{2}/2}{|\mathcal{S}|/p_{\mathsf{obs}}+t/(3p_{\mathsf{obs}})}\Big{)}.$$

$$\displaystyle\Pr\left\{\bigg{|}\sum_{(i,j)\in\mathcal{S}}Z^{(1)}_{i,j}\bigg{|}\geq 4\sqrt{\frac{|\mathcal{S}|n_{1}n_{2}}{N}\log(n_{1}n_{2})}+6\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})\right\}\leq(n_{1}n_{2})^{-4}.$$

$$\displaystyle\Pr\left\{\bigg{|}\sum_{(i,j)\in\mathcal{S}}Z^{(2)}_{i,j}\bigg{|}\geq\sqrt{\frac{16\zeta^{2}|\mathcal{S}|n_{1}n_{2}}{N}t}+\frac{4\zeta n_{1}n_{2}}{N}t\right\}\leq 2\exp(-t),$$

$$\displaystyle\Pr\left\{\bigg{|}\sum_{(i,j)\in\mathcal{S}}Z^{(2)}_{i,j}\bigg{|}\geq 8\zeta\sqrt{\frac{|\mathcal{S}|n_{1}n_{2}}{N}\log(n_{1}n_{2})}+16\zeta\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})\right\}\leq(n_{1}n_{2})^{-4}.$$

Full text

1 Introduction

Structured††Accepted for presentation at Conference on Learning Theory (COLT) 2018 matrices with entries in the range $[0,1]$ and unknown permutations acting on their rows and columns arise in multiple applications, including estimation from pairwise comparisons [BT52, SBGW17] and crowd-labeling [DS79, SBW16b]. Traditional parametric models [BT52, Luc59, Thu27, DS79] assume that these matrices are obtained from rank-one matrices via a known link function. Aided by tools such as maximum likelihood estimation and spectral methods, researchers have made significant progress in studying both statistical and computational aspects of these parametric models [HOX14, RA14, SBB*+*16, NOS16, ZCZJ16, GZ13, GLZ16, KOS11b, LPI12, DDKR13, GKM11] and their low-rank generalizations [RA16, NOTX17, KOS11a].

There has been evidence from empirical studies (e.g., [ML65, BW97]) that real-world data is not always well-captured by such parametric models. With the goal of increasing model flexibility, a recent line of work has studied the class of permutation-based models [Cha15, SBGW17, SBW16b]. Rather than imposing parametric conditions on the matrix entries, these models impose only shape constraints on the matrix, such as monotonicity, before unknown permutations act on the its rows and columns. This more flexible class reduces modeling bias compared to its parametric counterparts while, perhaps surprisingly, producing models that can be estimated at rates that differ only by logarithmic factors from parametric models. On the negative side, these advantages of permutation-based models are accompanied by significant computational challenges. The unknown permutations make the parameter space highly non-convex, so that efficient maximum likelihood estimation is unlikely. Moreover, spectral methods are often suboptimal in approximating shape-constrained sets of matrices [Cha15, SBGW17]. Consequently, results from many recent papers show a non-trivial statistical-computational gap in estimation rates for models with latent permutations [SBGW17, CM16, SBW16b, FMR16, PWC17].

Related work.

While the main motivation of our work comes from nonparametric methods for aggregating pairwise comparisons, we begin by discussing a few other lines of related work. The current paper lies at the intersection of shape-constrained estimation and latent permutation learning. Shape-constrained estimation has long been a major topic in nonparametric statistics, and of particular relevance to our work is the estimation of a bivariate isotonic matrix without latent permutations [CGS18]. There, it was shown that the minimax rate of estimating an $n\times n$ matrix from noisy observations of all its entries is $\widetilde{\Theta}(n^{-1})$. The upper bound is achieved by the least squares estimator, which is efficiently computable due to the convexity of the parameter space.

Shape-constrained matrices with permuted rows or columns also arise in applications such as seriation [FJBd13, FMR16] and feature matching [CD16]. In particular, the monotone subclass of the statistical seriation model [FMR16] contains $n\times n$ matrices that have increasing columns, and an unknown row permutation. The authors established the minimax rate $\widetilde{\Theta}(n^{-2/3})$ for estimating matrices in this class and proposed a computationally efficient algorithm with rate $\widetilde{\mathcal{O}}(n^{-1/2})$. For the subclass of such matrices where in addition, the rows are also monotone, the results of the current paper improve the two rates to $\widetilde{\mathcal{O}}(n^{-1})$ and $\widetilde{\mathcal{O}}(n^{-3/4})$ respectively.

Another related model is that of noisy sorting [BM08], which involves a latent permutation but no shape-constraint. In this prototype of a permutation-based ranking model, we have an unknown, $n\times n$ matrix with constant upper and lower triangular portions whose rows and columns are acted upon by an unknown permutation. The hardness of recovering any such matrix in noise lies in estimating the unknown permutation. As it turns out, this class of matrices can be estimated efficiently at minimax optimal rate $\widetilde{\Theta}(n^{-1})$ by multiple procedures: the original work by Braverman and Mossel [BM08] proposed an algorithm with time complexity $\mathcal{O}(n^{c})$ for some unknown and large constant $c$, and recently, an $\widetilde{\mathcal{O}}(n^{2})$-time algorithm was proposed by Mao et al. [MWR17]. These algorithms, however, do not generalize beyond the noisy sorting class, which constitutes a small subclass of an interesting class of matrices that we describe next.

The most relevant body of work to the current paper is that on estimating matrices satisfying the strong stochastic transitivity condition, or SST for short. This class of matrices contains all $n\times n$ bivariate isotonic matrices with unknown permutations acting on their rows and columns, with an additional skew-symmetry constraint. The first theoretical study of these matrices was carried out by Chatterjee [Cha15], who showed that a spectral algorithm achieved the rate $\widetilde{\mathcal{O}}(n^{-1/4})$ in the normalized Frobenius norm. Shah et al. [SBGW17] then showed that the minimax rate of estimation is given by $\widetilde{\Theta}(n^{-1})$, and also improved the analysis of the spectral estimator of Chatterjee [Cha15] to obtain the computationally efficient rate $\widetilde{\mathcal{O}}(n^{-1/2})$. In follow-up work [SBW16a], they also showed a second $\mathsf{CRL}$ estimator based on the Borda count that achieved the same rate, but in near-linear time. In related work, Chatterjee and Mukherjee [CM16] analyzed a variant of the $\mathsf{CRL}$ estimator, showing that for subclasses of SST matrices, it achieved rates that were faster than $\mathcal{O}(n^{-1/2})$. In a complementary direction, a superset of the current authors [PMM*+*17] analyzed the estimation problem under an observation model with structured missing data, and showed that for many observation patterns, a variant of the $\mathsf{CRL}$ estimator was minimax optimal.

Shah et al. [SBW16a] also showed that conditioned on the planted clique conjecture, it is impossible to improve upon a certain notion of adaptivity of the $\mathsf{CRL}$ estimator in polynomial time. Such results have prompted various authors [FMR16, SBW16a] to conjecture that a similar statistical-computational gap also exists when estimating SST matrices in the Frobenius norm.

Our contributions.

Our main contribution in the current work is to tighten the aforementioned statistical-computational gap. More precisely, we study the problem of estimating a bivariate isotonic matrix with unknown permutations acting on its rows and columns, given noisy, partial observations of its entries; this matrix class strictly contains the SST model [Cha15, SBGW17] for ranking from pairwise comparisons. As a corollary of our results, we show that when the underlying matrix has dimension $n\times n$ and $\Theta(n^{2})$ noisy entries are observed, our polynomial-time, two-dimensional sorting algorithm provably achieves the rate of estimation $\widetilde{}\mathcal{O}(n^{-3/4})$ in the normalized Frobenius norm; thus, this result breaks the previously mentioned $\widetilde{}\mathcal{O}(n^{-1/2})$ barrier [SBGW17, CM16]. Although the rate $\widetilde{}\mathcal{O}(n^{-3/4})$ still differs from the minimax optimal rate $\widetilde{\Theta}(n^{-1})$, our algorithm is, to the best of our knowledge, the first efficient procedure to obtain a rate faster than $\widetilde{}\mathcal{O}(n^{-1/2})$ uniformly over the SST class. This guarantee, which is stated in slightly more technical terms below, can be significant in practice (see Figure 1).

Main theorem (informal)

There is an estimator $\widehat{M}$ computable in time $\mathcal{O}(n^{2.5})$ such that for any $n\times n$ SST matrix $M^{*}$, given $\Theta(n^{2})$ Bernoulli observations of its entries, we have

[TABLE]

Our algorithm is novel in the sense that it is neither spectral in nature, nor simple variations of the Borda count estimator that was previously employed. Our algorithm takes advantage of the fine monotonicity structure of the underlying matrix along both dimensions, and this allows us to prove tighter bounds than before. In addition to making algorithmic contributions, we also briefly revisit the minimax rates of estimation.

Organization.

In Section 2, we formally introduce our estimation problem. Section 3 contains statements and discussions of our main results, and in Section 4, we describe in detail how the estimation problem that we study is connected to applications in crowd-labeling and ranking from pairwise comparisons. We provide the proofs of our main results in Section 5.

Notation.

For a positive integer $n$, let $[n]:\,=\{1,2,\ldots,n\}$. For a finite set $S$, we use $|S|$ to denote its cardinality. For two sequences $\{a_{n}\}_{n=1}^{\infty}$ and $\{b_{n}\}_{n=1}^{\infty}$, we write $a_{n}\lesssim b_{n}$ if there is a universal constant $C$ such that $a_{n}\leq Cb_{n}$ for all $n\geq 1$. The relation $a_{n}\gtrsim b_{n}$ is defined analogously. We use $c,C,c_{1},c_{2},\dots$ to denote universal constants that may change from line to line. We use $\mbox{\sf Ber}(p)$ to denote the Bernoulli distribution with success probability $p$, the notation $\mbox{\sf Bin}(n,p)$ to denote the binomial distribution with $n$ trials and success probability $p$, and the notation $\operatorname{\mathsf{Poi}}(\lambda)$ to denote the Poisson distribution with parameter $\lambda$. Given a matrix $M\in{}^{n_{1}\times n_{2}}$, its $i$-th row is denoted by $M_{i}$. For a vector $v\in\mathbb{R}^{n}$, define its variation as $\operatorname{\mathsf{var}}(v)=\max_{i}v_{i}-\min_{i}v_{i}$. Let $\mathfrak{S}_{n}$ denote the set of all permutations $\pi:[n]\to[n]$. Let $\mathsf{id}$ denote the identity permutation, where the dimension can be inferred from context.

2 Background and problem setup

In this section, we present the relevant technical background and notation on permutation-based models, and introduce the observation model of interest.

2.1 Matrix models

Our main focus is on designing efficient algorithms for estimating a bivariate isotonic matrix with unknown permutations acting on its rows and columns. Formally, we define $\mathbb{C}_{\mathsf{BISO}}$ to be the class of matrices in $[0,1]^{n_{1}\times n_{2}}$ with nondecreasing rows and nondecreasing columns. For readability, we assume throughout that $n_{1}\geq n_{2}$ unless otherwise stated; our results can be straightforwardly extended to the other case. Given a matrix $M\in{}^{n_{1}\times n_{2}}$ and permutations $\pi\in\mathfrak{S}_{n_{1}}$ and $\sigma\in\mathfrak{S}_{n_{2}}$, we define the matrix $M(\pi,\sigma)\in{}^{n_{1}\times n_{2}}$ by specifying its entries as

[TABLE]

Also define the class $\mathbb{C}_{\mathsf{BISO}}(\pi,\sigma):\,=\{M(\pi,\sigma):M\in\mathbb{C}_{\mathsf{BISO}}\}$ as the set of matrices that are bivariate isotonic when viewed along the row permutation $\pi$ and column permutation $\sigma$, respectively.

The class of matrices that we are interested in estimating contains bivariate isotonic matrices whose rows and columns are both acted upon by unknown permutations:

[TABLE]

2.2 Observation model

In order to study estimation from noisy observations of a matrix $M^{*}$ in the class $\mathbb{C}_{\mathsf{Perm}}$, we suppose that $N$ noisy entries are sampled independently and uniformly with replacement from all entries of $M^{*}$. This sampling model is popular in the matrix completion literature, and is a special case of the trace regression model [NW12, KLT11]. It has also been used in the context of permutation models by Mao et al. [MWR17] to study the noisy sorting class.

More precisely, let $E^{(i,j)}$ denote the $n_{1}\times n_{2}$ matrix with $1$ in the $(i,j)$-th entry and [math] elsewhere, and suppose that $X_{\ell}$ is a random matrix sampled independently and uniformly from the set $\{E^{(i,j)}:i\in[n_{1}],\,j\in[n_{2}]\}$. We observe $N\leq n_{1}n_{2}$ independent pairs $\{(X_{\ell},y_{\ell})\}_{\ell=1}^{N}$ from the model

[TABLE]

where the observations are contaminated by independent, centered, sub-Gaussian noise $z_{\ell}$ with variance parameter $\zeta^{2}$. Of particular interest is the noise model considered in applications such as crowd-labeling and ranking from pairwise comparisons. Here our samples take the form

[TABLE]

and consequently, the sub-Gaussian parameter $\zeta^{2}$ is bounded; for a discussion of other regimes of noise in a related matrix model, see Gao [Gao17].

For analytical convenience, we employ the standard trick of Poissonization, whereby we assume throughout the paper that $N^{\prime}=\operatorname{\mathsf{Poi}}(N)$ random samples are drawn according to the trace regression model (1). Upper and lower bounds derived under this model carry over with loss of constant factors to the model with exactly $N$ samples; for a detailed discussion, see Appendix B.

For notational convenience, denote the probability that an entry of the matrix is observed under Poissonized sampling by $p_{\mathsf{obs}}=1-\exp(-N/n_{1}n_{2})$. Since we assume throughout that $N\leq n_{1}n_{2}$, it can be verified that $\frac{N}{2n_{1}n_{2}}\leq p_{\mathsf{obs}}\leq\frac{N}{n_{1}n_{2}}$.

Now given $N^{\prime}=\operatorname{\mathsf{Poi}}(N)$ observations $\{(X_{\ell},y_{\ell})\}_{\ell=1}^{N^{\prime}}$, let us define the matrix of observations $Y=Y\left(\{(X_{\ell},y_{\ell})\}_{\ell=1}^{N^{\prime}}\right)$, with entry $(i,j)$ given by

[TABLE]

In words, the rescaled entry $p_{\mathsf{obs}}Y_{i,j}$ is the average of all the noisy realizations of $M^{*}_{i,j}$ that we have observed, or zero if the entry goes unobserved. Note that $\mathbb{E}[Y_{i,j}]=\frac{1}{p_{\mathsf{obs}}}M^{*}_{i,j}\cdot p_{\mathsf{obs}}=M^{*}_{i,j},$ so that ${\mathbb{E}}[Y]=M^{*}$. Moreover, we may write the model in the linearized form $Y=M^{*}+W$, where $W$ is a matrix of additive noise having independent, zero-mean, sub-Gaussian entries.

3 Main results

In this section, we present our main results—we begin by briefly revisiting the fundamental limits of estimation, and then introduce our algorithms in Section 3.2. We assume throughout this section that as per the setup, we have $n_{1}\geq n_{2}$ and $N\in[n_{1}n_{2}]$.

3.1 Statistical limits of estimation

We begin by characterizing the fundamental limits of estimation under the trace regression observation model (1) with $N^{\prime}=\operatorname{\mathsf{Poi}}(N)$ observations. We define the least squares estimator over the class of matrices $\mathbb{C}_{\mathsf{Perm}}$ as the projection

[TABLE]

The projection is a non-convex problem, and is unlikely to be computable exactly in polynomial time. However, studying this estimator allows us to establish a baseline that characterizes the best achievable statistical rate. The following theorem characterizes its risk up to a logarithmic factor in the dimension; recall the shorthand $Y=Y\left(\{X_{\ell},y_{\ell}\}_{\ell=1}^{N^{\prime}}\right)$.

Theorem 1**.**

For any matrix $M^{*}\in\mathbb{C}_{\mathsf{Perm}}$, we have

[TABLE]

*Additionally, under the Bernoulli observation model (2), any estimator $\widehat{M}$ satisfies *

[TABLE]

The factor $(\zeta^{2}\lor 1)$ appears in the upper bound instead of the noise variance $\zeta^{2}$ because even if the noise is zero, there are missing entries. The theorem characterizes the minimax rate of estimation for the class $\mathbb{C}_{\mathsf{Perm}}$ up to a logarithmic factor.

3.2 Efficient algorithms

Next, we propose polynomial-time algorithms for estimating the permutations $(\pi,\sigma)$ and the matrix $M^{*}$. Our main algorithm relies on two distinct steps: first, we estimate the unknown permutations, and then project onto the class of matrices that are bivariate isotonic when viewed along the estimated permutations. The formal meta-algorithm is described below.

Algorithm 1 (meta-algorithm)

• •

Step 0: Split the observations into two disjoint parts, each containing $N^{\prime}/2$ observations, and construct the matrices $Y^{(1)}=Y\left(\{X_{\ell},y_{\ell}\}_{\ell=1}^{N^{\prime}/2}\right)$ and $Y^{(2)}=Y\left(\{X_{\ell},y_{\ell}\}_{\ell=N^{\prime}/2+1}^{N^{\prime}}\right)$.

• •

Step 1: Use $Y^{(1)}$ to obtain the permutation estimates $(\widehat{\pi},\widehat{\sigma})$.

• •

Step 2: Return the matrix estimate $\widehat{M}(\widehat{\pi},\widehat{\sigma}):\,=\arg\min_{M\in\mathbb{C}_{\mathsf{BISO}}(\widehat{\pi},\widehat{\sigma})}\|Y^{(2)}-M\|_{F}^{2}.$

Owing to the convexity of the set $\mathbb{C}_{\mathsf{BISO}}(\widehat{\pi},\widehat{\sigma})$, the projection operation in Step 2 of the algorithm can be computed in near linear time [BDPR84, KRS15]. The following result, a slight variant of Proposition 4.2 of Chatterjee and Mukherjee [CM16], allows us to characterize the error rate of any such meta-algorithm as a function of the permutation estimates $(\widehat{\pi},\widehat{\sigma})$.

Proposition 1**.**

Suppose that $M^{*}\in\mathbb{C}_{\mathsf{BISO}}(\pi,\sigma)$ where $\pi$ and $\sigma$ are unknown permutations in $\mathfrak{S}_{n_{1}}$ and $\mathfrak{S}_{n_{2}}$ respectively. Then with probability at least $1-(n_{1}n_{2})^{-3}$, we have

[TABLE]

The first term on the right hand side of the bound (5) corresponds to an estimation error, if the true permutations $\pi$ and $\sigma$ were known a priori, and the latter two terms correspond to an approximation error that we incur as a result of having to estimate these permutations from data. Comparing the bound (5) to the minimax lower bound (4b), we see that up to a logarithmic factor, the first term of the bound (5) is unavoidable, and so we can restrict our attention to obtaining good permutation estimates $(\widehat{\pi},\widehat{\sigma})$. We now present our main permutation estimation procedure that can be plugged into Step 1 of this meta-algorithm.

3.2.1 Two-dimensional sorting

Since the underlying matrix of interest is individually monotonic along each dimension, the row and column sums provide noisy information about the respective unknown permutations. Consequently, variants of such procedures are popular in the literature [CM16, FMR16]. However, such a procedure does not take simultaneous advantage of the fact that the underlying matrix is monotonic in both dimensions. To improve upon simply sorting row (resp. column) sums, we propose an algorithm that first sorts the columns (resp. rows) of the matrix approximately, and then exploits this approximate ordering to sort the rows (resp. columns) of the matrix.

We need more notation to facilitate the description of the algorithm. For a partition

$\mathsf{bl}=(\mathsf{bl}_{1},\dots,\mathsf{bl}_{K})$ of the set $[n_{2}]$111$\mathsf{bl}$ is a partition of $[n_{2}]$ if $[n_{2}]=\cup_{k=1}^{K}\mathsf{bl}_{k}$ and $\mathsf{bl}_{j}\cap\mathsf{bl}_{k}=\emptyset$ for $j\neq k$, we group the columns of a matrix $Y\in{}^{n_{1}\times n_{2}}$ into $K$ blocks according to their indices in $\mathsf{bl}$, and refer to $\mathsf{bl}$ as a partition or blocking of the columns of $Y$.

Given a data matrix $Y\in{}^{n_{1}\times n_{2}}$, the following blocking subroutine returns a column partition $\mathsf{BL}(Y)$. In the main algorithm, partial row sums are computed on indices contained in each block.

Subroutine 1 (blocking)

• •

Step 1: Compute the column sums $\{C(j)\}_{j=1}^{n_{2}}$ of the matrix $Y$ as

[TABLE]

Let $\widehat{\sigma}_{\sf pre}$ be the permutation along which the sequence $\{C(\widehat{\sigma}_{\sf pre}(j))\}_{j=1}^{n_{2}}$ is nondecreasing.

• •

Step 2: Set \tau=16(\zeta+1)\Big{(}\sqrt{\frac{n_{1}^{2}n_{2}}{N}\log(n_{1}n_{2})}+\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})\Big{)} and $K=\lceil n_{2}/\tau\rceil$. Partition the columns of $Y$ into $K$ blocks by defining

[TABLE]

Note that each block is contiguous when the columns are permuted by $\widehat{\sigma}_{\sf pre}$.

• •

Step 3 (aggregation): Set $\beta=n_{2}\sqrt{\frac{n_{1}}{N}\log(n_{1}n_{2})}$. Call a block $\mathsf{bl}_{k}$ “large” if $|\mathsf{bl}_{k}|\geq\beta$ and “small” otherwise. Aggregate small blocks in $\mathsf{bl}$ while leaving the large blocks as they are, to obtain the final partition $\mathsf{BL}$.

More precisely, consider the matrix $Y^{\prime}=Y(\mathsf{id},\widehat{\sigma}_{\sf pre})$ having nondecreasing column sums and contiguous blocks. Call two small blocks “adjacent” if there is no other small block between them. Take unions of adjacent small blocks to ensure that the size of each resulting block is in the range $[\frac{1}{2}\beta,2\beta]$. If the union of all small blocks is smaller than $\frac{1}{2}\beta$, aggregate them all.

Return the resulting partition $\mathsf{BL}(Y)=\mathsf{BL}$.

The threshold $\tau$ is a chosen to be a high probability bound on the perturbation of any column sum. In particular, this ensures that we obtain blocks containing columns that are close when the matrix is ordered according to the correct permutation. Computing partial row sums within each block then provides more refined information about the underlying row permutation than simply computing full row sums, and this intuition underlies the two-dimensional sorting algorithm to follow. As a technical detail, it is important to note that Step 3 aggregates small blocks into large enough ones to reduce noise in these partial row sums. We are now in a position to describe the two-dimensional sorting algorithm.

Algorithm 2 (two-dimensional sorting)

• •

Step 0: Split the observations into two independent subsamples of equal size, and form the corresponding matrices $Y^{(1)}$ and $Y^{(2)}$ according to equation (3).

• •

Step 1: Apply Subroutine 1 to the matrix $Y^{(1)}$ to obtain a partition $\mathsf{BL}=\mathsf{BL}(Y^{(1)})$ of the columns. Let $K$ be the number of blocks in $\mathsf{BL}$.

• •

Step 2: Using the second sample $Y^{(2)}$, compute the row sums

[TABLE]

and the partial row sums within each block

[TABLE]

Create a directed graph $G$ with vertex set $[n_{1}]$, where an edge $u\to v$ is present if either

[TABLE]

• •

Step 3: Compute a topological sort $\widehat{\pi}_{\sf tds}$ of the graph $G$; if none exists, set $\widehat{\pi}_{\sf tds}=\mathsf{id}$.

• •

Step 4: Repeat Steps 1–3 with $(Y^{(i)})^{\top}$ replacing $Y^{(i)}$ for $i=1,2$, the roles of $n_{1}$ and $n_{2}$ switched, and the roles of $\pi$ and $\sigma$ switched, to compute the permutation estimate $\widehat{\sigma}_{\sf tds}$.

• •

Step 5: Return the permutation estimates $(\widehat{\pi}_{\sf tds},\widehat{\sigma}_{\sf tds})$.

Recall that a permutation $\pi$ is called a topological sort of $G$ if $\pi(u)<\pi(v)$ for every directed edge $u\to v$. The construction of the graph $G$ in Step 2 dominates the computational complexity, and takes time $\mathcal{O}(n_{1}^{2}n_{2}/\beta)=\mathcal{O}(n_{1}^{2}n_{2}^{1/2})$. We have the following guarantee for the two-dimensional sorting algorithm.

Theorem 2**.**

For any matrix $M^{*}\in\mathbb{C}_{\mathsf{Perm}}$, we have

[TABLE]

with probability at least $1-9(n_{1}n_{2})^{-3}$.

In particular, setting $N=n_{1}n_{2}$, we have proved that our efficient estimator enjoys the rate

[TABLE]

which is the main theoretical guarantee established in this paper for permutation-based models.

4 Applications

We now discuss in detail how the matrix models studied in this paper arise in practice. The class $\mathbb{C}_{\mathsf{Perm}}$ was studied as a permutation-based model for crowd-labeling [SBW16b] in the case of binary questions, and was proposed as a strict generalization of the classical Dawid-Skene model [DS79, KOS11b, LPI12, DDKR13, GKM11]. Here there is a set of $n_{2}$ questions of a binary nature; the true answer to these questions can be represented by a vector $x^{*}\in\{0,1\}^{n_{2}}$, and our goal is to estimate this vector by asking these questions to $n_{1}$ workers on a crowdsourcing platform. A key to this problem is being able to model the probabilities with which workers answer questions correctly, and we do so by collecting these probabilities within a matrix $M^{*}\in[0,1]^{n_{1}\times n_{2}}$. Assuming that workers have a strict ordering $\pi$ of their abilities, and that questions have a strict ordering $\sigma$ of their difficulties, the matrix $M^{*}$ is bivariate isotonic when the rows are ordered in increasing order of worker ability, and columns are ordered in decreasing order of question difficulty. However, since worker abilities and question difficulties are unknown a priori, the matrix of probabilities obeys the inclusion $M^{*}\in\mathbb{C}_{\mathsf{Perm}}$.

In the calibration problem, we would like to ask questions whose answers we know a priori, so that we can estimate worker abilities and question difficulties, or more generally, the entries of the matrix $M^{*}$. This corresponds to estimating matrices in the class $\mathbb{C}_{\mathsf{Perm}}$ from noisy observations of their entries, whose rate of estimation is our main result.

A subclass of $\mathbb{C}_{\mathsf{Perm}}$ specializes to the case $n_{1}=n_{2}=n$, and also imposes an additional skew symmetry constraint. More precisely, define $\mathbb{C}_{\mathsf{BISO}}^{\prime}$ analogously to the class $\mathbb{C}_{\mathsf{BISO}}$, except with matrices having columns that are nonincreasing instead of nondecreasing. Also define the class $\mathbb{C}_{\mathsf{skew}}(n):\,=\{M\in[0,1]^{n_{1}\times n_{2}}:M+M^{\top}=11^{\top}\}$, and the strong stochastic transitivity class

[TABLE]

The class $\mathbb{C}_{\mathsf{SST}}(n)$ is useful as a model for estimation from pairwise comparisons [Cha15, SBGW17], and was proposed as a strict generalization of parametric models for this problem [BT52, NOS16, RA14]. In particular, given $n$ items obeying some unknown underlying ranking $\pi$, entry $(i,j)$ of a matrix $M^{*}\in\mathbb{C}_{\mathsf{SST}}(n)$ represents the probability $\Pr(i\succ j)$ with which item $i$ beats item $j$ in a pairwise comparison between them. The shape constraint encodes the transitivity condition that for all triples $(i,j,k)$ obeying $\pi(i)<\pi(j)<\pi(k)$, we must have

[TABLE]

For a more classical introduction to these models, see the papers [Fis73, ML65, BW97] and the references therein. Our task is to estimate the underlying ranking from results of passively chosen pairwise comparisons222Such a passive, simultaneous setting should be contrasted with the active case (e.g., [HSRW16, FOPS17, AAAK17]), where we may sequentially choose pairs of items to compare depending on the results of previous comparisons. between the $n$ items, or more generally, to estimate the underlying probabilities $M^{*}$ that govern these comparisons333Accurate, proper estimates of $M^{*}$ translate to accurate estimates of the ranking $\pi$ (see Shah et al. [SBGW17]).. All the results we obtain in this work clearly extend to the class $\mathbb{C}_{\mathsf{SST}}(n)$ with minimal modifications; for example, either of the two estimates $\widehat{\pi}_{\sf tds}$ or $\widehat{\sigma}_{\sf tds}$ may be returned as an estimate of the permutation $\pi$. Consequently, the informal theorem stated in the introduction is an immediate corollary of Theorem 2 once these modifications are made to the algorithm.

5 Proofs

Throughout the proofs, we assume without loss of generality that $M^{*}\in\mathbb{C}_{\mathsf{BISO}}(\mathsf{id},\mathsf{id})=\mathbb{C}_{\mathsf{BISO}}$. Because we are interested in rates of estimation up to universal constants, we assume that each independent subsample contains $N^{\prime}=\operatorname{\mathsf{Poi}}(N)$ observations (instead of $\operatorname{\mathsf{Poi}}(N)/2$ or $\operatorname{\mathsf{Poi}}(N)/4$). We use the shorthand $Y=Y\left(\{(X_{\ell},y_{\ell})\}_{\ell=1}^{N^{\prime}}\right)$, throughout.

5.1 Some preliminary lemmas

Before turning to the proof of Theorems 1 and 2, we provide three lemmas that underlie many of our arguments. The first lemma can be readily distilled from the proof of Theorem 5 of Shah et al. [SBGW17] with slight modifications. It is worth mentioning that similar lemmas characterizing the estimation error of a bivariate isotonic matrix were also proved by [CGS18, CM16].

Lemma 1** ([SBGW17]).**

Let $n_{1}\geq n_{2}$, and let $M^{*}\in\mathbb{C}_{\mathsf{Perm}}$. Assume that our observation model takes the form $Y=M^{*}+W$, where the noise matrix $W$ satisfies the properties

1. (a)

the entries $W_{i,j}$ are independent, centered, $\frac{c_{1}}{p_{\mathsf{obs}}}(\zeta\lor 1)$-sub-Gaussian random variables; 2. (b)

the second moments are bounded as ${\mathbb{E}}[|W_{i,j}|^{2}]\leq\frac{c_{2}}{p_{\mathsf{obs}}}(\zeta^{2}\lor 1)$ for all $i\in[n_{1}],j\in[n_{2}]$.

Then the least squares estimator $\widehat{M}_{{\sf LS}}(Y)$ satisfies

[TABLE]

Moreover, the same result holds if the class $\mathbb{C}_{\mathsf{Perm}}$ is replaced by the class $\mathbb{C}_{\mathsf{BISO}}$.

The proof closely follows that of Shah et al. [SBGW17, Theorem 5]; consequently, we postpone it to Appendix A. The next lemma establishes concentration of sums of our observations around their means.

Lemma 2**.**

For any nonempty subset $\mathcal{S}\subset[n_{1}]\times[n_{2}]$, it holds that

[TABLE]

Proof.

According to definitions (1) and (3), we have

[TABLE]

where $W^{\prime}$ is a $\zeta$-sub-Gaussian noise matrix with independent entries. Consequently, we can express the noise on each entry as $W_{i,j}=Z^{(1)}_{i,j}+Z^{(2)}_{i,j}$ where $\{Z^{(1)}_{i,j}\}_{i\in[n_{1}],j\in[n_{2}]}$ are independent, zero-mean random variables given by

[TABLE]

and $\{Z^{(2)}_{i,j}\}_{i\in[n_{1}],j\in[n_{2}]}$ are independent, zero-mean random variables such that

[TABLE]

We control the two separately. First, we have $|Z^{(1)}_{i,j}|\leq 1/p_{\mathsf{obs}}$ and the variance of each $Z^{(1)}_{i,j}$ is bounded by $(1-p_{\mathsf{obs}})^{2}/p_{\mathsf{obs}}+(1-p_{\mathsf{obs}})\leq 1/p_{\mathsf{obs}}$. Hence Bernstein’s inequality for bounded noise yields

[TABLE]

Taking $t=4\sqrt{\frac{|\mathcal{S}|n_{1}n_{2}}{N}\log(n_{1}n_{2})}+6\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})$ and recalling that $p_{\mathsf{obs}}\geq\frac{N}{2n_{1}n_{2}}$, we obtain

[TABLE]

In order to control the deviation of the sum of $Z^{(2)}_{i,j}$, we note that the $q$-th moment of $Z^{(2)}_{i,j}$ is bounded by $\frac{N}{n_{1}n_{2}}(\frac{2\zeta}{p_{\mathsf{obs}}}\sqrt{q})^{q}\leq\frac{q!}{2}\frac{8\zeta^{2}n_{1}n_{2}}{N}(\frac{4\zeta n_{1}n_{2}}{N})^{q-2}.$ Then another version of Bernstein’s inequality [BLM13] yields

[TABLE]

and setting $t=4\log(n_{1}n_{2})$ gives

[TABLE]

Combining the above two deviation bounds completes the proof. ∎

The last lemma is a deterministic result.

Lemma 3**.**

Let $\{a_{i}\}_{i=1}^{n}$ be a nondecreasing sequence of real numbers. If $\pi$ is a permutation in $\mathfrak{S}_{n}$ such that $\pi(i)<\pi(j)$ whenever $a_{j}-a_{i}>\tau$ where $\tau>0$, then $|a_{\pi(i)}-a_{i}|\leq\tau$ for all $i\in[n]$.

Proof.

Suppose that $a_{j}-a_{\pi(j)}>\tau$ for some index $j\in[n]$. Since $\pi$ is a bijection, there must exist an index $i\leq\pi(j)$ such that $\pi(i)>\pi(j)$. However, we then have $a_{j}-a_{i}\geq a_{j}-a_{\pi(j)}>\tau$, which contradicts the assumption. A similar argument shows that $a_{\pi(j)}-a_{j}>\tau$ also leads to a contradiction. Therefore, we obtain that $|a_{\pi(j)}-a_{j}|\leq\tau$ for every $j\in[n]$. ∎

With these lemmas in hand, we are now ready to prove our main theorems.

5.2 Proof of Theorem 1

We split the proof into two parts by proving the upper and lower bounds separately.

5.2.1 Proof of upper bound

The upper bound follows from Lemma 1 once we check the conditions on the noise for our model. We have seen in the proof of Lemma 2 that the noise on each entry can be written as $W_{i,j}=Z^{(1)}_{i,j}+Z^{(2)}_{i,j}$. Again, $Z^{(1)}_{i,j}$ and $Z^{(2)}_{i,j}$ are $\frac{c}{p_{\mathsf{obs}}}$-sub-Gaussian and $\frac{c\,\zeta}{p_{\mathsf{obs}}}$-sub-Gaussian respectively, and have variances bounded by $\frac{1}{p_{\mathsf{obs}}}$ and $\frac{c\,\zeta^{2}}{p_{\mathsf{obs}}}$ respectively. Hence the conditions on $W$ in Lemma 1 are satisfied. Then we can apply the lemma, recall the relation $p_{\mathsf{obs}}\geq\frac{N}{2n_{1}n_{2}}$ and normalize the bound by $\frac{1}{n_{1}n_{2}}$ to complete the proof.

5.2.2 Proof of lower bound

The lower bound follows from an application of Fano’s lemma. The technique is standard, and we briefly review it here. Suppose we wish to estimate a parameter $\theta$ over an indexed class of distributions $\mathcal{P}=\{\mathbb{P}_{\theta}\,\mid\,\theta\in\Theta\}$ in the square of a (pseudo-)metric $\rho$. We refer to a subset of parameters $\{\theta^{1},\theta^{2},\ldots,\theta^{K}\}$ as a local $(\delta,\epsilon)$-packing set if

[TABLE]

Note that this set is a $\delta$-packing in the metric $\rho$ with the average KL-divergence bounded by $\epsilon$. The following result is a straightforward consequence of Fano’s inequality:

Lemma 4** (Local packing Fano lower bound).**

For any $(\delta,\epsilon)$-packing set of cardinality $K$, we have

[TABLE]

In addition, the Gilbert-Varshamov bound [Gil52, Var57] guarantees the existence of binary vectors $\{v^{1},v^{2},\ldots,v^{K}\}\subseteq\{0,1\}^{n_{1}}$ such that

[TABLE]

for some fixed tuple of constants $(c_{1},c_{2})$. We use this guarantee to design a packing of matrices in the class $\mathbb{C}_{\mathsf{Perm}}$. For each $i\in[K]$, fix some $\delta\in[0,1/4]$ to be precisely set later, and define the matrix $M^{i}$ having identical columns, with entries given by

[TABLE]

Clearly, each of these matrices $\{M^{i}\}_{i=1}^{K}$ is a member of the class $\mathbb{C}_{\mathsf{Perm}}$, and each distinct pair of matrices $(M^{i},M^{j})$ satisfies the inequality $\|M^{i}-M^{j}\|_{F}^{2}\geq c_{2}n_{1}n_{2}\delta^{2}$.

Let $\mathbb{P}_{M}$ denote the probability distribution of the observations in the model (2) with underlying matrix $M\in\mathbb{C}_{\mathsf{Perm}}$. Our observations are independent across entries of the matrix, and so the KL divergence tensorizes to yield

[TABLE]

Let us now examine one term of this sum. We observe $T_{k,\ell}=\operatorname{\mathsf{Poi}}(\frac{N}{n_{1}n_{2}})$ samples of entry $(k,\ell)$; conditioned on the event $T_{k,\ell}=m$, we have the distributions

[TABLE]

Consequently, the KL divergence conditioned on $T_{k,\ell}=m$ is given by

[TABLE]

where we have used $D(p\|q)=p\log(\frac{p}{q})+(1-p)\log(\frac{1-p}{1-q})$ to denote the KL divergence between the Bernoulli random variables $\mbox{\sf Ber}(p)$ and $\mbox{\sf Ber}(q)$.

Note that for $p,q\in[1/2,3/4]$, we have

[TABLE]

Here, step ${\sf(i)}$ follows from the inequality $\log x\leq x-1$, and step ${\sf(ii)}$ from the assumption $q\in[\frac{1}{2},\frac{3}{4}]$. Taking the expectation with respect to $T_{k,\ell}$, we have

[TABLE]

Summing over $k\in[n_{1}],\ell\in[n_{2}]$ yields $D(\mathbb{P}_{M^{i}}\|\mathbb{P}_{M^{j}})\leq\frac{16}{3}N\delta^{2}.$

Substituting into the Fano’s inequality (7), we have

[TABLE]

Finally, choosing $\delta^{2}=c\frac{n_{1}}{N}$ and normalizing by $n_{1}n_{2}$ yields the claim.

5.3 Proof of Proposition 1

Recall the definition of $\widehat{M}(\widehat{\pi},\widehat{\sigma})$ in the meta-algorithm, and additionally, define the projection of any matrix $M\in\mathbb{R}^{n_{1}\times n_{2}}$, as

[TABLE]

and letting $W=Y^{(2)}-M^{*}$, we have

[TABLE]

where step ${\sf(ii)}$ follows from the non-expansiveness of a projection onto a convex set, and steps ${\sf(i)}$ and ${\sf(iii)}$ from the triangle inequality.

The first term in (10) is the estimation error of a bivariate isotonic matrix with known permutations. Since the sample used to obtain $(\widehat{\pi},\widehat{\sigma})$ is independent from the sample used in the projection step, it is equivalent to control the error $\|\mathcal{P}_{\mathsf{id},\mathsf{id}}(M^{*}+W)-M^{*}\|_{F}^{2}$. As before, the noise matrix $W$ satisfies the conditions of Lemma 1. Therefore, applying Lemma 1 in the case $M^{*}\in\mathbb{C}_{\mathsf{BISO}}$ with $p_{\mathsf{obs}}\geq\frac{N}{2n_{1}n_{2}}$ yields the desired bound of order $(\zeta^{2}\lor 1)\frac{n_{1}\log^{2}n_{1}}{N}$.

It remains to bound the second term of (10), the approximation error of the permutation estimates. Note that the approximation error can be split into two components: one along the rows of the matrix, and the other along the columns. More explicitly, we have

[TABLE]

Recall that we assumed without loss of generality that the true permutations are identity permutations, so this completes the proof of Proposition 1. The proof readily extends to the general case by precomposing $\widehat{\pi}$ and $\widehat{\sigma}$ with $\pi^{-1}$ and $\sigma^{-1}$ respectively.

5.4 Proof of Theorem 2

Recall that according to Proposition 1, it suffices to bound the approximation error of our permutation estimate $\|M^{*}-M^{*}(\widehat{\pi}_{\sf tds},\mathsf{id})\|_{F}^{2}$. To ease the notation, we use the shorthand

[TABLE]

and for each block $\mathsf{BL}_{k}$ in Algorithm 2 where $k\in[K]$, we use the shorthand

[TABLE]

throughout the proof. Applying Lemma 2 with $\mathcal{S}=\{i\}\times[n_{2}]$ and then with $\mathcal{S}=\{i\}\times\mathsf{BL}_{k}$ for each $i\in[n_{1}],k\in[K]$, we obtain that

[TABLE]

Note that $K\leq n_{2}/\beta\leq n_{2}^{1/2}$, so a union bound over all $n_{1}(K+1)$ events in inequalities (11a) and (11b) yields that $\Pr\{\mathcal{E}\}\geq 1-2(n_{1}n_{2})^{-3}$, where we define the event

[TABLE]

We now condition on event $\mathcal{E}$. Applying the triangle inequality yields that if

[TABLE]

then we have

[TABLE]

It follows that $u<v$ since $M^{*}$ has nondecreasing columns. Thus, by the choice of thresholds $\eta$ and $\eta_{k}$ in inequalities (6a) and (6b), we have guaranteed that every edge $u\to v$ in the graph $G$ is consistent with the underlying permutation $\mathsf{id}$, so a topological sort exists on event $\mathcal{E}$.

Conversely, if we have

[TABLE]

then the triangle inequality implies that

[TABLE]

Hence the edge $u\to v$ is present in the graph $G$, so the topological sort $\widehat{\pi}_{\sf tds}(u)$ satisfies the relation $\widehat{\pi}_{\sf tds}(u)<\widehat{\pi}_{\sf tds}(v)$. Claim that this allows us to obtain the following bounds on event $\mathcal{E}$:

[TABLE]

We now prove inequality (12b). The proof of inequality (12a) follows in the same fashion. We split the proof into two cases.

Case 1.

First, suppose that $|\mathsf{BL}_{k}|\geq\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})$. Applying Lemma 3 with $a_{i}=\sum_{\ell\in\mathsf{BL}_{k}}M^{*}_{i,\ell}$, $\pi=\widehat{\pi}_{\sf tds}$ and $\tau=2\eta_{k}$, we see that for all $i\in[n_{1}]$,

[TABLE]

Case 2.

Otherwise, we have $|\mathsf{BL}_{k}|\leq\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})$. It then follows that

[TABLE]

where we have used the fact that $M\in[0,1]^{n_{1}\times n_{2}}$.

Next, we consider concentration of the column sums of $Y^{(1)}$. Applying Lemma 2 again with $\mathcal{S}=[n_{1}]\times\{j\}$, we obtain that

[TABLE]

for all $j\in[n_{2}]$ with probability at least $1-2(n_{1}n_{2})^{-3}$. We carry out the remainder of the proof conditioned on the event of probability at least $1-4(n_{1}n_{2})^{-3}$ that inequalities (12a), (12b) and (13) hold.

Having stated the necessary bounds, we now split the remainder of the proof into two parts for convenience. In order to do so, we first split the set $\mathsf{BL}$ into two disjoint sets of blocks, depending on whether a block comes from an originally large block (of size larger than $\beta=n_{2}\sqrt{\frac{n_{1}}{N}\log(n_{1}n_{2})}$ as in Step 3 of Subroutine 1) or from an aggregation of small blocks. More formally, define the sets

[TABLE]

For a set of blocks $\mathsf{B}$, define the shorthand $\cup\mathsf{B}=\bigcup_{B\in\mathsf{B}}B$ for convenience. We begin by focusing on the blocks $\mathsf{BL}^{\mathbb{L}}$.

5.4.1 Error on columns indexed by $\cup\mathsf{BL}^{\mathbb{L}}$

Recall that when the columns of the matrix are ordered according to $\widehat{\sigma}_{\sf pre}$, the blocks in $\mathsf{BL}^{\mathbb{L}}$ are contiguous and thus have an intrinsic ordering. We index the blocks according to this ordering as $B_{1},B_{2},\ldots,B_{\ell}$ where $\ell=|\mathsf{BL}^{\mathbb{L}}|$. Now define the disjoint sets

[TABLE]

Let $\ell_{t}=|\mathsf{BL}^{(t)}|$ for each $t=1,2$.

Recall that each block $B_{k}$ in $\mathsf{BL}^{\mathbb{L}}$ remains unchanged after aggregation, and that the threshold we used to block the columns is \tau=16(\zeta+1)\big{(}\sqrt{\frac{n_{1}^{2}n_{2}}{N}\log(n_{1}n_{2})}+2\frac{n_{1}n_{2}}{N}\log(n_{1}n_{2})\big{)}. Hence, applying the concentration bound (13) together with the definition of blocks in Step 2 of Subroutine 1 yields

[TABLE]

where we again used the argument leading to claim (12b) to combine the two terms. Moreover, since the threshold is twice the concentration bound, it holds that under the true ordering $\mathsf{id}$, every index in $B_{k}$ precedes every index in $B_{k+2}$ for any $k\in[K-2]$. By definition, we have thus ensured that the blocks in $\mathsf{BL}^{(t)}$ do not “mix” with each other.

The rest of the argument hinges on the following lemma, which is proved in Section 5.4.3.

Lemma 5**.**

For $m\in\mathbb{Z}_{+}$, let $J_{1}\sqcup\cdots\sqcup J_{\ell}$ be a partition of $[m]$ such that each $J_{k}$ is contiguous and $J_{k}$ precedes $J_{k+1}$. Let $a_{k}=\min J_{k}$, $b_{k}=\max J_{k}$ and $m_{k}=|J_{k}|$. Let $A$ be a matrix in $[0,1]^{n\times m}$ with nondecreasing rows and nondecreasing columns. Suppose that

[TABLE]

Additionally, suppose that there are positive reals $\rho,\rho_{1},\rho_{2},\ldots,\rho_{\ell}$, and a permutation $\pi$ such that for any $i\in[n]$, we have ${\sf(i)}$ $\sum_{j=1}^{m}|A_{\pi(i),j}-A_{i,j}|\leq\rho$, and ${\sf(ii)}$ $\sum_{j\in J_{k}}|A_{\pi(i),j}-A_{i,j}|\leq\rho_{k}$ for each $k\in[\ell]$. Then it holds that

[TABLE]

We apply the lemma as follows. For $t=1,2$, let the matrix $M^{(t)}$ be the submatrix of $M^{*}$ restricted to the columns indexed by the indices in $\cup\mathsf{BL}^{(t)}$. The matrix $M^{(t)}$ has nondecreasing rows and columns by assumption. We have shown that the blocks in $\mathsf{BL}^{(t)}$ do not mix with each other, so they are contiguous and correctly ordered in $M^{(t)}$. Moreover, the inequality assumptions of the lemma correspond to (15), (12a) and (12b) respectively, with the substitutions

[TABLE]

and setting $J_{1},\dots,J_{\ell}$ to be the blocks in $\mathsf{BL}^{(t)}$. Therefore, applying Lemma 5 yields

[TABLE]

where step ${\sf(i)}$ follows from the Cauchy-Schwarz inequality, and step ${\sf(ii)}$ follows from the fact that $\min_{B\in\mathsf{BL}^{(t)}}|B|\geq\beta=n_{2}\sqrt{\frac{n_{1}}{N}\log(n_{1}n_{2})}$ so that $\ell_{t}\leq n_{2}/\beta$. Substituting for $\beta$ and normalizing by $n_{1}n_{2}$ yields

[TABLE]

This proves the required result for the set of blocks $\mathsf{BL}^{(t)}$. Summing over $t=1,2$ then yields a bound of twice the size for columns of the matrix indexed by $\cup\mathsf{BL}^{\mathbb{L}}$.

5.4.2 Error on columns indexed by $\cup\mathsf{BL}^{\mathbb{S}}$

Next we bound the approximation error of each row of the matrix with column indices restricted to the union of all small blocks. In the easy case where $\mathsf{BL}^{\mathbb{S}}$ contains a single block of size less than $\frac{1}{2}n_{2}\sqrt{\frac{n_{1}}{N}\log(n_{1}n_{2})}$, we have

[TABLE]

where step ${\sf(i)}$ follows from the Hölder’s inequality and the fact that $M^{*}\in[0,1]^{n_{1}\times n_{2}}$, step ${\sf(ii)}$ from the monotonicity of the columns of $M^{*}$, and step ${\sf(iii)}$ from equation (12a).

Now we aim to prove a bound of the same order for the general case. Critical to our analysis is the following lemma:

Lemma 6**.**

For a vector $v\in\mathbb{R}^{n}$, define its variation as $\operatorname{\mathsf{var}}(v)=\max_{i}v_{i}-\min_{i}v_{i}$. Then we have

[TABLE]

See Section 5.4.4 for the proof of this claim.

For each $i\in[n_{1}]$, define $\Delta^{i}$ to be the restriction of the $i$-th row difference $M^{*}_{\widehat{\pi}_{{\sf tds}}(i)}-M^{*}_{i}$ to the union of blocks $\cup\mathsf{BL}^{\mathbb{S}}$. For each block $B\in\mathsf{BL}^{\mathbb{S}}$, denote the restriction of $\Delta^{i}$ to $B$ by $\Delta^{i}_{B}$. Lemma 6 applied with $v=\Delta^{i}_{B}$ yields

[TABLE]

We now analyze the quantities in inequality (17). By the aggregation step of Subroutine 1, we have $\frac{1}{2}\beta\leq|B|\leq 2\beta$, where $\beta=n_{2}\sqrt{\frac{n_{1}}{N}\log(n_{1}n_{2})}$. Additionally, the bounds (12a) and (12b) imply that

[TABLE]

Moreover, to bound the quantity $\sum_{B\in\mathsf{BL}^{\mathbb{S}}}\operatorname{\mathsf{var}}\left(\Delta^{i}_{B}\right)$, we proceed as in the proof for the large blocks in $\mathsf{BL}^{\mathbb{L}}$. Recall that if we permute the columns by $\widehat{\sigma}_{\sf pre}$ according to the column sums, then the blocks in $\mathsf{BL}^{\mathbb{S}}$ have an intrinsic ordering, even after adjacent small blocks are aggregated. Let us index the blocks in $\mathsf{BL}^{\mathbb{S}}$ by $B_{1},B_{2},\ldots,B_{m}$ according to this ordering, where $m=|\mathsf{BL}^{\mathbb{S}}|$. As before, the odd-indexed (or even-indexed) blocks do not mix with each other under the true ordering $\mathsf{id}$, because the threshold used to define the blocks is larger than twice the column sum perturbation. We thus have

[TABLE]

where inequality ${\sf(i)}$ holds because the odd (or even) blocks do not mix, and inequality ${\sf(ii)}$ holds because $M^{*}$ has monotone rows in $[0,1]^{n_{2}}$.

Finally, putting together all the pieces, we can substitute for $\beta$, sum over the indices $i\in n_{1}$, and normalize by $n_{1}n_{2}$ to obtain

[TABLE]

and so the error on columns indexed by the set $\cup\mathsf{BL}^{\mathbb{S}}$ is bounded as desired.

Combining the bounds (16) and (18), we conclude that

[TABLE]

with probability at least $1-4(n_{1}n_{2})^{-3}$. The same proof works with the roles of $n_{1}$ and $n_{2}$ switched and all the matrices transposed, so it holds with the same probability that

[TABLE]

Consequently,

[TABLE]

with probability at least $1-8(n_{1}n_{2})^{-3}$, where we have used the relation $n_{1}\geq n_{2}$. Applying Proposition 1 completes the proof.

5.4.3 Proof of Lemma 5

Since $A$ has increasing rows, for any $i,i_{2}\in[n]$ with $i\leq i_{2}$ and any $j,j_{2}\in J_{k}$, we have

[TABLE]

Choosing $j_{2}=\arg\min_{r\in J_{k}}(A_{i_{2},r}-A_{i,r})$, we obtain

[TABLE]

Together with the assumption on $\pi$, this implies that

[TABLE]

Hence it follows that

[TABLE]

According to the assumptions, we have

$\sum_{k=1}^{\ell}x_{i,k}\leq 1$ and $\sum_{i=1}^{n}x_{i,k}\leq\tau$ for any $i\in[n],k\in[\ell]$; 2. 2.

$\sum_{k=1}^{\ell}y_{i,k}\leq 1$ and $\sum_{i=1}^{n}y_{i,k}\leq\tau$ for any $i\in[n],k\in[\ell]$; 3. 3.

$z_{i,k}\leq\rho_{k}$ and $\sum_{k=1}^{\ell}z_{i,k}\leq\rho$ for any $i\in[n],k\in[\ell]$.

Consequently, the following bounds hold:

$\sum_{i=1}^{n}\sum_{k=1}^{\ell}z_{i,k}x_{i,k}\leq\sum_{i=1}^{n}\sum_{k=1}^{\ell}\rho_{k}x_{i,k}\leq\tau\sum_{k=1}^{\ell}\rho_{k}$; 2. 2.

$\sum_{i=1}^{n}\sum_{k=1}^{\ell}z_{i,k}y_{i,k}\leq\sum_{i=1}^{n}\sum_{k=1}^{\ell}\rho_{k}y_{i,k}\leq\tau\sum_{k=1}^{\ell}\rho_{k}$; 3. 3.

$\sum_{i=1}^{n}\sum_{k=1}^{\ell}z_{i,k}^{2}/m_{k}\leq\sum_{i=1}^{n}\sum_{k=1}^{\ell}z_{i,k}\cdot\max_{k\in[\ell]}(\rho_{k}/m_{k})\leq n\rho\max_{k\in[\ell]}(\rho_{k}/m_{k})$.

Combining these inequalities yields the claim.

5.4.4 Proof of Lemma 6

Let $a=\min_{i\in[n]}v_{i}$ and $b=\max_{i\in[n]}v_{i}=a+\operatorname{\mathsf{var}}(v)$. Since the quantities in the inequality remain the same if we replace $v$ by $-v$, we assume without loss of generality that $b\geq 0$. If $a\leq 0$, then $\|v\|_{\infty}\leq b-a=\operatorname{\mathsf{var}}(v)$. If $a>0$, then $a\leq\|v\|_{1}/n$ and $\|v\|_{\infty}=b\leq\|v\|_{1}/n+\operatorname{\mathsf{var}}(v)$. Hence in any case we have $\|v\|_{2}^{2}\leq\|v\|_{\infty}\|v\|_{1}\leq[\|v\|_{1}/n+\operatorname{\mathsf{var}}(v)]\|v\|_{1}$.

6 Discussion

While the current paper narrows the statistical-computational gap for estimation in permutation-based models with monotonicity constraints, several intriguing questions remain:

• •

Can Algorithm 2 be recursed so as to improve the rate of estimation, until we eventually achieve the statistically optimal rate (up to lower-order terms) in polynomial time?

• •

If not, does there exist a statistical-computational gap in this problem, and if so, what is the fastest rate achievable by computationally efficient estimators?

• •

Can the techniques from here be used to narrow statistical-computational gaps in other permutation-based models [SBW16b, FMR16, PWC17]?

As a partial answer to the first question, it can be shown that when our two-dimensional sorting algorithm is recursed in the natural way and applied to the noisy sorting subclass of the SST model, it yields another minimax optimal estimator for noisy sorting, similar to the multistage algorithm of Mao et al. [MWR17]. However, showing that this same guarantee is preserved for the larger class of SST matrices seems out of the reach of techniques introduced in this paper. In fact, we conjecture that any algorithm that only exploits partial row and column sums cannot achieve a rate faster than $O(n^{-3/4})$ for the SST class.

It is also worth noting that the model (1) allowed us to perform multiple sample-splitting steps while preserving the independence across observations. While our proofs also hold for the observation model where we have exactly $3$ independent samples per entry of the matrix, handling the weak dependence of the sampling model with one observation per entry is an interesting technical challenge that may also involve its own statistical-computational tradeoffs [Mon15].

Acknowledgments

CM thanks Philippe Rigollet for helpful discussions. The work of CM was supported in part by grants NSF CAREER DMS-1541099, NSF DMS-1541100 and ONR N00014-16-S-BA10, and the work of AP and MJW was supported in part by grants NSF-DMS-1612948 and DOD ONR-N00014. We thank Jingyan Wang for pointing out an error in an earlier version of the paper.

Appendix A Proof of Lemma 1

The proof parallels that of Shah et al. [SBGW17, Theorem 5(a)], so we only emphasize the differences and sketch the remaining argument. We may assume that $p_{\mathsf{obs}}\geq\frac{1}{n_{2}}$, since otherwise the bound is trivial.

We first employ a truncation argument. Consider the event

[TABLE]

If the universal constant $c_{3}$ is chosen to be sufficiently large, then it follows from the sub-Gaussianity of $W_{i,j}$ and a union bound over all index pairs $(i,j)\in[n_{1}]\times[n_{2}]$ that $\Pr\{\mathcal{E}\}\geq 1-(n_{1}n_{2})^{-4}$. Now define the truncation operator

[TABLE]

With the choice $\lambda=\frac{c_{3}}{p_{\mathsf{obs}}}(\zeta\lor 1)\sqrt{\log(n_{1}n_{2})}$, define the random variables $W^{(1)}_{i,j}=T_{\lambda}\left(W_{i,j}\right)$ for each pair of indices $(i,j)\in[n_{1}]\times[n_{2}]$. Consider the model where we observe $M^{*}+W^{(1)}$ instead of $Y=M^{*}+W$. Then the new model and the original one are coupled so that they coincide on the event $\mathcal{E}$. Therefore, it suffices to prove a high probability bound assuming that the noise is given by $W^{(1)}$.

Let us define $\mu=\mathbb{E}[W^{(1)}]$ and $\widetilde{W}=W^{(1)}-\mu$. We claim that for any $i\in[n_{1}],j\in[n_{2}]$, the following relations hold:

$|\mu_{i,j}|\leq\frac{c}{p_{\mathsf{obs}}}(\zeta\lor 1)(n_{1}n_{2})^{-4}$; 2. 2.

$\widetilde{W}_{i,j}$ are independent, centered and $\frac{c}{p_{\mathsf{obs}}}(\zeta\lor 1)$-sub-Gaussian; 3. 3.

$|\widetilde{W}_{i,j}|\leq\frac{c}{p_{\mathsf{obs}}}(\zeta\lor 1)\sqrt{\log(n_{1}n_{2})}$; 4. 4.

$\mathbb{E}[|\widetilde{W}_{i,j}|^{2}]\leq\frac{c}{p_{\mathsf{obs}}}(\zeta^{2}\lor 1)$.

Taking these claims as given for the moment, we turn to the main argument assuming that our observations take the form $Y=M^{*}+\widetilde{W}+\mu$.

For any permutations $\pi\in\mathfrak{S}_{n_{1}},\sigma\in\mathfrak{S}_{n_{2}}$, let $M_{\pi,\sigma}=\widehat{M}_{{\sf LS}}(Y)$. We claim that for any fixed pair $(\pi,\sigma)$ such that $\|Y-M_{\pi,\sigma}\|_{F}^{2}\leq\|Y-M^{*}\|_{F}^{2}$, we have

[TABLE]

Treating claim (20) as true for the moment, we see that since the least squares estimator $\widehat{M}$ is equal to $M_{\pi,\sigma}$ for some pair $(\pi,\sigma)$, a union bound over $\pi\in\mathfrak{S}_{n_{1}},\sigma\in\mathfrak{S}_{n_{2}}$ yields

[TABLE]

which completes the proof. Thus, to prove our result, it suffices to prove claim (20).

Let $\Delta_{\pi,\sigma}=M_{\pi,\sigma}-M^{*}$. The condition $\|Y-M_{\pi,\sigma}\|_{F}^{2}\leq\|Y-M^{*}\|_{F}^{2}$ yields the basic inequality

[TABLE]

Since $\Delta_{\pi,\sigma}\in[-1,1]^{n_{1}\times n_{2}}$, we have $\langle\Delta_{\pi,\sigma},\mu\rangle\leq\|\mu\|_{1}\leq\frac{c}{p_{\mathsf{obs}}}(\zeta\lor 1)n_{1}^{-6}$ by claim 1. If it holds that $\|\Delta_{\pi,\sigma}\|_{F}^{2}\leq\frac{4c}{p_{\mathsf{obs}}}(\zeta\lor 1)n_{1}^{-6}$, then the proof is immediate. Thus, we may assume the opposite, from which it follows that

[TABLE]

Consider the set of matrices

[TABLE]

Additionally, for every $t>0$, define the random variable

[TABLE]

For every $t>0$, define the event

[TABLE]

For $t\geq\delta_{n}$, either we already have $\|\Delta_{\pi,\sigma}\|_{F}^{2}\leq t\delta_{n}$, or we have $\|\Delta_{\pi,\sigma}\|_{F}>\sqrt{t\delta_{n}}$. In the latter case, on the complement of $\mathcal{A}_{t}$, we must have $\langle\Delta_{\pi,\sigma},\widetilde{W}\rangle\leq 4\|\Delta_{\pi,\sigma}\|_{F}\sqrt{t\delta_{n}}$. Combining this with inequality (21) then yields $\|\Delta_{\pi,\sigma}\|_{F}^{2}\leq ct\delta_{n}$. It thus remains to bound the probability $\Pr\{\mathcal{A}_{t}\}$.

Using the star-shaped nature of the set $\mathbb{C}_{\mathsf{DIFF}}(\pi,\sigma)$, a rescaling argument yields

[TABLE]

The following lemma bounds the tail behavior of the random variable $Z_{\pi,\sigma}(\delta_{n})$, and its proof is postponed to Section A.2.

Lemma 7**.**

For any $\delta>0$ and $u>0$, we have

[TABLE]

Taking the lemma as given and setting $\delta_{n}^{2}=\frac{c_{2}}{p_{\mathsf{obs}}}(\zeta^{2}\lor 1)n_{1}\log^{2}n_{1}$ and $u=c_{3}(\zeta\lor 1)n_{1}\log^{1.5}n_{1}$, we see that for any $t\geq\delta_{n}$, we have

[TABLE]

In particular, for $t=\delta_{n}$, on the complement of $\mathcal{A}_{t}$, we have

[TABLE]

which completes the proof. Note that the original proof sacrificed a logarithmic factor in proving the equivalent of equation (22), and this is why we recover the same logarithmic factors as in the bounded case in spite of the sub-Gaussian truncation argument.

In the setting where we know that $M^{*}\in\mathbb{C}_{\mathsf{BISO}}$, the same proof clearly works, except that we do not even need to take a union bound over $\pi\in\mathfrak{S}_{n_{1}},\sigma\in\mathfrak{S}_{n_{2}}$ as the columns and rows are ordered.

A.1 Proof of claims 1–4

We assume throughout that the constant $c_{3}$ is chosen to be sufficiently large. Claim 1 follows as a result of the following argument; we have

[TABLE]

By definition, the random variables $W^{(1)}_{i,j}-\mu_{i,j}$ are independent and zero-mean, and applying Lemma 8 (see Appendix C) yields that they are also sub-Gaussian with the claimed variance parameter, thus yielding claim 2. The triangle inequality together with the definition of $\widetilde{W}_{i,j}$ then yields claim 3.

Finally, since $|T(x)|\leq|x|$, we have

[TABLE]

yielding claim 4.

A.2 Proof of Lemma 7

The chaining argument from the proof of Shah et al. [SBGW17, Lemma 10] can applied to show that

[TABLE]

as $\widetilde{W}_{i,j}$ is $\frac{c}{p_{\mathsf{obs}}}(\zeta\lor 1)$-sub-Gaussian by claim 2. Note that although we are considering a set of rectangular matrices $\mathbb{C}_{\mathsf{DIFF}}(\pi,\sigma)\subset[-1,1]^{n_{1}\times n_{2}}$ instead of square matrices as in [SBGW17], we can augment each matrix by zeros to obtain an $n_{1}\times n_{1}$ matrix, and so $\mathbb{C}_{\mathsf{DIFF}}(\pi,\sigma)$ can be viewed as a subset of its counterpart consisting of $n_{1}\times n_{1}$ matrices. Hence the entropy bound depending on $n_{1}$ can be employed so that the chaining argument indeed goes through.

In order to obtain the deviation bound, we apply Lemma 11 of Shah et al. [SBGW17] (i.e., Theorem 1.1(c) of Klein and Rio [KR05]) with $\mathcal{V}=\mathbb{C}_{\mathsf{DIFF}}(\pi,\sigma)\cap\mathcal{B}_{\delta}$, $m=n_{1}n_{2}$, $X=\frac{p_{\mathsf{obs}}}{c(\zeta\lor 1)\sqrt{\log n_{1}}}\widetilde{W}$ and $X^{\dagger}=\frac{p_{\mathsf{obs}}}{c(\zeta\lor 1)\sqrt{\log n_{1}}}Z_{\pi,\sigma}(\delta)$. Claim 3 guarantees that $|X|$ is uniformly bounded by $1$. We also have $\mathbb{E}[\langle D,\widetilde{W}\rangle^{2}]\leq\frac{c}{p_{\mathsf{obs}}}(\zeta^{2}\lor 1)\delta^{2}$ by claim 4 for $\|D\|_{F}^{2}\leq\delta^{2}$. Therefore, we conclude that

[TABLE]

Combining the expectation and the deviation bounds completes the proof.

Appendix B Poissonization reduction

In this section, we show that Poissonization only affects the rates of estimation up to a constant factor. Note that we may assume that $N\geq 4\log(n_{1}n_{2})$, since otherwise, all the bounds in the theorems hold trivially.

Let us first show that an estimator designed for a Poisson number of samples may be employed for estimation with a fixed number of samples. Assume that $N$ is fixed, and we have an estimator $\widehat{M}_{\operatorname{\mathsf{Poi}}}(N)$, which is designed under $N^{\prime}=\operatorname{\mathsf{Poi}}(N)$ observations $\{y_{\ell}\}_{\ell=1}^{N^{\prime}}$. Now, given exactly $N$ observations $\{y_{\ell}\}_{\ell=1}^{N}$ from the model (1), choose an integer $\widetilde{N}=\operatorname{\mathsf{Poi}}(N/2)$, and output the estimator

[TABLE]

Recalling the assumption $N\geq 4\log(n_{1}n_{2})$, we have

[TABLE]

Thus, the error of the estimator $\widehat{M}(N)$, which always uses at most $N$ samples, is bounded by $\frac{1}{n_{1}n_{2}}\|\widehat{M}_{\operatorname{\mathsf{Poi}}}(N/2)-M^{*}\|_{F}^{2}$ with probability greater than $1-(n_{1}n_{2})^{-2}$, and moreover, we have

[TABLE]

We now show the reverse, that an estimator $\widehat{M}(N)$ designed using exactly $N$ samples may be used to estimate $M^{*}$ under a Poissonized observation model. Given $\widetilde{N}=\operatorname{\mathsf{Poi}}(2N)$ samples, define the estimator

[TABLE]

where in the former case, $\widehat{M}(N)$ is computed by discarding $\widetilde{N}-N$ samples at random.

Again, using the fact that $N\geq 4\log(n_{1}n_{2})$ yields

[TABLE]

and so once again, the error of the estimator $\widehat{M}_{\operatorname{\mathsf{Poi}}}(2N)$ is bounded by $\frac{1}{n_{1}n_{2}}\|\widehat{M}(N)-M^{*}\|_{F}^{2}$ with probability greater than $1-(n_{1}n_{2})^{-4}$. A similar guarantee also holds in expectation.

Appendix C Truncation preserves sub-Gaussianity

In this appendix, we show that truncating a sub-Gaussian random variable preserves its sub-Gaussianity to within a constant factor.

Lemma 8**.**

Let $X$ be a (not necessarily centered) $\sigma$-sub-Gaussian random variable, and for some choice $\lambda\geq 0$, let $T_{\lambda}(X)$ denote its truncation according to equation (19). Then $T_{\lambda}(X)$ is $\sqrt{2}\sigma$-sub-Gaussian.

Proof.

The proof follows a symmetrization argument. Let $X^{\prime}$ denote an i.i.d. copy of $X$, and use the shorthand $Y=T_{\lambda}(X)$ and $Y^{\prime}=T_{\lambda}(X^{\prime})$. Let $\varepsilon$ denote a Rademacher random variable that is independent of everything else. Then $Y$ and $Y^{\prime}$ are i.i.d., and $\varepsilon(Y-Y^{\prime})\stackrel{{\scriptstyle d}}{{=}}Y-Y^{\prime}$. Hence we have

[TABLE]

Using the Taylor expansion of $e^{x}$, we have

[TABLE]

since only the even moments remain. Finally, since the map $T_{\lambda}:\mathbb{R}\to\mathbb{R}$ is $1$-Lipschitz, we have $|Y-Y^{\prime}|\leq|X-X^{\prime}|$, and combining this with the fact that $X-X^{\prime}$ has odd moments equal to zero yields

[TABLE]

where the last step follows since the random variable $X-X^{\prime}$ is zero-mean and $\sqrt{2}\sigma$-sub-Gaussian. ∎