# Optimal Estimator for Linear Regression
with Shuffled Labels

Hang Zhang, Ping Li
Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004, USA
{zhanghanghitomi<EMAIL_ADDRESS>
###### Abstract

111Preliminary results appeared in Proceedings of the 37th International
Conference on Machine Learning (ICML’20).

This paper considers the task of linear regression with shuffled labels, i.e.,
$\mathbf{Y}=\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural}+\mathbf{W}$,
where $\mathbf{Y}\in\mathbb{R}^{n\times
m},\bm{\Pi}^{\natural}\in\mathbb{R}^{n\times
n},\mathbf{X}\in\mathbb{R}^{n\times
p},\mathbf{B}^{\natural}\in\mathbb{R}^{p\times m}$, and
$\mathbf{W}\in\mathbb{R}^{n\times m}$, respectively, represent the sensing
results, (unknown or missing) corresponding information, sensing matrix,
signal of interest, and additive sensing noise. Given the observation
$\mathbf{Y}$ and sensing matrix $\mathbf{X}$, we propose a one-step estimator
to reconstruct $(\bm{\Pi}^{\natural},\mathbf{B}^{\natural})$. From the
computational perspective, our estimator’s complexity is $O(n^{3}+np^{2}m)$,
which is no greater than the maximum complexity of a linear assignment
algorithm (e.g., $O(n^{3})$) and a least square algorithm (e.g.,
$O(np^{2}m)$). From the statistical perspective, we divide the minimum
$\mathsf{SNR}$ requirement into four regimes, e.g., unknown, hard, medium, and
easy regimes; and present sufficient conditions for the correct permutation
recovery under each regime: $(i)$ $\mathsf{SNR}\geq\Omega(1)$ in the easy
regime; $(ii)$ $\mathsf{SNR}\geq\Omega(\log n)$ in the medium regime; and
$(iii)$ $\mathsf{SNR}\geq\Omega((\log n)^{c_{0}}\cdot
n^{\nicefrac{{c_{1}}}{{\textup{srank}(\mathbf{B}^{\natural})}}})$ in the hard
regime ($c_{0},c_{1}$ are some positive constants and
$\textup{srank}(\mathbf{B}^{\natural})$ denotes the stable rank of
$\mathbf{B}^{\natural}$). In the end, we also provide numerical experiments to
confirm the above claims.

##  1 Introduction

The recent years have witnessed a renaissance of permuted linear regression,
or “unlabeled sensing", due to its broad spectrum of applications ranging from
database merging, to privacy, to communications, to computer vision, to
robotics, to sensor networks, etc (Pananjady et al., 2018; Unnikrishnan et
al., 2015; Slawski et al., 2020; Slawski and Ben-David, 2019; Pananjady et
al., 2017; Zhang et al., 2022). Here, we briefly review several use cases of
permuted linear regression, including record linkage, data de-anonymization,
and header-free communications. For more information about other applications,
we refer interested readers to the references thereof.

  * •

Record linkage. Given multiple databases containing information about the same
entities, the objective is to merge them into one comprehensive database.
However, these data may not be well-aligned due to the data formatting or data
quality issues. Modeling the mismatches as a permutation is investigated as a
mitigation strategy.

  * •

Data de-anonymization. This task can be regarded as the opposite side of
privacy protection. The intruders aim to infer the hidden identities/labels in
certain private networks with public information. One commonly used method is
to compare the correlation between this information and to find matching pairs
with maximum correlation sum. Here, the permuted linear regression arises as a
natural generalization.

  * •

Header-free communication. Another potential application of permuted linear
regression is in the _internet of things_ (IOT) network, where communications
happen frequently while the transmitted messages are usually in short length.
To improve the bandwidth efficiency, the sensor identity is omitted during
transmission and thus the signal decoding involves first restoring the
correspondence.

In this paper, we consider the canonical setting for a permuted linear
sensing, which is written as

$\displaystyle\mathbf{Y}=\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural}+\mathbf{W},$
(1)

where $\bm{\Pi}^{\natural}\in\\{0,1\\}^{n\times n}$ denotes the unknown
permutation matrix, $\mathbf{X}\in\mathbb{R}^{n\times p}$ represents the
design (sensing) matrix, $\mathbf{B}^{\natural}\in\mathbb{R}^{p\times m}$ is
the signal of interest, and $\mathbf{W}\in\mathbb{R}^{n\times m}$ denotes the
additive noise. Compared with the canonical model of linear regression with
well-aligned data, our task is to infer both the signal
$\mathbf{B}^{\natural}$ and the missing correspondence information
$\bm{\Pi}^{\natural}$ from the pair $(\mathbf{X},\mathbf{Y})$.

###  1.1 Related work

The study on permuted linear regression has a long history and can be at least
traced back to 1970s under the name of “broken sample problem” (DeGroot and
Goel, 1976, 1980; Goel, 1975; Bai and Hsing, 2005). In recent years, we have
witnessed a revival of its study and can broadly divide these research works
into two categories: $(i)$ single observation model and $(ii)$ multiple
observations model.

We first discuss the work on the single observation model. In Unnikrishnan et
al. (2015), they focus on the noiseless single observation model, namely,
$\mathbf{W}=\bm{0}$ and $m=1$. Assuming the entries in the sensing matrix
$\mathbf{X}$ are drawn from a continuous distribution and
$\mathbf{B}^{\natural}\in\mathbb{R}^{p}$ is an arbitrary vector residing
within a linear space with dimension $p$, Unnikrishnan et al. (2015) establish
the necessary condition $n\geq 2p$ for the correct recovery. Similar results
have been discovered by Dokmanic (2019) but with different approaches. In
Pananjady et al. (2018), they investigate the noisy case under the single
observation model, i.e., $m=1$. After obtaining the statistical limit of the
minimum $\mathsf{SNR}$ required for permutation recovery, e.g.,
$\Omega(n^{c})$ ($c>0$ is some positive constant), they analyze the _maximum
likelihood_ (ML) estimator and show its performance matches the order of the
statistical limits thereof. However, the ML estimator is NP-hard in general
except for the special case when $p=1$, in other words,
$\mathbf{B}^{\natural}$ is a scalar. Later, an approximation algorithm for
permutation recovery is presented in Hsu et al. (2017); and in Slawski and
Ben-David (2019), the permutation recovery is studied from the viewpoint of
denoising. We will give a more technical discussion of these works (Pananjady
et al., 2018; Hsu et al., 2017; Slawski and Ben-David, 2019), as their
settings have a huge overlap with ours. The content is deferred to Section 3.3
until we have collected the required facts. An independent line of research
can be found in Tsakiris and Peng (2019); Peng et al. (2021), where the
permuted linear regression is studied from the viewpoint of algebraic
geometry.

Next, we survey the work on the multiple observations model (Pananjady et al.,
2017; Slawski et al., 2020; Zhang et al., 2022; Zhang and Li, 2023a, b). In
Pananjady et al. (2017), the authors study the problem of reconstructing the
product $\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural}$. In Zhang et al.
(2022), the focus is shifted to reconstruct the individual values of
$\bm{\Pi}^{\natural}$ and $\mathbf{B}^{\natural}$. Statistical limits are
presented as well as an investigation of the ML estimator. In Slawski et al.
(2020), they take a similar viewpoint of Slawski and Ben-David (2019) and use
a denoising-based method for correspondence recovery. In Zhang and Li (2023b),
a follow-up of this work, they consider the sparse matrix setting (i.e., each
column of $\mathbf{B}^{\natural}$ is $k$-sparse) and design the estimator
based on Zhang and Li (2020). A more detailed comparison between our works and
the existing literature can be found in Section 4.3.

In addition to the works mentioned above, other research, e.g., Haghighatshoar
and Caire (2018); Tang et al. (2021); Jeong et al. (2020); Slawski and Sen
(2022); Fang and Li (2023), are also related to our works. A detailed
discussion is omitted due to their relatively loose connection.

###  1.2 Contributions

We define the _signal-to-noise-ratio_ ($\mathsf{SNR}$) as
${|\kern-1.72218pt|\kern-1.72218pt|\mathbf{B}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}^{2}/(m\cdot\sigma^{2})$
before describing our contributions.

  * •

We propose a one-step estimator for permuted linear regression, which consists
of two sub-parts: one for the permutation recovery and the other for the
signal recovery. We show the first sub-part is with computational complexity
$O(n^{3})$, which is the same as the oracle estimator (i.e., a linear
assignment algorithm) where $\mathbf{B}^{\natural}$ is given a prior. For the
second sub-part, our estimator’s computational cost is $O(np^{2}m)$, which is
the same as the least square algorithm generally associated with classical
linear regression.

  * •

We prove that our estimator almost reaches the statistical limits for
permutation recovery. First, we study the single observation model (i.e.,
$m=1$) and show our estimator can yield the ground truth permutation matrix
when $\log\mathsf{SNR}\gtrsim\log n$ and signal length $p$ is one, in other
words, $\mathbf{B}^{\natural}$ is a scalar. This bound matches the statistical
limit up to some positive constant. Second, we investigate the multiple
observations model, namely, $m>1$. We divide the required $\mathsf{SNR}$ into
four regimes, i.e., unknown, hard, medium, and easy regimes. In the medium and
hard regimes, we show that our $\mathsf{SNR}$ requirement for correct
permutation recovery matches the statistical limits up to a multiplicative
polynomial of $\log n$; while in the easy regime, our $\mathsf{SNR}$
requirement matches the limits up to a multiplicative factor of certain
positive constant. Compared with previous works, our work has a lead in both
the computational and statistical perspectives. A detailed summary is put in
Table 1.

###  1.3 Notations

Denote $c$, $c^{{}^{\prime}}$, $c_{i}$ as some positive constants, whose
values are not necessarily the same even for those with the same notations. We
denote $a\lesssim b$ if there exists some positive constants $c_{0}>0$ such
that $a\leq c_{0}b$. Similarly, we define $a\gtrsim b$ provided $a\geq c_{0}b$
for some positive constant $c_{0}$. We write $a\asymp b$ when $a\lesssim b$
and $a\gtrsim b$ hold simultaneously.

We call $z$ a sub-gaussian _random variable_ (RV) with $\|z\|_{\psi_{2}}\leq
K$ if it satisfies $\mathbb{E}e^{\nicefrac{{z^{2}}}{{K^{2}}}}\leq 2$ (Section
$2.5.2$ in Vershynin (2018)). A centered isotropic random vector $\bm{x}$ is
defined such that $\mathbb{E}\bm{x}=\bm{0}$ and
$\mathbb{E}\bm{x}\bm{x}^{\top}=\mathbf{I}$ (Definition $3.2.1$ in Vershynin
(2018)).

For an arbitrary matrix $\mathbf{M}$, we denote $\mathbf{M}_{i,:}$ as its
$i$th row, $\mathbf{M}_{:,i}$ as its $i$th column, and $\mathbf{M}_{ij}$ as
its $(i,j)$th element. The Frobenius norm of $\mathbf{M}$ is defined as
${|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}$
while the operator norm is denoted as
${|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}$,
whose definitions can be found in Section $2.3$ of Golub and Van Loan (2012).
Its stable rank is defined as
$\textup{srank}(\cdot)\triangleq{|\kern-1.72218pt|\kern-1.72218pt|\cdot|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}^{2}/{|\kern-1.72218pt|\kern-1.72218pt|\cdot|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}^{2}$
(Section $2.1.15$ in Tropp (2015)). The inner product between matrices are
denoted as
$\left\langle\kern-2.75554pt\left\langle\cdot,\cdot\right\rangle\kern-2.75554pt\right\rangle$;
while the inner product between vectors are denoted as
$\left\langle\cdot,\cdot\right\rangle$.

Associate with each permutation matrix $\bm{\Pi}$, we define the operator
$\pi(\cdot)$ that transforms index $i$ to $\pi(i)$ under $\bm{\Pi}$. The
Hamming distance $\mathsf{d_{H}}(\bm{\Pi}_{1},\bm{\Pi}_{2})$ between
permutation matrix $\bm{\Pi}_{1}$ and $\bm{\Pi}_{2}$ is defined as
$\mathsf{d_{H}}\left(\bm{\Pi}_{1},\bm{\Pi}_{2}\right)=\sum_{i=1}^{n}\mathbbm{1}\left(\pi_{1}(i)\neq\pi_{2}(i)\right)$.
The _signal-to-noise-ratio_ ($\mathsf{SNR}$) is defined as
$\mathsf{SNR}={|\kern-1.72218pt|\kern-1.72218pt|\mathbf{B}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}^{2}/(m\cdot\sigma^{2})$.
More notations are in the supplementary material.

Table 1: Comparison with the prior art. All results are presented in their
best orders, which may only hold true in certain regimes. The computational
cost refers to the number of iterations. Moreover, the logarithmic term is
omitted in the notation $\widetilde{\Omega}(\cdot)$ and
$\widetilde{O}(\cdot)$. Notation $\textup{rank}(\cdot)$ denotes the rank of
the corresponding matrix, $n_{\textup{min}}$ denotes the minimum sample
number, and $h_{\textup{max}}$ denotes the maximum allowed number of permuted
rows. Notation $\checkmark$ means the requirement is met; $\bigtimes$ means
the requirement is not met; and N/A means not applied.

| Comput. Optim. |  | Statis. Optim. |  | ${n_{\textup{min}}}/{p}~{}(\geq)$ |  | ${h_{\textup{max}}}/{n}~{}(\leq)$
---|---|---|---|---|---|---|---
| $m=1$ | $m\gg 1$ |  | $m=1$ | $m\gg 1$ |  | $m=1$ | $m\gg 1$ |  | $m=1$ | $m\gg 1$
(Pananjady et al., 2018) | $\checkmark$ | N/A |  | $\checkmark$ | N/A |  | $\widetilde{\Omega}(1)$ | N/A |  | $\widetilde{O}(1)$ | N/A
(Hsu et al., 2017) | $\bigtimes$ | N/A |  | $\bigtimes$ | N/A |  | $\widetilde{\Omega}(1)$ | N/A |  | $\widetilde{O}(1)$ | N/A
(Slawski and Ben-David, 2019) | $\bigtimes$ | N/A |  | $\checkmark$ | N/A |  | $\widetilde{\Omega}(1)$ | N/A |  | $\widetilde{O}(\log^{-1}n)$ | N/A
(Zhang et al., 2022) | N/A | $\bigtimes$ |  | N/A | $\checkmark$ |  | N/A | $\widetilde{\Omega}(1)$ |  | N/A | $\widetilde{O}\left(\log^{-1}\textup{rank}(\mathbf{B}^{\natural})\right)$
(Slawski et al., 2020) | N/A | $\bigtimes$ |  | N/A | $\checkmark$ |  | N/A | $\widetilde{\Omega}(p)$ |  | N/A | $\widetilde{O}(\log^{-1}n)$
This work | $\checkmark$ | $\checkmark$ |  | $\checkmark$ | $\checkmark$ |  | $\widetilde{\Omega}(1)$ | $\widetilde{\Omega}(1)$ |  | $\widetilde{O}(1)$ | $\widetilde{O}(1)$

###  1.4 Road map

The organization of this paper is as follows. In Section 2, we formally state
our problem setting, present our estimator and its design insight, and review
the minimax lower bounds. Then, we separately investigate our estimator’s
statistical properties under the single observation model ($m=1$) and multiple
observations model ($m>1$). Corresponding discussions are put in Section 3 and
Section 4, respectively. Simulation results are presented in Section 5 and the
conclusions are drawn in Section 6. The technical details are deferred to the
Appendix.

##  2 Problem Setting

We start the discussion with a formal restatement of the sensing model

$\displaystyle\mathbf{Y}=\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural}+\mathbf{W},\vspace{-0.05in}$
(2)

where $\mathbf{Y}\in\mathbb{R}^{n\times m}$ denotes the observation,
$\bm{\Pi}^{\natural}\in\\{0,1\\}^{n\times n}$ is the unknown permutation
matrix such that
$\sum_{i}\bm{\Pi}^{\natural}_{i,j}=\sum_{j}\bm{\Pi}^{\natural}_{i,j}=1$,
$\mathbf{X}\in\mathbb{R}^{n\times p}$ denotes the sensing matrix such that
each entry $\mathbf{X}_{ij}$ are i.i.d. centered isotropic sub-gaussian RV
with ${\left\|\mathbf{X}_{ij}\right\|}_{\psi_{2}}\lesssim 1$, i.e.,
$\mathbb{E}\mathbf{X}_{ij}=0$ and $\mathbb{E}\mathbf{X}_{ij}^{2}=1$,
$\mathbf{B}^{\natural}\in\mathbb{R}^{p\times m}$ denotes the signal of
interests, and $\mathbf{W}\in\mathbb{R}^{n\times m}$ represents the additive
Gaussian noise with each entry $\mathbf{W}_{ij}$ being Gaussian RV with zero
mean and $\sigma^{2}$ variance, namely,
$\mathbf{W}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\mathsf{N}(0,\sigma^{2})$.
222We call $n$ sample number, $p$ length of signal, and $m$ measurement
number.

Our goal is to reconstruct the pair
$(\bm{\Pi}^{\natural},\mathbf{B}^{\natural})$ from observations $\mathbf{Y}$
and the sensing matrix $\mathbf{X}$. In the following context, we put our
major focus on the permutation recovery: on one hand, this problem reduces to
the classical setting of linear regression once ground truth permutation is
obtained; on the other hand, no meaningful bound on
${|\kern-1.72218pt|\kern-1.72218pt|\widehat{\mathbf{B}}-\mathbf{B}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}$
can be obtained with incorrect correspondence information (e.g.,
$\bm{\Pi}^{\natural}$).

###  2.1 Estimator and its design insight

Algorithm 1 One-Step Estimator.

Input: observation $\mathbf{Y}$ and sensing matrix $\mathbf{X}$.

Output: pair $(\bm{\Pi}^{\textup{opt}},~{}\mathbf{B}^{\textup{opt}})$, which
is written as  $\displaystyle\bm{\Pi}^{\textup{opt}}$
$\displaystyle=\textup{argmax}_{\bm{\Pi}\in{\cal P}_{n}}\
\left\langle\kern-2.75554pt\left\langle\bm{\Pi},\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}\right\rangle\kern-2.75554pt\right\rangle,$
(3) $\displaystyle\mathbf{B}^{\textup{opt}}$
$\displaystyle=\mathbf{X}^{\dagger}\widehat{\bm{\Pi}}^{\top}\mathbf{Y},$ (4)
where
$\mathbf{X}^{\dagger}=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}$
denotes the pseudo-inverse of $\mathbf{X}$ (Section $5.5.2$ in Golub and Van
Loan (2012)) and ${\cal P}_{n}$ is the set of all possible permutation
matrices.

We propose a one-step estimator, whose details are summarized in Algorithm 1.
Before a thorough investigation of our estimator’s properties, we first
present its underlying design insight, which is quite straightforward.
Considering the oracle situation where $\mathbf{B}^{\natural}$ is given a
prior, we can reconstruct the permutation matrix $\bm{\Pi}^{\natural}$ via

$\displaystyle\bm{\Pi}^{\textup{opt}}=\textup{argmax}_{\bm{\Pi}\in{\cal
P}_{n}}\left\langle\kern-2.75554pt\left\langle\bm{\Pi},\mathbf{Y}(\nicefrac{{\mathbf{B}^{\natural\top}}}{{\alpha}})\mathbf{X}^{\top}\right\rangle\kern-2.75554pt\right\rangle,$
(5)

where $\alpha>0$ is an arbitrary scaling constant. Back to our case, we can
see that the major difficulty comes from the lack of knowledge about
$\mathbf{B}^{\natural}$, to put more precisely, the direction of
$\mathbf{B}^{\natural}$, since the solution to (5) remains the same up to some
positive scaling factor.

Inspired by the recent progress in non-convex optimization Candès and Plan
(2010); Chi et al. (2019); Balakrishnan et al. (2017), we would like to
approximate $\mathbf{B}^{\natural}$’s direction with $\mathbf{X}$ and
$\mathbf{Y}$. Note that
$\mathbb{E}\mathbf{X}^{\top}\mathbf{Y}=(n-h)\mathbf{B}^{\natural}$, which is
parallel to $\mathbf{B}^{\natural}$ if $h<n$ ($h$ denotes the number of
permuted rows by $\bm{\Pi}^{\natural}$, e.g.,
$h\triangleq\mathsf{d_{H}}(\mathbf{I},\bm{\Pi}^{\natural})$). Assuming that
$\mathbf{X}^{\top}\mathbf{Y}$ is close to
$\mathbb{E}\mathbf{X}^{\top}\mathbf{Y}$, we design our estimator by
substituting $\mathbf{B}^{\natural}$ in (5) with
$\mathbf{X}^{\top}\mathbf{Y}$, which then leads to the permutation estimation
in (3). Once the permutation is obtained, we can restore (2) to the classical
setting of linear regression and then estimate $\mathbf{B}^{\natural}$ with a
_least-square_ estimator (Golub and Van Loan, 2012).

Although at first glance this idea looks simple, if not naive, in the
following context, we will show that our estimator can reach statistical
optimality in a broad regime. In addition, our estimator is tuning-free, to
put more specifically, no estimation of the noise variance $\sigma^{2}$ nor
the number of permuted rows $h$ is required.

###  2.2 Computational cost

This subsection concerns our estimator’s computational cost. We consider two
types of oracle estimators as benchmarks.

  * •

Oracle estimator I. We consider the oracle scenario where
$\mathbf{B}^{\natural}$ is given a prior. Then, we can recover
$\bm{\Pi}^{\textup{opt}}$ as in (5), which is with computational cost
$O(n^{3})$.

  * •

Oracle estimator II. We consider the oracle scenario when
$\bm{\Pi}^{\natural}$ is known in advance. The sensing relation in (2) reduces
to the classical multivariate linear regression and least-square algorithm for
$\mathbf{B}^{\textup{opt}}$ takes up to $O(np^{2}m)$ time.

Then we study our estimator’s computational cost. In the first step (3), we
only sacrifice one matrix multiplication, i.e., replacing
$\mathbf{B}^{\natural}$ by the product $\mathbf{X}^{\top}\mathbf{Y}$. Since
the computational bottleneck lies in solving the linear assignment problem
(Kuhn, 1955; Bertsekas and Castañón, 1992), one additional matrix
multiplication does not change the computational complexity, which is also of
order $O(n^{3})$. Similarly, in the second step (4) our estimator sacrifices
another matrix multiplication, whose cost is negligible when compared with the
total cost of Oracle Estimator II.

With the relation $n^{3}+np^{2}m\leq 2(n^{3}\vee np^{2}m)$, we conclude our
estimator’s computational cost $O(n^{3}+np^{2}m)$ is in the same order as the
maximum computational costs of the above-mentioned two oracle estimators.

###  2.3 Mini-max lower bounds

Before studying the statistical properties of our algorithm, we review the
mini-max lower bounds on the permutation recovery. First, we consider the
single observation model, i.e., $m=1$. We have

######  Theorem 1 (Theorem $2$ in Pananjady et al. (2018)).

For any estimator $\widehat{\bm{\Pi}}$, we have the error probability
$\mathbb{P}(\widehat{\bm{\Pi}}\neq\bm{\Pi}^{\natural})$ exceed $1-c_{0}\cdot
e^{-c_{1}n\delta}$ provided that $2+\log(1+\mathsf{SNR})\leq(2-\delta)\log
n~{},0<\delta<2$.

This theorem suggests that we need $\mathsf{SNR}$ to be at least of order
$\Omega(n^{c})$ to avoid construction failure of $\bm{\Pi}^{\natural}$. Then
we move on to the multiple observations model, i.e., $m>1$. The corresponding
lower bound is summarized as

######  Theorem 2 (Theorem $1$ in Zhang et al. (2022)).

For any estimator $\widehat{\bm{\Pi}}$, we have the error probability
$\mathbb{P}(\widehat{\bm{\Pi}}\neq\bm{\Pi}^{\natural})$ exceed $1/2$, provided
that
$\log\det\left(\mathbf{I}+\nicefrac{{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}}{{\sigma^{2}}}\right)<\frac{\log
n!-2}{n}$.

######  Remark 1.

Theorem 1 can be regarded as a special case of Theorem 2, since the former
theorem can be obtained from the latter theorem by letting $m=1$.

Denote $\lambda_{i}(\cdot)$ as the $i$th singular value, we can rewrite
$\log\det\left(\mathbf{I}+\nicefrac{{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}}{{\sigma^{2}}}\right)$
as
$\sum_{i}\log[1+\nicefrac{{\lambda_{i}^{2}(\mathbf{B}^{\natural})}}{{\sigma^{2}}}]$
and approximately
$\log\det\left(\mathbf{I}+\nicefrac{{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}}{{\sigma^{2}}}\right)$
as $\textup{srank}(\mathbf{B}^{\natural})\cdot\log(1+\mathsf{SNR})$. Then, we
conclude that no estimator can reliably recover $\bm{\Pi}^{\natural}$ if
$\log\mathsf{SNR}\lesssim\frac{\log
n}{\textup{srank}(\mathbf{B}^{\natural})}$.

The following context studies our algorithm’s statistical properties under the
single observation model ($m=1$) and multiple observations model ($m>1$),
respectively. We will show that the required $\mathsf{SNR}$s for correct
permutation recovery are close to the statistical limits thereof.

##  3 Single Observation Model

This section considers the single observation model, namely, $m=1$. To
distinguish this case with multiple observations model, i.e., $m>1$, we
rewrite the sensing relation in (2) as

$\displaystyle\bm{y}=\bm{\Pi}^{\natural}\mathbf{X}\bm{\beta}^{\natural}+\bm{w},$
(6)

where $\bm{y},~{}\bm{w}\in\mathbb{R}^{n}$ and
$\bm{\beta}^{\natural}\in\mathbb{R}^{p}$ all reduce to vectors. Based on
whether $p\geq 2$, we find that our estimator exhibits vastly different
behaviors.

###  3.1 A warm-up example: $p=1$

First, we consider a warm-up example where $m=1$ and $p=1$. In such situation,
we have $\mathbf{X}$ to be a vector of length $n$ and $\bm{\beta}^{\natural}$
be a scalar. Then we have

######  Theorem 3.

Assume $\bm{x}$ be an isotropic log-concave random vector with zero mean and
${\left\|\bm{x}\right\|}_{\psi_{2}}\lesssim 1$. Consider the large-system
limit where $n$ is sufficiently large and assume $(i)$
$h=\mathsf{d_{H}}(\mathbf{I},\bm{\Pi}^{\natural})\leq\nicefrac{{n}}{{4}}$, and
$(ii)$ $n\geq 2p$. Our estimator in (3) can return the correct permutation
matrix, i.e., $\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural}$, with probability
at least $1-c_{1}n^{-1}$ provided that $\log\mathsf{SNR}\geq c\cdot\log n$.

###### Proof outline.

To emphasize the fact that sensing matrix $\mathbf{X}$ reduces to a vector
when $p=1$, we denote it as $\bm{x}$. To start with, we consider the noiseless
case, where $\mathsf{SNR}$ is infinite. With some simple algebraic
manipulations, we can expand the inner product
$\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle$ in (3) as

$\displaystyle\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle=(\beta^{\natural})^{2}\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle\cdot\langle\bm{\Pi}\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle.$
(7)

Since the first term
$(\beta^{\natural})^{2}\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle$ is
positive with high probability, we conclude that the maximum is reached when
$\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural}$. Back to the noisy case, we
interpret the observation $\bm{y}$ as a perturbed version of product
$\bm{\Pi}^{\natural}\bm{x}\beta^{\natural}$. Provided the perturbation is
significantly small, our estimator can obtain the correct permutation matrix
$\bm{\Pi}^{\natural}$ with high likelihood. For the technical details, we
refer the interested readers to the appendix. ∎

First, we comment on the tightness of Theorem 3. According to Theorem $2$ in
Pananjady et al. (2018), which is restated as Theorem 1, correct permutation
recovery requires $\log\mathsf{SNR}\gtrsim\log n$ at least. Easily, we can see
that our estimator matches this statistical limit with a difference up to some
multiplicative constant and hence concludes its tightness.

In addition, we would like to mention the general case where $p\geq 2$ can be
transformed to the warm-up case where $p=1$ if the direction
$\bm{e}\triangleq\bm{\beta}^{\natural}/{|\kern-1.72218pt|\kern-1.72218pt|\bm{\beta}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}$
is known. The detailed construction method comes as follows. With the Gram-
Schmidt process (Section $5.2.7$ in Golub and Van Loan (2012)), we can
construct an orthonormal matrix $\mathbf{Q}\in\mathbb{R}^{p\times p}$ whose
first column is $\bm{e}$. Then we rewrite (6) as

$\displaystyle\bm{y}={|\kern-1.72218pt|\kern-1.72218pt|\bm{\beta}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}\bm{\Pi}^{\natural}\left(\mathbf{X}\mathbf{Q}\right)_{:,1}+\bm{w}.\vspace{1mm}$

Easily, we can show $\left(\mathbf{X}\mathbf{Q}\right)_{:,1}$ has similar
statistical properties as $\mathbf{X}_{:,1}$ and can restore it to the model
thereof, to which Theorem 3 also applies.

At last, we should stress that our proof of Theorem 3 relies heavily on the
expansion in (7), which only holds true when $p=1$ or
$\bm{e}\triangleq\bm{\beta}^{\natural}/{|\kern-1.72218pt|\kern-1.72218pt|\bm{\beta}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}$
is known. In fact, we will see shortly that our estimator fails to obtain
correct permutation matrix–even with infinite $\mathsf{SNR}$–when $p$ exceeds
over $2$.

###  3.2 The general case: $p\geq 2$

This subsection studies a more general case where $m=1$ and $p\geq 2$. We
consider the noiseless case ($\sigma=0$) and set $\bm{\beta}^{\natural}$ as
$[1000;~{}1000]^{\top}$. The sample number $n$ is picked to be $1000$. We
evaluate the permutation recovery by Hamming distance
$\mathsf{d_{H}}(\cdot,\bm{\Pi}^{\natural})$.

First, we obtain $(\widehat{\bm{\Pi}}^{(0)},\bm{\beta}^{(0)})$ with Algorithm
1. The error $\mathsf{d_{H}}(\widehat{\bm{\Pi}}^{(0)},\bm{\Pi}^{\natural})$ is
approximately $900$, which means almost all correspondence
$\pi^{\natural}(\cdot)$ are incorrectly detected. To refine the permutation
recovery, we carry out the following alternative minimization

$\displaystyle\widehat{\bm{\Pi}}^{(t+1)}$
$\displaystyle=\textup{argmax}_{\bm{\Pi}\in{\cal
P}_{n}}\left\langle\kern-2.75554pt\left\langle\bm{\Pi},\mathbf{Y}\widehat{\mathbf{B}}^{(t)\top}\mathbf{X}^{\top}\right\rangle\kern-2.75554pt\right\rangle;$
$\displaystyle\widehat{\bm{\beta}}^{(t+1)}$
$\displaystyle=\mathbf{X}^{\dagger}\widehat{\bm{\Pi}}^{(t+1)\top}\mathbf{Y},\vspace{1mm}$

where $\widehat{\bm{\Pi}}^{(t)}$ and $\widehat{\bm{\beta}}^{(t)}$ denote the
reconstructed value of $\bm{\Pi}$ and $\bm{\beta}$ in the $t$th iteration,
respectively. Numerical experiment suggests that the correct permutation
matrix remains out of reach even after $100$ iterations. An illustration is
put in Figure 1.

Figure 1: Hamming distance
$\mathsf{d_{H}}(\widehat{\bm{\Pi}}^{(t)},\bm{\Pi}^{\natural})$ when $n=1000$,
$p=2$, $m=1$, $\bm{\beta}^{\natural}=[1000;~{}1000]^{\top}$, and $\sigma=0$.
We only plot the behavior of
$\mathsf{d_{H}}(\widehat{\bm{\Pi}}^{(t)},\bm{\Pi}^{\natural})$ in the first
$100$ iterations since it remains almost the same in the following $400$
iterations.

The underlying reason is the low stable rank
$\textup{srank}(\mathbf{B}^{\natural})$, which is one in the single
observation model. In the next section, we will show that the ground truth
$\mathbf{B}^{\natural}$ can be obtained in an almost effortless way once
$\textup{srank}(\mathbf{B}^{\natural})$ exceeds certain threshold.

###  3.3 Discussion of related work

This subsection discusses the relation between our estimator and prior work
focusing on the single observation model (e.g., $m=1$) (Slawski and Ben-David,
2019; Pananjady et al., 2018; Hsu et al., 2017; Abid et al., 2017). In
Pananjady et al. (2018), the ML estimator is investigated, which is only
computable for the special case $(m,p)=(1,1)$ and NP-hard for the rest cases.
Their estimator gets the same $\mathsf{SNR}$ requirement as ours, namely,
$\log\mathsf{SNR}\gtrsim\log n$.

To handle the computational issue of the ML estimator, Hsu et al. (2017)
propose an approximation algorithm with polynomial complexity. Their
$\mathsf{SNR}$ requirement is $\mathsf{SNR}\geq c\min\left(1,p/\log\log
n\right)$, which has a gap with the mini-max lower bound. In addition, they
focus on the recovery of $\mathbf{B}^{\natural}$ rather than the permutation
matrix $\bm{\Pi}^{\natural}$.

Later, Slawski and Ben-David (2019) study the problem for the viewpoint of
denoising. By putting a sparse constraint on $h$ (e.g., the number of permuted
rows), they view the term
$(\mathbf{I}-\bm{\Pi}^{\natural})\mathbf{X}\bm{\beta}^{\natural}$ as an
additive sparse outlier. After obtaining an estimate of
$\bm{\beta}^{\natural}$, they restore the permutation information with a
linear assignment algorithm. Compared with Algorithm 1, the estimator in
Slawski and Ben-David (2019) has a more stringent requirement on $h$ (i.e.,
$h\ll\nicefrac{{n}}{{\log n}}$) but at the same time enjoys a broader class
for $\bm{\beta}^{\natural}$, to put more specific, they allow $p>1$ under the
single observation model.

A parallel line of work can be found in Abid et al. (2017), where a consistent
estimator is proposed based on the method-of-moments. However, their analysis
is only for the special case $(m,p)=(1,1)$ and focuses on
${|\kern-1.72218pt|\kern-1.72218pt|\widehat{\mathbf{B}}-\mathbf{B}^{\natural}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}$.
Direct comparison with ours can be difficult.

##  4 Multiple Observations Model

The previous section considers the single observation model ($m=1$) and this
section extends the discussion to the multiple observations model ($m>1$). We
$(i)$ discover a much richer behavior inherent in our estimator and $(ii)$
show that correct permutation can be obtained with $\mathsf{SNR}$ being a
positive constant.

Figure 2: Illustration of different regimes of the permutation recovery. In
regime hard, we require $\mathsf{SNR}\gtrsim
n^{\frac{c}{\textup{srank}(\mathbf{B}^{\natural})}}$; in regime medium
difficulty, we require $\mathsf{SNR}\gtrsim\log^{c}n$; and in regime easy, we
require $\mathsf{SNR}\geq c$.

Based on the requirement on $\mathsf{SNR}$ for correct permutation recovery,
we divide the estimator’s performance into four regimes in terms of
$\textup{srank}(\mathbf{B}^{\natural})$, which is illustrated in Figure 2.

  * •

Unknown regime: $\textup{srank}(\mathbf{B}^{\natural})\in[1,c_{0})$. In this
regime, our estimator’s behavior remains a mystery. To the best of our
knowledge, no estimator’s performance has been thoroughly studied in this
regime.

  * •

Hard regime: $\textup{srank}(\mathbf{B}^{\natural})\in[c_{0},c_{1}\log n)$. We
assume $\mathbf{X}_{ij}$ to be i.i.d. sub-gaussian random variables with log-
concavity. For the correct permutation recovery, we require $\mathsf{SNR}$
satisfying $\mathsf{SNR}\gtrsim
n^{\frac{c}{\textup{srank}(\mathbf{B}^{\natural})}}$.

  * •

Medium regime: $\textup{srank}(\mathbf{B}^{\natural})\in[c_{1}\log
n,c_{2}\log^{4}n)$. We only need $\mathbf{X}_{ij}$ to be i.i.d. sub-gaussian
random variables without enforcing log-concavity. In addition, we require
$\mathsf{SNR}$ satisfying $\mathsf{SNR}\gtrsim\log^{c}n$ to obtain the ground
truth permutation matrix.

  * •

Easy regime:
$\textup{srank}(\mathbf{B}^{\natural})\in[c_{2}\log^{4}n,\infty)$. Still, we
only assume $\mathbf{X}_{ij}$ to be i.i.d. sub-gaussian random variables.
Here, we can relax the requirement on $\mathsf{SNR}$ to be above some positive
constant for correct permutation recovery.

The formal statement of results is put in Theorem 4 and Theorem 5.

######  Theorem 4.

Consider the sensing matrix $\mathbf{X}$ with its entries $\mathbf{X}_{i,j}$
being sub-gaussian RV with zero mean and unit variance $(1\leq i\leq n,1\leq
j\leq p)$. Assuming that $(i)$ $n\gg p\cdot\log^{3}n\cdot\log^{2}(n^{2}p^{3})$
and $(ii)$ $h\leq c_{0}\cdot n$, we can reliably obtain the permutation matrix
with Algorithm 1, namely,
$\mathbb{P}(\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural})\geq
1-c_{1}n^{-c_{2}}-c_{3}p^{-c_{4}}$, in the following situations:

  * •

Easy regime $(\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n)$: we have
$\mathsf{SNR}\geq c$;

  * •

Medium regime $(\log n\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log^{4}n)$:
we have $\log\mathsf{SNR}\gtrsim\log\log n$.

Note that the above result only applies to the easy and medium regime, i.e.,
$\textup{srank}(\mathbf{B}^{\natural})\gg\log n$. By enforcing additional
constraints, to put more specifically, $\mathbf{X}_{ij}$ is log-concave, we
can generalize the above result to the hard regime, which is formally stated
as

######  Theorem 5.

Consider the sensing matrix $\mathbf{X}$ with its entries $\mathbf{X}_{i,j}$
being log-concave sub-gaussian RV with zero mean and unit variance $(1\leq
i\leq n,1\leq j\leq p)$. Assuming that $(i)$ $n\gg
p^{1+\varepsilon}\cdot\log^{3(1+\varepsilon)}n\cdot\log^{2(1+\varepsilon)}(n^{2}p^{3})$
and $(ii)$ $h\leq c_{0}\cdot n$, we can reliably obtain the permutation matrix
with Algorithm 1, namely,
$\mathbb{P}(\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural})\geq
1-c_{1}n^{-c_{2}}-c_{3}p^{-c_{4}}$, in the following situations:

  * •

Easy regime $(\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n)$: we have
$\mathsf{SNR}\geq c$;

  * •

Medium regime $(\log n\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log^{4}n)$:
we have $\log\mathsf{SNR}\gtrsim\log\log n$;

  * •

Hard regime
$(c(\varepsilon)\lesssim\textup{srank}(\mathbf{B}^{\natural})\ll\log n)$: we
have

$\displaystyle\log\mathsf{SNR}\gtrsim\frac{\log
n}{\textup{srank}(\mathbf{B}^{\natural})}+\log\log n,$ (8)

where $\varepsilon>0$ is an arbitrary positive constant and $c(\varepsilon)$
is a positive constant depending on $\varepsilon$.

Family of log-concave sub-gaussian distributions include standard Gaussian
distribution and uniform distribution among $[-c,~{}c]$, where $c$ denotes a
certain positive constant. However, not all sub-gaussian RVs are with log-
concavity. One example is that $\mathbf{X}_{ij}$ is a Rademacher RV, i.e.,
$\mathbb{P}(\mathbf{X}_{ij}=\pm 1)=\nicefrac{{1}}{{2}}$. We can see that
Theorem 5 requires log-concavity and can allow a much broader range of
$\textup{srank}(\mathbf{B}^{\natural})$ than Theorem 4, to put more specific,
a positive constant $\textup{srank}(\mathbf{B}^{\natural})$ is permitted. In
the following, we will see some numerical results implying that the log-
concavity may be inseparable from the positive constant
$\textup{srank}(\mathbf{B}^{\natural})$, in other words,
$\textup{srank}(\mathbf{B}^{\natural})$ cannot be $O(1)$ if the log-concavity
assumption is violated.

###  4.1 Results discussion

#### Comparison with single observation model.

The most noticeable implications of our theorems are that the $\mathsf{SNR}$
requirement for correct permutation recovery can be greatly reduced by making
multiple measurements, i.e., from $\Omega(n^{c})$ to
$\Omega(n^{\frac{c\cdot\textup{srank}(\mathbf{B}^{\natural})}{n}}\vee\log^{c}n)$
or even $\Omega(1)$.

Additionally, we notice that the constraint on signal length $p$ has been
removed when multiple measurements are made. According to Theorem 3 which
concerns the single observation model $(m=1)$, our estimator only works when
$p=1$. While in Theorem 4 and Theorem 5 where
$\textup{srank}(\mathbf{B}^{\natural})\gg 1$, a much broader range of $p$ is
permitted.

In summary, the take-home message is that diversity, i.e., large
$\textup{srank}(\mathbf{B}^{\natural})$, helps in the permutation recovery.

#### Statistical optimality.

First, we argue that our estimator in Algorithm 1 is almost mini-max optimal
w.r.t. $\mathsf{SNR}$ requirement. For an easy comparison, we consider the
special case where $\textup{srank}(\mathbf{B}^{\natural})$ is same as
$\textup{rank}(\mathbf{B}^{\natural})$, which corresponds to the situation
where $\mathbf{B}^{\natural}$’s signal strength is uniformly spread among all
its eigenvalues. Theorem 2 suggests that wrong permutation matrix
$\widehat{\bm{\Pi}}$ will be returned with a high probability if

$\displaystyle\textup{srank}(\mathbf{B}^{\natural})\cdot\log\mathsf{SNR}\approx\
\textup{srank}(\mathbf{B}^{\natural})\cdot\log(1+\mathsf{SNR})\lesssim\log
n,\vspace{1mm}$

which holds regardless of the estimator form. Comparing with Theorem 5, we
conclude that our estimator has optimal $\mathsf{SNR}$ requirement in the
following two regimes: $(i)$
$\textup{srank}(\mathbf{B}^{\natural})\lesssim\frac{\log n}{\log\log n}$ and
$(ii)$ $\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$. Only in the regime
$\nicefrac{{\log n}}{{\log\log
n}}\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log^{4}n$ our estimator
experiences a loss in the $\mathsf{SNR}$ requirement, which is up to
$O(\log^{c}n)$.

Second, our requirement on sample number $n$ is almost optimal: Unnikrishnan
et al. (2015) claims that $n\gtrsim p$ is required for correct permutation
recovery while our estimator only needs $n\gg
p\cdot\log^{3}n\cdot\log^{2}(n^{2}p^{3})$. In addition, we have almost minimal
constraint on permuted rows’ number, i.e., we allow $h_{\textup{max}}\asymp
n$, where $h_{\textup{max}}$ is the maximum allowed number of permuted rows.

###  4.2 Proof outline

To make the proof more digestible, we first outline the proof strategies and
technical challenges before delving into the technical details. The rigorous
proof is attached in the supplementary material, including all supporting
Lemmas. Here, we would like to explain the main technical challenges in the
proof of Theorem 4 and Theorem 5, which lies in the proof that

$\displaystyle\left\\{\left\langle\bm{\Pi},\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}\right\rangle\geq\langle\bm{\Pi}^{\natural},\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}\rangle,~{}~{}\exists~{}\bm{\Pi}\neq\bm{\Pi}^{\natural}\right\\}$
(9)

holds with probability near zero given the assumptions in Theorem 4 and
Theorem 5. Two noticeable challenges in analyzing the above events are $(i)$
combinatorial nature of the optimization problem; and $(ii)$ high-order
moments of $\mathbf{X}$ in the product
$\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}$. To address these
challenges, we construct our solutions with two building blocks: $(i)$
relaxations of error event and $(ii)$ modified leave-one-out techniques. The
following context presents a more detailed explanation.

To begin with, we define $\widehat{\mathbf{B}}$ and $\widetilde{\mathbf{B}}$
respectively as

$\displaystyle\widetilde{\mathbf{B}}$
$\displaystyle=\left(n-h\right)^{-1}\mathbf{X}^{\top}\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural},$
$\displaystyle\widehat{\mathbf{B}}$
$\displaystyle=\left(n-h\right)^{-1}\mathbf{X}^{\top}\mathbf{Y}=\widetilde{\mathbf{B}}+\
\left(n-h\right)^{-1}\mathbf{X}^{\top}\mathbf{W},\vspace{1mm}$

where $h$ is denoted as the Hamming distance between identity matrix
$\mathbf{I}$ and the ground truth permutation matrix $\bm{\Pi}^{\natural}$,
i.e., $h=\mathsf{d_{H}}(\mathbf{I},\bm{\Pi}^{\natural})$.

#### Stage I. Relaxation of error event.

To combat the combinatorial nature of (9), we first relax the error event to
make it amenable for analysis. Associated with different regimes comes
different forms of relaxations.

  * •

Easy regime. We relax the error event
$\\{\bm{\Pi}^{\textup{opt}}\neq\bm{\Pi}^{\natural}\\}$ as

$\displaystyle\left\\{\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\rangle\leq\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{j,:}\right\rangle,~{}\exists~{}1\leq\pi^{\natural}(i)\neq
j\leq n\right\\}.$ (10)

With this relaxation method, we find the constraint
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{2}n$ to be inevitable. However,
as compensation, we can show $\mathsf{SNR}\geq c$ is sufficient for the
correct permutation reconstruction.

  * •

Medium $\&$ hard regime. To get rid of the constraint on
$\textup{srank}(\mathbf{B}^{\natural})$, we first exploit the energy-
preserving property permutation matrix, i.e.,
${\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\bm{\Pi}\mathbf{M}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}={\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{M}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$
where $\mathbf{M}\in\mathbb{R}^{n\times(\cdot)}$ is an arbitrary matrix. Then
we adopt the relaxation

$\displaystyle\left\\{\big{\|}\mathbf{Y}_{i,:}-\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\big{\|}_{2}^{2}\geq\big{\|}\mathbf{Y}_{i,:}-\widehat{\mathbf{B}}^{\top}\mathbf{X}_{j,:}\big{\|}_{2}^{2},~{}~{}\exists~{}i,j\right\\}.$
(11)

In this way, we can relax the constraint on
$\textup{srank}(\mathbf{B}^{\natural})$ from
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{2}n$ to
$\textup{srank}(\mathbf{B}^{\natural})\gg\log n$ if $\mathbf{X}_{ij}$ is
i.i.d. isotropic sub-gaussian RV. If the log-concavity assumption is further
put on $\mathbf{X}_{ij}$’s distribution, we can relax the constraint to
$\textup{srank}(\mathbf{B}^{\natural})\gg c$.

However, compared with the above relaxation method, this method will
experience some loss in the $\mathsf{SNR}$ requirement. For instance, in the
Easy Regime, this approach will require $\mathsf{SNR}$ to satisfy
$\mathsf{SNR}\gtrsim\log^{c}n$ for the correct permutation recovery while the
previous approach only needs $\mathsf{SNR}\geq c$.

Bonuses of the above relaxations include $(i)$ we reduce the moments of
$\mathbf{X}$ from fourth order in (9) to third order, see (10) and (11); and
$(ii)$ the dependences between $\widehat{\mathbf{B}}$ and other terms in (10)
and (11) are restricted to the rows $\mathbf{X}_{\pi^{\natural}(i),:}$ and
$\mathbf{X}_{j,:}$. The latter bonus is essential for the leave-one-out
technique.

#### Stage II. Dependence decoupling with leave-one-out technique.

To further reduce the moments of $\mathbf{X}$ involved in the analysis, we
modify the _leave-one-out_ technique (Karoui, 2013; Karoui et al., 2013;
Karoui, 2018; Chen et al., 2020; Sur et al., 2019; Zhang and Li, 2020).

As mentioned above, $\widehat{\mathbf{B}}$ is weakly correlated with rows
$\mathbf{X}_{\pi^{\natural}(i),:},\mathbf{X}_{j,:}$, $\mathbf{Y}_{i,:}$, and
$\mathbf{Y}_{j,:}$. The basic idea of the leave-one-out technique is to
replace the correlated rows in $\widehat{\mathbf{B}}$ with their i.i.d.
copies, by which the dependence is decoupled. Since only a limited number of
rows are replaced, we expect its statistical properties should remain almost
identical. A detailed explanation comes as follows.

First, we draw an independent copy $\mathbf{X}_{s,:}^{{}^{\prime}}$ for each
row $\mathbf{X}_{s,:}$ ($s$th row of the sensing matrix $\mathbf{X}$). With
these independent copies, we construct leave-one-out samples
$\widetilde{\mathbf{B}}_{\setminus(s)}$ $(1\leq s\leq n)$ by replacing the
$s$th row in $\widetilde{\mathbf{B}}$ with its independent copy
$\mathbf{X}_{s,:}^{{}^{\prime}}$. In formulae:

$\displaystyle\widetilde{\mathbf{B}}_{\setminus(s)}=(n-h)^{-1}\bigg{(}\sum_{\begin{subarray}{c}k\neq
s\\\ \pi^{\natural}(k)\neq
s\end{subarray}}\mathbf{X}_{\pi^{\natural}(k),:}\mathbf{X}_{k,:}^{\top}+\sum_{\begin{subarray}{c}k=s\textup{
or}\\\
\pi^{\natural}(k)=s\end{subarray}}{\mathbf{X}}^{{}^{\prime}}_{\pi^{\natural}(k),:}{\mathbf{X}}^{{}^{\prime}\top}_{k,:}\bigg{)}\mathbf{B}^{\natural}.\vspace{1mm}$

Easily we can verify that $\widetilde{\mathbf{B}}_{\setminus(s)}$ is
independent of $\mathbf{X}_{s,:}$. Similarly, we construct matrices
$\\{\widetilde{\mathbf{B}}_{\setminus(s,t)}\\}_{1\leq s\neq t\leq n}$ as

$\displaystyle\widetilde{\mathbf{B}}_{\setminus(s,t)}=(n-h)^{-1}\bigg{(}\sum_{\begin{subarray}{c}k\neq
s,t\\\ \pi^{\natural}(k)\neq
s,t\end{subarray}}\mathbf{X}_{\pi^{\natural}(k),:}\mathbf{X}_{k,:}^{\top}+\sum_{\begin{subarray}{c}k=s\textup{
or }k=t\textup{ or}\\\ \pi^{\natural}(k)=s\textup{ or
}\pi^{\natural}(k)=t\end{subarray}}{\mathbf{X}}^{{}^{\prime}}_{\pi^{\natural}(k),:}{\mathbf{X}}^{{}^{\prime}\top}_{k,:}\bigg{)}\mathbf{B}^{\natural},\vspace{1mm}$

and can verify the independence between
$\widetilde{\mathbf{B}}_{\setminus(s,t)}$ and the rows $\mathbf{X}_{s,:}$ and
$\mathbf{X}_{t,:}$.

###  4.3 Discussion of related work

This subsection compares our results with the previous works on multiple
observations model ($m>1$) (Pananjady et al., 2017; Zhang et al., 2022;
Slawski et al., 2020). Pananjady et al. (2017) consider a similar setting as
ours while their focus is on the product
$\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural}$ rather than the
individual values of $\bm{\Pi}^{\natural}$ and $\mathbf{B}^{\natural}$.

In Zhang et al. (2022), statistical limits w.r.t. the $\mathsf{SNR}$ are
presented and the ML estimator is revisited under the setting of multiple
observations model. It is suggested that the ML estimator can reach the
statistical limits under certain regimes. In addition, a projected gradient
descent-based algorithm with monotonic descent property is proposed for
practical use. However, whether this algorithm will yield the ground truth
$\bm{\Pi}^{\natural}$ still remains a mystery.

In Slawski et al. (2020), they follow a similar idea of Slawski and Ben-David
(2019) and take the viewpoint of denoising for permutation recovery. Assuming
only a limited proportion of rows are permuted, they first obtain an estimate
of $\mathbf{B}^{\natural}$ and plug it into the ML estimator in (3) to
reconstruct $\bm{\Pi}^{\natural}$. Compared with our estimator in Algorithm 1,
their estimator can be applied to a broader family of matrices (i.e.,
selection matrices); but $(i)$ allow a smaller number of permuted rows and
$(ii)$ require a larger sample number $n$. A detailed comparison is referred
to Table 1.

##  5 Simulation Results

This section presents the numerical results. Since our estimator cannot
guarantee the correct permutation recovery under the single observation model,
our simulations focus on the multiple observations model, i.e., $m>1$.

We investigate the impact of the $\nicefrac{{p}}{{n}}$ ratio and number of
permuted rows $h$ on the permutation recovery when $\mathbf{X}_{ij}$ is $(i)$
Gaussian distributed, $(ii)$ uniformly distributed within $[-1,1]$, and
$(iii)$ Rademacher distributed such that $\mathbb{P}(\mathbf{X}_{ij}=\pm
1)=\nicefrac{{1}}{{2}}$.

First, we present our experiment settings. We set the $i$th column
$\mathbf{B}^{\natural}_{:,i}$ $(1\leq i\leq\min(m,p))$ to be the $i$th
canonical basis, which has $1$ on the $i$th entry and $0$ elsewhere. One
benefit of this setting is that the stable rank
$\textup{srank}(\mathbf{B}^{\natural})$ can be easily calculated, i.e.,
$m\wedge p$. In the following context, we will use
$\textup{srank}(\mathbf{B}^{\natural})$ and $m$ $(m<p)$ interchangeably.

Figure 3: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$, with $n=500$,
$p\in\\{50,100,150\\}$, and $h\in\\{50,125\\}$, versus $\mathsf{SNR}$ (upper
panels) and
$\frac{\log\det\left(\mathbf{I}+\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}/\sigma^{2}\right)}{\log
n}$ (lower panels).

###  5.1 Gaussian distribution

We assume $\mathbf{X}_{ij}$ to be i.i.d. standard normal distribution
$\mathsf{N}(0,1)$. We fix $\nicefrac{{p}}{{n}}$ to be
$\left\\{0.1,0.2,0.3\right\\}$ and $\nicefrac{{h}}{{n}}$ to be
$\left\\{0.1,0.25\right\\}$. We vary the sample number $n$ to be
$\left\\{500,1000,2000,4000\right\\}$. The corresponding results are put in
Figures 3, 4, 5, and 6, respectively. Apart from evaluating the permutation
recovery w.r.t. $\mathsf{SNR}$, we also evaluate it w.r.t.
$\nicefrac{{\log\det\left(\mathbf{I}+\frac{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}{\sigma^{2}}\right)}}{{\log
n}}$. This ratio appears in the statistical lower bounds in Theorem 2, which
claims

$\displaystyle\log\det\bigg{(}\mathbf{I}+\frac{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}{\sigma^{2}}\bigg{)}\gtrsim\log
n,\vspace{1mm}$

should hold for correct permutation recovery.

First, we notice the numerical results are well aligned with our theoretical
results, say Theorem 5, which concludes the optimality of our estimator.
Moreover, we find that smaller $\nicefrac{{p}}{{n}}$ and $h$ facilitate the
permutation recovery. The former factor leads to more concentrated behavior of
$\mathbf{X}^{\top}\mathbf{Y}$ around its means
$\mathbb{E}\mathbf{X}^{\top}\mathbf{Y}=(n-h)\mathbf{B}^{\natural}\parallel\mathbf{B}^{\natural}$;
while the latter factor brings more energy in
$\mathbb{E}\mathbf{X}^{\top}\mathbf{Y}$.

Figure 4: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$, with $n=1000$,
$p\in\\{100,200,300\\}$, and $h\in\\{100,250\\}$, versus $\mathsf{SNR}$ (upper
panels) and
$\frac{\log\det\left(\mathbf{I}+\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}/\sigma^{2}\right)}{\log
n}$ (lower panels).

Figure 5: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$, with $n=2000$,
$p\in\\{200,400,600\\}$, and $h\in\\{200,500\\}$, versus $\mathsf{SNR}$ (upper
panels) and
$\frac{\log\det\left(\mathbf{I}+\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}/\sigma^{2}\right)}{\log
n}$ (lower panels).

Figure 6: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$, with $n=4000$,
$p\in\\{400,800,1200\\}$, and $h\in\\{400,1000\\}$, versus $\mathsf{SNR}$
(upper panels) and
$\frac{\log\det\left(\mathbf{I}+\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}/\sigma^{2}\right)}{\log
n}$ (lower panels).

Figure 7: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$ with $n=500$,
$p\in\left\\{50,100,150\right\\}$, $h\in\left\\{50,125\right\\}$, and
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Unif}[-1,1]$,
with respect to $\mathsf{SNR}$.

Figure 8: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$ with $n=1000$,
$p\in\left\\{100,200,300\right\\}$, $h\in\left\\{100,250\right\\}$, and
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Unif}[-1,1]$,
with respect to $\mathsf{SNR}$.

###  5.2 Uniform distribution

This subsection investigates the recovery performance when
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Unif}[-1,1]$.
Similar to above, we fix $\nicefrac{{p}}{{n}}$ to be
$\left\\{0.1,0.2,0.3\right\\}$ and $\nicefrac{{h}}{{n}}$ to be
$\left\\{0.1,0.25\right\\}$. Since the performance does not change drastically
with increasing $n$, we limit sample number $n$ to be
$\left\\{500,1000\right\\}$ and put the corresponding results in Figure 7 and
Figure 8, respectively. For an easy comparison, we also put the results with
Gaussian distributed $\mathbf{X}$ in Figure 7 and Figure 8, which are shown in
dashed curves. One noteworthy fact is that uniform distribution seems to be
more friendly for permutation recovery. A plausible reason is that uniform
distribution has a smaller variance, namely,
$\textup{Var}(\cdot)=(\nicefrac{{1}}{{2}})\int_{-1}^{1}z^{2}dz=\nicefrac{{1}}{{3}}<1$.

Apart from that, we see its behavior is very similar to that of the Gaussian
distribution. For the conciseness of presentation, we omit the plot of
recovery rate in terms of
$\frac{\log\det(\mathbf{I}+\nicefrac{{\mathbf{B}^{\natural\top}\mathbf{B}^{\natural}}}{{\sigma^{2}}})}{\log
n}$.

###  5.3 Rademacher distribution

This subsection considers the Rademacher distribution, i.e.,
$\mathbb{P}(\mathbf{X}_{ij}=\pm 1)=\nicefrac{{1}}{{2}}$. To begin with, we
briefly comment on the $\textup{srank}(\mathbf{B}^{\natural})$ requirement. In
Theorem 5, we claim that log-concavity is required to remove the requirement
$\textup{srank}(\mathbf{B}^{\natural})\gtrsim\log n$. Here we give an example
to illustrate its necessity. For details see Figure 9. In the left panel, we
assume
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Unif}[-1,1]$,
which is a log-concave sub-gaussian RV; while in the right panel, we assume
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Rademacher}$,
which is sub-gaussian but not log-concave. When the stable rank
$\textup{srank}(\mathbf{B}^{\natural})$ is not sufficiently large, to put more
specifically, $\textup{srank}(\mathbf{B}^{\natural})\gtrsim\log n$, we
conclude the correct permutation cannot be reconstructed for Rademacher
distribution even when $\mathsf{SNR}$ is sufficiently large.

Figure 9: Comparison of permutation recovery performance between Uniform
distribution among $[-1,1]$ (left panel) and Rademacher distribution (right
panel). The dashed line corresponds to the performance when
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\mathsf{N}(0,1)$.

In addition, we investigate the impact of $\nicefrac{{p}}{{n}}$ and number of
permuted rows $h$. Similar to above, we fix $\nicefrac{{p}}{{n}}$ to be
$\left\\{0.1,0.2,0.3\right\\}$ and $\nicefrac{{h}}{{n}}$ to be
$\left\\{0.1,0.25\right\\}$. In Figure 10, we fix $n$ to be $500$; while in
Figure 11, we increase it to be $1000$. One noticeable phenomenon is that the
permutation recovery performance under the Rademacher setting is very similar
to that of Gaussian distribution and Uniform distribution provided that
$\textup{srank}(\mathbf{B}^{\natural})$ is large enough. This again suggests
that log-concavity assumption is unavoidable if the large stable rank
requirement is to be mitigated.

Figure 10: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$ with $n=500$,
$p\in\left\\{50,100,150\right\\}$, $h\in\left\\{50,125\right\\}$, and
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Rademacher}$,
with respect to $\mathsf{SNR}$. The dashed line corresponds to the performance
when
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\mathsf{N}(0,1)$.

Figure 11: Simulated recovery rate
$\mathbb{P}(\widehat{\bm{\Pi}}=\bm{\Pi}^{\natural})$ with $n=1000$,
$p\in\left\\{100,200,300\right\\}$, $h\in\left\\{100,250\right\\}$, and
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\textup{Rademacher}$,
with respect to $\mathsf{SNR}$. The dashed line corresponds to the performance
when
$\mathbf{X}_{ij}\stackrel{{\scriptstyle\textup{i.i.d}}}{{\sim}}\mathsf{N}(0,1)$.

##  6 Concluding Remarks

This paper considers unlabeled linear regression and proposes a one-step
estimator, which is optimal in both computational and statistical
perspectives. First, we show our estimator has the same computational
complexity as that of oracle estimators. Afterwards, we separately investigate
its statistical properties under the single observation model ($m=1$) and
multiple observations model ($m>1$). For the single observation model, our
estimator can obtain the ground truth permutation matrix when the signal is
with length one, i.e., $p=1$. Moreover, its $\mathsf{SNR}$ requirement matches
the minimax lower bound, to put more specifically,
$\log\mathsf{SNR}\gtrsim\log n$. For the multiple observations model, our
estimator exhibits much richer behavior, which can be broadly divided into
three regimes. In the Easy Regime
($\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$), our estimator is
minimax optimal and can obtain the correct permutation once
$\mathsf{SNR}\geq\Omega(1)$; in the Medium Regime ($\log
n\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log^{4}n$), our estimator
experiences a performance loss up to a multiplicative polynomial of $\log n$
as it requires $\mathsf{SNR}\gtrsim\log n$ rather than
$\mathsf{SNR}\gtrsim\Omega(1)$ for correct permutation recovery; and in the
hard regime ($1\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log n$), our
estimator needs $\mathsf{SNR}\gtrsim(\log n)^{c_{0}}\cdot
n^{\nicefrac{{c_{1}}}{{\textup{srank}(\mathbf{B}^{\natural})}}}$, which
matches the statistical limits when
$\textup{srank}(\mathbf{B}^{\natural})\ll\frac{\log n}{\log\log n}$ and
experiences a loss up to a multiplicative polynomial of $\log n$ when
$\frac{\log n}{\log\log n}\ll\textup{srank}(\mathbf{B}^{\natural})\ll\log n$.
Comparing with the single observation model, we conclude that high diversity,
namely, large $\textup{srank}(\mathbf{B}^{\natural})$, can greatly facilitate
the permutation recovery. Ultimately, we provide numerical experiments to
support all claims thereof.

## References

  * Abid et al. (2017) Abubakar Abid, Ada Poon, and James Zou.  Linear regression with shuffled labels.  _arXiv preprint arXiv:1705.01342_ , 2017.
  * Bai and Hsing (2005) Zhidong Bai and Tailen Hsing.  The broken sample problem.  _Probability Theory and Related Fields_ , 131(4):528–552, 2005.
  * Balakrishnan et al. (2017) Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu.  Statistical guarantees for the EM algorithm: From population to sample-based analysis.  _The Annals of Statistics_ , 45(1):77–120, 02 2017.
  * Bertsekas and Castañón (1992) Dimitri P. Bertsekas and David A. Castañón.  A forward/reverse auction algorithm for asymmetric assignment problems.  _Comput. Optim. Appl._ , 1(3):277–297, 1992.
  * Candès and Plan (2010) Emmanuel J. Candès and Yaniv Plan.  Matrix completion with noise.  _Proc. IEEE_ , 98(6):925–936, 2010.
  * Chen et al. (2020) Yuxin Chen, Yuejie Chi, Jianqing Fan, Cong Ma, and Yuling Yan.  Noisy matrix completion: Understanding statistical guarantees for convex relaxation via nonconvex optimization.  _SIAM J. Optim._ , 30(4):3098–3121, 2020.
  * Chi et al. (2019) Yuejie Chi, Yue M. Lu, and Yuxin Chen.  Nonconvex optimization meets low-rank matrix factorization: An overview.  _IEEE Trans. Signal Process._ , 67(20):5239–5269, 2019.
  * DeGroot and Goel (1976) Morris H. DeGroot and Prem K. Goel.  The matching problem for multivariate normal data.  _Sankhyā: The Indian Journal of Statistics, Series B (1960-2002)_ , 38(1):14–29, 1976.
  * DeGroot and Goel (1980) Morris H. DeGroot and Prem K. Goel.  Estimation of the correlation coefficient from a broken random sample.  _The Annals of Statistics_ , 8(2):264–278, 03 1980.
  * Dokmanic (2019) Ivan Dokmanic.  Permutations unlabeled beyond sampling unknown.  _IEEE Signal Process. Lett._ , 26(6):823–827, 2019.
  * Fang and Li (2023) Guanhua Fang and Ping Li.  Regression with label permutation in generalized linear model.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , pages 9716–9760, Honolulu, HI, 2023.
  * Goel (1975) Prem K. Goel.  On re-pairing observations in a broken random sample.  _The Annals of Statistics_ , 3(6):1364–1369, 11 1975.
  * Golub and Van Loan (2012) Gene H Golub and Charles F Van Loan.  _Matrix computations_ , volume 3.  JHU press, 2012.
  * Haghighatshoar and Caire (2018) Saeid Haghighatshoar and Giuseppe Caire.  Signal recovery from unlabeled samples.  _IEEE Trans. Signal Process._ , 66(5):1242–1257, 2018.
  * Hsu et al. (2017) Daniel J. Hsu, Kevin Shi, and Xiaorui Sun.  Linear regression without correspondence.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 1531–1540, Long Beach, CA, 2017.
  * Jeong et al. (2020) Min-Oh Jeong, Alex Dytso, Martina Cardone, and H. Vincent Poor.  Recovering data permutations from noisy observations: The linear regime.  _IEEE J. Sel. Areas Inf. Theory_ , 1(3):854–869, 2020.
  * Karoui (2013) Noureddine El Karoui.  Asymptotic behavior of unregularized and ridge-regularized high-dimensional robust regression estimators: rigorous results.  _arXiv preprint arXiv:1311.2445_ , 2013.
  * Karoui (2018) Noureddine EL Karoui.  On the impact of predictor geometry on the performance on high-dimensional ridge-regularized generalized robust regression estimators.  _Probability Theory and Related Fields_ , 170(1-2):95–175, 2018.
  * Karoui et al. (2013) Noureddine El Karoui, Derek Bean, Peter J Bickel, Chinghway Lim, and Bin Yu.  On robust regression with high-dimensional predictors.  _Proceedings of the National Academy of Sciences_ , 110(36):14557–14562, 2013.
  * Kuhn (1955) Harold W Kuhn.  The Hungarian method for the assignment problem.  _Naval research logistics quarterly_ , 2(1-2):83–97, 1955.
  * Latala et al. (2007) Rafal Latala, Piotr Mankiewicz, Krzysztof Oleszkiewicz, and Nicole Tomczak-Jaegermann.  Banach-mazur distances and projections on random subgaussian polytopes.  _Discret. Comput. Geom._ , 38(1):29–50, 2007\.
  * Pananjady et al. (2017) Ashwin Pananjady, Martin J Wainwright, and Thomas A Courtade.  Denoising linear models with permuted data.  In _Proceedings of the 2017 IEEE International Symposium on Information Theory (ISIT)_ , pages 446–450, Aachen, Germany, 2017.
  * Pananjady et al. (2018) Ashwin Pananjady, Martin J. Wainwright, and Thomas A. Courtade.  Linear regression with shuffled data: Statistical and computational limits of permutation recovery.  _IEEE Trans. Inf. Theory_ , 64(5):3286–3300, 2018.
  * Paouris (2012) Grigoris Paouris.  Small ball probability estimates for log-concave measures.  _Transactions of the American Mathematical Society_ , 364(1):287–308, 2012.
  * Peng et al. (2021) Liangzu Peng, Boshi Wang, and Manolis C. Tsakiris.  Homomorphic sensing: Sparsity and noise.  In _Proceedings of the 38th International Conference on Machine Learning (ICML)_ , pages 8464–8475, Virtual Event, 2021.
  * Slawski and Ben-David (2019) Martin Slawski and Emanuel Ben-David.  Linear regression with sparsely permuted data.  _Electronic Journal of Statistics_ , 1:1–36, 2019.
  * Slawski and Sen (2022) Martin Slawski and Bodhisattva Sen.  Permuted and unlinked monotone regression in $\mathbb{R}^{d}$: an approach based on mixture modeling and optimal transport.  _arXiv preprint arXiv:2201.03528_ , 2022.
  * Slawski et al. (2020) Martin Slawski, Emanuel Ben-David, and Ping Li.  Two-stage approach to multivariate linear regression with sparsely mismatched data.  _J. Mach. Learn. Res._ , 21:204:1–204:42, 2020.
  * Sur et al. (2019) Pragya Sur, Yuxin Chen, and Emmanuel J Candès.  The likelihood ratio test in high-dimensional logistic regression is asymptotically a rescaled chi-square.  _Probability Theory and Related Fields_ , 175(1):487–558, 2019.
  * Tang et al. (2021) Zhiwei Tang, Tsung-Hui Chang, Xiaojing Ye, and Hongyuan Zha.  Low-rank matrix recovery with unknown correspondence.  _arXiv preprint arXiv:2110.07959_ , 2021.
  * Tropp (2015) Joel A. Tropp.  An introduction to matrix concentration inequalities.  _Found. Trends Mach. Learn._ , 8(1-2):1–230, 2015\.
  * Tsakiris and Peng (2019) Manolis C. Tsakiris and Liangzu Peng.  Homomorphic sensing.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning (ICML)_ , pages 6335–6344, Long Beach, CA, 2019.
  * Unnikrishnan et al. (2015) Jayakrishnan Unnikrishnan, Saeid Haghighatshoar, and Martin Vetterli.  Unlabeled sensing: Solving a linear system with unordered measurements.  In _Proceedings of the 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 786–793, Monticello, IL, 2015.
  * Vershynin (2018) Roman Vershynin.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Zhang and Li (2020) Hang Zhang and Ping Li.  Optimal estimator for unlabeled linear regression.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 11153–11162, Virtual Event, 2020.
  * Zhang and Li (2023a) Hang Zhang and Ping Li.  Greed is good: correspondence recovery for unlabeled linear regression.  In _Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 2509–2518, Pittsburgh, PA, 2023a.
  * Zhang and Li (2023b) Hang Zhang and Ping Li.  One-step estimator for permuted sparse recovery.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , pages 41244–41267, Honolulu, HI, 2023b.
  * Zhang et al. (2022) Hang Zhang, Martin Slawski, and Ping Li.  The benefits of diversity: Permutation recovery in unlabeled sensing from multiple measurement vectors.  _IEEE Trans. Inf. Theory_ , 68(4):2509–2529, 2022.

##  Appendix A Appendix for Section 3

This section focuses on the special case where $p=1$ and $m=1$. Consider
$\bm{x}\in\mathbb{R}^{n}$ to be an isotropic log-concave random vector with
zero mean and ${\left\|\bm{x}\right\|}_{\psi_{2}}\lesssim 1$. Additionally, we
assume that the permutation matrix $\bm{\Pi}^{\natural}$ satisfies
$\mathsf{d_{H}}(\mathbf{I},\bm{\Pi}^{\natural})=h\leq\nicefrac{{n}}{{4}}$.

###  A.1 Notations: single observation model

First, we define the following events ${\cal E}_{i}$, $(1\leq i\leq 4)$,

$\displaystyle{\cal E}_{1}\triangleq$
$\displaystyle\left\\{\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle\geq
c_{0}n\right\\};$ $\displaystyle{\cal E}_{2}\triangleq$
$\displaystyle\left\\{\bm{w}^{\top}\bm{x}\bm{x}^{\top}(\bm{\Pi}^{\natural}-\bm{\Pi})\bm{w}\lesssim\sigma^{2}n^{2}\log
n,~{}\forall~{}\bm{\Pi}\neq\bm{\Pi}^{\natural}\right\\};$ $\displaystyle{\cal
E}_{3}\triangleq$
$\displaystyle\left\\{\left\lvert\langle\bm{w},\bm{x}\rangle\langle\bm{\Pi}^{\natural}\bm{x},(\bm{\Pi}^{\natural}-\bm{\Pi})^{\top}\bm{x}\rangle+\
\langle\bm{w},(\bm{\Pi}^{\natural}-\bm{\Pi})^{\top}\bm{x}\rangle\langle\bm{\Pi}^{\natural}\bm{x},\bm{x}\rangle\right\rvert\lesssim\sigma
n^{2}\sqrt{\log n},~{}\forall~{}\bm{\Pi}\neq\bm{\Pi}^{\natural}\right\\};$
$\displaystyle{\cal E}_{4}\triangleq$
$\displaystyle\left\\{{\left\|\bm{x}-\bm{\Pi}\bm{x}\right\|}_{2}^{2}\gtrsim
n^{-20},~{}~{}\forall~{}\bm{\Pi}\neq\bm{\Pi}^{\natural}\right\\}.$

###  A.2 Proof of Theorem 3

###### Proof.

Under the assumptions in Theorem 3, we will prove that the ground truth
permutation matrix $\bm{\Pi}^{\natural}$ will be returned with high
probability with Algorithm 1. Ahead of the technical details, we give an
outline of the proof strategy, which can be divided into two stages.

  * •

Stage I. We show the intersection of events $\bigcap_{i=1}^{5}{\cal E}_{i}$ is
a subset of the event $\\{\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural}\\}$
under the assumptions of Theorem 3.

  * •

Stage II. With the union bound, we can upper-bound the error probability
$\mathbb{P}(\bm{\Pi}^{\textup{opt}}\neq\bm{\Pi}^{\natural})$ by
$\sum_{\ell=1}^{4}\mathbb{P}(\overline{{\cal E}}_{\ell})$. The proof is then
completed by studying probability $\mathbb{P}(\overline{{\cal E}}_{(\cdot)})$,
respectively.

Technical details come as follows.

Stage I. We begin the proof by showing $\bigcap_{\ell=1}^{4}{\cal
E}_{\ell}\subseteq\\{\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural}\\}$ under the
assumptions of Theorem 3. First, we expand
$\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle$ as

$\displaystyle\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle=(\beta^{\natural})^{2}{\cal
T}_{2}(\bm{\Pi})+\beta^{\natural}{\cal T}_{1}(\bm{\Pi})+{\cal
T}_{0}(\bm{\Pi}),\vspace{1mm}$

where $\bm{\Pi}$ is an arbitrary permutation matrix, and ${\cal
T}_{i}(\bm{\Pi})$ $(0\leq i\leq 2)$ are defined as

$\displaystyle{\cal T}_{2}(\bm{\Pi})$
$\displaystyle=\langle\bm{\Pi}^{\natural}\bm{x},\bm{\Pi}\bm{x}\rangle\langle\bm{\Pi}^{\natural}\bm{x},\bm{x}\rangle;$
$\displaystyle{\cal T}_{1}(\bm{\Pi})$
$\displaystyle=\left\langle\bm{w},\bm{x}\right\rangle\langle\bm{\Pi}^{\natural}\bm{x},\bm{\Pi}^{\top}\bm{x}\rangle+\
\langle\bm{w},\bm{\Pi}^{\top}\bm{x}\rangle\langle\bm{\Pi}^{\natural}\bm{x},\bm{x}\rangle;$
$\displaystyle{\cal T}_{0}(\bm{\Pi})$
$\displaystyle=\langle\bm{w},\bm{\Pi}^{\top}\bm{x}\rangle\langle\bm{x},\bm{w}\rangle.\vspace{1mm}$

Then we can express the difference
$\langle\bm{\Pi}^{\natural},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle-\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle$
as

$\displaystyle\langle\bm{\Pi}^{\natural},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle-\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle$
$\displaystyle=~{}$ $\displaystyle\ (\beta^{\natural})^{2}({\cal
T}_{2}(\bm{\Pi}^{\natural})-{\cal T}_{2}(\bm{\Pi}))+\beta^{\natural}({\cal
T}_{1}(\bm{\Pi}^{\natural})-{\cal T}_{1}(\bm{\Pi}))+{\cal
T}_{0}(\bm{\Pi}^{\natural})-{\cal T}_{0}(\bm{\Pi})$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{=}}~{}$
$\displaystyle\frac{(\beta^{\natural})^{2}}{2}\langle\bm{\Pi}^{\natural}\bm{x},\bm{x}\rangle\big{\|}\bm{x}-\bm{\Pi}^{\natural\top}\bm{\Pi}\bm{x}\big{\|}_{2}^{2}+\beta^{\natural}({\cal
T}_{1}(\bm{\Pi}^{\natural})-{\cal T}_{1}(\bm{\Pi}))+{\cal
T}_{0}(\bm{\Pi}^{\natural})-{\cal T}_{0}(\bm{\Pi})$

where in \raisebox{-.8pt}{1}⃝ we exploit the fact such that $\bm{x}$ is a
vector and hence

$\displaystyle{\left\|\bm{x}\right\|}_{2}^{2}-\langle\bm{\Pi}^{\natural}\bm{x},\bm{\Pi}\bm{x}\rangle=\left({\left\|\bm{x}\right\|}_{2}^{2}+\big{\|}\bm{\Pi}^{\natural\top}\bm{\Pi}\bm{x}\big{\|}_{2}^{2}-2\langle\bm{\Pi}^{\natural}\bm{x},\bm{\Pi}\bm{x}\rangle\right)/2=\
\big{\|}\bm{x}-\bm{\Pi}^{\natural\top}\bm{\Pi}\bm{x}\big{\|}_{2}^{2}/2.\vspace{1mm}$

Note that this relation generally does not hold except for the vector case.
Conditioning on ${\cal E}_{1}\bigcap{\cal E}_{4}$, we have the relation ${\cal
T}_{2}(\bm{\Pi}^{\natural})-{\cal T}_{2}(\bm{\Pi})\gtrsim\ell\cdot n^{-19}$
and hence

$\displaystyle\langle\bm{\Pi}^{\natural},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle-\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle\geq~{}$
$\displaystyle\frac{(\beta^{\natural})^{2}}{2}c_{0}nn^{-20}-\
\beta^{\natural}|{\cal T}_{1}(\bm{\Pi}^{\natural})-{\cal T}_{1}(\bm{\Pi})|-\
|{\cal T}_{0}(\bm{\Pi}^{\natural})-{\cal T}_{0}(\bm{\Pi})|$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\gtrsim}}~{}$
$\displaystyle\frac{c_{0}(\beta^{\natural})^{2}}{n^{19}}-c_{1}\beta^{\natural}\sigma
n^{2}\sqrt{\log n}-c_{2}\sigma^{2}n^{2}\log n,\vspace{1mm}$

where in \raisebox{-.8pt}{2}⃝ we condition on ${\cal E}_{2}\bigcap{\cal
E}_{3}$. Under the assumptions in Theorem 3, i.e.,
$\log\mathsf{SNR}\gtrsim\log n$, we have

$\displaystyle\langle\bm{\Pi}^{\natural},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle>\langle\bm{\Pi},\bm{y}\bm{y}^{\top}\bm{x}\bm{x}^{\top}\rangle,~{}~{}\forall~{}~{}\bm{\Pi}\neq\bm{\Pi}^{\natural},\vspace{1mm}$

which suggests the correct permutation $\bm{\Pi}^{\natural}$ can always be
obtained and we can upper bound the error probability
$\mathbb{P}(\bm{\Pi}^{\textup{opt}}\neq\bm{\Pi}^{\natural})$ by
$\sum_{\ell=1}^{4}\mathbb{P}(\overline{{\cal E}}_{\ell})$.

Stage II. We upper bound the error probability
$\mathbb{P}(\bm{\Pi}^{\textup{opt}}\neq\bm{\Pi}^{\natural})$ by
$\sum_{\ell=1}^{4}\mathbb{P}(\overline{{\cal E}}_{\ell})$ and complete the
proof with Lemma 1, Lemma 2, Lemma 3, and Lemma 4. ∎

###  A.3 Supporting lemmas for Theorem 3

This subsection collects the supporting lemmas for the proof of Theorem 3.

######  Lemma 1.

We have $\mathbb{P}\left({\cal E}_{1}\right)\geq 1-c_{0}e^{-c_{1}n}$ when $n$
is sufficiently large, where $c_{0},c_{1}>0$ are some positive constants.

###### Proof.

W.l.o.g, we assume the first $h$ entries are permuted and expand the inner
product $\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle$ as

$\displaystyle\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle=\sum_{i=1}^{h}x_{i}x_{\pi^{\natural}(i)}+\sum_{i=h+1}^{n}x_{i}^{2}.\vspace{1mm}$

With union bound, we can upper bound
$\mathbb{P}\left(\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle\leq
c_{0}n\right)$ as

$\displaystyle\mathbb{P}\left(\langle\bm{x},\bm{\Pi}^{\natural}\bm{x}\rangle\leq
c_{0}n\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}\
\underbrace{\mathbb{P}\left(\sum_{i=1}^{h}x_{i}x_{\pi^{\natural}(i)}\leq-\frac{c_{0}n}{2}\right)}_{\triangleq~{}\zeta_{1}}+\underbrace{\mathbb{P}\left(\sum_{i=h+1}^{n}x_{i}^{2}\leq\frac{3c_{0}n}{2}\right)}_{\triangleq~{}\zeta_{2}}.$
(12)

We finish the proof by separately upper-bounding $\zeta_{1}$ and $\zeta_{2}$
as $\zeta_{1},\zeta_{2}\lesssim e^{-c_{1}n}$. The detailed computation comes
as follows.

Analysis of $\zeta_{1}$. The technical difficulties stem from the correlation
between the terms $x_{i}x_{\pi^{\natural}(i)}$. According to Lemma 8 in
Pananjady et al. (2018) (restated as Lemma 21), we can divide the index set
$\left\\{i:i\neq\pi^{\natural}(i)\right\\}$ into $3$ disjoint categories
${\cal I}_{\ell}$, $(1\leq\ell\leq 3)$, such that (i) indices $i$ and
$\pi^{\natural}(i)$ belongs to different categories; (ii) the cardinality
$h_{\ell}$ of ${\cal I}_{\ell}$ satisfies $h_{\ell}\geq\lfloor h/5\rfloor$.
Then we make the decomposition

$\displaystyle\sum_{i=1}^{h}x_{i}x_{\pi^{\natural}(i)}=\sum_{\ell=1}^{3}\sum_{i\in{\cal
I}_{\ell}}x_{i}x_{\pi^{\natural}(i)},\vspace{1mm}$

and again, using the union bound, we have

$\displaystyle\zeta_{1}\leq~{}$
$\displaystyle\sum_{\ell=1}^{3}\mathbb{P}\bigg{(}\sum_{i\in{\cal
I}_{\ell}}x_{i}x_{\pi^{\natural}(i)}\leq\frac{-c_{0}n}{6}\bigg{)}\leq\sum_{\ell=1}^{3}\mathbb{P}\bigg{(}\bigg{|}\sum_{i\in{\cal
I}_{\ell}}x_{i}x_{\pi^{\natural}(i)}\bigg{|}\geq c_{1}n\bigg{)}.$ (13)

Recalling the fact that $i$ and $\pi^{\natural}(i)$ belong to different
categories, we have $\sum_{i\in{\cal I}_{\ell}}x_{i}x_{\pi^{\natural}(i)}$ to
be identically distributed as $\langle\bm{z}_{1},\bm{z}_{2}\rangle$, where
$\bm{z}_{1},\bm{z}_{2}\in\mathbb{R}^{h_{\ell}}$ are i.i.d sub-gaussian random
vectors. Then we obtain

$\displaystyle\mathbb{P}\bigg{(}\bigg{|}\sum_{i\in{\cal
I}_{\ell}}x_{i}x_{\pi^{\natural}(i)}\bigg{|}\geq c_{1}n\bigg{)}\leq~{}$
$\displaystyle\mathbb{P}\big{(}{\left\|\bm{z}_{2}\right\|}_{2}\gtrsim\sqrt{h_{\ell}\log
n}\big{)}+\mathbb{P}\big{(}|\langle\bm{z}_{1},\bm{z}_{2}\rangle|\geq
c_{1}n,~{}{\left\|\bm{z}_{2}\right\|}_{2}\lesssim\sqrt{h_{\ell}\log n}\big{)}$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\leq}}~{}$
$\displaystyle e^{-cn}+2\exp\left(-\frac{c_{1}^{2}n^{2}}{h_{\ell}\log
n}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{3}⃝}}}}{{\lesssim}}e^{-cn},\vspace{1mm}$

where \raisebox{-.8pt}{2}⃝ is due to the tail bound of sub-gaussian RVs and
\raisebox{-.8pt}{3}⃝ is because $h_{\ell}\leq h\lesssim n$. Thus we show
$\zeta_{1}\lesssim e^{-cn}$.

Analysis of $\zeta_{2}$. Its analysis is a direct consequence of Lemma 20 as
$\sum_{i=h+1}^{n}x_{i}^{2}$ can be viewed as a $\chi^{2}$ RV with freedom
$n-h$. In formulae:

$\displaystyle\zeta_{2}\leq~{}$
$\displaystyle\mathbb{P}\left(\left\lvert\sum_{i=h+1}^{n}x_{i}^{2}-(n-h)\right\rvert\geq\frac{n-h}{2}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{4}⃝}}}}{{\leq}}\exp\left[-c_{0}\left(\frac{n-h}{{|\kern-1.72218pt|\kern-1.72218pt|\mathbf{I}_{(n-h)\times(n-h)}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}}\wedge\frac{(n-h)^{2}}{{|\kern-1.72218pt|\kern-1.72218pt|\mathbf{I}_{(n-h)\times(n-h)}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}^{2}}\right)\right]$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{5}⃝}}}}{{=}}~{}$
$\displaystyle e^{-cn},$ (14)

where in \raisebox{-.8pt}{4}⃝ we use the Hanson-Wright inequality (Theorem
$6.2.1$ in Vershynin (2018)) and in \raisebox{-.8pt}{5}⃝ we use the fact
$h\lesssim n$. The proof is thus completed by combining (12), (13) and (A.3).
∎

######  Lemma 2.

We have $\mathbb{P}\left({\cal E}_{2}\right)\geq
1-c_{0}e^{-c_{1}n}-c_{2}n^{-n}$, where $c_{0},c_{1}$ and $c_{2}$ are some
positive constants.

###### Proof.

We begin the proof with the union bound, which proceeds as

$\displaystyle\mathbb{P}(\overline{{\cal
E}}_{2})\leq\mathbb{P}(\|\bm{x}\|_{2}\geq\sqrt{2n})+\mathbb{P}\left(\overline{{\cal
E}}_{2},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)\leq
c_{0}e^{-c_{1}n}+\mathbb{P}\left(\overline{{\cal
E}}_{2},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right).$ (15)

To upper-bound $\mathbb{P}\left(\overline{{\cal
E}}_{2},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)$, we first consider a fixed
permutation matrix $\bm{\Pi}_{0}$ such that
$\bm{\Pi}_{0}\neq\bm{\Pi}^{\natural}$. In addition, we define $\mathbf{M}$ as
$\bm{x}\bm{x}^{\top}(\bm{\Pi}^{\natural}-\bm{\Pi}_{0})$. Due to the
independence of the $\bm{x}$ and $\bm{w}$, we have

$\displaystyle\mathbb{P}\left(\overline{{\cal
E}}_{2},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}~{}$
$\displaystyle\mathbb{P}\left(|\bm{w}^{\top}\mathbf{M}\bm{w}-\mathbb{E}\bm{w}^{\top}\mathbf{M}\bm{w}|\geq
c\sigma^{2}n^{2}\log n,~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right),\vspace{1mm}$

where in \raisebox{-.8pt}{1}⃝ we condition on $\|\bm{x}\|_{2}\leq\sqrt{2n}$
and use the fact

$\displaystyle\mathbb{E}\bm{w}^{\top}\mathbf{M}\bm{w}+c\sigma^{2}n^{2}\log n=\
\sigma^{2}\operatorname{Tr}(\mathbf{M})+c\sigma^{2}n^{2}\log
n\lesssim\sigma^{2}\|\bm{x}\|^{2}_{2}+c\sigma^{2}n^{2}\log
n\lesssim\sigma^{2}n^{2}\log n.$

Using Hanson-Wright inequality (Theorem $6.2.1$ in Vershynin (2018)), we
obtain

$\displaystyle\mathbb{P}\left(|\bm{w}^{\top}\mathbf{M}\bm{w}-\mathbb{E}\bm{w}^{\top}\mathbf{M}\bm{w}|\geq
c\sigma^{2}n^{2}\log n,~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)$
$\displaystyle\leq~{}$ $\displaystyle
2\mathbb{E}\left\\{\exp\left[-\left(\frac{c_{0}n^{4}\log^{2}n}{{|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}^{2}}\wedge\frac{c_{1}n^{2}\log
n}{{|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}}\right)\right]\mathbbm{1}(\|\bm{x}\|_{2}\leq\sqrt{2n})\right\\}.\vspace{1mm}$

Since $\mathbf{M}$ is a rank-$1$ matrix, we have
${|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{F}}={|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}$.
Conditioning on the event $\|\bm{x}\|_{2}\leq\sqrt{2n}$, we have
${|\kern-1.72218pt|\kern-1.72218pt|\mathbf{M}|\kern-1.72218pt|\kern-1.72218pt|}_{\textup{OP}}\lesssim\|\bm{x}\|^{2}_{2}\lesssim
n$ and hence
$\mathbb{P}\left(\bm{w}^{\top}\mathbf{M}\bm{w}\gtrsim\sigma^{2}n^{2}\log
n\right)\lesssim n^{-3n}$ for a fixed $\bm{\Pi}_{0}$. Iterating over all
possible $\bm{\Pi}\neq\bm{\Pi}^{\natural}$, we use (15) and complete the proof
as

$\displaystyle\mathbb{P}(\overline{{\cal E}}_{2})\leq c_{0}e^{-c_{1}n}+|{\cal
P}_{n}|\cdot\mathbb{P}(\overline{{\cal
E}}_{2},\|\bm{x}\|_{2}\leq\sqrt{2n})\leq c_{0}e^{-c_{1}n}+c_{2}\cdot n!\cdot
n^{-3n}\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\lesssim}}e^{-cn}+n^{-n},$
(16)

where we use the Stirling’s approximation, namely, $n!\sim\sqrt{2\pi
n}(n/e)^{n}$, in \raisebox{-.8pt}{2}⃝. ∎

######  Lemma 3.

We have $\mathbb{P}\left({\cal E}_{3}\right)\geq
1-c_{0}e^{-c_{1}n}-c_{2}n^{-n}$, where $c_{0},c_{1}$, and $c_{2}$ are some
positive constants.

###### Proof.

The proof follows a similar strategy as in that of Lemma 2. Defining
$\bm{u}_{\bm{\Pi}}\in\mathbb{R}^{n}$ as

$\displaystyle\bm{u}_{\bm{\Pi}}\triangleq\langle\bm{\Pi}^{\natural}\bm{x},(\bm{\Pi}^{\natural}-\bm{\Pi})^{\top}\bm{x}\rangle\cdot\bm{x}+\langle\bm{\Pi}^{\natural}\bm{x},\bm{x}\rangle(\bm{\Pi}^{\natural}-\bm{\Pi})^{\top}\bm{x},\vspace{1mm}$

we first rewrite the $\mathbb{P}\left(\overline{{\cal E}}_{3}\right)$ as

$\displaystyle\mathbb{P}\left(\overline{{\cal
E}}_{3}\right)=\mathbb{P}\left(|\left\langle\bm{u}_{\bm{\Pi}},\bm{w}\right\rangle|\gtrsim\sigma
n^{2}\sqrt{\log
n},~{}\exists~{}\bm{\Pi}\neq\bm{\Pi}^{\natural}\right).\vspace{1mm}$

With the union bound, we can obtain

$\displaystyle\mathbb{P}\left(\overline{{\cal E}}_{3}\right)\leq
e^{-cn}+\mathbb{P}\left(\overline{{\cal
E}}_{3},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)\leq e^{-cn}+|{\cal
P}_{n}|\cdot\mathbb{P}\left(|\left\langle\bm{u}_{\bm{\Pi}},\bm{w}\right\rangle|\gtrsim\sigma
n^{2}\sqrt{\log n},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right).\vspace{1mm}$

Conditioning on $\\{\|\bm{x}\|_{2}\leq\sqrt{2n}\\}$, we have the relation
$\|\bm{u}_{\bm{\Pi}}\|_{2}^{2}\lesssim\sigma^{2}(4{\left\|\bm{x}\right\|}_{2}^{3})^{2}=c\sigma^{2}n^{3}$
and hence

$\displaystyle\mathbb{P}\left(|\left\langle\bm{u}_{\bm{\Pi}},\bm{w}\right\rangle|\gtrsim\sigma
n^{2}\sqrt{\log n},~{}\|\bm{x}\|_{2}\leq\sqrt{2n}\right)\leq~{}$
$\displaystyle
2\cdot\mathbb{E}_{\bm{u}_{\bm{\Pi}}}\left[\exp\left(-\frac{c\sigma^{2}n^{4}\log
n}{2\|\bm{u}_{\bm{\Pi}}\|^{2}_{2}}\right)\mathbbm{1}(\|\bm{u}_{\bm{\Pi}}\|^{2}_{2}\lesssim\sigma^{2}n^{3})\right]$
$\displaystyle\leq~{}$ $\displaystyle 2\exp\left(-\frac{c\sigma^{2}n^{4}\log
n}{2\sigma^{2}n^{3}}\right)=2\cdot n^{-3n}.\vspace{1mm}$

Following the same logic as (16), we complete the proof.

∎

######  Lemma 4.

We have $\mathbb{P}\left({\cal E}_{4}\right)\geq 1-c_{0}\cdot n^{-c_{1}}$,
where $c_{0},c_{1}>0$ are some positive constants.

###### Proof.

To begin with, we notice the relation

$\displaystyle\mathbb{P}(\overline{{\cal
E}}_{4})\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}~{}$
$\displaystyle\mathbb{P}\left(\min_{s\neq
t}|\bm{x}_{s}-\bm{x}_{t}|^{2}\lesssim n^{-20}\right)\leq
n^{2}\cdot\mathbb{P}\left(|\bm{x}_{s}-\bm{x}_{t}|\lesssim n^{-10}\right)$
$\displaystyle\leq~{}$ $\displaystyle
n^{2}\cdot\left[\mathbb{P}(|\bm{x}_{s}|\lesssim
n^{-10})+\mathbb{P}(|\bm{x}_{t}|\lesssim
n^{-10})\right]=2n^{2}\cdot\mathbb{P}(|\bm{x}_{s}|\lesssim n^{-20}),$ (17)

where \raisebox{-.8pt}{1}⃝ is due to the relation
${\left\|\bm{x}-\bm{\Pi}\bm{x}\right\|}_{2}^{2}\geq\min_{s\neq
t}|\bm{x}_{s}-\bm{x}_{t}|^{2}$ as long as $\bm{\Pi}\neq\mathbf{I}$. Invoking
(23), we have

$\displaystyle\mathbb{P}(|\bm{x}_{s}|\lesssim n^{-20})\lesssim\exp\left(-c\log
n\right)=n^{-c},\vspace{1mm}$

and complete the proof when combining with (A.3).

∎

##  Appendix B Appendix for Section 4

This section provides a theoretical analysis for the multiple observations
model, i.e., $m>1$. Without specification, only the centered sub-gaussian
assumption is put on $\mathbf{X}_{ij}$ and the log-concavity is not assumed.

###  B.1 Notations: multiple observations model

We begin the proof by defining the notations. First, we define
$\widetilde{\mathbf{B}}$ and $\widehat{\mathbf{B}}$ as

$\displaystyle\widetilde{\mathbf{B}}\triangleq~{}$
$\displaystyle(n-h)^{-1}\mathbf{X}^{\top}\bm{\Pi}^{\natural}\mathbf{X}\mathbf{B}^{\natural};$
$\displaystyle\widehat{\mathbf{B}}\triangleq~{}$
$\displaystyle(n-h)^{-1}\mathbf{X}^{\top}\mathbf{Y}=\widetilde{\mathbf{B}}+(n-h)^{-1}\mathbf{X}^{\top}\mathbf{W}.\vspace{1mm}$

In addition, we define the following events ${\cal F}_{\ell}$ $(1\leq\ell\leq
8)$ as

$\displaystyle{\cal F}_{1}(\mathbf{M})$
$\displaystyle\triangleq\left\\{{\left\|\mathbf{M}^{\top}\mathbf{X}_{s,:}\right\|}_{2}\lesssim\sqrt{\log
n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{M}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\textup{
and
}\big{\|}\mathbf{M}^{\top}\mathbf{X}^{{}^{\prime}}_{s,:}\big{\|}_{2}\lesssim\sqrt{\log
n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{M}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}~{}~{}\forall~{}1\leq
s\leq n\right\\};$ $\displaystyle{\cal F}_{2,1}$
$\displaystyle\triangleq\left\\{\langle\mathbf{X}_{s,:},\mathbf{X}_{t,:}^{{}^{\prime}}\rangle\lesssim(\log
n)\sqrt{p},~{}~{}1\leq s,t\leq n\right\\};$ $\displaystyle{\cal F}_{2,2}$
$\displaystyle\triangleq\left\\{\langle\mathbf{X}_{s,:},\mathbf{X}_{t,:}\rangle\lesssim(\log
n)\sqrt{p},~{}~{}1\leq s\neq t\leq n\right\\};$ $\displaystyle{\cal F}_{2,3}$
$\displaystyle\triangleq\left\\{\langle\mathbf{X}^{{}^{\prime}}_{s,:},\mathbf{X}^{{}^{\prime}}_{t,:}\rangle\lesssim(\log
n)\sqrt{p},~{}~{}1\leq s\neq t\leq n\right\\};$ $\displaystyle{\cal F}_{2}$
$\displaystyle={\cal F}_{2,1}\bigcap{\cal F}_{2,2}\bigcap{\cal F}_{2,3};$
$\displaystyle{\cal F}_{3}$
$\displaystyle=\left\\{{\left\|\mathbf{X}_{s,:}\right\|}_{2}\leq\sqrt{p\log
n}\textup{ and }\|\mathbf{X}^{{}^{\prime}}_{s,:}\|_{2}\leq\sqrt{p\log
n},~{}~{}\forall~{}1\leq s\leq n\right\\};$ $\displaystyle{\cal F}_{4}$
$\displaystyle=\left\\{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{X}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\leq\sqrt{2np}~{}\textup{
and
}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{X}_{\setminus(s)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\leq\sqrt{2np},~{}~{}\forall~{}1\leq
s\leq n\right\\};$ $\displaystyle{\cal F}_{5}$
$\displaystyle=\left\\{{\left\|\mathbf{X}\mathbf{X}_{s,:}\right\|}_{2}\lesssim(\log
n)\sqrt{np},~{}~{}\forall~{}1\leq s\leq n\right\\};$ $\displaystyle{\cal
F}_{6,\textup{original}}$
$\displaystyle=\left\\{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}-\widetilde{\mathbf{B}}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\lesssim\frac{(\log
n)(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right\\};$
$\displaystyle{\cal F}_{6,\textup{single}}$
$\displaystyle=\left\\{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}-\widetilde{\mathbf{B}}_{\setminus(s)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\lesssim\frac{(\log
n)(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}~{}\forall~{}1\leq
s\leq n\right\\};$ $\displaystyle{\cal F}_{6,\textup{double}}$
$\displaystyle=\left\\{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}-\widetilde{\mathbf{B}}_{\setminus(s,t)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\lesssim\frac{(\log
n)(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}~{}\forall~{}1\leq
s\neq t\leq n\right\\};$ $\displaystyle{\cal F}_{6}$ $\displaystyle={\cal
F}_{6,\textup{original}}\bigcap{\cal F}_{6,\textup{single}}\bigcap{\cal
F}_{6,\textup{double}};$ $\displaystyle{\cal F}_{7,\textup{single}}$
$\displaystyle=\
\left\\{{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(s)})^{\top}\mathbf{X}_{s,:}\right\|}_{2}\lesssim\frac{p\log^{\nicefrac{{3}}{{2}}}n}{n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}\forall~{}1\leq
s\leq n\right\\};$ $\displaystyle{\cal F}_{7,\textup{double}}$
$\displaystyle=\left\\{{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(s,t)})^{\top}\mathbf{X}_{s,:}\right\|}_{2}\lesssim\frac{p\log^{\nicefrac{{3}}{{2}}}n}{n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}\forall~{}1\leq
s\neq t\leq n\right\\};$ $\displaystyle{\cal F}_{7}$ $\displaystyle={\cal
F}_{7,\textup{single}}\bigcap{\cal F}_{7,\textup{double}};$
$\displaystyle{\cal F}_{8}$
$\displaystyle=\left\\{\big{\|}(\widetilde{\mathbf{B}}-\mathbf{B}^{\natural})^{\top}\mathbf{X}_{s,:}\big{\|}_{2}\lesssim\frac{(\log
n)^{\nicefrac{{3}}{{2}}}(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}\forall~{}1\leq
s\leq n\right\\},\vspace{1mm}$

where $\mathbf{M}$ in ${\cal F}_{1}(\mathbf{M})$ is an arbitrary matrix
independent of $\mathbf{X}_{i,:}$. In addition, we define the quantities
$\Delta_{1}$, $\Delta_{2}$, and $\Delta_{3}$ as

$\displaystyle\Delta_{1}$
$\displaystyle=c_{0}\sigma(\log^{\nicefrac{{5}}{{2}}}n)\sqrt{\frac{p}{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}};$
(18) $\displaystyle\Delta_{2}$
$\displaystyle=c_{1}\sigma(\log^{2}n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}};$
(19) $\displaystyle\Delta_{3}$
$\displaystyle=c_{2}\left(\frac{mp\sigma^{2}(\log^{2}n)}{n}+\sigma^{2}(\log^{2}n)\sqrt{\frac{mp}{n}}\right),$
(20)

respectively. Besides, we denote $\Delta_{(\textup{multi, 1})}$ as
$\Delta_{1}+\Delta_{2}+\Delta_{3}$ and $\Delta_{(\textup{multi, 2})}$ as
$\Delta_{2}+\Delta_{3}$.

The following context presents the analysis corresponding to each regime,
which is organized in ascending order of difficulties. To facilitate
understanding, we put a diagram illustrating the dependence among lemmas in
Figure 12.

Easy Regimesub-gaussian
$\&$$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$Medium Regimesub-
gaussian $\&$$\textup{srank}(\mathbf{B}^{\natural})\gg\log n$Hard Regimelog-
concave sub-gaussian$\&$ $\textup{srank}(\mathbf{B}^{\natural})\geq c$Lemma
5Lemma 6Lemma 7Lemma 8Lemma 9Lemma 10Lemma 11Lemma 12, 13, 14, 15, 16, 17, 18,
and 19 Figure 12: Dependence diagram of lemmas.

Before proceeding to the technical details, we first restate our construction
method of $\mathbf{X}_{\setminus(s)}$,
$\widetilde{\mathbf{B}}_{\setminus(s)}$, and
$\widetilde{\mathbf{B}}_{\setminus(s,t)}$, respectively. To begin with, we
construct the leave-one-out sensing matrices. For each row $\mathbf{X}_{s,:}$
of $\mathbf{X}$ $(1\leq s\leq n)$, we draw an independent copy denoted as
$\mathbf{X}_{s,:}^{{}^{\prime}}$. Matrix $\mathbf{X}_{\setminus(s)}$ is then
constructed by $(i)$ replacing the $s$th row in the sensing matrix
$\mathbf{X}$ with $\mathbf{X}_{s,:}^{{}^{\prime}}$ and $(ii)$ copying the rest
rows.

In addition, we construct leave-one-out samples
$\widetilde{\mathbf{B}}_{\setminus(s)}$ ($1\leq s\leq n$). This is done by
replacing the role of $\mathbf{X}_{s,:}$ in $\widetilde{\mathbf{B}}$ with its
independent copy $\mathbf{X}_{s,:}^{{}^{\prime}}$. In formulae:

$\displaystyle\widetilde{\mathbf{B}}_{\setminus(s)}=(n-h)^{-1}\bigg{(}\sum_{\begin{subarray}{c}k\neq
s\\\ \pi^{\natural}(k)\neq
s\end{subarray}}\mathbf{X}_{\pi^{\natural}(k),:}\mathbf{X}_{k,:}^{\top}+\sum_{\begin{subarray}{c}k=s\textup{
or}\\\
\pi^{\natural}(k)=s\end{subarray}}{\mathbf{X}}^{{}^{\prime}}_{\pi^{\natural}(k),:}{\mathbf{X}}^{{}^{\prime}\top}_{k,:}\bigg{)}\mathbf{B}^{\natural}.\vspace{1mm}$

Easily, we can verify that $\widetilde{\mathbf{B}}_{\setminus(s)}$ is
independent of $\mathbf{X}_{s,:}$. Similarly, we construct matrices
$\\{\widetilde{\mathbf{B}}_{\setminus(s,t)}\\}_{1\leq s\neq t\leq n}$ as

$\displaystyle\widetilde{\mathbf{B}}_{\setminus(s,t)}=(n-h)^{-1}\bigg{(}\sum_{\begin{subarray}{c}k\neq
s,t\\\ \pi^{\natural}(k)\neq
s,t\end{subarray}}\mathbf{X}_{\pi^{\natural}(k),:}\mathbf{X}_{k,:}^{\top}+\sum_{\begin{subarray}{c}k=s\textup{
or }k=t\textup{ or}\\\ \pi^{\natural}(k)=s\textup{ or
}\pi^{\natural}(k)=t\end{subarray}}{\mathbf{X}}^{{}^{\prime}}_{\pi^{\natural}(k),:}{\mathbf{X}}^{{}^{\prime}\top}_{k,:}\bigg{)}\mathbf{B}^{\natural}.\vspace{1mm}$

Moreover, we can show that $\widetilde{\mathbf{B}}_{\setminus(s,t)}$ is
independent from the rows $\mathbf{X}_{s,:}$ and $\mathbf{X}_{t,:}$ $(1\leq
s\neq t\leq n)$.

###  B.2 Easy regime: proof of Theorem 4

###### Proof.

To begin with, we notice that the construction error occurs when there exists
some permutation matrix $\bm{\Pi}\neq\bm{\Pi}^{\natural}$ such that
$\langle\mathbf{Y},\bm{\Pi}\mathbf{X}\widehat{\mathbf{B}}\rangle\geq\langle\mathbf{Y},\bm{\Pi}^{\natural}\mathbf{X}\widehat{\mathbf{B}}\rangle$.

Step I. Defining event ${\cal F}_{\textup{err-relax}}^{(\textup{multi},1)}$ as

$\displaystyle{\cal F}_{\textup{err-
relax}}^{(\textup{multi},1)}\triangleq\left\\{\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\rangle\leq\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{j,:}\right\rangle,~{}\exists~{}1\leq\pi^{\natural}(i)\neq
j\leq n\right\\},\vspace{1mm}$

we first show
$\\{\bm{\Pi}^{\textup{opt}}\neq\bm{\Pi}^{\natural}\\}\subseteq{\cal
F}_{\textup{err-relax}}^{(\textup{multi},1)}$. The detailed reasoning is as
follows. Conditional on $\overline{{\cal F}}_{\textup{err-
relax}}^{(\textup{multi},1)}$, we have

$\displaystyle\langle\mathbf{Y},\bm{\Pi}^{\natural}\mathbf{X}\widehat{\mathbf{B}}\rangle=\sum_{i}\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\rangle\geq\sum_{i}\left\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi(i),:}\right\rangle=\langle\mathbf{Y},\bm{\Pi}\mathbf{X}\widehat{\mathbf{B}}\rangle,\vspace{1mm}$

which means ground truth $\bm{\Pi}^{\natural}$ will be returned by our
estimator in Algorithm 1, in other words, $\overline{{\cal F}}_{\textup{err-
relax}}^{(\textup{multi},1)}\subseteq\\{\bm{\Pi}^{\textup{opt}}=\bm{\Pi}^{\natural}\\}$.

Step II. Then we will upper bound the error probability with $\mathbb{P}({\cal
F}_{\textup{err-relax}}^{(\textup{multi},1)})$. With the relation
$\mathbf{Y}_{i,:}=\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}+\mathbf{W}_{i,:}$
and
$\widehat{\mathbf{B}}=\widetilde{\mathbf{B}}+(n-h)^{-1}\mathbf{X}^{\top}\mathbf{W}$,
we obtain an equivalent form of condition
$\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\rangle\leq\langle\mathbf{Y}_{i,:},\widehat{\mathbf{B}}^{\top}\mathbf{X}_{j,:}\rangle$
reading as

$\displaystyle\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}+\mathbf{W}_{i,:},\left(\widetilde{\mathbf{B}}+(n-h)^{-1}\mathbf{X}^{\top}\mathbf{W}\right)^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\rangle$
$\displaystyle\leq~{}$
$\displaystyle\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}+\mathbf{W}_{i,:},\left(\widetilde{\mathbf{B}}+(n-h)^{-1}\mathbf{X}^{\top}\mathbf{W}\right)^{\top}\mathbf{X}_{j,:}\right\rangle.$
(21)

For the notation conciseness, we define terms $\mathsf{Term}_{i}$ ($1\leq
i\leq 4$) as

$\displaystyle\mathsf{Term}_{\textup{tot}}$
$\displaystyle=\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{\pi^{\natural}(i),:}-\mathbf{X}_{j,:}\right)\right\rangle;$
(22) $\displaystyle\mathsf{Term}_{1}$
$\displaystyle=(n-h)^{-1}\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\mathbf{W}^{\top}\mathbf{X}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\rangle;$
(23) $\displaystyle\mathsf{Term}_{2}$
$\displaystyle=\left\langle\mathbf{W}_{i,:},\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\rangle;$
(24) $\displaystyle\mathsf{Term}_{3}$
$\displaystyle=(n-h)^{-1}\left\langle\mathbf{W}\mathbf{W}_{i,:},\mathbf{X}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\rangle.$
(25)

Then (B.2) can be rewritten in a concise form, i.e.,
$\mathsf{Term}_{\textup{tot}}\leq\mathsf{Term}_{1}+\mathsf{Term}_{2}+\mathsf{Term}_{3}$.
With the union bound, we obtain

$\displaystyle\mathbb{P}(\pi^{\natural}(i)\neq\widehat{\pi}(i),\exists~{}i)=~{}$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{\textup{tot}}\leq\mathsf{Term}_{1}+\mathsf{Term}_{2}+\mathsf{Term}_{3},~{}\exists~{}i,j\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]+\sum_{a=1}^{8}\mathbb{P}\left(\overline{{\cal
F}}_{a}\right)$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}~{}$
$\displaystyle
n^{2}\cdot\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{\textup{tot}}\leq\mathsf{Term}_{1}+\mathsf{Term}_{2}+\mathsf{Term}_{3}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]+c_{0}p^{-c_{1}}+c_{2}n^{-c_{3}},$ (26)

where in \raisebox{-.8pt}{1}⃝ we invoke Lemma 12, Lemma 13, Lemma 14, Lemma
15, Lemma 16, Lemma 17, Lemma 18, and Lemma 19.

Regarding the term
$\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{\textup{tot}}\leq\mathsf{Term}_{1}+\mathsf{Term}_{2}+\mathsf{Term}_{3},~{}\exists~{}i,j\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]$, we perform the following decomposition

$\displaystyle\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{\textup{tot}}\leq\mathsf{Term}_{1}+\mathsf{Term}_{2}+\mathsf{Term}_{3}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]$ $\displaystyle\leq~{}$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{\textup{tot}}\leq\Delta_{(\textup{multi,
1})}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]+\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{1}\geq\Delta_{1}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]$ $\displaystyle+~{}$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{2}\geq\Delta_{2}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right]+\mathbb{E}\left[\mathbbm{1}\left(\mathsf{Term}_{3}\geq\Delta_{3}\right)\mathbbm{1}\left(\bigcap_{a=1}^{8}{\cal
F}_{a}\right)\right],$ (27)

where the definitions of $\Delta_{1}$, $\Delta_{2}$, $\Delta_{3}$, and
$\Delta_{(\textup{multi, 1})}$ are referred to Subsection B.1. The proof is
then completed by combining (B.2) and (B.2) and invoking Lemma 5, Lemma 6,
Lemma 7, and Lemma 8.

∎

######  Lemma 5.

Consider fixed indices $\pi^{\natural}(i)$ and $j$ such that
$j\neq\pi^{\natural}(i)$. Assuming that $(i)$
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$, $(ii)$ $n\gtrsim
p\log^{6}n$, $(iii)$ intersection of events ${\cal
F}_{1}(\mathbf{B}^{\natural})\bigcap{{\cal
F}_{1}\big{(}\mathbf{B}^{\natural}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\big{)}}\bigcap{\cal
F}_{6}\bigcap{\cal F}_{7}\bigcap{\cal F}_{8}$ holds, and $(iv)$
$\mathsf{SNR}\geq c$, we have
$\\{\mathsf{Term}_{\textup{tot}}\geq\Delta_{(\textup{multi, 1})}\\}$ hold with
probability $1-n^{-c}$ when $n$ and $p$ are sufficiently large. Here
$\mathsf{Term}_{\textup{tot}}$ and $\Delta_{(\textup{multi, 1})}$ are defined
in (22) and Section B.1, respectively.

###### Proof.

We begin the proof with the decomposition

$\displaystyle\mathsf{Term}_{\textup{tot}}={\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}^{2}+\underbrace{\big{\langle}\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\left(\widetilde{\mathbf{B}}-\mathbf{B}^{\natural}\right)^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\big{\rangle}}_{\triangleq\mathsf{Term}_{\textup{tot},1}}-\underbrace{\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\widetilde{\mathbf{B}}^{\top}\mathbf{X}_{j,:}\right\rangle}_{\triangleq\mathsf{Term}_{\textup{tot},2}}.\vspace{1mm}$

Then we obtain

$\displaystyle\mathbb{P}\left(\mathsf{Term}_{\textup{tot}}\leq\Delta_{(\textup{multi,
1})}\right)=~{}$ $\displaystyle\mathbb{P}\left(\frac{\Delta_{(\textup{multi,
1})}}{{\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}^{2}}-\frac{\mathsf{Term}_{\textup{tot},1}}{{\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}^{2}}+\frac{\mathsf{Term}_{\textup{tot},2}}{{\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}^{2}}\geq
1\right)$ $\displaystyle\leq~{}$
$\displaystyle\underbrace{\mathbb{P}\left({\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\leq\delta\right)}_{\triangleq~{}\zeta_{1}}+\underbrace{\mathbb{P}\left(\frac{\Delta_{(\textup{multi,
1})}}{\delta^{2}}+\frac{\left\lvert\mathsf{Term}_{\textup{tot},1}\right\rvert}{\delta^{2}}+\frac{\left\lvert\mathsf{Term}_{\textup{tot},2}\right\rvert}{\delta^{2}}\geq
1\right)}_{\triangleq~{}\zeta_{2}}.$ (28)

Setting $\delta$ as
$1/2{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$,
we separately bound the probabilities $\zeta_{1}$ and $\zeta_{2}$. For the
term $\zeta_{1}$, we invoke the small ball probability (Lemma 22) and conclude

$\displaystyle\mathbb{P}\left({\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\leq\frac{1}{2}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)\leq
e^{-c_{0}\cdot\textup{srank}(\mathbf{B}^{\natural})}.$ (29)

For probability $\zeta_{2}$, we will show it is upper-bounded by $c_{1}\cdot
n^{-c_{2}}$ provided $\mathsf{SNR}\geq C$. A detailed explanation comes as
follows.

Phase I. First, we consider $\mathsf{Term}_{\textup{tot},1}$. Conditional on
the intersection of events ${\cal F}_{1}(\mathbf{B}^{\natural})\bigcap{\cal
F}_{8}$, we have

$\displaystyle\left\lvert\mathsf{Term}_{\textup{tot},1}\right\rvert\leq~{}$
$\displaystyle{\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{i,:}\right\|}_{2}\cdot{\left\|(\widetilde{\mathbf{B}}-\mathbf{B}^{\natural})^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\lesssim\sqrt{\log
n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\frac{(\log
n)^{\nicefrac{{3}}{{2}}}(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$
$\displaystyle=~{}$ $\displaystyle(\log^{2}n)(\log
n^{2}p^{3})\sqrt{\frac{p}{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2}.\vspace{1mm}$

Phase II. Then we turn to $\mathsf{Term}_{\textup{tot},2}$. Adopting the
leave-out-out trick, we can expand it as

$\displaystyle\mathsf{Term}_{\textup{tot},2}=~{}\underbrace{\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{j,:}\right\rangle}_{\mathsf{Term}_{\textup{tot},2,1}}+\underbrace{\left\langle\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\mathbf{X}_{j,:}\right\rangle}_{\mathsf{Term}_{\textup{tot},2,2}}.\vspace{1mm}$

For $\mathsf{Term}_{\textup{tot},2,1}$, we have

$\displaystyle\mathsf{Term}_{\textup{tot},2,1}\leq~{}$
$\displaystyle{\left\|\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\big{\|}(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{j,:}\big{\|}_{2}\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\lesssim}}\sqrt{\log
n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\frac{p\log^{\nicefrac{{3}}{{2}}}n}{n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$
$\displaystyle=~{}$
$\displaystyle\frac{p\log^{2}n}{n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2},\vspace{1mm}$

where in \raisebox{-.8pt}{1}⃝ we condition on event ${\cal
F}_{7,\textup{double}}$. Regarding $\mathsf{Term}_{2,2,2}$, we notice that
$\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}$ is independent of
the rows $\mathbf{X}_{\pi^{\natural}(i),:}$ and $\mathbf{X}_{j,:}$ due to its
construction method. Hence we can bound $\mathsf{Term}_{2,2,2}$ by
conditioning on rows $\\{\mathbf{X}_{s,:}\\}_{s\neq\pi^{\natural}}$ and only
viewing $\mathbf{X}_{\pi^{\natural}(i),:}$ as the RV, which yields

$\displaystyle\mathsf{Term}_{\textup{tot},2,2}\lesssim\sqrt{\log
n}{\left\|\mathbf{B}^{\natural}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\mathbf{X}_{j,:}\right\|}_{2},$
(30)

holds with probability $1-n^{-c}$. Conditional on event ${\cal
F}_{1}(\mathbf{B}^{\natural}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top})$,
we have

$\displaystyle\mathsf{Term}_{\textup{tot},2,2}\lesssim~{}$ $\displaystyle\log
n{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\lesssim\log
n{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{OP}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\leq}}~{}$
$\displaystyle(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{OP}}\left[{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}-\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}+{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right]\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{3}⃝}}}}{{\lesssim}}\frac{(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2}}{\sqrt{\textup{srank}(\mathbf{B}^{\natural})}},\vspace{1mm}$

where \raisebox{-.8pt}{2}⃝ uses the definition of stable rank, and
\raisebox{-.8pt}{3}⃝ is conditional on event ${\cal F}_{6}$, $n\geq p$, and
$n\gtrsim p\log^{6}n$.

Phase III. Conditioning on (30), we expand the sum
$\nicefrac{{\Delta_{(\textup{multi,
1})}}}{{\delta^{2}}}+\mathsf{Term}_{\textup{tot},1}/\delta^{2}+\mathsf{Term}_{\textup{tot},2}/\delta^{2}$
as

$\displaystyle\frac{\Delta_{(\textup{multi,
1})}}{\delta^{2}}+\frac{\mathsf{Term}_{\textup{tot},1}}{\delta^{2}}+\frac{\mathsf{Term}_{\textup{tot},2}}{\delta^{2}}$
$\displaystyle=~{}$ $\displaystyle c_{0}\sigma(\log
n)^{\nicefrac{{5}}{{2}}}\sqrt{\frac{p}{n}}\frac{1}{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}}+\frac{c_{1}\sigma(\log^{2}n)}{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}}+c_{2}\left(\frac{pm}{n}+\sqrt{\frac{mp}{n}}\right)\frac{\sigma^{2}(\log^{2}n)}{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2}}$
$\displaystyle+~{}$ $\displaystyle\frac{c_{3}(\log^{2}n)(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}+\frac{c_{4}p\log^{2}n}{n}+\frac{c_{5}\log
n}{\sqrt{\textup{srank}(\mathbf{B}^{\natural})}}$ $\displaystyle\asymp~{}$
$\displaystyle c_{0}\sqrt{\frac{p}{nm}}\frac{\left(\log
n\right)^{\nicefrac{{5}}{{2}}}}{\sqrt{\mathsf{SNR}}}+\frac{c_{1}\log^{2}n}{\sqrt{m\cdot\mathsf{SNR}}}+\frac{c_{2}p(\log^{2}n)}{n\cdot\mathsf{SNR}}+c_{2}\sqrt{\frac{p}{mn}}\frac{\log^{2}n}{\mathsf{SNR}}$
$\displaystyle+~{}$ $\displaystyle\frac{c_{3}(\log^{2}n)(\log
n^{2}p^{3})\sqrt{p}}{\sqrt{n}}+\frac{c_{4}p\log^{2}n}{n}+\frac{c_{5}\log
n}{\sqrt{\textup{srank}(\mathbf{B}^{\natural})}}.\vspace{1mm}$

Provided that $\mathsf{SNR}\geq c$,
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$ and $n\gtrsim p\log^{6}n$,
we can verify the sum $\nicefrac{{\Delta_{(\textup{multi,
1})}}}{{\delta^{2}}}+\mathsf{Term}_{\textup{tot},1}/\delta^{2}+\mathsf{Term}_{\textup{tot},2}/\delta^{2}\ll
1$ when $n$ and $p$ are sufficiently large. Hence we can conclude

$\displaystyle\zeta_{2}\leq\mathbb{P}\left(\mathsf{Term}_{\textup{tot},2,2}\gtrsim\sqrt{\log
n}{\left\|\mathbf{B}^{\natural}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\mathbf{X}_{j,:}\right\|}_{2}\right)\leq
n^{-c}.\vspace{1mm}$

Hence the proof is completed by combining (B.2) and (29). ∎

######  Remark 2.

If we strength the requirement on $\mathsf{SNR}$ from $\mathsf{SNR}\geq c$ to
$\mathsf{SNR}\gtrsim\log^{2}n$, we can relax the requirement on the stable
rank $\textup{srank}(\mathbf{B}^{\natural})$ from
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{4}n$ to
$\textup{srank}(\mathbf{B}^{\natural})\gg\log^{2}n$.

######  Lemma 6.

Conditional on the intersection of events ${\cal
F}_{1}(\mathbf{B}^{\natural})\bigcap{\cal F}_{5}$ and fixing the indices
$\pi^{\natural}(i)$ and $j$ $(j\neq\pi^{\natural}(i))$, we have

$\displaystyle\mathsf{Term}_{1}\lesssim\sigma(\log
n)^{\nicefrac{{5}}{{2}}}\sqrt{\frac{p}{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}.\vspace{1mm}$

hold with probability at least $1-n^{-c}$, where $\mathsf{Term}_{1}$ is
defined in (23).

###### Proof.

Define vectors $\bm{u}_{\mathbf{X}}$ and $\bm{v}_{\mathbf{X}}^{\top}$ as

$\displaystyle\bm{u}_{\mathbf{X}}$
$\displaystyle=\mathbf{X}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right),$
$\displaystyle\bm{v}_{\mathbf{X}}$
$\displaystyle=\mathbf{B}^{\natural\top}\mathbf{X}_{\pi^{\natural}(i),:},\vspace{1mm}$

respectively. We can rewrite $\mathsf{Term}_{1}$ as

$\displaystyle\mathsf{Term}_{1}=~{}$
$\displaystyle(n-h)^{-1}\operatorname{Tr}\left[\mathbf{X}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\mathbf{X}^{\top}_{\pi^{\natural}(i),:}\mathbf{B}^{\natural}\mathbf{W}^{\top}\right]=(n-h)^{-1}\bm{u}_{\mathbf{X}}^{\top}\mathbf{W}\bm{v}_{\mathbf{X}}.\vspace{1mm}$

Invoking the union bound, we conclude

$\displaystyle\mathbb{P}\left(\mathsf{Term}_{1}\gtrsim\sigma(\log
n)^{\nicefrac{{5}}{{2}}}\sqrt{\frac{p}{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle\leq~{}$
$\displaystyle\mathbb{P}\left(\mathsf{Term}_{1}\gtrsim\sigma(\log
n)^{\nicefrac{{5}}{{2}}}\sqrt{\frac{p}{n}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}{\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}\lesssim(\log
n)^{\nicefrac{{3}}{{2}}}\sqrt{np}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle+~{}$
$\displaystyle\mathbb{P}\left({\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}\gtrsim(\log
n)^{\nicefrac{{3}}{{2}}}\sqrt{np}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle\leq~{}$
$\displaystyle\underbrace{\mathbb{P}\left(\mathsf{Term}_{1}\gtrsim\frac{\sigma(\log
n){\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}}{n-h}\right)}_{\triangleq~{}\zeta_{1}}+\underbrace{\mathbb{P}\left({\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}\gtrsim(\log
n)^{\nicefrac{{3}}{{2}}}\sqrt{np}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)}_{\triangleq~{}\zeta_{2}}.$
(31)

Then we separately bound the probabilities $\zeta_{1}$ and $\zeta_{2}$.

Phase I. For probability $\zeta_{1}$, we exploit the independence between
$\mathbf{X}$ and $\mathbf{W}$. We can view $\mathsf{Term}_{1}$ as a Gaussian
RV conditional on $\mathbf{X}$, since it is a linear combination of Gaussian
RVs $\left\\{\mathbf{W}_{i,j}\right\\}_{1\leq i\leq n,1\leq j\leq m}$. Easily
we can calculate its mean as zero and its variance as

$\displaystyle\mathbb{E}_{\mathbf{W}}(\mathsf{Term}_{1})^{2}=\frac{\sigma^{2}}{(n-h)^{2}}{\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}^{2}.\vspace{1mm}$

Thus we can upper-bound $\zeta_{1}$ as

$\displaystyle\zeta_{1}=\mathbb{E}_{\mathbf{X}}\mathbb{E}_{\mathbf{W}}\mathbbm{1}\left(\mathsf{Term}_{1}\gtrsim\frac{\sigma(\log
n){\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}}{n-h}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}\mathbb{E}_{\mathbf{X}}\exp\left(-c_{0}\log
n\right)=n^{-c},$ (32)

where \raisebox{-.8pt}{1}⃝ is due to the bound on the tail-probability of
Gaussian RV.

Phase II. As for $\zeta_{2}$, easily we can verify it to be zero conditional
on the intersection of events ${\cal F}_{1}(\mathbf{B}^{\natural})\bigcap{\cal
F}_{5}$ since

$\displaystyle{\left\|\bm{u}_{\mathbf{X}}\right\|}_{2}{\left\|\bm{v}_{\mathbf{X}}\right\|}_{2}\lesssim\sqrt{\log
n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\cdot\left({\left\|\mathbf{X}\mathbf{X}_{j,:}\right\|}_{2}+{\left\|\mathbf{X}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\right)\lesssim\left(\log
n\right)^{\nicefrac{{3}}{{2}}}\sqrt{np}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}.\vspace{1mm}$

The proof is then completed by combining (B.2) and (32). ∎

######  Lemma 7.

Conditional on the intersection of events ${\cal F}_{6}\bigcap{\cal F}_{7}$
and fixing the indices $\pi^{\natural}(i)$ and $j$ $(j\neq\pi^{\natural}(i))$,
we know
$\mathsf{Term}_{2}\lesssim\sigma(\log^{2}n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}$
hold with probability at least $1-n^{-c}$, where $\mathsf{Term}_{2}$ is
defined in (24).

###### Proof.

Following a similar proof strategy as in Lemma 7, we first invoke the union
bound and obtain

$\displaystyle\mathbb{P}\left(\mathsf{Term}_{2}\gtrsim\sigma(\log
n)^{2}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle\leq~{}$
$\displaystyle\mathbb{P}\left(\mathsf{Term}_{2}\gtrsim\sigma(\log
n)^{2}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},~{}{\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\lesssim(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle+~{}$
$\displaystyle\mathbb{P}\left({\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\gtrsim(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)$
$\displaystyle\leq~{}$
$\displaystyle\underbrace{\mathbb{P}\left(\mathsf{Term}_{2}\gtrsim\sigma(\log
n){\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\right)}_{\triangleq~{}\zeta_{1}}+\underbrace{\mathbb{P}\left({\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\gtrsim(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)}_{\triangleq~{}\zeta_{2}}.$
(33)

The following analysis separately investigates the two probabilities
$\zeta_{1}$ and $\zeta_{2}$.

Phase I. Exploiting the independence between $\mathbf{X}$ and $\mathbf{W}$, we
can bound $\zeta_{1}$ as

$\displaystyle\zeta_{1}=\mathbb{E}_{\mathbf{X}}\mathbb{E}_{\mathbf{W}}\mathbbm{1}\left(\mathsf{Term}_{2}\gtrsim\sigma(\log
n){\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}\mathbb{E}_{\mathbf{X}}\exp\left(-c_{0}\log
n\right)=n^{-c_{0}},$ (34)

where in \raisebox{-.8pt}{1}⃝ we use the fact that $\mathsf{Term}_{2}$ is a
Gaussian RV
$\mathsf{N}(0,\|\widetilde{\mathbf{B}}^{\top}(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:})\|_{2})$
conditional on $\mathbf{X}$.

Phase II. Then we bound term $\zeta_{2}$. Note that

$\displaystyle{\left\|\widetilde{\mathbf{B}}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\leq~{}$
$\displaystyle{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}+{\left\|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}$
$\displaystyle\leq~{}$
$\displaystyle{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{j,:}\right\|}_{2}+{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}$
$\displaystyle+~{}$
$\displaystyle{\left\|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2},\vspace{1mm}$

we conclude

$\displaystyle\zeta_{2}\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\leq}}~{}$
$\displaystyle\underbrace{\mathbb{P}\left({\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{j,:}\right\|}_{2}+{\left\|(\widetilde{\mathbf{B}}-\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)})^{\top}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}\gtrsim\frac{p\log^{\nicefrac{{3}}{{2}}}n}{n}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)}_{\triangleq~{}\zeta_{2,1}}$
$\displaystyle+~{}$
$\displaystyle\underbrace{\mathbb{P}\left({\left\|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\left(\mathbf{X}_{j,:}-\mathbf{X}_{\pi^{\natural}(i),:}\right)\right\|}_{2}\gtrsim(\log
n){\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\right)}_{\triangleq~{}\zeta_{2,2}},$
(35)

where in \raisebox{-.8pt}{2}⃝ we use the fact $n\gtrsim p\log^{6}n$. Recalling
the definition of ${\cal F}_{7}$ then yields $\zeta_{2,1}=0$. For term
$\zeta_{2,2}$, we exploit the independence between
$\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}$ and
$\mathbf{X}_{j,:}$, $\mathbf{X}_{\pi^{\natural}(i),:}$. Via the Hanson-Wright
inequality (Theorem $6.2.1$ in Vershynin (2018)), we have

$\displaystyle\zeta_{2,2}\leq\exp\left[-c_{0}\left(\frac{(\log
n)^{2}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2}}{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{OP}}}\wedge\frac{(\log
n)^{4}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{4}}{{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}^{\top}\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}^{2}}\right)\right]\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{3}⃝}}}}{{\leq}}n^{-c},$
(36)

where \raisebox{-.8pt}{3}⃝ is due to the fact

$\displaystyle{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\leq{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}+{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\widetilde{\mathbf{B}}_{\setminus(\pi^{\natural}(i),j)}-\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}}\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{4}⃝}}}}{{\lesssim}}{\left|\kern-1.72218pt\left|\kern-1.72218pt\left|\mathbf{B}^{\natural}\right|\kern-1.72218pt\right|\kern-1.72218pt\right|}_{\textup{F}},\vspace{1mm}$

and \raisebox{-.8pt}{4}⃝ is conditional on event ${\cal F}_{6}$ and $n\gtrsim
p\log^{6}n$. Combining (B.2), (34), (B.2), and (36) then completes the proof.
∎

######  Lemma 8.

Conditional on event ${\cal F}_{2}\bigcap{\cal F}_{5}$ and fixing the indices
$\pi^{\natural}(i)$ and $j$, we have
$\mathsf{Term}_{3}\lesssim\frac{mp(\log^{2}n)\sigma^{2}}{n}+(\log^{2}n)\sigma^{2}\sqrt{\frac{mp}{n}}$
hold with probability exceeding $1-c_{0}n^{-c_{1}}$, where $\mathsf{Term}_{3}$
is defined in (25), and $c_{0},c_{1}>0$ are some fixed positive constants.

###### Proof.

For the benefits of presentation, we first define
$\bm{\Xi}^{\pi^{\natural}(i),j}$ as

$\displaystyle\bm{\Xi}^{\pi^{\natural}(i),j}=\mathbf{X}\left(\mathbf{X}_{\pi^{\natural}(i),:}-\mathbf{X}_{j,:}\right).\vspace{1mm}$

Then we can rewrite $\mathsf{Term}_{3}$ as
$(n-h)^{-1}\mathbf{W}_{i,:}^{\top}\mathbf{W}^{\top}\bm{\Omega}^{\pi^{\natural}(i),j}$
and expand it as

$\displaystyle\left\lvert\mathsf{Term}_{3}\right\rvert=~{}$
$\displaystyle\left(n-h\right)^{-1}\left\lvert\Xi^{\pi^{\natural}(i),j}_{i}\mathbf{W}_{i,:}^{\top}\mathbf{W}_{i,:}+\mathbf{W}_{i,:}^{\top}\bigg{(}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\bigg{)}\right\rvert$
$\displaystyle\leq~{}$
$\displaystyle\frac{1}{n-h}\left\lvert\Xi^{\pi^{\natural}(i),j}_{i}\right\rvert\cdot{\left\|\mathbf{W}_{i,:}\right\|}_{2}^{2}+\frac{1}{n-h}\left\lvert\left\langle\mathbf{W}_{i,:},\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\rangle\right\rvert$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{1}⃝}}}}{{\leq}}~{}$
$\displaystyle\frac{p\log
n}{n-h}{\left\|\mathbf{W}_{i,:}\right\|}_{2}^{2}+\frac{1}{n-h}\left\lvert\left\langle\mathbf{W}_{i,:},\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\rangle\right\rvert,\vspace{1mm}$

where in \raisebox{-.8pt}{1}⃝ we condition on event ${\cal F}_{2}$ and have
$\left\lvert\Xi^{\pi^{\natural}(i),j}_{i}\right\rvert\leq{\left\|\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}^{2}+{\left\|\mathbf{X}_{j,:}\right\|}_{2}^{2}\lesssim
p\log n$. With the union bound, we obtain

$\displaystyle\mathbb{P}\left(\mathsf{Term}_{3}\gtrsim\frac{mp(\log^{2}n)\sigma^{2}}{n}+\sigma^{2}(\log^{2}n)\sqrt{\frac{mp}{n}}\right)$
$\displaystyle\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{2}⃝}}}}{{\leq}}~{}$
$\displaystyle\underbrace{\mathbb{P}\left(\frac{p\log
n}{n-h}{\left\|\mathbf{W}_{i,:}\right\|}_{2}^{2}\gtrsim\frac{mp(\log^{2}n)\sigma^{2}}{n}\right)}_{\triangleq\zeta_{1}}+\underbrace{\mathbb{P}\left(\frac{1}{n-h}\left\lvert\left\langle\mathbf{W}_{i,:},\sum_{k\neq
i}\Omega^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\rangle\right\rvert\gtrsim\sigma^{2}(\log^{2}n)\sqrt{\frac{mp}{n}}\right)}_{\triangleq\zeta_{2}}.$
(37)

Then we separately bound the two terms $\zeta_{1}$ and $\zeta_{2}$.

Phase I. For term $\zeta_{1}$, we have

$\displaystyle\zeta_{1}\leq\mathbb{P}\left({\left\|\mathbf{W}_{i,:}\right\|}_{2}^{2}\gtrsim
m(\log
n)\sigma^{2}\right)\stackrel{{\scriptstyle\small{\raisebox{0.6pt}{\raisebox{-.8pt}{3}⃝}}}}{{=}}e^{-c_{0}\log
n}=n^{-c_{0}},$ (38)

where in \raisebox{-.8pt}{3}⃝ we use the fact that
${\left\|\mathbf{W}_{i,:}\right\|}_{2}^{2}/\sigma^{2}$ is a $\chi^{2}$-RV with
freedom $m$ and invoke Lemma 20.

Phase II. Then we upper-bound $\zeta_{2}$ as

$\displaystyle\zeta_{2}\leq~{}$
$\displaystyle\underbrace{\mathbb{P}\left(\frac{1}{n-h}\left\lvert\left\langle\mathbf{W}_{i,:},\sum_{k\neq
i}\Omega^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\rangle\right\rvert\gtrsim\frac{\sigma\sqrt{\log
n}}{n}\bigg{\|}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\bigg{\|}_{2}\right)}_{\triangleq~{}\zeta_{2,1}}$
$\displaystyle+~{}$
$\displaystyle\underbrace{\mathbb{P}\left(\bigg{\|}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\bigg{\|}_{2}^{2}\gtrsim
mnp(\log n)^{3}\sigma^{2}\right)}_{\triangleq~{}\zeta_{2,2}}.$ (39)

For term $\zeta_{2,1}$, we exploit the independence across the rows of the
matrix $\mathbf{W}$. Conditional on $\left\\{\mathbf{W}_{k,:}\right\\}_{k\neq
i}$, we conclude the inner-product $\left\langle\mathbf{W}_{i,:},\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\rangle$ to be a Gaussian
RV with zero mean and $\big{\|}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\big{\|}_{2}^{2}$ variance,
which yields $\zeta_{2,1}\leq n^{-c}$. For term $\zeta_{2,2}$, we have

$\displaystyle\zeta_{2,2}\leq~{}$ $\displaystyle\
\underbrace{\mathbb{P}\left(\bigg{\|}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\bigg{\|}_{2}^{2}\gtrsim m(\log
n)\sigma^{2}\bigg{[}\sum_{k\neq
i}(\Xi^{\pi^{\natural}(i),j}_{k})^{2}\bigg{]},~{}\sum_{k\neq
i}(\Xi^{\pi^{\natural}(i),j}_{k})^{2}\lesssim(\log^{2}n)np\right)}_{\triangleq~{}\zeta_{2,2,1}}$
$\displaystyle+~{}$ $\displaystyle\underbrace{\mathbb{P}\left(\sum_{k\neq
i}(\Xi^{\pi^{\natural}(i),j}_{k})^{2}\gtrsim(\log^{2}n)np\right)}_{\triangleq~{}\zeta_{2,2,2}}.$
(40)

Due to the independence across $\mathbf{X}$ and $\mathbf{W}$, we can verify
${\left\|\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\right\|}_{2}^{2}/[\sigma^{2}\sum_{k\neq
i}(\Xi^{\pi^{\natural}(i),j}_{k})^{2}]$ to be a $\chi^{2}$-RV with freedom $m$
when conditional on $\mathbf{X}$. Invoking Lemma 20, we can upper-bound
$\xi_{2,2,1}$ as

$\displaystyle\zeta_{2,2,1}\leq\mathbb{P}\left(\bigg{\|}\sum_{k\neq
i}\Xi^{\pi^{\natural}(i),j}_{k}\mathbf{W}_{k,:}\bigg{\|}_{2}^{2}\gtrsim m(\log
n)\sigma^{2}\bigg{[}\sum_{k\neq
i}(\Xi^{\pi^{\natural}(i),j}_{k})^{2}\bigg{]}\right)\leq n^{-c}.$ (41)

As for $\xi_{2,2,2}$, we condition on event ${\cal F}_{5}$ and have

$\displaystyle\zeta_{2,2,2}\leq\mathbb{P}\left({\left\|\mathbf{X}\mathbf{X}_{\pi^{\natural}(i),:}\right\|}_{2}+{\left\|\mathbf{X}\mathbf{X}_{j,:}\right\|}_{2}\gtrsim(\log
n)\sqrt{np}\right)=0.$ (42)

Then the proof is completed by combining (B.2), (38), (B.2), (B.2), (41), and
(42). ∎