Generalized Shuffled Linear Regression
Feiran Li1

Kent Fujiwara2
1

Fumio Okura1
2

Osaka University

Abstract
We consider the shuffled linear regression problem where
the correspondences between covariates and responses are
unknown. While the existing formulation assumes an ideal
underlying bijection in which all pieces of data should match,such an assumption barely holds in real-world applications
due to either missing data or outliers. Therefore, in this work,
we generalize the formulation of shuffled linear regression
to a broader range of conditions where only part of thedata should correspond. Moreover, we present a remarkably
simple yet effective optimization algorithm with guaranteed
global convergence. Distinct tasks validate the effectivenessof the proposed method.1

Yasuyuki Matsushita1

LINE Corporation

All these issues lead to a recently raised variant of linear
regression, namely, the shuffled linear regression (SLR) problem. It aims at simultaneously recovering both correspondences between data and the regression variable x [28, 45].
Formally, its objective function is written as
  \begin {aligned} \min _{\V {P} \in \Pi ^\star , \ \V {x} \in \mathcal {X}} \norm {\V {PAx-b}}_2^2, \end {aligned} \label {SLR_obj} (2)where P is an n-dimensional square permutation matrix
from the discrete and finite set Π⋆ satisfying
  \begin {aligned} \V {P1}=\V {P}^T\V {1}=\V {1}, \ \V {P}_{ij}\in \{0,1\}, \end {aligned} \label {constraint_on_p_star} 

(3)where 1 is the all-one vector, and Pij denotes the ij th entry of P. Intuitively, constraints (3) bring some plausible
properties on the relations between each covariate Ai andresponse bj :1. Introduction• Pij ∈ {0, 1} ensures Ai and bj are either absolutely
matched or completely independent from each other.

Linear regression has always been a powerful tool forparameter estimation problems in science and engineering.
In the conventional form, it is designed to estimate the parameters of a linear system given pairs of matched data. Its
objective function is written as2• Together with the binary property, P1 = PT 1 = 1 further guarantees bijectivity, i.e., each Ai and bj possess
exactly one match.

  \min _{\V {x}\in \mathcal {X}} \norm {\V {Ax-b}}_2^2, (1)

where A ∈ Rn×d (n ≫ d) and b ∈ Rn denote the collections of covariates and responses, respectively, and x is theregression variable belonging to set X .
There also exist cases where the correspondences between
A and b are unknown. For example, in computer vision tasks
such as point cloud registration, shape matching, perspectiven-points, and image registration, the regression variable x
stands for the rigid transformation, spectral mapping, camera projection, and transformation matrices, respectively. Insuch cases, correspondences among points, vertices, or pixels are often unknown. Similar situations also appear inother fields, such as data de-anonymization [25] in data science, artifacts dating [12] in archaeology, and sampling withjitter [28] in signal processing.
1 Source code can be found at https://github.com/SILI1994/

Generalized-Shuffled-Linear-Regression.
2 Without loss of generality, we assume the Euclidean norm as a metric.

This vanilla formulation of SLR demands exactly the same
cardinalities of {Ai } and {bj }, and all of them to be inliers.
However, such requirements are seldom satisfied in practicedue to missing data and outliers. Therefore, it is desired toextend SLR to work well for real-world applications.
In this work, we generalize SLR to a broader range ofcases where the cardinalities of covariates {Ai } and responses {bj } can be different, and only part of them shouldmatch. We denote this generalized setup of SLR as Generalized Shuffled Linear Regression (GSLR). Our contributions
are summarized as follows:• We present GSLR, a generalized formulation of SLR,
making it applicable to practical situations, where only
parts of the data should correspond.
• We propose a remarkably simple yet effective algorithmwith a detailed theoretical analysis for optimization.
• We employ distinct examples to demonstrate how
GSLR can benefit computer vision tasks in achieving
state-of-the-art accuracies.

!! !#"

Ground- !!
truth: !"#""(!!)
!"!
!"#

!!
!#
""(!#)

!! !#

NN:

!"!
!"#

Figure 1: NN-based corresponding strategy may violate
the one-to-one correspondences by either assigning a single
point to multiple queries or leaving points unmatched.

2. Related work
In this section, we review the prior-arts related to (G)SLR
and discuss about relevant techniques and problems.
Shuffled linear regression Some pioneering works of
SLR focus on solvability analysis. For example, Unnikrishnan et al. [37] prove that for a given set of covariatesA ∈ Rn×d , n ≥ 2d is necessary for obtaining a uniquesolution of the regression variable x in noiseless case. Pananjady et al. [28] further provide conditions for unique recovery
under noisy data. There also exists a work [19] that provides
theoretical guidance on the design of global algorithms withpolynomial time complexity.
However, they are still imprac
ticable due to the O nd time complexity.
Other literature focuses on practically solving the SLR
problem. For example, Abid et al. [1, 2] establish a MonteCarlo expectation-maximization framework by treating correspondences as the hidden parameter. Zhang et al. [45]reformulate objective (2) by eliminating x with maximum
likelihood criterion, which in fact leads to a quadratic assignment problem (QAP).
There also exist works that concentrate on simplified
versions of objective (2). For example, Ashwin et al. [27]estimate PT Ax as a whole in the manner of the denoising problem. Slawski and Ben-David [34] assume only asmall amount of data is disordered and propose to model thesparsely-shuffled property via L1-norm. Based on the same
assumption, they later suggest rejecting the mismatched ones
as outliers via robust estimation [35].
In contrast to these approaches that assume an ideal underlying bijection, our GSLR formulation aims at generalcases where only part of the data should correspond.

such as iterative closet point (ICP) [3] and its variants, areinfeasible algorithms for the SLR problem. In detail, as
shown in Fig. 1, their nearest neighbors (NN)-based matching strategy may either link a point to multiple queries orleave it unmatched, which consequently violates the one-toone matching constraint of permutation matrices. As we willshow later in Sec. 5, it is this violation that prevents themfrom achieving stably accurate results. For a solution, someliterature proposes to strengthen the bijectivity. For example,
both Jing et al. [30] and Szymon [32] propose to symmetrize
the registration process; and Gold et al. [14] convert theassignment matrix to a doubly stochastic one to balance theweights of correspondences. However, contrary to our GSLR
formulation, these attempts still cannot guarantee strict oneto-one correspondences. There also exist works [38, 39] thatcast the matching problem into a linear assignment problem(LAP) by maximizing kernel densities in product spaces,
whereas their generalizability is limited for assuming thecorrespondences between the covariates {Ai } and responses{bj } to be partial-to-all instead of partial-to-partial.
Robust estimation Robust estimation aims at alleviating
the effects of outliers on regression problems. A popular tool
to realize so are the M-estimators [36], which propose toeither down-weight or completely reject suspected outliers.
The ability of a robust estimator is typically assessed by two
criteria, namely, breakdown point and efficiency [48]. Specifically, the breakdown point, theoretically upper-bounded by
0.5, demonstrates the proportion of outliers an estimator cantolerate before providing an incorrect result. The efficiency,
calculated as the ratio of the theoretically minimal varianceprovided by the Cramer-Rao bound w.r.t. the actual one of
the estimator, stands for the quality of an estimator.
Relation to QAP QAP (e.g., second-order graph matching [46]) is another popular problem subject to permutational
constraints. However, different from QAP, (G)SLR additionally requires to estimate the regression variable x. Thus,
converting (G)SLR to QAP is hardly possible unless the
optimal x can be explicitly expressed by A, b, and P [45](e.g., when x ∈ Rd , we can eliminate x from (G)SLR by
−1 T Tx = AT PT PA
A P b). However, in general cases
where x ∈ X ≠ Rd , such an elimination can only present
coarse results as mentioned in [22].

3. Formulation and optimization of GSLR
Similar problems in applications Although it is only recently that SLR was formally formulated [28], the problem
itself has already been approximately handled in various
applications. For example, the rigid point cloud registration
problem, which aims at recovering the spatial transformation
T with unknown point-wise correspondences, can perfectly
fit in the (G)SLR formulation. However, existing approaches,

Following the formulation of SLR shown in objective (2),the GSLR problem can be similarly written as
  \begin {aligned} \min _{\V {P} \in \Pi , \ \V {x} \in \mathcal {X}} \norm {\V {PAx-b}}_2^2, \end {aligned} \label {obj_mat_form} (4)where A ∈ Rm×d , b ∈ Rn (m ̸= n in general), and Π isthe set of generalized permutation matrices.

3.1. Properties of ΠGiven objective (4), we first explore the properties thata generalized permutation matrix P ∈ Π should satisfy.
Specifically, apart from the binary and one-to-one matchingproperties inherited from the ordinary counterpart Π⋆ , wealso desire it to be capable of characterizing outliers, as• Pij = 1 if Ai and bj should match otherwise 0.
• Pij = 0 if either Ai or bj is an outlier.
PP•i Pij ,
j Pij ∈ {0, 1} to guarantee that each Ai
and bj can have one correspondence at most.
These requirements can be summarized as  \begin {aligned} \V {P1} \leq \V {1}, \ \V {P}^T\V {1} \leq \V {1}, \ \V {P}_{ij}\in \{0,1\}. \end {aligned} \label {constraint_1}   \begin {aligned} \min _{\V {P} \in \Pi , \ \V {x} \in \mathcal {X}} \sum \V {P}_{ji} \norm {\V {A}_i\V {x} - \V {b}_j}_2^2. \end {aligned} \label {final_obj_func} 

(8)

If we apply alternative optimization to objective (8), bytreating the correspondences Pji as known, updating theregression variable x leads to a weighted linear least square
problem, which should be globally solvable in most applications. On the other hand, if x is fixed, the objective function
for updating P is in the form of
  \min _{\V {P} \in \Pi } \V {D}_{ij} \V {P}_{ji}, \label {klap_obj_func} (5)However, conditions (5) cannot be directly used for optimization. If done so, the global minimum of objective (4)would be anchored to 0 with the associated trivial solution
P = 0. Therefore, it is necessary to introduce some otherconstraints to ensure the existence of matches.
A simple yet effective solution is to manually designate
the number of
we introduce an
h matches [40,41]. Specifically,iP
integer k ∈ min(m,n)
,
min
(m,
n)
and
force
Pij = k,2
indicating that there must exactly exist k inlier matches
between the sets of covariates and responses. Consequently,
the trivial solution can be avoided, and the constraints overthe set of generalized permutation matrices Π is  \V {P1} \leq \V {1}, \ \V {P}^T\V {1} \leq \V {1}, \ \V {P}_{ij}\in \{0,1\}, \ \sum \V {P}_{ij} = k. \label {final_constraint} 

The key observation of our proposal is that slight reformulation of objective (7) can result in a bi-globally optimizable
form. Specifically, per the definitions of linear regression,
we can facilely rewrite it into the summation form of(6)

For notational simplicity, we hereafter keep using Π to denote the set of all the generalized permutation matrices under
all possible k. Similar to the ordinary counterpart Π⋆ , constraints (6) imply that Π is discrete and combinatorial, anddegenerates to Π⋆ when k = m = n. Now, we can obtainthe formal formulation of GSLR by joining (4) and (6):  \begin {aligned} &\min _{\V {P, \ x \in \mathcal {X}}} \norm {\V {PAx-b}}_2^2 \qquad \\ \mathrm {s.t.} \ \V {P1} \leq \V {1}, \ &\V {P}^T\V {1} \leq \V {1}, \ \V {P}_{ij}\in \{0,1\}, \ \sum \V {P}_{ij} = k. \label {final_obj} \end {aligned} (7)3.2. A bi-global optimization algorithm
In this section, we consider the optimization method for
objective (7). Its difficulties are threefold. First, the discrete
property of the permutation P hinders optimizers designed
for smooth functions. Second, its high non-convexity requires careful initialization. Last, the number of matches k
is unknown. As a solution, we propose a simple and effective
algorithm, which can effectively maintain the discretenessand alleviate the initialization problem. For clarification,
we hereby temporally assume that k is known and put its
estimation in Sec. 3.3.

(9)2

where D is the cost matrix formulated by ∥Ai x − bj ∥2 .
When the number of matches k = min (m, n), problem(9) is a standard LAP, which can be globally optimized via
efficient approaches such as the Hungarian algorithm [20].
For an LAP with arbitrary k, it forms a k-cardinality linear
assignment problem (k-LAP), which can be transformed intoa standard one via paddling the cost matrix [40]. Therefore,
we can conclude that the GSLR problem is facilely solvable
via alternative optimization.

3.3. k as an M-estimator
It is obvious that the number of matches k characterizes
the ability of GSLR in rejecting outliers. In detail, it acts similarly to the trimming operator in robust estimation, whichsets the weights of a specific proportion of suspected outliers
to zero. In fact, such a technique has already been broadly
used in various similar tasks. For example, the Trim-ICP
algorithm [8] for point cloud registration suggests removing
a portion of paired points with large residuals; and the partial graph matching algorithm [41] assumes only half of thenodes as inliers, which leads to k = 21 min (m, n).
However, we hesitate to directly use the trimming operator for its well-studied limitation on simultaneously maintaining a high breakdown point and high efficiency. Specifically,
although a breakdown point of 0.5 can be facilely achievedby ignoring 50% of the data, the associated statistical efficiency would only be 7% [11]. Moreover, it is also preferable
for some applications [26, 31] to preserve as many matches
as possible while maintaining high robustness.
In this work, we employ another M-estimator to determine k, namely, the Huber-skip estimator [16]. Apart frominheriting the complete-rejection (i.e., set weights as binary)and easy-to-calculate properties of the trimming operator, italso achieves much higher efficiency while being capable ofmaintaining a breakdown point of 0.5. Specifically, given
some data (r1 , r2 , ..., rn ), it rejects ri that violates
  \abs {\frac {r_i - \mathrm {median}\left (r\right )}{\mathrm {MAD}\left (r\right )}} \leq \text {threshold}, \label {huber_skip} (10)

Algorithm 1: A bi-global algorithm for GSLR
Input: A ∈ Rm×d , b ∈ Rn
Initialization: x = x0 ∈ X
while not converged do• Calculate P̆t by solving the standard LAP
with kt = min (m, n) and xt−1 .
• Compute residuals under P̆t and xt−1 andapply the Huber-skip estimator to get kt .

since the linear least square problem on the regression variable x is globally optimizable, it is easy to see that  E\left (\V {P}_t, \V {x}_{t+1}\right ) \leq E\left (\V {P}_t, \V {x}_{t}\right ), \label {mono_eq1} (11)where the equality holds if and only if xt is already a
global optimum. Furthermore, we define a temporary permutation matrix Ptmp by randomly removing kt+1 − kt
non-zero assignments from Pt . Since the cost matrix
2
Dij = (Ai x − bj ) is non-negative, we can obtain that

• Set kt = min (kt , kt−1 )• Calculate P̂t by solving the k-LAP with theupdated kt and xt−1 .
• Get xt ∈ X by solving the weighted linear
regression problem under P̂t .
Output: x ∈ X , P ∈ Π

where MAD is the median absolute deviation, and the threshold is set to 3.5 throughout the paper as suggested by [16].
More details are provided in the supplementary material.
To embed the Huber-skip estimator into our aforementioned optimization scheme, we directly follow the formulation of least trimmed squares (LTS) [18] and compute k as
the number of pairs whose residuals satisfy Eq. (10). However, such a naive formulation would break the monotonicityof the alternative optimization scheme and make its convergence non-guaranteed. As a solution, in the tth iteration, we
set kˆt = min (kt−1 , kt ) to ensure a convergence. The entire
algorithm is summarized in Algorithm 1.

4. Algorithm analysis
In this section, we analyze the time and space complexities of Algorithm 1. Moreover, we also prove that it monotonically decreases the objective function and converges
globally to a critical point.

4.1. Complexity
We analyze the space and time complexities of Algorithm 1 under the assumption that #{A1 , . . . , Am } ≥#{b1 , . . . , bn } (i.e., m ≥ n) without loss of generality.
The main space complexity comes from the cost matrix D,
which is O (mn) theoretically and O (mnd) in our implementation to avoid loops. As for time complexity in each
iteration,
the Hungarian algorithm for solving the k-LAP is2
O n (m + n − k) , and the Huber-skipping proceduresis with O (n log n) time complexity.

4.2. Monotonicity
We now show that the total energy E (P, x) of objective(8) monotonically decreases with Algorithm 1. Specifically,  E\left (\V {P}_{\mathrm {tmp}}, \V {x}_{t+1}\right ) \leq E\left (\V {P}_{t}, \V {x}_{t+1}\right ). \label {mono_eq2} (12)Since Ptmp satisfies the constraints of the k-LAP formulated
by xt+1 and kt+1 , it is also within the solution space. On
the other hand, as Pt+1 is a global optimum of this objectivefunction, we can further derive that  E\left (\V {P}_{t+1}, \V {x}_{t+1}\right ) \leq E\left (\V {P}_{\mathrm {tmp}}, \V {x}_{t+1}\right ), \label {mono_eq3} (13)where the equality holds if and only if Ptmp is already within
the set of global minima. Moreover, we can obtain thefollowing relationship by summarizing Eqs. (11)-(13):  E\left (\V {P}_{t+1}, \V {x}_{t+1}\right ) \leq E\left (\V {P}_{\mathrm {tmp}}, \V {x}_{t+1}\right ) \leq E\left (\V {P}_{t}, \V {x}_{t}\right ), \label {monotone}  (14)where the equality holds if and only if (Pt , xt ) is an critical point of the objective function. From Eq. (14), we canconclude that Algorithm 1 has a monotonic property.

4.3. Convergence
Our proposal can globally converge to an critical point
owing to the monotonicity. To be precise, we first make two
assumptions about the domain of the regression variable x:• The set X for x ∈ X is compact.
• f (x) = Ax is continuous on X .
These assumptions are rather weak and can be satisfied inmost practical applications formulated in the GSLR manner.
For example, the orthogonal groups raised in geometric computer vision and geometric processing tasks can satisfy bothconditions. For the general real coordinate space Rn , despite
not being compact by itself, we can reasonably assume thatx lies in a closed and bounded subset of it, since infinities
are barely meaningful solutions in practical applications.
Similar to the proof of the global convergence of the
expectation-maximization algorithm [42] and the generalized alternative optimization [15], the global convergenceproperty of GSLR is a direct instantiation of the convergence
theorem of Zangwill [44]:Theorem 1 (Zangwill’s global convergence theorem). Let
H be a set-valued map on S that, given s0 ∈ S, generates
a sequence {sk }∞k via iterating sk+1 ∈ H (sk ). Also, let asolution set Γ ⊂ S be given. Suppose that′′
• The sequence {sk }∞k ⊂ S for S ⊂ S is compact.

• There exists a continuous function Ψ on S that
– if s ∈/ Γ: Ψ (s′ ) < Ψ (s) for all s′ ∈ H (s).

We now prove that the map (Pt+1 , xt+1 ) ∈ H (Pt , xt )is closed on Π × X . Specifically, for any (Pt , xt ), we canalways obtain some δ1 > 0 and δ2 > 0 that satisfy
  \begin {gathered} \forall \V {P}', \ \V {P}' \to \V {P} \coloneqq \{\V {P}' : \V {P}' \in \mathcal {N}_{\delta _1}(\V {P}_t)\} = \{\V {P}_t\}, \\ \forall \V {x}', \ \V {x}' \to \V {x} \coloneqq \abs {\V {x'}, \V {x}_t}_\V {x} < \delta _2, \end {gathered} \label {continuity_condition} – if s ∈ Γ: Ψ (s′ ) ≤ Ψ (s) for all s′ ∈ H (s).
• The set-valued map H is closed at all s for s ∈ S \ Γ.
Then, the limit of any convergent subsequence of {sk }∞k is inthe solution set Γ. Furthermore, limk→∞ Ψ (sk ) = Ψ (s⋆ )for all limit points s⋆ .
An iterative optimization algorithm can be described asa triplet {Ψ, H, Γ} in the terminology of Theorem 1. For
Algorithm 1, we can straightforwardly set the objective function of GSLR as Ψ, and H : Π×X → Π×X to one iterationof it. Also, we can formulate the solution set with all thecritical points that satisfy
  \begin {aligned} \Gamma = \{ \left (\V {P, x}\right ) \in \Pi \times \mathcal {X} : \V {P} &\in \argmin _\V {P'}\Psi \left (\V {P', x}\right ) \mathrm {and} \\ \V {x} &\in \argmin _\V {x'} \Psi \left (\V {P, x'}\right ) \}. \end {aligned} \label {set_of_extreme_point} 
(15)

We show in the following that all the three requests in Theorem 1 are satisfied by this GSLR triplet.
For the compactness of the feasible region Π × X , weonly need to verify the compactness of Π and X respectively,
since unions of compact sets are still compact. Specifically,
X is compact under our assumption; and for Π, its finite anddiscrete properties imply its compactness.
For the continuity and monotonicity of the objective func2
tion Ψ = ∥PAx − b∥2 , since it is already demonstrated tobe monotonic in Sec. 4.2, we hereby only focus on the continuity. For proof, since compositions of continuous functions
are still continuous, we hereby focus on the continuity off (P, x) = PAx. Given the fact that such a function is abilinear form w.r.t. (P, x), we only need to show that it isseparately continuous w.r.t. both P and x. For proof, per
the assumption mentioned above, we have f (x) = Ax
continuous w.r.t. x ∈ X . Moreover, since the set of all possible generalized permutations Π is finitely discrete, we havef (P) = PA continuous w.r.t. P ∈ Π. Therefore, we obtainthat f (P, x) = PAx is continuous and hence so is Ψ.
Definition 1 (Closedness of set-valued map). A set-valued
map H : U → V is said to be closed on a point u ∈ Uprovided that the following two characters(16)

where Nδ1 (Pt ) is the neighborhood of Pt . Therefore, we
can derive the differences of mapped x as  \begin {split} &\mathcal {H}\left (\V {P', x'}\right ) \setminus \mathcal {H}\left (\V {P}_t, \V {x}_t\right ) = \mathcal {H}\left (\V {P}_t, \V {x}_t\right ) \setminus \mathcal {H}\left (\V {P', x'}\right ) \\ =& \mathcal {H}\left (\V {P}_t, \V {x}'\right ) \setminus \mathcal {H}\left (\V {P}_t, \V {x}_t\right ) \\ =&\argmin _{\V {x} \in \mathcal {X}}\Psi \left (\V {P}_t, \V {x}\right ) \setminus \argmin _{\V {x} \in \mathcal {X}}\Psi \left (\V {P}_t, \V {x}\right ) = \emptyset . \end {split} \raisetag {3\normalbaselineskip } \label {continuity_x} (17)

Since Eq. (17) implies that the cost matrix for the k-LAP
generated by the regression variable x ∈ H (P, x) remains
unchanged when the input of H switches from (Pt , xt )to (P′ , x′ ), so will the mapped P. Therefore, we have
H (Pt , xt ) closed by further joining Eqs. (16) and (17).
As all the requirements of the Zangwill’s convergence
theorem are satisfied as presented above, we can now conclude that our algorithm is globally convergent to a criticalpoint as defined by the solution set Γ.

5. Experiments and applications
We conduct experiments to demonstrate the differences
between SLR and GSLR in terms of robust estimation, as
well as to explore how GSLR performs on suitable applications. We limit the discussion within the computer vision
field, although GSLR itself is a much more general problem.

5.1. Task 1: Image registration
We first carry out experiments on image registration todemonstrate the effectiveness of GSLR over SLR in termsof robust estimation. Given a pair of images portraying thesame plane, or a non-planar scene subject to 0 translation interms of the camera pose, homography describes the relationof p′ = αHp, where p′ and p are matched pixels on the two
images in homogeneous coordinates, α is a scaling factor,
and H is the homography matrix. Furthermore, if images arecaptured with long focal lengths, H can be assumed as affine
transformations [17]. Since such an assumption implies thatthe scaling α = 1, we can formulate the problem as  \min _{\V {H} \ \V {P} \in \Pi } \norm {\V {PHA-B}}_F^2, \label {homography} (18)• uk → u, where uk ∈ U,
• vk → v, where v, vk ∈ V,
• vk ∈ F(uk ),imply v ∈ H (u). Based on this, H is called closed if it isclosed at any u ∈ U.

where A and B stand for the collections of pixels from thetwo images in homogeneous coordinates. This problem isregarded as a GSLR problem rather than an SLR one due tothe inherent existence of outliers.
In the implementation, we use the SIFT feature [21] toretrieve the sets of pixel locations A and B. The initial

Source
SourceTargetTargetFeature matching

Ground-truth
w/ RANSACSLRSLRGSLR (ours)
GSLR(ours)

Figure 2: Results of image registration experiments. GSLR
achieves visually the same results w.r.t. the ground-truths
whereas SLR presents obvious distortions on all the trials.
On the source and target pairs, outliers detected by GSLR
are marked as red and inliers as green.
transformation H of all the trials is set to the identity matrix. For optimizing SLR, we fix the number of inliers kof Algorithm 1 to the number of feature points. We use
the mature descriptor and RANSAC-based method [17] toapproximate the ground-truths. Fig. 2 shows the results.
As shown, GSLR can effectively estimate the transformation matrices and provide visually the same results as thedescriptor-based approach (denoted as Ground-truth) even
without using the descriptors, where SLR results are obviously affected by outliers on all the trials.

Outliers

OriginPartial-overlap

MeanMedianMeanMedian

ICP

1.05

1.04

1.70e-1

1.19e-1

Trim-ICP

1.07

1.09

8.61e-1

8.56e-1

SVR

1.05

1.00

2.16e-1

1.71e-1CPD

1.58e-2

1.52e-2

9.03e-2

9.09e-2

FilterReg

2.20e-2

1.70e-2

1.12e-1

9.95e-2

FGR

8.55e-2

7.22e-2

1.07

1.12

GSLR (ours)

6.24e-3

5.65e-3

6.52e-2

3.78e-2ICP

CPD

FilterReg

GSLR (ours)

Figure 3: Results on the synthetic data. Upper left: Robustness of GSLR under the outlier-contaminated setup. Identified outliers are marked as green. Upper right: Summary of
the mean and median transformational errors. Bottom: An
example of the partially overlapped data. Trim-ICP, SVR,
and FGR present similar visualizations and thus omitted.

scaled Geman-McClure estimator to remove outliers frommatchings obtained via FPFH features [33], and CPD and FilterReg cast the registration problem into a distribution-fitting
one and use a uniform distribution to cope with outliers.
Other details are given in the supplementary material.

5.2. Task 2: Point cloud registration
Point cloud registration works toward estimating the optimal rigid transformation T between two point clouds. Its
objective function can be written as  \min _{\V {T} \in \mathrm {SE}\left (3\right ), \ \V {P} \in \Pi } \norm {\V {PTA - B}}_F^2, \label {regis_obj} 

(19)

where A and B represent the two point clouds in homogeneous coordinates. Usually, they would consist of numerous
outliers that should never be used for estimating T.
Experimental setup, metrics and peer methods We conduct experiments on both synthetic and real-world datasets.
All pairs of point clouds are registered back and forth dueto the unknown alignment order. For space limitation,
we only report the transformational errors with the metricTrans. err = Tgt T−1
esti − I F in the main paper and put
the respective rotational and translational counterparts in the
supplementary material. For competing methods, we use
ICP [3], Trim-ICP [8], SVR [7], CPD [24], FGR [47], andFilterReg [13] since they represent different strategies on handling outliers. Specifically, ICP works as the basic baseline,
Trim-ICP shares similar outlier-rejection techniques with
GSLR but employs the NN-based matching strategy, SVR
treats both point clouds as mixture distributions and minimizes the robust L2-divergence between them; FGR uses a

Registration with outliers and partial overlaps We use
the synthetic bunny, dragon, and armadillo point clouds [9]for tests. For the outlier-contaminated cases, we first downsample the point clouds to approximately 2500 points andcorrupt them with Gaussian noises. Then outliers are addedwith ratios sampled from (0.25, 1), representing that 20% to
50% points are outliers. We repeat this procedure 10 times
to obtain 30 trials in total. For the partial-overlapping setup,
we down-size the point clouds to around 3000 points each,and randomly crop some parts from them. We retrieve 4
cropped point clouds from each shape and conduct 36 trials
in total. For all the tests, we set the ground-truth rotation to30◦ in Euler angles along each axis together with a randomtranslation within [0, 1]. As shown in Fig. 3, SVR fails
to provide accurate results as dense clusters of outliers areoften similarly weighted as inliers. Moreover, comparedto CPD and FilterReg, which can only down-weight theoutliers, GSLR presents more accurate results by completely
rejecting them.
Why ICP-style algorithms suffer from outliers Considering the fact that Trim-ICP and GSLR are equipped withsimilar robust estimation techniques (i.e., binarily weighting
the pairs), it is worthy of analyzing the substandard performances of Trim-ICP presented in Fig. 3.

because of repeated patterns and a minimum overlap ratioof 63.9%. To let ICP-style algorithms perform reasonably,
we roughly pre-align the point clouds by fixing the groundtruth rotations to 15◦ in Euler angles. All point clouds are
down-sampled to around 3500 points. Figure 5 reports theresults. As shown, our method exhibits the highest accuracy compared to other state-of-the-art methods. On theother hand, since the local features are extremely unreliable
when repeated patterns exist, FGR fails to present reasonable
results.
Figure 4: Sensitivity of different approaches w.r.t. initializations on outlier-contaminated point clouds.
MeanMedianICP

5.82e-1

1.65e-1

Trim-ICP

1.64

1.53

SVR

3.03e-1

2.00e-1CPD

4.97e-1

1.22e-1

FilterReg

7.59e-1

6.62e-1

FGR

1.23

1.06

GSLR (ours)

1.40e-1

1.16e-1

Figure 5: Registration results on the real-world dataset. Left:
An example achieved by our GSLR formulation. Right:
Mean and median transformational errors of all the trials.

It is the matching policy to blame. In detail, we observethat the NN-based matching strategy of ICP and its variants makes them extremely sensitive to initializations. For
demonstration, we conduct a test with different initializations
on the aforementioned bunny point clouds with 1/3 points
as outliers. The ground-truth rotations in Euler angles are
split into 6 ranges, and the registration is repeated 10 times
within each of them. Fig. 4 records the results. As shown,
ICP and Trim-ICP can also present reasonable results given
plausible initializations. In such a condition, Trim-ICP caneven provide more accurate estimations compared to the others owing to its large trimming ratio. However, they would
fail drastically when the initializations are coarsened. On thecontrary, GSLR can remarkably alleviate this problem via
the one-to-one correspondences, demonstrated by the stableperformances even under large initialization errors.
Real-world registration We also examine the performances of GSLR on real-world point clouds. For experimental setup, we establish 7 groups consisting of 5 consecutive frames from the ETH Hauptgebaude dataset [29], andconduct pairwise registration within each group, resulting
in 140 trials altogether. This dataset is known challenging

Other results and experiments We have also conducted
other experiments to analyze statistical efficiencies, elapsed
time, as well as sensitivities w.r.t. different scales of noises
and ratios of outliers. The results are provided in the supplementary material for interested readers.

5.3. Task 3: Isometric shape matching
In contrast to point cloud registration that concentrates on
the regression variable, isometric shape matching looks forvertex-wise correspondences between meshes. Following
functional map [26], its objective function can be written as  \min _{\V {C} \in \mathrm {O}\left (s\right ), \ \V {P} \in \Pi } \norm {\V {PCF - G}}_F^2, \label {sm_obj} 

(20)

where O (·) denotes the orthogonal group, and F and G
stand for the vertex-wise Laplacian eigenfunctions of thetwo meshes. In general, they include outliers since the mappings can seldom be ideally isometric. Therefore, robust
estimation is still required although the underlying groundtruth is a full-rank matching [30].
Experimental setup, metrics and peer methods We
carry out experiments on commonly used datasets. Specifically, we conduct both intra-class and inter-class matchingon the FAUST dataset [5], and intra-class matching on the
TOSCA one [6]. Inter-class tests on TOSCA are out of scope
due to their high non-isometry. For each task, we randomly
select 30 pairs of shapes and report the mean values as the
final results. For the vertices that are labeled as outliers
by GSLR, we simply map them to their nearest neighborswith the estimated C, although more complicated methods(e.g., searching for correspondences within the connectedcomponents) can be employed to improve smoothness. For
evaluation, we employ the commonly used correspondencesgeodesic metric to quantify the performances. Moreover,
we also follow the definition of coverage presented in [30]to measure the bijectivity, which is defined as the ratio ofvertices that possess at least one correspondence to their total
amount. We select ICP [26], BCICP [30] and ZoomOut [23]as the competing methods. For setup, we use 50-dimensional
eigenfunctions of the Laplace-Beltrami operator to formulate the sets of functions F and G, and initialize all the(a) Intra-class matching(b) Inter-class matching

Figure 6: Matching results on the FAUST dataset. Within each sub-figure, left shows an example achieved by GSLR, middle is
the geodesic errors, and right shows the coverage rate. For the labels, FMAP is short for functional map, and FMAP+ZO-50
and FMAP+ZO-100 respectively stand for the ZoomOut algorithm with 50 and 100-dimensional eigenfunctions.
0

Source

FMAP

ICP

BCICP

ZoomOut-100> 0.1

GSLR (ours)

Figure 7: An example of geodesic error maps on the FAUST
dataset. Comparing to the others, more correspondences
estimated by GSLR coincide with the ground-truth.

Figure 8: Matching results on the TOSCA dataset. Left:
An example matching achieved by GSLR. Right: Geodesic
errors and coverages.
37.1%

52.6%

83.4%

62.0%

FMAP

ICP

BCICP

ZoomOut-100

94.4%

algorithms with the wave kernel signature [4] and functional
map under the direct operator [30]. More details are givenin the supplementary material.
SourceResults Figure 6 shows the overall results on the FAUST
dataset, and Fig. 7 depicts the estimation errors. Comparedto the others, GSLR can provide more accurate results by
estimating obviously higher percentages of matches thatcoincide with the ground-truth. Moreover, it can also effectively maintain the one-to-one correspondences by covering
more than 80% vertices on all the pairs. For the TOSCA
dataset, we decimate the meshes to approximately 3400 vertices each and obtain the error curves via NN-composition.
Fig. 8 records the results. As presented, the GSLR formulation can still champion both the matching accuracy andthe bijectivity. We emphasize that ZoomOut-100 doubles
the dimensions of input eigenfunctions although it achieves
similar results as GSLR in terms of accuracy. An example
of coverages is given in Fig. 9.

6. Discussion and conclusion
In this work, we propose the GSLR formulation together
with an optimization algorithm to generalize SLR to broaderranging conditions where only part of the data should match.
Distinct applications demonstrate its effectiveness.
The main limitation of our current algorithm lies in time
efficiency. I.e., while the k-LAP is tractable for medium-

GSLR (ours)

Figure 9: An example of coverages on the TOSCA dataset.
Vertices violating bijection are marked as orange. Upper
percentages indicate the coverage rates.

sized problems (e.g., the sizes in our experiments vary from3k to 8k), the elapsed time would increase significantly when
dealing with giant ones. Fortunately, this problem can beeffectively solved with off-the-shelf CUDA-based Hungarianalgorithms [10, 43], which are hundreds of times faster.
For future work, GSLR requires more theoretical analysis. For example, we are interested in studying whether theconditions of the unique solution proposed in [28] for SLR
can generalize to GSLR, as well as the convergence rate ofAlgorithm 1. For applications, we wish to work on variantsof GSLR for other tasks, such as partial shape matching andregistration with point-to-plane metrics.

Acknowledgment
This work was supported by NII CRIS collaborative research program operated by NII CRIS and LINE Corporation.

References[1] Abubakar Abid, Ada Poon, and James Zou. Linear regression
with shuffled labels. arXiv preprint arXiv:1705.01342, 2017.
2[2] Abubakar Abid and James Zou. A stochastic expectationmaximization approach to shuffled linear regression. In Proceedings of Annual Allerton Conference on Communication,
Control, and Computing, 2018. 2[3] K Somani Arun, Thomas S Huang, and Steven D Blostein.
Least-squares fitting of two 3-d point sets. Transactions on
Pattern Analysis and Machine Intelligence (PAMI), 9(5):698–
700, 1987. 2, 6
[4] Mathieu Aubry, Ulrich Schlickewei, and Daniel Cremers.
The wave kernel signature: A quantum mechanical approachto shape analysis. In International conference on computer
vision workshops (ICCV workshops), 2011. 8[5] Federica Bogo, Javier Romero, Matthew Loper, and Michael J.
Black. FAUST: Dataset and evaluation for 3D mesh registration. In Proceedings of Conference on Computer Vision and
Pattern Recognition (CVPR), 2014. 7[6] Alexander M Bronstein, Michael M Bronstein, and Ron Kimmel. Numerical geometry of non-rigid shapes. 2008. 7[7] Dylan Campbell and Lars Petersson. An adaptive data representation for robust point-set registration and merging. In
Proceedings of International Conference on Computer Vision(ICCV), 2015. 6[8] Dmitry Chetverikov, Dmitry Svirko, Dmitry Stepanov, andPavel Krsek. The trimmed iterative closest point algorithm.
Object recognition supported by user interaction for service
robots, 3:545–548, 2002. 3, 6
[9] Brian Curless and Marc Levoy. A volumetric method forbuilding complex models from range images. In Proceedingsof the 23rd annual conference on Computer graphics andinteractive techniques, 1996. 6[10] Ketan Date and Rakesh Nagi. Gpu-accelerated hungarianalgorithms for the linear assignment problem. Parallel Computing, 57:52–72, 2016. 8[11] Jurgen A Doornik. Robust estimation using least trimmed
squares. Institute for Economic Modelling, Oxford Martin
School, and Economics Department, University of Oxford,
UK, 2011. 3[12] Fajwel Fogel, Rodolphe Jenatton, Francis Bach, and Alexandre d’Aspremont. Convex relaxations for permutation problems. In Proceedings of Conference on Neural Information
Processing Systems (NIPS), 2013. 1[13] Wei Gao and Russ Tedrake. Filterreg: Robust and efficient
probabilistic point-set registration using gaussian filter andtwist parameterization. In Proceedings of Conference on
Computer Vision and Pattern Recognition (CVPR), 2019. 6[14] Steven Gold, Anand Rangarajan, Chien-Ping Lu, Suguna
Pappu, and Eric Mjolsness. New algorithms for 2d and 3d
point matching: pose estimation and correspondence. Pattern
Recognition, 31(8):1019–1031, 1998. 2[15] Asela Gunawardana and William Byrne. Convergence theorems for generalized alternating minimization procedures.
Journal of Machine Learning Research, 6:2049–2073, 2005.
4[16] Frank R Hampel. The breakdown points of the mean combined with some rejection rules. Technometrics, 27(2):95–107,
1985. 3, 4[17] Richard Hartley and Andrew Zisserman. Multiple view geometry in computer vision. 2003. 5, 6
[18] Douglas M Hawkins. The feasible solution algorithm for least
trimmed squares regression. Computational statistics & data
analysis, 17(2):185–196, 1994. 4[19] Daniel J Hsu, Kevin Shi, and Xiaorui Sun. Linear regression
without correspondence. In Proceedings of Conference on
Neural Information Processing Systems (NIPS), 2017. 2[20] Harold W Kuhn. The hungarian method for the assignment
problem. Naval research logistics quarterly, 2(1-2):83–97,
1955. 3[21] David G Lowe. Object recognition from local scale-invariant
features. In Proceedings of International Conference on Computer Vision (ICCV), 1999. 5[22] Joao Maciel and Joao Paulo Costeira. A global solution tosparse correspondence problems. Transactions on PatternAnalysis and Machine Intelligence (PAMI), 25(2):187–199,
2003. 2[23] Simone Melzi, Jing Ren, Emanuele Rodolà, Abhishek
Sharma, Peter Wonka, and Maks Ovsjanikov. Zoomout: Spectral upsampling for efficient shape correspondence. Transactions on Graphics (TOG), 38(6), 2019. 7[24] Andriy Myronenko and Xubo Song. Point set registration:
Coherent point drift. Transactions on Pattern Analysis andMachine Intelligence (PAMI), 32(12):2262–2275, 2010. 6[25] Arvind Narayanan and Vitaly Shmatikov. Robust deanonymization of large sparse datasets. In Proceedings ofSymposium on Security and Privacy, 2008. 1[26] Maks Ovsjanikov, Mirela Ben-Chen, Justin Solomon, Adrian
Butscher, and Leonidas Guibas. Functional maps: a flexible
representation of maps between shapes. Transactions on
Graphics (TOG), 31(4):1–11, 2012. 3, 7
[27] Ashwin Pananjady, Martin J Wainwright, and Thomas A
Courtade. Denoising linear models with permuted data. In
Proceedings of International Symposium on Information Theory (ISIT), 2017. 2[28] Ashwin Pananjady, Martin J Wainwright, and Thomas A
Courtade. Linear regression with shuffled data: Statistical andcomputational limits of permutation recovery. Transactions
on Information Theory, 64(5):3286–3300, 2017. 1, 2, 8
[29] François Pomerleau, M. Liu, Francis Colas, and Roland
Siegwart. Challenging data sets for point cloud registration algorithms. International Journal of Robotics Research,
31(14):1705–1711, Dec. 2012. 7[30] Jing Ren, Adrien Poulenard, Peter Wonka, and Maks Ovsjanikov. Continuous and orientation-preserving correspondences via functional maps. Transactions on Graphics (TOG),
37(6):1–16, 2018. 2, 7, 8[31] Emanuele Rodolà, Luca Cosmo, Michael M Bronstein, Andrea Torsello, and Daniel Cremers. Partial functional correspondence. In Computer Graphics Forum, 2017. 3[32] Szymon Rusinkiewicz. A symmetric objective function for
icp. Transactions on Graphics (TOG), 38(4):1–7, 2019. 2

[33] Radu Bogdan Rusu, Nico Blodow, and Michael Beetz. Fast
point feature histograms (fpfh) for 3d registration. In Proceedings of International Conference on Robotics and Automation(ICRA), 2009. 6[34] Martin Slawski, Emanuel Ben-David, et al. Linear regression
with sparsely permuted data. Electronic Journal of Statistics,
13(1):1–36, 2019. 2[35] Martin Slawski, Emanuel Ben-David, and Ping Li. A twostage approach to multivariate linear regression with sparsely
mismatched data. Journal of Machine Learning Research,2019. 2[36] Robert G Staudte and Simon J Sheather. Robust estimationand testing, volume 918. 2011. 2[37] Jayakrishnan Unnikrishnan, Saeid Haghighatshoar, and Martin Vetterli. Unlabeled sensing with random linear measurements. Transactions on Information Theory, 64(5):3237–
3253, 2018. 2[38] Matthias Vestner, Zorah Lähner, Amit Boyarski, Or Litany,
Ron Slossberg, Tal Remez, Emanuele Rodola, Alex Bronstein,Michael Bronstein, Ron Kimmel, et al. Efficient deformableshape correspondence via kernel matching. In Proceedings ofInternational Conference on 3D Vision (3DV), 2017. 2[39] Matthias Vestner, Roee Litman, Emanuele Rodolà, Alex Bronstein, and Daniel Cremers. Product manifold filter: Non-rigid
shape correspondence via kernel density estimation in the
product space. In Proceedings of Conference on Computer
Vision and Pattern Recognition (CVPR), 2017. 2[40] Anton Volgenant. Solving the k-cardinality assignment problem by transformation. European Journal of Operational
Research, 157(2):322–331, 2004. 3[41] Fudong Wang, Nan Xue, Jin-Gang Yu, and Gui-Song Xia.
Zero-assignment constraint for graph matching with outliers.
In Proceedings of Conference on Computer Vision and Pattern
Recognition (CVPR), 2020. 3[42] CF Jeff Wu. On the convergence properties of the em algorithm. The Annals of statistics, pages 95–103, 1983. 4[43] Satyendra Singh Yadav, Paulo Alexandre Crisóstomo Lopes,
Aleksandar Ilic, and Sarat Kumar Patra. Hungarian algorithm
for subcarrier assignment problem using gpu and cuda. International Journal of Communication Systems, 32(4):e3884,
2019. 8[44] Willard I Zangwill. Nonlinear programming: a unified approach, volume 52. 1969. 4[45] Hang Zhang, Martin Slawski, and Ping Li. Permutation recovery from multiple measurement vectors in unlabeled sensing.
In Proceedings of International Symposium on Information
Theory (ISIT), 2019. 1, 2
[46] Ruonan Zhang and Wenmin Wang. Second-and high-ordergraph matching for correspondence problems. Transactions
on Circuits and Systems for Video Technology, 28(10):2978–
2992, 2017. 2[47] Qian-Yi Zhou, Jaesik Park, and Vladlen Koltun. Fast global
registration. In Proceedings of European Conference on Computer Vision (ECCV), 2016. 6[48] Abdelhak M Zoubir, Visa Koivunen, Yacine Chakhchoukh,
and Michael Muma. Robust estimation in signal processing:
A tutorial-style treatment of fundamental concepts. Signal
Processing Magazine, 29(4):61–80, 2012. 2