Estimation of Monge Matrices

TL;DR

This paper investigates the statistical estimation of Monge and pre-Monge matrices under noise, establishing minimax rates, proposing efficient estimators, and validating results through numerical experiments.

Contribution

It introduces the first minimax analysis for Monge matrix estimation and proposes computationally efficient estimators for pre-Monge matrices.

Findings

Established minimax rates for Monge and pre-Monge matrix estimation.

Proposed two efficient estimators with proven convergence rates.

Validated theoretical results with numerical experiments.

Abstract

Monge matrices and their permuted versions known as pre-Monge matrices naturally appear in many domains across science and engineering. While the rich structural properties of such matrices have long been leveraged for algorithmic purposes, little is known about their impact on statistical estimation. In this work, we propose to view this structure as a shape constraint and study the problem of estimating a Monge matrix subject to additive random noise. More specifically, we establish the minimax rates of estimation of Monge and pre-Monge matrices. In the case of pre-Monge matrices, the minimax-optimal least-squares estimator is not efficiently computable, and we propose two efficient estimators and establish their rates of convergence. Our theoretical findings are supported by numerical experiments.

Equations

$${\theta }_{i,j}+{\theta }_{k,\ell }\le {\theta }_{i,\ell }+{\theta }_{k,j},\text{ for all }1\le i\le k\le n_1,1\le j\le \ell \le n_2.$$

$$\theta ({\pi }_1,{\pi }_2{)}_{i,j}=\theta ({\pi }_1(i),{\pi }_2(j)),\text{ for all }(i,j)\in [n_1]\times [n_2],$$

$${\theta }_{i,j}{\theta }_{k,\ell }\ge {\theta }_{i,\ell }{\theta }_{k,j},\text{ for all }1\le i\le k\le n_1,1\le j\le \ell \le n_2.$$

$$D=\left[\begin{array}{cccccc}-1& 1& 0& \dots & 0& 0\\ 0& -1& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& & ⋮& ⋮\\ 0& 0& 0& \dots & -1& 1\end{array}\right],$$

$$\mathcal{M}={\mathcal{M}}^{n_1,n_2}:=\{\theta \in \text{mathbbm}{R}^{n_1\times n_2}:D\theta {\tilde{D}}^{\top }\ge 0\},$$

$$V(\theta ):={\theta }_{1,1}+{\theta }_{n_1,n_2}-{\theta }_{n_1,1}-{\theta }_{1,n_2}=\Vert D\theta {\tilde{D}}^{\top }{\Vert }_1,$$

$$\text{mathbbm}E[\exp (\mathrm{Tr}(M^{\top }\epsilon )]\le \exp ({\sigma }^2\Vert M{\Vert }_F^2/2).$$

$$y={\theta }^{\ast }+\epsilon ,$$

$${\hat{\theta }}^{\mathsf{ls}}:={\mathrm{argmin}}_{\theta \in \mathcal{M}}\Vert \theta -y{\Vert }_F^2,$$

$$\frac{1}{n_1{n}_2}\Vert \hat{\theta }-{\theta }^{\ast }{\Vert }_F^2.$$

$$\frac{1}{n_1{n}_2}\Vert {\hat{\theta }}^{\mathsf{ls}}-{\theta }^{\ast }{\Vert }_F^2\lesssim \left[\frac{{\sigma }^2}{n_2}+{\left(\frac{{\sigma }^{2}V({\theta }^{\ast })}{n_1{n}_2}\right)}^{2/3}\log (n_1{)}^{1/3}\log (n_2{)}^{2/3}\right]\wedge {\sigma }^2$$

$${\mathcal{M}}_{V_0}={\mathcal{M}}_{V_0}^{n_1,n_2}:=\{\theta \in {\mathcal{M}}^{n_1,n_2}:V(\theta )\le V_0\}.$$

$$\underset{\tilde{\theta}}{\inf}\underset{\begin{subarray}{c}\theta^{*}\in\mathcal{M}_{V_{0}}\end{subarray}}{\sup}\mathbbm{E}\Big{[}\frac{1}{n_{1}n_{2}}\|\tilde{\theta}-\theta^{*}\|_{F}^{2}\Big{]}\gtrsim\Big{[}\frac{\sigma^{2}}{n_{2}}+\left(\frac{\sigma^{2}V_{0}}{n_{1}n_{2}}\right)^{2/3}\Big{]}\land\sigma^{2},$$

$$\mathcal{M}({\pi }_1,{\pi }_2):=\{\theta ({\pi }_1,{\pi }_2):\theta \in \mathcal{M}\}\text{ and }{\mathcal{M}}_{V_0}({\pi }_1,{\pi }_2):=\{\theta ({\pi }_1,{\pi }_2):\theta \in {\mathcal{M}}_{V_0}\}$$

$$y={\theta }^{\ast }({\pi }_1^{\ast },{\pi }_2^{\ast })+\epsilon ,$$

$$\Vert \tilde{\theta }-{\theta }^{\ast }({\pi }_1^{\ast },{\pi }_2^{\ast }){\Vert }_F$$

$$\tilde{\theta }\in \bar{\mathcal{M}}:=\bigcup _{{\pi }_1\in {\mathcal{S}}_{n_1},{\pi }_2\in {\mathcal{S}}_{n_2}}\mathcal{M}({\pi }_1,{\pi }_2),$$

$${\bar{\mathcal{M}}}_{V_0}:=\bar{\mathcal{M}}\cap \{\theta \in \text{mathbbm}{R}^{n_1\times n_2}:V(\theta )\le V_0\}.$$

$${\hat{\theta }}^{\mathsf{gls}}\in {\mathrm{argmin}}_{\theta \in {\bar{\mathcal{M}}}_{V_0}}\Vert \theta -y{\Vert }_F^2.$$

$$\frac{1}{n_1{n}_2}\Vert {\hat{\theta }}^{\mathsf{gls}}-{\theta }^{\ast }({\pi }_1^{\ast },{\pi }_2^{\ast }){\Vert }_F^2\lesssim \left[\frac{{\sigma }^2\log (n_1)}{n_2}+{\left(\frac{{\sigma }^2{V}_0}{n_1{n}_2}\right)}^{2/3}\log (n_1{)}^{1/3}\log (n_2{)}^{2/3}\right]\wedge {\sigma }^2$$

$${\epsilon }_{i,j}\sim \mathsf{subG}(C{\sigma }^2)\text{ and }\mathrm{Var}[{\epsilon }_{i,j}]={\sigma }^2.$$

$$\sum _{k=1}^{n_2}{\left[{\theta }_{i,k}^{\ast }-{\theta }_{1,k}^{\ast }-\frac{1}{n_2}\sum _{\ell =1}^{n_2}({\theta }_{i,\ell }^{\ast }-{\theta }_{1,\ell }^{\ast })\right]}^2,$$

$$\displaystyle\xi{(i,j)}=\sum_{k=1}^{n_{2}}\Big{[}y_{i,k}-y_{j,k}-\frac{1}{n_{2}}\sum_{\ell=1}^{n_{2}}(y_{i,\ell}-y_{j,\ell})\Big{]}^{2},$$

$$(i_0,j_0)={\mathrm{argmax}}_{(i,j)\in [n_1{]}^2,i<j}\xi (i,j).$$

$$\underset{\theta \in {\mathcal{M}}_{V_0}}{\min }\Vert \theta ({\hat{\pi }}_1,{\hat{\pi }}_2)-y{\Vert }_F^2\le \underset{\theta \in {\mathcal{M}}_{V_0}}{\min }\Vert \theta ({\pi }_1^{\mathsf{r}}\circ {\hat{\pi }}_1,{\hat{\pi }}_2)-y{\Vert }_F^2,$$

$$\hat{\theta }:={\mathrm{argmin}}_{\theta \in {\mathcal{M}}_{V_0}({\hat{\pi }}_1^{\prime },{\hat{\pi }}_2)}\Vert \theta -y{\Vert }_F^2.$$

$${\mathcal{M}}_{V_0}({\hat{\pi }}_1,{\hat{\pi }}_2)\cup {\mathcal{M}}_{V_0}({\pi }_1^{\mathsf{r}}\circ {\hat{\pi }}_1,{\hat{\pi }}_2)$$

$$={\mathcal{M}}_{V_0}({\hat{\pi }}_1,{\hat{\pi }}_2)\cup {\mathcal{M}}_{V_0}({\pi }_1^{\mathsf{r}}\circ {\hat{\pi }}_1,{\hat{\pi }}_2)\cup {\mathcal{M}}_{V_0}({\hat{\pi }}_1,{\pi }_2^{\mathsf{r}}\circ {\hat{\pi }}_2)\cup {\mathcal{M}}_{V_0}({\pi }_1^{\mathsf{r}}\circ {\hat{\pi }}_1,{\pi }_2^{\mathsf{r}}\circ {\hat{\pi }}_2).$$

$$\displaystyle\frac{1}{n_{1}n_{2}}\|\hat{\theta}-\theta^{*}(\pi^{\ast}_{1},\pi^{\ast}_{2})\|_{F}^{2}\lesssim\big{(}\sigma^{2}+\sigma V_{0}\big{)}\Big{(}\frac{\log n_{1}}{n_{2}}\Big{)}^{1/2}.$$

$$y=\sum _{i=1}^{n_2}{\lambda }_i{u}_i{v}_i^{\top }.$$

$${\hat{\theta }}^{\mathsf{svt}}:=\sum _{i=1}^{n_2}\text{mathbbm}1\{{\lambda }_i>\rho \}{\lambda }_i{u}_i{v}_i^{\top },$$

Full text

1 Introduction

A matrix $\theta\in\mathbbm{R}^{n_{1}\times n_{2}}$ is called a Monge matrix [38] or a submodular matrix [57], if

[TABLE]

In addition, a matrix $\theta\in\mathbbm{R}^{n_{1}\times n_{2}}$ is called an anti-Monge matrix or a supermodular matrix if $-\theta$ is a Monge matrix. The Monge property dates back to Gaspard Monge’s work on optimal transport [51]. Since then, it has been widely used and studied in optimization, discrete mathematics and computer science [38, 1, 13, 60, 11, 8] as it allows for simple and fast algorithms in a variety of instances [38, 11, 56, 55, 9].

Many of these problems turn out to be invariant under relabeling of the rows and columns of the Monge matrix. Consequently, we introduce the following definition. A matrix $\theta\in\mathbbm{R}^{n_{1}\times n_{2}}$ is called pre-Monge if there exist permutations $\pi_{1}:[n_{1}]\to[n_{1}]$ and $\pi_{2}:[n_{2}]\to[n_{2}]$ such that the matrix $\theta(\pi_{1},\pi_{2})$ defined by

[TABLE]

is Monge. Note that the terminology permuted Monge has also been used to define the same object [11]. A pre-anti-Monge matrix is defined analogously. Like Monge matrices, pre-Monge matrices have also been studied in the context of optimization [10, 14] where the latent permutation yields new computational challenges. For example, even checking that a matrix is pre-Monge is a nontrivial algorithmic task [45, 22].

However, previous work on pre-Monge matrices has focused on the noiseless setting and algorithms typically fail when the pre-Monge matrix is contaminated by random noise. This motivates us to take a statistical approach and study estimation of a pre-Monge matrix under random noise.

1.1 Geometric interpretation and spectral ordering.

The Monge property has strong ties with geometric properties of certain datasets, starting with the seminal work of Monge on optimal transport [51]. In this subsection, we demonstrate how the Monge property arises in the context of seriation [43, 44, 3, 31, 32, 30] where the goal is to recover the latent ordering of objects based on pairwise distances or inner products.

Let $X\in\mathbbm{R}^{n\times d}$ be a data matrix with rows $x_{1}^{\top},\dots,x_{n}^{\top}\in\mathbbm{R}^{d}$. Moreover, suppose that for all $i,j\in[n-1]$, we have that $(x_{i+1}-x_{i})^{\top}(x_{j+1}-x_{j})\geq 0$. In other words, the differences between consecutive points form an acute angle so that the points $x_{1},\ldots,x_{n}$ may be ordered along a common direction. In this case, it is easy to check that the Gram matrix $\theta=XX^{\top}$ is an anti-Monge matrix, and the distance matrix $D$, defined by $D_{i,j}=\|x_{i}-x_{j}\|_{2}^{2}$ for $i,j\in[n]$, is a Monge matrix.

Suppose that we do not know the order of the points, or equivalently, we observe $x_{\pi(1)},\dots,x_{\pi(n)}$, where there is an unknown permutation $\pi:[n]\to[n]$. How can we reorder the points in order to recover the above geometric structure (i.e. so that the differences between consecutive points form an acute angle)? Intuitively, assuming that such a reordering exists suggests that the $n$ points should approximately lie along a hidden direction. Therefore, we can apply principal component analysis as follows. Let us assume without loss of generality that the points are centered so that $\sum_{i=1}^{n}x_{i}=0$ and thus $\sum_{i=1}^{n}\theta_{i,j}=\sum_{j=1}^{n}\theta_{i,j}=0$. Then the leading right singular vector of $X$ gives the hidden direction, and the leading left singular vector $v$ of $X$, i.e., the leading eigenvector of the Gram matrix $\theta$ is the first principal component of the data points. The entries of $v$ then give a one-dimensional embedding of the points, from which we easily recover the original order.

Indeed, this intuition can be made rigorous using Corollary 2.6 of [28], which is a variant of the Perron-Frobenius theorem and states that the leading eigenvector of a doubly centered anti-Monge matrix (i.e. having row and column sums equal to zero) is monotone. Hence the leading eigenvector $v$ of the Gram matrix $\theta$ is monotone. If the unknown permutation $\pi$ relabels the points, then the leading eigenvector of the Gram matrix becomes $v_{\pi}$, defined by $(v_{\pi})_{i}=v_{\pi(i)}$. As a result, sorting the entries of $v_{\pi}$ recovers the permutation $\pi$ and, therefore, the latent order of the points. The above method for spectral ordering is similar to the one for seriation proposed in [3].

1.2 Our contribution

In this work, we study the estimation of pre-(anti-)Monge matrices under additive sub-Gaussian noise. Statistically, we establish the minimax rates of estimation (up to logarithmic factors) for both Monge and pre-Monge matrices in Sections 2 and 3.1 respectively, where the upper bounds are achieved by the least-squares estimators.

Algorithmically, for estimating pre-Monge matrices, we further introduce two efficient estimators and study their rates of convergence. The Variance Sorting estimator introduced in Section 3.2, as the name suggests, employs second-order information to estimate the latent permutation. In Section 3.3, we study the singular value thresholding estimator based on our result (Proposition 6) on the approximation of pre-Monge matrices by low-rank ones.

Furthermore, we provide various numerical experiments in Section 4 to corroborate the theoretically established rates of estimation. Using Dykstra’s projection algorithm, we give a detailed implementation of the least-squares estimator for (anti-)Monge matrices, which is of practical interest.

The proofs of all theorems and auxiliary lemmas can be found in Section 5.

1.3 Related work

This work connects to several lines of research that are described below.

Total positivity.

The Monge property is closely related to the notion of total positivity [41]. An entrywise positive matrix $\theta\in\mathbbm{R}^{n_{1}\times n_{2}}$ is called totally positive (of order $2$), if

[TABLE]

Therefore, an entrywise positive matrix $\theta$ is totally positive if and only if $\log(\theta)$ is anti-Monge, where $\log(\cdot)$ is applied to each entry of $\theta$ individually. As a result, total positivity is also known as log-supermodularity. Total positivity plays an essential role in statistical physics via the FKG inequality [33] and appears frequently in many other areas of probability and statistics [41, 42]. More recently, there have been new developments in studying totally positive distributions and related estimation problems [26, 47, 58]. In a companion paper [39], we study minimax estimation of a totally positive distribution by employing mathematical tools that are closely related to those in the current paper.

Latent permutation learning.

Estimating a pre-Monge matrix from its noisy version falls into the category of matrix learning with latent permutations, which has recently observed a surge of interest. Models involving latent permutations include noisy sorting [7], the strong stochastic transitivity model [17, 62], feature matching [20], crowd labeling [64], statistical seriation [30] and graph matching [48, 24], to name a few. Many of the previous approaches for learning latent permutations under such models are based on sorting row or column sums of the observed matrix (or equivalently, degrees of vertices) [63, 19, 54] or certain refinements [50, 49]. However, since adding a constant to all entries in a row or column of a Monge matrix does not change its Monge property, first-order information such as row sums is uninformative for the Monge structure, and thus cannot be used to identify the latent permutation. Instead, we propose a new algorithm based on variance sorting. We show in Section 3.2 that this novel use of second order information is decisive when estimating pre-Monge matrices.

Graphon estimation.

Another related, substantial body of literature is that on graphon estimation [70, 15, 34, 5], where the goal is to estimate a bivariate function $f:[0,1]^{2}\to\mathbbm{R}$ from samples $\{f(X_{i},Y_{j}):1\leq i\leq n_{1},\,1\leq j\leq n_{2}\}$. Unlike regression, the design points $(X_{i},Y_{j})$ are not observed in graphon estimation, so the observations can be viewed as an $n_{1}\times n_{2}$ matrix with latent permutations acting on its rows and columns. There have been extensive studies on graphon estimation with various structures, including block models [2], smoothness [46] and low-rank structure [61]. While our setup is not about recovering an underlying function $f$, the current work can be viewed as a study of denoising observations in graphon estimation with the Monge structure.

Shape-constrained estimation.

Estimation of a Monge matrix, which we study in Section 2, falls in the scope of shape-constrained estimation. Closest to the present work is the estimation of a bivariate isotonic matrix under Gaussian noise [18]. In fact, every anti-Monge matrix can be written as the sum of a rank-two matrix and a bivariate isotonic matrix (Lemma 7). However, our results suggest that the set of Monge matrices is in fact qualitatively different from the set of bivariate isotonic matrices. Particularly, the minimax rate of estimation in Theorem 1 is different from that given by Theorem 2.1 of [18], and the low-rank approximation rate in Proposition 6 is different from that given by Lemma 4 of [62].

Shortly before completing the current work, we became aware of a concurrent work by Fang, Guntuboyina and Sen [27] that studies multivariate extensions of isotonic regression. The two-dimensional version almost coincides with the anti-Monge structure (without permutations) that we study, and the rate achieved by the least-squares estimator specialized to dimension two, as expected, coincides with the main term of the rate given by Theorem 1 in our current paper. However, it is worth noting that the two proofs follow drastically different paths. While the proof in [27] relies on metric entropy estimates from [4, 35], our proof is based on spectral decomposition of the difference operator $D$ defined in (2.2), a technique which has been used for example to study the performance of total variation regularization [69, 40]. Moreover, assuming $n=n_{1}=n_{2}$, our upper bound given in Theorem 1 contains a log factor of order $\log(n)$, while the one in Theorem 4.1 of [27] potentially scales like $\log(n)^{3}$, a minor improvement which nonetheless shows the potential merits of our proof technique.

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}$ of real numbers, 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$ and $C$ (possibly with subscripts) to denote universal constants that may change from line to line. Let $\land$ and $\lor$ denote the min and the max operators between two real numbers respectively. Given a matrix $M\in\mathbbm{R}^{n_{1}\times n_{2}}$, we denote its $i$-th row by $M_{i,\cdot}$ and its $j$-th column by $M_{\cdot,j}$. We denote by $\|M\|_{F}$ and $\|M\|$ the Frobenius norm and the operator norm of $M$, and by $\|M\|_{1}$ and $\|M\|_{\infty}$ the $\ell^{1}$ and $\ell^{\infty}$-norm of $M$ when viewed as a vector in $\mathbbm{R}^{n_{1}n_{2}}$, respectively. We write $M^{{\dagger}}$ for the Moore-Penrose pseudoinverse of $M$. Finally, let $\mathcal{S}_{n}$ denote the set of permutations $\pi:[n]\to[n]$.

2 Anti-Monge matrix estimation

We start with estimation of a Monge matrix under sub-Gaussian noise, without latent permutations. It is mathematically equivalent to study estimation of an anti-Monge matrix $\theta^{*}\in\mathbbm{R}^{n_{1}\times n_{2}}$, which we find more convenient for the presentation. Throughout this work, without loss of generality, we also assume that $n_{1}\geq n_{2}$.

Consider the difference operator $D\in\mathbbm{R}^{(n_{1}-1)\times n_{1}}$ defined by

[TABLE]

and we define $\tilde{D}\in\mathbbm{R}^{(n_{2}-1)\times n_{2}}$ in the same way.

Using a telescoping sum argument, it is easy to check that the set of anti-Monge matrices $\theta$ such that $-\theta$ satisfies (1.1) can be expressed as

[TABLE]

where the symbol $\geq$ denotes entrywise inequality.

For each $\theta\in\mathcal{M}$, we define the quantity

[TABLE]

where the last equality follows from a telescoping sum. We remark that $V(\theta)$ is a global seminorm of $\theta$, and turns out to play a role in the rate of estimation.

In this work, we consider additive sub-Gaussian noise. Namely, for a zero-mean random matrix $\varepsilon\in\mathbbm{R}^{n_{1}\times n_{2}}$, we say that $\varepsilon$ is sub-Gaussian with variance proxy $\sigma^{2}$, or simply $\varepsilon\sim\mathsf{subG}_{n_{1}\times n_{2}}(\sigma^{2})$, if for any matrix $M\in\mathbbm{R}^{n_{1}\times n_{2}}$, it holds that

[TABLE]

Suppose that we observe

[TABLE]

where $\varepsilon\sim\mathsf{subG}_{n_{1}\times n_{2}}(\sigma^{2})$. We study the performance of the least-squares estimator

[TABLE]

in terms of the mean squared error

[TABLE]

Our upper bound is stated in the following theorem.

Theorem 1**.**

Let $\theta^{*}\in\mathcal{M}^{n_{1},n_{2}}$ be an anti-Monge matrix, and suppose that we observe $y=\theta^{*}+\varepsilon$ where $\varepsilon\sim\mathsf{subG}_{n_{1}\times n_{2}}(\sigma^{2}).$ Then, the least-squares estimator $\hat{\theta}^{\sf ls}$ achieves the rate

[TABLE]

with probability at least $1-\exp(-n_{1})$. Moreover, the same bound holds in expectation.

Assuming Gaussian noise, the following theorem provides a lower bound that matches the above upper bound up to a logarithmic factor. For $V_{0}\geq 0$, let us define

[TABLE]

Theorem 2**.**

*Consider the model $y=\theta^{*}+\varepsilon,$ where $\theta^{*}\in\mathcal{M}_{V_{0}}^{n_{1},n_{2}}$ and $\varepsilon$ has i.i.d. $\mathcal{N}(0,\sigma^{2})$ entries. For any $V_{0}\geq 0$, it holds that *

[TABLE]

where the infimum is taken over all estimators measurable with respect to the observation $y$.

3 Pre-anti-Monge matrix estimation

In this section, we move on to study the estimation of a pre-anti-Monge matrix, that is, an anti-Monge matrix whose rows and columns have been shuffled by latent permutations. Let $\mathcal{S}_{n}$ denote the set of permutations $\pi:[n]\to[n]$. For any matrix $\theta\in\mathbbm{R}^{n_{1}\times n_{2}}$ and permutations $\pi_{1}\in\mathcal{S}_{n_{1}},\pi_{2}\in\mathcal{S}_{n_{2}}$, recall that $\theta(\pi_{1},\pi_{2})$ denotes the matrix defined by $\theta(\pi_{1},\pi_{2})_{i,j}=\theta(\pi_{1}(i),\pi_{2}(j)).$ Define the sets

[TABLE]

of anti-Monge matrices shuffled by fixed permutations.

Suppose that we observe

[TABLE]

where $(\pi_{1}^{*},\pi_{2}^{*},\theta^{*})\in\mathcal{S}_{n_{1}}\times\mathcal{S}_{n_{2}}\times\mathcal{M}$ and $\varepsilon\sim\mathsf{subG}_{n_{1}\times n_{2}}(\sigma^{2})$. Our goal is to estimate the pre-anti-Monge matrix $\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})$.

If two rows (or columns) of $\theta^{*}$ differ by a constant vector, then the matrix we obtain from switching these two rows is still anti-Monge. Therefore, even if the noise $\varepsilon$ is zero, neither the pair of permutations $(\pi_{1}^{*},\pi_{2}^{*})$ nor the matrix $\theta^{*}$ can be inferred from $y$. As a result, measures of permutation and estimation errors such as $\|\theta^{*}(\hat{\pi}_{1},\hat{\pi}_{2})-\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})\|_{F}$ and $\|\hat{\theta}-\theta^{*}\|_{F},$ may be not be pertinent. This is why, instead of studying identifiability of the permutations and the anti-Monge matrix, we focus on the denoising error

[TABLE]

for any estimator $\tilde{\theta}$ of the pre-anti-Monge matrix.

Depending on the application, it might be important to differentiate between proper and improper estimators $\tilde{\theta}$. In this context, a proper estimator is an estimator

[TABLE]

that is, an estimator that needs to be a pre-anti-Monge matrix itself. By contrast, an improper estimator can be any matrix $\tilde{\theta}\in\mathbbm{R}^{n_{1}\times n_{2}}$.

The rest of this section is organized as follows. We first establish the minimax rate for estimating a pre-anti-Monge matrix in Section 3.1. It is achieved by the global least-squares estimator, which is proper by nature, but is likely to be computationally infeasible. Next, we give a computationally feasible proper estimator in Section 3.2 under additional assumptions. Finally, in Section 3.3, we present another computationally feasible estimator based on singular value thresholding that yields a better rate than the one in Section 3.2, but may be improper. This presents a shortcoming if one wants to leverage the Monge structure for downstream numerical computations.

3.1 Minimax rates of estimation

We work under the technical assumption that $\theta^{*}\in\mathcal{M}_{V_{0}}$ where $V_{0}$ is known. Define

[TABLE]

Our upper bound is achieved (up to a logarithmic factor) by the global least-squares estimator over the entire parameter space

[TABLE]

If the minimizer is not unique, an arbitrary one is chosen.

Theorem 3**.**

Suppose that we have $y=\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})+\varepsilon$, where $\theta^{*}\in\mathcal{M}_{V_{0}}^{n_{1},n_{2}}$ and $\varepsilon\sim\mathsf{subG}(\sigma^{2})$. Then the global least-squares estimator (3.13) achieves the rate

[TABLE]

with probability at least $1-n_{1}^{-n_{1}}$. Moreover, the same bound holds in expectation.

Note that this rate is the same (up to a logarithmic factor in the first term) as that for estimating an anti-Monge matrix without latent permutations in view of Theorem 1. Therefore, the lower bound of Theorem 2 for the smaller class implies minimax optimality of the above upper bound (up to a logarithmic factor).

We conjecture that a similar bound holds true for a version of the least-squares estimator where the projection onto $\bar{\mathcal{M}}_{V_{0}}$ is replaced by the unrestricted version $\bar{\mathcal{M}}$, but our current proof technique does not allow us to conclude this.

3.2 Efficient estimation via variance sorting

While the global least-squares estimator retains the minimax rate even in the presence of latent permutations, solving the optimization problem (3.13) is unlikely to be computationally efficient. Thus we now discuss polynomial-time estimators. In this subsection, we assume that the noise matrix $\varepsilon$ is homoscedastic with independent sub-Gaussian entries, i.e.,

[TABLE]

As in the previous section, the estimator is based on projecting a permuted version of the observations onto $\mathcal{M}_{V_{0}}$, but we use an efficient method to find estimators of the permutations with respect to which we project on. Let us first focus on estimating the row permutation $\pi_{1}$. Since adding a constant to all entries in a row of the underlying matrix does not change its anti-Monge property, there is no first-order information that helps distinguish between the rows of $y$. Instead, we exploit second-order information, namely, the variance of row differences of $y$.

The intuition behind the following algorithm is that if we knew the index $\pi_{1}^{-1}(1)$ corresponding to the first row of $\theta^{*}$, the anti-Monge property would imply that the variances between any other row $i\in[n_{1}]$ and row $1$ in the unpermuted matrix $\theta^{\ast}$,

[TABLE]

are monotonically increasing in $i$. Hence, given $\pi_{1}(1)$, we could estimate these variances and sort the rows accordingly. The precise method is given in the following Variance Sorting Subroutine.

Note that in Algorithm 1, the pair $(i_{0},j_{0})$ is an estimator for the extremal rows $\pi_{1}^{-1}(1)$ and $\pi_{1}^{-1}(n_{1})$, but the choice of which index corresponds to $\pi_{1}^{-1}(1)$ is broken arbitrarily by the constraint $i_{0}<j_{0}$. In turn, the resulting estimator $\hat{\pi}_{1}$ can only be reliable up to a global flip of the coordinates. In order to obtain denoising rates, this indeterminacy can be overcome by projecting $y$ onto the set of anti-Monge matrices under both possible orientations and picking the best fit.

To facilitate our presentation, we define the reversal permutation $\pi_{1}^{\mathsf{r}}\in\mathcal{S}_{n_{1}}$ by $\pi_{1}^{\mathsf{r}}(i)=n_{1}-i+1$ for $i\in[n_{1}]$, and define similarly $\pi_{2}^{\mathsf{r}}\in\mathcal{S}_{n_{2}}$ by $\pi_{2}^{\mathsf{r}}(i)=n_{2}-i+1$ for $i\in[n_{2}]$. In short, Algorithm 2 below applies the Variance Sorting subroutine twice to estimate both row and column permutations, and then estimates $\theta$ by the (computationally efficient) least-squares estimator in the convex set of anti-Monge matrices along these estimated permutations.

Note that we only allowed a potential flip $\pi_{1}^{\mathsf{r}}$ for $\hat{\pi}_{1}$, although there is also such an ambiguity for $\hat{\pi}_{2}$. This suffices because if $\theta\in\mathcal{M}_{V_{0}}$, then $\theta(\pi_{1}^{\mathsf{r}},\pi_{2}^{\mathsf{r}})\in\mathcal{M}_{V_{0}}$, and as a result

[TABLE]

The estimator computed by the Main Algorithm achieves the following rate of estimation.

Theorem 4**.**

Suppose that $y=\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})+\varepsilon$, where $\theta^{*}\in\mathcal{M}_{V_{0}}^{n_{1},n_{2}}$ and $\varepsilon$ has independent $\mathsf{subG}(C\sigma^{2})$ entries with variance $\sigma^{2}.$ Let the estimator $\hat{\theta}$ be given by the Main Algorithm. Then it holds with probability at least $1-n_{1}^{-n_{1}}$ that

[TABLE]

Moreover, the same bound holds in expectation.

This rate achieved by our efficient estimator is consistent, but it is suboptimal in view of the minimax rate given by Theorem 3.

3.3 Denoising via singular value thresholding

While the Variance Sorting algorithm above yields efficient estimators of the latent permutations, the rate of convergence it achieves is suboptimal. We now aim for the easier task of denoising the pre-anti-Monge matrix without learning the latent permutations, in the hope of obtaining an efficient estimator with a faster rate of convergence. More precisely, under model (3.10), we look for a possibly improper estimator $\tilde{\theta}\in\mathbbm{R}^{n_{1}\times n_{2}}$ so that \big{\|}\tilde{\theta}-\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})\big{\|}_{F}^{2} is small.

To this end, we consider the well-studied singular value thresholding (SVT) estimator [17, 36]. Let the singular value decomposition of $y$ be

[TABLE]

Then the SVT (hard-thresholding) estimator is defined as

[TABLE]

where we choose the threshold to be $\rho:=C\sigma\sqrt{n_{1}}$ for a sufficiently large constant $C>0$. The rate of estimation achieved by the SVT estimator is given in the following theorem.

Theorem 5**.**

Suppose that we have $y=\theta^{*}(\pi_{1}^{*},\pi_{2}^{*})+\varepsilon$, where $\theta^{*}\in\mathcal{M}^{n_{1},n_{2}}$ and $\varepsilon\sim\mathsf{subG}(\sigma^{2})$. The singular value thresholding estimator $\hat{\theta}^{\sf svt}$ achieves the rate

[TABLE]

with probability at least $1-\exp(-n_{1})$. Moreover, the same bound holds in expectation.

This rate sits between the minimax rate given by Theorem 3, and the rate for the Variance Sorting estimator given by Theorem 4. Note that for this result, the noise $\varepsilon$ needs not be homoscedastic, and moreover, no knowledge of $V_{0}$ is required, i.e., the SVT estimator adapts to the quantity $V(\theta^{*})$.

The proof technique leading to upper bounds for the SVT estimator is well developed [17, 62]. Our contribution mainly lies in the following low-rank approximation result for an anti-Monge matrix, which is of independent interest.

Proposition 6**.**

For any $\theta\in\mathcal{M}^{n_{1},n_{2}}$ and positive integer $r$, there exists a rank-$(3r+3)$ matrix $\tilde{\theta}\in\mathbbm{R}^{n_{1}\times n_{2}}$ such that

[TABLE]

Note that using a similar proof, the same rate as in Theorem 5 can be obtained for a soft-thresholding estimator as well, that is, for

[TABLE]

with a similar scaling for $\rho$.

As the rate given in Theorem 5 does not match the minimax rate, it is natural to ask whether this suboptimality is an artifact of the proof or a true weakness of the SVT estimator. In Appendix B, we present a worst-case anti-Monge matrix which cannot be approximated by any low-rank matrix at a rate better than that given by Proposition 6. This in turn gives evidence that the rate of convergence for the SVT estimator in Theorem 5 might be the best achievable by this method.

4 Numerical experiments

In order to compare our theoretical guarantees with the empirical performance of the proposed estimators, we conducted experiments on synthetic data, using Dykstra’s algorithm to project onto the cone of anti-Monge matrices.

We first present this projection algorithm in Section 4.1. We then show the experimental results of the projection onto the cone of anti-Monge matrices in Section 4.2 and of the two efficient strategies for denoising pre-anti-Monge matrices in Section 4.3.

4.1 Dykstra’s algorithm for projecting onto the set of anti-Monge matrices

Since the set $\mathcal{M}$ is a convex cone specified by $O(n_{1}n_{2})$ constraints, the least-squares estimator (2.6) can be calculated by a general purpose convex optimization software such as SCS [52, 53] or EOCS [25]. The most computationally intensive subroutine of these methods is usually solving linear systems associated with the constraints specifying $\mathcal{M}$. Using direct methods to find these solutions results in a runtime that scales like $(n_{1}n_{2})^{3}$, rendering calculations relatively slow even for moderate values of $n_{1}$ and $n_{2}$. Hence, we chose to implement a specialized algorithm to calculate $\theta$ based on Dykstra’s projection algorithm [6, 21].

In its general form (see Algorithm 3), this algorithm is designed to calculate the projection of a vector $y\in\mathbbm{R}^{d}$ onto the intersection of $m$ convex sets $\mathcal{M}_{1},\dots,\mathcal{M}_{m}$ by iteratively projecting carefully chosen points to each individual set. This is similar to alternate projections of a point to each of the sets $\mathcal{M}_{1},\dots,\mathcal{M}_{m}$, but when initialized with $y\in\mathbbm{R}^{d}$, Dykstra’s algorithm not only finds a point in the intersection $\bigcap_{j\in[m]}\mathcal{M}_{j}$, but its iterates actually converge to the projection of $y$ onto $\bigcap_{j\in[m]}\mathcal{M}_{j}$.

To apply Dysktra’s algorithm to the problem of projecting onto the cone of anti-Monge matrices, note that we can write $\mathcal{M}=\bigcap_{i_{1}=1}^{n_{1}-1}\bigcap_{i_{2}=1}^{n_{2}-1}\mathcal{M}_{i_{1},i_{2}}$ with

[TABLE]

because a matrix is anti-Monge if and only if each contiguous $2\times 2$ submatrix is anti-Monge. The projection of $y$ onto $\mathcal{M}_{i_{1},i_{2}}$ can be explicitly calculated to be the matrix with entries

[TABLE]

for $j_{1},j_{2}\in\{0,1\}$, and

[TABLE]

if $(\ell_{1},\ell_{2})\notin(i_{1}+\{0,1\})\times(i_{2}+\{0,1\})$.

This leads to Algorithm 4 for projecting a matrix $y\in\mathbbm{R}^{n_{1}\times n_{2}}$ onto $\mathcal{M}$.

The rate of convergence of Dykstra’s method can be shown to be linearly exponential in the iterations [23], that is, if we denote by $\theta^{(k)}$ the $k$th iterate of $\theta$ in Algorithm 4 and by $\theta^{\ast}=\Pi_{\mathcal{M}}(y)$, then $\|\theta^{(k)}-\theta^{\ast}\|_{2}\lesssim c^{k}$ for a constant $c<1$. However, note that the constant $c$ may get closer to one with increasing $n_{1}$ and $n_{2}$, which is the case for isotonic regression as shown in [23] and matches our experience: simulations for larger values of $n_{1}$ and $n_{2}$ require more iterates before convergence.

In practice, convergence in Algorithm 4 can be checked by evaluating a measure of feasibility such as $\|D\theta\tilde{D}^{\top}\|_{\infty}$, or by checking when the distance between two successive iterates is small.

4.2 Experiments for anti-Monge matrices

In the following two sections, we assume $n=n_{1}=n_{2}$ for simplicity.

For the estimation of anti-Monge matrices, we consider the following family of ground truth signals, motivated by the construction of the lower bounds in the proof of Theorem 2. First, for $n\in\mathbb{N}$ and $V,\sigma>0$, define $\theta_{1,(n)}\in\mathbbm{R}^{n\times n}$ as

[TABLE]

where $k=(Vn/\sigma)^{1/3}$. The ground truth $\theta^{\ast}_{1,(n)}$ is obtained by centering $\theta_{1,(n)}$ to have zero column and row sums. Finally, we set $y=\theta^{\ast}_{1,(n)}+\varepsilon$ where $\varepsilon_{i,j}\overset{\text{\small}\rm{iid}}{\sim}\mathcal{N}(0,\sigma)$ and report the average denoising error $\|\theta^{\mathsf{ls}}-\theta^{\ast}_{1,(n)}\|_{F}^{2}/n^{2}$ over 20 repetitions.

Our simulations recover the three regimes for $n$ that appear in Theorem 1, although at different signal-to-noise ratios governed by $V/\sigma$. Namely, on the one hand, for $V=\sigma=1$, we see in Figure 1(a) an error decay of $n^{-1.02}\approx n^{-1}$ for $n$ between $10$ and $160$, obtained by linearly regressing the logarithm of the errors onto the logarithm of the $n$ values. On the other hand, for $V=2\cdot 10^{6}$, we can see both a plateau when the trivial $\sigma^{2}$ error bound in Theorem 1 is active, as well as a decay of $n^{-1.34}\approx n^{-4/3}$ at the beginning of the decay becoming effective, where the slope in the doubly logarithmic plot is read off between two consecutive points as indicated in Figure 1(b).

Similarly, fixing $n=200$, $\sigma=1$, and varying $V$ between $10^{-2}$ and $10^{7}$, we can observe a $V^{0.65}\approx V^{2/3}$ scaling in Figure 2(a). The overall curve is shallower, plateauing both at the far low and high end of $V$, corresponding to the $\sigma^{2}/n$ and $\sigma^{2}$ rates becoming active, respectively.

Finally, in Figure 2(b), when setting $n=300$, $V=1$, and varying $\sigma$ between $10^{-7}$ and $1$, we obtain slopes of $\sigma^{2.01}$ and $\sigma^{1.84}$ on the low and high end, while the lowest slope between consecutive points in the curve is $\sigma^{1.34}$, which matches the theoretical rates of $\sigma^{2}$, $\sigma^{2}/n$ and $(V\sigma^{2}/n^{2})^{2/3}$, respectively.

4.3 Experiments for pre-anti-Monge matrices

To illustrate the practical performance of the efficient methods presented for denoising a pre-anti-Monge matrix, Variance Sorting and singular value thresholding (see Sections 3.2 and 3.3 respectively), we further perform experiments by using both methods on the following family of ground truth matrices:

[TABLE]

These were chosen because the singular value decay we proved in Proposition 6 is tight for these matrices (see Lemma 23). By contrast, each ground truth example in the previous subsection, $\theta^{\ast}_{1,(n)}$, is a rank-one matrix, and hence should lead to an overall better performance of singular value thresholding that is independent of $n$.

For the Variance Sorting algorithm, we set $V=1$, $\sigma=0.5$ and report the approximation error induced by the estimated permutations, i.e.,

[TABLE]

for $\theta^{\ast}=\theta^{\ast}_{2,(n)}$, averaged over $256$ repetitions. This measure of the approximation quality of the estimated permutations corresponds to the upper bound used in the proof of Proposition 14 (see (5.98)) and is applicable since by construction, $\theta^{\ast}$ has row and column sums equal to zero. It is the dominating part in the error analysis, leading to the rate reported in Theorem 4, and it allows us to study a larger range of $n$, avoiding the need for subsequent projection of the permuted $y$ matrix.

In Figure 3(a), we observe that while for smaller $n$, we see a slower decay than predicted, for larger $n$, the decay scales like $n^{-0.47}\approx n^{-1/2}$, close to the predicted rate.

Finally, we perform singular value thresholding on the same set of ground truth matrices, this time setting $V=1$, $\sigma=0.1$, and varying $n$ between $20$ and $500$. For this experiment, in Figure 3(b), we plotted the full denoising error,

[TABLE]

averaged over $64$ repetitions. As in the other experiments, we can see an error decay that is close to our theoretical guarantees, that is, $n^{-0.73}\approx n^{-3/4}$.

5 Proofs

In this section, we provide the proofs of our results. Recall that $D$ is defined in (2.2) and $\tilde{D}$ is defined analogously for dimension $n_{2}$. In the sequel, whenever we introduce notation in dimension $n_{1}$, the analogous object in dimension $n_{2}$ is denoted by the same symbol with a tilde.

5.1 A structural lemma

To gain insight into the set of anti-Monge matrices, we first state and prove the following structural lemma, which says that any anti-Monge matrix can be written as the sum of a constant-row matrix, a constant-column matrix, and an anti-Monge, bivariate isotonic matrix. Several structural decompositions of this type are known, e.g., Lemma 2.1 of [11] and Lemma 2.5 of [29]. Our lemma is stated in a form that is convenient for application, and we provide a self-contained proof which also facilitates understanding the structure of an anti-Monge matrix.

Lemma 7**.**

For each $\theta\in\mathcal{M}^{n_{1},n_{2}}$, there exists a unique triple $(R,S,B)$ of matrices in $\mathbbm{R}^{n_{1}\times n_{2}}$ such that

• •

$\theta=R+S+B$;

• •

$R\tilde{D}^{\top}=0$* and $DS=0$, i.e., $R$ has constant rows and $S$ has constant columns;*

• •

$S_{i,1}=B_{i,1}=B_{1,j}=0$* for $i\in[n_{1}],\,j\in[n_{2}].$*

Moreover, we have that

• •

$B$* is anti-Monge;*

• •

$B$* is bivariate isotonic, i.e., $B$ has nondecreasing rows and columns;*

• •

$\max\{\|R\|_{\infty},\|S\|_{\infty},\|B\|_{\infty}\}\leq 4\|\theta\|_{\infty}$;

• •

$\|B\|_{\infty}=B_{n_{1},n_{2}}=V(\theta)$.

Proof.

By the condition $S_{i,1}=B_{i,1}=0$ for $i\in[n_{1}]$, the first column of $R$ is equal to that of $\theta$. Since $R$ has constant rows, it is uniquely defined by $R=\theta_{\cdot,1}\tilde{\mathbf{1}}^{\top}$ where $\tilde{\mathbf{1}}$ denotes the all-ones vector in $\mathbbm{R}^{n_{2}}$. Moreover, by the condition $B_{1,j}=0$ for $j\in[n_{2}]$, the first row of $S$ is equal to that of $\theta-R$. Since $S$ has constant columns, it is uniquely defined by $S=\mathbf{1}(\theta-R)_{1,\cdot}$ where $\mathbf{1}$ denotes the all-ones vector in $\mathbbm{R}^{n_{1}}$. Finally, define $B=\theta-R-S.$

Since $DS=R\tilde{D}^{\top}=0$, we have $DB\tilde{D}^{\top}=D(\theta-R-S)\tilde{D}^{\top}=D\theta\tilde{D}^{\top}\geq 0$. so that $B$ is anti-Monge. Moreover, this implies

[TABLE]

The first column and first row of $B$ are equal to [math] by construction. Consequently, we obtain that $(DB)_{\cdot,1}=0$, and by (5.27), this gives $DB\geq 0$ and similarly $\theta\tilde{D}^{\top}\geq 0$ by (5.28) . Together, these facts yield that $B$ is bivariate isotonic.

Additionally, from the triangle inequality and the way we constructed $R,S,$ and $B$, we get

[TABLE]

Since $B$ also inherits the variation of $\theta$ and is nonnegative, $\|B\|_{\infty}=V(B)=V(\theta)$, which completes the proof. ∎

5.2 Proof of Theorem 1

To control the performance of a least-squares estimator, we employ Chatterjee’s variational formula [16] that we recall below. See, e.g., Lemma 6.1 of [30] for this deterministic form.

Lemma 8** (Chatterjee’s variational formula).**

Let $\mathcal{M}$ be a closed subset of ${\mathbbm{R}}^{d}$. Suppose that $y=\theta^{*}+\varepsilon$ where $\theta^{*}\in\mathcal{M}$ and $\varepsilon\in{\mathbbm{R}}^{d}$. Let $\hat{\theta}\in\operatorname*{argmin}_{\theta\in\mathcal{M}}\|y-\theta\|_{2}^{2}$ be a projection of $y$ onto $\mathcal{M}$. Define the function $f_{\theta^{*}}:{\mathbbm{R}}_{+}\to{\mathbbm{R}}$ by

[TABLE]

Then we have

[TABLE]

Moreover, if there exists $t^{*}>0$ such that $f_{\theta^{*}}(t)<0$ for all $t\geq t^{*}$, then $\|\hat{\theta}-\theta^{*}\|_{2}\leq t^{*}.$

To control the supremum in Lemma 8, note that it suffices to consider Gaussian noise here, since the generalization to sub-Gaussian noise is taken care of by Theorem 20.

Proposition 9**.**

*Fix an anti-Monge matrix $\theta^{*}\in\mathcal{M}$, and suppose that $Z\in\mathbbm{R}^{n_{1}\times n_{2}}$ has i.i.d. $\mathcal{N}(0,1)$ entries. Then for any integer $k\in[n_{1}n_{2}]$ and any $t>0$, we have *

[TABLE]

To show Theorem 1 taking Proposition 9 as given, let $t>0$ and $1\leq k\leq n_{1}n_{2}$ to be chosen later. Note that by Theorem 19 and Proposition 9, we obtain

[TABLE]

where $\gamma_{2}$ denotes Talagrand’s $\gamma_{2}$ functional. Therefore, Theorem 20 yields that with probability $1-4\exp(-s^{2})$,

[TABLE]

Let us define

[TABLE]

we obtain that for any

[TABLE]

where $C$ is a sufficiently large constant, it holds with probability at least $1-4\exp(-s^{2})$ that

[TABLE]

Therefore by Lemma 8, we obtain that with probability at least $1-4\exp(-s^{2})$,

[TABLE]

We now choose $s=2\sqrt{n_{1}}$. Balancing the terms that depend on $k$ leads to the choice

[TABLE]

in addition to possibly rounding $k^{\ast}$ to an integer which we omitted to simplify the presentation. Therefore, we obtain that with probability $1-\exp(-n_{1})$, if $1\leq k^{\ast}\leq n_{1}n_{2}$, then

[TABLE]

If $k^{\ast}<1$, we replace it by $1$, increasing the $k\log(n_{2})/(n_{1}n_{2})$ term while decreasing the ones with $1/\sqrt{k}$ in (5.35), hence leading to the same rate as in (5.37). If $k^{\ast}>n_{1}n_{2}$, note that the $k/(n_{1}n_{2})$ term is already of the order $\sigma^{2}$, so a basic bound of $\sigma^{2}$ on the empirical process term yields the rate

[TABLE]

Combined, this yields that with probability at least $1-\exp(-n_{1})$,

[TABLE]

To obtain the bound in expectation, we can first integrate the exponentially decaying tale of (5.35), and then choose the optimal $k$ in the same way.

5.3 Proof of Proposition 9

Our main strategy consists in decomposing the noise matrix $Z$ into three terms according to the spectral decomposition of the linear map $\mathcal{D}$, defined by $\mathcal{D}(A)=DA\tilde{D}^{\top}$ for $A\in\mathbbm{R}^{n_{1}\times n_{2}}$.

Spectral decomposition of the difference operator.

Denote the (reduced) singular value decomposition of $D$ by $D=U\Sigma W^{\top}$, where we order the singular values in $\Sigma$ in ascending magnitude. In addition, we write W=\left[\begin{array}[]{@{}c|c|c@{}}w_{1}&\cdots&w_{n_{1}-1}\end{array}\right].

First, let $\Pi_{1}$ denote the projection onto $\ker\mathcal{D}$. Moreover, let $J=\{(l,r)\in[n_{1}]\times[n_{2}]:lr\leq k\}$ and $J^{c}=[n_{1}]\times[n_{2}]\setminus J$. Define the projection $\Pi_{2}$ by

[TABLE]

With these two projections, we decompose

[TABLE]

We now bound the three terms in (5.42) separately.

Bounding the first term in (5.42).

Recall that $\Pi_{1}$ be the projection onto $\ker\mathcal{D}$. We claim that $\dim(\ker\mathcal{D})=n_{1}+n_{2}-1$. Given a matrix $\theta\in\ker\mathcal{D}$, i.e., $D\theta\tilde{D}^{\top}=0$, we apply Lemma 7 to obtain the unique decomposition $\theta=R+S+B$, where $R\tilde{D}^{\top}=0$ and $DS=0$. It follows that $DB\tilde{D}^{\top}=0$. Since the first column and the first row of $B$ are both identically zero, it is easy to see from an inductive argument that $B=0$ so that $\ker\mathcal{D}$ contains only matrices of the form $\theta=R+S$. The set of constant-row matrices $R$ has dimension $n_{1}$; the set of constant-column matrices $S$ with $S_{i,1}=0$ for $i\in[n_{1}]$ has dimension $n_{2}-1$. Thus $\dim(\ker\mathcal{D})=n_{1}+n_{2}-1$. Consequently, we have

[TABLE]

Bounding the second term in (5.42).

Similarly, it suffices to compute the rank of $\Pi_{2}$, which is bounded as follows

[TABLE]

Therefore, we obtain

[TABLE]

Bounding the third term in (5.42).

Note that $I-\Pi_{1}$ is the projection onto the image of the linear map $\mathcal{D}^{\top}$, defined by $\mathcal{D}^{\top}(A)=D^{\top}A\tilde{D}$. Hence we have

[TABLE]

Since $Z$ has mean zero, it is sufficient to control

[TABLE]

To bound this quantity, we need the following lemma, whose proof is deferred to Section 5.4.

Lemma 10**.**

For any $i\in[n_{1}],\,j\in[n_{2}]$, the quantity $[(D^{\dagger})^{\top}(I-\Pi_{2})(Z)(\tilde{D}^{\dagger})]_{i,j}$ is sub-Gaussian with variance proxy

[TABLE]

Let us define $\Phi\in\mathbbm{R}^{(n_{1}-1)\times(n_{2}-1)}$ by

[TABLE]

and let $\oslash$ denote element-wise division. Lemma 10, together with a union bound readily yields

[TABLE]

In addition, it holds for every $\theta$ that

[TABLE]

by Hölder’s inequality. We therefore obtain

[TABLE]

It remains to bound this supremum. For $\theta\in\mathcal{M}$, because $D\theta\tilde{D}^{\top}\geq 0$ we can write

[TABLE]

The second term in (5.53) can be bounded by

[TABLE]

For the first term in (5.53), we need the following lemma, whose proof is deferred to Section 5.5.

Lemma 11**.**

We have the estimate

[TABLE]

If $\|\theta-\theta^{*}\|_{F}\leq t$, then the above lemma together with the Cauchy-Schwarz inequality yields

[TABLE]

Combining (5.47)–(5.55), we conclude that

[TABLE]

The bounds on the three terms of (5.42) together yield the desired result.

5.4 Proof of Lemma 10

By definition of $\Pi_{2}$, it holds that

[TABLE]

We now study the sub-Gaussianity of the $(i,j)$-th entry of this quantity. Since $Z$ has i.i.d. $\mathcal{N}(0,1)$ entries, it holds for each $\lambda>0$ that

[TABLE]

Note that $\|\tilde{w}_{r}w_{\ell}^{\top}\|_{F}=1$, and $\langle\tilde{w}_{r}w_{\ell}^{\top},\tilde{w}_{r^{\prime}}w_{\ell^{\prime}}^{\top}\rangle=0$ for any pairs $(r,\ell)\neq(r^{\prime},\ell^{\prime})$, so we have

[TABLE]

It remains to bound this quantity. Without loss of generality, assume that $n_{1}$ is odd, so $n_{1}-1$ is even. The matrix $D$ has the same left-singular vectors as

[TABLE]

which are known [65] to be

[TABLE]

Moreover, the matrix $D$ has (non-zero) singular values

[TABLE]

Note that because of the symmetry

[TABLE]

it is enough to consider $i=1,\dots,\frac{n_{1}-1}{2}$. We make use of the following inequalities to control the $\sin$ terms involved:

[TABLE]

Plugging in the entries of $U$ and $\Sigma$ yields

[TABLE]

where we used (5.69) on the numerator and (5.71) on the denominator in $(i)$ and the bound $\sum_{r=k+1}^{\infty}\frac{1}{r^{2}}\leq\frac{1}{k}$ for any $k\geq 1$ in $(ii)$ and $(iii)$.

On the other hand, even without using the constraint $(\ell,r)\in J^{c}$, we have

[TABLE]

where in $(i)$, we used (5.70) for the numerator and (5.71) for the denominator as well as (5.76) with $k$ replaced by $\frac{n_{1}n_{2}}{ij}$, and $(ii)$ follows from counting integer points in the set $\{(\ell,r):\ell r\leq\frac{n_{1}n_{2}}{ij}\}$ as in (5.44).

A similar argument yields bounds with $i$ replaced by $n_{1}-i$, or $j$ replaced by $n_{2}-j$. Combining this observation with (5.61), (5.63), (5.64), (5.76) and (5.77) completes the proof.

5.5 Proof of Lemma 11

Define $\Phi^{\prime}=\frac{\Phi}{\sqrt{\log(n_{1})\log(n_{2})}}$, i.e.,

[TABLE]

To simplify the notation, let $\Phi^{\prime}_{i,0}=\Phi^{\prime}_{i,n_{2}}=\Phi^{\prime}_{0,j}=\Phi^{\prime}_{n_{1},j}=0$ for all $i\in[n_{1}],\,j\in[n_{2}]$. We need to bound

[TABLE]

By symmetry, it suffices to consider summation over $i\in[\frac{n_{1}-1}{2}],\,j\in[\frac{n_{2}-1}{2}]$ where $\Phi^{\prime}=\sqrt{ij}\land\sqrt{\frac{n_{1}n_{2}}{k}}$. Moreover, note that the summand vanishes if $(i-1)(j-1)\geq\frac{n_{1}n_{2}}{k}$. Hence we can further split the sum into two parts:

over $\{i\in[\frac{n_{1}-1}{2}],\,j\in[\frac{n_{2}-1}{2}]:(i-1)(j-1)<\frac{n_{1}n_{2}}{k}<ij\}$, and 2. 2.

over $\{i\in[\frac{n_{1}-1}{2}],\,j\in[\frac{n_{2}-1}{2}]:ij\leq\frac{n_{1}n_{2}}{k}\}$.

To bound the first part of the sum, first consider the case $i\leq j$, where we have

[TABLE]

Adjusting signs to ensure both differences in the following expression are positive, it is easily checked that

[TABLE]

By the conditions $(i-1)(j-1)<\frac{n_{1}n_{2}}{k}<ij$ and $i\leq j$, we obtain

[TABLE]

Therefore, it remains to bound

[TABLE]

An analogous argument yields the same bound for the case $i>j$.

Next, we consider the sum of \big{(}\Phi^{\prime}_{i-1,j-1}+\Phi^{\prime}_{i,j}-\Phi^{\prime}_{i-1,j}-\Phi^{\prime}_{i,j-1}\big{)}^{2} over $(i,j)$ such that $ij\leq\frac{n_{1}n_{2}}{k}$, where we have

[TABLE]

Now even summing over all indices $i\in[n_{1}],\,j\in[n_{2}]$ yields

[TABLE]

Combining all the pieces, we obtain

[TABLE]

5.6 Proof of Theorem 2

First, consider the set $\mathcal{C}$ of constant-row matrices, which is a subset of $\mathcal{M}$. Note that $V(\theta^{*})=0$ for any $\theta^{*}\in\mathcal{C}$. Since there are $n_{1}$ rows, it is not hard to see that the minimax rate of estimation over $\mathcal{C}$ is $\sigma^{2}\,n_{1}$ in the squared Frobenius norm. The lower bound of order $\sigma^{2}/n_{2}$ then follows from normalization by $n_{1}n_{2}$.

Next, we turn to the second term in the lower bound, which is based on Assouad’s Lemma (Theorem 21). To this end, we construct an embedding of the hypercube into $\mathcal{M}_{V_{0}}$.

Consider integers $k_{1}\in[n_{1}]$ and $k_{2}\in[n_{2}]$, and let $m_{1}=n_{1}/k_{1}$ and $m_{2}=n_{2}/k_{2}$. Assume without loss of generality that $m_{1}$ and $m_{2}$ are integer. Denote the elements of the hypercube $\{-1,1\}^{k_{1}\times k_{2}}$ by $(\tau_{u,v}:(u,v)\in[k_{1}]\times[k_{2}])$. For each $\tau\in\{-1,1\}^{k_{1}\times k_{2}}$, define $\theta^{\tau}\in\mathcal{M}$ in the following way. For $i\in[n_{1}]$ and $j\in[n_{2}]$, first identify the unique $u\in[k_{1}],v\in[k_{2}]$ for which $(u-1)m_{1}<i\leq um_{1}$ and $(v-1)m_{2}<j\leq vm_{2}$, and then take

[TABLE]

First, we check that $\theta^{\tau}\in\mathcal{M}$ and that $V(\theta^{\tau})\leq V_{0}$. For the former, it is enough to check that for each $1\leq i<n_{1},1\leq j<n_{2}$, we have that $\theta^{\tau}_{i,j}+\theta^{\tau}_{i+1,j+1}-\theta^{\tau}_{i,j+1}-\theta^{\tau}_{i+1,j}\geq 0$. We distinguish the following two cases:

There exists $u\in[k_{1}]$ such that $(u-1)m_{1}<i<i+1\leq um_{1}$, or there exists $v\in[k_{2}]$ such that $(v-1)m_{2}<j<j+1\leq vm_{2}$. Then, either $\theta^{\tau}_{i,j}=\theta^{\tau}_{i+1,j},\,\theta^{\tau}_{i,j+1}=\theta^{\tau}_{i+1,j+1}$ or $\theta^{\tau}_{i,j}=\theta^{\tau}_{i,j+1},\,\theta^{\tau}_{i+1,j}=\theta^{\tau}_{i+1,j+1}$ respectively. In both cases, the difference above is [math]. 2. 2.

There exist $u\in[k_{1}],v\in[k_{2}]$ such that $(u-1)m_{1}<i\leq um_{1}<i+1\leq(u+1)m_{1}$ and $(v-1)m_{2}<j\leq vm_{2}<j+1\leq(v+1)m_{2}$. In this case, for any $\tau$, we have

[TABLE]

Thus, $\theta^{\tau}\in\mathcal{M}$. We now check that $V(\theta^{\tau})\leq V_{0}$. Note that $V(\theta^{\tau})$ can be written as the sum $\sum_{i,j}(\theta^{\tau}_{i,j}+\theta^{\tau}_{i+1,j+1}-\theta^{\tau}_{i,j+1}-\theta^{\tau}_{i+1,j})$. As we have seen above, this sum is nonzero in only $(k_{1}-1)(k_{2}-1)$ cases, and it equals exactly

[TABLE]

We now proceed to show our lower bound using Theorem 21, which states that

[TABLE]

where $d_{H}$ denotes the Hamming distance and $\ell^{2}(\theta,\theta^{\prime})=\frac{1}{n_{1}n_{2}}\|\theta-\theta^{\prime}\|_{F}^{2}$. Note that

[TABLE]

Thus, we have

[TABLE]

To bound $\|\mathbb{P}_{\theta^{\tau}}-\mathbb{P}_{\theta^{\tau^{\prime}}}\|_{TV}$, we use Pinsker’s inequality:

[TABLE]

It follows that

[TABLE]

Putting things together in Assouad’s lemma, we obtain

[TABLE]

If $\frac{4\sigma}{\sqrt{n_{1}n_{2}}}\leq V_{0}\leq 4\sigma n_{1}n_{2}$, then we can choose $k_{1}$ and $k_{2}$ such that $k_{1}k_{2}$ is of order \big{(}\frac{V_{0}\sqrt{n_{1}n_{2}}}{4\sigma}\big{)}^{\frac{2}{3}}, which yields that

[TABLE]

If $V_{0}\leq\frac{4\sigma}{\sqrt{n_{1}n_{2}}}$, then \big{(}\frac{\sigma^{2}V_{0}}{n_{1}n_{2}}\big{)}^{2/3}\lesssim\frac{\sigma^{2}}{n_{1}n_{2}}, so the second term in the statement of Theorem 2 is dominated by the first term. Finally, if $V_{0}\geq 4\sigma n_{1}n_{2}$, then \big{(}\frac{\sigma^{2}V_{0}}{n_{1}n_{2}}\big{)}^{2/3}\geq\sigma^{2}, so the rate becomes trivial. This completes the proof.

5.7 Proof of Theorem 3

Similar to the proof of Theorem 1, the main technical step in this proof is to bound the supremum of an empirical process. This is dealt with in the following proposition, whose proof is deferred to Section 5.8. Note that both the statement and the proof are very similar to Proposition 9, and that we can restrict our attention to noise matrices that are Gaussian.

Proposition 12**.**

Fix any matrix $\theta^{*}\in\mathbbm{R}^{n_{1}\times n_{2}}$ and permutations $\pi_{1}\in\mathcal{S}_{n_{1}}$ and $\pi_{2}\in\mathcal{S}_{n_{2}}$. Suppose that $Z\in\mathbbm{R}^{n_{1}\times n_{2}}$ has i.i.d. $\mathcal{N}(0,1)$ entries. Then for any integer $k\in[n_{1}n_{2}]$ and any $t>0$, we have

[TABLE]

where we use the convention that the supremum over the empty set is $-\infty$.

We assume without loss of generality that the underlying true permutations $\pi_{1}^{*}$ and $\pi_{2}^{*}$ are the identities throughout the proof. For fixed permutations $\pi_{1}\in\mathcal{S}_{n_{1}}$ and $\pi_{2}\in\mathcal{S}_{n_{2}}$, we obtain from Theorem 19 and Proposition 12 that

[TABLE]

Therefore, Theorem 20 yields that with probability $1-4\exp(-s^{2})$,

[TABLE]

Taking $s=2\sqrt{n_{1}\log(n_{1})}$ and applying a union bound over all $(\pi_{1},\pi_{2})\in\mathcal{S}_{n_{1}}\times\mathcal{S}_{n_{2}}$ (which has log-cardinality $\log(n_{1}!n_{2}!)\leq 2n_{1}\log n_{1}$), we get that with probability at least $1-n_{1}^{-n_{1}}$,

[TABLE]

Let us define

[TABLE]

Then for any

[TABLE]

where $C$ is a sufficiently large constant, it holds with probability at least $1-n_{1}^{-n_{1}}$ that

[TABLE]

Therefore by Lemma 8, we obtain

[TABLE]

noting that by assumption, $\theta^{\ast}\in\bigcup_{\begin{subarray}{c}\scriptscriptstyle\pi_{1}\in\mathcal{S}_{n_{1}}\\ \scriptscriptstyle\pi_{2}\in\mathcal{S}_{n_{2}}\end{subarray}}\mathcal{M}_{V_{0}}(\pi_{1},\pi_{2})$.

Balancing the terms that depend on $k$ leads to the choice

[TABLE]

and therefore we obtain that with probability $1-n_{1}^{-n_{1}}$,

[TABLE]

if $1\leq k^{\ast}\leq n_{1}n_{2}$. We conclude by arguing similar to the proof of Theorem 1 to handle the other possible cases of $k^{\ast}$ and to get an error bound in expectation instead of with high probability.

5.8 Proof of Proposition 12

Note that Proposition 12 is very similar to Proposition 9, and so are their proofs. The difference is that Proposition 12 has extra complications arising from the presence of permutations, while at the same time it is simpler because we restrict $\theta$ to $\mathcal{M}_{V_{0}}\subset\mathcal{M}$. Hence, we focus only the differences to the proof of Proposition 9 here.

Since $Z$ is equal in distribution to $Z(\pi_{1}^{-1},\pi_{2}^{-1})$, we have

[TABLE]

in distribution for any permutations $\pi_{1}\in\mathcal{S}_{n_{1}},\pi_{2}\in\mathcal{S}_{n_{2}}$. Therefore, by replacing $\theta^{*}(\pi_{1}^{-1},\pi_{2}^{-1})$ with $\theta^{*}$, it suffices to prove that for any matrix $\theta^{*}\in\mathbbm{R}^{n_{1}\times n_{2}}$, it holds that

[TABLE]

Note that this supremum is very similar to that studied in Proposition 9, with the main differences being that $\theta^{*}$ can be any matrix in $\mathbbm{R}^{n_{1}\times n_{2}}$ while $\theta$ is restricted to $\mathcal{M}_{V_{0}}$. As in the proof of Proposition 9, (5.42), the supremum can be split into three terms:

[TABLE]

The first two terms can be bounded exactly as before, because we only need the condition $\|\theta-\theta^{*}\|_{F}\leq t$ but not any other property of $\theta^{*}$. Up to a constant factor, the third term can be bounded by (recall (5.47) and (5.52) in the proof of Proposition 9)

[TABLE]

Thanks to the constraint $\theta\in\mathcal{M}_{V_{0}}$ in this case, we immediately obtain

[TABLE]

Hence, it holds

[TABLE]

Combining the bounds on the three terms completes the proof.

5.9 Proof of Theorem 4

We assume without loss of generality that the underlying true permutations $\pi_{1}^{*}$ and $\pi_{2}^{*}$ are the identities throughout the proof, except where the notations $\pi_{1}^{\ast}$ and $\pi_{2}^{\ast}$ are explicitly used. Recall that we defined the reversal permutation $\pi_{1}^{\mathsf{r}}\in\mathcal{S}_{n_{1}}$ by $\pi_{1}^{\mathsf{r}}(i)=n_{1}-i+1$ for $i\in[n_{1}]$. Given estimators $\hat{\pi}_{1},\hat{\pi}_{2}$, let us define

[TABLE]

The theorem follows from the next two propositions combined. The first proposition says that the final denoising error can be bounded by the sum of the minimax rate (the error rate incurred by the projection step of the algorithm), and the error incurred by the permutation estimators. The second proposition controls the error incurred by the permutation estimators.

Proposition 13**.**

Suppose that we have $y=\theta^{*}+\varepsilon$, where the noise matrix $\varepsilon$ has independent $\mathsf{subG}(C\sigma^{2})$ entries with $\mathrm{Var}[\varepsilon_{i,j}]=\sigma^{2}.$ Let the estimators $(\hat{\pi}_{1},\hat{\pi}_{1}^{\prime},\hat{\pi}_{2},\hat{\theta})$ be given by Algorithm 2, and define $\tilde{\theta}$ according to (5.88). Then it holds with probability at least $1-n_{1}^{-n_{1}}$ that

[TABLE]

Proposition 14**.**

Suppose that we have $y=\theta^{*}+\varepsilon$, where $\theta^{*}\in\mathcal{M}_{V_{0}}^{n_{1},n_{2}}$ and the noise matrix $\varepsilon$ has independent $\mathsf{subG}(C\sigma^{2})$ entries. For the permutation estimators $\hat{\pi}_{1}$ and $\hat{\pi}_{2}$ given by the Variance Sorting subroutine, Algorithm 1, let $\tilde{\theta}$ be defined by (5.88). Then it holds with probability at least $1-n_{1}^{-9}$ that

[TABLE]

The error bound in Proposition 14 is clearly dominating, and combining the two propositions yields the statement of the theorem in probability. Taking into account that we can do the same analysis keeping track of the failure probability independently of $n_{1}$ and integrate, we obtain bounds in expectation instead of probability as well.

5.10 Proof of Proposition 13

To employ the variational formula in Lemma 8, let us view $\tilde{\theta}$ defined by (5.88) as the ground truth, and view $y-\tilde{\theta}=\varepsilon+\theta^{*}-\tilde{\theta}$ as the noise. Correspondingly, we define

[TABLE]

To facilitate our analysis, for each pair of permutations $(\pi_{1},\pi_{2})\in\mathcal{S}_{n_{1}}\times\mathcal{S}_{n_{2}}$, we define

[TABLE]

and note that we have that either $\tilde{\theta}=\tilde{\theta}^{\hat{\pi}_{1},\hat{\pi}_{2}}$ or $\tilde{\theta}=\tilde{\theta}^{\pi_{1}^{\mathsf{r}}\circ\hat{\pi}_{1},\hat{\pi}_{2}}$, so $\tilde{\theta}\in\mathcal{M}_{V_{0}}(\hat{\pi}_{1},\hat{\pi}_{2})\cup\mathcal{M}_{V_{0}}(\pi_{1}^{\mathsf{r}}\circ\hat{\pi}_{1},\hat{\pi}_{2})$ and Lemma 8 is applicable.

We further estimate the supremum in (5.90) by

[TABLE]

Note that the random variables $\hat{\pi}_{1},\hat{\pi}_{2}$ and $\tilde{\theta}$ depend on $\varepsilon$, so it is not clear how to control the first supremum in (5.91). Instead, in (5.92) we take a supremum over all $(\pi_{1},\pi_{1}^{\prime},\pi_{2},\pi_{2}^{\prime},\theta)\in\mathcal{S}_{n_{1}}^{2}\times\mathcal{S}_{n_{2}}^{2}\times\mathcal{M}_{V_{0}}$, where each individual quantity $\theta(\pi_{1}^{\prime},\pi_{2}^{\prime})-\tilde{\theta}^{\pi_{1},\pi_{2}}$ is deterministic.

For fixed permutations $\pi_{1},\pi_{1}^{\prime}\in\mathcal{S}_{n_{1}}$ and $\pi_{2},\pi_{2}^{\prime}\in\mathcal{S}_{n_{2}}$, we obtain from Theorem 19 and Proposition 12 that

[TABLE]

Therefore, Theorem 20 yields that with probability $1-4\exp(-s^{2})$,

[TABLE]

Taking $s=3\sqrt{n_{1}\log(n_{1})}$ and applying a union bound over all $\pi_{1},\pi_{1}^{\prime}\in\mathcal{S}_{n_{1}}$ and $\pi_{2},\pi_{2}^{\prime}\in\mathcal{S}_{n_{2}}$, we see that with probability at least $1-n_{1}^{-n_{1}}$,

[TABLE]

This together with inequality (5.92) yields that for any

[TABLE]

where $C$ is a sufficiently large constant, it holds with probability at least $1-n_{1}^{-n_{1}}$ that

[TABLE]

Therefore, by Lemma 8 applied to the set $\mathcal{M}_{V_{0}}(\hat{\pi}_{1},\hat{\pi}_{2})\cup\mathcal{M}_{V_{0}}(\pi_{1}^{\mathsf{r}}\circ\hat{\pi}_{1},\hat{\pi}_{2})$, we obtain

[TABLE]

Balancing the terms that depend on $k$ yields that with probability $1-n_{1}^{-n_{1}}$,

[TABLE]

if for the optimal $k^{\ast}$, we have $1\leq k^{\ast}\leq n_{1}n_{2}$. The other cases can be handled as in the proof of Theorem 1.

5.11 Proof of Proposition 14

5.11.1 Reduction to the row/column-centered case.

First, we reduce the problem to the case where the underlying matrix $\theta^{*}$ has centered rows and columns. If $\theta^{*}$ is not centered, we may choose a matrix $R\in\mathbbm{R}^{n_{1}\times n_{2}}$ with constant rows and a matrix $S\in\mathbbm{R}^{n_{1}\times n_{2}}$ with constant columns, so that if $\bar{\theta}:=\theta^{*}-R-S$, then for all $i\in[n_{1}],\,j\in[n_{2}]$,

[TABLE]

Since $R$ and $S$ have constant rows and columns respectively, we have $V(\bar{\theta})=V(\theta^{*})$ by definition (2.4), and thus $\bar{\theta}\in\mathcal{M}_{V_{0}}$. More importantly, according to the definition of $\xi(i,j)$ in (3.16), its value does not change if we replace $y$ by $y-R-S$. Therefore, we may assume without loss of generality that

[TABLE]

which does not change the estimators $\hat{\pi}_{1}$, $\hat{\pi}_{1}^{\prime}$ and $\hat{\pi}_{2}$ output by the algorithm.

Furthermore, we have

[TABLE]

since if $\tilde{\theta}$ minimizes the left-hand side, then $\tilde{\theta}-R-S$ minimizes the right-hand side. Hence it suffices to show that the right-hand side of (5.95) is bounded by the desired rate.

Note that by symmetry of the anti-Monge constraint,

[TABLE]

In light of this, it is sufficient to show

[TABLE]

which we do in the sequel. This gives an upper bound on (5.95) and thus completes the proof.

5.11.2 Preliminaries.

Before proceeding to proving (5.98), we start with some lemmas.

Lemma 15**.**

For $\gamma\in\mathbbm{R}^{n}$ such that $\sum_{k=1}^{n}\gamma_{k}=0$, it holds that

[TABLE]

The proof follows by inspection.

Lemma 16**.**

Let $f:[n]\times[n]\to\mathbbm{R}$ be a symmetric bivariate function such that

[TABLE]

Let $\pi:[n]\to[n]$ be a permutation and $\tau\in\mathbbm{R}$. Suppose that

[TABLE]

Then we have $f(\pi(i),i)\leq\tau$ for all $i\in[n]$.

Proof.

Suppose that $f(\pi(j),j)>\tau$ and $\pi(j)<j$ for some index $j\in[n]$. Since $\pi$ is a bijection, there must exist an index $i\leq\pi(j)<j$ such that $\pi(i)>\pi(j)$. However, by (5.99) we then have

[TABLE]

which contradicts assumption (5.100). A similar argument yields a contradiction in the case that $f(\pi(j),j)>\tau$ and $\pi(j)>j$. Therefore, we obtain that $f(\pi(j),j)\leq\tau$ for all $j\in[n]$. ∎

Next, we study the quantity $\xi(i,j)$ used in the algorithm defined by (3.16). Throughout the rest of the proof, we use the notation

[TABLE]

Lemma 17**.**

Suppose that $y=\bar{\theta}+\varepsilon$, where $\bar{\theta}$ has centered rows and columns, and $\varepsilon$ has independent $\mathsf{subG}(C\sigma^{2})$ entries with $\mathrm{Var}[\varepsilon_{i,j}]=\sigma^{2}$. Then it holds that for all distinct $i,j\in[n_{1}]$,

[TABLE]

and that with probability $1-n_{1}^{-10}$, for all $i,j\in[n_{1}]$,

[TABLE]

where $C>0$ is a universal constant and $V(\bar{\theta})$ is defined as in (2.4).

The proof of the lemma is deferred to Section 5.12.

Next, we study properties of the expectation $\mathbbm{E}[\xi(i,j)]$, or equivalently $f(i,j)$, thanks to (5.102).

Lemma 18**.**

It holds for all $1\leq i\leq m\leq j\leq n_{1}$ that $f(i,m)+f(m,j)\leq f(i,j).$

Proof.

Fix indices $1\leq i\leq m\leq j\leq n_{1}$. Since all the row sums of $\bar{\theta}$ are zero by assumption, it follows from Lemma 15 that

[TABLE]

Note that we have

[TABLE]

where each of the three bracketed terms is nonnegative because $\bar{\theta}$ is anti-Monge. Therefore, we obtain

[TABLE]

This together with (5.104) completes the proof. ∎

5.11.3 Proof of main bound (5.98).

We condition on the event of probability $1-n_{1}^{-10}$ that (5.103) holds. Consider the permutation estimator $\hat{\pi}_{1}$ defined as in the algorithm so that $\{\xi(i_{0},\hat{\pi}_{1}^{-1}(i))\}_{i=1}^{n_{1}}$ is nondecreasing. Note that we made the decision of whether to consider $i_{0}$ or $j_{0}$ as an estimator for $\hat{\pi}_{1}^{-1}(1)$ arbitrarily by demanding $i_{0}<j_{0}$. In the following, we make use of the fact that this orientation aligns with the assumption of $\pi_{1}^{\ast}=\mathsf{id}$, that is, we use $\pi_{1}^{\ast}(i_{0})<\pi_{1}^{\ast}(j_{0})$. If the reverse inequality holds, then we may repeat the same proof with $\hat{\pi}_{1}$ replaced by $\pi_{1}^{\mathsf{r}}\circ\hat{\pi}_{1}$, to obtain permutation guarantees for the reversed permutation instead.

We claim that for $f$ defined in (5.101),

[TABLE]

To establish the claim, we first consider any pair $(i,j)$ for which $i_{0}<i<j$ and $\hat{\pi}_{1}(i)>\hat{\pi}_{1}(j)$. Thus by the definition of $\hat{\pi}_{1}$, we have $\xi(i_{0},i)>\xi(i_{0},j)$. Then it follows from Lemma 18, (5.102) and (5.103) that

[TABLE]

Next, consider any pair $(i,j)$ where $i\leq i_{0}$ and $\hat{\pi}_{1}(i)>\hat{\pi}_{1}(j)$. By the definition of $(i_{0},j_{0})$ in (3.17), we have that $i_{0}<j_{0}$ and $\xi(i_{0},j_{0})\geq\xi(i,j_{0})$. Together with Lemma 18, (5.102) and (5.103), this implies that

[TABLE]

Moreover, we have $\xi(i,i_{0})\geq\xi(i_{0},j)$ since $\hat{\pi}_{1}(i)>\hat{\pi}_{1}(j)$. Therefore, by (5.102) and (5.103) it holds

[TABLE]

Combining (5.106) and (5.107), we obtain

[TABLE]

Therefore, claim (5.105) is established.

Note that assumption (5.99) of Lemma 16 holds in view of Lemma 18 and the fact that $f(i,j)\geq 0$. Hence claim (5.105) and Lemma 16 together yield that for all $i\in[n_{1}]$,

[TABLE]

Summing over $i$, we conclude that

[TABLE]

contingent on the assumption that for the ground truth permutation $\pi_{1}^{\ast}$, we have $(\pi_{1}^{\ast})^{-1}(i_{0})<(\pi_{1}^{\ast})^{-1}(j_{0})$. If not, repeat the same proof with $\hat{\pi}_{1}$ replaced by $\pi_{1}^{\mathsf{r}}\circ\hat{\pi}_{1}$ to obtain

[TABLE]

Finally, the proof remains valid if we replace $\bar{\theta}$ and $y$ by their transposes, and switch the roles of row and column indices. Hence it also holds with probability $1-n_{1}^{-10}$ that

[TABLE]

We then complete the proof of (5.98) by using the triangle inequality to see that for some choice of $\bar{\pi}_{1}\in\{\mathsf{id},\pi_{1}^{\mathsf{r}}\}$, $\bar{\pi}_{2}\in\{\mathsf{id},\pi_{2}^{\mathsf{r}}\}$, it holds that

[TABLE]

5.12 Proof of Lemma 17

Recall that $y=\bar{\theta}+\varepsilon$ where $\mathrm{Var}[\varepsilon_{i,j}]=\sigma^{2}$, and recall the notation

[TABLE]

For the statement about its expectation, (5.102), we need to prove that for distinct $i,j\in[n_{1}]$,

[TABLE]

For the statement about its deviation, (5.103), we claim that it suffices to prove that with probability $1-n_{1}^{-12}$,

[TABLE]

To see this, note that Lemma 15 gives

[TABLE]

where we used that $|\bar{\theta}_{i,k}+\bar{\theta}_{j,\ell}-\bar{\theta}_{j,k}-\bar{\theta}_{i,\ell}|\leq V_{0}$. Plugging this bound into (5.114) and applying a union bound over all $i,j\in[n_{1}]$ then completes the proof.

The claims (5.113) and (5.114) can be simplified as follows. For distinct $i,j\in[n_{1}]$, we let

[TABLE]

Note that $\delta$ has independent $\mathsf{subG}(C\sigma^{2})$ entries and $\mathbbm{E}[\delta_{k}^{2}]=2\sigma^{2}$. We need to prove that

[TABLE]

and that with probability $1-n_{1}^{-12}$,

[TABLE]

Recall that all the rows and columns of $\bar{\theta}$ are centered by assumption, so we have $\sum_{k=1}^{n}\gamma_{k}=0$. Using this, we get the following Hoeffding decomposition

[TABLE]

In particular, it follows that

[TABLE]

Moreover, we have that with probability $1-n_{1}^{-12}$,

[TABLE]

where concentration of the first two terms is due to Lemma 22, and the last term is $\mathsf{subG}(C\sigma^{2}\|\gamma\|_{2}^{2})$. This completes the proof.

5.13 Proof of Theorem 5

This proof technique was developed by [17, 62], and our presentation follows that of Theorem 2 of [62]. We assume without loss of generality that the underlying true permutations $\pi_{1}^{*}$ and $\pi_{2}^{*}$ are the identities throughout the proof. Let $\|\cdot\|$ denotes the operator norm of a matrix. It is well known (see Theorem 4.4.5 of [68]) that $\|\varepsilon\|\leq C_{1}\sigma\sqrt{n_{1}}$ with probability at least $1-\exp(-n_{1})$. We condition on this event in the sequel.

Recall that the singular value decomposition of $y$ is

[TABLE]

where $\lambda_{1},\dots,\lambda_{n_{2}}$ are ordered non-increasingly, and the SVT estimator is defined as

[TABLE]

where we choose $\rho=2C_{1}\sigma\sqrt{n_{1}}$. Moreover, write the singular value decomposition of $\theta^{*}$ as

[TABLE]

where $\lambda^{*}_{1},\dots,\lambda^{*}_{n_{2}}$ are ordered non-increasingly. Let $s$ be the number of singular values of $\theta^{*}$ that are larger than $C_{1}\sigma\sqrt{n_{1}}$, and define

[TABLE]

Note that the for each $i>s$, by Weyl’s inequality, we have

[TABLE]

Therefore, $\hat{\theta}^{\sf svt}$ has rank at most $s$, and so has $\theta_{s}$. It follows that

[TABLE]

Moreover, it holds that

[TABLE]

Plugging this bound into the previous one, we obtain

[TABLE]

For any integer $r\geq 6$, Proposition 6 yields a rank-$r$ matrix $\tilde{\theta}\in\mathbbm{R}^{n_{1}\times n_{2}}$ such that

[TABLE]

Since $\theta_{r}=\sum_{i=1}^{r}\lambda^{*}_{i}u^{*}_{i}(v^{*}_{i})^{\top}$ is by definition the best rank-$r$ approximation of $\theta^{*}$ in the Frobenius norm, we see that

[TABLE]

Hence it follows from (5.122) that

[TABLE]

Choosing the optimal $r^{\ast}$ and considering the boundary cases $r^{\ast}<1$ and $r^{\ast}>n_{2}$ then yields

[TABLE]

By repeating the same proof keeping track of the failure probability of the statement, we can obtain bounds in expectation as well.

5.14 Proof of Proposition 6

By rescaling, we may assume that $V(\theta)=1$ without loss of generality. Lemma 7 yields that $\theta=R+S+B$, where $R$ and $S$ are rank-one matrices, and $B$ is anti-Monge, bivariate isotonic (i.e., $B$ has nondecreasing rows and columns). Additionally, we have $B_{i,1}=B_{1,j}=0$ for $i\in[n_{1}],j\in[n_{2}]$, and $B_{n_{1},n_{2}}=1$. It suffices to find a low-rank approximation of the matrix $B$.

5.14.1 Subdivision.

We claim that there exist two increasing sequences of indices $\{i_{k}\}_{k=1}^{r+1}$ and $\{j_{\ell}\}_{\ell=1}^{2r}$ such that

• •

$0\leq i_{k}-i_{k-1}\leq n_{1}/r$ for $k\in[r+1]$;

• •

$0\leq j_{\ell}-j_{\ell-1}\leq n_{2}/r$ and $B_{n_{1},j_{\ell}}-B_{n_{1},j_{\ell-1}+1}\leq 1/r$ for $\ell\in[2r]$.

For $\{i_{k}\}_{k=1}^{r+1}$, it suffices to choose $i_{0}=0$, $i_{k}=i_{k-1}+\lfloor n_{1}/r\rfloor$ for $k\in[r]$ and $i_{r+1}=n_{1}$. For $\{j_{\ell}\}_{\ell=1}^{2r}$, since $B_{n_{1},1}=0$, $B_{n_{1},n_{2}}=1$ and $B$ has nondecreasing rows, there is an increasing sequence of indices $\{j_{\ell}^{\prime}\}_{\ell=1}^{r}$ such that $B_{n_{1},j_{\ell}^{\prime}}-B_{n_{1},j_{\ell-1}^{\prime}+1}\leq 1/r$ for all $j\in[r]$. Moreover, by inserting (at most) another $r$ indices between the indices $j_{\ell}^{\prime}$ to obtain a new sequence $\{j_{\ell}\}_{\ell=1}^{2r}$, we can guarantee that not only $B_{n_{1},j_{\ell}}-B_{n_{1},j_{\ell-1}+1}\leq 1/r$, but also $j_{\ell}-j_{\ell-1}\leq n_{2}/r$.

5.14.2 Low-rank approximation.

Let $\{i_{k}\}_{k=1}^{r+1}$ and $\{j_{\ell}\}_{\ell=1}^{2r}$ be chosen so that the above conditions are satisfied. We define a matrix $X\in\mathbbm{R}^{n_{1}\times n_{2}}$ by setting $X_{i,j}=B_{i,j_{\ell-1}+1}$ for all $i\in[n_{1}]$ and $j_{\ell-1}<j\leq j_{\ell}$ where $\ell\in[2r]$. By definition, all columns of $X$ with indices in $(j_{\ell-1},j_{\ell}]$ are the same, so $X$ has rank at most $2r$.

Furthermore, we define a matrix $Y\in\mathbbm{R}^{n_{1}\times n_{2}}$ by setting $Y_{i,j}=(B-X)_{i_{k-1}+1,j}$ for all $j\in[n_{2}]$ and $i_{k-1}<i\leq i_{k}$ where $k\in[r+1]$. Similarly, all rows of $Y$ with indices in $(i_{k-1},i_{k}]$ are the same, so $Y$ has rank at most $r+1$.

It remains to bound $\|\Delta\|_{F}^{2}$ where $\Delta=X+Y-B$. First, let us focus on a block with double indices in $(i_{k-1},i_{k}]\times(j_{\ell-1},j_{\ell}]$. On each of these blocks, by definition it holds that:

• •

$X$ has constant rows, and $Y$ has constant columns;

• •

the first column of $Y$ is zero, and the first column and first row of $\Delta$ is zero.

Then by Lemma 7 (applied with the corresponding blocks of $(B,X,Y,\Delta)$ in place of $(\theta,R,S,B)$), we see that $\Delta$ is bivariate isotonic on each of these block, and

[TABLE]

Therefore, it follows that

[TABLE]

where the second inequality holds thanks to the above choice of indices.

Summing over all the blocks, we obtain

[TABLE]

5.14.3 Bounding the sum of variations.

It remains to bound the above sum. Recall that a telescoping sum gives

[TABLE]

where $(i)$ holds because the first row and first column of $B$ are zero, and $(ii)$ holds because of our choice of $\{j_{\ell}\}_{\ell=1}^{2r}$. As a result, it holds by Hölder’s inequality that

[TABLE]

We therefore obtain $\|\Delta\|_{F}^{2}\leq\frac{2n_{1}n_{2}}{r^{3}}.$ The proof is complete since $\Delta=X+Y-B=R+S+X+Y-\theta$, where the matrix $R+S+X+Y$ has rank at most $3r+3$.

Acknowledgment

We thank Adityanand Guntuboyina for discussing their concurrent work with us. PR was supported by NSF awards IIS-1838071, DMS-1712596 and DMS-TRIPODS-1740751; ONR grant N00014-17- 1-2147 and grant 2018-182642 from the Chan Zuckerberg Initiative DAF. ER was supported by an NSF MSPRF DMS-1703821.

Appendix A Existing results

In this appendix, we state some existing results that are used in our proofs.

The following result is Talagrand’s majorizing measure theorem [66]. See also Theorem 8.6.1 of [68].

Theorem 19** (Talagrand’s majorizing measure theorem).**

For any set $\mathcal{M}\subset\mathbbm{R}^{d}$ and a Gaussian random vector $\varepsilon\sim\mathcal{N}(0,I_{d})$, we have

[TABLE]

where $\gamma_{2}(\cdot)$ denotes Talagrand’s $\gamma_{2}$ functional.

The following theorem gives a tail bound on the supremum of a sub-Gaussian process [66, 68].

Theorem 20** (Generic chaining tail bound).**

Consider a sub-Gaussian vector $\varepsilon\sim\mathsf{subG}_{d}(\sigma^{2})$. For any set $\mathcal{M}\subset\mathbbm{R}^{d}$ and $s>0$, it holds with probability at least $1-2\exp(-s^{2})$ that

[TABLE]

Proof.

By Theorem 8.5.5 of [68], we have that for any $\theta^{*}\in\mathcal{M}$,

[TABLE]

with probability at least $1-2\exp(-s^{2})$, where $\mathsf{diam}(\cdot)$ denotes the diameter of a set. Moreover, we have with probability $1-2\exp(-s^{2})$ that

[TABLE]

It follows that with probability at least $1-4\exp(-s^{2})$,

[TABLE]

Since $\mathsf{diam}(\mathcal{M})+\|\theta^{*}\|_{2}\lesssim\sup_{\theta\in\mathcal{M}}\|\theta\|_{2}$, the proof is complete. ∎

Assouad’s lemma is used to prove the lower bounds (see [67, Lemma 24.3]).

Theorem 21** (Assouad’s Lemma).**

Consider a parameter space $\mathcal{M}$. Let $\mathbbm{P}_{\theta}$ denote the distribution of the observation given that the true parameter is $\theta\in\mathcal{M}$. Let $\mathbbm{E}_{\theta}$ denote the corresponding expectation. Suppose that for each $\tau\in\{-1,1\}^{d}$, there is an associated $\theta^{\tau}\in\mathcal{M}$. Then it holds that

[TABLE]

where $\ell$ denotes any distance function on $\mathcal{M}$, $d_{H}$ denotes the Hamming distance, $\|\cdot\|_{TV}$ denotes the total variation distance, and the infimum is taken over all estimators $\tilde{\theta}$ measurable with respect to the observation $y$.

The following lemma can be proven with basic concentration theory, or follows as a special instance from the Hanson-Wright inequality [37, 59].

Lemma 22**.**

Suppose that $\varepsilon\in\mathbbm{R}^{n}$ is a random vector with independent centered $\mathsf{subG}(\sigma^{2})$ entries. Then it holds that for all $t\geq 0$,

[TABLE]

Appendix B Lower bounds for the SVT estimator

Let us focus on the special case $n_{1}=n_{2}=n$, and study a worst-case matrix in $\mathcal{M}=\mathcal{M}^{n,n}$, for which the approximation rate given by Proposition 6 is tight.

Lemma 23**.**

There exists a matrix $\theta\in\mathcal{M}^{n,n}$ such that

[TABLE]

Proof.

Consider the matrix

[TABLE]

for $V_{0}>0$. It is anti-Monge because

[TABLE]

It also follows from (2.4) that

[TABLE]

Moreover, the singular values of $D$ are given in (5.67), so the eigenvalues of $\theta$ are

[TABLE]

Using the fact that $x/2\leq\sin x\leq x$ for $x\in[0,\pi/2]$, we obtain

[TABLE]

If the spectral decomposition of $\theta$ is $\theta=\sum_{i=1}^{n}\mu_{i}w_{i}w_{i}^{\top}$, then the best rank-$r$ approximation of $\theta$ in the Frobenius norm is $\theta_{r}=\sum_{i=1}^{r}\mu_{i}w_{i}w_{i}^{\top}$, and

[TABLE]

which completes the proof. ∎

Since the anti-Monge matrix $\theta$ in the above proof cannot be approximated by a low-rank matrix at a better rate, we conjecture that for this choice of $\theta$, the rate of convergence given by Theorem 5 for the SVT estimator (3.19) is tight. Intuitively, if we set threshold $\rho$ in definition (3.19) to be larger, then the resulting estimator has a lower rank, thus incurring a larger bias according to the above lemma. On the other hand, setting threshold $\rho$ to be smaller incurs a larger variance due to the noise $\varepsilon$.

More precisely, under the model $y=\theta+\varepsilon$ where $\varepsilon$ has i.i.d. $\mathcal{N}(0,\sigma^{2})$ entries, we conjecture that for any choice of threshold $\rho$ in the estimator (3.19), it holds with constant probability that

[TABLE]

This is because we believe that the bias-variance trade-off in the proof of Theorem 5 is optimal. However, we are unable to prove a lower bound based on a similar argument in [62], since there, the authors are able to exploit a varying signal-to-noise ratio within the classes of matrices they consider. This allows them to employ a triangle inequality argument instead of a explicit bias-variance decomposition that holds with equality. Potential other approaches include analyzing an explicit unbiased estimate of the risk for SVT [12], and studying the exact asymptotic optimal choice of threshold $\rho$ as in [36]. However, since any of these approaches require asymptotic random matrix theory, we consider them beyond the scope of the current work.