\section{Introduction}

\subsection{The Statistical Problem of Feature Matching}

In this paper, we present a rigorous statistical analysis of the problem of permutation estimation and multiple feature
matching from noisy observations. More precisely, let $\{X_1,\ldots,X_n\}$ and $\{X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m\}$ be
two sets of vectors from $\mathbb R^d$, hereafter referred to as noisy features, containing many matching elements. That is,
for many $X_i$'s there is a $X^{\text{\tt\#}}_j$ such that $X_i$ and $X^{\text{\tt\#}}_j$ coincide up to an observation noise (or measurement
error). Our goal is to estimate an application $\pi^* : \{1,\ldots,n\}\to \{1,\ldots,m\}$ for which each $X_i$ matches
with $X^{\text{\tt\#}}_{\pi^*(i)}$ and to provide tight conditions which make it possible to accurately recover $\pi^*$ from data.

In order to define a statistical framework allowing us to compare different estimators of $\pi^*$, we
confine\footnote{These assumptions are imposed for the purpose of getting transparent theoretical results and are in no
way necessary for the validity of the considered estimation procedures, as discussed later in the paper.} our attention to
the case $n=m$, that is when the two sets of noisy features have equal sizes.  Furthermore, we assume that there exists a
unique permutation of $\{1,\ldots,n\}$, denoted $\pi^*$, leading to pairs of features $(X_i,X^{\text{\tt\#}}_{\pi^*(i)})$ that match
up to a measurement error. In such a situation, it is clearly impossible to recover the true permutation $\pi^*$ if
some features within the set $\{X_1,\ldots,X_n\}$ are too close. Based on this observation, we propose to measure
the quality of a procedure of permutation estimation by the minimal distance between pairs of different features for which
the given procedure is still consistent. This quantity will be called \textit{separation distance} and will be the
main concept of interest in the present study. In this respect, the approach we adopted is close in spirit to
the minimax theory of hypotheses testing (see, for instance, \cite{Spokoiny96, Ingster}).

\subsection{A Motivating Example: Feature Matching in Computer Vision}

Many tasks of computer vision, such as object recognition, motion tracking or structure from motion, are currently
carried out using algorithms that contain a step of feature matching, cf.~\cite{Szeliski, Hartley_Zisserman}.
The features are usually local descriptors that serve to summarize the images. The most famous examples of such
features are perhaps SIFT \citep{Lowe2004} and SURF \citep{surf}. Once the features have been computed for each image,
an algorithm is applied to match features of one image to those of another one. The matching pairs are then used
for estimating the deformation of the object, for detecting the new position of the followed object, for creating a
panorama, etc. In this paper, we are interested in simultaneous matching of a large number of features. The main focus
is on the case when the two sets of features are extracted from the images that represent the same scene with a large
overlap, and therefore the sets of features are (nearly) of the same size and every feature in the first image is
also present in the second one. This problem is made more difficult by the presence of noise in the images, and
thus in the features as well. Typically, due to the high resolution of most images, the number of features is large
and their dimension is relatively large as well (128 for SIFT and 64 for SURF). It is therefore important to
characterize the behavior of various matching procedures as a function of the number of features, the dimension and
the noise level.

\subsection{Main Contributions}

We consider four procedures of permutation estimation that naturally arise in this context. The first one is a greedy
procedure that sequentially assigns to each feature $X_i$ the closest feature $X^{\text{\tt\#}}_j$ among those features that have
not been assigned at an earlier step. The three other estimators are defined as minimizers of the (profiled-)log-likelihood under
three different modeling assumptions. These three modeling assumptions are that the
noise level is constant across all the features (homoscedastic noise), that the noise level is variable (heteroscedastic noise) but known
and that the noise level is variable and unknown. The corresponding estimators are respectively called least sum of squares (LSS) estimator,
least sum of normalized squares (LSNS) estimator and least sum of logarithms (LSL) estimator.

We first consider the homoscedastic setting and show that all the considered estimators are consistent under similar
conditions on the minimal distance between distinct features $\kappa$. These conditions state that $\kappa$ is larger than some function of
the noise level $\sigma$, the sample size $n$ and the dimension $d$. This function is the same for the four aforementioned procedures
and is given, up to a multiplicative factor, by
\begin{align}\label{kappastar}
\kappa^*(\sigma,n,d) = \sigma \max((\log n)^{1/2}, (d\log n)^{1/4}).
\end{align}
Then, we prove that this expression provides the optimal rate of the separation distance in the sense that for some
absolute constant $c$  if $\kappa \le c\kappa^*(\sigma,n,d)$ then there is no procedure capable of consistently estimating $\pi^*$.

In the heteroscedastic case, we provide an upper bound on the identifiability threshold ensuring the consistency of the LSNS and
LSL estimators. Up to a proper normalization by the noise level, this bound is of the same form as (\ref{kappastar}) and, therefore,
the ignorance of the noise level does not seriously affect the quality of estimation. Furthermore, the LSL estimator is easy
to adapt to the case $n\not = m$ and is robust to the presence of outliers in the features. We carried out a small experimental
evaluation that confirms that in the heteroscedastic setting the LSL estimator is as good as the LSNS (pseudo-) estimator and
that they outperform the two other estimators: the greedy estimator and the least sum of squares. We also argue that the three
estimators stemming from the maximum likelihood methodology are efficiently computable either by linear programming or by the
Hungarian algorithm.

Note that different loss functions may be used for measuring the distance between an estimated permutation and the true one. Most results
of this paper are established for the 0\,-1 loss, which equals one if the estimator and the true permutation differ at least at one location
and equals 0 otherwise. However, it is of interest to analyze the situation with the Hamming distance as well, since it amounts to controlling
the proportion of the mismatched features and, hence, offers a more graduated evaluation of the quality of estimation.  We show that in the
case of the Hamming distance, in the regime of moderately large dimension (\textit{i.e., } $d\ge c\log n$ for some constant $c>0$) the rate of separation
is exactly the same as in the case of the 0\,-1 distance. The picture is more complex in the regime of small dimensionality $d=o(\log(n))$,
in which we get the same upper bound as for the 0\,-1 loss but the lower bound is expressed in terms of the logarithm of the packing number
of an $\ell_2$ ball of the symmetric group. We conjecture that this quantity is of the order of $n\log(n)$ and check this conjecture for
relatively small values of $n$. If this conjecture is correct, our lower bound coincides up to a multiplicative factor with the upper bound.

Finally, let us mention that some of the results of the present work have been presented in the AI-STATS 2013
conference and published
in the proceedings \citep{jmlr_CD13}.

\subsection{Plan of the Paper}

We introduce in Section~\ref{section:model} a model for the problem of matching two sets of features and of estimation of a permutation.
The estimating procedures analyzed in this work are presented in Section~\ref{section:estimators}, whereas their performances in terms
of rates of separation distance are described in Section~\ref{section:performance}. In Section~\ref{section:computationalaspects}, computational aspects of the
estimating procedures are discussed while Section~\ref{section:extensions} is devoted to the
statement of some extensions of our results. We report in Section~\ref{section:simulations} the results of some numerical experiments.
The proofs of the theorems and of the lemmas are postponed to Sections~\ref{section:theorems} and~\ref{section:lemmas}, respectively.

\section{Notation and Problem Formulation}\label{section:model}

We begin with formalizing the problem of matching two sets of features $\{X_1,\ldots,X_n\}$ and $\{X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m\}$ with $n,m\ge 2$.
In what follows we assume that the observed features are randomly generated from the model
\begin{equation}\label{model}
\begin{cases}
X_i = \theta_i + \sigma_i\xi_i \ , \\
X^{\text{\tt\#}}_j = {\theta}^\text{\tt\#}_j + \s^{\text{\tt\#}}_j\xi^\text{\tt\#}_j,
\end{cases}\quad i=1,\ldots,n \text{ and } j=1,\ldots,m
\end{equation}
where
\begin{itemize}\setlength{\itemsep}{-2pt}
\item $\boldsymbol{\theta}=\{\theta_1,\ldots,\theta_n\}$ and $\boldsymbol{\theta}^\text{\tt\#}=\{{\theta}^\text{\tt\#}_1,\ldots,{\theta}^\text{\tt\#}_m\}$ are two collections of vectors from $\mathbb R^d$, corresponding to
the original features, which are unavailable,
\item $\sigma_1,\ldots,\sigma_n,\s^{\text{\tt\#}}_1,\ldots,\s^{\text{\tt\#}}_m$ are positive real numbers corresponding to the levels of noise contaminating each feature,
\item $\xi_1,\ldots,\xi_n$ and $\xi^\text{\tt\#}_1,\ldots,\xi^\text{\tt\#}_m$ are two independent sets of i.i.d.\ random vectors drawn from
    the Gaussian distribution with zero mean and identity covariance matrix.
\end{itemize}

The task of feature matching consists in finding a bijection $\pi^*$ between the largest possible subsets $S_1$ and $S_2$ of $\{1,\ldots,n\}$ and $\{1,\ldots,m\}$ respectively, such that
\begin{equation}\label{pistar}
\forall\ i\in S_2,\quad {\theta}^\text{\tt\#}_i \equiv \theta_{\pi^*(i)},
\end{equation}
where $\equiv$ is an equivalence relation that we call \textit{matching criterion}. The features that do not belong to $S_1$ or $S_2$ are called \textit{outliers}. To ease presentation, we mainly focus on the case where the matching criterion is the equality of two vectors. However, as discussed in Section~\ref{subsection5.2}, most results carry over the equivalence corresponding to equality of two vectors transformed by a given linear transformation.
Furthermore, it turns out that statistical inference for the matching problem is already quite involved when no outlier is present in the data. Therefore,
we make the following assumption
\begin{equation}\label{m=n}
m=n \qquad\text{and}\qquad S_1=S_2=\{1,\ldots,n\}.
\end{equation}
Note however that the procedures we consider below admit natural counterparts in the setting with outliers. We will also restrict ourselves to
noise levels satisfying some constraints. The two types of constraints we consider, referred to as homoscedastic and heteroscedastic setting, correspond to
the relations $\sigma_i=\s^{\text{\tt\#}}_i=\sigma$, $\forall i=1,\ldots,n$, and $\sigma_{\pi^*(i)}=\s^{\text{\tt\#}}_i$, $\forall i=1,\ldots,n$.

In this formulation, the data generating distribution is defined by the (unknown) parameters $\boldsymbol{\theta}$,
$\boldsymbol{\sigma} = (\sigma_i,\ldots,\sigma_n)$ and $\pi^*$.
In the problem of matching, we focus our attention on the problem of estimating the parameter $\pi^*$ only,
considering $\boldsymbol{\theta}$ and $\boldsymbol{\sigma}$ as nuisance parameters. In what follows, we denote by  $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}$ the probability distribution
of the vector $(X_1,\ldots,X_n,X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_n)$ defined by (\ref{model}) under the conditions
(\ref{pistar}) and (\ref{m=n}). We write $\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}$ for the expectation
with respect to $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}$. The symmetric group, \textit{i.e., } the set of all permutations
of $\{1,\ldots,n\}$, will be denoted by $\mathfrak{S}_n$.

We will use two measures of quality for quantifying the error of an estimator $\hat\pi$ of the permutation $\pi^*$. These
errors are defined as the 0\,-1 distance and the normalized Hamming distance between $\hat\pi$ and $\pi^*$, given by
\begin{equation}\label{dist}
\delta_{\text{0\,-1}}(\hat\pi,\pi^*)\triangleq\mathds{1}_{\{\hat\pi\not=\pi^*\}},\qquad
\delta_H(\hat\pi,\pi^*)\triangleq\frac1n\sum_{k=1}^n \mathds{1}_{\{\hat\pi(k)\not=\pi^*(k)\}}.
\end{equation}
Our ultimate goal is to design estimators that have an expected error smaller than a prescribed level $\alpha$
under the weakest possible conditions on the nuisance parameter $\boldsymbol{\theta}$. The estimation of the permutation or, equivalently,
the problem of matching is more difficult when the features are hardly distinguishable. To quantify this phenomenon, we introduce
the separation distance $\kappa(\boldsymbol{\theta})$ and the relative separation distance $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})$, which measure the minimal distance
between distinct features and the minimal distance-to-noise ratio, respectively. The precise definitions are
\begin{equation}
\kappa(\boldsymbol{\theta}) \triangleq \min_{i\neq j} \|\theta_i-\theta_j\|,\qquad
\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma}) \triangleq \min_{i\neq j} \frac{\|\theta_i-\theta_j\|}{(\sigma_i^2+\sigma_j^2)^{1/2}}. \label{kappa}
\end{equation}
We will see that in the heteroscedastic case the last quantity is more suitable for characterizing the behavior of the estimators
than the first one.

Clearly, if $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})=0$ and $\sigma_i$'s are all equal, then the parameter $\pi^*$ is nonidentifiable, in the sense that there exist
two different permutations $\pi^*_1$ and $\pi^*_2$ such that the distributions $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*_1}$
and $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*_2}$ coincide. Therefore, the condition $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})>0$ is necessary
for the existence of consistent estimators of $\pi^*$. Furthermore, good estimators are those consistently estimating $\pi^*$
even if $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})$ is small. To give a precise sense to these considerations, let $\alpha\in(0,1)$ be a prescribed
tolerance level and let us call \textit{perceivable separation distance} of a given estimation procedure $\hat\pi$ the quantity
\begin{equation*}
\bar\kappa_\alpha(\hat\pi) \triangleq \inf\Big\{\kappa>0 \,\big|\, \max_{\pi\in\mathfrak{S}_n} \sup_{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})>\kappa} \mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi} (\hat\pi\neq\pi)\le\alpha\Big\},
\end{equation*}
where we skip the dependence on $n$, $d$ and $\boldsymbol{\sigma}$. Here, the perceivable separation distance
is defined with respect to the 0\,-1 distance, the corresponding definition for the Hamming distance is obtained
by replacing $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi} (\hat\pi\neq\pi)$ by $\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi} [\delta_H(\hat\pi,\pi)]$.
Finally, we call \textit{minimax separation distance} the smallest possible {perceivable separation distance} achieved by an estimator
$\hat\pi$, \textit{i.e., }
\begin{equation}\label{mmx}
\bar\kappa_\alpha \triangleq \inf_{\hat\pi}~\bar\kappa_\alpha(\hat\pi),
\end{equation}
where the infimum is taken over all possible estimators of $\pi^*$. In the following sections, we establish nonasymptotic upper and lower bounds on
the minimax separation distance which coincide up to a multiplicative constant independent of $n$, $d$ and $\boldsymbol{\sigma}$. We also
show that a suitable version of the maximum profiled likelihood estimator is minimax-separation-rate-optimal both in the homoscedastic and heteroscedastic
settings.

\section{Estimation Procedures}\label{section:estimators}

As already mentioned, we will consider four estimators. The simplest one, called greedy algorithm and
denoted by $\pi^{\rm gr}$ is defined as follows: $\pi^{\rm gr}(1) = \arg\min_{j\in\{1,\ldots,n\}} \|X_j-X^{\text{\tt\#}}_1\|$
and, for every $i\in\{ 2,\ldots,n\}$, recursively define
\begin{align}\label{greedy}
\pi^{\rm gr}(i)\triangleq \arg\min_{j\not\in\{\pi^{\rm gr}(1),\ldots,\pi^{\rm gr}(i-1)\}} \|X_j-X^{\text{\tt\#}}_i\|.
\end{align}
A drawback of this estimator is that it is not symmetric: the resulting permutation depends on the initial numbering of the
features. However, we will show that this estimator is minimax-separation-rate-optimal in the homoscedastic setting.

A common approach for avoiding incremental estimation and taking into consideration all the observations at the same time
consists in defining the estimator $\hat\pi$ as a maximizer of the profiled likelihood. In the homoscedastic case, which is
the first setting we studied, the computations lead to the estimator
\begin{equation}
\pi^{\textrm{LSS}}  \triangleq \arg\min_{\pi\in\mathfrak{S}_n} \sum_{i=1}^n \|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2. \label{LSS}
\end{equation}
which will be referred to as the \textit{Least Sum of Squares} (LSS) estimator.

The LSS estimator takes into account the distance between the observations irrespectively of the noise levels.
The fact of neglecting the noise levels, while harmless in the homoscedastic setting, turns out to cause serious
loss of efficiency in terms of the perceivable distance of separation in the setting of heteroscedastic noise.
Yet, in the latter setting, the distance between the observations is not as relevant as the signal-to-noise ratio
${\|\theta_i-\theta_j\|^2}/({\sigma_i^2+\sigma_j^2})$. Indeed, when the noise levels are small, two noisy but distinct vectors
are easier to distinguish than when the noise levels are large.

The computation of the maximum likelihood estimator in the heteroscedastic case with known noise levels also suggests that
the signal-to-noise ratio should be taken into account. In this setting, the likelihood maximization leads to the
\textit{Least Sum of Normalized Squares} (LSNS) estimator
\begin{equation}
\pi^{\textrm{LSNS}} \triangleq \arg\min_{\pi\in\mathfrak{S}_n}\sum_{i=1}^n \frac{\|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2}{\sigma_{\pi(i)}^2+{\s^{\text{\tt\#}}_i}^2}. \label{LSNS}
\end{equation}
We will often call the LSNS a pseudo-estimator, to underline the fact that it requires the knowledge of the noise levels $\sigma_i$ which
are generally unavailable.

In the general setting, when no information on the noise levels is available, the likelihood is maximized over all nuisance parameters (features and noise levels). But this problem is underconstrained, and the result of this maximization is ${+\infty}$. This can be circumvented by assuming a proper relation between the noise levels. As mentioned earlier, we chose the  assumption
\begin{equation}\label{sigma_h}
\forall\ i\in\{1,\ldots,n\},\quad \s^{\text{\tt\#}}_i = \sigma_{\pi^*(i)}.
\end{equation}
The maximum likelihood estimator under this constraint is the \textit{Least Sum of Logarithms} (LSL) defined as
\begin{equation}
\pi^{\textrm{LSL}}  \triangleq \arg\min_{\pi\in\mathfrak{S}_n} \sum_{i=1}^n \log \|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2. \label{LSL}
\end{equation}
We will prove that this estimator is minimax-rate-optimal both in the homoscedastic and the heteroscedastic cases.

\section{Performance of the Estimators}\label{section:performance}

The purpose of this section is to assess the quality of the aforementioned procedures. To this end, we present
conditions for the consistency of these estimators in the form of upper bounds on their perceivable separation
distance. Furthermore, to compare this bounds with the minimax separation distance, we establish lower bounds on the latter
and prove that it coincides up to a constant factor with the perceivable separation  distance of the LSL estimator.

\subsection{Homoscedastic Setup}
We start by considering the homoscedastic case, in which upper and lower bounds matching up to a constant are obtained for all the estimators
introduced in the previous section.

\begin{theorem}\label{upperbound1}
Let $\alpha\in (0,1)$ be a tolerance level and $\sigma_j=\s^{\text{\tt\#}}_j=\sigma$ for all $j\in\{1,\ldots,n\}$.
If $\hat\pi$ denotes any one of the estimators (\ref{greedy})-(\ref{LSL}), then
\begin{equation*}
\bar\kappa_\alpha(\hat\pi)\le 4 \max\Big\{\Big(2\log \frac{8n^2}{\alpha}\Big)^{1/2}, \Big(d\log\frac{4n^2}{\alpha}\Big)^{1/4} \Big\}.
\end{equation*}
\end{theorem}

An equivalent way of stating this result is that if
\begin{equation*}
\kappa = 4\sigma \max\Big\{\Big(4\log \frac{8n^2}{\alpha}\Big)^{1/2}, \Big(4d\log\frac{4n^2}{\alpha}\Big)^{1/4} \Big\}
\end{equation*}
and $\Theta_\kappa$ is the set of all $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ such that $\kappa(\boldsymbol{\theta})\ge \kappa$, then
\begin{equation*}
\max_{\pi^*\in\mathfrak{S}_n}\sup_{\boldsymbol{\theta}\in\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}(\hat\pi\neq\pi^*) \le \alpha
\end{equation*}
for all the estimators defined in Section~\ref{section:estimators}.
Note that this result is nonasymptotic. Furthermore, it tells us that the perceivable separation distance of the
procedures under consideration is at most of the order of
\begin{equation}
\max\Big\{(\log n)^{1/2}, (d \log n)^{1/4}\Big\}.
\end{equation}
It is interesting to observe that there are two regimes in this rate, the boundary of which corresponds to the case where $d$ is of the order of $\log n$.
For dimensions that are significantly smaller than $\log n$, the perceivable distance of separation is dimension free. On the other hand,
when $d$ is larger than $c\log n$ for some absolute constant $c>0$, the perceivable distance of separation deteriorates with increasing
$d$ at the polynomial rate $d^{1/4}$. However, this result does not allow us to deduce any hierarchy between the four estimators, since it
provides the same upper bound for all of them. Moreover, as stated in the next theorem, this bound is optimal up to a multiplicative constant.

\begin{theorem}\label{lowerbound1}
Assume that $n\ge6$, $\Theta_\kappa$ is the set of all $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ such that $\kappa(\boldsymbol{\theta})\ge \kappa$.
Then there exist two absolute constants $c,C>0$ such that for
\begin{equation*}
\kappa=2^{-5/2}{\sigma}\max\Big\{(\log n)^{1/2}, (c d \log n)^{1/4}\Big\},
\end{equation*}
the following lower bound holds
\begin{equation*}
\inf_{\hat\pi} \max_{\pi^*\in\mathfrak{S}_n} \sup_{\boldsymbol{\theta}\in\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}\big(\hat\pi\neq\pi^*\big) >  C,
\end{equation*}
where the infimum is taken over all permutation estimators.
\end{theorem}

An equivalent way of stating this result is to say that, for any $\alpha\le C$, the minimax distance of separation satisfies the inequality
$$
\bar\kappa_\alpha \ge \frac18\max\big\{(\log n)^{1/2}, (c d \log n)^{1/4}\big\}.
$$
Combined with Theorem~\ref{upperbound1}, this implies that the minimax rate of separation is given by the expression
$\max\big\{(\log n)^{1/2}, (d \log n)^{1/4}\big\}$.

In order to avoid any possible confusion, we emphasize that the rate obtained in this and subsequent sections concerns the
speed of decay of the separation distance and not the estimation risk measured by $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}\big(\hat\pi\neq\pi^*\big)$.
For the latter, considering $\kappa$ as fixed, one readily derives from Theorem~\ref{upperbound1} that
\begin{align}\label{rem:a}
\max_{\pi^*\in\mathfrak{S}_n} \sup_{\boldsymbol{\theta}\in\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}\big(\hat\pi\neq\pi^*\big)
\le \max\Big\{8n^2\exp\Big(-\frac{\kappa^2}{2^{6}\sigma^2}\Big),4n^2\exp\Big(-\frac{\kappa^4}{2^{10}d\sigma^4}\Big)\Big\}
\end{align}
for the four estimators $\hat\pi$ defined in the previous section by equations (\ref{greedy})-(\ref{LSL}). We do not know whether the right-hand side of this
inequality is the correct rate for the minimax risk $R^{\rm mmx} =\inf_{\hat\pi}  \sup_{(\boldsymbol{\theta},\pi^*)\in\Theta_\kappa\times\mathfrak{S}_n}
\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}\big(\hat\pi\neq\pi^*\big)$. In fact, one can adapt the proof of Theorem~\ref{lowerbound1}
to get a lower bound on  $R^{\rm mmx}$ which is of the same form as the right-hand side
in (\ref{rem:a}), but with constants $2^6$ and $2^{10}$ replaced by smaller ones. The ratio of such a lower bound and the upper bound in
(\ref{rem:a}) tends to 0 and, therefore, does not provide the minimax rate of estimation, in the most common sense of the
term. However, one may note that the minimax rate of separation established in this work is the analogue of the
minimax rate of estimation of the Bahadur risk, see \cite{Bahadur,Korostelev,KorostelevSpok}.

\subsection{Heteroscedastic Setup}

We switch now to the heteroscedastic setting, which allows us to discriminate between the four procedures.
Note that the greedy algorithm, the LSS and the LSL have a serious advantage over the LSNS since they can
be computed without knowing the noise levels $\boldsymbol{\sigma}$.

\begin{theorem}\label{upperbound2}
Let $\alpha\in(0,1)$ and condition (\ref{sigma_h}) be fulfilled. If $\hat\pi$ is either $\pi^{\textrm{LSNS}}$ (if the noise
levels $\sigma_i$ are known) or $\pi^{\textrm{LSL}}$ (when the noise levels are unknown),
then
\begin{equation*}
\bar\kappa_\alpha(\hat\pi) \leq 4\max\Big\{\Big(2\log \frac{8n^2}{\alpha}\Big)^{1/2}, \Big(d\log\frac{4n^2}{\alpha}\Big)^{1/4} \Big\}.
\end{equation*}
\end{theorem}

This result tells us that the performance of the LSNS and LSL estimators, measured in terms of the order of magnitude of the
separation distance, is not affected by the heteroscedasticity of the noise levels. Two natural questions arise: 1) is this performance the best possible
over all the estimators of $\pi^*$ and 2) is the performance of the LSS and the greedy estimator as good as that of
the LSNS  and LSL?

To answer the first question, we start by discarding the degenerate situations where faster rates can be achieved by estimators that are
heavily based on the knowledge of the noise levels $\boldsymbol\sigma$. In fact, although considering $\boldsymbol{\sigma}$ as known limits considerably
the scope of applications of the methods, the definition (\ref{mmx}) involves a minimum over all possible estimators $\hat\pi$ which are
allowed to depend on $\boldsymbol{\sigma}$. For some specific noise levels $\boldsymbol\sigma$, from a purely
theoretical point of view, knowing $\boldsymbol{\sigma}$ may lead to a substantially better minimax rate and even
to a separation equal to zero, which corresponds to an estimator that has no real practical
interest.  Indeed, under condition~(\ref{sigma_h}), it is possible to estimate the permutation $\pi^*$ by fitting the
noise levels. For instance, we can define an estimator $\hat\pi$ as follows:
$\hat\pi(1) = \arg\min_{i=1,\ldots,n} \big| \frac1{2d}{\|X_i-X^{\text{\tt\#}}_1\|^2} - \sigma_i^2 \big|$
and, recursively, for every $j\in\{ 2,\ldots,n\}$,
\begin{equation*}
\hat\pi(j) = \arg\min_{i\not\in\{\hat\pi(1),\ldots,\hat\pi(j-1)\}} \Big| \frac1{2d}{\|X_i-X^{\text{\tt\#}}_j\|^2} - \sigma_i^2 \Big|.
\end{equation*}
This provides an accurate estimator of $\pi^*$---for vectors $\{\theta_j\}$ that are very close---as soon as the noise
levels are different enough from each other.
In particular, the estimated permutation $\hat\pi$ coincides with the true permutation $\pi^*$ on the event
\begin{equation*}
\forall i\in\{1,\ldots,n\},\qquad \Big|\frac1{2d}{\|X_{\pi^*(i)}-X^{\text{\tt\#}}_i\|^2} - \sigma_{\pi^*(i)}^2 \Big| < \min_{j\neq \pi^*(i)} \Big|  \frac{\|X_j-X^{\text{\tt\#}}_i\|^2}{2d} - \sigma_j^2\Big|,
\end{equation*}
which includes the event:
\begin{equation}\label{eq:14bis}
\Big|\frac1{2d}{\|X_{\pi^*(i)}-X^{\text{\tt\#}}_i\|^2} - \sigma_{\pi^*(i)}^2 \Big|+
\Big|  \frac{\|X_j-X^{\text{\tt\#}}_i\|^2}{2d} -  \frac{\sigma_j^2+\sigma_{\pi^*(i)}^2}2\Big|<  \frac{|\sigma_j^2-\sigma_{\pi^*(i)}^2|}2
\end{equation}
for all $(i,j)\in\{1,\ldots,n\}^2$ such that $j\neq \pi^*(i)$.
Using standard bounds on the tails of the $\chi^2$ distribution recalled in Lemma~\ref{concentration} of Section~\ref{section:theorems},
in the case when all the vectors $\theta_j$ are equal, one can check that the left-hand side in (\ref{eq:14bis}) is of the order
of $(\sigma_{\pi^*(i)}^2+\sigma_{j}^2)\max\big\{\sqrt{(\log n)/d},(\log n)/d\big\}$. This implies that we can consistently identify the
permutation when
\begin{equation*}
\forall (i,j)\in\{1,\ldots,n\}^2,
\quad i\not= j,\qquad \Big|\frac{\sigma_j^2}{\sigma_i^2}-1\Big| \ \gg \ \max\Big\{\sqrt{\frac{\log n}{d}},\frac{\log n}{d}\Big\},
\end{equation*}
even if the separation distance is equal to zero. In order to discard such kind of estimators from the competition in the procedure of determining the minimax rates, we restrict our attention to the noise levels for which the values $|\frac{\sigma_j^2}{\sigma_i^2}-1|$ are not
larger than $C\max\big\{\sqrt{(\log n)/d},(\log n)/d\big\}$ for $j\not=i$.

\begin{theorem}\label{lowerbound2}
Assume that $n\ge6$, $\bar\Theta_\kappa$ is the set of all $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ such that $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\ge \kappa$ and
\begin{equation*}
\frac{\max_i\sigma^2_i}{\min_i\sigma^2_i} -1 \le  \ \max\Big\{\sqrt{\frac{\log n}{16d}},\frac{\log n}{16d}\Big\}.
\end{equation*}
Then there exist two constants $c,C>0$ such that
$\kappa<(1/8) \max\{ (\log n)^{1/2},  (cd\log n)^{1/4} \}$,
implies that
\begin{equation*}
\inf_{\hat\pi} \max_{\pi^*\in\mathfrak{S}_n} \sup_{\boldsymbol{\theta}\in\bar\Theta_{\kappa}} \mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*} (\hat\pi\neq\pi^*) > C,
\end{equation*}
where the infimum is taken over all permutation estimators.
\end{theorem}

It is clear that the constants $c$ and $C$ of the previous theorem are closely related. The inspection of the proof shows that, for instance,
if $c\le 1/20$ then  $C$ is larger than $17\%$.

Let us discuss now the second question raised earlier in this section and concerning the theoretical properties
of the greedy algorithm and the LSS under heteroscedasticity. In fact, the perceivable distances of separation
of these two procedures are significantly worse than those of the LSNS and the LSL especially for large
dimensions $d$.  We state the corresponding result for the greedy algorithm, a similar conclusion
being true for the LSS as well. The superiority of the LSNS and LSL is also confirmed by the numerical
simulations presented in Section~\ref{section:simulations} below. In the next theorem and in the sequel of the paper,
we denote by $\textit{id}$ the identity permutation defined by $\textit{id}(i) = i$ for all $i\in\{1,\ldots,n\}$.

\begin{theorem}\label{lowerbound3}
Assume that $d\ge 225\log 6$, $n=2$, $\sigma_1^2 =3$ and $\sigma_2^2=1$.
Then the condition  $\kappa < 0.1(2d)^{1/2}$ implies that
$$
\sup_{\boldsymbol{\theta}\in\bar\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\textit{id}}(\pi^{\textrm{gr}}\neq\textit{id}) \ge 1/2.
$$
\end{theorem}

This theorem shows that if $d$ is large, the necessary condition for $\pi^{\textrm{gr}}$ to be
consistent is much stronger than the one obtained for $\pi^{\textrm{LSL}}$ in Theorem~\ref{upperbound2}. Indeed,
for the consistency of $\pi^{\textrm{gr}}$, $\kappa$ needs to be at least of the order of $d^{1/2}$, whereas
$d^{1/4}$ is sufficient for the consistency of $\pi^{\textrm{LSL}}$. Hence, the maximum likelihood estimators LSNS and LSL
that take into account noise heteroscedasticity are, as expected, more interesting than the simple greedy estimator\footnote{It
should be noted that the conditions of Theorem~\ref{lowerbound3} are not compatible with those of Theorem~\ref{lowerbound2}.
Hence, strictly speaking, the former does not imply that the greedy estimator is not minimax under the conditions of the latter.}.

\section{Computational Aspects}\label{section:computationalaspects}

At first sight, the computation of the estimators (\ref{LSS})-(\ref{LSL}) requires to
perform an exhaustive search over the set of all possible permutations, the number of which, $n!$, is
prohibitively large. This is in practice impossible to do on a standard PC as soon as $n \geq 20$.
In this section, we show how to compute these (maximum likelihood) estimators in polynomial time using, for instance,
algorithms of linear programming\footnote{The idea of reducing the problem of permutation estimation to a linear program
has been already used in the literature, however, without sufficient theoretical justification: see, for instance,
\citet{Jebara2003}.}.

To explain the argument, let us consider the LSS estimator
\begin{equation*}
\pi^{{\textrm{LSS}}} = \arg\min_{\pi\in\mathfrak{S}_n} \sum_{i=1}^n \|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2.
\end{equation*}
For every permutation $\pi$, we denote by $P^\pi$ the $n\times n$ permutation matrix with coefficients
$P^\pi_{ij} = \mathds{1}_{\{ j=\pi(i) \}}$. Then we can give the equivalent formulation
\begin{equation}\label{discrete}
\pi^{\textrm{LSS}} = \arg\min_{\pi \in\mathfrak{S}_n} \mathbf{tr} \big(M P^\pi\big),
\end{equation}
where $M$ is the matrix with coefficient $\|X_i-X^{\text{\tt\#}}_j\|^2$ at the $i^{th}$ row and $j^{th}$ column. The cornerstone of our next argument is the Birkhoff-von Neumann theorem stated below, which can be found for example in \citep{BudishCheKojimaMilgrom2009}.

\begin{theorem}[Birkhoff-von Neumann Theorem]\label{birkhoff-von-neumann}
Assume that $\mathcal{P}$ is the set of all doubly stochastic matrices of size $n$, \textit{i.e., } the matrices whose
entries are nonnegative and sum up to $1$ in every row and every column. Then every matrix in $\mathcal{P}$
is a convex combination of matrices $\{P^\pi : \pi\in\mathfrak{S}_n\}$. Furthermore, permutation matrices are the vertices
of the simplex $\mathcal P$.
\end{theorem}

In view of this result, the combinatorial optimization problem (\ref{discrete}) is equivalent to
the following problem of continuous optimization:
\begin{equation}\label{cont}
P^{\textrm{LSS}} = \arg\min_{P \in\mathcal P} \mathbf{tr} \big(M P\big),
\end{equation}
in the sense that $\pi$ is a solution to (\ref{discrete}) if and only if $P^\pi$ is a solution to (\ref{cont}).
To prove this claim, let us remark that for every $P\in\mathcal{P}$, there exist coefficients $\alpha_1,\ldots,\alpha_{n!}\in[0,1]$ such that
$P = \sum_{i=1}^{n!} \alpha_i P^{\pi_i}$ and $\sum_{i=1}^{n!} \alpha_i = 1$.
Therefore, we have $\mathbf{tr}\big(MP\big) = \sum_{i=1}^{n!} \alpha_i \mathbf{tr}\big(MP^{\pi_i}\big) \ge \min_{\pi\in\mathfrak{S}_n} \mathbf{tr}\big(MP^{\pi}\big)$ and $\mathbf{tr}\big(MP^{\textrm{LSS}}\big)\ge \mathbf{tr}\big(MP^{\pi^{\textrm{LSS}}}\big).$

The great advantage of (\ref{cont}) is that it concerns the minimization of a linear function under linear constraints and, therefore,
is a problem of linear programming that can be efficiently solved even for large values of $n$. The same arguments apply to the estimators $\pi^{\textrm{LSNS}}$ and $\pi^{\textrm{LSL}}$, only the matrix $M$ needs to be changed.

There is a second way to compute the estimators LSS, LSNS and LSL efficiently. Indeed, the computation of the aforementioned maximum likelihood estimators is a particular case of the assignment problem, which consists in finding a minimum weight matching in a weighted bipartite graph, where the matrix of the costs is the matrix $M$ from above. This means that the cost of assigning the $i^{\text{th}}$ feature of the first image to the $j^{\text{th}}$ feature of the second image is either
\begin{itemize}\setlength{\itemsep}{0pt}
\item the squared distance  ${\|X_i-X^{\text{\tt\#}}_j\|^2}$,
\item or the normalized squared distance $\|X_i-X^{\text{\tt\#}}_j\|^2/(\sigma_i^2+(\s^{\text{\tt\#}}_j)^2)$,
\item or the logarithm of the squared distance $\log \|X_i-X^{\text{\tt\#}}_j\|^2$.
\end{itemize}
The so-called Hungarian algorithm presented in \citet*{Kuhn1955} solves the assignment problem in time $O(n^3)$.

\section{Extensions}\label{section:extensions}

In this section, we briefly discuss possible extensions of the foregoing results to
other distances, more general matching criteria and to the estimation of an arrangement.

\subsection{Minimax Rates for the Hamming Distance}
In the previous sections, the minimax rates were obtained for the error of estimation measured by the risk $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}(\hat\pi\neq\pi^*)
=\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}[\delta_{\text{0\,-1}}(\hat\pi,\pi^*)]$, which may be considered as too restrictive. Indeed, one could find acceptable an estimate
having a small number of mismatches, if it makes it possible to significantly reduce the perceivable distance of separation. These considerations lead
to investigating the behavior of the estimators in terms of the Hamming loss, \textit{i.e.}, to studying the risk
\begin{equation*}
\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}[\delta_{H}(\hat\pi,\pi^*)]=\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}
\Big[\frac1n \sum\limits_{i=1}^n \mathds{1}_{\{\hat\pi(i)\neq\pi^*(i)\}}\Big]
\end{equation*}
corresponding to the expected average number of mismatched features. Another advantage of studying the Hamming loss instead of
the 0\,-1 loss is that
the former sharpens the difference between the performances of various estimators. Note, however, that thanks to the inequality
$\delta_{H}(\hat\pi,\pi^*)\le \delta_{\text{0\,-1}}(\hat\pi,\pi^*)$, all the upper bounds established for the minimax rate of separation
under the 0-1 loss directly carry over to the case of the Hamming loss.
This translates into the following theorem.

\begin{theorem}\label{upperbound4}
Let $\alpha\in(0,1)$ and condition (\ref{sigma_h}) be fulfilled. If $\hat\pi$ is either $\pi^{\textrm{LSNS}}$
or $\pi^{\textrm{LSL}}$,
then $\bar\kappa_\alpha(\hat\pi) \leq 4\max\big\{\big(2\log \frac{8n^2}{\alpha}\big)^{1/2}, \big(d\log\frac{4n^2}{\alpha}\big)^{1/4} \big\}$. That is,
if
$$\kappa=4\max\Big\{\Big(2\log \frac{8n^2}{\alpha}\Big)^{1/2}, \Big(d\log\frac{4n^2}{\alpha}\Big)^{1/4} \Big\}
$$
and $\bar\Theta_\kappa$ is the set of all $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ such that $\bar\kappa(\boldsymbol{\theta},\sigma)\ge \kappa$, then
\begin{equation*}
\max_{\pi^*\in\mathfrak{S}_n}\sup_{\boldsymbol{\theta}\in\bar\Theta_\kappa}\mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}[\delta_{H}(\hat\pi,\pi^*)] \le \alpha.
\end{equation*}
\end{theorem}

\begin{table}
\begin{center}
\begin{tabular}{r|rrrrrrrrr}
\toprule
$n=$ & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12\\
$M_n \ge $ & 19 & 57 & 179 & 594 & 1939 & 3441 & 11680 & 39520 & 86575
\\
$\frac{\log M_n}{n\log n}\ge$  & 0.53 & 0.50 & 0.48 & 0.47 & 0.455 & 0.412 & 0.407 & 0.401 & 0.381
\\
$\frac{\log M_n}{n\log n}\le$  & 0.53 & 0.50 & 0.48 & 0.47 & 0.455 & 0.445 & 0.436 & 0.427 & 0.420
\\
\bottomrule
\end{tabular}
\caption{The values of $M_n=\mathcal M(1/4, B_{2,n}(2),\delta_H)$ for $n\in\{4,\ldots,12\}$.
The lower bound is just the cardinality of one $\epsilon$-packing, not necessarily the largest one.
The upper bound is merely the cardinality of the $\ell_2$-ball.
\label{tab:1}}
\end{center}
\end{table}
While this upper bound is an immediate consequence of Theorem~\ref{upperbound2}, getting lower bounds for the Hamming loss
appears to be more difficult. To state the corresponding result, let is consider the case of homoscedastic noise and introduce some
notation. We denote by $\delta_2(\cdot,\cdot)$ the normalized $\ell^2$-distance on the space of permutations $\mathfrak{S}_n$:
$\delta_2(\pi,\pi')^2=\frac1n\sum_{k=1}^n\big(\pi(k)-\pi'(k)\big)^2$. Let $B_{2,n}(R)$ be the ball of $(\mathfrak{S}_n,\delta_2)$
with radius $R$ centered at $\textit{id}$. As usual, we denote by $\mathcal M(\epsilon, B_{2,n}(R),\delta_H)$ the $\epsilon$-packing number
of the $\ell_2$-ball $B_{2,n}(R)$ in the metric $\delta_H$. This means that $\mathcal M(\epsilon, B_{2,n}(R),\delta_H)$ is the
largest integer $M$  such that there exist permutations $\pi_1,\ldots,\pi_M\in B_{2,n}(R)$ satisfying $\delta_H(\pi_i,\pi_j)\ge \epsilon$
for every $i\not=j$. One can easily check that replacing $B_{2,n}(R)$ by any other ball of radius $R$ leaves the packing number
$\mathcal M(\epsilon, B_{2,n}(R),\delta_H)$ unchanged. We set $M_n=\mathcal M(1/4, B_{2,n}(2),\delta_H)$.

\begin{theorem}\label{lowerbound4}
Let $\sigma_k=\sigma$ for all $k\in\{1,\ldots,n\}$ and $\bar\Theta_\kappa$ is the set of all $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ such that $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\ge \kappa$. Furthermore, assume that one of the following two conditions is fulfilled:
\begin{itemize}
\item $n\ge 3$ and $\kappa = (1/4)\big(\frac{\log M_n}{n}\big)^{1/2}$,
\item $n\ge 26$, $d\ge 24\log n$ and $\kappa\le (1/8)(d\log n)^{1/4}$.
\end{itemize}
Then
$\inf_{\hat\pi} \max_{\pi^*\in\mathfrak{S}_n} \sup_{\boldsymbol{\theta}\in\bar\Theta_\kappa} \mathbf E_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*} [\delta_H(\hat\pi,\pi^*)] >  2.15\%$.
\end{theorem}

This result implies that in the regime of moderately large dimension, $d\ge 24\log n$, the minimax rate of the separation is the same as
the one under the 0-1 loss and it is achieved by the LSL estimator. The picture is less clear in the regime of small dimensions, $d=o(\log n)$.
If one proves that for some $c>0$, the inequality $\log M_n\ge cn\log n$ holds for every $n\ge 3$, then the lower bound of the last theorem
matches the upper bound of Theorem~\ref{upperbound4} up to constant factors and leads to the minimax rate of separation
$\max\{(\log n)^{1/2},(d\log n)^{1/4}\}$. Unfortunately, we were unable to find any result on the order of magnitude of $\log M_n$,
therefore,  we cannot claim that there is no gap between our lower bound and the upper one. However, we did a small experiment for
evaluating $M_n$ for small values of $n$. The result is reported in Table~\ref{tab:1}.

\subsection{More General Matching Criteria}\label{subsection5.2}

In the previous sections, we were considering two vectors $\theta_i$ and ${\theta}^\text{\tt\#}_j$ as matching if $\theta_i\equiv{\theta}^\text{\tt\#}_j$, and $\equiv$
was the usual equality. In this part, we show that our results can be extended to more general matching criteria, defined as follows.
Let $A$, ${A}^{\texttt{\#}}$ be two known $p\times d$ matrices with some $p\in\mathbb N$ and $b$, ${b}^{\texttt{\#}}$ be two known vectors from $\mathbb R^{d}$.
We write $\theta\equiv_{A,b} {\theta}^\text{\tt\#}$, if
\begin{equation}\label{A,b}
A(\theta-b)={A}^{\texttt{\#}}({\theta}^\text{\tt\#} -{b}^{\texttt{\#}}).
\end{equation}
Note that the case of equality studied in previous sections is obtained for $A ={A}^{\texttt{\#}}=\mathbf I_d$ and $b={b}^{\texttt{\#}}=0$, where $\mathbf I_d$ is
the identity matrix of size $d$. Let us first note that without loss of generality, by a simple transformation of the features, one can replace
(\ref{A,b}) by the simpler relation
\begin{equation}\label{A'}
\bar\theta=B \bar\t^\text{\tt\#},
\end{equation}
where $\bar\theta\in\mathbb R^{d_1}$ for $d_1 = \text{rank}(A)$, $\bar\t^\text{\tt\#}\in\mathbb R^{d_2}$ for $d_2 = \text{rank}({A}^{\texttt{\#}})$  and
$B$ is a $d_1\times d_2$ known matrix. Indeed, let $A = U^\top\Lambda V$ (resp.\ ${A}^{\texttt{\#}} = \tilde U^\top \tilde\Lambda\tilde V$)
be the singular value decomposition of $A$ (resp.\ ${A}^{\texttt{\#}}$), with orthogonal matrices $U\in\mathbb R^{d_1\times p}$, $V\in\mathbb R^{d_1\times d}$
and a diagonal matrix $\Lambda\in\mathbb R^{d_1\times d_1}$ with positive entries. Then, one can deduce (\ref{A'}) from (\ref{A,b}) by setting
$\bar\theta = V(\theta-b)$, $\bar\t^\text{\tt\#}=\tilde V({\theta}^\text{\tt\#}-{b}^{\texttt{\#}})$ and $B = \Lambda^{-1}U\tilde U^\top\tilde\Lambda$. Of course, the same
transformation should be applied to the observed noisy features, which leads to $\bar X_i= V(X_i-b)$ and
$\bar X^\text{\tt\#}_i=\tilde V(X^{\text{\tt\#}}-{b}^{\texttt{\#}})$. Since $V$ and $\tilde V$ are orthogonal matrices, \textit{i.e.}, satisfy the relations
$VV^\top = \mathbf I_{d_1}$ and $\tilde V\tilde V^\top = \mathbf I_{d_2}$, the noise component in the transformed noisy features is still white
Gaussian.

All the four estimators introduced in Section~\ref{section:estimators} can be adapted to deal with such type of
criterion. For example, denoting by $M$ the matrix $B(B^\top B)^+B^\top+BB^\top$ where $M^+$ is the
Moore-Penrose pseudoinverse of the matrix $M$, the LSL estimator should be modified as follows
\begin{equation}
\pi^{\textrm{LSL}}  \triangleq \arg\min_{\pi\in\mathfrak{S}_n} \sum_{i=1}^n \log \|M^+(\bar X_{\pi(i)}-B\bar X^\text{\tt\#}_i)\|^2. \label{LSL1}
\end{equation}
All the results presented before can be readily extended to this case. In particular, if we assume that all the nonzero singular values of
$B$ are bounded and bounded away from 0 and $\text{rank}(B)=q$, then the minimax rate of separation is given by
$\max\big\{(\log n)^{1/2}, (q \log n)^{1/4}\big\}$. This rate is achieved, for instance, by the LSL estimator.

Let us briefly mention two situations in which this kind of general affine criterion may be useful. First, if each feature $\theta$ (resp.\
${\theta}^\text{\tt\#}$) corresponds to a patch in an image $I$ (resp.\ $I^{\text{\tt\#}}$), then for detecting pairs of patches that match each other
it is often useful to neglect the changes in illumination. This may be achieved by means of the criterion $A\theta=A{\theta}^\text{\tt\#}$
with $A =  \mathbf I_d-\frac1{d}\mathbf 1_d\mathbf 1_d^\top$. Indeed, the multiplication of the vector $\theta$ by $A$ corresponds
to removing from pixel intensities the mean pixel intensity of the patch. This makes the feature invariant by change of illumination.
The method described above applies to this case and the corresponding rate is $\max\big\{(\log n)^{1/2}, ((d-1)\log n)^{1/4}\big\}$
since the matrix $A$ is of rank $d-1$. Second, consider the case when each feature
combines the local descriptor of an image and the location in the image at which this descriptor is computed. If we have at our
disposal an estimator of the transformation that links the two images and if this transformation is linear, then we are in the
aforementioned framework. For instance, let each $\theta$ be composed of a local descriptor $\mathbf d\in\mathbb R^{d-2}$ and its location
$\mathbf x\in\mathbb R^2$. Assume that the first image $I$ from which the features $\theta_i=[\mathbf d_i,\mathbf x_i]$ are extracted is
obtained from the image $I^{\text{\tt\#}}$ by a rotation of the plane. Let $R$ be an estimator of this rotation and
${\theta}^\text{\tt\#}_i=[\mathbf{d}^\text{\tt\#}_i,\mathbf{x}^\text{\tt\#}_i]$ be the features extracted from the image $I^{\text{\tt\#}}$. Then, the aim is to find the permutation $\pi$
such that $\mathbf{d}^\text{\tt\#}_i=\mathbf d_{\pi(i)}$ and $\mathbf{x}^\text{\tt\#}_i=R\mathbf x_{\pi(i)}$ for every $i=1,\ldots,n$. This corresponds to taking in (\ref{A'})
the matrix $B$ given by
$$
B =\begin{pmatrix}
\mathbf I_{d-2} & 0\\
0 & R
\end{pmatrix}.
$$
This matrix $B$ being orthogonal, the resulting minimax rate of separation is exactly the same as when $B = \mathbf I_d$.

\begin{remark}
An interesting avenue for future research concerns the determination of the minimax rates in the case when the
equivalence of two features is understood under some transformation $A$ which is not completely determined.
For instance, one may consider that the features $\theta$ and $\theta'$ match if there is a matrix $A$
in a given parametric family $\{A_\tau:\tau\in\mathbb R\}\subset \mathbb R^{d\times d}$ for which $\theta = A\theta'$. In
other terms, $\theta\equiv\theta'$ is understood as $\inf_\tau \|\theta-A_\tau\theta'\|=0$. Such a criterion of
matching may be useful for handling various types of invariance (see \cite{CollierDalalyan2015,Collier2012} for
invariance by translation).
\end{remark}

\subsection{Estimation of an Arrangement}\label{sec:arrangement}

An interesting extension concerns the case of the estimation of a general arrangement, \textit{i.e., } the case when $m$
and $n$ are not necessarily identical. In such a situation, without loss of generality, one can assume that
$n\le m$ and look for an injective function $$\pi^*:\{1,\ldots,n\}\to\{1,\ldots,m\}.$$ All the estimators presented
in Section~\ref{section:estimators} admit natural counterparts in this rectangular setting. Furthermore, the
computational method using the Birkhoff-von Neumann theorem is still valid in this setting, and is justified
by the extension of the Birkhoff-von Neumann theorem recently proved by \citet{BudishCheKojimaMilgrom2009}.
In this case, the minimization should be carried out over the set of all matrices $P$ of size $(n,m)$ such
that $P_{i,j}\ge 0,$ and
\begin{equation*}
\begin{cases}
\ \sum_{i=1}^n P_{i,j}\le 1 \\
\ \sum_{j=1}^m P_{i,j}= 1
\end{cases},\quad (i,j)\in\{1,\ldots,n\}\times\{1,\ldots,m\}.
\end{equation*}
From a practical point of view, it is also important to consider the issue of
robustness with respect to the presence of outliers, \textit{i.e., } when for some $i$ there is no
$X^{\text{\tt\#}}_j$ matching with $X_i$. The detailed exploration of this problem being out
of scope of the present paper, let us just underline that the LSL-estimator seems
to be well suited for such a situation because of the robustness of the
logarithmic function. Indeed, the correct matches are strongly rewarded because
$\log(0)=-\infty$ and the outliers do not interfere too much with the estimation of
the arrangement thanks to the slow growth of $\log$ in ${+\infty}$.

In the rectangular case of $n \le m$ we modify the estimator $\pi^{\text{LSNS}}$ from \eqref{LSNS} as follows
\begin{equation}
\bar\pi^{\textrm{LSNS}} \triangleq \arg\min_{\pi:S_1\to S_2}\sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2},
\end{equation}
where the minimum is taken over all injective functions $\pi$ with $S_1 = \{1, \dots, n\}$ and $S_2 = \{1, \dots, m\}$. The analogue assumption of \eqref{sigma_h} reads as

\begin{align}\label{noise-match}
\forall i \in \{1, \dots, n\}, ~~ \sigma_i = \s^{\text{\tt\#}}_{\pi^*(i)}.
\end{align}

\begin{theorem}\label{upperbound_outliers}
Let $\alpha \in (0,1)$ and condition \eqref{noise-match} be fulfilled. If $\hat{\pi} \equiv \bar\pi^{\text{LSNS}}$, then
\begin{align*}
\bar\kappa_\alpha(\hat\pi) \leq 4\max\Big\{\Big(2\log \frac{8nm}{\alpha}\Big)^{1/2}, \Big(d\log\frac{4nm}{\alpha}\Big)^{1/4} \Big\}.
\end{align*}
\end{theorem}

\section{Proofs of the Theorems}\label{section:theorems}

In this section we collect the proofs of the theorems. We start with the proof of Theorem~\ref{upperbound2}, since it concerns the more general setting and the proof of Theorem~\ref{upperbound1} can be deduced from that of Theorem~\ref{upperbound2} by simple arguments. We then prove the other theorems in the usual order and postpone the proofs of some technical lemmas to the next section.

\begin{proofof}[Proof of Theorem~\ref{upperbound2}]
To ease notation and without loss of generality, we assume that $\pi^*$ is the identity permutation denoted by $\textit{id}$. Furthermore, since there is no risk of confusion, we write $\mathbf P$ instead of $\mathbf P_{\boldsymbol{\theta},\boldsymbol{\sigma},\pi^*}$. We wish to bound the probability of the event $\Omega = \{\hat\pi\neq \textit{id}\}$.

Let us first denote by $\hat\pi$ the maximum likelihood estimator $\pi^{\textrm{LSL}}$ defined by (\ref{LSL}). We have
\begin{equation*}
\Omega \subset \bigcup_{\pi\not=\textit{id}} \Omega_\pi,
\end{equation*}
where
\begin{align*}
\Omega_\pi &
		= \Big\{\sum_{i=1}^n \log \frac{\|X_i-X^{\text{\tt\#}}_i\|^2}{\|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2}\ge 0\Big\}
		= \Big\{\sum_{i : \pi(i)\not= i} \log \frac{\|X_i-X^{\text{\tt\#}}_i\|^2}{\|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2}\ge 0\Big\}.
\end{align*}
On the one hand, for every permutation $\pi$,
\begin{align*}
\sum_{\pi(i)\not =i} \log\Big({\textstyle\frac{2\sigma_i^2}{\sigma_i^2+\sigma_{\pi(i)}^2}}\Big)
 &= \sum_{i=1}^n \big(\log(2\sigma_i^2) - \log(\sigma_i^2+\sigma_{\pi(i)}^2) \big)\\
 &=  \sum_{i=1}^n \frac{\log(2\sigma_i^2)+\log(2\sigma_{\pi(i)}^2)}2-\log(\sigma_i^2+\sigma_{\pi(i)}^2)
\end{align*}
so, using the concavity of the logarithm, this quantity is nonpositive. Therefore,
\begin{align*}
\Omega_\pi
    &\subset  \Big\{\sum_{i : \pi(i)\not= i}\!\! \log \textstyle\frac{\|X_i-X^{\text{\tt\#}}_i\|^2/(2\sigma_i^2)}
        {\|X_{\pi(i)}-X^{\text{\tt\#}}_i\|^2/(\sigma_i^2+\sigma_{\pi(i)}^2)}\ge 0\Big\}\\
    &\subset \bigcup_{i=1}^n \bigcup_{j\not =i}\Big\{\frac{\|X_i-X^{\text{\tt\#}}_i\|^2}{2\sigma_i^2}\ge\frac{\|X_{j}-X^{\text{\tt\#}}_i\|^2}{\sigma_j^2+\sigma_{i}^2}\Big\}.
\end{align*}
This readily yields $\Omega\subset  \bar\Omega$, where
\begin{align}\label{inclusion}
\bar\Omega =\bigcup_{i=1}^n \bigcup_{j\not =i}\Big\{\frac{\|X_i-X^{\text{\tt\#}}_i\|^2}{2\sigma_i^2}\ge \frac{\|X_{j}-X^{\text{\tt\#}}_i\|^2}{\sigma_j^2+\sigma_{i}^2}\Big\}.
\end{align}
Furthermore, the same inclusion is true for the LSNS estimator as well. Therefore, the rest of the proof
is common for the estimators LSNS and LSL.

We set $\sigma_{i,j}=(\sigma_i^2+\sigma_j^2)^{1/2}$ and
\begin{align*}
\zeta_1 &= \max_{i\not =j} \bigg|
\frac{(\theta_{i}-\theta_j)^\top(\sigma_i\xi_i-\sigma_j\xi^\text{\tt\#}_j)}{\|\theta_i-\theta_j\|\sigma_{i,j}}\bigg|, \quad
\zeta_2 = d^{-1/2} \max_{i,j }  \bigg|\,\bigg\|\frac{\sigma_i\xi_i-\sigma_j\xi^\text{\tt\#}_j}{\sigma_{i,j}}\bigg\|^2 - d\bigg|.
\end{align*}
Since $\pi^*=\textit{id}$, it holds that for every $i\in\{1,\ldots,n\}$,
\begin{align*}
\|X_i-X^{\text{\tt\#}}_i\|^2 &=  \sigma_i^2 \|\xi_i-\xi^\text{\tt\#}_i\|^2 \le 2\sigma_i^2(d+\sqrt{d}\zeta_2).
\end{align*}
Similarly, for every $j\not =i$,
\begin{align*}
\|X_{j}-X^{\text{\tt\#}}_i\|^2
    &=  \|\theta_{j}-\theta_i\|^2+ \|\sigma_{j}\xi_{j}-\sigma_i\xi^\text{\tt\#}_i\|^2+ 2(\theta_{j}-\theta_i)^\top(\sigma_{j}\xi_{j}-\sigma_i\xi^\text{\tt\#}_i).
\end{align*}
Therefore,
\begin{align*}
\|X_{j}-X^{\text{\tt\#}}_i\|^2
    &\ge  \|\theta_{j}-\theta_i\|^2+\sigma_{i,j}^2 (d-\sqrt{d}\zeta_2) - 2\|\theta_{j}-\theta_i\|\sigma_{i,j}\zeta_1.
\end{align*}
This implies that on the event $\Omega_1 = \{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\ge \zeta_1\}$  it holds that
\begin{align*}
\frac{\|X_{j}-X^{\text{\tt\#}}_i\|^2}{\sigma_{i,j}^2}
    \ge  \bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2 - 2\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\zeta_1+d-\sqrt{d}\zeta_2.
\end{align*}
Combining these bounds, we get that
$$
\Omega\cap \Omega_1\subset  \Big\{d+\sqrt{d}\zeta_2\ge\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2 - 2\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\zeta_1+d-\sqrt{d}\zeta_2\Big\},
$$
which implies that
\begin{align}
\mathbf P(\Omega)
        &\le \mathbf P(\Omega_1^\complement) + \mathbf P\big(\Omega\cap\Omega_1\big)\nonumber\\
        &\le \mathbf P\big(\zeta_1\ge \bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
            +\mathbf P(2\sqrt{d}\zeta_2+2\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\zeta_1\ge\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2)\nonumber\\
        &\le 2\mathbf P\big(\zeta_1\ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})}4\big)+
            \mathbf P\big(\zeta_2\ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}{4\sqrt{d}}\big).\label{ineq:2}
\end{align}
Finally, one easily checks that for suitably chosen random variables $\zeta_{i,j}$ drawn from the standard Gaussian distribution,
it holds that $\zeta_1 = \max_{i\not=j}|\zeta_{i,j}|$. Therefore, using the well-known tail bound for the standard Gaussian
distribution in conjunction with the union bound, we get
\begin{align}
\mathbf P\big(\zeta_1\ge \textstyle\frac14\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
    &\le \sum_{i\not =j}\nolimits \mathbf P\big(|\zeta_{i,j}|\ge \textstyle\frac14\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
    \le 2 n^2 e^{-\frac1{32}\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}.\label{ineq:3}
\end{align}
To bound the large deviations of the random variable $\zeta_2$, we rely on the following result.
\begin{lemma}[\cite{LaurentMassart2000}, Eq. (4.3) and (4.4)]\label{concentration}
If $Y$ is drawn from the chi-squared distribution  $\chi^2(D)$, where $D\in\mathbb N^*$, then, for every $x>0$,
\begin{equation*}
\begin{cases}
\ \mathbf P\big(Y-D \le -2\sqrt{Dx} \big) \le e^{-x}, \phantom{\Big()}\\
\ \mathbf P\big(Y-D \ge 2\sqrt{Dx} + 2x \big) \le e^{-x}. \phantom{\Big()}
\end{cases}
\end{equation*}
As a consequence, $\forall y>0$,
$\mathbf P\big(D^{-1/2}|Y-D| \ge y\big) \le 2\exp\big\{-\textstyle\frac18 y(y\wedge \sqrt{D})\big\}$.
\end{lemma}
This inequality, combined with the union bound, yields
\begin{align}
\mathbf P\Big(\zeta_2 \ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}{4\sqrt{d}}\Big)
        &\le 2n^2 \exp\Big\{-\frac{(\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})/16)^2}{d}(\bar\kappa^2(\boldsymbol{\theta},\boldsymbol{\sigma})\wedge 8d)\Big\}.\label{ineq:4}
\end{align}
Combining inequalities (\ref{ineq:2})-(\ref{ineq:4}), we obtain that as soon as
\begin{align*}
\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\ge 4\Big(\sqrt{2\log(8n^2/\alpha)} \vee \big(d\log(4n^2/\alpha)\big)^{1/4}\Big),
\end{align*}
we have $\mathbf P(\hat\pi\not = \pi^*)=\mathbf P(\Omega)\le \alpha$.
\vspace{-20pt}

~
\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{upperbound_outliers}]
The proof is a simple extension of Theorem \ref{upperbound2}. We start by assuming that $\pi^* \equiv id$, where $\forall i \in \{1, \dots, n\}, ~ id(i)=i$. We wish to bound the probability of the event $\Omega = \{ \hat\pi \neq id \}$, where $\hat\pi = \bar\pi^{\text{LSNS}}$. It is evident that

\begin{align*}
\Omega \subset \bigcup_{\pi\neq\textit{id}} \Omega_{\pi},
\end{align*}
where
\begin{align*}
\Omega_{\pi} &= \bigg\{ \sum_{i=1}^n \frac{\| X_i - X^{\text{\tt\#}}_i \|^2}{2\sigma_i^2} \ge \sum_{i=1}^n \frac{\| X_i - X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_i^2 + (\s^{\text{\tt\#}}_{\pi(i)})^2 }\bigg\} \\ &\subset
\bigcup_{i=1}^n \bigcup_{j \in S_2 \setminus \{ i \}} \bigg\{ \frac{\| X_i - X^{\text{\tt\#}}_i \|^2}{2\sigma_i^2} \ge \frac{\| X_i - X^{\text{\tt\#}}_{j}\|^2}{\sigma_i^2 + (\s^{\text{\tt\#}}_{j})^2 }\bigg\}.
\end{align*}
Denote $\sigma_{i, j} = (\sigma_i^2 + (\s^{\text{\tt\#}}_j)^2)^{1/2}$ and
\begin{align*}
\zeta_1 = \max_{i \neq j} \bigg| \frac{(\theta_i - {\theta}^\text{\tt\#}_j)^\top (\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j)}{\|\theta_i - {\theta}^\text{\tt\#}_j\|\sigma_{i, j}}\bigg|, \quad \zeta_2 = d^{-1/2} \max_{i, j} \bigg| \bigg\| \frac{\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j}{\sigma_{i,j}} \bigg\|^2 - d \bigg|.
\end{align*}
Since $\pi^* = \textit{id}$, it holds for every $i \in \{1, \dots, n\}$,
\begin{align*}
\| X_i - X^{\text{\tt\#}}_i \|^2 = \sigma_i^2 \| \xi_i - \xi^\text{\tt\#}_i\|^2 \le 2\sigma_i^2 (d + \sqrt{d} \zeta_2).
\end{align*}
Similarly, for every $j \in \{1, \dots, m\}$ and $j \neq i$,
\begin{align*}
\| X_i - X^{\text{\tt\#}}_j \|^2 &= \| \theta_i - {\theta}^\text{\tt\#}_j \|^2 + \| \sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j \|^2 + 2(\theta_i - {\theta}^\text{\tt\#}_j)^\top (\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j) \\ &\ge
\| \theta_i - {\theta}^\text{\tt\#}_j \|^2 + \sigma_{i,j}^2 (d - \sqrt{d}\zeta_2) - 2 \| \theta_i - {\theta}^\text{\tt\#}_j \| \sigma_{i,j} \zeta_1.
\end{align*}
On the event $\Omega_1 = \{ \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \ge \zeta_1\}$ from the previous inequality we have

\begin{align*}
\frac{\| X_i - X^{\text{\tt\#}}_j \|^2}{\sigma_{i,j}^2} \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 - 2\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \zeta_1 + d - \sqrt{d} \zeta_2,
\end{align*}
where $\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) = \min_{i \neq j} \frac{\| \theta_i - {\theta}^\text{\tt\#}_j\|}{\sigma_{i,j}}$

\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{upperbound1}]
Without loss of generality, we assume $\pi^* =\textit{id}$. It holds that, on the event
$$
\mathcal A= \bigcap_{i=1}^n\bigcap_{j\not =i}\Big\{ \|X_{i}-X^{\text{\tt\#}}_{i}\|< \|X_{j}-X^{\text{\tt\#}}_{i}\|\Big\},
$$
all the four estimators coincide with the true permutation $\textit{id}$. Therefore, we have
$$
\{\hat \pi\not=\textit{id}\} \subseteq \bigcup_{i=1}^n\bigcup_{j\not =i}\Big\{ \|X_{i}-X^{\text{\tt\#}}_{i}\|\geq  \|X_{j}-X^{\text{\tt\#}}_{i}\|\Big\}.
$$
The latter event is included in $\bar\Omega$ at the right-hand side of (\ref{inclusion}), the probability of which has been
already shown to be small in the previous proof.
\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{lowerbound1}]
We refer the reader to the proof of Theorem~\ref{lowerbound2} below, which concerns the more general situation. Indeed,
when all the variances $\sigma_j$ are equal, Theorem~\ref{lowerbound2} boils down to
Theorem~\ref{lowerbound1}.
\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{lowerbound2}]
To establish lower bounds for various types of risks we will  use the following lemma:
\begin{lemma}[\citet{Tsybakov2009}, Theorem 2.5]\label{Tsybakov}
Assume that for some integer $M\ge 2$ there exist distinct permutations $\pi_0,\ldots,\pi_M\in\mathfrak{S}_n$ and mutually absolutely continuous
probability measures $\mathbf Q_0,\ldots,\mathbf Q_M$ defined on a common probability space $(\mathcal Z,\mathscr Z)$ such that
\begin{equation*}
 \frac{1}{M} \sum_{j=1}^M K(\mathbf Q_j,\mathbf Q_0) \le \frac{1}{8} \log M.
\end{equation*}
Then, for every measurable mapping $\tilde\pi:\mathcal Z\to\mathfrak{S}_n$,
\begin{equation*}
\max_{j=0,\ldots,M} \mathbf Q_j(\tilde\pi\neq\pi_j) \ge \frac{\sqrt{M}}{\sqrt{M}+1} \Big( \frac{3}{4}-\frac{1}{2\sqrt{\log M}} \Big).
\end{equation*}
\end{lemma}
To prove Theorem~\ref{lowerbound2}, we split the inequality $8\kappa\le \max\{(\log n)^{1/2},(cd\log n)^{1/4}\}$ into two cases
\begin{align*}
&\textbf{Case 1:}\qquad 8\kappa \le  (\log n)^{1/2},\\
&\textbf{Case 2:}\qquad (\log n)^{1/2}\le 8\kappa \le (cd\log n)^{1/4}.
\end{align*}
\paragraph{\textbf{Case 1:} We assume that $8\kappa \le (\log n)^{1/2}$. }
Denote by $m$ the largest integer such that $2m\le n$. We assume without loss of generality that the noise levels are ranked in increasing order: $\sigma_1\leq\ldots\leq\sigma_n.$
Then, we construct a least favorable set of vectors for the estimation of the permutation. To ease notation, we set $\sigma_{i,j}=(\sigma_i^2+\sigma_j^2)^{1/2}$.
\begin{lemma}\label{construction_btheta}
Assume that $m$ is the largest integer such that $2m\le n$. Then there is a set of vectors $\boldsymbol{\theta}$ such that
\begin{equation*}
\frac{\|\theta_1-\theta_2\|}{\sigma_{1,2}} = \ldots = \frac{\|\theta_{2m-1}-\theta_{2m}\|}{\sigma_{2m-1,2m}} = \kappa,
\end{equation*}
and for every pair $\{i,j\}$ different from the pairs $\{1,2\},\ldots,\{2m-1,2m\} $ we have
\begin{equation*}
\frac{\|\theta_i-\theta_j\|}{\sigma_{i,j}} > \kappa\bigg(1+\frac{\max_{1\le\ell\le n} \sigma_\ell}{\min_{1\le\ell\le n} \sigma_\ell}\bigg).
\end{equation*}
\end{lemma}

Let $\boldsymbol{\theta}^0$ be the set constructed in Lemma~\ref{construction_btheta}. The latter implies that $\boldsymbol{\theta}^0$ belongs to $\bar\Theta_\kappa$, so that,
denoting for every $k\in\{1,\ldots,m\}, \pi_k=(2k-1~2k)$ the transposition of $\mathfrak{S}_n$ that only permutes $2k-1$ and $2k$, and $\pi_0=\textit{id},$
we get the following lower bound for the risk:
\begin{equation*}
\inf_{\hat\pi} \sup_{(\pi,\boldsymbol{\theta})\in\mathfrak{S}_n\times\bar\Theta_\kappa} \mathbf P_{\boldsymbol{\theta},\pi} (\hat\pi\neq\pi) \ge \inf_{\hat\pi} \max_{k=0,\ldots,m} \mathbf P_{\boldsymbol{\theta}^0,\pi_k}(\hat\pi\neq\pi_k).
\end{equation*}
In order to use Lemma~\ref{Tsybakov} with $\mathbf Q_j=\mathbf P_{\boldsymbol{\theta}^0,\pi_j}$, we compute for every $k\in\{1,\ldots,m\}$
\begin{align*}
K(\mathbf P_{\boldsymbol{\theta}^0,\pi_k},\mathbf P_{\boldsymbol{\theta}^0,\pi_0}) = &\|\theta_{2k-1}^0-\theta_{2k}^0\|^2\Big(\frac{1}{2\sigma^2_{2k}} + \frac{1}{2\sigma^2_{2k-1}} \Big)+
\frac{d}2\Big(\frac{\sigma^2_{2k-1}}{\sigma^2_{2k}}+\frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-2\Big) .
\end{align*}
Using the fact that $\|\theta_{2k-1}^0-\theta_{2k}^0\|^2 = \kappa^2(\sigma_{2k-1}^2+\sigma_{2k}^2)$, we get
\begin{align*}
 K(\mathbf P_{\boldsymbol{\theta}^k,\pi_k},\mathbf P_{\boldsymbol{\theta}^0,\pi_0})
=  \frac{\kappa^2}2\big(2+\sigma_{2k-1}^2\sigma_{2k}^{-2}+\sigma_{2k}^2\sigma_{2k-1}^{-2}\big) + \frac{d}2\Big(\frac{\sigma^2_{2k-1}}{\sigma^2_{2k}}+\frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-2\Big) .
\end{align*}
The inequality $\sigma_{2k-1}^2\le \sigma_{2k}^2$ implies the following two relations:
\begin{align*}
\frac{\sigma^2_{2k-1}}{\sigma^2_{2k}}+\frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-2 &\le \frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-1\\
\frac{\sigma^2_{2k-1}}{\sigma^2_{2k}}+\frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-2 &=\frac{\sigma_{2k-1}^2}{\sigma_{2k}^2}\Big(\frac{\sigma_{2k}^2}{\sigma_{2k-1}^2}-1\Big)^2
\le \Big(\frac{\sigma_{2k}^2}{\sigma_{2k-1}^2}-1\Big)^2.
\end{align*}
This readily yields $\frac{\sigma^2_{2k-1}}{\sigma^2_{2k}}+\frac{\sigma^2_{2k}}{\sigma^2_{2k-1}}-2 \le \frac{\log n}{16d}$
and, therefore,
\begin{align*}
K(\mathbf P_{\boldsymbol{\theta}^k,\pi_k},\mathbf P_{\boldsymbol{\theta}^0,\pi_0})\le  \kappa^2\Big(\frac32+\frac{\sigma_{2k}^2}{2\sigma_{2k-1}^2}\Big) + \frac{\log n}{32}.
\end{align*}
Next, we apply the following result.
\begin{lemma}
Let $a_1,a_2,\ldots,a_m$ be real numbers larger than one such that $\prod_{k=1}^m a_k\le A$. Then,
$\sum_{k=1}^m a_k \le m+\log A\max_k a_k$.
\end{lemma}
\begin{proof}
We use the simple inequality $e^x\le 1+xe^x$ for all $x\ge 0$. Replacing $x$ by $\log a_k$ and summing over
$k=1,\ldots,m$ we get
$$
\sum_{k=1}^m a_k \le \sum_{k=1}^m 1+ a_k\log a_k\le m+ \max_{k=1,\ldots,m} a_k \log \prod_{k=1}^m a_k \le m+ \max_{k=1,\ldots,m} a_k \log A.
$$
This completes the proof of the lemma.
\end{proof}
We apply this lemma to $a_k = \sigma_{2k}^2/\sigma_{2k-1}^2$. Since the variances are sorted in increasing order, we have
$\prod_{k=1}^m a_k \le \prod_{i=1}^{n-1} {\sigma_{i+1}^2}/{\sigma_i^2} = \sigma_n^2/\sigma_1^2\le 1+\gamma_{n,d}$ with
$\gamma_{n,d}=\max\big((\frac{\log n}{16d})^{1/2},\frac{\log n}{16d}\big)$. In conjunction
with the inequality $\log(1+x)\le x$, this entails that $\sum_{k=1}^m \sigma_{2k}^2/\sigma_{2k-1}^2\le m+ \gamma_{n,d}(1+\gamma_{n,d})$.
Then, since $\log n\le 1.8\log m$ for $n\ge 6$ and $\gamma_{n,d}\le  0.2\log n\le 0.36\log m$, we get
\begin{align*}
\frac1m\sum_{k=1}^m K(\mathbf P_{\boldsymbol{\theta}^k,\pi_k},\mathbf P_{\boldsymbol{\theta}^0,\pi_0})
&\le\kappa^2\bigg(\frac{3}{2}+\frac{1}{2m}\sum_{k=1}^m\frac{\sigma_{2k}^2}{\sigma_{2k-1}^2}\bigg) + \frac{9\log m}{160}\\
&\le\kappa^2\Big(2+\frac{\gamma_{n,d}}{2m}\big\{1+\gamma_{n,d}\big\}\Big) + \frac{9\log m}{160}\\
&\le\kappa^2\Big(2+\frac{0.18\log m}{m}\big\{1+0.36\log m\big\}\Big) + \frac{9\log m}{160}.
\end{align*}
Finally, using the fact that $m\ge 3$, we get $\frac1m\sum_{k=1}^m K(\mathbf P_{\boldsymbol{\theta}^k,\pi_k},\mathbf P_{\boldsymbol{\theta}^0,\pi_0}) \le
2.1\kappa^2+\frac{9\log m}{160} \le  \frac{\log m}{8}$ since $\kappa^2\le \frac1{64}\log n \le \frac{9}{320}\log m$.

We conclude by Lemma~\ref{Tsybakov} and by the monotonicity of the function $m\mapsto \frac{\sqrt{m}}{1+\sqrt{m}}(\frac34-\frac1{2\sqrt{\log m}})$ that
\begin{equation*}
\inf_{\hat\pi} \sup_{(\pi,\boldsymbol{\theta})\in\mathfrak{S}_n\times\bar\Theta_\kappa} \mathbf P_{\boldsymbol{\theta},\pi} (\hat\pi\neq\pi) \ge \frac{\sqrt{3}}{\sqrt{3}+1} \Big( \frac{3}{4} - \frac{1}{2\sqrt{\log3}} \Big)\ge 0.17.
\end{equation*}

\paragraph{\textbf{Case 2:} We assume that $(\log n)^{1/2}\le 8\kappa \le (cd\log n)^{1/4}$ with $c\le 1/20$.}
In this case, we have $d\ge\frac{1}{c}\log n$.  To get the desired result, we use Lemma~\ref{Tsybakov} for a
properly chosen family of probability measures described below.

\begin{lemma}\label{kullback_leibler_divergence_heteroscedastic}
Let $\epsilon_1,\ldots,\epsilon_n$ be real numbers defined by
$$
\epsilon_{{k}} = \sqrt{2/d}\;\kappa\sigma_{{k}},\qquad \forall{k}\in\{1,\ldots,n\},
$$
and let $\mu$ be the uniform distribution on $\mathcal E=\{\pm\epsilon_1\}^d\times\ldots\times\{\pm\epsilon_{n}\}^d$. We denote by $\mathbf P_{\mu,\pi}$
the probability measure on $\mathbb R^{d\times n}$ defined by $\mathbf P_{\mu,\pi}(A) = \int_{\mathcal E} \mathbf P_{\boldsymbol{\theta},\pi}(A)\,\mu(d\boldsymbol{\theta})$.
Assume that $\sigma_1\le\ldots\le \sigma_n$. For two positive integers $k<k'\le n$, set $\gamma = \frac{\sigma_{k'}^2}{\sigma_{k}^2}$ and let $\pi=({k}~{k}')$
be the transposition that only permutes ${k}$ and ${k}'$. Then
\begin{align*}
K(\mathbf P_{\mu,\pi},\mathbf P_{\mu,\textit{id}}) &\le  {4\kappa^2}\big(1-\gamma^{-1}\big)+\frac{8\kappa^4}{d}\big(2+\big(1+(2/d)\kappa^2\big)^2\gamma^{2}\big
+ \frac12\big(d+2\kappa^2\big)(\gamma-1)^2
\end{align*}
and $\mu(\mathcal E\backslash\bar\Theta_\kappa) \le ({n(n-1)}/{2})\ e^{-d/8}$.
\end{lemma}

The assumption on the noise levels entails that, for any integer $k\in\{1,\ldots,k'\}$,
$1\le \frac{\sigma_{k'}^2}{\sigma_{k}^2} \le 1+\frac14(\frac{\log n}{d})^{1/2}$,
and consequently, $(\gamma-1)^2=\big(\frac{\sigma_{k'}^2}{\sigma_{k}^2} -1\big)^2\le 4^{-2}\big(\frac{\log n}{d}\big)$.
Furthermore, $\frac{\kappa^2}{d}\le \frac{c}{64}\le \frac1{64}$ provided that $c\le 1$.
Finally, for the Kullback-Leibler divergence between $\mathbf P_{\mu,\pi_{k,k'}}$ and $\mathbf P_{\mu,\textit{id}}$, where
$\pi_{k,k'}=(k~k')$ is the transposition from $\mathfrak{S}_n$ permuting only $k$ and $k'$, it holds
\begin{align*}
K(\mathbf P_{\mu,\pi_{k,k'}},\mathbf P_{\mu,\textit{id}}) &\le {\kappa^2}\sqrt{\frac{\log n}{d}}+\frac{8\kappa^4}{d}\Big(2+\frac{33^2}{32^2}\big(1+0.25\big)^2\Big)+
\frac{33d}{64}\times\frac{\log n}{16d}\\
&\le \frac{\log n}{8}\le \frac{\log n(n-1)/2}{8},
\end{align*}
where we have used once again the facts that $c\le 1$ and $n\ge 3$.
Applying Lemma~\ref{Tsybakov} with
$M = n(n-1)/2$, $\mathbf Q_0 = \mathbf P_{\mu,\textit{id}}$ and $\{\mathbf Q_j\}_{j=1,\ldots,M}=\{\mathbf P_{\mu,\pi_{k,k'}}\}_{k\not= k'}$,
we obtain
\begin{align*}
\max_{\pi^*\in\mathfrak{S}_n}\sup_{\boldsymbol{\theta}\in\bar\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\pi^*}\big(\hat\pi\neq\pi^*\big) &\ge \max_{\pi^*\in\{\textit{id}\}\cup\{\pi_{k,k'}\}}
\int_{\bar\Theta_\kappa} \mathbf P_{\boldsymbol{\theta},\pi^*} \big( \hat\pi \neq \pi^* \big) \frac{\mu(d\boldsymbol{\theta})}{\mu(\bar\Theta_\kappa)} \\
&\ge \max_{\pi^*\in\{\textit{id}\}\cup\{\pi_{k,k'}\}} \mathbf P_{\mu,\pi^*} \big( \hat\pi\neq\pi^*\big) - \mu(\mathcal E\backslash\bar\Theta_\kappa) \\
&\ge \frac{\sqrt{15}}{\sqrt{15}+1} \Big( \frac{3}{4} - \frac{1}{2\sqrt{\log15}} \Big) - \frac{n(n-1)}2 e^{-d/8}.
\end{align*}
In view of the inequalities $d\ge (1/c)\log n$, $c\le 1/20$ and $n\ge 6$, we get
the inequality $\max_{\pi^*\in\mathfrak{S}_n}\sup_{\boldsymbol{\theta}\in\bar\Theta_\kappa}\mathbf P_{\boldsymbol{\theta},\pi^*}\big(\hat\pi\neq\pi^*\big)\ge 22.4\%$.
\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{lowerbound3}]
Since the event $\{\pi^{\textrm{gr}}\not =\textit{id}\}$  includes the event
$$
\Omega_2=\{\|X_{1}-X^{\text{\tt\#}}_1\|^2>\|X_{2}-X^{\text{\tt\#}}_1\|^2\},
$$
it is sufficient to bound from below the probability of $\Omega_2$.
To this end, we choose any $\boldsymbol{\theta}\in\mathbb R^{n\times d}$ satisfying $\|\theta_1-\theta_2\|=2\kappa$. This readily implies that
$\boldsymbol{\theta}$ belongs to $\bar\Theta_\kappa$.  Furthermore, for suitably chosen random variables $\eta_1\sim \chi^2_d$,
$\eta_2\sim \chi^2_d$ and $\eta_3\sim\mathcal N(0,1)$, it holds that
$\|X_{1}-X^{\text{\tt\#}}_1\|^2-\|X_{2}-X^{\text{\tt\#}}_1\|^2 = 6\eta_1 - 4\kappa^2-8\kappa \eta_3 - 4\eta_2$.
The random terms in the last sum can be controlled using Lemma~\ref{concentration}.
More precisely, for every $x>0$, each one of the following three inequalities holds true
with probability at least $1-e^{-x^2}$:
\begin{align*}
\eta_1 \ge  d- 2\sqrt{d} x,\quad
\eta_2 \le  d+ 2\sqrt{d} x + 2x^2,\quad
\eta_3 \le \sqrt{2} x.
\end{align*}
This implies that with probability at least $1-3e^{-x^2}$, we have
\begin{equation*}
\|X_{1}-X^{\text{\tt\#}}_1\|^2-\|X_{2}-X^{\text{\tt\#}}_1\|^2 \ge  2d-20\sqrt{d} x-4(\kappa + \sqrt2 x)^2.
\end{equation*}
If $x=\sqrt{\log 6}$, then the conditions imposed on $\kappa$  and $d$ ensure that the right-hand side of the last inequality is positive. Therefore,
$\mathbf P(\Omega_2)\ge 1-3e^{-x^2}=1/2$.
\end{proofof}

\begin{proofof}[Proof of Theorem~\ref{lowerbound4}]
The proof is split into two parts. In the first part, we consider the case
$\kappa\le \frac14\sqrt{\frac{\log M_n}{n}}$, while in the second part the case $\kappa\le \frac18\big(\frac{\log n}{d}\big)^{1/4}$ with
$d\ge 24\log n$ and is analyzed. In both cases, the main tool we use is the following result.

\begin{lemma}[\citet{Tsybakov2009}, Theorem 2.5]\label{Tsybakov2}
Assume that for some integer $M\ge 2$ there exist distinct permutations $\pi_0,\ldots,\pi_M\in\mathfrak{S}_n$ and mutually absolutely continuous
probability measures $\mathbf Q_0,\ldots,\mathbf Q_M$ defined on a common probability space $(\mathcal Z,\mathscr Z)$ such that
\begin{equation*}
\begin{cases}
\ \exists\ s>0,\ \forall\ i\neq j, \ \delta(\pi_i,\pi_j)\ge 2s, \\
\  \frac{1}{M} \sum_{j=1}^M K(\mathbf Q_j,\mathbf Q_0) \le \frac{1}{8} \log M.
\end{cases}
\end{equation*}
Then, for every measurable mapping $\tilde\pi:\mathcal Z\to\mathfrak{S}_n$,
\begin{equation*}
\max_{j=0,\ldots,M} \mathbf Q_j\big(\delta(\hat\pi,\pi_j)\ge s\big) \ge \frac{\sqrt{M}}{\sqrt{M}+1} \Big( \frac{3}{4}-\frac{1}{2\sqrt{\log M}} \Big).
\end{equation*}
\end{lemma}

We now have to choose $M$ and $\pi_0,\ldots,\pi_M$ in a suitable manner, which will be done differently according to the relationship between
$n$ and $d$.

\paragraph{\textbf{Case 1:} We assume that $\kappa\le \frac14\sqrt{\frac{\log M_n}{n}}$}
Let $M = M_n$ with $M_n= \mathcal M(1/4, B_{2,n}(2),\delta_H)$ and let $\boldsymbol{\theta} = (\theta_1,\ldots,\theta_n)$ be the set of vectors
$\theta_k = k\kappa\sigma(1,0,\ldots,0)\in\mathbb R^d$. Clearly, $\boldsymbol{\theta}$ belongs to $\bar\Theta_\kappa$. By definition of the packing number, there exist $\pi_1,\ldots,\pi_{M_n}$, permutations from $\mathfrak{S}_n$,
such that
$$
\delta_2(\pi_j,\textit{id})\le 2,\qquad \delta_H(\pi_i,\pi_j)\ge \frac14;\qquad\forall i,j\in\{1,\ldots,M_n\}, \ i\not=j.
$$
Defining $\mathbf Q_j=\mathbf P_{\boldsymbol{\theta},\pi_j}$ for $j=1,\ldots,M_n$ and $\mathbf Q_0=\mathbf P_{\boldsymbol{\theta},\textit{id}}$ we get
\begin{align*}
K(\mathbf Q_j,\mathbf Q_0)
	&= \frac1{2\sigma^2} \sum_{k=1}^n \|\theta_{\pi_j(k)}-\theta_{k}\|^2 = \frac{\kappa^2}{2}\sum_{k=1}^n (\pi_j(k)-k)^2\\
	&= \frac{n\kappa^2}{2}\;\delta_2(\pi_j,\textit{id})^2\le 2n\kappa^2.
\end{align*}
Therefore, using Lemma~\ref{Tsybakov2} with $s=1/8$ we infer from $\kappa\le \frac14\big(\frac{\log M_n}{n}\big)^{1/2}$ that
$$
\min_{\hat\pi}\max_{j=0,\ldots,M_n} \mathbf P_{\boldsymbol{\theta},\pi_j}(\delta_H(\hat\pi,\pi_j)\ge 1/8) \ge \frac{\sqrt{3}}{\sqrt{3}+1}\Big(\frac34-\frac1{2\sqrt{\log 3}}\Big) \approx 17.31\%.
$$
As a consequence, we obtain that
\begin{align*}
\min_{\hat\pi}\max_{(\pi,\boldsymbol{\theta})\in\mathfrak{S}_n\times\bar\Theta_\kappa} \mathbf E_{\boldsymbol{\theta},\pi}[\delta_H(\hat\pi,\pi)]
    &\ge \min_{\hat\pi}\max_{j=0,\ldots,M} \mathbf E_{\boldsymbol{\theta},\pi}[\delta_H(\hat\pi,\pi)\mathds{1}_{\{\delta_H(\hat\pi,\pi)\ge 1/8\}}]\\
    &\ge \frac18 \min_{\hat\pi}\max_{j=0,\ldots,M} \mathbf E_{\boldsymbol{\theta},\pi}[\mathds{1}_{\{\delta_H(\hat\pi,\pi)\ge 1/8\}}]\\
    &\ge 2.15\%.
\end{align*}
This completes the proof of the first case.

\paragraph{\textbf{Case 2:} We assume that $d\ge 24\log n$ and $\kappa \le \frac18 (d\log n)^{1/4}$ }
Let $\mu$ be the uniform distribution on $ \{\pm\epsilon\}^{m\times d}$ with $\epsilon={\sqrt{2/d}\;\sigma\kappa}$, as in
Lemma~\ref{kullback_leibler_divergence_heteroscedastic}. For any set of permutations $\{\pi_0,\ldots,\pi_M\}\subset\mathfrak{S}_n$,
in view of Markov's inequality,
\begin{align*}
\sup_{(\pi,\boldsymbol{\theta})\in\mathfrak{S}_n\times\bar\Theta_\kappa} \mathbf E_{\boldsymbol{\theta},\pi} \big[\delta_H(\hat\pi,\pi)\big] &\ge
\frac{3}{16} \Big( \max_{i=0,\ldots,M} \mathbf P_{\mu,\pi_i} \Big(\delta_H(\hat\pi,\pi_i) \ge \frac{3}{16} \Big) - \mu(\bar\Theta_\kappa^\complement)\Big).
\end{align*}
We choose $M$ and $\pi_0,\ldots,\pi_M$ as in the following lemma.
\begin{lemma}\label{existence2}
For any integer $n\ge 4$ there exist permutations $\pi_0,\ldots,\pi_M\in\mathfrak{S}_n$ such that
$$
\pi_0=\textit{id},\qquad\ M\ge(n/24)^{n/6},
$$
each $\pi_i$  is a composition of at most $n/2$ transpositions with disjoint supports,
and for every distinct pair of indices $i,j\in\{0,\ldots,M\}$ we have
\begin{equation*}
\delta_H(\pi_i, \pi_j) \ge {3}/{8}.
\end{equation*}
\end{lemma}
As $\pi_i$ is a product of transpositions, the Kullback-Leibler divergence between $\mathbf P_{\mu,\pi_i}$ and $\mathbf P_{\mu,\pi_0}$ can be computed by independence thanks to Lemma~\ref{kullback_leibler_divergence_heteroscedastic}:
\begin{equation*}
\frac1M\sum_{i=1}^M K(\mathbf P_{\mu,\pi_i},\mathbf P_{\mu,\pi_0}) \le \frac{n}{2} \times  \frac{8\kappa^4}{d}\Big(2+\Big[1+\frac{2\kappa^2}{d}\Big]^2\Big) = \frac{16n\kappa^4}{d},
\end{equation*}
where the last inequality follows from the bound $\kappa\le 0.45 d^{1/2}$. For $n\ge 26$, it holds that $M\ge 2$ and
$$
\log M\ge \frac{n(\log n -\log 24)}{6}\ge \frac{\log n}{512}.
$$
Consequently, $\frac{16n\kappa^4}{d}\le \frac18\log M$ which allows us to apply Lemma~\ref{Tsybakov2}. This yields
\begin{align*}
\inf_{\hat\pi} \sup_{(\pi,\boldsymbol{\theta})\in\mathfrak{S}_n\times\bar\Theta_\kappa} \mathbf E_{\boldsymbol{\theta},\pi} [\delta_H(\hat\pi,\pi)]
&\ge \frac{3}{16}\Big[\frac{\sqrt{2}}{\sqrt{2}+1} \Big(\frac{3}{4}-\frac{1}{2\sqrt{\log 2}} \Big) - \frac{n^2}{2}e^{-d/8}\Big]\\
&\ge \frac{3}{16}\Big[0.077-\frac{n^2}{2}e^{-24\log(n)/8}\Big]\\
&\ge \frac{3}{16}\Big[0.077-\frac{1}{2n}\Big]\ge 5.81\%.
\end{align*}
\vspace{-28pt}

~
\end{proofof}

\section{Proofs of the Lemmas}\label{section:lemmas}

\begin{proofof}[Proof of Lemma~\ref{construction_btheta}]
Let us denote $r_\sigma =\max_{1\le\ell\le n} \sigma_\ell/\min_{1\le\ell\le n} \sigma_\ell$.  It suffices to set $\theta_1=0\in\mathbb R^d$,
\begin{align*}
\theta_{2k+1}&=\kappa\Big(\sigma_{1,2}+\ldots+\sigma_{2k-1,2k}+k\big(1+r_\sigma\big),0,\ldots,0\Big)\in\mathbb R^d,\\
\theta_{2k}&=\kappa\Big(\sigma_{1,2}+\ldots+\sigma_{2k-1,2k}+(k-1)\big(1+r_\sigma\big),0,\ldots,0\Big)\in\mathbb R^d
\end{align*}
for all $k=1,\ldots,m-1$. If $n$ is impair, one can set $\theta_n=\theta_{n-1}+\kappa \big(1+r_\sigma\big) (1,0,\ldots,0)$.
One readily checks that these vectors satisfy the desired conditions.
\end{proofof}

\begin{proofof}[Proof of Lemma~\ref{kullback_leibler_divergence_heteroscedastic}]
Without loss of generality, we assume hereafter that $\pi\in\mathfrak{S}_n$ is the transposition permuting $1$ and $2$.
Recall that the uniform distribution on $\{\pm\epsilon_1\}^d\times\ldots\times\{\pm\epsilon_{n}\}^d$ can also be written as the product
$\mu = \bigotimes_{\ell=1}^n\mu_\ell$, where $\mu_\ell$ is the uniform distribution on $\{\pm\epsilon_{\ell}\}^d$.
Let us introduce an auxiliary probability distribution $\tilde\mu$ on $\mathbb R^{n\times d}$ defined as $\tilde\mu = \delta_{\mathbf 0}\otimes\delta_{\mathbf 0}
\otimes \mu_2\otimes\ldots\otimes\mu_m$ with $\delta_{\mathbf 0}$ being the Dirac delta measure at $\mathbf 0\in\mathbb R^d$.
 We set $\mathbf P_{\tilde\mu,\textit{id}}(\cdot)= \int_{\Theta} \mathbf P_{\boldsymbol{\theta},\textit{id}}(\cdot)\tilde\mu(d\boldsymbol{\theta})$.

We first compute the density of $\mathbf P_{\mu,\pi}$ w.r.t.\ $\mathbf P_{\mu,\textit{id}}$, which can be written as
\begin{align*}
\frac{d\mathbf P_{\mu,\pi}}{d\mathbf P_{\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}) = \frac{d\mathbf P_{\mu,\pi}}{d\mathbf P_{\tilde\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}) \bigg/ \bigg(\frac{d\mathbf P_{\mu,\textit{id}}}{d\mathbf P_{\tilde\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}})\bigg),\quad  \boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}\in \mathbb R^{n\times d}.
\end{align*}
For every $\theta_i\in\mathbb R^d$ we denote by $\mathbf P_{\theta_i,\sigma_i}$ the probability distribution of $X_i$ from (\ref{model}), given by
\begin{equation*}
\frac{d\mathbf P_{\theta_i,\sigma_i}}{d\mathbf P_{0,\sigma_j}}(x) = \exp\Big\{-\frac{\|\theta_i\|^2}{2\sigma_i^2} + \frac{1}{\sigma_i^2}(x,\theta_i)-\frac{\|x\|^2}{2}(\sigma_i^{-2}-\sigma_j^{-2})\Big\},\quad \forall x\in\mathbb R^d.
\end{equation*}
With this notation, we have
\begin{align*}
&\frac{d\mathbf P_{\mu,\pi}}{d\mathbf P_{\tilde\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}) = \mathbf E_{\mu} \Bigg[\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X_1)\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X_2)\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_1}}(X^{\text{\tt\#}}_1)\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_2}}(X^{\text{\tt\#}}_2)\Bigg] \\
&= \mathbf E_{\mu} \Bigg[\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X_1)\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_2}}(X^{\text{\tt\#}}_2)\Bigg] \times \mathbf E_{\mu}\Bigg[\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X_2)\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_1}}(X^{\text{\tt\#}}_1)\Bigg]\\
&= \prod_{k=1}^d\cosh\Big(\frac{\epsilon_1}{\sigma_1^2}(X_{1,k}+X^{\text{\tt\#}}_{2,k})\Big)\cosh\Big(\frac{\epsilon_2}{\sigma_2^2}(X_{2,k}+X^{\text{\tt\#}}_{1,k})\Big)\\
&\qquad\times\exp\Big\{-\frac12(\|X^{\text{\tt\#}}_{1}\|^2-\|X^{\text{\tt\#}}_{2}\|^2)(\sigma_2^{-2}-\sigma_1^{-2})\Big\}.
\end{align*}
Similarly,
\begin{align*}
\frac{d\mathbf P_{\mu,\textit{id}}}{d\mathbf P_{\tilde\mu,\textit{id}}}&(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}) =\mathbf E_{\mu} \Bigg[\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X_1)\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X_2)\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X^{\text{\tt\#}}_1)
\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X^{\text{\tt\#}}_2)\Bigg] \\
&= \mathbf E_{\mu} \Bigg[\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X_1)\frac{d\mathbf P_{\theta_1,\sigma_1}}{d\mathbf P_{0,\sigma_1}}(X^{\text{\tt\#}}_1)\Bigg] \times \mathbf E_{\mu}\Bigg[\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X_2)\frac{d\mathbf P_{\theta_2,\sigma_2}}{d\mathbf P_{0,\sigma_2}}(X^{\text{\tt\#}}_2)\Bigg]\\
&= \prod_{k=1}^d\cosh\Big(\frac{\epsilon_1}{\sigma_1^2}(X_{1,k}+X^{\text{\tt\#}}_{1,k})\Big)\cosh\Big(\frac{\epsilon_2}{\sigma_2^2}(X_{2,k}+X^{\text{\tt\#}}_{2,k})\Big).
\end{align*}
Thus, we get that
\begin{align*}
\frac{d\mathbf P_{\mu,\pi}}{d\mathbf P_{\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}}) &= \prod_{k=1}^d \frac{\cosh\big(\frac{\epsilon_1}{\sigma_1^2}(X_{1,k}+X^{\text{\tt\#}}_{2,k})\big)}{\cosh\big(\frac{\epsilon_1}{\sigma^2_1}(X_{1,k}+X^{\text{\tt\#}}_{1,k})\big)}
\times \prod_{k=1}^d \frac{\cosh\big(\frac{\epsilon_2}{\sigma_2^2}(X_{2,k}+X^{\text{\tt\#}}_{1,k})\big)}{\cosh\big(\frac{\epsilon_2}{\sigma^2_2}(X_{2,k}+X^{\text{\tt\#}}_{2,k})\big)}\\
&\qquad\times\exp\Big\{-\frac12(\|X^{\text{\tt\#}}_{1}\|^2-\|X^{\text{\tt\#}}_{2}\|^2)(\sigma_2^{-2}-\sigma_1^{-2})\Big\}.
\end{align*}
Then, we compute the Kullback-Leibler divergence,
\begin{align}
K&(\mathbf P_{\mu,\pi},\mathbf P_{\mu,\textit{id}})
	 = \int \log \Big(\frac{d\mathbf P_{\mu,\pi}}{d\mathbf P_{\mu,\textit{id}}}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}})\Big) d\mathbf P_{\mu,\pi}(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}})\nonumber\\
	& =  \sum_{k=1}^d\sum_{j=1}^2 \bigg\{\mathbf E_{\mu} \bigg[ \int \log\cosh\Big[ \frac{\epsilon_j}{\sigma^2_j}(2\theta_{j,k}+\sigma_j\sqrt{2}x)\Big] \varphi(x)dx \bigg] \nonumber\\
	& \qquad -  \mathbf E_{\mu}\bigg[ \int \log\cosh \Big[\frac{\epsilon_j}{\sigma^2_j}(\theta_{1,k}+\theta_{2,k}+\sigma_{12}x)\Big]\varphi(x) dx \bigg]\bigg\}\nonumber\\
	&\qquad +\frac{d}2\mathbf E_{\mu}\bigg[ \int_\mathbb R ((\theta_{1,1}+\sigma_1 x)^2-(\theta_{2,1}+\sigma_2 x)^2)\varphi(x)\,dx\bigg](\sigma_2^{-2}-\sigma_1^{-2}),\nonumber
\end{align}
where $\varphi$ is the density function of the standard Gaussian distribution. We evaluate the first two terms of the last display using the following inequalities:
\begin{equation}\label{cosh}
\forall u\in\mathbb R,\qquad
\frac{u^2}{2}-\frac{u^4}{12} \le \log\cosh(u) \le \frac{u^2}{2},
\end{equation}
while for the third term the exact computation yields:
\begin{align*}
\mathbf E_{\mu}\bigg[ \int_\mathbb R ((\theta_{1,1}+\sigma_1 x)^2-(\theta_{2,1}+\sigma_2 x)^2)\varphi(x)\,dx\bigg]&=
\epsilon_1^2+\sigma_1^2-\epsilon_2^2-\sigma_2^2\\
&= (\sigma_1^2-\sigma_2^2)(1+(2/d)\kappa^2).
\end{align*}
In conjunction with the facts that $\epsilon_1/\sigma_1=\epsilon_2/\sigma_2$, $\sigma_1\le\sigma_2$ and $\epsilon_1\le\epsilon_2$, this leads to
\begin{align*}
&\mathbf E_{\mu} \bigg[ \int \log\cosh\Big[ \frac{\epsilon_j}{\sigma_j^2}(2\theta_{j,k}+\sigma_j\sqrt{2}x)\Big] \varphi(x)dx \le
\frac{\epsilon_j^2}{\sigma_j^2} + 2\frac{\epsilon_j^4}{\sigma_j^4}=\frac{\epsilon_2^2}{\sigma_2^2} + 2\frac{\epsilon_2^4}{\sigma_2^4},\\
&\mathbf E_{\mu}\bigg[ \int \log\cosh \Big[\frac{\epsilon_j}{\sigma_j^2}(\theta_{1,k}+\theta_{2,k}+\sigma_{1,2}x)\Big]\varphi(x) dx \bigg]\\
&\qquad\ge\frac{\epsilon^2_j}{2\sigma^4_j}\Big(\epsilon_1^2+\epsilon_2^2+\sigma_1^2+\sigma_2^2\Big)\\
&\qquad\qquad- \frac{\epsilon^4_j}{12\sigma^8_j}\Big(\epsilon_1^4+\epsilon_2^4+3(\sigma_1^2+\sigma_2^2)^2+6\epsilon_1^2\epsilon_2^2+6(\sigma_1^2+\sigma_2^2)(\epsilon_1^2+\epsilon_2^2)\Big)\\
&\qquad\ge\frac{\epsilon^2_2(\epsilon_1^2+\sigma_1^2)}{\sigma^4_2}- \frac{\epsilon^4_1(\epsilon_2^2+\sigma_2^2)^2}{\sigma^8_1}.
\end{align*}
Thus, we get that
\begin{align*}
(1/d)K(\mathbf P_{\mu,\pi},\mathbf P_{\mu,\textit{id}}) &\le \frac{2\epsilon_2^2}{\sigma_2^2} + \frac{4\epsilon_2^4}{\sigma_2^4}-\frac{2\epsilon^2_2(\epsilon_1^2+\sigma_1^2)}{\sigma^4_2}+\frac{2\epsilon^4_1(\epsilon_2^2+\sigma_2^2)^2}{\sigma^8_1}\\
&\qquad+ \frac12\big(1+(2/d)\kappa^2)(\sigma_1^2-\sigma_2^2)(\sigma_2^{-2}-\sigma_1^{-2})\\
&\le \frac{4\kappa^2}{d}\Big(1-\frac{\sigma_1^2}{\sigma_2^2}\Big)+\frac{16\kappa^4}{d^2}+\frac{8\kappa^4}{d^2}\big(1+(2/d)\kappa^2\big)^2\frac{\sigma_2^4}{\sigma_1^4}\\
&\qquad+ \frac12\big(1+(2/d)\kappa^2\big)\Big(\frac{\sigma_2^2}{\sigma_1^2}-1\Big)^2.
\end{align*}
To complete the proof, we need to evaluate $\mu(\mathcal E\setminus\bar\Theta_\kappa)$. We note that in view of the union bound,
\begin{align*}
\mu(\mathcal E\setminus\bar\Theta_\kappa) & = \mu\Big(\bigcup_{k=1}^n\bigcup_{k'\not=k}\{\boldsymbol{\theta} : \|\theta_{k}-\theta_{k'}\| < \kappa\sigma_{k,k'} \}\Big)\\
&\le \frac{n(n-1)}{2} \max_{k\not=k'}\mu\big(\{\boldsymbol{\theta} : \|\theta_{k}-\theta_{k'}\|^2 < \kappa^2\sigma_{k,k'}\}\big)\\
&= \frac{n(n-1)}{2}  \max_{k\not=k'}\mathbf P\big(d\epsilon_{k}^2+d\epsilon_{k'}^2-2d\epsilon_{k}\epsilon_{k'}\bar\zeta < \kappa^2\sigma_{k,k'}^2 \big),
\end{align*}
where $\bar\zeta=\frac1d\sum_{j=1}^d \zeta_j$ with $\zeta_1,\ldots,\zeta_d$ being i.i.d.\ Rademacher random variables (\textit{i.e.}, random variables taking the values $+1$ and $-1$
with probability $1/2$). One easily checks that
$$
\frac{d\epsilon_{k}^2+d\epsilon_{k'}^2-\kappa^2\sigma_{k,k'}^2}{2d\epsilon_{k}\epsilon_{k'}} = \frac{2\sigma_{k}^2+2\sigma_{k'}^2-(\sigma_{k}^2+\sigma_{k'}^2)}{4\sigma_{k}\sigma_{k'}}\ge \frac12.
$$
Therefore, using the Hoeffding inequality, we get $\mu(\mathcal E\setminus\bar\Theta_\kappa)  \le \frac12n(n-1) \mathbf P\big(\bar\zeta > 1/2 \big)\le \frac12n(n-1) e^{-d/8}$.
\end{proofof}

\begin{proofof}[Proof of Lemma~\ref{existence2}]
We first prove an auxiliary result.

\begin{lemma}\label{existence}
For any integer $n\ge 2$ there exist permutations $\pi_0,\pi_1,\ldots,\pi_M$ in $\mathfrak{S}_n$ such that
$$
\pi_0=\textit{id},\qquad\ M \ge (n/8)^{n/2}
$$
and for any pair $i,j\in\{0,\ldots,M\}$ of distinct indices we have
$\delta_H(\pi_i, \pi_j) \ge \frac{1}{2}$.
\end{lemma}

\begin{proof}
When $n\le 8$, the claim of this lemma is trivial since one can always find at least one permutation that
differs from the identity at all the positions and thus $M\ge 1\ge (n/8)^{n/2}$. Let us consider the case
$n>8$. For every $\pi\in\mathfrak{S}_n$, denote
\begin{equation*}
E_\pi \triangleq \Big\{ \pi'\in\mathfrak{S}_n \,\big|\, \frac1n \sum_{i=1}^n \mathds{1}_{\{\pi(i)\neq \pi'(i)\}} \ge 1/2 \Big\}.
\end{equation*}
We first notice that for every $\pi\in\mathfrak{S}_n$, there is a one-to-one correspondence between $E_{id}$ and $E_\pi$ through the bijection
\begin{equation*}
\fonction{\phi}{E_{id}}{E_\pi\phantom{\Big()}}{\pi'}{\pi\circ\pi'},
\end{equation*}
so that $\# E_\pi = \# E_{id}$. The following lemma, proved later on in this section, gives a bound for this number.
\begin{lemma}\label{cardinalEid}
Let $n\ge 2$ be an integer and $m$ be the smallest integer such that $2m\ge n$.
Then
\begin{equation*}
\# E_{id}^\complement \le \frac{4n!}{m!}.
\end{equation*}
\end{lemma}
Now we denote $\pi_0=id$ and choose $\pi_1$ in $E_{id}$. Then, it is sufficient to choose $\pi_2$ as any element from $E_{id}\cap E_{\pi_1}$,
the latter set being nonempty since
\begin{align*}
\#\big(E_{\pi_0}\cap E_{\pi_1}\big) &\ge \# \mathfrak{S}_n - \# E_{\pi_0}^\complement - \# E_{\pi_1}^\complement \\
&\ge n!\times\Big(1-\frac{8}{m!}\Big) > 0.
\end{align*}
We can continue the construction until $\pi_i$ if
\begin{equation*}
1-\frac{4i}{m!} > 0\qquad \ssi\qquad i < \frac{m!}{4}.
\end{equation*}
To conclude, we observe that
\begin{equation*}
\frac{m!}{4} > \frac14 \Big( \frac{n}{2e} \Big)^{n/2} \ge \Big( \frac{n}{8} \Big)^{n/2}.
\end{equation*}
\vspace{-20pt}

~
\end{proof}

Let us denote by $m$ the largest integer such that $2m\le n$, and choose
$$
\tilde\pi_0,\tilde\pi_1,\ldots,\tilde\pi_M\in\mathfrak{S}_m \quad \text{ with } \quad M\ge(m/8)^{m/2}
$$
as in Lemma~\ref{existence}, so that for every $i\neq j\in\{0,\ldots,M\}$,
$\delta_H(\tilde\pi_i,\tilde\pi_j) \ge \frac{1}{2}$.
We use each permutation $\tilde\pi_i\in\mathfrak{S}_m$ to construct a permutation $\pi_i\in\mathfrak{S}_n$.
The idea of the construction is as follows: the permutation $\pi_i$ is a product of $m$ transpositions of
distinct supports, and each transposition permutes an even integer with an odd one.
We set $\pi_0=id$ and for every $i$ in $\{1,\ldots,M\}$,
\begin{equation*}
\pi_i = \big(1 \ 2\tilde\pi_i(1)\big) \circ \big(3 \ 2\tilde\pi_i(2)\big) \circ \ldots \circ \big(2m-1 \ 2\tilde\pi_i(m)\big) \in \mathfrak{S}_n.
\end{equation*}
With these choices, the number of differences between $\pi_i$ and $\pi_j$ is exactly twice as much as the number of differences
between $\tilde\pi_i$ and $\tilde\pi_j$. To sum up, for every pair of distinct indices $i,j\in\{0,\ldots,M\}$,
\begin{equation*}
\frac1n \sum_{k=1}^n \mathds{1}_{\{\pi_i(k)\neq \pi_j(k)\}} = \frac{2m}{n} \times \frac1m \sum_{k=1}^m \mathds{1}_{\{\tilde\pi_i(k)\neq \tilde\pi_j(k)\}} \ge \frac{m}{n}\ge\frac{3}{8},\qquad \forall n\ge 4.
\end{equation*}
To complete the proof, we note that $m\ge n/3$.
\end{proofof}

\begin{proofof}[Proof of Lemma~\ref{cardinalEid}]
For every $\ell\in\{m,\ldots,n\}$, counting all the permutations $\pi$ such that $\sum_{k=1}^n \mathds{1}_{\{\pi(k)\neq k\}}=\ell$, we get
\begin{equation*}
\# E_{id} = !n + !(n-1)\binom{n}{1}  + \ldots + !(n-m)\binom{n}{m} ,
\end{equation*}
where $!\ell$ is the number of derangements, the permutations such that none of the elements appear in their original position, in $\mathfrak{S}_\ell$
for $\ell\ge 1$. We know that
\begin{equation*}
\forall \ \ell\ge1, \qquad !\ell = \ell! \times \sum_{j=0}^\ell \frac{(-1)^j}{j!},
\end{equation*}
which, using the alternating series test, yields
\begin{equation*}
\forall \ \ell\ge1, \qquad !\ell \ge \ell !\times\Big(e^{-1}-\frac{1}{(\ell+1)!}\Big).
\end{equation*}
It follows that
\begin{align*}
\# E_{id} &\ge n!\times\Big(e^{-1}-\frac{1}{(n-m+1)!}\Big)\times\Big( 1 + \frac{1}{1!} + \ldots + \frac{1}{m!} \Big) \\
&\ge n!\times\Big(e^{-1}-\frac{1}{(n-m+1)!}\Big)\times\Big(e-\frac{e}{(m+1)!}\Big)\\
&\ge n!\times\Big(1-\frac{e}{(n-m+1)!} - \frac{1}{(m+1)!} \Big).
\end{align*}
Therefore,
\begin{equation*}
\# E_{id}^\complement \le n!\times\Big( \frac{e}{(n-m+1)!} + \frac{1}{(m+1)!} \Big) \le \frac{4n!}{m!}.
\end{equation*}
\vspace{-20pt}

~
\end{proofof}

\section*{Acknowledgments}
This work was partially supported by the grants Investissements d'Avenir (ANR-11-IDEX-0003/Labex Ecodec/ANR-11-LABX-0047) and
CALLISTO. The authors thank the Reviewers for many valuable suggestions. Special thanks to an anonymous Referee for pointing a
mistake in the proof of Theorem~\ref{lowerbound2}.

\section{Introduction}

\section{Problem Formulation}\label{section:model}

We begin with formalizing the problem of matching two sets of features $\{X_1,\ldots,X_n\}$ and $\{X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m\}$ with different sizes $n$ and $m$ such that $m \ge n \ge 2$.
In what follows we assume that the observed features are randomly generated from the model
\begin{equation}\label{model}
\begin{cases}
X_i = \theta_i + \sigma_i\xi_i \ , \\
X^{\text{\tt\#}}_j = {\theta}^\text{\tt\#}_j + \s^{\text{\tt\#}}_j\xi^\text{\tt\#}_j,
\end{cases}\quad i=1,\ldots,n \text{ and } j=1,\ldots,m
\end{equation}
where
\begin{itemize}\setlength{\itemsep}{-2pt}
\item $\boldsymbol{\theta}=\{\theta_1,\ldots,\theta_n\}$ and $\boldsymbol{\theta}^\text{\tt\#}=\{{\theta}^\text{\tt\#}_1,\ldots,{\theta}^\text{\tt\#}_m\}$ are two collections of vectors from $\mathbb R^d$, corresponding to
the original features, which are unavailable,
\item $\boldsymbol{\sigma}=(\sigma_1,\ldots,\sigma_n), \boldsymbol{\sigma}^\text{\tt\#} = (\s^{\text{\tt\#}}_1,\ldots,\s^{\text{\tt\#}}_m)$ are positive real numbers corresponding to the levels of noise contaminating each feature,
\item $\xi_1,\ldots,\xi_n$ and $\xi^\text{\tt\#}_1,\ldots,\xi^\text{\tt\#}_m$ are two independent sets of i.i.d.\ random vectors drawn from
    the Gaussian distribution with zero mean and identity covariance matrix.
\end{itemize}

In the most generic setting both $\boldsymbol{X} = \{X_1, \dots, X_n\}$ and $\boldsymbol{X}^{\text{\tt\#}} = \{X^{\text{\tt\#}}_1, \dots, X^{\text{\tt\#}}_m\}$ may contain outliers, i.e. features that don't have their corresponding pair. In other words, the task is to find two sets of features $S_k$ and $S^{\text{\tt\#}}_{\pi(k)}$ of cardinality $k < n$ such that $S_k \subset \boldsymbol{X}$ and $S^{\text{\tt\#}}_{\pi(k)} \subset \boldsymbol{X}^{\text{\tt\#}}$, where the parameter $\pi(\cdot)$ gives a correspondence between these sets. In this case there are $n + m - 2k$ outliers and the value of $k$ is unknown. Thus, for simplicity we only consider the case where $k=n$, hence only the larger set of features, namely $\boldsymbol{X}^{\text{\tt\#}}$, contains outliers.

To formalize the task of feature matching with outliers we aim to find an injection $\pi^*: \{1,
\dots, n\} \to \{1, \dots, m\}$ such that
\begin{equation}\label{pistar}
\forall\ i\in \{1,
\dots, n\},\quad \theta_i \equiv {\theta}^\text{\tt\#}_{\pi^*(i)},
\end{equation}
where $\equiv$ is an equivalence relation that we call \textit{matching criterion}. We assume that there exists an injection $\pi^*$ such that $\forall i \in \{1, \dots, n\}$ we have $\theta_i = {\theta}^\text{\tt\#}_{\pi^*(i)}$ and $\sigma_i = \s^{\text{\tt\#}}_{\pi^*(i)}$.

In what follows we call a feature $X^{\text{\tt\#}}_j$ an {\it{outlier}}, if $j \notin J_{\pi^*} \triangleq\text{Im}(\pi^*)$, where we have that $|J_{\pi^*}| = n$. Let us also define $O_{\pi^*} \triangleq \{1, \dots, m\} \setminus J_{\pi^*}$, for which we have $|O_{\pi^*}| = m - n$. Consequently, in addition to the features $\boldsymbol{X}$, the rest of the features from $\boldsymbol{X}^{\text{\tt\#}}$ are called {\it{inliers}}.

In this formulation, the data generating distribution is defined by the (unknown) parameters\footnote{We omit the set of parameters $\boldsymbol{\theta}$ and $\boldsymbol{\sigma}$, since they are automatically contained in $\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$, respectively due to the fact that $\pi^*$ is included in the set of parameters.} $\boldsymbol{\theta}^\texttt{\#}$,
$\boldsymbol{\sigma}^\text{\tt\#}$ and $\pi^*$.
In the problem of matching, we focus our attention on the problem of estimating the parameter $\pi^*$ only,
considering $\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ as nuisance parameters. In what follows, we denote by  $\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*}$ the probability distribution
of the vector $(X_1,\ldots,X_n,X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m)$ defined by (\ref{model}) under the condition
(\ref{pistar}). We write $\mathbf E_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*}$ for the expectation
with respect to $\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*}$.

We use the binary loss function, which equals $1$, if $\hat\pi \neq \pi^*$ and $0$ otherwise, i.e. 0\,-1 distance between $\hat\pi$ and $\pi^*$, given by
\begin{equation}\label{dist}
\delta_{\text{0\,-1}}(\hat\pi,\pi^*)\triangleq\mathds{1}_{\{\hat\pi\not=\pi^*\}}.
\end{equation}

Our ultimate goal is to design estimators that have an expected error smaller than a prescribed level $\alpha$ under the weakest possible conditions on the nuisance parameter $\boldsymbol{\theta}^\texttt{\#}$ and noise level $\boldsymbol{\sigma}^\text{\tt\#}$. The problem of matching becomes more difficult with hardly distinguishable features. To quantify this phenomenon, we introduce the relative separation distance $\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})$ and the relative outlier separation distance $\bar\kappa_{\text{out}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$, which measure the minimal distance-to-noise ratio between inlier features and the minimal distance-to-noise ratio between inlier and outlier features, respectively. Denoting $[n] \triangleq \{1, \dots, n\}$ the precise definitions read as
\begin{equation}
\bar\kappa(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*) \triangleq \min_{\substack{\pi^*(i) \neq j, \\ i\in[n], j \in J_{\pi^*}}} \frac{\|\theta_i-{\theta}^\text{\tt\#}_j\|}{(\sigma_i^2+(\s^{\text{\tt\#}}_j)^2)^{1/2}},\qquad
\bar\kappa_{\text{out}}(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*) \triangleq \min_{\substack{ i\in[n], \\ j \in O_{\pi^*}}} \frac{\|\theta_i-{\theta}^\text{\tt\#}_j\|}{(\sigma_i^2+(\s^{\text{\tt\#}}_j)^2)^{1/2}}. \label{kappa}
\end{equation}
Notice that $\bar\kappa(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ can be rewritten as
\begin{align*}
    \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \triangleq \min_{\substack{i, j \in [n] \\ i \neq j}} \frac{\| \theta_i - \theta_j\|}{(\sigma_i^2 + \sigma_j^2)^{1/2}}
\end{align*}
and the dependence of $\bar\kappa(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ on $\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}$ and $\pi^*$ is redundant. Hence, we will refer to $\bar\kappa(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*) \equiv \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})$ as $\bar\kappa$ to ease the notation.

Clearly, if $\bar\kappa(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)=0$ (or, $\bar\kappa_{\text{out}}(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)=0$) and $\sigma_i$'s are all equal, then the parameter $\pi^*$ is nonidentifiable, in the sense that there exist two different permutations $\pi^*_1$ and $\pi^*_2$ such that the distributions $\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_1}$ and $\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_2}$ coincide. Therefore, the conditions $\bar\kappa(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)>0, \bar\kappa_{\text{out}}(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)>0$ are necessary for the existence of consistent estimators of $\pi^*$. Furthermore, good estimators are those consistently estimating $\pi^*$ even if $\bar\kappa(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ and $\bar\kappa_{\text{out}}(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ are small.

In the setting of permutation estimation, i.e. when there are no outliers in both $\boldsymbol{X}$ and $\boldsymbol{X}^{\text{\tt\#}}$, \citep{collier2016minimax} established the optimal minimax rate of relative separation distance $\bar\kappa$. The minimax rate is defined as
\begin{align*}
    \bar\kappa_{\alpha} =  \inf_{\hat\pi} \inf \big\{ \kappa > 0  ~|~ \, \max_{\pi \in \mathfrak{S}_n} \, \sup_{\bar\kappa > \kappa} \mathbf P_{\boldsymbol{\theta}, \boldsymbol{\sigma}, \pi}(\hat\pi \neq \pi) \le \alpha\big\},
\end{align*}
where $\mathfrak{S}_n$ is the symmetric group over the set $[n]$. The minimax rate is proved to be proportional to $(\log n)^{1/2} \vee (d\log n )^{1/4}$. Since this rate is optimal we establish a threshold $\bar\kappa_{\textup{out}}^{\alpha}$ such that for all  $\bar\kappa_{\textup{out}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*) > \bar\kappa_{\textup{out}}^{\alpha}$ and given estimator $\hat\pi$ we have
\begin{align*}
    \bar\kappa_{\textup{out}}^{\alpha}(\hat\pi) = \inf\big\{\kappa > 0 ~|~ \max_{\pi} \sup_{\bar\kappa_{\textup{out}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*) > \kappa} \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi}(\hat\pi \neq \pi) \le \alpha \big\}.
\end{align*}
The minimax outlier separation distance is the smallest possible outlier separation distance achieved by an estimator $\hat\pi$, i.e.
\begin{equation}\label{mmx}
\bar\kappa_{\textup{out}}^\alpha \triangleq \inf_{\hat\pi}~\bar\kappa_{\textup{out}}^\alpha(\hat\pi),
\end{equation}
where the infimum is taken over all possible estimators of $\pi^*$.

In what follows we devote our main attention to the case of  heteroscedastic noise scenario when the $\sigma_i$'s can be different.

\section{Main results}
The general approach here is based on maximum likelihood estimation. Writing the log-likelihood for the vector $(X_1, \dots, X_n, X^{\text{\tt\#}}_1, \dots, X^{\text{\tt\#}}_m)$ yields

\begin{align*}
    -l_n(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi, \boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}) &= \sum_{i=1}^n \bigg(\frac{\| X_i - \theta_i\|_2^2}{2\sigma_i^2} + \frac 1 2 \log(\sigma_i^2)\bigg) \\ &+  \sum_{j=1}^m \bigg(\frac{\| X^{\text{\tt\#}}_j - {\theta}^\text{\tt\#}_j\|_2^2}{2(\s^{\text{\tt\#}}_j)^2} + \frac 1 2 \log((\s^{\text{\tt\#}}_j)^2)\bigg)
\end{align*}

\subsection{Upper bound for $\bar\kappa_{\textup{out}}$}
In this section we obtain an upper bound for $\bar\kappa_{\textup{out}}$

Further, we consider the case of known noise variances, i.e. we assume that all $\s^{\text{\tt\#}}_i$'s are known.

\subsubsection{Known variance}

In the case of known $\sigma_i$'s the estimator of $\pi^*$ reads as
\begin{equation}\label{ml-sigmas}
\bar\pi^{\textrm{LSNS}} \triangleq \arg\min_{\pi:[n]\to [m]}\sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2},
\end{equation}
where the minimum is taken over all injective functions $\pi$ mapping from $[n]$ to $[m]$.

The next theorem states that if both $\bar\kappa_{\alpha}$ and $\bar\kappa^{\textup{out}}_{\alpha}$ are large enough then the probability of finding the correct matching using LSNS method is at least $1 - \alpha$.

\begin{theorem}[Upper bound for LSNS]\label{version1:upperbound}
Let $\alpha \in (0,1)$ and condition \eqref{pistar} be fulfilled. If $\bar\kappa_{\alpha}(\boldsymbol{\theta}, \boldsymbol{\sigma}) \ge \dots$ and $\bar\kappa_{\textup{out}}^{\alpha}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*) \ge \dots $ then
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi \neq \pi^*) \le \alpha.
\end{align}

\end{theorem}

\begin{proofof}[Proof of Theorem~\ref{version1:upperbound}.]

We prove the upper bound for $\bar\kappa$ in the presence of outliers. Without loss of generality we can assume that  $ \pi^*(i)=i, ~ ~ \forall i \in \{1, \dots, n\}$. We wish to bound the probability of the event $\Omega = \{ \hat\pi \neq \pi^* \}$, where $\hat\pi = \bar\pi^{\text{LSNS}}$. It is evident that

\begin{align*}
\Omega \subset \bigcup_{\pi\neq\pi^*} \Omega_{\pi},
\end{align*}
where
\begin{align*}
\Omega_{\pi} &= \bigg\{ \sum_{i=1}^n \frac{\| X_i - X^{\text{\tt\#}}_i \|^2}{2\sigma_i^2} \ge \sum_{i=1}^n \frac{\| X_i - X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_i^2 + (\s^{\text{\tt\#}}_{\pi(i)})^2 }\bigg\} \\ &\subset
\bigcup_{i=1}^n \bigcup_{j \in [m] \setminus \{ i \}} \bigg\{ \frac{\| X_i - X^{\text{\tt\#}}_i \|^2}{2\sigma_i^2} \ge \frac{\| X_i - X^{\text{\tt\#}}_{j}\|^2}{\sigma_i^2 + (\s^{\text{\tt\#}}_{j})^2 }\bigg\}.
\end{align*}
Denote $\sigma_{i, j} = (\sigma_i^2 + (\s^{\text{\tt\#}}_j)^2)^{1/2}$ and
\begin{align}\label{zetas}
\zeta_1 = \max_{i \neq j} \bigg| \frac{(\theta_i - {\theta}^\text{\tt\#}_j)^\top (\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j)}{\|\theta_i - {\theta}^\text{\tt\#}_j\|\sigma_{i, j}}\bigg|, \quad \zeta_2 = d^{-1/2} \max_{i, j} \bigg| \bigg\| \frac{\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j}{\sigma_{i,j}} \bigg\|^2 - d \bigg|.
\end{align}
Since $\pi^*(i) = i$ for every $i \in \{1, \dots, n\}$, then
\begin{align}\label{eq:i=j}
\| X_i - X^{\text{\tt\#}}_i \|^2 = \sigma_i^2 \| \xi_i - \xi^\text{\tt\#}_i\|^2 \le 2\sigma_i^2 (d + \sqrt{d} \zeta_2).
\end{align}
Similarly, for every $j \in \{1, \dots, m\}$ and $j \neq i$,
\begin{align*}
\| X_i - X^{\text{\tt\#}}_j \|^2 &= \| \theta_i - {\theta}^\text{\tt\#}_j \|^2 + \| \sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j \|^2 + 2(\theta_i - {\theta}^\text{\tt\#}_j)^\top (\sigma_i \xi_i - \s^{\text{\tt\#}}_j \xi^\text{\tt\#}_j) \\ &\ge
\| \theta_i - {\theta}^\text{\tt\#}_j \|^2 + \sigma_{i,j}^2 (d - \sqrt{d}\zeta_2) - 2 \| \theta_i - {\theta}^\text{\tt\#}_j \| \sigma_{i,j} \zeta_1.
\end{align*}
On the event $\Omega_1 = \{ \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \ge \zeta_1\}$ from the previous inequality we have

\begin{align}\label{eq:ineqj}
\frac{\| X_i - X^{\text{\tt\#}}_j \|^2}{\sigma_{i,j}^2} \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 - 2\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \zeta_1 + d - \sqrt{d} \zeta_2,
\end{align}
where $$\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) = \min_{\substack{i\in[n], j \in [m] \\ i \neq j}} \frac{\| \theta_i - {\theta}^\text{\tt\#}_j\|}{\sigma_{i,j}}.$$
Hence, combining obtained bounds \eqref{eq:i=j} and \eqref{eq:ineqj} we get that
\begin{align}\label{eq:omega_omega1}
\Omega \cap \Omega_1 \subset \bigg\{ d + \sqrt{d} \zeta_2 \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 - 2\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})\zeta_1 + d - \sqrt{d} \zeta_2\bigg\}.
\end{align}
Using \eqref{eq:omega_omega1} we can show, that
\begin{align}
\mathbf P(\Omega)
        &\le \mathbf P(\Omega_1^\complement) + \mathbf P\big(\Omega\cap\Omega_1\big)\nonumber\\
        &\le \mathbf P\big(\zeta_1\ge \bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
            +\mathbf P(2\sqrt{d}\zeta_2+2\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\zeta_1\ge\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2)\nonumber\\
        &\le 2\mathbf P\big(\zeta_1\ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})}4\big)+
            \mathbf P\big(\zeta_2\ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}{4\sqrt{d}}\big).\label{ineq:2}
\end{align}
For suitably chosen standard Gaussian random variables $\zeta_{i,j}$
it holds that $\zeta_1 = \max_{i\not=j}|\zeta_{i,j}|$. Therefore, using the tail bound for the standard Gaussian
distribution and the union bound, we get
\begin{align}
\mathbf P\big(\zeta_1\ge \textstyle\frac14\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
    &\le \sum_{i\not =j}\nolimits \mathbf P\big(|\zeta_{i,j}|\ge \textstyle\frac14\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\big)
    \le 2 n m e^{-\frac1{32}\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}.\label{ineq:3}
\end{align}
To bound the random variable $\zeta_2$, we use the following result.
\begin{lemma}[\cite{LaurentMassart2000}, Eq. (4.3) and (4.4)]\label{concentration}
If $Y$ is drawn from the chi-squared distribution  $\chi^2(D)$, where $D\in\mathbb N^*$, then, for every $x>0$,
\begin{equation*}
\begin{cases}
\ \mathbf P\big(Y-D \le -2\sqrt{Dx} \big) \le e^{-x}, \phantom{\Big()}\\
\ \mathbf P\big(Y-D \ge 2\sqrt{Dx} + 2x \big) \le e^{-x}. \phantom{\Big()}
\end{cases}
\end{equation*}
As a consequence, $\forall y>0$,
$\mathbf P\big(D^{-1/2}|Y-D| \ge y\big) \le 2\exp\big\{-\textstyle\frac18 y(y\wedge \sqrt{D})\big\}$.
\end{lemma}
This inequality, combined with the union bound, yields
\begin{align}
\mathbf P\Big(\zeta_2 \ge \textstyle\frac{\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})^2}{4\sqrt{d}}\Big)
        &\le 2 n m \exp\Big\{-\frac{(\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})/16)^2}{d}(\bar\kappa^2(\boldsymbol{\theta},\boldsymbol{\sigma})\wedge 8d)\Big\}.\label{ineq:4}
\end{align}
Using inequalities (\ref{ineq:2})-(\ref{ineq:4}), we get that
\begin{align*}
\bar\kappa(\boldsymbol{\theta},\boldsymbol{\sigma})\ge 4\Big(\sqrt{2\log(8n m/\alpha)} \vee \big(d\log(4n m/\alpha)\big)^{1/4}\Big),
\end{align*}
which implies $\mathbf P(\hat\pi\neq \pi^*)=\mathbf P(\Omega)\le \alpha.$

\end{proofof}

\subsubsection{Unknown variance}
In the setup when the observation variances are unknown the minimization problem displayed in \eqref{ml-sigmas} can be further minimized with respect to parameters $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^\text{\tt\#}$, taking into account the constraint \eqref{pistar}. This readily yields
\begin{align}
    \bar\pi^{\textup{LSL}} \triangleq \arg\min_{\pi :[n] \to [m]} \sum_{i=1}^n \log \| X_i - X^{\text{\tt\#}}_{\pi(i)} \|^2.
\end{align}

\begin{theorem}[Upper bound for LSL]\label{thm:upper-bound-LSL}
Let $\alpha \in (0,1)$ then for $\hat\pi = \bar\pi^{\textup{LSL}}$ we have
\begin{align}
    \bar\kappa^{\alpha}_{\textup{out}}(\hat\pi)^2 \le d/2 + \Big(\sqrt{2d \log \frac{8nm}{\alpha}} \vee 32\log\frac{4nm}{\alpha}\Big).
\end{align}
\end{theorem}

\begin{proofof}[Proof of Theorem \ref{thm:upper-bound-LSL}]
In the same spirit as in the proof of Theorem \ref{version1:upperbound} we have the event of wrong matching denoted by $\Omega = \{\hat\pi \neq \pi^*\}$ with $\hat\pi = \bar\pi^{\textup{LSL}}$. Notice that
\begin{align}
    \Omega \in \bigcup_{\pi \neq \pi^*} \Omega_{\pi},
\end{align}
where
\begin{align}
    \Omega_{\pi} &= \Big\{\sum_{i=1}^n \log \| X_i - X^{\text{\tt\#}}_i\|^2 \ge \sum_{i=1}^n \log \| X_i - X^{\text{\tt\#}}_{\pi(i)}\|^2\Big\} \\
    &\subset \bigcup_{i=1}^n \bigcup_{j \in [m] \setminus \{i\}} \Big\{ \log \| X_i - X^{\text{\tt\#}}_i \|^2 \ge \log \| X_i -X^{\text{\tt\#}}_j\|^2\Big\}
\end{align}
For the same random variables $\zeta_1$ and $\zeta_2$ defined in \eqref{zetas} we can upper bound the term $\log \| X_i - X^{\text{\tt\#}}_i \|^2$ and lower bound $\log \| X_i -X^{\text{\tt\#}}_j\|^2$ as follows
\begin{align*}
    \| X_i - X^{\text{\tt\#}}_i \|^2 \le 2 \sigma_i^2(d + \sqrt{d} \zeta_2)
\end{align*}
and
\begin{align*}
    \| X_i - X^{\text{\tt\#}}_j \|^2 \ge \| \theta_i - {\theta}^\text{\tt\#}_j \|^2 + \sigma_{i, j}^2 (d - \sqrt{d} \zeta_2) - 2\zeta_1 \| \theta_i - {\theta}^\text{\tt\#}_j \| \sigma_{i,j}.
\end{align*}
Further, on the event $\Omega_1 = \{ \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \ge \zeta_1 \}$ from the last display we get that
\begin{align*}
    \frac{\| X_i - X^{\text{\tt\#}}_j \|^2}{\sigma_{i, j}^2} \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 - 2\zeta_1 \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) + d - \sqrt{d} \zeta_2
\end{align*}
with $\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})$ defined as
\begin{align*}
  \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) = \min_{\substack{i\in[n], j \in [m] \\ i \neq j}} \frac{\| \theta_i - {\theta}^\text{\tt\#}_j\|}{\sigma_{i,j}}.
\end{align*}
Therefore, combining these bounds we get that
\begin{align*}
    \Omega \cap \Omega_1 \subset \Big\{ 2\sigma_i^2 (d+\sqrt{d}\zeta_2) \ge \sigma_{i, j}^2 (\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 - 2\zeta_1 \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) + d - \sqrt{d} \zeta_2) \Big\},
\end{align*}
which implies
\begin{align*}
    \mathbf P(\Omega) &\le \mathbf P(\Omega_1^c) + \mathbf P(\Omega \cap \Omega_1) \\
    &\le \mathbf P(\zeta_1 \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})) + \mathbf P(\sqrt{d}\zeta_2(3\sigma_i^2 + \sigma_j^2) + 2\sigma_{i, j}^2 \zeta_1 \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) \ge \sigma_{i, j}^2 \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 + d(\sigma_j^2 - \sigma_i^2)) \\
    \label{prob-error}
    &\le 2\mathbf P(\zeta_1 \ge \frac 1 4 \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})) + \mathbf P(\zeta_2 \ge \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma}) A(\boldsymbol{\sigma}) + \sqrt{d} |B(\boldsymbol{\sigma})|),
\end{align*}
where
\begin{align*}
    A(\boldsymbol{\sigma}) = \frac{\sigma_i^2 + \sigma_j^2}{2(3\sigma_i^2 + \sigma_j^2)}, \qquad B(\boldsymbol{\sigma}) = \frac{\sigma_j^2 - \sigma_i^2}{3\sigma_i^2 + \sigma_j^2}.
\end{align*}
Notice that the term $A(\boldsymbol{\sigma})$ is bounded, namely it is easy to show that $A(\boldsymbol{\sigma}) \in (1/6, 1/2)$ for any combination of $\sigma_i$ and $\sigma_j$ involved in the definition of $A(\boldsymbol{\sigma})$. Similarly, it can be verified that $|B(\boldsymbol{\sigma})| \in [0, 1)$.
Hence, combining both bounds obtained from \eqref{prob-error} and taking $\bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})$ to be greater than $\bar\kappa$ we get that if
\begin{align*}
    \bar\kappa(\boldsymbol{\theta}, \boldsymbol{\sigma})^2 \ge \bar\kappa^2 \triangleq d/2 + \Big(\sqrt{2d \log \frac{8nm}{\alpha}} \vee 32\log\frac{4nm}{\alpha}\Big),
\end{align*}
then $\mathbf P(\hat\pi \neq \pi^*) \le \alpha$.
\end{proofof}

\subsubsection{Outlier detection}
In this section we show that the outlier detection accuracy depends only on the quantity $\bar\kappa_{\textup{out}}$. In Theorem \ref{version1:upperbound} we showed that if both $\bar\kappa_{\textup{in}}$ and $\bar\kappa_{\textup{out}}$ are large enough then the LSNS procedure finds the correct matching with high probability.

\begin{lemma}\label{lem:small-kappa}
Assume that for two distinct inliers $i \neq j \in [n]$ the probability of the mismatch between two pairs $(X_i, X^{\text{\tt\#}}_i)$ and $(X_j, X^{\text{\tt\#}}_j)$ is larger than some universal constant $c$. Considering the LSS procedure the last statement formally reads as
\begin{align}
    \mathbf P(\|X_i - X^{\text{\tt\#}}_j\|^2 + \| X_j - X^{\text{\tt\#}}_i \|^2 \le \| X_i - X^{\text{\tt\#}}_i \|^2 + \| X_j - X^{\text{\tt\#}}_j\|^2) \ge c.
\end{align}
Then,
\begin{align}
    \kappa_{i,j} \triangleq \frac{\|\theta_i - \theta_j\|}{\sqrt{\sigma_i^2 + \sigma_j^2}} \le 2\sqrt{\log \frac 2 c} + 1.5.
\end{align}
\end{lemma}

\todo{put $c=1/4$}
\begin{proofof}[Proof of Lemma \ref{lem:small-kappa}]
First notice that the condition
\begin{align*}
    \|X_i - X^{\text{\tt\#}}_j\|^2 + \| X_j - X^{\text{\tt\#}}_i \|^2 \le \| X_i - X^{\text{\tt\#}}_i \|^2 + \| X_j - X^{\text{\tt\#}}_j\|^2
\end{align*}
can be rewritten as follows
\begin{align}\label{diff-regrouped}
    (X_i - X_j)^\top (X^{\text{\tt\#}}_i - X^{\text{\tt\#}}_j) \le 0.
\end{align}
Under the assumption that $\pi^*(i) = i$ for all $i \in [n]$ the inequality from \eqref{diff-regrouped} reads as
\begin{align*}
    D_{ij} \triangleq \| \theta_i - \theta_j\|^2 + (\eta_{ij} + \eta^{\text{\tt\#}}_{ij})^\top(\theta_i - \theta_j) + \eta_{ij}^\top \eta^{\text{\tt\#}}_{ij} \le 0,
\end{align*}
where $\eta_{ij} = \sigma_i \xi_i - \sigma_j \xi_j$ and $\eta^{\text{\tt\#}}_{ij} = \sigma_i \xi^\text{\tt\#}_i - \sigma_j\xi^\text{\tt\#}_j$. We introduce the normalized version of $D_{ij}$, denoted by $\bar{D}_{ij}$:
\begin{align}\label{diff-normalized}
    \bar{D}_{ij} = \frac{D_{ij}}{\sigma_i^2 + \sigma_j^2} = \kappa_{ij}^2 + \frac{(\eta_{ij} + \eta^{\text{\tt\#}}_{ij})^\top(\theta_i - \theta_j)}{\sigma^2_{i, j}} + \frac{\eta_{ij}^\top \eta^{\text{\tt\#}}_{ij}}{\sigma_{i, j}^2},
\end{align}
where $\kappa_{ij} = \frac{\| \theta_i - \theta_j\|}{\sigma_{i, j}}$. Due to the independence of $\eta_{ij}$ and $\eta^{\text{\tt\#}}_{ij}$ we can decompose the sum in \eqref{diff-normalized} in the following way
\begin{align}
    \bar{D}_{ij} = \zeta_1 + \frac{\eta_{ij}^\top \eta^{\text{\tt\#}}_{ij}}{\sigma_{i, j}^2},
\end{align}
where $\zeta_1 \triangleq \kappa_{ij}^2 + \frac{(\eta_{ij} + \eta^{\text{\tt\#}}_{ij})^\top(\theta_i - \theta_j)}{\sigma^2_{i, j}} \sim \mathcal{N}(\kappa_{ij}^2, 2\kappa_{ij}^2)$. The second term from the last display can be dealt as follows
\begin{align*}
    \zeta_2 &\triangleq \frac{\eta_{ij}^\top \eta^{\text{\tt\#}}_{ij}}{\sigma_{i, j}^2} = \frac{1}{4\sigma_{i,j}^2} \Big[(\eta_{ij} + \eta^{\text{\tt\#}}_{ij})^2 - (\eta_{ij} - \eta^{\text{\tt\#}}_{ij})^2\Big] \\
    &\sim \frac 1 2 (Q - Q^{\texttt{\#}}),
\end{align*}
where $Q$ and $Q^{\texttt{\#}}$ are independent and $Q, Q^{\texttt{\#}} \sim \chi^2_1$. The independence of $Q$ and $Q^{\texttt{\#}}$ follows from the fact that $\mathbb{E}\eta_{12}^2 = \mathbb{E}[\eta^{\text{\tt\#}}_{12}]^2$.
Thus, on the event $\Omega \triangleq \Big\{ \zeta_2 \le \log \frac 1 \delta + 2 \sqrt{\log \frac 1 {\delta}}\Big\}$ we have
\begin{align}
    \mathbf P(\bar{D}_{12} \le 0) &\ge \mathbf P\Big(\zeta_1 + \log \frac 1 \delta + 2\sqrt{\log \frac 1 {\delta}} \le 0\Big) \\
    &= \mathbf P\Big(Z_0 \ge \frac{\kappa_{ij}^2 + \log \frac 1 \delta + 2\sqrt{\log \frac 1 {\delta}}}{\kappa_{ij}\sqrt{2}}\Big) \ge c,

\end{align}
where $Z_0$ is a standard normal random variable.
From the last insertion we derive the upper bound for $\kappa_{ij}$, by using the lower bound for the complementary error function.\footnote{proof or citation needed.}
\begin{align}
    e^{} \ge c.
\end{align}

\end{proofof}

\begin{remark}
Taking
\end{remark}

\subsection{Lower bounds for $\bar\kappa_{\textup{out}}$}

\section{Numerical results}

\begin{figure}
\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{outlier_results.png}
\caption{Results of the numerical experiments. Every point is a distinct experiment.}
\label{fig:results}
\end{figure}

We performed several numerical experiments to corroborate our theoretical results. All experiments were implemented using python (numpy). For linear sum assignment problem we used the generalized Hungarian algorithm
implemented in scipy library. The experiment was carried out in the following manner. First we generate $\pi^*$ (the original matching) and $\boldsymbol{\theta}^\texttt{\#}$ (note that $\boldsymbol{\theta}$ can be derived using $\boldsymbol{\theta}^\texttt{\#}$ and $\pi^*$). $\boldsymbol{\theta}^\texttt{\#}$ is sampled from Gaussian distribution with 0 mean and $\tau$ variance (which is also sampled randomly from uniform distribution on $[0, 2]$). Additionally, for every $\theta_i \in \boldsymbol{\theta}^\texttt{\#}$ for which $i \notin J_{\pi^*}$ ($\theta_i$ is an outlier) we increment every coordinate of $\theta_i$ by $i$. $\boldsymbol{\sigma}^\text{\tt\#}$ (also $\boldsymbol{\sigma}$) is sampled from uniform distribution over $[0.5, 2]$. Afterwards we generate $\boldsymbol{X}$ and $\boldsymbol{X}^{\text{\tt\#}}$ according to the procedure described in section \ref{section:model}. We try to solve the problem with four aforementioned algorithms (Greedy, LSS, LSNS, LSL). The results are summed up in Figure \ref{fig:results}. In all experiments $n = 100, m = 130$ and $d = 50$.
\todo{Add explanation of obtained results.}\todo{Average the plots over 50 independent trials.}

\section*{Acknowledgments}
This work was partially supported by the grants Investissements d'Avenir (ANR-11-IDEX-0003/Labex Ecodec/ANR-11-LABX-0047) and
CALLISTO. The authors thank the Reviewers for many valuable suggestions.

\part{}

\maketitle
\vspace{-5pt}
\begin{abstract}

    We consider the problem of finding the matching map between two sets of $d$-dimensional vectors from noisy observations, where the second set contains outliers.

    The matching map is then an injection, which can be
    consistently estimated only if the vectors of
    the second set are well separated.
    The main result shows that,

    in the high-dimensional setting,
    a detection region of unknown injection may be characterized by the sets
    of vectors for which the inlier-inlier distance is
    of order at least $d^{1/4}$ and the inlier-outlier
    distance is of order at least $d^{1/2}$.

    These rates are achieved using the estimated matching minimizing the sum of logarithms of distances between matched pairs of points.

    We also prove lower bounds establishing optimality of these rates.
    Finally, we report the results of numerical
    experiments on both synthetic and real world data that illustrate our theoretical results
    and provide further insight into the properties
    of the estimators studied in this work.
\end{abstract}

\vspace{-5pt}
\section{Introduction}\label{sec:intro}

Finding the best match between two clouds of
points is a problem encountered in many real
problems. In computer vision, one can look for
correspondences between two sets of local
descriptors extracted from two images. In
text analysis, one can be interested in matching
vector representations of the words of two
similar texts, potentially in two different
languages. The goal of the present work is to
gain theoretical understanding of the statistical
limits of the matching problem.

In the sequel, we use the notation $[n] = \{1,\ldots,
n\}$ for any integer $n$, and define $\|\cdot\|$
as the Euclidean norm in $\mathbb{R}^d$. Assume
that two independent sequences $\boldsymbol{X} = (X_i;
i\in[n])$ and $\boldsymbol{Y} =(Y_i; i\in[n])$ of independent
vectors are generated such that $X_i$ and $Y_i$
are drawn from the same distribution $P_i$
on $\mathbb{R}^d$, for every $i\in[n]$. The
statistician observes the sequence $\boldsymbol{X}$
and a shuffled version
$\boldsymbol{X}^{\text{\tt\#}}$ of the sequence $\boldsymbol{Y}$. More precisely,
$\boldsymbol{X}^{\text{\tt\#}}$ is such that $X^{\text{\tt\#}}_i = Y_{\pi^*(i)}$
for some unobserved permutation $\pi^*$. The
goal of matching is to estimate the permutation
$\pi^*$ from data $(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}})$. In the case
of Gaussian distributions $P_i$, this problem
has been studied in \citep{jmlr_CD13, collier2016minimax}. Clearly, consistent estimation
of the matching map $\pi^*$ is impossible if
there are two data generating distributions
$P_i$ and $P_j$ that are very close. In
\citep{jmlr_CD13, collier2016minimax}, a
precise quantification of the separation
between these distributions is given that
enables consistent estimation of $\pi^*$.
Furthermore, it is shown that the permutation
estimator minimizing the sum of logarithms
of pairwise distances between the elements of
$\boldsymbol{X}$ and the elements of the shuffled version $\boldsymbol{X}^{\text{\tt\#}}$ is an optimal estimator of $\pi^*$.

In this paper, we extend the model studied in \citep{collier2016minimax} to the case when the
set $\boldsymbol{X}^{\text{\tt\#}}$ is contaminated by outliers. The number of outliers is supposed to be known and is equal to $m - n$, where $n = |\boldsymbol{X}|$ and $m = |\boldsymbol{X}^{\text{\tt\#}}|$ are the sizes of considered
two sequences, however the indices of the
outliers are unknown. All the distributions
are assumed in this paper to be
spherical Gaussian, although all the probabilistic
tools used in the proofs have their sub-Gaussian
counterparts. Thus, we consider that two spherical
Gaussian distributions
\footnote{We use the notation $\mathbf I_d$ for the $d\times d$ identity matrix}
$P_1 = \mathcal N_d(\mu_1, \sigma_1^2\mathbf I_d)$ and $P_2 = \mathcal
N_d(\mu_2, \sigma_2^2\mathbf I_d)$ are well
separated if the ``distance to noise ratio''
$\kappa(P_1,P_2) = \|\mu_1-\mu_2\|/\sqrt{
\sigma_1^2+\sigma_2^2}$ is large. Main findings
of \citep{collier2016minimax}, in terms of
smallest separation distance $\bar\kappa =
\min_{i\neq j}\kappa(P_i,P_j)$ are summarized
in the second columns of \Cref{tab:1}.
Likewise, the last column of the table
provides a summary of the contributions of
the present paper in terms of $\bar\k_{\textup{in-in}} =
\min_{i\neq j} \kappa(P_i,P_j)$ and $\bar\k_{\textup{in-out}} =
\min_{i,j} \kappa(P_i,Q_j)$, where $Q_1,\ldots,
Q_{m-n}$ are the distributions of the outliers.

An unexpected finding of this work is that the
``degree'' of heteroscedasticity of the model
has a strong impact on the separation distances
and the detection regions (sets of values of
$(\bar\k_{\textup{in-in}},\bar\k_{\textup{in-out}})$ for which it is possible to
detect the feature map $\pi^*$). This is in sharp
contrast with the outlier-free case, where consistent
estimation requires $\bar\kappa$ to be
at least of order $(d\log n)^{1/4}$ irrespective
from the behaviour of variances of $P_i$.
We prove in this work that in the high dimensional
regime $d\ge c\log n$, which is arguably more
appealing than the low dimensional regime $d\le c\log n$,
the following statements are true:
\begin{itemize}
    \item  If there is no heteroscedasticity, \textit{i.e.}, when all the variances are equal, consistent estimation
    of $\pi^*$ is possible if and only if $\bar\kappa
    = \bar\k_{\textup{in-in}}\wedge\bar\k_{\textup{in-out}} $
    is at least of order $(d\log(nm))^{1/4}$. This is
    the same rate as in the outlier-free case.
    \item If the heteroscedasticity is mild, \textit{i.e.},
    all the variances are of the same order, the
    condition
    $\bar\k_{\textup{in-in}}\gtrsim (d\log(nm))^{1/4}$ is the same as in the
    previous item, but the stronger condition $\bar\k_{\textup{in-out}}
    \gtrsim d^{1/2}$ is needed for the inlier-outlier
    separation distance.
    \item Finally, in the general heteroscedastic setting
    both $\bar\k_{\textup{in-in}}$ and $\bar\k_{\textup{in-out}}$ should be at least of
    order $d^{1/2}$. Furthermore, in all these cases
    consistent estimation is performed by the same
    estimator: the Least Sum of Logarithms (LSL).
\end{itemize}
Note also that the empirical evaluation done in this
paper shows that LSL is interesting not only from the
theoretical but also from the practical point of view.

\begin{table}
    \centering
    \begin{tabular}{ c|c|c }
        \toprule
         & Without outliers  & With outliers \\
         &\citep{collier2016minimax}& this paper \\
        \toprule
        \shortstack{known $\boldsymbol{\sigma}^\text{\tt\#}$ or \\
        all equal $\boldsymbol{\sigma}^\text{\tt\#}$s} & \shortstack{
        LSNS is optimal \\ $\bar\kappa \gtrsim (d\log n)^{1/4}$}
        & \shortstack{ LSNS is optimal [Thm. 1] \\
        $\bar\k_{\textup{in-in}} \wedge \bar\k_{\textup{in-out}}\gtrsim (d
        \log(nm))^{1/4}$} \\
        \midrule
        \begin{tabular}{c}
        \shortstack{unknown $\boldsymbol{\sigma}^\text{\tt\#}$ \\
        $\nicefrac{\sigma_{\max}}{\sigma_{\min}} \le C$} \\
        \midrule
        \shortstack{unknown $\boldsymbol{\sigma}^\text{\tt\#}$ \\
        arbitrary}
        \end{tabular}
        & \shortstack{LSL is optimal \\[5pt] $\bar\kappa
        \gtrsim (d\log n)^{1/4}$} &
        \begin{tabular}{c}
            \shortstack{LSL is optimal [Thm. 4] \\
            $\bar\k_{\textup{in-in}} \gtrsim (d\log(nm))^{1/4}$  \&
            $\bar\k_{\textup{in-out}}\gtrsim d^{1/2}$}  \\
            \midrule
            \shortstack{ LSL is optimal
            [Thm. 2, 3] \\ $\bar\k_{\textup{in-in}} \wedge
            \bar\k_{\textup{in-out}} \gtrsim d^{1/2}$}
        \end{tabular}
        \\
        \bottomrule
    \end{tabular}
    \vspace{-0.1cm}
    \caption{A brief overview of the contributions
    in the high-dimensional regime $d\ge c\log n$.
    The table provides the condition on the normalized
    inlier-inlier distance $\bar\k_{\textup{in-in}}$ and inlier-outlier
    distance $\bar\k_{\textup{in-out}}$, making it possible to consistently
    detect the matching map between two sets of
    $d$-dimensional vectors. LSL and LSNS refer to
    least sum of logarithms and least sum of normalized
    squares, respectively.
    }
    \vspace{-19pt}
    \label{tab:1}
\end{table}

\paragraph{Agenda} \Cref{sec:model} describes the
framework of the vector matching problem and
introduces the ter\-minology
used throughout this paper. Precise statements
of the main theoretical results are gathered in \Cref{sec:main}. The prior work is briefly discussed in
\Cref{sec:prior}. \Cref{sec:numerical} contains numerical experiments carried out both for synthetic and real data. A brief summary and some concluding remarks are presented in \Cref{sec:conc}. Proofs of all theoretical claims are deferred to the supplemental material.

\section{Problem Formulation}\label{sec:model}

We begin with formalizing the problem of matching two
sequences of feature vectors $(X_1,\ldots,X_n)$ and
$(X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m)$ with different sizes $n$
and $m$ such that $m \ge n \ge 2$. In what follows, we
assume that the observed feature vectors are randomly
generated from the model
\begin{equation}\label{model}
\begin{cases}
X_i = \theta_i + \sigma_i\xi_i \ , \\
X^{\text{\tt\#}}_j = {\theta}^\text{\tt\#}_j + \s^{\text{\tt\#}}_j\xi^\text{\tt\#}_j,
\end{cases}\quad i=1,\ldots,n \text{ and }
j=1,\ldots,m.
\end{equation}
In this model, illustrated in \Cref{fig:ani}, it
is assumed that
\vspace{-5pt}
\begin{itemize}\setlength{\itemsep}{-2pt}
\item $\boldsymbol{\theta}=(\theta_1,\ldots,\theta_n)$ and $\boldsymbol{\theta}^\text{\tt\#} =
({\theta}^\text{\tt\#}_1,\ldots,{\theta}^\text{\tt\#}_m)$ are two sequences of
vectors from $\mathbb R^d$, corresponding to the original
features, which are unavailable,
\item $\boldsymbol{\sigma}=(\sigma_1,\ldots,\sigma_n)^\top, \boldsymbol{\sigma}^\text{\tt\#}
= (\s^{\text{\tt\#}}_1,\ldots,\s^{\text{\tt\#}}_m)^\top$ are positive real
numbers corresponding to the magnitudes of the noise
contaminating each feature,
\item $\xi_1, \ldots, \xi_n$ and $\xi^\text{\tt\#}_1, \ldots,
\xi^\text{\tt\#}_m$ are two independent sequences of i.i.d.\
random vectors drawn from the Gaussian distribution
with zero mean and identity covariance matrix.
\end{itemize}

The simplest special case of \eqref{model},
considered in \citep{collier2016minimax}, corresponds
to the situation where a perfect matching exists
between the two sequences $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^\text{\tt\#}$.
This means that $m = n$ and, for some bijective
mapping $\pi^*:[n]\to[n]$, $\theta_i = {\theta}^\text{\tt\#}_{
\pi(i)}$ for all $i\in[n]$. In the general
case, both $\boldsymbol{X} = (X_1, \dots, X_n)$ and
$\boldsymbol{X}^{\text{\tt\#}} = (X^{\text{\tt\#}}_1, \dots, X^{\text{\tt\#}}_m)$ may
contain outliers, \textit{i.e.} feature vectors
that have no corresponding pair. In such a
situation, it is merely assumed that there exists
a set $S\subset [n]$ and an  injective mapping
$\pi^*:S\to [m]$ such that
\begin{align}\label{pistar}
    \theta_i = {\theta}^\text{\tt\#}_{\pi^*(i)} \quad\text
    {and}\quad \sigma_i = \s^{\text{\tt\#}}_{\pi^*(i)},
    \qquad \forall\,
    i\in S.
\end{align}
In this case we say that the vectors $\{X_i:i\in[n]\setminus S\}$ and $\{X^{\text{\tt\#}}_j:j\in[m]\setminus\pi^*(S)\}$ are
outliers. The ultimate goal is to estimate the
feature matching map $\pi^*$.

In this work we consider the case when $S = [n]$ and $m>n$. This means
that only the larger set of feature vectors, namely
$\boldsymbol{X}^{\text{\tt\#}}$, contains outliers. Let us also define the
set $O_{\pi^*} \triangleq [m] \setminus
\text{Im}(\pi^*)$,  which contains the indices of
outliers and satisfies  $|O_{\pi^*}| = m - n$.
Naturally, the feature vectors contained in $\boldsymbol{X}$,
as well as those vectors from $\boldsymbol{X}^{\text{\tt\#}}$ that are
not outliers, are called {\it{inliers}}.

In this formulation, the data generating distribution is defined by the parameters
$\boldsymbol{\theta}^\texttt{\#}$,
$\boldsymbol{\sigma}^\text{\tt\#}$ and $\pi^*$. We omit the set of parameters $\boldsymbol{\theta}$ and $\boldsymbol{\sigma}$, since they are automatically identified using $\pi^*$, $\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ by the formula $(\theta_i, \sigma_i) = ({\theta}^\text{\tt\#}_{\pi^*(i)},\s^{\text{\tt\#}}_{\pi^*(i)})$ for $i\in[n]$. Since our goal is to match the
feature vectors, we focus our attention on the problem of
estimating the parameter $\pi^*$ only, considering
$\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ as nuisance parameters.
In what follows, we denote by $\mathbf P_{ \boldsymbol{\theta}^\texttt{\#},
\boldsymbol{\sigma}^\text{\tt\#}, \pi^*}$ the probability distribution of the
sequence $(X_1,\ldots,X_n,X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m)$
defined by (\ref{model}) under condition (\ref{pistar})
with $S= [n]$.

We are interested in designing estimators that have an expected error smaller than a prescribed level $\alpha$ under the weakest possible conditions on the nuisance parameter $\boldsymbol{\theta}^\texttt{\#}$ and noise level $\boldsymbol{\sigma}^\text{\tt\#}$. Clearly, the problem of matching becomes more difficult with hardly distinguishable features. To quantify this phenomenon, we introduce the normalized separation distance $\bar\k_{\textup{in-in}} = \bar\k_{\textup{in-in}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ and the normalized outlier separation distance $\bar\k_{\textup{in-out}} = \bar\k_{\textup{in-out}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$, which measure the minimal distance-to-noise ratio between inliers and the minimal distance-to-noise ratio between inliers and outliers, respectively. The precise
definitions read as
\begin{equation}
\bar\k_{\textup{in-in}} \triangleq \min_{\substack{ i,j\not\in O_{\pi^*}, \\ j \not=i}} \frac{\|{\theta}^\text{\tt\#}_i-{\theta}^\text{\tt\#}_j\|}{(\s^{\text{\tt\#}}_i{}^2+\s^{\text{\tt\#}}_j{}^2)^{1/2}},\qquad
\bar\k_{\textup{in-out}} \triangleq \min_{\substack{ i\not\in O_{\pi^*}, \\ j \in O_{\pi^*}}} \frac{\|{\theta}^\text{\tt\#}_i-{\theta}^\text{\tt\#}_j\|}{(\s^{\text{\tt\#}}_i{}^2+\s^{\text{\tt\#}}_j{}^2)^{1/2}}. \label{kappa}
\end{equation}
Notice that $\bar\k_{\textup{in-in}}$ can be rewritten as
\begin{align}
    \bar\k_{\textup{in-in}} = \min_{\substack{i, j \in [n] \\ i \neq j}} \frac{\| \theta_i - \theta_j\|}{(\sigma_i^2 + \sigma_j^2)^{1/2}}.
\end{align}

\begin{figure}
\centering
\includegraphics[height=0.32\textwidth]{
ani_match_crop2.pdf}
\hspace{20pt}
\vspace{-11pt}
\includegraphics[height=0.32\textwidth]{
ani_match_crop1.pdf}
\caption{Illustration of the considered framework described in \eqref{model}.
We wish to match a set of 7 patches extracted from
the first image to the 9 patches from the second
image. The picture on the left shows the locations
of patches as well as the true matching map $\pi^*$
(the yellow lines).
}
\vspace{-13pt}
\label{fig:ani}
\end{figure}

Clearly, if $\bar\k_{\textup{in-in}}=0$ or, $\bar\k_{\textup{in-out}}=0$, there are
two identical feature vectors in $\boldsymbol{X}^{\text{\tt\#}}$. In such
a situation, assuming $\sigma_i$'s are all equal,
the parameter $\pi^*$ is nonidentifiable, in the
sense that there exist two different permutations
$\pi^*_1$ and $\pi^*_2$ such that the distributions
$\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_1}$ and
$\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_2}$
coincide.  Therefore, to ensure the existence of
consistent estimators of $\pi^*$ it is necessary
to impose the conditions $\bar\k_{\textup{in-in}}>0$ and
$\bar\k_{\textup{in-out}}>0$. Moreover, good estimators are those
consistently estimating $\pi^*$ even if either
$\bar\k_{\textup{in-in}}$ or $\bar\k_{\textup{in-out}}$ are small. We are interested
here in finding the detection boundary in terms
of the order of magnitude of $(\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}})
$. More precisely, for any given $\alpha\in(0,1)$
we wish to find a region $\mathcal R_{n,m,d}^\alpha$
in $\mathbb R^2$ such that:
\vspace{-8pt}
\setlength{\leftmargini}{2em}
\begin{itemize}
    \item There is an estimator $\hat\pi_{n,m}$ of
    $\pi^*$ satisfying

        $\mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi
        \neq \pi^*) \le \alpha$

    for every $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$
    lying in the detection region, \textit{i.e.}, for
    which $(\bar\k_{\textup{in-in}}, \bar\kappa_{
    \textup{in-out}}) \in \mathcal R_{n,m,d}^\alpha$.
    \item There is a constant $C<1$ such that for any
    estimator $\bar\pi_{n,m}$ of $\pi^*$, we can find
    a parameter value $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$
    in the region $\{(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*) :
    (\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}})
    \in C\mathcal R_{n,m,d}^\alpha\}$
    such that $\bar\pi$ fails to detect $\pi^*$ with a
    probability larger than $\alpha$.
\end{itemize}
\vspace{-4pt}
Let us make two remarks. First, note that in the
outlier-free case considered
in \citep{collier2016minimax}, $\bar\k_{\textup{in-out}}$
is meaningless and, therefore, the detection region  is one-dimensional $\mathcal R_{n,m,d}^\alpha$. Thus, it
is necessarily a half-line and is proven to be of the
form $\bar\k_{\textup{in-in}}\ge C(\log \nicefrac{
n}{\alpha})^{1/2} \vee (d\log \nicefrac{n}{\alpha} )^{1/4}$
for some universal constant $C$. Second, the aforementioned definition of the detection region $\mathcal R_{n,m,d}^\alpha$ does not guarantee
its uniqueness (even up to a scaling by a universal
constant). This is in contrast with the
outlier-free case. To overcome this difficulty, we
look for $\mathcal R_{n,m,d}^\alpha$ of
the form $[t_{\textup{in-in}},+\infty)\times [t_{\textup{in-out}},+\infty)$ with the smallest possible
threshold $t_{\textup{in-out}}$ for the normalized inlier-outlier distance $\bar\k_{\textup{in-out}}$.

\def\boldsymbol{\Xi}{\boldsymbol{\Xi}}
\section{Main theoretical results}\label{sec:main}

In this section, we have collected the main theoretical findings of the paper.  When the noise is homoscedastic,
\textit{i.e.},  when all $\boldsymbol{\sigma}$'s are equal,
the results obtained by \cite{collier2016minimax} in the outlier-free setting can be easily extended to the setting
with outliers. Therefore, in the present paper, we focus
on the heteroscedastic case. For the sake of clarity of
exposition, we will present the results in the case of known
variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$ prior to investigating
the more interesting case of unknown variances.

The detection regions we study below are based on the
profile maximum likelihood estimator. The model presented in \eqref{model} has the parameter $\boldsymbol{\Xi} = (\boldsymbol{\theta}^\text{\tt\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi)$, while the observations are the sequences
of feature vectors $\boldsymbol{X}$ and $\boldsymbol{X}^{\text{\tt\#}}$. The negative
log-likelihood of this model is given by
\begin{align}
    \ell_n(\boldsymbol{\Xi}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\}) &= \sum_{i=1}^n \bigg(\frac{\| X_i - {\theta}^\text{\tt\#}_{\pi(i)}\|_2^2}{2\s^{\text{\tt\#}2}_{\pi(i)}} + \frac 1 2 \log(\s^{\text{\tt\#}2}_{\pi(i)})\bigg) +  \sum_{j=1}^m \bigg(\frac{\| X^{\text{\tt\#}}_j - {\theta}^\text{\tt\#}_j\|_2^2}{2\s^{\text{\tt\#}2}_j} + \frac 1 2 \log(\s^{\text{\tt\#}2}_j)\bigg).
\end{align}
The profile negative log-likelihood is then defined as
the minimum
with respect to $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#})$ of
the log-likelihood $\ell_n(\boldsymbol{\Xi}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\})$.

\subsection{Warming up: known variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

One can check
that the minimization with respect to $\boldsymbol{\theta}^\texttt{\#}$ leads
to the variance-dependent cost function
\begin{align}\label{liklihood-sigma}
    \ell_n(\pi,\boldsymbol{\sigma}^\text{\tt\#}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\}) =  \sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2} + \sum_{i=1}^n \frac{1}{2}\log(\s^{\text{\tt\#}2}_{\pi(i)}) + \sum_{j=1}^m \frac{1}{2} \log(\s^{\text{\tt\#}2}_j).
\end{align}
When $m=n$ and there is no outlier, the last two sums of the
last display do not depend on $\pi$ and, therefore, the
maximum profile likelihood estimator of $\pi^*$ is obtained
by the Least Sum of Normalized Squares (LSNS) criterion
\begin{align}\label{LSNS}
    \hat\pi^{\textrm{LSNS}}_{n,m}\in \text{arg}\min_{\pi}
    \sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2},
\end{align}
where the minimum is over all injective mappings
$\pi:[n]\to[m]$. This, and the other estimators defined
in this work, can be efficiently computed using suitable
versions of the Hungarian algorithm \citep{Kuhn,Kuhn2010,
Munkres}. As shows the next theorem, it
turns out that even when $m>n$, the estimator
$\hat\pi^{\textrm{LSNS}}_{ n,m } $  defined above
leads to an optimal detection region.

\begin{theorem}[Upper bound for LSNS]\label{thm:upperLSNS}
Let $\alpha \in (0,1)$ and condition \eqref{pistar}
be fulfilled. If the separation distances $\bar\k_{\textup{in-in}}$
and $\bar\k_{\textup{in-out}}$ corresponding to $(\boldsymbol{\theta}^\texttt{\#},
\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$ and defined in \eqref{kappa}
satisfy the condition
\begin{align}\label{thresh:1}
    \min\{\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}}\} \ge 4 \Big\{ \big(
    d\log(\nicefrac{4nm}{\alpha})\big)^{1/4} \vee
    \big( 2\log(\nicefrac{8nm}{ \alpha})\big)^{1/2}
    \Big\}
\end{align}
then the LSNS estimator defined in \eqref{LSNS} detects
the true matching map $\pi^*$ with probability at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi^{
    \textup{LSNS}}_{n,m} = \pi^*) \ge 1-\alpha.
\end{align}
\vspace{-15pt}
\end{theorem}

The similarity---both its statement and its proof---
of this result to its counterpart
in the outlier-free setting might suggest that
the presence of outliers does not make the problem
any harder from a statistical point of view.
However, this is not true in the more
appealing setting of unknown variances.

\subsection{Detection of $\pi^*$ for unknown and
arbitrary variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

The LSNS procedure analyzed in \Cref{thm:upperLSNS}
exploits the values of known noise variances to normalize
the squares of distances between vectors $X_i$ and
$X^{\text{\tt\#}}_{\pi(i)}$. Therefore, LSNS is inapplicable
in the case of unknown noise variances.
Instead, we consider the Least Sum of
Logarithms (LSL) estimator
\begin{align}\label{LSL}
    \hat\pi^{\textup{LSL}}_{n,m} \triangleq \arg\min_{\pi:[n]\to [m]}\sum_{i=1}^n {\log \|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2},
\end{align}
where the minimum is over all injective maps $\pi
: [n] \to [m]$. This estimator can be seen as the
minimizer of a criterion defined as the minimum of
the cost function from \eqref{liklihood-sigma} with
respect to $\boldsymbol{\sigma}^\text{\tt\#}$ under the constraint
$\min_{j\not\in\Im(\pi)} \s^{\text{\tt\#}}_{j} \ge
\sigma_{\min}$, for some fixed (but unknown)
constant  $\sigma_{\min}>0$.

To provide a quick overview of what follows, let us
stick in the remaining of this paragraph to the
case $\log(nm) = O(d)$ so that the right hand side
of \eqref{thresh:1} is of order $\big( d\log(nm)
\big){}^{1/4}$. Recall that in the outlier-free
case, the LSL estimator has been shown to perform
as well as the LSNS while having the advantage of
not requiring the knowledge of variances
$\boldsymbol{\sigma}^\text{\tt\#}$ \citep{collier2016minimax}.
Somewhat unexpectedly, the situation is significantly
different in the presence of outliers.  Indeed, the
best we manage to prove in the presence of
outliers is that the detection of the matching
map by LSL is possible whenever $\min\{ \bar\kappa_{
\textup{in-in}}, \bar\k_{\textup{in-out}}\}\ge
C\sqrt{d}$ for some sufficiently large constant
$C$. The precise statement being given in the next
theorem, let us mention right away that the
discrepancy between this rate $\sqrt{d}$ and the
rate $\big( d\log(nm) \big){}^{1/4}$ in
\eqref{thresh:1} is due to the inherent hardness
of the setting and not merely an artefact of the
proof. This will be made clear below.

\begin{theorem}[Upper bound for LSL]
\label{thm:upper-bound-LSL}
Let $\alpha \in (0,1/2)$ and condition
\eqref{pistar} be fulfilled. If the separation
distances $\bar\k_{\textup{in-in}}$ and $\bar\k_{\textup{in-out}}$ corresponding
to $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$ and
defined by \eqref{kappa} satisfy
\begin{align}\label{detect}
    \min\{\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}}\} \ge \sqrt{2d} + 4
    \Big\{\big(2d\log (\nicefrac{4nm}{ \alpha})
    \big)^{1/4} \vee \big(3 \log(\nicefrac{8nm}{
    \alpha}\big)^{1/2}\Big\}
\end{align}
then the LSL estimator \eqref{LSL} detects
the matching map $\pi^*$ with probability at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}
    (\hat\pi^{\textup{LSL}}_{n,m} = \pi^*)
    \ge 1 - \alpha.
\end{align}
\vspace{-15pt}
\end{theorem}

This result is disappointing since it requires
the distance between different feature vectors to
be larger than $\sqrt{2d}$ in order to be able
to consistently estimate the matching map $\pi^*$.
As we show below, without any further condition
(for instance, on the noise variances), this
rate cannot be improved. Furthermore, the
rate $\sqrt{d}$ is optimal not only for LSL but
also for the larger class of so called
\textit{distance based $M$-estimators.}

We will say that an estimator $\hat\pi_n$ of $\pi^*$
is a distance based $M$-estimator, if for a sequence of non-decreasing functions $\rho_i:\mathbb{R}_+
\to \mathbb{R}$, $i=1,\ldots,n$, the following is correct
\begin{align}
    \hat\pi_n \in \arg\min_{\pi :[n] \to [m]} \sum_{i=1}^n \rho_i\big(\| X_i - X^{\text{\tt\#}}_{\pi(i)}\|\big),
\end{align}
where the minimum is over all injective mappings
$\pi:[n]\to[m]$. We denote by $\mathcal M$ the set
of all distance based $M$-estimators.
We show that there is indeed a setup where $\bar\k_{\textup{in-in}} \wedge \bar\kappa_{ \textup{in-out}}$ is as large as $0.2\sqrt{d}$ but
any estimator from $\mathcal M$ fails to detect
$\pi^*$ with probability at least $1/4$.  The next theorem formalizes the described result.
\begin{theorem}[Lower bound]\label{thm:lower-bound}
Assume that $m > n \ge 4$ and $d \ge 422\log(4n)$. There exists a triplet $(\boldsymbol{\sigma}^\text{\tt\#}, \boldsymbol{\theta}^\text{\tt\#}, \pi^*)$ such that  $\bar\k_{\textup{in-in}} = \bar\k_{\textup{in-out}} = \sqrt{d/20}$ and
    \begin{align}\label{lower:1}
        \inf_{\hat\pi \in \mathcal{M}}

        \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi
        \neq \pi^*) > \nicefrac 1 4.
    \end{align}
    \vspace{-15pt}
\end{theorem}

The proof of this theorem, postponed to the appendix,
is constructive. This means that we exhibit a triplet
$(\boldsymbol{\sigma}^\text{\tt\#}, \boldsymbol{\theta}^\text{\tt\#}, \pi^*)$ satisfying
\eqref{lower:1}. Careful inspection shows that in the
case $d = O(\log(nm))$ the same triplet satisfies
$\bar\k_{\textup{in-in}} \wedge \bar\kappa_{ \textup{in-out}} \asymp
\sqrt{\log(nm)}$ and \eqref{lower:1} is still true.
This implies that the order of magnitude of the right
hand side of \eqref{detect} is optimal both in the
high-dimensional regime $d\ge 422\log(4n) $ and in the
low-dimensional regime $d< 422\log(4n)$. This shows the optimality of LSL among all estimators
from $\mathcal{M}$. Note that the estimator $\hat\pi^{
\textup{LSNS}}_{n,m}$ does not belong to the family of
\textit{distance based M-estimators}. Furthermore, in the
low dimensional regime $d=O(\log(nm))$, the separation
rate of the LSL, $\sqrt{\log(nm)}$, is the same as that of the LSNS.

In the next section we show that under some mild
conditions on $\boldsymbol{\sigma}^\text{\tt\#}$ it is indeed possible to
obtain different rates for $\bar\k_{\textup{in-in}}$
and $\bar\kappa_{ \textup{in-out}}$, namely we show that
if $\bar\k_{\textup{in-in}} \gtrsim d^{1/4}$ and
$\bar\k_{\textup{in-out}} \gtrsim d^{1/2}$ then the
LSL estimator detects correct matching with high
probability.

\subsection{Detection of $\pi^*$ for unknown and mildly
varying variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

The results of the last two theorems are disappointing,
since they indicate that the features should be very
different from one another for detection of the matching
map to be possible. An interesting finding, presented
below, is that strong constraint can be significantly
alleviated in the context of mild heteroscedasticity.
By mild heteroscedasticity we understand here the
situation in which all variances $\s^{\text{\tt\#}}_i$ are of
the same order of magnitude.

\begin{theorem}[Upper bound under mild
heteroscedasticity]\label{thm:mild}
Let $r_\sigma = \max_{i,j}(\s^{\text{\tt\#}}_i/\s^{\text{\tt\#}}_j)
<\infty$. If the separation distances $\bar\k_{\textup{in-in}}$
and $\bar\k_{\textup{in-out}}$ defined in \eqref{kappa} satisfy
\begin{align}
    \bar\k_{\textup{in-in}} &\ge
    2\big(4d\log(\nicefrac{4nm}{\alpha})
    \big)^{1/4} + 2\big(2\log(\nicefrac{4nm}{
    \alpha})\big)^{1/2}\\
    \bar\k_{\textup{in-out}}
        &\ge\  \sqrt{2(r_\sigma-1)d} + 2\big(4
        r_\sigma^2 d\log(\nicefrac{4nm}{\alpha}
        ) \big)^{1/4} + 2\big(2 r_\sigma \log(
        \nicefrac{4nm}{\alpha}) \big)^{1/2} ,
\end{align}
then the LSL estimator \eqref{LSL}
detects the matching map $\pi^*$ with probability
at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}
    (\hat\pi^{\textup{LSL}}_{n,m} = \pi^*) \ge
    1 - \alpha.
\end{align}
\vspace{-15pt}
\end{theorem}
Note that a lower bound similar to that of
\Cref{thm:lower-bound} can be proved in the case
of mild heterescodestacity as well, showing that
there is an example for which $\bar\k_{\textup{in-in}}$ is of order
$d^{1/4}$, $\bar\k_{\textup{in-out}}$ is of order $d^{1/2}$ and
any estimator from $\mathcal M$ fails to detect
$\pi^*$ with probability at least $1/4$.

We complete this section by summarizing the
joint contribution of
\Cref{thm:upperLSNS,thm:upper-bound-LSL,thm:lower-bound,thm:mild}.
To simplify this discussion, we consider two
cases: high-dimensional case refers to $d\ge \log(4nm/
\alpha)$ (presented in \Cref{tab:1}) and low-dimensional case refers to
the condition $d< \log(4nm/\alpha)$. In the
high dimensional setting with arbitrary noise
variances, the detection region for the LSL
estimator is given by $\{\bar\k_{\textup{in-in}}\wedge \bar\k_{\textup{in-out}}
\ge 15 \sqrt{d}\}$, which is much worse than
the detection region for LSNS, $\{\bar\k_{\textup{in-in}}
\wedge \bar\k_{\textup{in-out}}\ge 8 (d\log(4nm/\alpha))^{1/4}
\}$, obtained in the known-variance scenario.
Somewhat surprisingly, in such a setting,
even a strong assumption on the outliers, such
as requiring them to be at least at a
distance $0.2\sqrt{d}$ of the inliers, is not
enough for relaxing the assumption on the
inlier-inlier separation distance. Finally, on a positive
note, in the intermediate case of mildly
varying variances, the detection region for
the LSL estimator is of the form
$\{\bar\k_{\textup{in-in}}\ge 7 (d\log(4nm/\alpha))^{1/4}; \
\bar\k_{\textup{in-out}} \ge 10 \sqrt{d}\}$. This means that
if the outliers are at a distance $\Omega(
\sqrt{d})$ of the inliers, then the LSL
recovers the true matching under the same
condition on $\bar\k_{\textup{in-in}}$ as in the outlier-free
setting.
}

\begin{figure
\begin{center}
\includegraphics[width=0.85\textwidth]{figure1c.png}
\vspace{-10pt}
\caption{The performance of the methods (Greedy, LSS,
LSL, LSNS) in the setup described in Exp. 1.
Plots show the error rate (percentage of repetitions
in which the estimated matching map differs from the
true one) as a function of separation distances. The
left plot illustrates that LSS, LSL and LSNS require
much lower value of $\bar\k_{\textup{in-in}}$ in order to find the
correct mapping. The left plot shows that there is
a clear improvement of the error when the inlier-inlier
separation distance increases, while the right plot shows
that the error rate might not be made small by
augmenting $\bar\k_{\textup{in-out}}$ only. }
\label{fig:exp1}
\end{center}
\vspace{-20pt}
\end{figure}

\section{Other related work}\label{sec:prior}

Measuring the quality of the various statistical
procedures of decision making by their minimal
separation rates became the standard in hypotheses
testing, see the seminal papers \citep{Burnashev,
Ingster82} and the monographs \citep{Ingster,JudNem}.
Currently this approach is widely adopted in machine
learning literature \citep{Xing20,Wolfer,Balchard,
RamdasISW16,Wei,Collier12a}. Beyond the classical
setting of two hypotheses, it can also be applied
to multiple hypotheses testing framework, for instance,
variable selection \citep{Ndaoud,azais2020multiple,
ComDal1} or the matching problem considered here.

On the other hand, feature matching is a well
studied problem in computer vision. In recent
years, a great deal of attention was devoted to
the acceleration of greedy matching algorithms,
based on approximate and fast methods of finding
nearest neighbors (e.g.
\cite{jiang2016, wang2018, wang2011, harwood2016,malkov2020}).
Another direction that helps to improve feature
matching problem is using alternative local
descriptors \citep{orb2011,piifd2010, brief2010}
for given keypoints. Naturally, the question of
how to chose keypoints arises, which is addressed,
for instance, in \citep{bai2020d3feat, Tian_2020_ACCV}.
For more complete overview of the field we refer
to \citep{ma2021} and references therein.

Finally, permutation estimation and related problems
have been recently investigated in different contexts
such as statistical seriation
\citep{flammarion2019optimal}, noisy sorting
\citep{mao2018minimax}, regression with shuffled data
\citep{pananjady2017linear,slawski2019linear},
isotonic regression and matrices
\citep{Mao18,pananjady2020isotonic,ma2020optimal},
crowd labeling \citep{ShahBW16a}, and recovery of
general discrete structure \citep{gao2019iterative}.

\section{Numerical results}\label{sec:numerical}
In this section, we report the results of some
numerical experiments carried out on simulated and real
data. We applied aforementioned methods LSNS and LSL
and computed different measures of their performance.
To get a more complete picture, we included in this
study the Least Sum of Squeres (LSS) estimator
and the greedy estimator. LSS is an unnormalized
version of LSNS, given by
\begin{align}\label{LSS}
    \hat\pi^{\textrm{LSS}}_{n,m}\in \text{arg}\min_{\pi}
    \sum_{i=1}^n \|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2.
\end{align}
It coincides with LSNS in the case of homoscedastaic
noise. The greedy estimator is obtained by sequentially
matching each vector from $\boldsymbol{X}$ to the nearest vector
from $\boldsymbol{X}^{\text{\tt\#}}$. Experiments were implemented using
python or matlab. For solving linear sum assignment
problems such as \eqref{LSL} or \eqref{LSS}, the
generalized and
improved versions of the Hungarian algorithm
were used \citep{Kuhn,Kuhn2010,Munkres,Duff}, implemented in SciPy library \citep{2020SciPy-NMeth}.

\paragraph{Experiment 1: Synthetic data with random features}
We first randomly generated $\pi^*$,
$\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ as follows. We
randomly sampled from uniform distribution on $[0, 2]$
independent variables $\tau_{ij}$, $i\in[m]$, $j\in[d]$.
Then, $({\theta}^\text{\tt\#}_i)_j$ are independently sampled from
the Gaussian distribution with 0 mean and variance
$\tau_{ij}$. Additionally, for every ${\theta}^\text{\tt\#}_i \in
\boldsymbol{\theta}^\texttt{\#}$ such that $i \notin J_{\pi^*}$
(${\theta}^\text{\tt\#}_i$ is an outlier), we incremented every
coordinate of ${\theta}^\text{\tt\#}_i$ by $i$. Entries of
$\boldsymbol{\sigma}^\text{\tt\#}$ were sampled from the uniform
distribution over $[0.5, 2]$. Sequences $\boldsymbol{X}$ and
$\boldsymbol{X}^{\text{\tt\#}}$ were generated according to \Cref{sec:model}
with $\pi^*(i) = i$ for $i\in[n]$. We applied to this
data the following matching algorithms: Greedy, LSS, LSNS and LSL.

We chose $n = 100$, $m = 130$ and $d = 50$, and
generated $N=50$ datasets according to the
foregoing process. For each dataset, we computed
the 0-1 error of the considered estimators and the
values of $(\bar\k_{\textup{in-in}},\bar\k_{\textup{in-out}})$.  We plotted in the
left (resp., the right) panel of \Cref{fig:exp1}
the error averaged over all datasets with a given
value of $\bar\k_{\textup{in-in}}$ (resp.\ $\bar\k_{\textup{in-out}}$). In this specific
example, we see that the error decreases fast with
$\bar\k_{\textup{in-in}}$, while there is no monotonicity of the error
rate as a function of $\bar\k_{\textup{in-out}}$.

\paragraph{Experiment 2: Synthetic data with deterministic features}
The second experiment is conducted on data generated
by features $\boldsymbol{\theta}^\texttt{\#}$ and variances $\boldsymbol{\sigma}^\text{\tt\#}$
inspired by the example constructed in the proof of
\Cref{thm:lower-bound}. More precisely, for some
real numbers $a$ and $b$ representing, respectively,
the scale of inlier-inlier distance $\bar\k_{\textup{in-in}}$ and
inlier-outlier distance $\bar\k_{\textup{in-out}}$, we set
${\theta}^\text{\tt\#}_k = [ka, 0,\ldots,0]^\top$ for
$k\in[n]$ and ${\theta}^\text{\tt\#}_{n+k} = [na + kb, 0,\ldots,
0]^\top$. We also used decreasing variances
$\s^{\text{\tt\#}}_k = 1/k^{3/2}$ for $k\in[m]$ and true identity
mapping $\pi^*(k)=k$ for $k\in [n]$. We chose
$n=100$, $m=120$ and dimension $d$ varying in
the set $\{10,20,40\}$. For each pair of values
$(a,b)$ in a suitably chosen grid, we repeated
$n_{\textup{rep}} = 400$ times the experiment
that consisted in generating data according to
\eqref{model} and computing estimators $\hat\pi_{n,m
}^{\textup{LSS}}$ and $\hat\pi_{n,m}^{\textup{LSL}}$
defined respectively by \eqref{LSS} and \eqref{LSL}.
We then computed, for each pair $(a,b)$ and for each
estimator LSS and LSL, the percentage of successful
detection among $n_{\textup{rep}}$ repetitions.

The obtained detection regions are depicted in
\Cref{fig:heatmap} in the form of heatmaps. This
visualisation allows us to grasp the forms of the
detection regions for the specific choice of
parameters considered in this example. The first
observation is that LSL is clearly superior to
LSS for all the considered values of the dimension.
Second, we clearly see the deterioration of the
detection region when the dimension $d$ becomes
larger. Third, the
values of $\bar\k_{\textup{in-out}}$ used in the plots are at least
one order of magnitude larger than those of
$\bar\k_{\textup{in-in}}$. This is in line with the claim of
\Cref{thm:mild}. We also observe in these pictures
that successful detection occurs when $\bar\k_{\textup{in-out}}$ is
larger than some threshold even if $\bar\k_{\textup{in-in}}$ is small.
This must be a nice feature of LSL and LSS in this
specific example, which unfortunately does not
generalize to other examples as shown by our
theoretical results.

\paragraph{Experiment 3:
Real data example}
This experiment is conducted on the IMC-PT 2020 dataset
from \citep{Jin2020} that consists of images of 16
different scenes with corresponding 3D point-clouds
of landmarks, which are used to obtain
(pseudo) ground-truth local keypoint matchings.
For a given scene, we sampled 1000 pairs of distinct
images of the same landmark with different camera
locations, angles, weather conditions etc. For each
image pair we generated 2D keypoints from original
set of 3D points (note that the same 3D point
appears in both of the images, so we have ground
truth keypoint matching between 2 images).
Subsequently, we computed SIFT descriptors
\citep{lowe2004distinctive} for every keypoint
in images using Python OpenCV interface \citep{opencv}.
Some pairs of images being more challenging than others,
we split the dataset into two sets of image pairs
in order to gain more understanding on the
behaviour of the algorithms. The challenging
pairs are those for which the OpenCV default
matching algorithm has accuracy less than $0.5$.
Then, for every image pair, we fixed randomly
chosen 100 keypoints in the first image (and
corresponding keypoints in the second image) and
added outliers to the second image from the
remaining keypoints. The outlier rate is chosen
to be between $0\%$ to $70\%$. Finally, 3
descriptor matching algorithms were applied
(OpenCV default matching algorithm, LSS and LSL).
{Note that $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^\text{\tt\#}$
from \eqref{model} are unknown, hence LSNS is
not applicable. One can consider using the
estimates $\hat\boldsymbol{\sigma}$ and $\boldsymbol{\hat\sigma}^\text{\tt\#}$
instead of $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^\text{\tt\#}$
in \eqref{LSNS}, but this is beyond the scope
of this paper.} Further details on this experiment
along with some additional results are deferred
to the supplemental material (\Cref{sec:numerical-details}).

The median estimation accuracy measured in the Hamming loss---for the image pairs from Temple Nara Japan scene---is plotted in \Cref{fig:hamming}. Similar results
for other scenes are presented in the supplemental material.
The error bars with borderlines corresponding to $75\%$ and $25\%$ percentiles are also displayed. The first observation is that LSS and LSL outperform the OpenCV matcher in terms of Hamming distance. Second, the task of feature matching becomes harder with the growth of outlier rate, with a deterioration that
seems to be linear in the rate of outliers. Notice that in this experiment the outliers can be very similar to the inliers and the separation condition imposed on $\bar\kappa_{\textup{in-out}}$, of \Cref{thm:upper-bound-LSL,thm:mild}, is violated. This implies that the larger the outlier rate the harder it is to find the correct match, causing a larger number of mistakes.

\begin{figure}
    \centering
    \includegraphics[scale=0.8, width=0.85\textwidth]{
    exp2_heatmaps}
    \vspace{-13pt}
    \caption{Heatmaps of the error rate
    of the LSL (top row) and the LSS (bottom row)
    estimators in Experiment 2. We chose $n=100$, $m=120$ and
    $d \in \{10,20,40\}$ from left to right. The
    parameter $a$ representing the scale of $\bar\k_{\textup{in-in}}$
    and corresponding to $y$-axis varies from 0.02 to
    0.08, whereas $b$ representing the scale of
    $\bar\k_{\textup{in-out}}$ and corresponding to $x$-axis varies
    from 0.3 to 10. Dark colour means that probability
    of successful detection is close to 1 (error rate
    close to zero).}
    \vspace{-10pt}
    \label{fig:heatmap}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.85\textwidth,
    height = 0.25\textwidth
    ]{accuracies_nara.pdf}
    \vspace{-10pt}

    \caption{The estimation accuracy measured in the Hamming loss of the estimated matching in Exp.~3 for different values of the outlier rate, $(m-n)/n$, varying from $0\%$ to $70\%$. The medians of estimation accuracy both for challenging pairs (right plot) and simple pairs (left plot) of images from Temple Nara  scene was computed using OpenCV, LSS and LSL matchers. The green region represents the interquartile range (lower and upper bounds being $25\%$ and $75\%$ percentiles, respectively).  }

    \vspace{-15pt}
    \label{fig:hamming}
\end{figure}

\section{Conclusion}\label{sec:conc}
We have investigated the detection regions in the
problem of estimation of the matching map between
two sequences of noisy vectors. We have shown that
the presence of outliers in one  of the two sequences
has a strong negative impact on the detection
region. Interestingly, this negative impact is
mitigated in the regime of mild heteroscedasticity,
\textit{i.e.}, when noise variances are of the same
order of magnitude. In the extremely favorable case
of homoscedastic noise (all variances are equal),
the presence of outliers does not make the problem
any harder, provided that the outliers are at least
as different from inliers as two distinct inliers are
different one another. Precise forms of the detection
region in these different cases can be found in
\Cref{tab:1}. The results of the numerical
experiments conducted on both synthetic and real data
confirm our findings and, furthermore, show the good
behaviour of the LSL estimator in terms of its
robustness to noise and to outliers, not only in
the problem of detection but also in the problem of
estimation. In the future, we plan to investigate the
case when both sequences contain outliers and to
obtain theoretical guarantees on the estimation
error measured by the Hamming distance.

\newpage

\part{}

\maketitle
\vspace{-5pt}
\begin{abstract}

    We consider the problem of finding the matching map between two sets of $d$-dimensional vectors from noisy observations, where the second set contains outliers.

    The matching map is then an injection, which can be
    consistently estimated only if the vectors of
    the second set are well separated.
    The main result shows that,

    in the high-dimensional setting,
    a detection region of unknown injection may be characterized by the sets
    of vectors for which the inlier-inlier distance is
    of order at least $d^{1/4}$ and the inlier-outlier
    distance is of order at least $d^{1/2}$.

    These rates are achieved using the estimated matching minimizing the sum of logarithms of distances between matched pairs of points.

    We also prove lower bounds establishing optimality of these rates.
    Finally, we report the results of numerical
    experiments on both synthetic and real world data that illustrate our theoretical results
    and provide further insight into the properties
    of the estimators studied in this work.
\end{abstract}

\vspace{-5pt}
\section{Introduction}\label{sec:intro}

Finding the best match between two clouds of
points is a problem encountered in many real
problems. In computer vision, one can look for
correspondences between two sets of local
descriptors extracted from two images. In
text analysis, one can be interested in matching
vector representations of the words of two
similar texts, potentially in two different
languages. The goal of the present work is to
gain theoretical understanding of the statistical
limits of the matching problem.

In the sequel, we use the notation $[n] = \{1,\ldots,
n\}$ for any integer $n$, and define $\|\cdot\|$
as the Euclidean norm in $\mathbb{R}^d$. Assume
that two independent sequences $\boldsymbol{X} = (X_i;
i\in[n])$ and $\boldsymbol{Y} =(Y_i; i\in[n])$ of independent
vectors are generated such that $X_i$ and $Y_i$
are drawn from the same distribution $P_i$
on $\mathbb{R}^d$, for every $i\in[n]$. The
statistician observes the sequence $\boldsymbol{X}$
and a shuffled version
$\boldsymbol{X}^{\text{\tt\#}}$ of the sequence $\boldsymbol{Y}$. More precisely,
$\boldsymbol{X}^{\text{\tt\#}}$ is such that $X^{\text{\tt\#}}_i = Y_{\pi^*(i)}$
for some unobserved permutation $\pi^*$. The
goal of matching is to estimate the permutation
$\pi^*$ from data $(\boldsymbol{X},\boldsymbol{X}^{\text{\tt\#}})$. In the case
of Gaussian distributions $P_i$, this problem
has been studied in \citep{jmlr_CD13, collier2016minimax}. Clearly, consistent estimation
of the matching map $\pi^*$ is impossible if
there are two data generating distributions
$P_i$ and $P_j$ that are very close. In
\citep{jmlr_CD13, collier2016minimax}, a
precise quantification of the separation
between these distributions is given that
enables consistent estimation of $\pi^*$.
Furthermore, it is shown that the permutation
estimator minimizing the sum of logarithms
of pairwise distances between the elements of
$\boldsymbol{X}$ and the elements of the shuffled version $\boldsymbol{X}^{\text{\tt\#}}$ is an optimal estimator of $\pi^*$.

In this paper, we extend the model studied in \citep{collier2016minimax} to the case when the
set $\boldsymbol{X}^{\text{\tt\#}}$ is contaminated by outliers. The number of outliers is supposed to be known and is equal to $m - n$, where $n = |\boldsymbol{X}|$ and $m = |\boldsymbol{X}^{\text{\tt\#}}|$ are the sizes of considered
two sequences, however the indices of the
outliers are unknown. All the distributions
are assumed in this paper to be
spherical Gaussian, although all the probabilistic
tools used in the proofs have their sub-Gaussian
counterparts. Thus, we consider that two spherical
Gaussian distributions
\footnote{We use the notation $\mathbf I_d$ for the $d\times d$ identity matrix}
$P_1 = \mathcal N_d(\mu_1, \sigma_1^2\mathbf I_d)$ and $P_2 = \mathcal
N_d(\mu_2, \sigma_2^2\mathbf I_d)$ are well
separated if the ``distance to noise ratio''
$\kappa(P_1,P_2) = \|\mu_1-\mu_2\|/\sqrt{
\sigma_1^2+\sigma_2^2}$ is large. Main findings
of \citep{collier2016minimax}, in terms of
smallest separation distance $\bar\kappa =
\min_{i\neq j}\kappa(P_i,P_j)$ are summarized
in the second columns of \Cref{tab:1}.
Likewise, the last column of the table
provides a summary of the contributions of
the present paper in terms of $\bar\k_{\textup{in-in}} =
\min_{i\neq j} \kappa(P_i,P_j)$ and $\bar\k_{\textup{in-out}} =
\min_{i,j} \kappa(P_i,Q_j)$, where $Q_1,\ldots,
Q_{m-n}$ are the distributions of the outliers.

An unexpected finding of this work is that the
``degree'' of heteroscedasticity of the model
has a strong impact on the separation distances
and the detection regions (sets of values of
$(\bar\k_{\textup{in-in}},\bar\k_{\textup{in-out}})$ for which it is possible to
detect the feature map $\pi^*$). This is in sharp
contrast with the outlier-free case, where consistent
estimation requires $\bar\kappa$ to be
at least of order $(d\log n)^{1/4}$ irrespective
from the behaviour of variances of $P_i$.
We prove in this work that in the high dimensional
regime $d\ge c\log n$, which is arguably more
appealing than the low dimensional regime $d\le c\log n$,
the following statements are true:
\begin{itemize}
    \item  If there is no heteroscedasticity, \textit{i.e.}, when all the variances are equal, consistent estimation
    of $\pi^*$ is possible if and only if $\bar\kappa
    = \bar\k_{\textup{in-in}}\wedge\bar\k_{\textup{in-out}} $
    is at least of order $(d\log(nm))^{1/4}$. This is
    the same rate as in the outlier-free case.
    \item If the heteroscedasticity is mild, \textit{i.e.},
    all the variances are of the same order, the
    condition
    $\bar\k_{\textup{in-in}}\gtrsim (d\log(nm))^{1/4}$ is the same as in the
    previous item, but the stronger condition $\bar\k_{\textup{in-out}}
    \gtrsim d^{1/2}$ is needed for the inlier-outlier
    separation distance.
    \item Finally, in the general heteroscedastic setting
    both $\bar\k_{\textup{in-in}}$ and $\bar\k_{\textup{in-out}}$ should be at least of
    order $d^{1/2}$. Furthermore, in all these cases
    consistent estimation is performed by the same
    estimator: the Least Sum of Logarithms (LSL).
\end{itemize}
Note also that the empirical evaluation done in this
paper shows that LSL is interesting not only from the
theoretical but also from the practical point of view.

\begin{table}
    \centering
    \begin{tabular}{ c|c|c }
        \toprule
         & Without outliers  & With outliers \\
         &\citep{collier2016minimax}& this paper \\
        \toprule
        \shortstack{known $\boldsymbol{\sigma}^\text{\tt\#}$ or \\
        all equal $\boldsymbol{\sigma}^\text{\tt\#}$s} & \shortstack{
        LSNS is optimal \\ $\bar\kappa \gtrsim (d\log n)^{1/4}$}
        & \shortstack{ LSNS is optimal [Thm. 1] \\
        $\bar\k_{\textup{in-in}} \wedge \bar\k_{\textup{in-out}}\gtrsim (d
        \log(nm))^{1/4}$} \\
        \midrule
        \begin{tabular}{c}
        \shortstack{unknown $\boldsymbol{\sigma}^\text{\tt\#}$ \\
        $\nicefrac{\sigma_{\max}}{\sigma_{\min}} \le C$} \\
        \midrule
        \shortstack{unknown $\boldsymbol{\sigma}^\text{\tt\#}$ \\
        arbitrary}
        \end{tabular}
        & \shortstack{LSL is optimal \\[5pt] $\bar\kappa
        \gtrsim (d\log n)^{1/4}$} &
        \begin{tabular}{c}
            \shortstack{LSL is optimal [Thm. 4] \\
            $\bar\k_{\textup{in-in}} \gtrsim (d\log(nm))^{1/4}$  \&
            $\bar\k_{\textup{in-out}}\gtrsim d^{1/2}$}  \\
            \midrule
            \shortstack{ LSL is optimal
            [Thm. 2, 3] \\ $\bar\k_{\textup{in-in}} \wedge
            \bar\k_{\textup{in-out}} \gtrsim d^{1/2}$}
        \end{tabular}
        \\
        \bottomrule
    \end{tabular}
    \vspace{-0.1cm}
    \caption{A brief overview of the contributions
    in the high-dimensional regime $d\ge c\log n$.
    The table provides the condition on the normalized
    inlier-inlier distance $\bar\k_{\textup{in-in}}$ and inlier-outlier
    distance $\bar\k_{\textup{in-out}}$, making it possible to consistently
    detect the matching map between two sets of
    $d$-dimensional vectors. LSL and LSNS refer to
    least sum of logarithms and least sum of normalized
    squares, respectively.
    }
    \vspace{-19pt}
    \label{tab:1}
\end{table}

\paragraph{Agenda} \Cref{sec:model} describes the
framework of the vector matching problem and
introduces the ter\-minology
used throughout this paper. Precise statements
of the main theoretical results are gathered in \Cref{sec:main}. The prior work is briefly discussed in
\Cref{sec:prior}. \Cref{sec:numerical} contains numerical experiments carried out both for synthetic and real data. A brief summary and some concluding remarks are presented in \Cref{sec:conc}. Proofs of all theoretical claims are deferred to the supplemental material.

\section{Problem Formulation}\label{sec:model}

We begin with formalizing the problem of matching two
sequences of feature vectors $(X_1,\ldots,X_n)$ and
$(X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m)$ with different sizes $n$
and $m$ such that $m \ge n \ge 2$. In what follows, we
assume that the observed feature vectors are randomly
generated from the model
\begin{equation}\label{model}
\begin{cases}
X_i = \theta_i + \sigma_i\xi_i \ , \\
X^{\text{\tt\#}}_j = {\theta}^\text{\tt\#}_j + \s^{\text{\tt\#}}_j\xi^\text{\tt\#}_j,
\end{cases}\quad i=1,\ldots,n \text{ and }
j=1,\ldots,m.
\end{equation}
In this model, illustrated in \Cref{fig:ani}, it
is assumed that
\vspace{-5pt}
\begin{itemize}\setlength{\itemsep}{-2pt}
\item $\boldsymbol{\theta}=(\theta_1,\ldots,\theta_n)$ and $\boldsymbol{\theta}^\text{\tt\#} =
({\theta}^\text{\tt\#}_1,\ldots,{\theta}^\text{\tt\#}_m)$ are two sequences of
vectors from $\mathbb R^d$, corresponding to the original
features, which are unavailable,
\item $\boldsymbol{\sigma}=(\sigma_1,\ldots,\sigma_n)^\top, \boldsymbol{\sigma}^\text{\tt\#}
= (\s^{\text{\tt\#}}_1,\ldots,\s^{\text{\tt\#}}_m)^\top$ are positive real
numbers corresponding to the magnitudes of the noise
contaminating each feature,
\item $\xi_1, \ldots, \xi_n$ and $\xi^\text{\tt\#}_1, \ldots,
\xi^\text{\tt\#}_m$ are two independent sequences of i.i.d.\
random vectors drawn from the Gaussian distribution
with zero mean and identity covariance matrix.
\end{itemize}

The simplest special case of \eqref{model},
considered in \citep{collier2016minimax}, corresponds
to the situation where a perfect matching exists
between the two sequences $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^\text{\tt\#}$.
This means that $m = n$ and, for some bijective
mapping $\pi^*:[n]\to[n]$, $\theta_i = {\theta}^\text{\tt\#}_{
\pi(i)}$ for all $i\in[n]$. In the general
case, both $\boldsymbol{X} = (X_1, \dots, X_n)$ and
$\boldsymbol{X}^{\text{\tt\#}} = (X^{\text{\tt\#}}_1, \dots, X^{\text{\tt\#}}_m)$ may
contain outliers, \textit{i.e.} feature vectors
that have no corresponding pair. In such a
situation, it is merely assumed that there exists
a set $S\subset [n]$ and an  injective mapping
$\pi^*:S\to [m]$ such that
\begin{align}\label{pistar}
    \theta_i = {\theta}^\text{\tt\#}_{\pi^*(i)} \quad\text
    {and}\quad \sigma_i = \s^{\text{\tt\#}}_{\pi^*(i)},
    \qquad \forall\,
    i\in S.
\end{align}
In this case we say that the vectors $\{X_i:i\in[n]\setminus S\}$ and $\{X^{\text{\tt\#}}_j:j\in[m]\setminus\pi^*(S)\}$ are
outliers. The ultimate goal is to estimate the
feature matching map $\pi^*$.

In this work we consider the case when $S = [n]$ and $m>n$. This means
that only the larger set of feature vectors, namely
$\boldsymbol{X}^{\text{\tt\#}}$, contains outliers. Let us also define the
set $O_{\pi^*} \triangleq [m] \setminus
\text{Im}(\pi^*)$,  which contains the indices of
outliers and satisfies  $|O_{\pi^*}| = m - n$.
Naturally, the feature vectors contained in $\boldsymbol{X}$,
as well as those vectors from $\boldsymbol{X}^{\text{\tt\#}}$ that are
not outliers, are called {\it{inliers}}.

In this formulation, the data generating distribution is defined by the parameters
$\boldsymbol{\theta}^\texttt{\#}$,
$\boldsymbol{\sigma}^\text{\tt\#}$ and $\pi^*$. We omit the set of parameters $\boldsymbol{\theta}$ and $\boldsymbol{\sigma}$, since they are automatically identified using $\pi^*$, $\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ by the formula $(\theta_i, \sigma_i) = ({\theta}^\text{\tt\#}_{\pi^*(i)},\s^{\text{\tt\#}}_{\pi^*(i)})$ for $i\in[n]$. Since our goal is to match the
feature vectors, we focus our attention on the problem of
estimating the parameter $\pi^*$ only, considering
$\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ as nuisance parameters.
In what follows, we denote by $\mathbf P_{ \boldsymbol{\theta}^\texttt{\#},
\boldsymbol{\sigma}^\text{\tt\#}, \pi^*}$ the probability distribution of the
sequence $(X_1,\ldots,X_n,X^{\text{\tt\#}}_1,\ldots,X^{\text{\tt\#}}_m)$
defined by (\ref{model}) under condition (\ref{pistar})
with $S= [n]$.

We are interested in designing estimators that have an expected error smaller than a prescribed level $\alpha$ under the weakest possible conditions on the nuisance parameter $\boldsymbol{\theta}^\texttt{\#}$ and noise level $\boldsymbol{\sigma}^\text{\tt\#}$. Clearly, the problem of matching becomes more difficult with hardly distinguishable features. To quantify this phenomenon, we introduce the normalized separation distance $\bar\k_{\textup{in-in}} = \bar\k_{\textup{in-in}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$ and the normalized outlier separation distance $\bar\k_{\textup{in-out}} = \bar\k_{\textup{in-out}}(\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*)$, which measure the minimal distance-to-noise ratio between inliers and the minimal distance-to-noise ratio between inliers and outliers, respectively. The precise
definitions read as
\begin{equation}
\bar\k_{\textup{in-in}} \triangleq \min_{\substack{ i,j\not\in O_{\pi^*}, \\ j \not=i}} \frac{\|{\theta}^\text{\tt\#}_i-{\theta}^\text{\tt\#}_j\|}{(\s^{\text{\tt\#}}_i{}^2+\s^{\text{\tt\#}}_j{}^2)^{1/2}},\qquad
\bar\k_{\textup{in-out}} \triangleq \min_{\substack{ i\not\in O_{\pi^*}, \\ j \in O_{\pi^*}}} \frac{\|{\theta}^\text{\tt\#}_i-{\theta}^\text{\tt\#}_j\|}{(\s^{\text{\tt\#}}_i{}^2+\s^{\text{\tt\#}}_j{}^2)^{1/2}}. \label{kappa}
\end{equation}
Notice that $\bar\k_{\textup{in-in}}$ can be rewritten as
\begin{align}
    \bar\k_{\textup{in-in}} = \min_{\substack{i, j \in [n] \\ i \neq j}} \frac{\| \theta_i - \theta_j\|}{(\sigma_i^2 + \sigma_j^2)^{1/2}}.
\end{align}

\begin{figure}
\centering
\includegraphics[height=0.32\textwidth]{
ani_match_crop2.pdf}
\hspace{20pt}
\vspace{-11pt}
\includegraphics[height=0.32\textwidth]{
ani_match_crop1.pdf}
\caption{Illustration of the considered framework described in \eqref{model}.
We wish to match a set of 7 patches extracted from
the first image to the 9 patches from the second
image. The picture on the left shows the locations
of patches as well as the true matching map $\pi^*$
(the yellow lines).
}
\vspace{-13pt}
\label{fig:ani}
\end{figure}

Clearly, if $\bar\k_{\textup{in-in}}=0$ or, $\bar\k_{\textup{in-out}}=0$, there are
two identical feature vectors in $\boldsymbol{X}^{\text{\tt\#}}$. In such
a situation, assuming $\sigma_i$'s are all equal,
the parameter $\pi^*$ is nonidentifiable, in the
sense that there exist two different permutations
$\pi^*_1$ and $\pi^*_2$ such that the distributions
$\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_1}$ and
$\mathbf P_{\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*_2}$
coincide.  Therefore, to ensure the existence of
consistent estimators of $\pi^*$ it is necessary
to impose the conditions $\bar\k_{\textup{in-in}}>0$ and
$\bar\k_{\textup{in-out}}>0$. Moreover, good estimators are those
consistently estimating $\pi^*$ even if either
$\bar\k_{\textup{in-in}}$ or $\bar\k_{\textup{in-out}}$ are small. We are interested
here in finding the detection boundary in terms
of the order of magnitude of $(\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}})
$. More precisely, for any given $\alpha\in(0,1)$
we wish to find a region $\mathcal R_{n,m,d}^\alpha$
in $\mathbb R^2$ such that:
\vspace{-8pt}
\setlength{\leftmargini}{2em}
\begin{itemize}
    \item There is an estimator $\hat\pi_{n,m}$ of
    $\pi^*$ satisfying

        $\mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi
        \neq \pi^*) \le \alpha$

    for every $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$
    lying in the detection region, \textit{i.e.}, for
    which $(\bar\k_{\textup{in-in}}, \bar\kappa_{
    \textup{in-out}}) \in \mathcal R_{n,m,d}^\alpha$.
    \item There is a constant $C<1$ such that for any
    estimator $\bar\pi_{n,m}$ of $\pi^*$, we can find
    a parameter value $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$
    in the region $\{(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*) :
    (\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}})
    \in C\mathcal R_{n,m,d}^\alpha\}$
    such that $\bar\pi$ fails to detect $\pi^*$ with a
    probability larger than $\alpha$.
\end{itemize}
\vspace{-4pt}
Let us make two remarks. First, note that in the
outlier-free case considered
in \citep{collier2016minimax}, $\bar\k_{\textup{in-out}}$
is meaningless and, therefore, the detection region  is one-dimensional $\mathcal R_{n,m,d}^\alpha$. Thus, it
is necessarily a half-line and is proven to be of the
form $\bar\k_{\textup{in-in}}\ge C(\log \nicefrac{
n}{\alpha})^{1/2} \vee (d\log \nicefrac{n}{\alpha} )^{1/4}$
for some universal constant $C$. Second, the aforementioned definition of the detection region $\mathcal R_{n,m,d}^\alpha$ does not guarantee
its uniqueness (even up to a scaling by a universal
constant). This is in contrast with the
outlier-free case. To overcome this difficulty, we
look for $\mathcal R_{n,m,d}^\alpha$ of
the form $[t_{\textup{in-in}},+\infty)\times [t_{\textup{in-out}},+\infty)$ with the smallest possible
threshold $t_{\textup{in-out}}$ for the normalized inlier-outlier distance $\bar\k_{\textup{in-out}}$.

\def\boldsymbol{\Xi}{\boldsymbol{\Xi}}
\section{Main theoretical results}\label{sec:main}

In this section, we have collected the main theoretical findings of the paper.  When the noise is homoscedastic,
\textit{i.e.},  when all $\boldsymbol{\sigma}$'s are equal,
the results obtained by \cite{collier2016minimax} in the outlier-free setting can be easily extended to the setting
with outliers. Therefore, in the present paper, we focus
on the heteroscedastic case. For the sake of clarity of
exposition, we will present the results in the case of known
variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$ prior to investigating
the more interesting case of unknown variances.

The detection regions we study below are based on the
profile maximum likelihood estimator. The model presented in \eqref{model} has the parameter $\boldsymbol{\Xi} = (\boldsymbol{\theta}^\text{\tt\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi)$, while the observations are the sequences
of feature vectors $\boldsymbol{X}$ and $\boldsymbol{X}^{\text{\tt\#}}$. The negative
log-likelihood of this model is given by
\begin{align}
    \ell_n(\boldsymbol{\Xi}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\}) &= \sum_{i=1}^n \bigg(\frac{\| X_i - {\theta}^\text{\tt\#}_{\pi(i)}\|_2^2}{2\s^{\text{\tt\#}2}_{\pi(i)}} + \frac 1 2 \log(\s^{\text{\tt\#}2}_{\pi(i)})\bigg) +  \sum_{j=1}^m \bigg(\frac{\| X^{\text{\tt\#}}_j - {\theta}^\text{\tt\#}_j\|_2^2}{2\s^{\text{\tt\#}2}_j} + \frac 1 2 \log(\s^{\text{\tt\#}2}_j)\bigg).
\end{align}
The profile negative log-likelihood is then defined as
the minimum
with respect to $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#})$ of
the log-likelihood $\ell_n(\boldsymbol{\Xi}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\})$.

\subsection{Warming up: known variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

One can check
that the minimization with respect to $\boldsymbol{\theta}^\texttt{\#}$ leads
to the variance-dependent cost function
\begin{align}\label{liklihood-sigma}
    \ell_n(\pi,\boldsymbol{\sigma}^\text{\tt\#}; \{\boldsymbol{X}, \boldsymbol{X}^{\text{\tt\#}}\}) =  \sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2} + \sum_{i=1}^n \frac{1}{2}\log(\s^{\text{\tt\#}2}_{\pi(i)}) + \sum_{j=1}^m \frac{1}{2} \log(\s^{\text{\tt\#}2}_j).
\end{align}
When $m=n$ and there is no outlier, the last two sums of the
last display do not depend on $\pi$ and, therefore, the
maximum profile likelihood estimator of $\pi^*$ is obtained
by the Least Sum of Normalized Squares (LSNS) criterion
\begin{align}\label{LSNS}
    \hat\pi^{\textrm{LSNS}}_{n,m}\in \text{arg}\min_{\pi}
    \sum_{i=1}^n \frac{\|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2}{\sigma_{i}^2+{\s^{\text{\tt\#}}_{\pi(i)}}^2},
\end{align}
where the minimum is over all injective mappings
$\pi:[n]\to[m]$. This, and the other estimators defined
in this work, can be efficiently computed using suitable
versions of the Hungarian algorithm \citep{Kuhn,Kuhn2010,
Munkres}. As shows the next theorem, it
turns out that even when $m>n$, the estimator
$\hat\pi^{\textrm{LSNS}}_{ n,m } $  defined above
leads to an optimal detection region.

\begin{theorem}[Upper bound for LSNS]\label{thm:upperLSNS}
Let $\alpha \in (0,1)$ and condition \eqref{pistar}
be fulfilled. If the separation distances $\bar\k_{\textup{in-in}}$
and $\bar\k_{\textup{in-out}}$ corresponding to $(\boldsymbol{\theta}^\texttt{\#},
\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$ and defined in \eqref{kappa}
satisfy the condition
\begin{align}\label{thresh:1}
    \min\{\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}}\} \ge 4 \Big\{ \big(
    d\log(\nicefrac{4nm}{\alpha})\big)^{1/4} \vee
    \big( 2\log(\nicefrac{8nm}{ \alpha})\big)^{1/2}
    \Big\}
\end{align}
then the LSNS estimator defined in \eqref{LSNS} detects
the true matching map $\pi^*$ with probability at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi^{
    \textup{LSNS}}_{n,m} = \pi^*) \ge 1-\alpha.
\end{align}
\vspace{-15pt}
\end{theorem}

The similarity---both its statement and its proof---
of this result to its counterpart
in the outlier-free setting might suggest that
the presence of outliers does not make the problem
any harder from a statistical point of view.
However, this is not true in the more
appealing setting of unknown variances.

\subsection{Detection of $\pi^*$ for unknown and
arbitrary variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

The LSNS procedure analyzed in \Cref{thm:upperLSNS}
exploits the values of known noise variances to normalize
the squares of distances between vectors $X_i$ and
$X^{\text{\tt\#}}_{\pi(i)}$. Therefore, LSNS is inapplicable
in the case of unknown noise variances.
Instead, we consider the Least Sum of
Logarithms (LSL) estimator
\begin{align}\label{LSL}
    \hat\pi^{\textup{LSL}}_{n,m} \triangleq \arg\min_{\pi:[n]\to [m]}\sum_{i=1}^n {\log \|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2},
\end{align}
where the minimum is over all injective maps $\pi
: [n] \to [m]$. This estimator can be seen as the
minimizer of a criterion defined as the minimum of
the cost function from \eqref{liklihood-sigma} with
respect to $\boldsymbol{\sigma}^\text{\tt\#}$ under the constraint
$\min_{j\not\in\Im(\pi)} \s^{\text{\tt\#}}_{j} \ge
\sigma_{\min}$, for some fixed (but unknown)
constant  $\sigma_{\min}>0$.

To provide a quick overview of what follows, let us
stick in the remaining of this paragraph to the
case $\log(nm) = O(d)$ so that the right hand side
of \eqref{thresh:1} is of order $\big( d\log(nm)
\big){}^{1/4}$. Recall that in the outlier-free
case, the LSL estimator has been shown to perform
as well as the LSNS while having the advantage of
not requiring the knowledge of variances
$\boldsymbol{\sigma}^\text{\tt\#}$ \citep{collier2016minimax}.
Somewhat unexpectedly, the situation is significantly
different in the presence of outliers.  Indeed, the
best we manage to prove in the presence of
outliers is that the detection of the matching
map by LSL is possible whenever $\min\{ \bar\kappa_{
\textup{in-in}}, \bar\k_{\textup{in-out}}\}\ge
C\sqrt{d}$ for some sufficiently large constant
$C$. The precise statement being given in the next
theorem, let us mention right away that the
discrepancy between this rate $\sqrt{d}$ and the
rate $\big( d\log(nm) \big){}^{1/4}$ in
\eqref{thresh:1} is due to the inherent hardness
of the setting and not merely an artefact of the
proof. This will be made clear below.

\begin{theorem}[Upper bound for LSL]
\label{thm:upper-bound-LSL}
Let $\alpha \in (0,1/2)$ and condition
\eqref{pistar} be fulfilled. If the separation
distances $\bar\k_{\textup{in-in}}$ and $\bar\k_{\textup{in-out}}$ corresponding
to $(\boldsymbol{\theta}^\texttt{\#},\boldsymbol{\sigma}^\text{\tt\#},\pi^*)$ and
defined by \eqref{kappa} satisfy
\begin{align}\label{detect}
    \min\{\bar\k_{\textup{in-in}}, \bar\k_{\textup{in-out}}\} \ge \sqrt{2d} + 4
    \Big\{\big(2d\log (\nicefrac{4nm}{ \alpha})
    \big)^{1/4} \vee \big(3 \log(\nicefrac{8nm}{
    \alpha}\big)^{1/2}\Big\}
\end{align}
then the LSL estimator \eqref{LSL} detects
the matching map $\pi^*$ with probability at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}
    (\hat\pi^{\textup{LSL}}_{n,m} = \pi^*)
    \ge 1 - \alpha.
\end{align}
\vspace{-15pt}
\end{theorem}

This result is disappointing since it requires
the distance between different feature vectors to
be larger than $\sqrt{2d}$ in order to be able
to consistently estimate the matching map $\pi^*$.
As we show below, without any further condition
(for instance, on the noise variances), this
rate cannot be improved. Furthermore, the
rate $\sqrt{d}$ is optimal not only for LSL but
also for the larger class of so called
\textit{distance based $M$-estimators.}

We will say that an estimator $\hat\pi_n$ of $\pi^*$
is a distance based $M$-estimator, if for a sequence of non-decreasing functions $\rho_i:\mathbb{R}_+
\to \mathbb{R}$, $i=1,\ldots,n$, the following is correct
\begin{align}
    \hat\pi_n \in \arg\min_{\pi :[n] \to [m]} \sum_{i=1}^n \rho_i\big(\| X_i - X^{\text{\tt\#}}_{\pi(i)}\|\big),
\end{align}
where the minimum is over all injective mappings
$\pi:[n]\to[m]$. We denote by $\mathcal M$ the set
of all distance based $M$-estimators.
We show that there is indeed a setup where $\bar\k_{\textup{in-in}} \wedge \bar\kappa_{ \textup{in-out}}$ is as large as $0.2\sqrt{d}$ but
any estimator from $\mathcal M$ fails to detect
$\pi^*$ with probability at least $1/4$.  The next theorem formalizes the described result.
\begin{theorem}[Lower bound]\label{thm:lower-bound}
Assume that $m > n \ge 4$ and $d \ge 422\log(4n)$. There exists a triplet $(\boldsymbol{\sigma}^\text{\tt\#}, \boldsymbol{\theta}^\text{\tt\#}, \pi^*)$ such that  $\bar\k_{\textup{in-in}} = \bar\k_{\textup{in-out}} = \sqrt{d/20}$ and
    \begin{align}\label{lower:1}
        \inf_{\hat\pi \in \mathcal{M}}

        \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}(\hat\pi
        \neq \pi^*) > \nicefrac 1 4.
    \end{align}
    \vspace{-15pt}
\end{theorem}

The proof of this theorem, postponed to the appendix,
is constructive. This means that we exhibit a triplet
$(\boldsymbol{\sigma}^\text{\tt\#}, \boldsymbol{\theta}^\text{\tt\#}, \pi^*)$ satisfying
\eqref{lower:1}. Careful inspection shows that in the
case $d = O(\log(nm))$ the same triplet satisfies
$\bar\k_{\textup{in-in}} \wedge \bar\kappa_{ \textup{in-out}} \asymp
\sqrt{\log(nm)}$ and \eqref{lower:1} is still true.
This implies that the order of magnitude of the right
hand side of \eqref{detect} is optimal both in the
high-dimensional regime $d\ge 422\log(4n) $ and in the
low-dimensional regime $d< 422\log(4n)$. This shows the optimality of LSL among all estimators
from $\mathcal{M}$. Note that the estimator $\hat\pi^{
\textup{LSNS}}_{n,m}$ does not belong to the family of
\textit{distance based M-estimators}. Furthermore, in the
low dimensional regime $d=O(\log(nm))$, the separation
rate of the LSL, $\sqrt{\log(nm)}$, is the same as that of the LSNS.

In the next section we show that under some mild
conditions on $\boldsymbol{\sigma}^\text{\tt\#}$ it is indeed possible to
obtain different rates for $\bar\k_{\textup{in-in}}$
and $\bar\kappa_{ \textup{in-out}}$, namely we show that
if $\bar\k_{\textup{in-in}} \gtrsim d^{1/4}$ and
$\bar\k_{\textup{in-out}} \gtrsim d^{1/2}$ then the
LSL estimator detects correct matching with high
probability.

\subsection{Detection of $\pi^*$ for unknown and mildly
varying variances $\boldsymbol{\sigma},\boldsymbol{\sigma}^\text{\tt\#}$}

The results of the last two theorems are disappointing,
since they indicate that the features should be very
different from one another for detection of the matching
map to be possible. An interesting finding, presented
below, is that strong constraint can be significantly
alleviated in the context of mild heteroscedasticity.
By mild heteroscedasticity we understand here the
situation in which all variances $\s^{\text{\tt\#}}_i$ are of
the same order of magnitude.

\begin{theorem}[Upper bound under mild
heteroscedasticity]\label{thm:mild}
Let $r_\sigma = \max_{i,j}(\s^{\text{\tt\#}}_i/\s^{\text{\tt\#}}_j)
<\infty$. If the separation distances $\bar\k_{\textup{in-in}}$
and $\bar\k_{\textup{in-out}}$ defined in \eqref{kappa} satisfy
\begin{align}
    \bar\k_{\textup{in-in}} &\ge
    2\big(4d\log(\nicefrac{4nm}{\alpha})
    \big)^{1/4} + 2\big(2\log(\nicefrac{4nm}{
    \alpha})\big)^{1/2}\\
    \bar\k_{\textup{in-out}}
        &\ge\  \sqrt{2(r_\sigma-1)d} + 2\big(4
        r_\sigma^2 d\log(\nicefrac{4nm}{\alpha}
        ) \big)^{1/4} + 2\big(2 r_\sigma \log(
        \nicefrac{4nm}{\alpha}) \big)^{1/2} ,
\end{align}
then the LSL estimator \eqref{LSL}
detects the matching map $\pi^*$ with probability
at least
$1-\alpha$, that is
\begin{align}
    \mathbf P_{\boldsymbol{\theta}^\texttt{\#}, \boldsymbol{\sigma}^\text{\tt\#}, \pi^*}
    (\hat\pi^{\textup{LSL}}_{n,m} = \pi^*) \ge
    1 - \alpha.
\end{align}
\vspace{-15pt}
\end{theorem}
Note that a lower bound similar to that of
\Cref{thm:lower-bound} can be proved in the case
of mild heterescodestacity as well, showing that
there is an example for which $\bar\k_{\textup{in-in}}$ is of order
$d^{1/4}$, $\bar\k_{\textup{in-out}}$ is of order $d^{1/2}$ and
any estimator from $\mathcal M$ fails to detect
$\pi^*$ with probability at least $1/4$.

We complete this section by summarizing the
joint contribution of
\Cref{thm:upperLSNS,thm:upper-bound-LSL,thm:lower-bound,thm:mild}.
To simplify this discussion, we consider two
cases: high-dimensional case refers to $d\ge \log(4nm/
\alpha)$ (presented in \Cref{tab:1}) and low-dimensional case refers to
the condition $d< \log(4nm/\alpha)$. In the
high dimensional setting with arbitrary noise
variances, the detection region for the LSL
estimator is given by $\{\bar\k_{\textup{in-in}}\wedge \bar\k_{\textup{in-out}}
\ge 15 \sqrt{d}\}$, which is much worse than
the detection region for LSNS, $\{\bar\k_{\textup{in-in}}
\wedge \bar\k_{\textup{in-out}}\ge 8 (d\log(4nm/\alpha))^{1/4}
\}$, obtained in the known-variance scenario.
Somewhat surprisingly, in such a setting,
even a strong assumption on the outliers, such
as requiring them to be at least at a
distance $0.2\sqrt{d}$ of the inliers, is not
enough for relaxing the assumption on the
inlier-inlier separation distance. Finally, on a positive
note, in the intermediate case of mildly
varying variances, the detection region for
the LSL estimator is of the form
$\{\bar\k_{\textup{in-in}}\ge 7 (d\log(4nm/\alpha))^{1/4}; \
\bar\k_{\textup{in-out}} \ge 10 \sqrt{d}\}$. This means that
if the outliers are at a distance $\Omega(
\sqrt{d})$ of the inliers, then the LSL
recovers the true matching under the same
condition on $\bar\k_{\textup{in-in}}$ as in the outlier-free
setting.
}

\begin{figure
\begin{center}
\includegraphics[width=0.85\textwidth]{figure1c.png}
\vspace{-10pt}
\caption{The performance of the methods (Greedy, LSS,
LSL, LSNS) in the setup described in Exp. 1.
Plots show the error rate (percentage of repetitions
in which the estimated matching map differs from the
true one) as a function of separation distances. The
left plot illustrates that LSS, LSL and LSNS require
much lower value of $\bar\k_{\textup{in-in}}$ in order to find the
correct mapping. The left plot shows that there is
a clear improvement of the error when the inlier-inlier
separation distance increases, while the right plot shows
that the error rate might not be made small by
augmenting $\bar\k_{\textup{in-out}}$ only. }
\label{fig:exp1}
\end{center}
\vspace{-20pt}
\end{figure}

\section{Other related work}\label{sec:prior}

Measuring the quality of the various statistical
procedures of decision making by their minimal
separation rates became the standard in hypotheses
testing, see the seminal papers \citep{Burnashev,
Ingster82} and the monographs \citep{Ingster,JudNem}.
Currently this approach is widely adopted in machine
learning literature \citep{Xing20,Wolfer,Balchard,
RamdasISW16,Wei,Collier12a}. Beyond the classical
setting of two hypotheses, it can also be applied
to multiple hypotheses testing framework, for instance,
variable selection \citep{Ndaoud,azais2020multiple,
ComDal1} or the matching problem considered here.

On the other hand, feature matching is a well
studied problem in computer vision. In recent
years, a great deal of attention was devoted to
the acceleration of greedy matching algorithms,
based on approximate and fast methods of finding
nearest neighbors (e.g.
\cite{jiang2016, wang2018, wang2011, harwood2016,malkov2020}).
Another direction that helps to improve feature
matching problem is using alternative local
descriptors \citep{orb2011,piifd2010, brief2010}
for given keypoints. Naturally, the question of
how to chose keypoints arises, which is addressed,
for instance, in \citep{bai2020d3feat, Tian_2020_ACCV}.
For more complete overview of the field we refer
to \citep{ma2021} and references therein.

Finally, permutation estimation and related problems
have been recently investigated in different contexts
such as statistical seriation
\citep{flammarion2019optimal}, noisy sorting
\citep{mao2018minimax}, regression with shuffled data
\citep{pananjady2017linear,slawski2019linear},
isotonic regression and matrices
\citep{Mao18,pananjady2020isotonic,ma2020optimal},
crowd labeling \citep{ShahBW16a}, and recovery of
general discrete structure \citep{gao2019iterative}.

\section{Numerical results}\label{sec:numerical}
In this section, we report the results of some
numerical experiments carried out on simulated and real
data. We applied aforementioned methods LSNS and LSL
and computed different measures of their performance.
To get a more complete picture, we included in this
study the Least Sum of Squeres (LSS) estimator
and the greedy estimator. LSS is an unnormalized
version of LSNS, given by
\begin{align}\label{LSS}
    \hat\pi^{\textrm{LSS}}_{n,m}\in \text{arg}\min_{\pi}
    \sum_{i=1}^n \|X_{i}-X^{\text{\tt\#}}_{\pi(i)}\|^2.
\end{align}
It coincides with LSNS in the case of homoscedastaic
noise. The greedy estimator is obtained by sequentially
matching each vector from $\boldsymbol{X}$ to the nearest vector
from $\boldsymbol{X}^{\text{\tt\#}}$. Experiments were implemented using
python or matlab. For solving linear sum assignment
problems such as \eqref{LSL} or \eqref{LSS}, the
generalized and
improved versions of the Hungarian algorithm
were used \citep{Kuhn,Kuhn2010,Munkres,Duff}, implemented in SciPy library \citep{2020SciPy-NMeth}.

\paragraph{Experiment 1: Synthetic data with random features}
We first randomly generated $\pi^*$,
$\boldsymbol{\theta}^\texttt{\#}$ and $\boldsymbol{\sigma}^\text{\tt\#}$ as follows. We
randomly sampled from uniform distribution on $[0, 2]$
independent variables $\tau_{ij}$, $i\in[m]$, $j\in[d]$.
Then, $({\theta}^\text{\tt\#}_i)_j$ are independently sampled from
the Gaussian distribution with 0 mean and variance
$\tau_{ij}$. Additionally, for every ${\theta}^\text{\tt\#}_i \in
\boldsymbol{\theta}^\texttt{\#}$ such that $i \notin J_{\pi^*}$
(${\theta}^\text{\tt\#}_i$ is an outlier), we incremented every
coordinate of ${\theta}^\text{\tt\#}_i$ by $i$. Entries of
$\boldsymbol{\sigma}^\text{\tt\#}$ were sampled from the uniform
distribution over $[0.5, 2]$. Sequences $\boldsymbol{X}$ and
$\boldsymbol{X}^{\text{\tt\#}}$ were generated according to \Cref{sec:model}
with $\pi^*(i) = i$ for $i\in[n]$. We applied to this
data the following matching algorithms: Greedy, LSS, LSNS and LSL.

We chose $n = 100$, $m = 130$ and $d = 50$, and
generated $N=50$ datasets according to the
foregoing process. For each dataset, we computed
the 0-1 error of the considered estimators and the
values of $(\bar\k_{\textup{in-in}},\bar\k_{\textup{in-out}})$.  We plotted in the
left (resp., the right) panel of \Cref{fig:exp1}
the error averaged over all datasets with a given
value of $\bar\k_{\textup{in-in}}$ (resp.\ $\bar\k_{\textup{in-out}}$). In this specific
example, we see that the error decreases fast with
$\bar\k_{\textup{in-in}}$, while there is no monotonicity of the error
rate as a function of $\bar\k_{\textup{in-out}}$.

\paragraph{Experiment 2: Synthetic data with deterministic features}
The second experiment is conducted on data generated
by features $\boldsymbol{\theta}^\texttt{\#}$ and variances $\boldsymbol{\sigma}^\text{\tt\#}$
inspired by the example constructed in the proof of
\Cref{thm:lower-bound}. More precisely, for some
real numbers $a$ and $b$ representing, respectively,
the scale of inlier-inlier distance $\bar\k_{\textup{in-in}}$ and
inlier-outlier distance $\bar\k_{\textup{in-out}}$, we set
${\theta}^\text{\tt\#}_k = [ka, 0,\ldots,0]^\top$ for
$k\in[n]$ and ${\theta}^\text{\tt\#}_{n+k} = [na + kb, 0,\ldots,
0]^\top$. We also used decreasing variances
$\s^{\text{\tt\#}}_k = 1/k^{3/2}$ for $k\in[m]$ and true identity
mapping $\pi^*(k)=k$ for $k\in [n]$. We chose
$n=100$, $m=120$ and dimension $d$ varying in
the set $\{10,20,40\}$. For each pair of values
$(a,b)$ in a suitably chosen grid, we repeated
$n_{\textup{rep}} = 400$ times the experiment
that consisted in generating data according to
\eqref{model} and computing estimators $\hat\pi_{n,m
}^{\textup{LSS}}$ and $\hat\pi_{n,m}^{\textup{LSL}}$
defined respectively by \eqref{LSS} and \eqref{LSL}.
We then computed, for each pair $(a,b)$ and for each
estimator LSS and LSL, the percentage of successful
detection among $n_{\textup{rep}}$ repetitions.

The obtained detection regions are depicted in
\Cref{fig:heatmap} in the form of heatmaps. This
visualisation allows us to grasp the forms of the
detection regions for the specific choice of
parameters considered in this example. The first
observation is that LSL is clearly superior to
LSS for all the considered values of the dimension.
Second, we clearly see the deterioration of the
detection region when the dimension $d$ becomes
larger. Third, the
values of $\bar\k_{\textup{in-out}}$ used in the plots are at least
one order of magnitude larger than those of
$\bar\k_{\textup{in-in}}$. This is in line with the claim of
\Cref{thm:mild}. We also observe in these pictures
that successful detection occurs when $\bar\k_{\textup{in-out}}$ is
larger than some threshold even if $\bar\k_{\textup{in-in}}$ is small.
This must be a nice feature of LSL and LSS in this
specific example, which unfortunately does not
generalize to other examples as shown by our
theoretical results.

\paragraph{Experiment 3:
Real data example}
This experiment is conducted on the IMC-PT 2020 dataset
from \citep{Jin2020} that consists of images of 16
different scenes with corresponding 3D point-clouds
of landmarks, which are used to obtain
(pseudo) ground-truth local keypoint matchings.
For a given scene, we sampled 1000 pairs of distinct
images of the same landmark with different camera
locations, angles, weather conditions etc. For each
image pair we generated 2D keypoints from original
set of 3D points (note that the same 3D point
appears in both of the images, so we have ground
truth keypoint matching between 2 images).
Subsequently, we computed SIFT descriptors
\citep{lowe2004distinctive} for every keypoint
in images using Python OpenCV interface \citep{opencv}.
Some pairs of images being more challenging than others,
we split the dataset into two sets of image pairs
in order to gain more understanding on the
behaviour of the algorithms. The challenging
pairs are those for which the OpenCV default
matching algorithm has accuracy less than $0.5$.
Then, for every image pair, we fixed randomly
chosen 100 keypoints in the first image (and
corresponding keypoints in the second image) and
added outliers to the second image from the
remaining keypoints. The outlier rate is chosen
to be between $0\%$ to $70\%$. Finally, 3
descriptor matching algorithms were applied
(OpenCV default matching algorithm, LSS and LSL).
{Note that $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^\text{\tt\#}$
from \eqref{model} are unknown, hence LSNS is
not applicable. One can consider using the
estimates $\hat\boldsymbol{\sigma}$ and $\boldsymbol{\hat\sigma}^\text{\tt\#}$
instead of $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}^\text{\tt\#}$
in \eqref{LSNS}, but this is beyond the scope
of this paper.} Further details on this experiment
along with some additional results are deferred
to the supplemental material (\Cref{sec:numerical-details}).

The median estimation accuracy measured in the Hamming loss---for the image pairs from Temple Nara Japan scene---is plotted in \Cref{fig:hamming}. Similar results
for other scenes are presented in the supplemental material.
The error bars with borderlines corresponding to $75\%$ and $25\%$ percentiles are also displayed. The first observation is that LSS and LSL outperform the OpenCV matcher in terms of Hamming distance. Second, the task of feature matching becomes harder with the growth of outlier rate, with a deterioration that
seems to be linear in the rate of outliers. Notice that in this experiment the outliers can be very similar to the inliers and the separation condition imposed on $\bar\kappa_{\textup{in-out}}$, of \Cref{thm:upper-bound-LSL,thm:mild}, is violated. This implies that the larger the outlier rate the harder it is to find the correct match, causing a larger number of mistakes.

\begin{figure}
    \centering
    \includegraphics[scale=0.8, width=0.85\textwidth]{
    exp2_heatmaps}
    \vspace{-13pt}
    \caption{Heatmaps of the error rate
    of the LSL (top row) and the LSS (bottom row)
    estimators in Experiment 2. We chose $n=100$, $m=120$ and
    $d \in \{10,20,40\}$ from left to right. The
    parameter $a$ representing the scale of $\bar\k_{\textup{in-in}}$
    and corresponding to $y$-axis varies from 0.02 to
    0.08, whereas $b$ representing the scale of
    $\bar\k_{\textup{in-out}}$ and corresponding to $x$-axis varies
    from 0.3 to 10. Dark colour means that probability
    of successful detection is close to 1 (error rate
    close to zero).}
    \vspace{-10pt}
    \label{fig:heatmap}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.85\textwidth,
    height = 0.25\textwidth
    ]{accuracies_nara.pdf}
    \vspace{-10pt}

    \caption{The estimation accuracy measured in the Hamming loss of the estimated matching in Exp.~3 for different values of the outlier rate, $(m-n)/n$, varying from $0\%$ to $70\%$. The medians of estimation accuracy both for challenging pairs (right plot) and simple pairs (left plot) of images from Temple Nara  scene was computed using OpenCV, LSS and LSL matchers. The green region represents the interquartile range (lower and upper bounds being $25\%$ and $75\%$ percentiles, respectively).  }

    \vspace{-15pt}
    \label{fig:hamming}
\end{figure}

\section{Conclusion}\label{sec:conc}
We have investigated the detection regions in the
problem of estimation of the matching map between
two sequences of noisy vectors. We have shown that
the presence of outliers in one  of the two sequences
has a strong negative impact on the detection
region. Interestingly, this negative impact is
mitigated in the regime of mild heteroscedasticity,
\textit{i.e.}, when noise variances are of the same
order of magnitude. In the extremely favorable case
of homoscedastic noise (all variances are equal),
the presence of outliers does not make the problem
any harder, provided that the outliers are at least
as different from inliers as two distinct inliers are
different one another. Precise forms of the detection
region in these different cases can be found in
\Cref{tab:1}. The results of the numerical
experiments conducted on both synthetic and real data
confirm our findings and, furthermore, show the good
behaviour of the LSL estimator in terms of its
robustness to noise and to outliers, not only in
the problem of detection but also in the problem of
estimation. In the future, we plan to investigate the
case when both sequences contain outliers and to
obtain theoretical guarantees on the estimation
error measured by the Hamming distance.

\newpage