Optimal detection of the feature matching map in presence of noise and outliers

Abstract:
We consider the problem of finding the matching map between two sets of d2 dimensional vectors from noisy observations, where the second set contains outliers. 3 The matching map is then an injection, which can be consistently estimated only 4 if the vectors of the second set are well separated. The main result shows that, 5 in the high-dimensional setting, a detection region of unknown injection may be 6 characterized by the sets of vectors for which the inlier-inlier distance is of order 7 at least d and the inlier-outlier distance is of order at least d. These rates 8 are achieved using the estimated matching minimizing the sum of logarithms of 9 distances between matched pairs of points. We also prove lower bounds establishing 10 optimality of these rates. Finally, we report the results of numerical experiments on 11 both synthetic and real world data that illustrate our theoretical results and provide 12 further insight into the properties of the estimators studied in this work. 13

text: We consider the problem of finding the matching map between two sets of d-2 dimensional vectors from noisy observations, where the second set contains outliers.3 The matching map is then an injection, which can be consistently estimated only4 if the vectors of the second set are well separated. The main result shows that,5 in the high-dimensional setting, a detection region of unknown injection may be6 characterized by the sets of vectors for which the inlier-inlier distance is of order7 at least d1/4 and the inlier-outlier distance is of order at least d1/2. These rates8 are achieved using the estimated matching minimizing the sum of logarithms of9 distances between matched pairs of points. We also prove lower bounds establishing10 optimality of these rates. Finally, we report the results of numerical experiments on11 both synthetic and real world data that illustrate our theoretical results and provide12 further insight into the properties of the estimators studied in this work.13
1 Introduction14
Finding the best match between two clouds of points is a problem encountered in many real problems.15 In computer vision, one can look for correspondences between two sets of local descriptors extracted16 from two images. In text analysis, one can be interested in matching vector representations of the17 words of two similar texts, potentially in two different languages. The goal of the present work is to18 gain theoretical understanding of the statistical limits of the matching problem.19
In the sequel, we use the notation [n] = {1, . . . , n} for any integer n, and define ‖ · ‖ as the Euclidean20 norm in Rd. Assume that two independent sequences X = (Xi; i ∈ [n]) and Y = (Yi; i ∈ [n]) of21 independent vectors are generated such that Xi and Yi are drawn from the same distribution Pi on Rd,22 for every i ∈ [n]. The statistician observes the sequenceX and a shuffled versionX# of the sequence23 Y . More precisely,X# is such that X#i = Yπ∗(i) for some unobserved permutation π
∗. The goal of24 matching is to estimate the permutation π∗ from data (X,X#). In the case of Gaussian distributions25 Pi, this problem has been studied in (Collier and Dalalyan, 2013, 2016). Clearly, consistent estimation26 of the matching map π∗ is impossible if there are two data generating distributions Pi and Pj that are27 very close. In (Collier and Dalalyan, 2013, 2016), a precise quantification of the separation between28 these distributions is given that enables consistent estimation of π∗. Furthermore, it is shown that the29 permutation estimator minimizing the sum of logarithms of pairwise distances between the elements30 ofX and the elements of the shuffled versionX# is an optimal estimator of π∗.31
In this paper, we extend the model studied in (Collier and Dalalyan, 2016) to the case when the setX#32 is contaminated by outliers. The number of outliers is supposed to be known and is equal to m− n,33 where n = |X| and m = |X#| are the sizes of considered two sequences, however the indices of the34
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
outliers are unknown. All the distributions are assumed in this paper to be spherical Gaussian, although35 all the probabilistic tools used in the proofs have their sub-Gaussian counterparts. Thus, we consider36 that two spherical Gaussian distributions 1 P1 = Nd(µ1, σ21Id) and P2 = Nd(µ2, σ22Id) are well37 separated if the “distance to noise ratio” κ(P1, P2) = ‖µ1 − µ2‖/ √ σ21 + σ 2 2 is large. Main findings38 of (Collier and Dalalyan, 2016), in terms of smallest separation distance κ̄ = mini 6=j κ(Pi, Pj) are39 summarized in the second columns of Table 1. Likewise, the last column of the table provides40 a summary of the contributions of the present paper in terms of κ̄in-in = mini6=j κ(Pi, Pj) and41 κ̄in-out = mini,j κ(Pi, Qj), where Q1, . . . , Qm−n are the distributions of the outliers.42
An unexpected finding of this work is that the “degree” of heteroscedasticity of the model has a strong43 impact on the separation distances and the detection regions (sets of values of (κ̄in-in, κ̄in-out) for which44 it is possible to detect the feature map π∗). This is in sharp contrast with the outlier-free case, where45 consistent estimation requires κ̄ to be at least of order (d log n)1/4 irrespective from the behaviour46 of variances of Pi. We prove in this work that in the high dimensional regime d ≥ c log n, which is47 arguably more appealing than the low dimensional regime d ≤ c log n, the following statements are48 true:49
• If there is no heteroscedasticity, i.e., when all the variances are equal, consistent estimation50 of π∗ is possible if and only if κ̄ = κ̄in-in ∧ κ̄in-out is at least of order (d log(nm))1/4. This is51 the same rate as in the outlier-free case.52
• If the heteroscedasticity is mild, i.e., all the variances are of the same order, the condition53 κ̄in-in & (d log(nm))1/4 is the same as in the previous item, but the stronger condition54 κ̄in-out & d1/2 is needed for the inlier-outlier separation distance.55
• Finally, in the general heteroscedastic setting both κ̄in-in and κ̄in-out should be at least of56 order d1/2. Furthermore, in all these cases consistent estimation is performed by the same57 estimator: the Least Sum of Logarithms (LSL).58
Note also that the empirical evaluation done in this paper shows that LSL is interesting not only from59 the theoretical but also from the practical point of view.60
Agenda Section 2 describes the framework of the vector matching problem and introduces the ter-61 minology used throughout this paper. Precise statements of the main theoretical results are gathered in62 Section 3. The prior work is briefly discussed in Section 4. Section 5 contains numerical experiments63 carried out both for synthetic and real data. A brief summary and some concluding remarks are64 presented in Section 6. Proofs of all theoretical claims are deferred to the supplemental material.65
2 Problem Formulation66
We begin with formalizing the problem of matching two sequences of feature vectors (X1, . . . , Xn)67 and (X#1 , . . . , X # m) with different sizes n and m such that m ≥ n ≥ 2. In what follows, we assume68
1We use the notation Id for the d× d identity matrix
that the observed feature vectors are randomly generated from the model69 { Xi = θi + σiξi ,
X#j = θ # j + σ # j ξ # j ,
i = 1, . . . , n and j = 1, . . . ,m. (1)
In this model, illustrated in Figure 1, it is assumed that70
• θ = (θ1, . . . , θn) and θ# = (θ#1 , . . . , θ # m) are two sequences of vectors from Rd, corre-71 sponding to the original features, which are unavailable,72 • σ = (σ1, . . . , σn)>,σ# = (σ#1 , . . . , σ # m) > are positive real numbers corresponding to the73 magnitudes of the noise contaminating each feature,74 • ξ1, . . . , ξn and ξ#1 , . . . , ξ # m are two independent sequences of i.i.d. random vectors drawn75 from the Gaussian distribution with zero mean and identity covariance matrix.76
The simplest special case of (1), considered in (Collier and Dalalyan, 2016), corresponds to the77 situation where a perfect matching exists between the two sequences θ and θ#. This means that78 m = n and, for some bijective mapping π∗ : [n]→ [n], θi = θ#π(i) for all i ∈ [n]. In the general case,79 both X = (X1, . . . , Xn) and X# = (X#1 , . . . , X # m) may contain outliers, i.e. feature vectors that80 have no corresponding pair. In such a situation, it is merely assumed that there exists a set S ⊂ [n]81 and an injective mapping π∗ : S → [m] such that82
θi = θ # π∗(i) and σi = σ # π∗(i), ∀ i ∈ S. (2)
In this case we say that the vectors {Xi : i ∈ [n] \ S} and {X#j : j ∈ [m] \ π∗(S)} are outliers. The83 ultimate goal is to estimate the feature matching map π∗.84
In this work we consider the case when S = [n] and m > n. This means that only the larger set of85 feature vectors, namelyX#, contains outliers. Let us also define the set Oπ∗ , [m] \ Im(π∗), which86 contains the indices of outliers and satisfies |Oπ∗ | = m− n. Naturally, the feature vectors contained87 inX , as well as those vectors fromX# that are not outliers, are called inliers.88
In this formulation, the data generating distribution is defined by the parameters θ#, σ# and π∗. We89 omit the set of parameters θ and σ, since they are automatically identified using π∗, θ# and σ# by90 the formula (θi, σi) = (θ#π∗(i), σ # π∗(i)) for i ∈ [n]. Since our goal is to match the feature vectors, we91 focus our attention on the problem of estimating the parameter π∗ only, considering θ# and σ# as92 nuisance parameters. In what follows, we denote by Pθ#,σ#,π∗ the probability distribution of the93 sequence (X1, . . . , Xn, X#1 , . . . , X # m) defined by (1) under condition (2) with S = [n].94
We are interested in designing estimators that have an expected error smaller than a prescribed level95 α under the weakest possible conditions on the nuisance parameter θ# and noise level σ#. Clearly,96 the problem of matching becomes more difficult with hardly distinguishable features. To quantify97 this phenomenon, we introduce the normalized separation distance κ̄in-in = κ̄in-in(θ#,σ#, π∗) and98 the normalized outlier separation distance κ̄in-out = κ̄in-out(θ#,σ#, π∗), which measure the minimal99 distance-to-noise ratio between inliers and the minimal distance-to-noise ratio between inliers and100 outliers, respectively. The precise definitions read as101
κ̄in-in , min i,j 6∈Oπ∗ , j 6=i
‖θ#i − θ#j ‖ (σ#i 2 + σ#j 2)1/2
, κ̄in-out , min i 6∈Oπ∗ , j∈Oπ∗
‖θ#i − θ#j ‖ (σ#i 2 + σ#j 2)1/2 . (3)
Notice that κ̄in-in can be rewritten as102
κ̄in-in = min i,j∈[n] i 6=j
‖θi − θj‖ (σ2i + σ 2 j ) 1/2 .
Clearly, if κ̄in-in = 0 or, κ̄in-out = 0, there are two identical feature vectors in X#. In such a103 situation, assuming σi’s are all equal, the parameter π∗ is nonidentifiable, in the sense that there104 exist two different permutations π∗1 and π ∗ 2 such that the distributions Pθ#,σ#,π∗1 and Pθ#,σ#,π∗2105 coincide. Therefore, to ensure the existence of consistent estimators of π∗ it is necessary to impose106 the conditions κ̄in-in > 0 and κ̄in-out > 0. Moreover, good estimators are those consistently estimating107 π∗ even if either κ̄in-in or κ̄in-out are small. We are interested here in finding the detection boundary in108 terms of the order of magnitude of (κ̄in-in, κ̄in-out). More precisely, for any given α ∈ (0, 1) we wish109 to find a regionRαn,m,d in R2 such that:110
• There is an estimator π̂n,m of π∗ satisfying Pθ#,σ#,π∗(π̂ 6= π∗) ≤ α for every (θ#,σ#, π∗)111 lying in the detection region, i.e., for which (κ̄in-in, κ̄in-out) ∈ Rαn,m,d.112 • There is a constant C < 1 such that for any estimator π̄n,m of π∗, we can find a parameter value113 (θ#,σ#, π∗) in the region {(θ#,σ#, π∗) : (κ̄in-in, κ̄in-out) ∈ CRαn,m,d} such that π̄ fails to detect114 π∗ with a probability larger than α.115
Let us make two remarks. First, note that in the outlier-free case considered in (Collier and Dalalyan,116 2016), κ̄in-out is meaningless and, therefore, the detection region is one-dimensionalRαn,m,d. Thus,117 it is necessarily a half-line and is proven to be of the form κ̄in-in ≥ C(log n/α)1/2 ∨ (d log n/α)1/4118 for some universal constant C. Second, the aforementioned definition of the detection region119 Rαn,m,d does not guarantee its uniqueness (even up to a scaling by a universal constant). This is120 in contrast with the outlier-free case. To overcome this difficulty, we look for Rαn,m,d of the form121 [tin-in,+∞)× [tin-out,+∞) with the smallest possible threshold tin-out for the normalized inlier-outlier122 distance κ̄in-out.123
3 Main theoretical results124
In this section, we have collected the main theoretical findings of the paper. When the noise is125 homoscedastic, i.e., when all σ’s are equal, the results obtained by Collier and Dalalyan (2016) in126 the outlier-free setting can be easily extended to the setting with outliers. Therefore, in the present127 paper, we focus on the heteroscedastic case. For the sake of clarity of exposition, we will present128 the results in the case of known variances σ,σ# prior to investigating the more interesting case of129 unknown variances.130
The detection regions we study below are based on the profile maximum likelihood estimator. The131 model presented in (1) has the parameter Ξ = (θ#,σ#, π), while the observations are the sequences132 of feature vectorsX andX#. The negative log-likelihood of this model is given by133
`n(Ξ; {X,X#}) = n∑ i=1 (‖Xi − θ#π(i)‖22 2σ#2π(i) + 1 2 log(σ#2π(i)) ) + m∑ j=1 (‖X#j − θ#j ‖22 2σ#2j + 1 2 log(σ#2j ) ) .
The profile negative log-likelihood is then defined as the minimum with respect to (θ#,σ#) of the134 log-likelihood `n(Ξ; {X,X#}).135

heading: 3.1 Warming up: known variances σ,σ#136
text: One can check that the minimization with respect to θ# leads to the variance-dependent cost function137
`n(π,σ #; {X,X#}) = n∑ i=1 ‖Xi −X#π(i)‖ 2 σ2i + σ # π(i) 2 + n∑ i=1 1 2 log(σ#2π(i)) + m∑ j=1 1 2 log(σ#2j ). (4)
When m = n and there is no outlier, the last two sums of the last display do not depend on π138 and, therefore, the maximum profile likelihood estimator of π∗ is obtained by the Least Sum of139
Normalized Squares (LSNS) criterion140
π̂LSNSn,m ∈ arg min π n∑ i=1 ‖Xi −X#π(i)‖ 2 σ2i + σ # π(i) 2 , (5)
where the minimum is over all injective mappings π : [n] → [m]. This, and the other estimators141 defined in this work, can be efficiently computed using suitable versions of the Hungarian algorithm142 (Kuhn, 1955, 2010; Munkres, 1957). As shows the next theorem, it turns out that even when m > n,143 the estimator π̂LSNSn,m defined above leads to an optimal detection region.144
Theorem 1 (Upper bound for LSNS) Let α ∈ (0, 1) and condition (2) be fulfilled. If the separation145 distances κ̄in-in and κ̄in-out corresponding to (θ #,σ#, π∗) and defined in (3) satisfy the condition146
min{κ̄in-in, κ̄in-out} ≥ 4 {( d log(4nm/α) )1/4 ∨ (2 log(8nm/α))1/2} (6) then the LSNS estimator defined in (5) detects the true matching map π∗ with probability at least147 1− α, that is148
Pθ#,σ#,π∗(π̂ LSNS n,m = π ∗) ≥ 1− α.
149 The similarity—both its statement and its proof— of this result is to its counterpart in the outlier-free150 setting might suggest that the presence of outliers does not make the problem any harder from a151 statistical point of view. However, this is not true in the more appealing setting of unknown variances.152

heading: 3.2 Detection of π∗ for unknown and arbitrary variances σ,σ#153
text: The LSNS procedure analyzed in Theorem 1 exploits the values of known noise variances to normalize154 the squares of distances between vectors Xi and X#π(i). Therefore, LSNS is inapplicable in the case155 of unknown noise variances. Instead, we consider the Least Sum of Logarithms (LSL) estimator156
π̂LSLn,m , arg min π:[n]→[m] n∑ i=1 log ‖Xi −X#π(i)‖ 2, (7)
where the minimum is over all injective maps π : [n] → [m]. This estimator can be seen as the157 minimizer of a criterion defined as the minimum of the cost function from (4) with respect to σ#158 under the constraint minj 6∈Im(π) σ#j ≥ σmin, for some fixed (but unknown) constant σmin > 0.159 To provide a quick overview of what follows, let us stick in the remaining of this paragraph to the160 case log(nm) = O(d) so that the right hand side of (6) is of order ( d log(nm) ) 1/4. Recall that in the161 outlier-free case, the LSL estimator has been shown to perform as well as the LSNS while having the162 advantage of not requiring the knowledge of variances σ# (Collier and Dalalyan, 2016). Somewhat163 unexpectedly, the situation is significantly different in the presence of outliers. Indeed, the best we164 manage to prove in the presence of outliers is that the detection of the matching map by LSL is165 possible whenever min{κ̄in-in, κ̄in-out} ≥ C √ d for some sufficiently large constant C. The precise166 statement being given in the next theorem, let us mention right away that the discrepancy between167 this rate √ d and the rate ( d log(nm) ) 1/4 in (6) is due to the inherent hardness of the setting and not168 merely an artefact of the proof. This will be made clear below.169
Theorem 2 (Upper bound for LSL) Let α ∈ (0, 1/2) and condition (2) be fulfilled. If the separa-170 tion distances κ̄in-in and κ̄in-out corresponding to (θ #,σ#, π∗) and defined by (3) satisfy171
min{κ̄in-in, κ̄in-out} ≥ √ 2d+ 4 {( 2d log(4nm/α) )1/4 ∨ (3 log(8nm/α)1/2} (8)
then the LSL estimator (7) detects the matching map π∗ with probability at least 1− α, that is172
Pθ#,σ#,π∗(π̂ LSL n,m = π ∗) ≥ 1− α.
173 This result is disappointing since it requires the distance between different feature vectors to be larger174 than √ 2d in order to be able to consistently estimate the matching map π∗. As we show below,175 without any further condition (for instance, on the noise variances), this rate cannot be improved.176
Furthermore, the rate √ d is optimal not only for LSL but also for the larger class of so called distance177 based M -estimators.178
We will say that an estimator π̂n of π∗ is a distance based M -estimator, if for a sequence of non-179 decreasing functions ρi : R+ → R, i = 1, . . . , n, the following is correct180
π̂n ∈ arg min π:[n]→[m] n∑ i=1 ρi ( ‖Xi −X#π(i)‖ ) ,
where the minimum is over all injective mappings π : [n] → [m]. We denote byM the set of all181 distance based M -estimators. We show that there is indeed a setup where κ̄in-in ∧ κ̄in-out is as large as182 0.2 √ d but any estimator fromM fails to detect π∗ with probability at least 1/4. The next theorem183 formalizes the described result.184
Theorem 3 (Lower bound) Assume that m > n ≥ 4 and d ≥ 422 log(4n). There exists a triplet185 (σ#,θ#, π∗) such that κ̄in-in = κ̄in-out = √ d/20 and186
inf π̂∈M
Pθ#,σ#,π∗(π̂ 6= π∗) > 1/4. (9)
187 The proof of this theorem, postponed to the appendix, is constructive. This means that we exhibit188 a triplet (σ#,θ#, π∗) satisfying (9). Careful inspection shows that in the case d = O(log(nm))189 the same triplet satisfies κ̄in-in ∧ κ̄in-out √ log(nm) and (9) is still true. This implies that the190 order of magnitude of the right hand side of (8) is optimal both in the high-dimensional regime191 d ≥ 422 log(4n) and in the low-dimensional regime d < 422 log(4n). This shows the optimality of192 LSL among all estimators fromM. Note that the estimator π̂LSNSn,m does not belong to the family of193 distance based M-estimators. Furthermore, in the low dimensional regime d = O(log(nm)), the194 separation rate of the LSL, √ log(nm), is the same as that of the LSNS.195
In the next section we show that under some mild conditions on σ# it is indeed possible to obtain196 different rates for κ̄in-in and κ̄in-out, namely we show that if κ̄in-in & d1/4 and κ̄in-out & d1/2 then the197 LSL estimator detects correct matching with high probability.198

heading: 3.3 Detection of π∗ for unknown and mildly varying variances σ,σ#199
text: The results of the last two theorems are disappointing, since they indicate that the features should200 be very different from one another for detection of the matching map to be possible. An interesting201 finding, presented below, is that strong constraint can be significantly alleviated in the context of mild202 heteroscedasticity. By mild heteroscedasticity we understand here the situation in which all variances203 σ#i are of the same order of magnitude.204
Theorem 4 (Upper bound under mild heteroscedasticity) Let rσ = maxi,j(σ#i /σ#j ) <∞. If the205 separation distances κ̄in-in and κ̄in-out defined in (3) satisfy206
κ̄in-in ≥ 2 ( 4d log(4nm/α) )1/4 + 2 ( 2 log(4nm/α) )1/2 κ̄in-out ≥ √ 2(rσ − 1)d+ 2 ( 4r2σd log(4nm/α) )1/4 + 2 ( 2rσ log(4nm/α) )1/2 ,
then the LSL estimator (7) detects the matching map π∗ with probability at least 1− α, that is207
Pθ#,σ#,π∗(π̂ LSL n,m = π ∗) ≥ 1− α.
208 Note that a lower bound similar to that of Theorem 3 can be proved in the case of mild heterescodestac-209 ity as well, showing that there is an example for which κ̄in-in is of order d1/4, κ̄in-out is of order d1/2210 and any estimator fromM fails to detect π∗ with probability at least 1/4.211 We complete this section by summarizing the joint contribution of Theorems 1 to 4. To simplify this212 discussion, we consider two cases: high-dimensional case refers to d ≥ log(4nm/α) (presented in213 Table 1) and low-dimensional case refers to the condition d < log(4nm/α). In the high dimensional214 setting with arbitrary noise variances, the detection region for the LSL estimator is given by {κ̄in-in ∧215 κ̄in-out ≥ 15 √ d}, which is much worse than the detection region for LSNS, {κ̄in-in ∧ κ̄in-out ≥216 8(d log(4nm/α))1/4}, obtained in the known-variance scenario. Somewhat surprisingly, in such a217
setting, even a strong assumption on the outliers, such as requiring them to be at least at a distance218 0.2 √ d of the inliers, is not enough for relaxing the assumption on the inlier-inlier separation distance.219 Finally, on a positive note, in the intermediate case of mildly varying variances, the detection region220 for the LSL estimator is of the form {κ̄in-in ≥ 7(d log(4nm/α))1/4; κ̄in-out ≥ 10 √ d}. This means221
that if the outliers are at a distance Ω( √ d) of the inliers, then the LSL recovers the true matching222 under the same condition on κ̄in-in as in the outlier-free setting.223
4 Other related work224
Measuring the quality of the various statistical procedures of decision making by their minimal225 separation rates became the standard in hypotheses testing, see the seminal papers (Burnashev, 1979;226 Ingster, 1982) and the monographs (Ingster and Suslina, 2003; Juditsky and Nemirovski, 2020).227 Currently this approach is widely adopted in machine learning literature (Xing et al., 2020; Wolfer228 and Kontorovich, 2020; Blanchard et al., 2018; Ramdas et al., 2016; Wei et al., 2019; Collier, 2012).229 Beyond the classical setting of two hypotheses, it can also be applied to multiple hypotheses testing230 framework, for instance, variable selection (Ndaoud and Tsybakov, 2020; Azaïs and de Castro, 2020;231 Comminges and Dalalyan, 2012) or the matching problem considered here.232
On the other hand, feature matching is a well studied problem in computer vision. In recent years,233 a great deal of attention was devoted to the acceleration of greedy matching algorithms, based on234 approximate and fast methods of finding nearest neighbors (e.g. Jiang et al. (2016); Wang et al.235 (2018); Wang (2011); Harwood and Drummond (2016); Malkov and Yashunin (2020)). Another236 direction that helps to improve feature matching problem is using alternative local descriptors (Rublee237 et al., 2011; Chen et al., 2010; Calonder et al., 2010) for given keypoints. Naturally, the question of238 how to chose keypoints arises, which is addressed, for instance, in (Bai et al., 2020; Tian et al., 2020).239 For more complete overview of the field we refer to (Ma et al., 2021) and references therein.240
Finally, permutation estimation and related problems have been recently investigated in different241 contexts such as statistical seriation (Flammarion et al., 2019), noisy sorting (Mao et al., 2018),242 regression with shuffled data (Pananjady et al., 2017; Slawski and Ben-David, 2019), isotonic243 regression and matrices (Mao et al., 2020; Pananjady and Samworth, 2020; Ma et al., 2020), crowd244 labeling (Shah et al., 2021), and recovery of general discrete structure (Gao and Zhang, 2019).245
5 Numerical results246
In this section, we report the results of some numerical experiments carried out on simulated and real247 data. We applied aforementioned methods LSNS and LSL and computed different measures of their248 performance. To get a more complete picture, we included in this study the Least Sum of Squeres249 (LSS) estimator and the greedy estimator. LSS is an unnormalized version of LSNS, given by250
π̂LSSn,m ∈ arg min π n∑ i=1 ‖Xi −X#π(i)‖ 2. (10)
It coincides with LSNS in the case of homoscedastaic noise. The greedy estimator is obtained251 by sequentially matching each vector from X to the nearest vector from X#. Experiments were252 implemented using python or matlab. For solving linear sum assignment problems such as (7) or (10),253 the generalized and improved versions of the Hungarian algorithm were used (Kuhn, 1955, 2010;254 Munkres, 1957; Duff and Koster, 2001), implemented in SciPy library (Virtanen et al., 2020).255
Experiment 1: Synthetic data with random features We first randomly generated π∗, θ# and256 σ# as follows. We randomly sampled from uniform distribution on [0, 2] independent variables τij ,257 i ∈ [m], j ∈ [d]. Then, (θ#i )j are independently sampled from the Gaussian distribution with 0 mean258 and variance τij . Additionally, for every θ#i ∈ θ
# such that i /∈ Jπ∗ (θ#i is an outlier), we incremented259 every coordinate of θ#i by i. Entries of σ
# were sampled from the uniform distribution over [0.5, 2].260 SequencesX andX# were generated according to Section 2 with π∗(i) = i for i ∈ [n]. We applied261 to this data the following matching algorithms: Greedy, LSS, LSNS and LSL.262
We chose n = 100, m = 130 and d = 50, and generated N = 50 datasets according to the foregoing263 process. For each dataset, we computed the 0-1 error of the considered estimators and the values264 of (κ̄in-in, κ̄in-out). We plotted in the left (resp., the right) panel of Figure 2 the error averaged over265 all datasets with a given value of κ̄in-in (resp. κ̄in-out). In this specific example, we see that the error266 decreases fast with κ̄in-in, while there is no monotonicity of the error rate as a function of κ̄in-out.267
Experiment 2: Synthetic data with deterministic features The second experiment is conducted268 on data generated by features θ# and variances σ# inspired by the example constructed in the proof269 of Theorem 3. More precisely, for some real numbers a and b representing, respectively, the scale of270 inlier-inlier distance κ̄in-in and inlier-outlier distance κ̄in-out, we set θ#k = [ka, 0, . . . , 0]
> for k ∈ [n]271 and θ#n+k = [na+ kb, 0, . . . , 0] >. We also used decreasing variances σ#k = 1/k 3/2 for k ∈ [m] and272 true identity mapping π∗(k) = k for k ∈ [n]. We chose n = 100, m = 120 and dimension d varying273 in the set {10, 20, 40}. For each pair of values (a, b) in a suitably chosen grid, we repeated nrep = 400274 times the experiment that consisted in generating data according to (1) and computing estimators275 π̂LSSn,m and π̂ LSL n,m defined respectively by (10) and (7). We then computed, for each pair (a, b) and for276 each estimator LSS and LSL, the percentage of successful detection among nrep repetitions.277
The obtained detection regions are depicted in Figure 3 in the form of heatmaps. This visualisation278 allows us to grasp the forms of the detection regions for the specific choice of parameters considered279 in this example. The first observation is that LSL is clearly superior to LSS for all the considered280 values of the dimension. Second, we clearly see the deterioration of the detection region when the281 dimension d becomes larger. Third, the values of κ̄in-out used in the plots are at least one order of282 magnitude larger than those of κ̄in-in. This is in line with the claim of Theorem 4. We also observe in283 these pictures that successful detection occurs when κ̄in-out is larger than some threshold even if κ̄in-in284 is small. This must be a nice feature of LSL and LSS in this specific example, which unfortunately285 does not generalize to other examples as shown by our theoretical results.286
Experiment 3: Real data example This experiment is conducted on the IMC-PT 2020 dataset287 from Jin et al. (2020) that consists of images of 16 different scenes with corresponding 3D point-288 clouds of landmarks, which are used to obtain (pseudo) ground-truth local keypoint matchings. For a289 given scene, we sampled 1000 pairs of distinct images of the same landmark with different camera290 locations, angles, weather conditions etc. For each image pair we generated 2D keypoints from291 original set of 3D points (note that the same 3D point appears in both of the images, so we have292 ground truth keypoint matching between 2 images). Subsequently, we computed SIFT descriptors293 (Lowe, 2004) for every keypoint in images using Python OpenCV interface (Itseez, 2015). Some294 pairs of images being more challenging than others, we split the dataset into two sets of image pairs in295 order to gain more understanding on the behaviour of the algorithms. The challenging pairs are those296 for which the OpenCV default matching algorithm has accuracy less than 0.5. Then, for every image297 pair, we fixed randomly chosen 100 keypoints in the first image (and corresponding keypoints in the298 second image) and added outliers to the second image from the remaining keypoints. The outlier299 rate is chosen to be between 0% to 70%. Finally, 3 descriptor matching algorithms were applied300 (OpenCV default matching algorithm, LSS and LSL). Note that σ and σ# from (1) are unknown,301 hence LSNS is not applicable. One can consider using the estimates σ̂ and σ̂# instead of σ and σ#302 in (5), but this is beyond the scope of this paper.. Further details on this experiment along with some303 additional results are deferred to the supplemental material (Appendix D).304
The median estimation accuracy measured in the Hamming loss—for the image pairs from Temple305 Nara Japan scene—is plotted in Figure 4. Similar results for other scenes are presented in the306 supplemental material. The error bars with borderlines corresponding to 75% and 25% percentiles307 are also displayed. The first observation is that LSS and LSL outperform the OpenCV matcher in308 terms of Hamming distance. Second, the task of feature matching becomes harder with the growth309 of outlier rate, with a deterioration that seems to be linear in the rate of outliers. Notice that in this310 experiment the outliers can be very similar to the inliers and the separation condition imposed on311 κ̄in-out, of Theorems 2 and 4, is violated. This implies that the larger the outlier rate the harder it is to312 find the correct match, causing a larger number of mistakes.313
6 Conclusion314
We have investigated the detection regions in the problem of estimation of the matching map between315 two sequences of noisy vectors. We have shown that the presence of outliers in one of the two316 sequences has a strong negative impact on the detection region. Interestingly, this negative impact is317 mitigated in the regime of mild heteroscedasticity, i.e., when noise variances are of the same order318 of magnitude. In the extremely favorable case of homoscedastic noise (all variances are equal), the319 presence of outliers does not make the problem any harder, provided that the outliers are at least as320 different from inliers as two distinct inliers are different one another. Precise forms of the detection321 region in these different cases can be found in Table 1. The results of the numerical experiments322 conducted on both synthetic and real data confirm our findings and, furthermore, show the good323 behaviour of the LSL estimator in terms of its robustness to noise and to outliers, not only in the324 problem of detection but also in the problem of estimation. In the future, we plan to investigate the325 case when both sequences contain outliers and to obtain theoretical guarantees on the estimation error326 measured by the Hamming distance.327
References328 Azaïs, J. M. and de Castro, Y. (2020). Multiple testing and variable selection along least angle329 regression’s path.330
Bai, X., Luo, Z., Zhou, L., Fu, H., Quan, L., and Tai, C.-L. (2020). D3feat: Joint learning of dense331 detection and description of 3d local features. In Proceedings of the IEEE/CVF Conference on332 Computer Vision and Pattern Recognition (CVPR).333
Blanchard, G., Carpentier, A., and Gutzeit, M. (2018). Minimax Euclidean separation rates for testing334 convex hypotheses in Rd. Electron. J. Stat., 12(2):3713–3735.335
Burnashev, M. V. (1979). On the minimax detection of an inaccurately known signal in a white336 Gaussian noise background. Theory Probab. Appl., 24:107–119.337
Calonder, M., Lepetit, V., Strecha, C., and Fua, P. (2010). Brief: Binary robust independent elementary338 features. In Daniilidis, K., Maragos, P., and Paragios, N., editors, Computer Vision – ECCV 2010,339 pages 778–792, Berlin, Heidelberg. Springer Berlin Heidelberg.340
Chen, J., Tian, J., Lee, N., Zheng, J., Smith, R. T., and Laine, A. F. (2010). A partial intensity invariant341 feature descriptor for multimodal retinal image registration. IEEE Transactions on Biomedical342 Engineering, 57(7):1707–1718.343
Collier, O. (2012). Minimax hypothesis testing for curve registration. In Lawrence, N. D. and344 Girolami, M. A., editors, Proceedings of the Fifteenth International Conference on Artificial345 Intelligence and Statistics, AISTATS 2012, La Palma, Canary Islands, Spain, April 21-23, 2012,346 volume 22 of JMLR Proceedings, pages 236–245. JMLR.org.347
Collier, O. and Dalalyan, A. S. (2013). Permutation estimation and minimax rates of identifiability.348 Journal of Machine Learning Research, W & CP 31 (AI-STATS 2013):10–19.349
Collier, O. and Dalalyan, A. S. (2016). Minimax rates in permutation estimation for feature matching.350 The Journal of Machine Learning Research, 17(1):162–192.351
Comminges, L. and Dalalyan, A. S. (2012). Tight conditions for consistency of variable selection in352 the context of high dimensionality. Ann. Stat., 40(5):2667–2696.353
Duff, I. S. and Koster, J. (2001). On algorithms for permuting large entries to the diagonal of a sparse354 matrix. SIAM Journal on Matrix Analysis and Applications, 22(4):973–996.355
Flammarion, N., Mao, C., and Rigollet, P. (2019). Optimal rates of statistical seriation. Bernoulli,356 25(1):623–653.357
Gao, C. and Zhang, A. Y. (2019). Iterative algorithm for discrete structure recovery. arXiv preprint358 arXiv:1911.01018.359
Harwood, B. and Drummond, T. (2016). Fanng: Fast approximate nearest neighbour graphs. In 2016360 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 5713–5722.361
Ingster, Y. I. (1982). Minimax nonparametric detection of signals in white Gaussian noise. Probl. Inf.362 Transm., 18:130–140.363
Ingster, Y. I. and Suslina, I. A. (2003). Nonparametric goodness-of-fit testing under Gaussian models,364 volume 169 of Lecture Notes in Statistics. Springer-Verlag, New York.365
Itseez (2015). Open source computer vision library. https://github.com/itseez/opencv.366
Jiang, Z., Xie, L., Deng, X., Xu, W., and Wang, J. (2016). Fast nearest neighbor search in the hamming367 space. In Tian, Q., Sebe, N., Qi, G.-J., Huet, B., Hong, R., and Liu, X., editors, MultiMedia368 Modeling, pages 325–336, Cham. Springer International Publishing.369
Jin, Y., Mishkin, D., Mishchuk, A., Matas, J., Fua, P., Yi, K. M., and Trulls, E. (2020). Image370 Matching across Wide Baselines: From Paper to Practice. International Journal of Computer371 Vision.372
Juditsky, A. and Nemirovski, A. (2020). Statistical inference via convex optimization. Princeton, NJ:373 Princeton University Press.374
Kuhn, H. W. (1955). The Hungarian Method for the Assignment Problem. Naval Research Logistics375 Quarterly, 2(1–2):83–97.376
Kuhn, H. W. (2010). The Hungarian Method for the Assignment Problem, pages 29–47. Springer377 Berlin Heidelberg, Berlin, Heidelberg.378
Laurent, B. and Massart, P. (2000). Adaptive estimation of a quadratic functional by model selection.379 Ann. Statist., 28(5):1302–1338.380
Lowe, D. G. (2004). Distinctive image features from scale-invariant keypoints. International journal381 of computer vision, 60(2):91–110.382
Ma, J., Jiang, X., Fan, A., Jiang, J., and Yan, J. (2021). Image matching from handcrafted to deep383 features: A survey. International Journal of Computer Vision, 129.384
Ma, R., Tony Cai, T., and Li, H. (2020). Optimal permutation recovery in permuted monotone matrix385 model. Journal of the American Statistical Association, pages 1–15.386
Malkov, Y. A. and Yashunin, D. A. (2020). Efficient and robust approximate nearest neighbor search387 using hierarchical navigable small world graphs. IEEE Transactions on Pattern Analysis and388 Machine Intelligence, 42(4):824–836.389
Mao, C., Pananjady, A., and Wainwright, M. J. (2020). Towards optimal estimation of bivariate390 isotonic matrices with unknown permutations. Annals of Statistics, 48(6):3183–3205.391
Mao, C., Weed, J., and Rigollet, P. (2018). Minimax rates and efficient algorithms for noisy sorting.392 In Algorithmic Learning Theory, pages 821–847. PMLR.393
Munkres, J. (1957). Algorithms for the assignment and transportation problems. Journal of the394 Society for Industrial and Applied Mathematics, 5(1):32–38.395
Ndaoud, M. and Tsybakov, A. B. (2020). Optimal variable selection and adaptive noisy compressed396 sensing. IEEE Trans. Inf. Theory, 66(4):2517–2532.397
Pananjady, A. and Samworth, R. J. (2020). Isotonic regression with unknown permutations: Statistics,398 computation, and adaptation. arXiv preprint arXiv:2009.02609.399
Pananjady, A., Wainwright, M. J., and Courtade, T. A. (2017). Linear regression with shuffled data:400 Statistical and computational limits of permutation recovery. IEEE Transactions on Information401 Theory, 64(5):3286–3300.402
Ramdas, A., Isenberg, D., Singh, A., and Wasserman, L. A. (2016). Minimax lower bounds for403 linear independence testing. In IEEE International Symposium on Information Theory, ISIT 2016,404 Barcelona, Spain, July 10-15, 2016, pages 965–969. IEEE.405
Rublee, E., Rabaud, V., Konolige, K., and Bradski, G. (2011). Orb: An efficient alternative to sift or406 surf. In 2011 International Conference on Computer Vision, pages 2564–2571.407
Shah, N. B., Balakrishnan, S., and Wainwright, M. J. (2021). A permutation-based model for408 crowd labeling: Optimal estimation and robustness. IEEE Transactions on Information Theory,409 67(6):4162–4184.410
Slawski, M. and Ben-David, E. (2019). Linear regression with sparsely permuted data. Electronic411 Journal of Statistics, 13(1):1–36.412
Tian, Y., Balntas, V., Ng, T., Barroso-Laguna, A., Demiris, Y., and Mikolajczyk, K. (2020). D2d:413 Keypoint extraction with describe to detect approach. In Proceedings of the Asian Conference on414 Computer Vision (ACCV).415
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E.,416 Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J.,417 Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y.,418 Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero,419 E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy420 1.0 Contributors (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.421 Nature Methods, 17:261–272.422
Wang, K., Zhu, N., Cheng, Y., Li, R., Zhou, T., and Long, X. (2018). Fast feature matching based on423 r -nearest k -means searching. CAAI Transactions on Intelligence Technology, 3(4):198–207.424
Wang, X. (2011). A fast exact k-nearest neighbors algorithm for high dimensional search using425 k-means clustering and triangle inequality. In The 2011 International Joint Conference on Neural426 Networks, pages 1293–1299.427
Wei, Y., Wainwright, M. J., and Guntuboyina, A. (2019). The geometry of hypothesis testing over428 convex cones: Generalized likelihood ratio tests and minimax radii. The Annals of Statistics,429 47(2):994 – 1024.430
Wolfer, G. and Kontorovich, A. (2020). Minimax testing of identity to a reference ergodic markov431 chain. In Chiappa, S. and Calandra, R., editors, The 23rd International Conference on Artificial432 Intelligence and Statistics, AISTATS 2020, 26-28 August 2020, Online [Palermo, Sicily, Italy],433 volume 108 of Proceedings of Machine Learning Research, pages 191–201. PMLR.434
Xing, X., Liu, M., Ma, P., and Zhong, W. (2020). Minimax nonparametric parallelism test. Journal435 of Machine Learning Research, 21(94):1–47.436
Checklist437
1. For all authors...438 (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s439 contributions and scope? [Yes] See Section 1.440 (b) Did you describe the limitations of your work? [Yes] See Section 1 and Section 2 for441 informal statements, precise assumptions/limitations are presented in Section 3.442 (c) Did you discuss any potential negative societal impacts of your work? [N/A]443 (d) Have you read the ethics review guidelines and ensured that your paper conforms to444 them? [Yes]445 2. If you are including theoretical results...446
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section 3.447 (b) Did you include complete proofs of all theoretical results? [Yes] See supplemental448 material (Appendix B).449 3. If you ran experiments...450
(a) Did you include the code, data, and instructions needed to reproduce the main experi-451 mental results (either in the supplemental material or as a URL)? [Yes] See supplemen-452 tal material.453
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they454 were chosen)? [Yes] See Section 5 and supplemental material (Appendix D).455
(c) Did you report error bars (e.g., with respect to the random seed after running experi-456 ments multiple times)? [Yes] See Section 5 and supplemental material (Appendix D).457
(d) Did you include the total amount of compute and the type of resources used (e.g., type458 of GPUs, internal cluster, or cloud provider)? [N/A]459
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...460 (a) If your work uses existing assets, did you cite the creators? [Yes]461 (b) Did you mention the license of the assets? [N/A]462 (c) Did you include any new assets either in the supplemental material or as a URL? [N/A]463
464
(d) Did you discuss whether and how consent was obtained from people whose data you’re465 using/curating? [N/A]466
(e) Did you discuss whether the data you are using/curating contains personally identifiable467 information or offensive content? [N/A]468
5. If you used crowdsourcing or conducted research with human subjects...469 (a) Did you include the full text of instructions given to participants and screenshots, if470 applicable? [N/A]471 (b) Did you describe any potential participant risks, with links to Institutional Review472 Board (IRB) approvals, if applicable? [N/A]473 (c) Did you include the estimated hourly wage paid to participants and the total amount474 spent on participant compensation? [N/A]475

Review summaries:
This paper studies the problem of finding matches between two observations, one of which has outliers. It provides theoretical analysis about the conditions when using the Least Sum of Logarithms estimator can lead to optical solution. It studies different situation with regards to different cases of observations, providing the upper bound and lower bound theoretically. Experiments are conducted on synthetic datasets and one real image matching dataset.

This paper considers the problem of matching two sets of vectors where the second set consists of a noisy version of the first set (inliers) together with outliers. Under a statistical framework, the paper develops theoretical bounds on the inlier-inlier and inlier-outlier separation that enable correct matching for a class of estimators including the Least-Sum-of-Logarithms (LSL). The approach follows Collier and Dalalyan, 2016, who considered the same problem without outliers. The LSL estimator is also evaluated using synthetic data, as well as in a semi-synthetic computer vision experiment using the IMC-PT-2020 real dataset for the task of feature matching towards obtaining point correspondences from two views.  

This paper provides a thorough investigation of how accurate detection regions can be mapped to each other.  The detection regions are typically represented by noisy vectors with outliers, and those problems can be formulated as how to finding the matching map between two sets of d-dimensional vectors. The investigation shows that the LSL estimator is robust to noise and outliers for both problems of detection and estimation.