# ARCS: Accurate Rotation and Correspondence Search

Liangzu Peng
Johns Hopkins University
<EMAIL_ADDRESS>
Manolis C. Tsakiris
ShanghaiTech University
<EMAIL_ADDRESS>René Vidal
Johns Hopkins University
<EMAIL_ADDRESS>

###### Abstract

This paper is about the old Wahba problem in its more general form, which we
call “simultaneous rotation and correspondence search”. In this generalization
we need to find a rotation that best aligns two partially overlapping $3$D
point sets, of sizes $m$ and $n$ respectively with $m\geq n$. We first propose
a solver, ARCS, that i) assumes noiseless point sets in general position, ii)
requires only $2$ inliers, iii) uses $O(m\log m)$ time and $O(m)$ space, and
iv) can successfully solve the problem even with, _e.g_., $m,n\approx 10^{6}$
in about $0.1$ seconds. We next robustify ARCS to noise, for which we
approximately solve consensus maximization problems using ideas from robust
subspace learning and interval stabbing. Thirdly, we refine the approximately
found consensus set by a Riemannian subgradient descent approach over the
space of unit quaternions, which we show converges globally to an
$\varepsilon$-stationary point in $O(\varepsilon^{-4})$ iterations, or locally
to the ground-truth at a linear rate in the absence of noise. We combine these
algorithms into ARCS+, to simultaneously search for rotations and
correspondences. Experiments show that ARCS+ achieves state-of-the-art
performance on large-scale datasets with more than $10^{6}$ points with a
$10^{4}$ time-speedup over alternative methods.
https://github.com/liangzu/ARCS

##  1 Introduction

> The villain Procrustes forced his victims to sleep on an iron bed; if they
> did not fit the bed he cut off or stretched their limbs to make them fit
> [32].
>
> Richard Everson

Modern sensors have brought the classic Wahba problem [86], or slightly
differently the Procrustes analysis problem [36], into greater generality that
has increasing importance to computer vision [57, 39], computer graphics [65],
and robotics [14]. We formalize this generalization as follows.

######  Problem 1 (simultaneous rotation and correspondence search).

Consider point sets
$\mathcal{Q}=\\{\bm{q}_{1},\dots,\bm{q}_{m}\\}\subset\mathbb{R}^{3}$ and
$\mathcal{P}=\\{\bm{p}_{1},\dots,\bm{p}_{n}\\}\subset\mathbb{R}^{3}$ with
$m\geq n$. Let $\mathcal{C}^{*}$ be a subset of
$[m]\times[n]:=\\{1,\dots,m\\}\times\\{1,\dots,n\\}$ of size $k^{*}$, called
the inlier correspondence set, such that all pairs $(i_{1},j_{1})$ and
$(i_{2},j_{2})$ of $\mathcal{C}^{*}$ satisfy $i_{1}\neq i_{2}$ and $j_{1}\neq
j_{2}$. Assume that

$\displaystyle\bm{q}_{i}$
$\displaystyle=\bm{R}^{*}\bm{p}_{j}+\bm{\epsilon}_{i,j},\ \ \ \textnormal{if}\
(i,j)\in\mathcal{C}^{*}$ (1)

where $\bm{\epsilon}_{i,j}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$ is noise,
$\bm{R}^{*}$ is an unknown $3$D rotation, and $(\bm{q}_{i},\bm{p}_{j})$ is
called an inlier. If $(i,j)\notin\mathcal{C}^{*}$ then
$(\bm{q}_{i},\bm{p}_{j})$ is arbitrary and is called an outlier. The goal of
the simultaneous rotation and correspondence search problem is to
simultaneously estimate the $3$D rotation $\bm{R}^{*}$ and the inlier
correspondence set $\mathcal{C}^{*}$ from point sets $\mathcal{Q}$ and
$\mathcal{P}$.

We focus on Problem 1 for two reasons. First, it already encompasses several
vision applications such as image stitching [19]. Second, the more general and
more important simultaneous pose and correspondence problem, which involves an
extra unknown translation in (1), reduces to Problem 1 by eliminating the
translation parameters (at the cost of squaring the number of measurements)
[91]. As surveyed in [43], whether accurate and fast algorithms exist for
solving the pose and correspondence search is largely an open question.
Therefore, solving the simpler Problem 1 efficiently is an important step for
moving forward.

For Problem 1 or its variants, there is a vast literature of algorithms that
are based on i) local optimization via iterative closest points
($\operatorname{\texttt{ICP}}$) [15, 23, 76] or graduated non-convexity
($\operatorname{\texttt{GNC}}$) [96, 87, 2] or others [68, 46, 26], ii) global
optimization by branch $\&$ bound [57, 24, 92, 73, 21, 80, 61, 62], iii)
outlier removal techniques [19, 72, 71, 91, 79], iv) semidefinite programming
[65, 88, 44, 90, 81], v) RANSAC [33, 58, 82, 59], vi) deep learning [25, 42,
10, 5], and vii) spherical Fourier transform [14]. But all these methods, if
able to accurately solve Problem 1 with the number $k^{*}$ of inliers
extremely small, take $\Omega(mn)$ time. Yet we have:

######  Theorem 1 ($\operatorname{\texttt{ARCS}}$).

Suppose there are at least two inliers, $k^{*}\geq 2$, and that the point sets
$\mathcal{Q}$ and $\mathcal{P}$ of Problem 1 are noiseless “in general
position”. Then there is an algorithm that solves Problem 1 in $O(m\log m)$
time and $O(m)$ space.

######  Remark 1 (general position assumption).

In Theorem 1, by “in general position” we mean that i) for any outlier
$(\bm{q}_{i},\bm{p}_{j})$, we have
$\big{\|}\bm{q}_{i}\big{\|}_{2}\neq\big{\|}\bm{p}_{j}\big{\|}_{2}$, ii) there
exists some inlier pairs $(\bm{q}_{i_{1}},\bm{p}_{j_{1}})$ and
$(\bm{q}_{i_{2}},\bm{p}_{j_{2}})$ such that $\bm{q}_{i_{1}}$ and
$\bm{q}_{i_{2}}$ are not parallel. If point sets $\mathcal{Q}$ and
$\mathcal{P}$ are randomly sampled from $\mathbb{R}^{3}$, these two conditions
hold true with probability $1$.

A numerical illustration of Theorem 1 is that our
$\operatorname{\texttt{ARCS}}$ solver, to be described in §3, can handle the
case where $m=10^{6},n=8\times 10^{5}$ and $k^{*}=2$, in about $0.1$ seconds
(cf. Table 1).111We run experiments on an Intel(R) i7-1165G7, $16$GB laptop.
In the paper we consider random instead of adversarial outliers. However, like
other correspondence-based minimal solvers for geometric vision [34, 70, 51,
52, 53], $\operatorname{\texttt{ARCS}}$ might be fragile to noise. That being
said, it can be extended to the noisy case, leading to a three-step algorithm
called $\operatorname{\texttt{ARCS+}}$, which we summarize next.

The first step $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ of
$\operatorname{\texttt{ARCS+}}$ extends $\operatorname{\texttt{ARCS}}$ by
establishing correspondences under noise.
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ outputs in $O(\ell+m\log m)$ time
a candidate correspondence set $\overline{\mathcal{C}}$ of size $\ell$ that
contains $\mathcal{C}^{*}$. Problem 1 then reduces to estimating $\bm{R}^{*}$
and $\mathcal{C}^{*}$ from $\mathcal{P},\mathcal{Q},$ and hypothetical
correspondences $\overline{\mathcal{C}}$, a simpler task of robust rotation
search [96, 19, 72, 88].

The second step $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ of
$\operatorname{\texttt{ARCS+}}$ is to remove outliers from the previous step
$1$. To do so we approximately maximize an appropriate consensus over
$\operatorname{SO}(3)$ (§4.2). Instead of mining inliers in
$\operatorname{SO}(3)$ [57, 39, 11, 73, 48], we show that the parameter space
of consensus maximization can be reduced from $\operatorname{SO}(3)$ to
$\mathbb{S}^{2}$ and further to $[0,\pi]$ (see [19] for a different
reduction). With this reduction, $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$
removes outliers via repeatedly solving in $O(\ell\log\ell)$ time a
computational geometry problem, interval stabbing [28] (§4.2.1). Note that
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ only repeats for $s\approx 90$
times to reach satisfactory accuracy. Therefore, conceptually, for $\ell\geq
10^{6}$, it is $10^{4}$ times faster than the most related outlier removal
method $\operatorname{\texttt{GORE}}$ [19], which uses $O(\ell^{2}\log\ell)$
time (Table 4).

The third and final step $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ of our
$\operatorname{\texttt{ARCS+}}$ pipeline is to accurately estimate the
rotation, using the consensus set from the second step (§4.3). In short,
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ is a Riemannian subgradient
descent method. Our novelty here is to descend in the space $\mathbb{S}^{3}$
of unit quaternions, not $\operatorname{SO}(3)$ [16]. This allows us to
derive, based on [60], that $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$
converges linearly though locally to the ground-truth unit quaternion, thus
obtaining the first to our knowledge convergence rate guarantee for robust
rotation search.

Numerical highlights are in order (§5).
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ is an outlier pruning procedure
for robust rotation search that can handle extremely small inlier ratios
$k^{*}/\ell=3000/10^{7}=0.03\%$ in $5$ minutes;
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}+\operatorname{\texttt{ARCS+}}_{\texttt{R}}$,
or $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ for short, accurately solves
the robust rotation search problem with $k^{*}/\ell=10^{3}/10^{6}$ in $23$
seconds (see Table 4).
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}+\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$,
that is $\operatorname{\texttt{ARCS+}}$, solves Problem 1 with
$m=10^{4},n=8000,k^{*}=2000$ in $90$ seconds (see Figure 2). To the best of
our knowledge, all these challenging cases have not been considered in prior
works. In fact, as we will review soon (§2), applying state-of-the-art methods
to those cases either gives wrong estimates of rotations, or takes too much
time ($\geq 8$ hours), or exhausts the memory (Table 4).

##  2 Prior Art: Accuracy Versus Scalability

Early efforts on Problem 1 have encountered an accuracy versus scalability
dilemma. The now classic $\operatorname{\texttt{ICP}}$ algorithm [15]
estimates the rotation and correspondences in an alternating fashion, running
in real time but requiring a high-quality and typically unavailable
initialization to avoid local and usually poor minima; the same is true for
its successors [23, 76, 26, 68, 46]. The $\operatorname{\texttt{GO-ICP}}$
method [93, 92] of the branch $\&$ bound type enumerates initializations fed
to $\operatorname{\texttt{ICP}}$ to reach a global minimum—in exponential
time; the same running time bound is true for its successors [21, 73, 62].

The above $\operatorname{\texttt{ICP}}$ versus $\operatorname{\texttt{GO-
ICP}}$ dilemma was somewhat alleviated by a two-step procedure: i) compute a
candidate correspondence set $\hat{\mathcal{C}}$, via hand-crafted [77] or
learned [35] feature descriptors, and ii) estimate the rotation from point
sets indexed by $\hat{\mathcal{C}}$. But, as observed in [91], due to the
quality of the feature descriptors, there could be fewer than $2$ inliers
remaining in $\hat{\mathcal{C}}$, from which the ground-truth rotation can
never be determined. An alternative and more conservative idea is to use all-
to-all correspondences $\hat{\mathcal{C}}:=[m]\times[n]$, although now the
inlier ratio becomes extremely small.

This justifies why researchers have recently focused on designing robust
rotation search algorithms for extreme outlier rates, _e.g_., $\geq 90$
outliers out of $100$. One such design is $\operatorname{\texttt{GORE}}$ [19],
a guaranteed outlier removal algorithm of $O(\ell^{2}\log\ell)$ time
complexity that heavily exploits the geometry of $\operatorname{SO}(3)$. The
other one is the semidefinite relaxation $\operatorname{\texttt{QUASAR}}$ of
[88], which involves sophisticated manipulation on unit quaternions;
$\ell\approx 1000$ constitutes the current limit on the number of points this
relaxation can handle. Yet another one is $\operatorname{\texttt{TEASER++}}$
[91]; its robustness to outliers comes mainly from finding via parallel branch
$\&$ bound [75] a maximum clique of the graph whose vertices represent point
pairs and whose edges indicate whether two point pairs can simultaneously be
inliers. This maximum clique formulation was also explored by [71] where it
was solved via a different branch $\&$ bound algorithm. Since finding a
maximum clique is in general NP-hard, their algorithms take exponential time
in the worst case; in addition, $\operatorname{\texttt{TEASER++}}$ was
implemented to trade $O(\ell^{2})$ space for speed. One should also note
though that if noise is small then the graph is sparse so that the otherwise
intractable branch $\&$ bound algorithm can be efficient. Since constructing
such a graph entails checking $\binom{\ell}{2}$ point pairs, recent follow-up
works [82, 81, 58, 79, 63] that use such a graph entail $O(\ell^{2})$ time
complexity. While all these methods are more accurate than scalable, the
following two are on the other side. $\operatorname{\texttt{FGR}}$ [96]
combines graduated non-convexity ($\operatorname{\texttt{GNC}}$) and
alternating minimization, while $\operatorname{\texttt{GNC-TLS}}$ [87]
combines truncated least squares, iteratively reweighted least-squares, and
$\operatorname{\texttt{GNC}}$. Both of them scale gracefully with $\ell$,
while being robust against up to $80/100=80\%$ outliers.

Is such accuracy versus scalability dilemma of an inherent nature of the
problems here, or can we escape from it?

##  3 ARCS: Accuracy & Scalability

Basic Idea. Although perhaps not explicitly mentioned in the literature, it
should be known that there is a simple algorithm that solves Problem 1 under
the assumptions of Theorem 1. This algorithm first computes the $\ell_{2}$
norm of each point in $\mathcal{Q}$ and $\mathcal{P}$ and the difference
$d_{i,j}:=\big{\|}\bm{q}_{i}\big{\|}_{2}-\big{\|}\bm{p}_{j}\big{\|}_{2}$.
Since $\mathcal{Q}$ and $\mathcal{P}$ are in general position (Remark 1), we
have that $(\bm{q}_{i},\bm{p}_{j})$ is an inlier pair if and only if
$d_{i,j}=0$. Based on the $d_{i,j}$’s, extract all such inlier pairs. Since
$k^{*}\geq 2$, and by the general position assumption (Remark 1), there exist
two inlier pairs say $(\bm{q}_{1},\bm{p}_{1}),(\bm{q}_{2},\bm{p}_{2})$ such
that $\bm{q}_{1}$ and $\bm{q}_{2}$ are not parallel. As a result and as it has
been well-known since the $1980$’s [40, 41, 64, 4], if not even earlier [86,
78], $\bm{R}^{*}$ can be determined from the two inlier pairs by SVD.

ARCS: Efficient Implementation. Not all the $d_{i,j}$’s should be computed in
order to find the correspondence set $\mathcal{C}^{*}$, meaning that the
otherwise $O(mn)$ time complexity can be reduced. Our
$\operatorname{\texttt{ARCS}}$ Algorithm 1 seeks all point pairs
$(\bm{q}_{i},\bm{p}_{j})$’s whose norms are close, _i.e_., they satisfy
$|d_{i,j}|\leq c$, for some sufficiently small $c\geq 0$. Here $c$ is provided
as an input of $\operatorname{\texttt{ARCS}}$ and set as $0$ in the current
context. It is clear that, under the general position assumption of Theorem 1,
the set $\overline{\mathcal{C}}$ returned by $\operatorname{\texttt{ARCS}}$ is
exactly the ground-truth correspondence set $\mathcal{C}^{*}$. It is also
clear that $\operatorname{\texttt{ARCS}}$ takes $O(m\log m)$ time and $O(m)$
space (recall $m\geq n\geq|\mathcal{C}^{*}|$).

1Input:
$\mathcal{Q}=\\{\bm{q}_{i}\\}_{i=1}^{m},\mathcal{P}=\\{\bm{p}_{j}\\}_{j=1}^{n},c\geq
0$;

2Sort $\mathcal{Q}$ so that (w.l.o.g.)
$\big{\|}\bm{q}_{1}\big{\|}_{2}\leq\cdots\leq\big{\|}\bm{q}_{m}\big{\|}_{2}$;

3Sort $\mathcal{P}$ so that (w.l.o.g.)
$\big{\|}\bm{p}_{1}\big{\|}_{2}\leq\cdots\leq\big{\|}\bm{p}_{n}\big{\|}_{2}$;

4$i=1$; $j=1$; $\overline{\mathcal{C}}=\varnothing$;

5while  _$i\leq m$ and $j\leq n$_ do

6
$d_{i,j}\leftarrow\big{\|}\bm{q}_{i}\big{\|}_{2}-\big{\|}\bm{p}_{j}\big{\|}_{2}$;

7 if  _$d_{i,j} >c$_ then

8 $j\leftarrow j+1$;

9 end if

10 if  _$d_{i,j} <-c$_ then

11 $i\leftarrow i+1$;

12 end if

13

14 if  _$-c\leq d_{i,j}\leq c$_ then

15 $\overline{\mathcal{C}}\leftarrow\overline{\mathcal{C}}\cup(i,j)$;
$(i,j)\leftarrow(i+1,j+1)$;

16 end if

17

18 end while

return $\overline{\mathcal{C}}$;

Algorithm 1 $\operatorname{\texttt{ARCS}}$ Table 1: Time (msec) of generating noiseless Gaussian point sets (G) and solving Problem 1 by $\operatorname{\texttt{ARCS}}$ ($100$ trials, $k^{*}=2$).  $m$ | $10^{4}$ | $10^{5}$ | $10^{6}$
---|---|---|---
$n$ | $8\times 10^{3}$ | $8\times 10^{4}$ | $8\times 10^{5}$
G | $5.9$ | $15.0$ | $212.8$
Brute Force | $73.8$ | $8304$ | $8380441.5$
$\operatorname{\texttt{ARCS}}$ | $1.51$ | $8.4$ | $121.1$

We proved Theorem 1. It is operating in the noiseless case that allows us to
show that Problem 1 can be solved accurately and at large scale. Indeed,
$\operatorname{\texttt{ARCS}}$ can handle more than $10^{6}$ points with
$k^{*}=2$ in about $0.1$ seconds, even though generating those points has
taken more than $0.2$ seconds, as shown in Table 1.222For experiments in
Tables 1 and 2 we generate data as per Section 5.1. Note that in the setting
of Table 1 we have only $k^{*}=2$ overlapping points, a situation where all
prior methods mentioned in §1 and §2, if directly applicable, in principle
break down. One reason is that they are not designed to handle the noiseless
case. The other reason is that the overlapping ratio $k^{*}/m$ of Table 1 is
the minimum possible. While the achievement in Table 1 is currently limited to
the noiseless case, it forms a strong motivation that urges us to robustify
$\operatorname{\texttt{ARCS}}$ to noise, while keeping as much of its accuracy
and scalability as possible. Such robustification is the main theme of the
next section.

##  4 ARCS+: Robustifying ARCS to Noise

Here we consider Problem 1 with noise $\bm{\epsilon}_{i,j}$. We will
illustrate our algorithmic ideas by assuming
$\bm{\epsilon}_{i,j}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$, although this is
not necessary for actual implementation. As indicated in §1,
$\operatorname{\texttt{ARCS+}}$ has three steps. We introduce them
respectively in the next three subsections.

###  4.1 Step 1: Finding Correspondences Under Noise

A simple probability fact is
$\big{\|}\bm{q}_{i}-\bm{R}^{*}\bm{p}_{j}\big{\|}_{2}\leq 5.54\sigma$ for any
inlier $(\bm{q}_{i},\bm{p}_{j})$, so $|d_{i,j}|\leq 5.54\sigma$ with
probability at least $1-10^{-6}$ (see, _e.g_., [91]). To establish
correspondences under noise, we need to modify333The details of this
modification can be found at:
https://github.com/liangzu/ARCS/blob/main/ARCSplus_N.m the while loop of
Algorithm 1, such that, in $O(\ell+m\log m)$ time, it returns the set
$\overline{\mathcal{C}}$ of all correspondences of size $\ell$ where each
$(i,j)\in\overline{\mathcal{C}}$ satisfies $|d_{i,j}|\leq c$, with $c$ now set
to $5.54\sigma$. Note that, to store the output correspondences, we need an
extra $O(\ell)$ time, which can not be simply ignored as $\ell$ is in general
larger than $m$ in the presence of noise (Table 2). We call this modified
version $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$.
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ gives a set
$\overline{\mathcal{C}}$ that contains all inlier correspondences
$\mathcal{C}^{*}$ with probability at least $(1-10^{-6})^{k^{*}}$. This
probability is larger than $99.9\%$ if $k^{*}\leq 10^{3}$, or larger than
$99\%$ if $k^{*}\leq 10^{4}$.

######  Remark 2 (feature matching versus all-to-all correspondences versus
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$).

Feature matching methods provide fewer than $n$ hypothetical correspondences
and thus speed up the subsequent computation, but they might give no inliers.
Using all-to-all correspondences preserves all inliers, but a naive
computation needs $O(mn)$ time and leads to a large-scale problem with extreme
outlier rates. $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ strikes a balance
by delivering in $O(\ell+m\log m)$ time a candidate correspondence set
$\overline{\mathcal{C}}$ of size $\ell$ containing all inliers with high
probability and with $\ell\ll mn$.

For illustration, Table 2 reports the number $\ell$ of correspondences that
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ typically yields. As shown, even
though $\ell/(mn)$ is usually smaller than $5\%$, yet $\ell$ itself could be
very large, and the inlier ratio $k^{*}/\ell$ is extremely small (_e.g_.,
$\leq 0.05\%$). This is perhaps the best we could do for the current stage,
because for now we only considered every point pair individually, while any
pair $(\bm{q}_{i},\bm{p}_{i})$ is a potential inlier if it satisfies the
necessary (but no longer sufficient) condition $|d_{i,j}|\leq c$. On the other
hand, collectively analyzing the remaining point pairs allows to further
remove outliers, and this is the major task of our next stage (§4.2).

Table 2: The number $\ell$ of candidate correspondences produced by $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ on synthetic noisy Gaussian point sets. A single trial. $m$ | $1000$ | $5000$ | $10000$
---|---|---|---
$n$ | $800$ | $4000$ | $8000$
$k^{*}$ | $200$ | $1000$ | $2000$
$\ell$ | $36622$ | $931208$ | $3762888$
$\ell/(mn)$ | $4.58\%$ | $4.66\%$ | $4.70\%$

###  4.2 Step 2: Outlier Removal

Let there be some correspondences given, by, _e.g_., either
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ or feature matching (cf. Remark
2). Then we arrive at an important special case of Problem 1, called robust
rotation search. For convenience we formalize it below:

######  Problem 2.

(robust rotation search) Consider $\ell$ pairs of $3$D points
$\\{(\bm{y}_{i},\bm{x}_{i})\\}_{i=1}^{\ell}$, with each pair satisfying

$\displaystyle\bm{y}_{i}=\bm{R}^{*}\bm{x}_{i}+\bm{o}_{i}+\bm{\epsilon}_{i}.$
(2)

Here $\bm{\epsilon}_{i}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$ is noise,
$\bm{o}_{i}=\bm{0}$ if $i\in\mathcal{I}^{*}$ where
$\mathcal{I}^{*}\subset[\ell]$ is of size $k^{*}$, and if
$i\notin\mathcal{I}^{*}$ then $\bm{o}_{i}$ is nonzero and arbitrary. The task
is to find $\bm{R}^{*}$ and $\mathcal{I}^{*}$.

The percentage of outliers in Problem 2 can be quite large (cf. Table 2), so
our second step $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ here is to remove
outliers. In §4.2.1, we shortly review the interval stabbing problem, on which
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ of §4.2.2 is based.

####  4.2.1 Preliminaries: Interval Stabbing

Consider a collection of subsets of $\mathbb{R}$,
$\\{\mathcal{J}_{i}\\}_{i=1}^{L}$, where each $\mathcal{J}_{i}$ is an interval
of the form $[a,b]$. In the interval stabbing problem, one needs to determine
a point $\omega\in\mathbb{R}$ and a subset $\mathcal{I}$ of
$\\{\mathcal{J}_{i}\\}_{i=1}^{L}$, so that $\mathcal{I}$ is a maximal subset
whose intervals overlap at $\omega$. Formally, we need to solve

$\displaystyle\max_{\mathcal{I}\subset[L],\omega\in\mathbb{R}}\ \ \ \ \ \ \
|\mathcal{I}|$ (3) s.t. $\displaystyle\ \ \omega\in\mathcal{J}_{i},\ \ \forall
i\in\mathcal{I}$

For this purpose, the following result is known.

######  Lemma 1 (interval stabbing).

Problem (3) can be solved in $O(L\log L)$ time and $O(L)$ space.

Actually, the interval stabbing problem can be solved using sophisticated data
structures such as interval tree [28] or interval skip list [37]. On the other
hand, it is a basic exercise to find an algorithm that solves Problem (3),
which, though also in $O(L\log L)$ time, involves only a sorting operation and
a for loop (details are omitted, see, _e.g_., [20]). Finally, note that the
use of interval stabbing for robust rotation search is not novel, and can be
found in $\operatorname{\texttt{GORE}}$ [19, 72]. However, as the reader might
realize after §4.2.2, our use of interval stabbing is quite different from
$\operatorname{\texttt{GORE}}$.

####  4.2.2 The Outlier Removal Algorithm

We now consider the following consensus maximization:

$\begin{split}&\max_{\mathcal{I}\subset[\ell],\bm{R}\in\operatorname{SO}(3)}\
\ \ \ \ \ \ |\mathcal{I}|\\\ \textnormal{s.t.}\ \ \ \ &\ \
\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}\leq c,\ \ \forall
i\in\mathcal{I}.\end{split}$ (4)

It has been shown in [84] that for the very related robust fitting problem,
such consensus maximization is in general NP-hard444Interestingly, consensus
maximization over $\operatorname{SO}(2)$, _i.e_., the $2$D version of
(LABEL:eq:R-3D), can be solved in $O(\ell\log\ell)$ time; see [20].. Thus it
seems only prudent to switch our computational goal from solving
(LABEL:eq:R-3D) exactly to approximately.

From $\operatorname{SO}(3)$ to $\mathbb{S}^{2}$. Towards this goal, we first
shift our attention to $\mathbb{S}^{2}$ where the rotation axis $\bm{b}^{*}$
of $\bm{R}^{*}$ lives. An interesting observation is that the axis
$\bm{b}^{*}$ has the following interplay with data, independent of the
rotation angle of $\bm{R}^{*}$.

######  Proposition 1.

Let $\bm{v}_{i}:=\bm{y}_{i}-\bm{x}_{i}$. Recall
$\bm{\epsilon}_{i}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$. If
$(\bm{y}_{i},\bm{x}_{i})$ is an inlier pair, then
$\bm{v}_{i}^{\top}\bm{b}^{*}\sim\mathcal{N}(0,\sigma^{2})$, and so
$|\bm{v}_{i}^{\top}\bm{b}^{*}|\leq 4.9\sigma$ with probability at least
$1-10^{-6}$.

Proposition 1 (cf. Appendix C) leads us to Problem (5):

$\begin{split}&\max_{\mathcal{I}\subset[\ell],\bm{b}\in\mathbb{S}^{2}}\ \ \ \
\ \ \ |\mathcal{I}|\\\ \textnormal{s.t.}\ \ \ \ &\ \
|\bm{v}_{i}^{\top}\bm{b}|\leq\bar{c},\ \ \forall i\in\mathcal{I}\\\ &\ \ \ \ \
\ \ \ b_{2}\geq 0.\end{split}$ (5)

In (5) the constraint on the second entry $b_{2}$ of $\bm{b}$ is to eliminate
the symmetry, and Proposition 1 suggests to set $\bar{c}:=4.9\sigma$. Problem
(5) is easier than (LABEL:eq:R-3D) as it has fewer degrees of freedom; see
also [19] where a different reduction to a 2 DoF (sub-)problem was derived for
$\operatorname{\texttt{GORE}}$.

Solving (5) is expected to yield an accurate estimate of $\bm{b^{*}}$, from
which the rotation angle can later be estimated. Problem (5) reads: find a
plane (defined by the normal $\bm{b}$) that approximately contains as much
points $\bm{v}_{i}$’s as possible. This is an instance of the robust subspace
learning problem [85, 98, 97, 54, 30, 31, 94], for which various scalable
algorithms with strong theoretical guarantees have been developed in more
tractable formulations (_e.g_., $\ell_{1}$ minimization) than consensus
maximization. Most notably, the so-called dual principal component pursuit
formulation [85] was proved in [98] to be able to tolerate
$O\big{(}(k^{*})^{2}\big{)}$ outliers. Still, all these methods can not handle
as many outliers as we currently have (cf. Table 2), even though they can
often minimize their objective functions to global optimality.

From $\mathbb{S}^{2}$ to $[0,\pi]$. We can further “reduce” the degrees of
freedom in (5) by $1$, through the following lens. Certainly
$\bm{b}\in\mathbb{S}^{2}$ in (5) is determined by two angles
$\theta\in[0,\pi]$, $\phi\in[0,\pi]$. Now consider the following problem:

$\begin{split}&\max_{\mathcal{I}\subset[\ell],\theta\in[0,\pi]}\ \ \ \ \ \ \
|\mathcal{I}|\\\ \textnormal{s.t.}\ \ \ \ &\ \
|\bm{v}_{i}^{\top}\bm{b}|\leq\bar{c},\ \ \forall i\in\mathcal{I}\\\ &\
\bm{b}=[\sin(\theta)\cos(\phi),\ \sin(\theta)\sin(\phi),\
\cos(\theta)]^{\top}.\end{split}$ (6)

Problem (6) is a simplified version of (5) with $\phi$ given. Clearly, to
solve (5) it suffices to minimize the function $f:[0,\pi]\to\mathbb{R}$ which
maps any $\phi_{0}\in[0,\pi]$ to the objective value of (6) with
$\phi=\phi_{0}$. Moreover, we have:

######  Proposition 2.

Problem (6) can be solved in $O(\ell\log\ell)$ time and $O(\ell)$ space via
interval stabbing.

Proposition 2 gives an $O(\ell\log\ell)$ time oracle to access the values of
$f$. Since computing the objective value of (5) given $\theta,\phi$ already
needs $O(\ell)$ time, the extra cost of the logarithmic factor in Proposition
2 is nearly negligible. Since $f$ has only one degree of freedom, its global
minimizer can be found by one-dimensional branch $\&$ bound [47]. But this
entails exponential time complexity in the worst case, a situation we wish to
sidestep. Alternatively, the search space $[0,\pi]$ is now so small that the
following algorithm $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ turns out to
be surprisingly efficient and robust: i) sampling from $[0,\pi]$, ii) stabbing
in $\mathbb{S}^{2}$, and iii) stabbing in $\operatorname{SO}(3)$.

Sampling from $[0,\pi]$. Take $s$ equally spaced points
$\phi_{j}=(2j-1)\pi/(2s)$, $\forall j\in[s]$, on $[0,\pi]$. The reader may
find this choice of $\phi_{j}$’s similar to the uniform grid approach [69]; in
the latter Nesterov commented that “the reason why it works here is related to
the dimension of the problem”.

Stabbing in $\mathbb{S}^{2}$. For each $j\in[s]$, solve (6) with
$\phi=\phi_{j}$ to get $s$ candidate consensus set $\mathcal{I}_{j}$’s and $s$
angles $\theta_{j}$’s. From each $\phi_{j}$ and $\theta_{j}$ we obtain a
candidate rotation axis $\bm{b}_{j}$.

Stabbing in $\operatorname{SO}(3)$. Since now we have estimates of rotation
axes, $\bm{b}_{j}$’s, there is one degree of freedom remaining, the rotation
angle $\omega$. For this we consider:

$\begin{split}&\max_{\mathcal{I}\subset[\ell],\omega\in[0,2\pi]}\ \ \ \ \ \ \
|\mathcal{I}|\\\ \textnormal{s.t.}\ \ \ \ &\ \
\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}\leq c,\ \ \forall
i\in\mathcal{I}\\\ \bm{R}&=\bm{b}\bm{b}^{\top}+[\ \bm{b}\
]_{\times}\sin(\omega)+(\bm{I}_{3}-\bm{b}\bm{b}^{\top})\cos(\omega)\end{split}$
(7)

Here $[\ \bm{b}\ ]_{\times}\in\mathbb{R}^{3\times 3}$ denotes the matrix
generating the cross product $\times$ by $\bm{b}$, that is $[\ \bm{b}\
]_{\times}\bm{a}=\bm{b}\times\bm{a}$ for all $\bm{a}\in\mathbb{R}^{3}$.
Similarly to Proposition 2, we have the following result:

######  Proposition 3.

Problem (7) can be solved in $O(\ell\log\ell)$ time and $O(\ell)$ space via
interval stabbing.

After solving (7) with $\bm{b}=\bm{b}_{j}$ for each $j\in[s]$, we obtain $s$
candidate consensus sets
$\tilde{\mathcal{I}}_{1},\dots,\tilde{\mathcal{I}}_{s}$, and we choose the one
with maximal cardinality as an approximate solution to (LABEL:eq:R-3D).
Finally, notice that the time complexity $O(s\ell\log\ell)$ of
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ depends on the hyper-parameter
$s$. We set $s=90$ as an invariant choice, as suggested by Figure 1.

This output consensus set $\tilde{\mathcal{I}}$ typically has very few
outliers; see Table 3. Thus it will be used next in
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$, our final step for accurately
estimating the rotation (§4.3).

Table 3: The output of $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ with inputs from Table 2. Input Inlier Ratio | $\frac{200}{36622}$ | $\frac{1000}{931208}$ | $\frac{2000}{3762888}$
---|---|---|---
Output Inlier Ratio | $\frac{199}{213}$ | $\frac{993}{1314}$ | $\frac{1951}{3184}$

$1$$10$$10$$30$$50$$70$$90$$s$Rotation
Error$k^{*}/\ell=10^{1}/10^{2}$$k^{*}/\ell=10^{2}/10^{4}$$k^{*}/\ell=10^{3}/10^{6}$
(a)

$0.1$$1$$10$$10$$30$$50$$70$$90$$s$ (b)

Figure 1: Rotation errors (in degrees) of steps 2 and 3 for robust rotation
search methods with $s$ varying ($500$ trials, $\sigma=0.01$).

###  4.3 Step 3: Rotation Estimation

The final step $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ of
$\operatorname{\texttt{ARCS+}}$ is a refinement procedure that performs robust
rotation search on the output correspondences $\tilde{\mathcal{I}}$ of
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$. Since $\tilde{\mathcal{I}}$
contains much fewer outlier correspondences than we previously had (cf. Table
2 and 3), in what follows we simplify the notations by focusing on the point
set $\\{(\bm{y}_{i},\bm{x}_{i})\\}_{i\in[\ell]}$, which we assume has few
outliers (say $\leq 50\%$). Then, a natural formulation is

$\displaystyle\min_{\bm{R}\in\operatorname{SO}(3)}\sum_{i=1}^{\ell}\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}.$
(8)

Problem (8) appears easier to solve than consensus maximization
(LABEL:eq:R-3D), as it has a convex objective function at least. Next we
present the $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ algorithm and its
theory.

Algorithm. We start with the following equivalence.

######  Proposition 4.

We have
$\bm{w}^{\top}\bm{D}_{i}\bm{w}=\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}^{2}$,
where $\bm{w}\in\mathbb{S}^{3}$ is a quaternion representation of $\bm{R}$ of
$\eqref{eq:L1-SO3}$, and $\bm{D}_{i}\in\mathbb{R}^{4\times 4}$ is a positive
semi-definite matrix whose entries depend on $\bm{x}_{i}$, $\bm{y}_{i}$. So
Problem (8) is equivalent to

$\displaystyle\min_{\bm{w}\in\mathbb{S}^{3}}h(\bm{w}),\ \ \ \
h(\bm{w})=\sum_{i=1}^{\ell}\sqrt{\bm{w}^{\top}\bm{D}_{i}\bm{w}}.$ (9)

The exact relation between unit quaternions and rotations is reviewed in
Appendix A, where Proposition 4 is proved and the expression of $\bm{D}_{i}$
is given. For what follows, it suffices to know that a unit quaternion is
simply a unit vector of $\mathbb{R}^{4}$, and that the space of unit
quaternions is $\mathbb{S}^{3}$.

Note that the objective $h$ of (9) is convex, while both problems (8) and (9)
are nonconvex (due to the constraint) and nonsmooth (due to the objective).
Though (8) and (9) are equivalent, the advantage of (9) will manifest itself
soon. Before that, we first introduce the
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ algorithm for solving (9).
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ falls into the general Riemannian
subgradient descent framework (see, _e.g_., [60]). It is initialized at some
unit quaternion $\bm{w}^{(0)}\in\mathbb{S}^{3}$ and proceeds by

$\displaystyle\bm{w}^{(t+1)}\leftarrow\operatorname{Proj}_{\mathbb{S}^{3}}\big{(}\bm{w}^{(t)}-\gamma^{(t)}\operatorname{\tilde{\nabla}_{\textnormal{s}}}h(\bm{w}^{(t)})\big{)},$
(10)

where $\operatorname{Proj}_{\mathbb{S}^{3}}(\cdot)$ projects a vector onto
$\mathbb{S}^{3}$, $\gamma^{(t)}$ is some stepsize,
$\operatorname{\tilde{\nabla}_{\textnormal{s}}}h(\bm{w}^{(t)})$ is a
Riemannian subgradient555We follow [60] where a Riemannian subgradient
$\operatorname{\tilde{\nabla}_{\textnormal{s}}}h(\bm{w})$ at
$\bm{w}\in\mathbb{S}^{3}$ is defined as the projection of some subgradient
$\operatorname{\nabla_{\textnormal{s}}}h(\bm{w})$ of $h$ at $\bm{w}$ onto the
tangent space of $\mathbb{S}^{3}$ at $\bm{w}$, _i.e_.,
$\operatorname{\tilde{\nabla}_{\textnormal{s}}}h(\bm{w}):=(\bm{I}_{4}-\bm{w}\bm{w}^{\top})\operatorname{\nabla_{\textnormal{s}}}h(\bm{w})$.
See [12] for how to compute a subgradient of some given function. of $h$ at
$\bm{w}^{(t)}$.

Theory. Now we are able to compare (8) and (9) from a theoretical perspective.
As proved in [16], for any fixed outlier ratio and $k^{*}>0$, Riemannian
subgradient descent when applied to (8) with proper initialization converges
to $\bm{R}^{*}$ in finite time, as long as i) $\ell$ is sufficiently large,
ii) all points $\bm{y}_{i}$’s and $\bm{x}_{i}$’s are uniformly distributed on
$\mathbb{S}^{2}$, iii) there is no noise. But in [16] no convergence rate is
given. One main challenge of establishing convergence rates there is that
projecting on $\operatorname{SO}(3)$ does not enjoy a certain kind of
nonexpansiveness property, which is important for convergence analysis (cf.
Lemma 1 of [60]). On the other hand, projection onto $\mathbb{S}^{3}$ of (9)
does satisfy such property. As a result, we are able to provide convergence
rate guarantees for $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$. For example,
it follows directly from Theorem 2 of [60] that
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ (10) converges to an
$\varepsilon$-stationary point in $O(\varepsilon^{-4})$ iterations, even if
initialized arbitrarily.

We next give conditions for $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ to
converge linearly to the ground-truth unit quaternion $\pm\bm{w}^{*}$ that
represents $\bm{R}^{*}$. Let the distance between a unit quaternion $\bm{w}$
and $\pm\bm{w}^{*}$ be

$\displaystyle\operatorname{dist}(\bm{w},\pm\bm{w}^{*}):=\min\big{\\{}\big{\|}\bm{w}-\bm{w}^{*}\big{\|}_{2},\big{\|}\bm{w}+\bm{w}^{*}\big{\|}_{2}\big{\\}}.$

If $\operatorname{dist}(\bm{w},\pm\bm{w}^{*})<\rho$ with $\rho>0$ then
$\bm{w}$ is called $\rho$-close to $\pm\bm{w}^{*}$. We need the following
notion of sharpness.

######  Definition 1 (sharpness [18, 56, 50, 60]).

We say that $\pm\bm{w}^{*}$ is an $\alpha$-sharp minimum of (9) if $\alpha>0$
and if there exists a number $\rho_{\alpha}>0$ such that any unit quaternion
$\bm{w}\in\mathbb{S}^{3}$ that is $\rho_{\alpha}$-close to $\pm\bm{w}^{*}$
satisfies the inequality

$\displaystyle
h(\bm{w})-h(\bm{w}^{*})\geq\alpha\operatorname{dist}(\bm{w},\pm\bm{w}^{*}).$
(11)

We provide a condition below for $\pm\bm{w}^{*}$ to be $\alpha^{*}$-sharp:

######  Proposition 5.

If $\alpha^{*}:=k^{*}\eta_{\min}/\sqrt{2}-(\ell-k^{*})\eta_{\max}>0$ and if
$\bm{\epsilon}_{i}=\bm{0}$ in Problem 2, then Problem (9) admits $\pm w^{*}$
as an $\alpha^{*}$-sharp minimum. Here $\eta_{\min}$, $\eta_{\max}$ are
respectively

$\displaystyle\eta_{\min}$
$\displaystyle:=\frac{1}{k^{*}}\min_{\bm{w}\in\mathcal{S}^{*}\cap\mathbb{S}^{3}}\sum_{i\in\mathcal{I}^{*}}\sqrt{\bm{w}^{\top}\bm{D}_{i}\bm{w}},\text{\
\ and }$ (12) $\displaystyle\eta_{\max}$
$\displaystyle:=\frac{1}{\ell-k^{*}}\max_{\bm{w}\in\mathbb{S}^{3}}\sum_{i\in[\ell]\backslash\mathcal{I}^{*}}\sqrt{\bm{w}^{\top}\bm{D}_{i}\bm{w}},$
(13)

where $\mathcal{S}^{*}$ is the hyperplane of $\mathbb{R}^{4}$ perpendicular to
$\pm\bm{w}^{*}$.

Proposition 5 is proved in Appendix B.1. The condition $\alpha^{*}>0$ defines
a relation between the number of inliers $(k^{*})$ and outliers
$(\ell-k^{*})$, and involves two quantities $\eta_{\min}$ and $\eta_{\max}$
whose values depend on how $\bm{D}_{i}$’s are distributed on the positive
semi-definite cone. We offer probabilistic interpretations for $\eta_{\min}$
and $\eta_{\max}$ in Appendix B.2.

Table 4: Average errors in degrees $|$ standard deviation $|$ running times in seconds of various algorithms on synthetic data ($20$ trials).  Inlier Ratio $\frac{k^{*}}{\ell}$ | $\frac{10^{3}}{10^{5}}=1\%$ | $\frac{10^{3}}{10^{6}}=0.1\%$ | $\frac{3\times 10^{3}}{5\times 10^{6}}=0.06\%$ | $\frac{3\times 10^{3}}{10^{7}}=0.03\%$ | $\frac{10^{3}}{10^{7}}=0.01\%$
---|---|---|---|---|---
$\operatorname{\texttt{TEASER++}}$ [91] | out-of-memory |  |  |  |
$\operatorname{\texttt{RANSAC}}$ |  $0.39$ $|$ $0.20$ $|$ $29.1$ |  $\geq 8.4$ hours |  |  |
$\operatorname{\texttt{GORE}}$ [19, 72] |  $3.43$ $|$ $2.10$ $|$ $1698$ |  $\geq 12$ hours |  |  |
$\operatorname{\texttt{FGR}}$ [96] |  $52.2$ $|$ $68.5$ $|$ $3.64$ |  $95.0$ $|$ $60.9$ $|$ $37.7$ |  $84.9$ $|$ $59.4$ $|$ $145$ |  $86.5$ $|$ $56.9$ $|$ $311$ |  $97.3$ $|$ $61.3$ $|$ $314$
$\operatorname{\texttt{GNC-TLS}}$ [87] |  $3.86$ $|$ $9.51$ $|$ $0.13$ |  $63.4$ $|$ $50.5$ $|$ $2.26$ |  $49.9$ $|$ $31.1$ $|$ $15.9$ |  $90.2$ $|$ $45.6$ $|$ $40.1$ |  $120$ $|$ $34.3$ $|$ $36.3$
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ |  $9.92$ $|$ $13.1$ $|$ $0.12$ |  $65.2$ $|$ $48.9$ $|$ $0.96$ |  $55.6$ $|$ $38.3$ $|$ $5.58$ |  $88.4$ $|$ $36.2$ $|$ $12.6$ |  $98.2$ $|$ $36.0$ $|$ $12.2$
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ |  $0.86$ $|$ $0.29$ $|$ $1.71$ |  $0.99$ $|$ $0.37$ $|$ $23.2$ |  $0.91$ $|$ $0.30$ $|$ $125$ |  $0.98$ $|$ $0.42$ $|$ $287$ |  $55.6$ $|$ $60.9$ $|$ $281$
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ |  $0.03$ $|$ $0.03$ $|$ $1.72$ |  $0.09$ $|$ $0.07$ $|$ $23.2$ |  $0.11$ $|$ $0.07$ $|$ $125$ |  $0.22$ $|$ $0.15$ $|$ $287$ |  $55.4$ $|$ $60.1$ $|$ $281$

With Theorem 4 of [60] and Proposition 5 we have that
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ (10), if initialized properly and
with suitable stepsizes, converges linearly to the ground-truth unit
quaternion $\pm\bm{w}^{*}$, as long as $\pm\bm{w}^{*}$ is $\alpha^{*}$-sharp.
A formal statement is:

######  Theorem 2.

Suppose $\alpha^{*}:=k^{*}\eta_{\min}/\sqrt{2}-(\ell-k^{*})\eta_{\max}>0$. Let
$L_{h}$ be a Lipschitz constant of $h$. Run Riemannian subgradient descent
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ (10) with initialization
$\bm{w}^{(0)}$ satisfying
$\operatorname{dist}(\bm{w}^{(0)},\pm\bm{w}^{*})\leq\min\\{\alpha^{*}/L_{h},\rho_{\alpha^{*}}\\}$
and with geometrically diminishing stepsizes
$\gamma^{(t)}=\beta^{t}\gamma^{(0)}$, where

$\displaystyle\gamma^{(0)}$
$\displaystyle<\min\Bigg{\\{}\frac{2e_{0}(\alpha^{*}-L_{h}e_{0})}{L_{h}^{2}},\frac{e_{0}}{2(\alpha^{*}-L_{h}e_{0})}\Bigg{\\}},$
$\displaystyle\beta^{2}$
$\displaystyle\in\Bigg{[}1+2\Big{(}L_{h}-\frac{\alpha^{*}}{e_{0}}\Big{)}\gamma^{(0)}+\frac{L_{h}^{2}(\gamma^{(0)})^{2}}{e_{0}^{2}},\
1\Bigg{)},$ $\displaystyle e_{0}$
$\displaystyle=\min\Bigg{\\{}\max\Big{\\{}\operatorname{dist}(\bm{w}^{(0)},\pm\bm{w}^{*}),\frac{\alpha^{*}}{2L_{h}}\Big{\\}},\rho_{\alpha^{*}}\Bigg{\\}}.$

In the noiseless case ($\bm{\epsilon}_{i}=0$) we have each $\bm{w}^{(t)}$
satisfying

$\displaystyle\operatorname{dist}(\bm{w}^{(t)},\pm\bm{w}^{*})\leq\beta^{t}e_{0}.$
(14)

######  Remark 3 (a posteriori optimality guarantees).

Theorem 2 endows $\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ (10) with
convergence guarantee. On the other hand, a posteriori optimality guarantees
can be obtained via semidefinite certification [6, 22, 89, 91].

##  5 Experiments

In this section we evaluate $\operatorname{\texttt{ARCS+}}$ via synthetic and
real experiments for Problem 1, simultaneous rotation and correspondence
search. We also evaluate its components, namely
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ (§4.2) and
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ (§4.3) for Problem 2, robust
rotation search, as it is a task of independent interest. For both of the two
problems we compare the following state-of-the-art methods (reviewed in §2):
$\operatorname{\texttt{FGR}}$ [96], $\operatorname{\texttt{GORE}}$ [19],
$\operatorname{\texttt{RANSAC}}$, $\operatorname{\texttt{GNC-TLS}}$ [87], and
$\operatorname{\texttt{TEASER++}}$ [91].

###  5.1 Experiments on Synthetic Point Clouds

Setup. We set $\sigma=0.01$, $\bar{c}=c=5.54\sigma$, $n=\lfloor 0.8m\rfloor$,
and $s=90$ unless otherwise specified. For all other methods we used default
or otherwise appropriate parameters. We implemented
$\operatorname{\texttt{ARCS+}}$ in MATLAB. No parallelization was explicitly
used and no special care was taken for speed.

Robust Rotation Search. From $\mathcal{N}(0,\bm{I}_{3})$ we randomly sampled
point pairs $\\{(\bm{y}_{i},\bm{x}_{i})\\}_{i=1}^{\ell}$ with $k^{*}$ inliers
and noise $\bm{\epsilon}_{i}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$.
Specifically, we generated the ground-truth rotation $\bm{R}^{*}$ from an axis
randomly sampled from $\mathbb{S}^{2}$ and an angle from $[0,2\pi]$, rotated
$k^{*}$ points randomly sampled from $\mathcal{N}(0,\bm{I}_{3})$ by
$\bm{R}^{*}$, and added noise to obtain $k^{*}$ inlier pairs. Every outlier
point $\bm{y}_{j}$ or $\bm{x}_{j}$ was randomly sampled from
$\mathcal{N}(0,\bm{I}_{3})$ with the constraint
$-c\leq\big{\|}\bm{y}_{j}\big{\|}_{2}-\big{\|}\bm{x}_{j}\big{\|}_{2}\leq c$;
otherwise $(\bm{y}_{j},\bm{x}_{j})$ might simply be detected and removed by
computing $\big{\|}\bm{y}_{j}\big{\|}_{2}-\big{\|}\bm{x}_{j}\big{\|}_{2}$.

We compared $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ and
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ and their combination
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ with prior works. The results
are in Table 4. We first numerically illustrate the accuracy versus
scalability dilemma in prior works (§2). On the one hand, we observed an
extreme where accuracy overcomes scalability: $\operatorname{\texttt{RANSAC}}$
performed well with error $0.39$ when $k^{*}/\ell=10^{3}/10^{5}$, but its
running time increased greatly with decreasing inlier ratio, from $29$ seconds
to more than $8.4$ hours. The other extreme is where scalability overcomes
accuracy: Both $\operatorname{\texttt{GNC-TLS}}$ and
$\operatorname{\texttt{FGR}}$ failed in presence of such many outliers—as
expected—even though their running time scales linearly with $\ell$.

Table 4 also depicted the performance of our proposals
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ and
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$. Our approximate consensus
strategy $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ reached a balance
between accuracy and scalability. In terms of accuracy, it made errors smaller
than $1$ degree, as long as there are more than $3\times 10^{3}/10^{7}=0.03\%$
inliers; this was further refined by Riemannian subgradient descent
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$, so that their combination
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ had even lower errors. In terms
of scalability, we observed that $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$
is uniformly faster than $\operatorname{\texttt{FGR}}$, and is at least $1800$
times faster than $\operatorname{\texttt{GORE}}$ for
$k^{*}/\ell=10^{3}/10^{6}=0.1\%$. But it had been harder to measure exactly
how faster $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ is than
$\operatorname{\texttt{GORE}}$ and $\operatorname{\texttt{RANSAC}}$ for even
larger point sets. Finally, $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$
failed at $k^{*}/\ell=10^{3}/10^{7}=0.01\%$.

Simultaneous Rotation and Correspondence Search. We randomly sampled point
sets $\mathcal{Q}$ and $\mathcal{P}$ from $\mathcal{N}(0,\bm{I}_{3})$ with
$k^{*}$ inlier pairs and noise
$\bm{\epsilon}_{i,j}\sim\mathcal{N}(0,\sigma^{2}\bm{I}_{3})$ (cf. Problem 1).
Each outlier point was randomly and independently drawn also from
$\mathcal{N}(0,\bm{I}_{3})$. Figure 2 shows that
$\operatorname{\texttt{ARCS+}}$ accurately estimated the rotations for
$k^{*}\geq 2000$ (in $90$ seconds), and broke down at $k^{*}=1000$, a
situation where there were $k^{*}/m=10\%$ overlapping points. We did not
compare methods like $\operatorname{\texttt{TEASER++}}$,
$\operatorname{\texttt{GORE}}$, $\operatorname{\texttt{RANSAC}}$ here, because
giving them correspondences from $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$
would result unsatisfactory running time or accuracy (recall Tables 2 and 4),
while feature matching methods like $\operatorname{\texttt{FPFH}}$ do not
perform well on random synthetic data.

051015$7000$$5000$$3000$$1000$# Inlier Pairs ($k^{*}$)Rotation Error Figure 2: Rotation errors of $\operatorname{\texttt{ARCS+}}$ on synthetic Gaussian point clouds. $20$ trials, $m=10^{4}$, $n=\lfloor 0.8m\rfloor$, $\sigma=0.01$.  Table 5: Success rates of methods run on the scene pairs of the 3DMatch dataset [95] for which the ground-truth rotation and translation are provided (rotation error smaller than $10$ degree means a success [91]; see also the first paragraph of Appendix E).  Scene Type | Kitchen | Home 1 | Home 2 | Hotel 1 | Hotel 2 | Hotel 3 | Study Room | MIT Lab | Overall
---|---|---|---|---|---|---|---|---|---
# Scene Pairs | $506$ | $156$ | $208$ | $226$ | $104$ | $54$ | $292$ | $77$ | $1623$
$\operatorname{\texttt{TEASER++}}$ | 99.0% | 98.1% | $94.7\%$ | $98.7\%$ | 99.0% | $98.1\%$ | $97.0\%$ | $94.8\%$ | $97.72\%$
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ | $98.4\%$ | $97.4\%$ | 95.7% | $98.7\%$ | $98.1\%$ | 100% | 97.3% | 96.1% | $97.72\%$

###  5.2 Experiments on 3DMatch

The 3DMatch666License info: https://3dmatch.cs.princeton.edu/ dataset [95]
contains more than $1000$ point clouds for testing, representing $8$ different
scenes (such as kitchen, hotel, etc.), while the number of point clouds for
each scene ranges from $77$ to $506$. Each point cloud has more than $10^{5}$
points, yet in [95] there are $5000$ keypoints for each cloud. We used the
pretrained model777https://github.com/zgojcic/3DSmoothNet of the 3DSmoothNet
[35] to extract descriptors from these key points, and matched them using the
Matlab function pcmatchfeatures, with its parameter MatchThreshold set to the
maximum $1$. We assume that the ground-truth translation $\bm{t}^{*}$ is
given, and run $\operatorname{\texttt{TEASER++}}$ and
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ on
$(\bm{y}_{i}-\bm{t}^{*},\bm{x}_{i})$’s; the performance is comparable. We did
not compare other methods here, as $\operatorname{\texttt{TEASER++}}$
currently has the best performance on 3DMatch (to the best of our knowledge);
see [91] for comparison with optimization-based methods, and see [25] for the
success rates (recall) of deep learning methods.

See the supplementary materials for more experiments.

##  6 Discussion and Future Work

Despite of the progress that we made for robust rotation search and
simultaneous rotation and correspondence search on large-scale point clouds,
our $\operatorname{\texttt{ARCS+}}$ pipeline has a few limitations, and we
discuss them next.

For small datasets (_e.g_., $\ell\leq 500$), as in homography fitting [19],
other methods, _e.g_., MAGSAC++ [7, 8, 9], VSAC [45],
$\operatorname{\texttt{TEASER++}}$ [91], and $\operatorname{\texttt{GORE}}$
[19] might be considered with higher priority; they come with efficient C++
implementations. For more points, _e.g_., $\ell\geq 10^{4}$, but with higher
inlier rates than in Table 4 (_e.g_., $\geq 15\%$),
$\operatorname{\texttt{GNC-TLS}}$ [87] and $\operatorname{\texttt{RANSAC}}$
are our recommendations for what to use.

Modern point clouds have more than $10^{5}$ points, and are naturally
correspondences-less (cf. [20]). $\operatorname{\texttt{ARCS}}$ operates at
that scale in the absence of noise (Table 1), while
$\operatorname{\texttt{ARCS+}}$ can handle $m,n\approx 10000$ (Figure 2) and
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ can handle $\ell\approx 10^{7}$
correspondences (Table 4); all these are limited to the rotation-only case. To
find rotation (and translation) from such point sets “in the wild”, it seems
inevitable to downsample them. An interesting future work is to theoretically
quantify the tradeoff between downsampling factors and the final registration
performance. Another tradeoff to quantify, as implied by Remark 2, is this:
Can we design a correspondence matching algorithm that better balances the
number of remaining points and the number of remaining inliers? In particular,
such matching should take specific pose into consideration (cf.
$\operatorname{\texttt{ARCS}}$); many methods did not.

Like $\operatorname{\texttt{TEASER++}}$, $\operatorname{\texttt{GORE}}$,
$\operatorname{\texttt{GNC-TLS}}$, $\operatorname{\texttt{RANSAC}}$, our
algorithm relies on an inlier threshold $c$. While how to set this hyper-
parameter suitably is known for Gaussian noise with given variance, in
practice the distance threshold is usually chosen empirically, as Hartley &
Zisserman wrote [38]. While mis-specification of $c$ could fail the
registration, certain heuristics have been developed to alleviate the
sensitivity to such mis-specification; see [7, 8, 9, 3]. Finally, our
experience is to set $c$ based on the scale of the point clouds.

Our outlier removal component $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$
presented good performance (Table 3), yet with no optimality guarantees. Note
that, with $s=90$ we have $|\phi_{j}-\phi^{*}|\leq 1$ for some $\phi_{j}$,
while Figure 1(a) shows that $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ gave
roughly $1$ degree error at $s=90$. Theoretically justifying this is left as
future work. Without guarantees, registration could fail, which might lead to
undesired consequences in safety-critical applications. On the other hand, we
believe that $\operatorname{\texttt{ARCS+}}$ is a good demonstration of
trading optimality guarantees for accuracy and scalability; enforcing all of
the three properties amounts to requiring solving NP hard problems efficiently
at large scale! In fact, since any solutions might get certified for
optimality (Remark 3), bold algorithmic design ideas can be taken towards
improving accuracy and scalability, while relying on other tools for
optimality certification.

Acknowledgments. The first author was supported by the MINDS PhD fellowship at
Johns Hopkins University. This work was supported by NSF Grants 1704458 and
1934979, and by the Northrop Grumman Mission Systems Research in Applications
for Learning Machines (REALM) initiative.

## References

  * [1] Simon L Altmann.  Hamilton, Rodrigues, and the quaternion scandal.  Mathematics Magazine, 62(5):291–308, 1989.
  * [2] Pasquale Antonante, Vasileios Tzoumas, Heng Yang, and Luca Carlone.  Outlier-robust estimation: Hardness, minimally-tuned algorithms, and applications.  Technical report, arXiv:2007.15109v2 [cs.CV], 2020.
  * [3] Pasquale Antonante, Vasileios Tzoumas, Heng Yang, and Luca Carlone.  Outlier-robust estimation: Hardness, minimally tuned algorithms, and applications.  IEEE Transactions on Robotics, 2021.
  * [4] K Somani Arun, Thomas S Huang, and Steven D Blostein.  Least-squares fitting of two 3D point sets.  IEEE Transactions on Pattern Analysis and Machine Intelligence, (5):698–700, 1987.
  * [5] Xuyang Bai, Zixin Luo, Lei Zhou, Hongkai Chen, Lei Li, Zeyu Hu, Hongbo Fu, and Chiew-Lan Tai.  Pointdsc: Robust point cloud registration using deep spatial consistency.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 15859–15869, 2021.
  * [6] Afonso S Bandeira.  A note on probably certifiably correct algorithms.  Comptes Rendus Mathematique, 354(3):329–333, 2016.
  * [7] Daniel Barath, Jiri Matas, and Jana Noskova.  Magsac: marginalizing sample consensus.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10197–10205, 2019.
  * [8] Daniel Barath, Jana Noskova, Maksym Ivashechkin, and Jiri Matas.  Magsac++, a fast, reliable and accurate robust estimator.  In IEEE/CVF conference on computer vision and pattern recognition, pages 1304–1312, 2020.
  * [9] Daniel Barath, Jana Noskova, and Jiri Matas.  Marginalizing sample consensus.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [10] Dominik Bauer, Timothy Patten, and Markus Vincze.  Reagent: Point cloud registration using imitation and reinforcement learning.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 14586–14594, 2021.
  * [11] Jean-Charles Bazin, Yongduek Seo, and Marc Pollefeys.  Globally optimal consensus set maximization through rotation search.  In Asian Conference on Computer Vision, pages 539–551, 2012.
  * [12] Amir Beck.  First-Order Methods in Optimization.  Society for Industrial and Applied Mathematics, 2017.
  * [13] Amir Beck, Petre Stoica, and Jian Li.  Exact and approximate solutions of source localization problems.  IEEE Transactions on signal processing, 56(5):1770–1778, 2008.
  * [14] Lukas Bernreiter, Lionel Ott, Juan Nieto, Roland Siegwart, and Cesar Cadena.  PHASER: A robust and correspondence-free global pointcloud registration.  IEEE Robotics and Automation Letters, 6(2):855–862, 2021.
  * [15] PJ Besl and Neil D McKay.  A method for registration of 3D shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(2):239–256, 1992.
  * [16] Cindy Orozco Bohorquez, Yuehaw Khoo, and Lexing Ying.  Maximizing robustness of point-set registration by leveraging non-convexity.  Technical report, arXiv:2004.08772v3 [math.OC], 2020.
  * [17] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration inequalities: A nonasymptotic theory of independence.  Oxford university press, 2013.
  * [18] James V Burke and Michael C Ferris.  Weak sharp minima in mathematical programming.  SIAM Journal on Control and Optimization, 31(5):1340–1359, 1993\.
  * [19] Álvaro Parra Bustos and Tat-Jun Chin.  Guaranteed outlier removal for rotation search.  In IEEE International Conference on Computer Vision, pages 2165–2173, 2015.
  * [20] Zhipeng Cai, Tat-Jun Chin, Alvaro Parra Bustos, and Konrad Schindler.  Practical optimal registration of terrestrial lidar scan pairs.  ISPRS Journal of Photogrammetry and Remote Sensing, 147:118–131, 2019.
  * [21] Dylan Campbell and Lars Petersson.  GOGMA: Globally-optimal gaussian mixture alignment.  In IEEE Conference on Computer Vision and Pattern Recognition, 2016\.
  * [22] Luca Carlone, Giuseppe C. Calafiore, Carlo Tommolillo, and Frank Dellaert.  Planar pose graph optimization: Duality, optimal solutions, and verification.  IEEE Transactions on Robotics, 32(3):545–565, 2016.
  * [23] Dmitry Chetverikov, Dmitry Svirko, Dmitry Stepanov, and Pavel Krsek.  The trimmed iterative closest point algorithm.  In Object recognition supported by user interaction for service robots, volume 3, pages 545–548, 2002.
  * [24] Tat-Jun Chin, Yang Heng Kee, Anders Eriksson, and Frank Neumann.  Guaranteed outlier removal with mixed integer linear programs.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 5858–5866, 2016.
  * [25] Christopher Choy, Wei Dong, and Vladlen Koltun.  Deep global registration.  In IEEE Conference on Computer Vision and Pattern Recognition, 2020\.
  * [26] Haili Chui and Anand Rangarajan.  A new point matching algorithm for non-rigid registration.  Computer Vision and Image Understanding, 89(2-3):114–141, 2003\.
  * [27] Brian Curless and Marc Levoy.  A volumetric method for building complex models from range images.  In Annual Conference on Computer Graphics and Interactive Techniques, pages 303–312, 1996.
  * [28] Mark De Berg, Marc Van Kreveld, Mark Overmars, and Otfried Schwarzkopf.  Computational Geometry.  Springer, 1997.
  * [29] Tianjiao Ding, Yunchen Yang, Zhihui Zhu, Daniel P Robinson, René Vidal, Laurent Kneip, and Manolis C Tsakiris.  Robust homography estimation via dual principal component pursuit.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 6080–6089, 2020.
  * [30] Tianyu Ding, Zhihui Zhu, Tianjiao Ding, Yunchen Yang, René Vidal, Manolis C. Tsakiris, and Daniel Robinson.  Noisy dual principal component pursuit.  In International Conference on Machine Learning, pages 1617–1625, 2019.
  * [31] Tianyu Ding, Zhihui Zhu, René Vidal, and Daniel P Robinson.  Dual principal component pursuit for robust subspace learning: Theory and algorithms for a holistic approach.  In International Conference on Machine Learning, pages 2739–2748, 2021.
  * [32] Richard Everson.  Orthogonal, but not orthonormal, Procrustes problems.  Advances in Computational Mathematics, 3(4), 1998.
  * [33] Martin A Fischler and Robert C Bolles.  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM, 24(6):381–395, 1981.
  * [34] Xiao-Shan Gao, Xiao-Rong Hou, Jianliang Tang, and Hang-Fei Cheng.  Complete solution classification for the perspective-three-point problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(8):930–943, 2003.
  * [35] Zan Gojcic, Caifa Zhou, Jan D. Wegner, and Andreas Wieser.  The perfect match: 3D point cloud matching with smoothed densities.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 5540–5549, 2019.
  * [36] John C Gower, Garmt B Dijksterhuis, et al.  Procrustes Problems, volume 30.  Oxford University Press on Demand, 2004.
  * [37] Eric N. Hanson.  The interval skip list: A data structure for finding all intervals that overlap a point.  In Algorithms and Data Structures, pages 153–164, 1991.
  * [38] Richard Hartley and Andrew Zisserman.  Multiple View Geometry in Computer Vision.  Cambridge University Press, 2004.
  * [39] Richard I Hartley and Fredrik Kahl.  Global optimization through rotation space search.  International Journal of Computer Vision, 82(1):64–79, 2009.
  * [40] Berthold KP Horn.  Closed-form solution of absolute orientation using unit quaternions.  Journal of the Optical Society of America A, 4(4):629–642, 1987\.
  * [41] Berthold KP Horn, Hugh M Hilden, and Shahriar Negahdaripour.  Closed-form solution of absolute orientation using orthonormal matrices.  Journal of the Optical Society of America A, 5(7):1127–1135, 1988\.
  * [42] Shengyu Huang, Zan Gojcic, Mikhail Usvyatsov, Andreas Wieser, and Konrad Schindler.  PREDATOR: Registration of 3D point clouds with low overlap.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 4267–4276, 2021.
  * [43] Xiaoshui Huang, Guofeng Mei, Jian Zhang, and Rana Abbas.  A comprehensive survey on point cloud registration.  Technical report, arXiv:2103.02690v2 [cs.CV], 2021.
  * [44] Jose Pedro Iglesias, Carl Olsson, and Fredrik Kahl.  Global optimality for point set registration using semidefinite programming.  In IEEE Conference on Computer Vision and Pattern Recognition, 2020\.
  * [45] Maksym Ivashechkin, Daniel Barath, and Jiří Matas.  Vsac: Efficient and accurate estimator for h and f.  In IEEE/CVF International Conference on Computer Vision, pages 15243–15252, 2021.
  * [46] Bing Jian and Baba C. Vemuri.  Robust point set registration using gaussian mixture models.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(8):1633–1645, 2011.
  * [47] Yanmei Jiao, Yue Wang, Bo Fu, Qimeng Tan, Lei Chen, Minhang Wang, Shoudong Huang, and Rong Xiong.  Globally optimal consensus maximization for robust visual inertial localization in point and line map.  In International Conference on Intelligent Robots and Systems, pages 4631–4638, 2020.
  * [48] Kyungdon Joo, Hongdong Li, Tae-Hyun Oh, and In So Kweon.  Robust and efficient estimation of relative pose for cameras on selfie sticks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [49] Sham Kakade.  Symmetrization and Rademacher averages.  Technical report, Lecture 11 of Stat 928: Statistical Learning Theory, 2011.
  * [50] Mohammad Mahdi Karkhaneei and Nezam Mahdavi-Amiri.  Nonconvex weak sharp minima on Riemannian manifolds.  Journal of Optimization Theory and Applications, 183(1):85–104, 2019.
  * [51] Zuzana Kukelova, Martin Bujnak, and Tomas Pajdla.  Automatic generator of minimal problem solvers.  In European Conference on Computer Vision, pages 302–315, 2008\.
  * [52] Viktor Larsson, Kalle Astrom, and Magnus Oskarsson.  Efficient solvers for minimal problems by syzygy-based reduction.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 2383–2392, 2017.
  * [53] Viktor Larsson, Magnus Oskarsson, Kalle Astrom, Alge Wallis, Zuzana Kukelova, and Tomas Pajdla.  Beyond grobner bases: Basis selection for minimal solvers.  In IEEE Conference on Computer Vision and Pattern Recognition, 2018\.
  * [54] Gilad Lerman and Tyler Maunu.  An overview of robust subspace recovery.  Proceedings of the IEEE, 106(8):1380–1410, 2018.
  * [55] Gilad Lerman, Michael B. McCoy, Joel A. Tropp, and Teng Zhang.  Robust computation of linear models by convex relaxation.  Foundations of Computational Mathematics, 15(2):363–410, 2015.
  * [56] Chong Li, Boris S. Mordukhovich, Jinhua Wang, and Jen-Chih Yao.  Weak sharp minima on Riemannian manifolds.  SIAM Journal on Optimization, 21(4):1523–1560, 2011.
  * [57] Hongdong Li and Richard Hartley.  The 3D-3D registration problem revisited.  In IEEE International Conference on Computer Vision, pages 1–8, 2007.
  * [58] Jiayuan Li, Qingwu Hu, and Mingyao Ai.  GESAC: Robust graph enhanced sample consensus for point cloud registration.  ISPRS Journal of Photogrammetry and Remote Sensing, 167:363–374, 2020.
  * [59] Jiayuan Li, Qingwu Hu, and Mingyao Ai.  Point cloud registration based on one-point RANSAC and scale-annealing biweight estimation.  IEEE Transactions on Geoscience and Remote Sensing, pages 1–14, 2021.
  * [60] Xiao Li, Shixiang Chen, Zengde Deng, Qing Qu, Zhihui Zhu, and Anthony Man-Cho So.  Weakly convex optimization over Stiefel manifold using Riemannian subgradient-type methods.  SIAM Journal on Optimization, 31(3):1605–1634, 2021.
  * [61] Wei Lian, Lei Zhang, and Ming-Hsuan Yang.  An efficient globally optimal algorithm for asymmetric point matching.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(7):1281–1293, 2017.
  * [62] Yinlong Liu, Chen Wang, Zhijian Song, and Manning Wang.  Efficient global point cloud registration by matching rotation invariant features through translation search.  In European Conference on Computer Vision, 2018.
  * [63] Parker C. Lusk, Kaveh Fathian, and Jonathan P. How.  CLIPPER: A graph-theoretic framework for robust data association.  Technical report, arXiv:2011.10202v2 [cs.RO], 2021.
  * [64] F Landis Markley.  Attitude determination using vector observations and the singular value decomposition.  Journal of the Astronautical Sciences, 36(3):245–258, 1988.
  * [65] Haggai Maron, Nadav Dym, Itay Kezurer, Shahar Kovalsky, and Yaron Lipman.  Point registration via efficient convex relaxation.  ACM Transactions on Graphics, 35(4), 2016.
  * [66] Andreas Maurer.  A vector-contraction inequality for Rademacher complexities.  In International Conference on Algorithmic Learning Theory, pages 3–17. Springer, 2016.
  * [67] Colin McDiarmid et al.  On the method of bounded differences.  Surveys in Combinatorics, 141(1):148–188, 1989.
  * [68] Andriy Myronenko and Xubo Song.  Point set registration: Coherent point drift.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(12):2262–2275, 2010.
  * [69] Yurii Nesterov.  Lectures on Convex Optimization.  Springer, 2018.
  * [70] David Nistér.  An efficient solution to the five-point relative pose problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(6):756–770, 2004.
  * [71] Álvaro Parra, Tat-Jun Chin, Frank Neumann, Tobias Friedrich, and Maximilian Katzmann.  A practical maximum clique algorithm for matching with pairwise constraints.  Technical report, arXiv:1902.01534v2 [cs.CV], 2020.
  * [72] Álvaro Parra Bustos and Tat-Jun Chin.  Guaranteed outlier removal for point cloud registration with correspondences.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(12):2868–2882, 2018.
  * [73] Álvaro Parra Bustos, Tat-Jun Chin, Anders Eriksson, Hongdong Li, and David Suter.  Fast rotation search with stereographic projections for 3D registration.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(11):2227–2240, 2016.
  * [74] Mikael Persson and Klas Nordberg.  Lambda twist: An accurate fast robust perspective three point (p3p) solver.  In European conference on computer vision, pages 318–332, 2018\.
  * [75] Ryan A. Rossi, David F. Gleich, and Assefaw H. Gebremedhin.  Parallel maximum clique algorithms with applications to network analysis.  SIAM Journal on Scientific Computing, 37(5):C589–C616, 2015.
  * [76] Szymon Rusinkiewicz and Marc Levoy.  Efficient variants of the ICP algorithm.  In International Conference on 3D Digital Imaging and Modeling, pages 145–152. IEEE, 2001.
  * [77] Radu Bogdan Rusu, Nico Blodow, and Michael Beetz.  Fast point feature histograms (FPFH) for 3D registration.  In IEEE International Conference on Robotics and Automation, pages 3212–3217, 2009.
  * [78] Peter H Schönemann.  A generalized solution of the orthogonal Procrustes problem.  Psychometrika, 31(1):1–10, 1966.
  * [79] Jingnan Shi, Heng Yang, and Luca Carlone.  ROBIN: a graph-theoretic approach to reject outliers in robust estimation using invariants.  Technical report, arXiv:2011.03659v2 [cs.CV], 2021.
  * [80] Julian Straub, Trevor Campbell, Jonathan P How, and John W Fisher.  Efficient global point cloud alignment using bayesian nonparametric mixtures.  In IEEE Conference on Computer Vision and Pattern Recognition, pages 2941–2950, 2017.
  * [81] Lei Sun.  IRON: Invariant-based highly robust point cloud registration.  Technical report, arXiv:2103.04357v2 [cs.CV], 2021.
  * [82] Lei Sun.  RANSIC: Fast and highly robust estimation for rotation search and point cloud registration using invariant compatibility.  Technical report, arXiv:2104.09133v3 [cs.CV], 2021.
  * [83] Richard Szeliski.  Computer Vision: Algorithms and Applications.  Springer Science & Business Media, 2010.
  * [84] Chin Tat-Jun, Cai Zhipeng, and Frank Neumann.  Robust fitting in computer vision: Easy or hard?  International Journal of Computer Vision, 128(3):575–587, 2020\.
  * [85] Manolis C. Tsakiris and René Vidal.  Dual principal component pursuit.  Journal of Machine Learning Research, 19(18):1–50, 2018.
  * [86] Grace Wahba.  A least squares estimate of satellite attitude.  SIAM Review, 7(3):409, 1965.
  * [87] Heng Yang, Pasquale Antonante, Vasileios Tzoumas, and Luca Carlone.  Graduated non-convexity for robust spatial perception: From non-minimal solvers to global outlier rejection.  IEEE Robotics and Automation Letters, 5(2):1127–1134, 2020.
  * [88] Heng Yang and Luca Carlone.  A quaternion-based certifiably optimal solution to the Wahba problem with outliers.  In IEEE International Conference on Computer Vision, pages 1665–1674, 2019.
  * [89] Heng Yang and Luca Carlone.  One ring to rule them all: Certifiably robust geometric perception with outliers.  In Advances in Neural Information Processing Systems, 2020.
  * [90] Heng Yang, Ling Liang, Kim-Chuan Toh, and Luca Carlone.  STRIDE along spectrahedral vertices for solving large-scale rank-one semidefinite relaxations.  Technical report, arXiv:2105.14033 [math.OC], 2021.
  * [91] Heng Yang, Jingnan Shi, and Luca Carlone.  TEASER: Fast and certifiable point cloud registration.  IEEE Transactions on Robotics, 37(2):314–333, 2021.
  * [92] Jiaolong Yang, Hongdong Li, Dylan Campbell, and Yunde Jia.  Go-ICP: A globally optimal solution to 3D ICP point-set registration.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(11):2241–2254, 2016.
  * [93] Jiaolong Yang, Hongdong Li, and Yunde Jia.  Go-ICP: Solving 3D registration efficiently and globally optimally.  In IEEE International Conference on Computer Vision, pages 1457–1464, 2013.
  * [94] Yunzhen Yao, Liangzu Peng, and Manolis Tsakiris.  Unlabeled principal component analysis.  Advances in Neural Information Processing Systems, 2021.
  * [95] Andy Zeng, Shuran Song, Matthias Nießner, Matthew Fisher, Jianxiong Xiao, and T Funkhouser.  3DMatch: Learning the matching of local 3D geometry in range scans.  In IEEE Conference on Computer Vision and Pattern Recognition, page 4, 2017.
  * [96] Qian-Yi Zhou, Jaesik Park, and Vladlen Koltun.  Fast global registration.  In European Conference on Computer Vision, pages 766–782, 2016\.
  * [97] Zhihui Zhu, Tianyu Ding, Daniel Robinson, Manolis Tsakiris, and René Vidal.  A linearly convergent method for non-smooth non-convex optimization on the grassmannian with applications to robust subspace and dictionary learning.  In Advances in Neural Information Processing Systems, 2019.
  * [98] Zhihui Zhu, Yifan Wang, Daniel Robinson, Daniel Naiman, René Vidal, and Manolis C. Tsakiris.  Dual principal component pursuit: Improved analysis and efficient algorithms.  In Advances in Neural Information Processing Systems, 2018.

##  Appendix A Proof of Proposition 4

We consider a stronger version of Proposition 4:

######  Proposition 6.

We have
$\bm{w}^{\top}\bm{D}_{i}\bm{w}=\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}^{2}$,
where $\bm{w}$ is a quaternion representation of $\bm{R}$ of
$\eqref{eq:L1-SO3}$, and $\bm{D}_{i}\in\mathbb{R}^{4\times 4}$ is a positive
semi-definite matrix whose entries depend on $\bm{x}_{i}$, $\bm{y}_{i}$. So
Problem (8) is equivalent to

$\displaystyle\min_{\bm{w}\in\mathbb{S}^{3}}h(\bm{w}),\ \ \ \
h(\bm{w})=\sum_{i=1}^{\ell}\sqrt{\bm{w}^{\top}\bm{D}_{i}\bm{w}}.$ (15)

Moreover, $\bm{D}_{i}$ has eigenvalues $4,4,0,0$ if $\bm{x}_{i}$ and
$\bm{y}_{i}$ are normalized (that is
$\big{\|}\bm{x}_{i}\big{\|}_{2}=\big{\|}\bm{y}_{i}\big{\|}_{2}=1$).

We first recall some basics about unit quaternions, an algebraic construction
invented by Hamilton in the $1840$’s, when the notion of vector does not
exist; see the beautiful account of [1]. In our current notation, each element
$\bm{w}$ of $\mathbb{S}^{3}$ is called a unit quaternion. The most crucial
fact is that $\operatorname{SO}(3)$ is isomorphic to the $3$-sphere
$\mathbb{S}^{3}$ up to sign, that is
$\operatorname{SO}(3)\equiv\mathbb{S}^{3}/\\{\pm 1\\}$. This implies a two-to-
one correspondence between unit quaternions and $3$D rotations. Algebraically,
any $\bm{R}\in\operatorname{SO}(3)$ can be written as a $3\times 3$ matrix

$\displaystyle\footnotesize{\begin{matrix}w_{1}^{2}+w_{2}^{2}-w_{3}^{2}-w_{4}^{2}&2(w_{2}w_{3}-w_{1}w_{4})&2(w_{2}w_{4}+w_{1}w_{3})\\\
2(w_{2}w_{3}+w_{1}w_{4})&w_{1}^{2}+w_{3}^{2}-w_{2}^{2}-w_{4}^{2}&2(w_{3}w_{4}-w_{1}w_{2})\\\
2(w_{2}w_{4}-w_{1}w_{3})&2(w_{3}w_{4}+w_{1}w_{2})&w_{1}^{2}+w_{4}^{2}-w_{2}^{2}-w_{3}^{2}\\\
\end{matrix}}$

where $\bm{w}=[w_{1},w_{2},w_{3},w_{4}]^{\top}\in\mathbb{S}^{3}$. We can now
write the three entries of $\bm{R}\bm{x}_{i}$ as quadratic forms
$\bm{w}^{\top}\bm{X}_{i,1}\bm{w}$, $\bm{w}^{\top}\bm{X}_{i,2}\bm{w}$, and
$\bm{w}^{\top}\bm{X}_{i,3}\bm{w}$, respectively. Here $\bm{X}_{i,1}$,
$\bm{X}_{i,2}$, and $\bm{X}_{i,3}$ are $4\times 4$ symmetric matrices, defined
as

$\displaystyle\bm{X}_{i,1}$
$\displaystyle=\begin{bmatrix}[\bm{x}_{i}]_{1}&0&[\bm{x}_{i}]_{3}&-[\bm{x}_{i}]_{2}\\\
0&[\bm{x}_{i}]_{1}&[\bm{x}_{i}]_{2}&[\bm{x}_{i}]_{3}\\\
[\bm{x}_{i}]_{3}&[\bm{x}_{i}]_{2}&-[\bm{x}_{i}]_{1}&0\\\
-[\bm{x}_{i}]_{2}&[\bm{x}_{i}]_{3}&0&-[\bm{x}_{i}]_{1}\end{bmatrix}$ (16)
$\displaystyle\bm{X}_{i,2}$
$\displaystyle=\begin{bmatrix}[\bm{x}_{i}]_{2}&-[\bm{x}_{i}]_{3}&0&[\bm{x}_{i}]_{1}\\\
-[\bm{x}_{i}]_{3}&-[\bm{x}_{i}]_{2}&[\bm{x}_{i}]_{1}&0\\\
0&[\bm{x}_{i}]_{1}&[\bm{x}_{i}]_{2}&[\bm{x}_{i}]_{3}\\\
[\bm{x}_{i}]_{1}&0&[\bm{x}_{i}]_{3}&-[\bm{x}_{i}]_{2}\end{bmatrix}$ (17)
$\displaystyle\bm{X}_{i,3}$
$\displaystyle=\begin{bmatrix}[\bm{x}_{i}]_{3}&[\bm{x}_{i}]_{2}&-[\bm{x}_{i}]_{1}&0\\\
[\bm{x}_{i}]_{2}&-[\bm{x}_{i}]_{3}&0&[\bm{x}_{i}]_{1}\\\
-[\bm{x}_{i}]_{1}&0&-[\bm{x}_{i}]_{3}&[\bm{x}_{i}]_{2}\\\
0&[\bm{x}_{i}]_{1}&[\bm{x}_{i}]_{2}&[\bm{x}_{i}]_{3}\end{bmatrix}$ (18)

Defining
$\bm{C}_{i}:=[\bm{y}_{i}]_{1}\bm{X}_{i,1}+[\bm{y}_{i}]_{2}\bm{X}_{i,2}+[\bm{y}_{i}]_{3}\bm{X}_{i,3}$,
we get that $\bm{y}_{i}^{\top}\bm{R}\bm{x}_{i}=\bm{w}^{\top}\bm{C}_{i}\bm{w}$.
And defining

$\displaystyle\bm{D}_{i}=(\big{\|}\bm{y}_{i}\big{\|}_{2}^{2}+\big{\|}\bm{x}_{i}\big{\|}_{2}^{2})\bm{I}_{4}-2\bm{C}_{i}$
(19)

with $\bm{I}_{4}$ the $4\times 4$ identity matrix, we obtain the equality

$\displaystyle\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}^{2}$
$\displaystyle=\big{\|}\bm{y}_{i}\big{\|}_{2}^{2}+\big{\|}\bm{x}_{i}\big{\|}_{2}^{2}-2\bm{y}_{i}^{\top}\bm{R}\bm{x}_{i}$
(20) $\displaystyle=\bm{w}^{\top}\bm{D}_{i}\bm{w}.$ (21)

Since $\bm{D}_{i}$ is symmetric and $\bm{w}^{\top}\bm{D}_{i}\bm{w}\geq 0$ for
any $\bm{w}\in\mathbb{S}^{3}$, we know that $\bm{D}_{i}\in\mathbb{R}^{4\times
4}$ is positive semi-definite.

Suppose $\big{\|}\bm{y}_{i}\big{\|}_{2}=\big{\|}\bm{x}_{i}\big{\|}_{2}=1$.
Then there is at least two different $3$D rotations $\bm{R}_{1}$ and
$\bm{R}_{2}$ satisfying
$\bm{y}_{i}=\bm{R}_{1}\bm{x}_{i}=\bm{R}_{2}\bm{x}_{i}$. Thus, with the
factorization $\bm{D}_{i}=\bm{Z}_{i}\bm{Z}_{i}^{\top}$, there are at least two
quaternions $\bm{w}_{1}$ and $\bm{w}_{2}$ with $\bm{w}_{1}\neq\pm\bm{w}_{2}$
satisfying that $\bm{Z}_{i}^{\top}\bm{w}_{1}=\bm{Z}_{i}^{\top}\bm{w}_{2}=0$.
So $\operatorname{rank}(\bm{D}_{i})=\operatorname{rank}(\bm{Z}_{i})\leq 2$.
Recalling $\bm{D}_{i}=2\bm{I}_{4}-2\bm{C}_{i}$, we see that $1$ is an
eigenvalue of $\bm{C}_{i}$ that has multiplicity at least $2$. Similarly, we
can derive that
$\big{\|}\bm{y}_{i}+\bm{R}\bm{x}_{i}\big{\|}_{2}^{2}=\bm{w}^{\top}\bm{D}_{i}^{\prime}\bm{w}$
where
$\bm{D}_{i}^{\prime}=(\big{\|}\bm{y}_{i}\big{\|}_{2}^{2}+\big{\|}\bm{x}_{i}\big{\|}_{2}^{2})\bm{I}_{4}+2\bm{C}_{i}=2\bm{I}_{4}+2\bm{C}_{i}$
is positive semi-definite of rank at most $2$. That is, $-1$ is an eigenvalue
of $\bm{C}_{i}$ of multiplicity at least $2$. Concluding, $\bm{C}_{i}$ has
eigenvalues $1,1,-1,-1$ and $\bm{D}_{i}$ has eigenvalues $4,4,0,0$.

##  Appendix B Proposition 5: Proof and Interpretation

Here we provide a proof (Appendix B.1) and probabilistic interpretation
(Appendix B.2) for Proposition 5. In this section, we use the notation
$\bm{D}_{i}=\bm{Z}_{i}\bm{Z}_{i}^{\top}$ from Appendix A where we decomposed
every positive semidefinite matrix $\bm{D}_{i}$ into the product of its root
$\bm{Z}_{i}$. Since we could always normalize the point sets $\bm{y}_{i}$ and
$\bm{x}_{i}$, and then normalize $\bm{D}_{i}$, we assume without loss of
generality that $\bm{D}_{i}$ has eigenvalues $1,1,0,0$ (cf. Proposition 6). In
this situation, we can now specify that $\bm{Z}_{i}$ is a matrix of size
$4\times 2$ and it has orthonormal columns, i.e.,
$\bm{Z}_{i}^{\top}\bm{Z}_{i}=\bm{I}_{2}$. Also, we see that the objective
function (9) of interest can be rewritten as

$\displaystyle\min_{\bm{w}\in\mathbb{S}^{3}}h(\bm{w}),\ \ \ \
h(\bm{w})=\sum_{i=1}^{\ell}\big{\|}\bm{Z}_{i}^{\top}\bm{w}\big{\|}_{2}.$ (22)

Note that, if $\bm{Z}_{i}$ had a single column, then (22) is exactly the
problem of dual principal component pursuit (DPCP) [85]. On the other hand,
one could think of (22) as a group version of DPCP, as
$\big{\|}\bm{Z}_{i}^{\top}\bm{w}\big{\|}_{2}$ here promotes group sparsity. A
similar group version of DPCP was considered by [29] in the context of
homography estimation. In [29], the authors provided conditions under which
any global minimizer of (22) coincides with the ground-truth normal vector,
or, in our context, the ground-truth unit quaternion $\pm\bm{w}^{*}$. Thus,
our contribution here, if viewed from the angle of group-DPCP, is to show
that, there is actually an efficient algorithm that exactly reaches the
guaranteed ground-truth normal. We present our contribution next.

###  B.1 Proof of Proposition 5

The proof follows from Proposition 4 of [60] with some simplification for
specializing arbitrary Stiefel manifolds to $\mathbb{S}^{3}$, and with some
modification to tighten a constant factor (from $2$ to $\sqrt{2}$). We also
note that $\eta_{\min}$ and $\eta_{\max}$ are motivated from their
corresponding definitions.

Write $\bm{w}:=c_{0}\bm{w}_{0}+c^{*}\bm{w}^{*}$ with $c_{0}^{2}+(c^{*})^{2}=1$
and $\bm{w}_{0}\in\mathcal{S}^{*}$. Without loss of generality assume
$c^{*}\geq 0$. Then

$\displaystyle\operatorname{dist}(\bm{w},\pm\bm{w}^{*})$
$\displaystyle=\min\big{\\{}\sqrt{2+2c^{*}},\sqrt{2-2c^{*}}\big{\\}}$
$\displaystyle=\sqrt{2-2c^{*}}$ (23)
$\displaystyle\leq\sqrt{2-2(c^{*})^{2}}=\sqrt{2}c_{0}.$

If $i\in\mathcal{I}^{*}$ then by Proposition 4 we have

$\displaystyle\big{\|}\bm{Z}_{i}^{\top}\bm{w}^{*}\big{\|}_{2}=\sqrt{\bm{w}^{\top}\bm{D}_{i}\bm{w}}=\big{\|}\bm{y}_{i}-\bm{R}^{*}\bm{x}_{i}\big{\|}_{2}=0.$
(24)

Hence the difference $h(\bm{w})-h(\bm{w}^{*})$ is equal to

$\displaystyle
c_{0}\sum_{i\in\mathcal{I}^{*}}\big{\|}\bm{Z}_{i}^{\top}\bm{w}_{0}\big{\|}_{2}+\sum_{i\in[\ell]\backslash\mathcal{I}^{*}}\Big{(}\big{\|}\bm{Z}_{i}^{\top}\bm{w}\big{\|}_{2}-\big{\|}\bm{Z}_{i}^{\top}\bm{w}^{*}\big{\|}_{2}\Big{)}.$

By (B.1) and the definition of $\eta_{\min}$ (12), we know that

$\displaystyle
c_{0}\sum_{i\in\mathcal{I}^{*}}\big{\|}\bm{Z}_{i}^{\top}\bm{w}_{0}\big{\|}_{2}\geq\frac{k^{*}\eta_{\min}\operatorname{dist}(\bm{w},\pm\bm{w}^{*})}{\sqrt{2}}.$
(25)

By triangle inequality the second summation in the above the difference
$h(\bm{w})-h(\bm{w}^{*})$ is smaller than or equal to
$\sum_{i\in[\ell]\backslash\mathcal{I}^{*}}\big{\|}\bm{Z}_{i}^{\top}(\bm{w}-\bm{w}^{*})\big{\|}_{2}$,
but this bound satisfies

$\displaystyle\sum_{i\in[\ell]\backslash\mathcal{I}^{*}}\big{\|}\bm{Z}_{i}^{\top}(\bm{w}-\bm{w}^{*})\big{\|}_{2}\leq(\ell-k^{*})\eta_{\max}\operatorname{dist}(\bm{w},\pm\bm{w}^{*}),$

where we used
$\operatorname{dist}(\bm{w},\pm\bm{w}^{*})=\sqrt{2-2c^{*}}=\big{\|}\bm{w}-\bm{w}^{*}\big{\|}_{2}$
(B.1) and the definition of (13). We finished the proof.

###  B.2 Probabilistic Interpretation of Proposition 5

####  B.2.1 Technical Assumptions

We assume there is no noise for two reasons. First, analysis for noisy data is
more challenging and requires a full different chapter to penetrate. Second,
analysis in the noiseless case typically serves as a starting point for and
sheds enough light on analysis for noise. For example, see the trajectory of
the development from the noiseless case [85] to the noisy case [30] in the
context of DPCP.

Next, we discuss probabilistic assumptions on inliers. For an inlier index
$i\in\mathcal{I}^{*}$, each column of $\bm{Z}_{i}$ lies in the ground-truth
hyperplane $\mathcal{S}^{*}\subset\mathbb{R}^{4}$ that is perpendicular to the
ground-truth unit quaternion $\pm\bm{w}^{*}$, and the two columns of
$\bm{Z}_{i}$ span a subspace $\mathcal{S}_{i}$ of dimension $2$ that is
contained in $\mathcal{S}^{*}$. Note that any
$\bm{Z}_{i}^{\prime}\in\mathbb{R}^{4\times 2}$ whose columns are in
$\mathcal{S}_{i}\cap\mathbb{S}^{3}$ are equivalent to $\bm{Z}_{i}$ in the
sense that
$\big{\|}\bm{Z}_{i}^{\top}\bm{w}^{*}\big{\|}_{2}=\big{\|}(\bm{Z}_{i}^{\prime})^{\top}\bm{w}^{*}\big{\|}_{2}=0$.
To impose randomness assumptions on $\bm{Z}_{i}$, one could simply replace
$\bm{Z}_{i}$ by a $4\times 2$ random matrix whose columns are independently
sampled uniformly at random from the intersection
$\mathcal{S}_{i}\cap\mathbb{S}^{3}$. In fact, we need a slightly stronger
assumption:

######  Assumption 1 (randomness on inliers).

For each $i\in\mathcal{I}^{*}$, every column of $\bm{Z}_{i}$ is independently
sampled uniformly at random from the intersection
$\mathcal{S}^{*}\cap\mathbb{S}^{3}$.

This assumption destroys some good property of $\bm{Z}_{i}$: it might not be
orthonormal in general. However, it is orthonormal in expectation, i.e., it
satisfies $\mathbb{E}[\bm{Z}_{i}^{\top}\bm{Z}_{i}]=\bm{I}_{2}$. This will
suffice for our later analysis.

On the other hand, Assumption 1 simplifies matters by a lot. This can be
appreciated in comparison with a “common” approach, where one makes
assumptions on the “source data”, which are point pairs
$(\bm{y}_{i},\bm{x}_{i})$’s in our case. Let us first recall the “data flow”
from $(\bm{y}_{i},\bm{x}_{i})$ to $\bm{Z}_{i}$:

$\displaystyle(\bm{y}_{i},\bm{x}_{i})\xmapsto{\text{Proposition
\ref{prop:SO3toDPCP1}}}\bm{D}_{i}\xmapsto{\text{factorizing
$\bm{D}_{i}$}}\bm{Z}_{i}$ (26)

In view of the above flow (or graphical model), one intuitively (not very
rigorously) feels that, if $(\bm{y}_{i},\bm{x}_{i})$’s are independent, then
$\bm{Z}_{i}$’s are independent; the latter is implied by Assumption 1. On the
other hand, it seems hard to know the distribution of $\bm{Z}_{i}$’s, even if
the distribution of $(\bm{y}_{i},\bm{x}_{i})$’s is given or assumed. It is via
Assumption 1 that this challenge is circumvented and that our theorems are
developed.

Finally, we need randomness on outliers:

######  Assumption 2 (randomness on outliers).

Each column of any outlier $\bm{Z}_{j}$, where
$j\in[\ell]\backslash\mathcal{I}^{*}$, is independently sampled uniformly at
random from $\mathbb{S}^{3}$.

Since an outlier $\bm{Z}_{j}$ could be distributed arbitrarily, this
assumption is the most natural, if not the most challenging, as the outliers
try their best to mimic the distribution of inliers. Assumptions 1 and 2
(together with the noiseless assumption) are all we need for the next section.

####  B.2.2 Probabilistic Interpretation

Recall that the quantities $\eta_{\min}$, $\eta_{\max}$ of interest are equal
to

$\displaystyle\eta_{\min}$
$\displaystyle=\frac{1}{k^{*}}\min_{\bm{w}\in\mathcal{S}^{*}\cap\mathbb{S}^{3}}\sum_{i\in\mathcal{I}^{*}}\big{\|}\bm{Z}_{i}^{\top}\bm{w}\big{\|}_{2},\text{\
\ and }$ (27) $\displaystyle\eta_{\max}$
$\displaystyle=\frac{1}{\ell-k^{*}}\max_{\bm{w}\in\mathbb{S}^{3}}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\big{\|}\bm{Z}_{j}^{\top}\bm{w}\big{\|}_{2}.$
(28)

The following proposition gives probabilistic upper and lower bounds for
$\eta_{\max}$ and $\eta_{\min}$ respectively.

######  Proposition 7.

Under the assumptions of §B.2.1, we have

  1. (i)

With probability at least $1-\exp(-\zeta^{2}/2)$ it holds that

$\displaystyle\eta_{\max}\leq\frac{1}{\sqrt{2}}+\frac{(4+\zeta)}{\sqrt{\ell-k^{*}}}.$
(29)

  2. (ii)

With probability at least $1-\exp(-\zeta^{2}/2)$ it holds that

$\displaystyle\eta_{\min}\geq\sqrt{\frac{2}{3}}-\frac{(4+\zeta)}{\sqrt{k^{*}}}$
(30)

To prove Proposition 7 (cf. Appendix B.2.3), we combine the proof strategies
of [60] and [98], where both sets of the authors found inspirations from [55].
We can now see that the condition of Proposition 5,
$\alpha^{*}:=k^{*}\eta_{\min}/\sqrt{2}-(\ell-k^{*})\eta_{\max}>0$, holds with
high probability as long as

$\displaystyle\sqrt{\frac{2}{3}}k^{*}-(4+\zeta)\sqrt{k^{*}}\geq\frac{1}{\sqrt{2}}(\ell-k^{*})+(4+\zeta)\sqrt{\ell-k^{*}}.$

Ignoring lower-order terms we get the condition

$\displaystyle
k^{*}\gtrsim\frac{\sqrt{3}}{2}(\ell-k^{*})\Leftrightarrow\frac{k^{*}}{\ell}\gtrsim\frac{\sqrt{3}}{\sqrt{3}+2},$
(31)

which holds true whenever there are sufficiently many inliers. This condition
ensures the $\alpha^{*}$-sharpness, from which local linear convergence to
$\pm\bm{w}^{*}$ from a good enough initialization with proper stepsize ensues.

####  B.2.3 Details: Proof of Proposition 7

We need the following simple result, with its proof omitted.

######  Lemma 2.

If $\bm{z}=[z_{1},z_{2},z_{3},z_{4}]^{\top}$ sampled uniformly at random from
$\mathbb{S}^{3}$, we have for any $\bm{w}\in\mathbb{S}^{3}$ that

$\displaystyle\mathbb{E}\big{[}(\bm{z}^{\top}\bm{w})^{2}\big{]}=\frac{1}{4}.$
(32)

On the other hand, if $\hat{\bm{z}}$ is sampled uniformly at random from
$\mathbb{S}^{3}\cap\mathcal{S}$ where $\mathcal{S}$ is a linear subspace of
$\mathbb{R}^{4}$ of dimension $3$, then we have for every
$\hat{\bm{w}}\in\mathbb{S}^{3}\cap\mathcal{S}$ that

$\displaystyle\mathbb{E}\big{[}(\hat{\bm{z}}^{\top}\hat{\bm{w}})^{2}\big{]}=\frac{1}{3}.$
(33)

Upper Bounding $\eta_{\max}$ (i). We first prove (i) of Proposition 7.
Consider matrix $\bm{Z}\in\mathbb{R}^{4\times 2}$ whose columns are sampled
independently and uniformly at random from the $3$-sphere $\mathbb{S}^{3}$. We
will give upper bounds respectively for

$\displaystyle(\ell-k^{*})\max_{\bm{w}\in\mathbb{S}^{3}}\mathbb{E}\Big{[}\big{\|}\bm{Z}^{\top}\bm{w}\big{\|}_{2}\Big{]}\text{
and}$ (34)
$\displaystyle\max_{\bm{w}\in\mathbb{S}^{3}}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\bigg{(}\big{\|}\bm{Z}_{j}^{\top}\bm{w}\big{\|}_{2}-\mathbb{E}\Big{[}\big{\|}\bm{Z}^{\top}\bm{w}\big{\|}_{2}\Big{]}\bigg{)},$
(35)

while summing the two bounds gives an upper bound for
$(\ell-k^{*})\eta_{\max}$. For (34), Jensen’s inequality gives

$\displaystyle\max_{\bm{w}\in\mathbb{S}^{3}}\mathbb{E}\Big{[}\big{\|}\bm{Z}^{\top}\bm{w}\big{\|}_{2}\Big{]}$
$\displaystyle\leq\max_{\bm{w}\in\mathbb{S}^{3}}\sqrt{\mathbb{E}\Big{[}\big{\|}\bm{Z}^{\top}\bm{w}\big{\|}_{2}^{2}\Big{]}}$
(36)
$\displaystyle=\max_{\bm{w}\in\mathbb{S}^{3}}\sqrt{2\cdot\frac{1}{4}}=\frac{1}{\sqrt{2}}.$
(37)

To obtain (37) we used (32) and the linearity of the expectation. The second
term (35) is harder to handle, and we first consider its expectation
$\mathbb{E}[\eqref{eq:37-second}]$. We know from a standard symmetrization
argument (cf. [49], Lemma 11.4 of [17]) that, since $\bm{Z}_{j}$’s are
independent (Assumption 2), the expectation $\mathbb{E}[\eqref{eq:37-second}]$
has the following bound:

$\displaystyle\mathbb{E}[\eqref{eq:37-second}]\leq 2\
\mathbb{E}\bigg{[}\max_{\bm{w}\in\mathbb{S}^{3}}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}r_{j}\big{\|}\bm{Z}_{j}^{\top}\bm{w}\big{\|}_{2}\bigg{]},$
(38)

where $r_{j}$’s are independent Radeamacher random variables which take values
$1$, $-1$ with probabilities $1/2$ each and independent of $\bm{Z}_{j}$’s. We
also know from the vector contraction inequality (cf. Corollary 1 of [66])
that the right-hand side of (38), and thus $\mathbb{E}[\eqref{eq:37-second}]$,
is has the following bound:

$\displaystyle\mathbb{E}[\eqref{eq:37-second}]\leq$ $\displaystyle\ 2\sqrt{2}\
\mathbb{E}\bigg{[}\max_{\bm{w}\in\mathbb{S}^{3}}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\Big{(}r_{j1}\bm{Z}_{j1}^{\top}\bm{w}+r_{j2}\bm{Z}_{j2}^{\top}\bm{w}\Big{)}\bigg{]}$
$\displaystyle=$ $\displaystyle\ 2\sqrt{2}\
\mathbb{E}\bigg{[}\Big{\|}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\big{(}r_{j1}\bm{Z}_{j1}+r_{j2}\bm{Z}_{j2}\big{)}\Big{\|}_{2}\bigg{]}$
(39)

where $\bm{Z}_{j1}$’s and $\bm{Z}_{j2}$’s are the first and second columns of
$\bm{Z}_{j}$ respectively, while $r_{j1}$’s and $r_{j2}$’s are independent
Radeamacher random variables that are also independent of entries of
$\bm{Z}_{j}$’s. Applying Jensen’s inequality to (39) we get

$\displaystyle\mathbb{E}[\eqref{eq:37-second}]\leq$ $\displaystyle\ 2\sqrt{2}\
\sqrt{\mathbb{E}\bigg{[}\Big{\|}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\big{(}r_{j1}\bm{Z}_{j1}+r_{j2}\bm{Z}_{j2}\big{)}\Big{\|}_{2}^{2}\bigg{]}}$
$\displaystyle=$ $\displaystyle\ 2\sqrt{2}\
\sqrt{\mathbb{E}\bigg{[}\sum_{j\in[\ell]\backslash\mathcal{I}^{*}}\big{(}r_{j1}^{2}\bm{Z}_{j1}^{\top}\bm{Z}_{j1}+r_{j2}^{2}\bm{Z}_{j2}^{\top}\bm{Z}_{j2}\big{)}\bigg{]}}$
$\displaystyle=$ $\displaystyle\ 4\sqrt{\ell-k^{*}}$

To summarize, we have $\mathbb{E}[\eqref{eq:37-second}]\leq
4\sqrt{\ell-k^{*}}$. Treat now (35) as a function of $\bm{Z}_{j}$’s. It is
straightforward to verify that this function has bounded difference $2$ (cf.
[67]). Since $\bm{Z}_{j}$’s are independent (Assumption 2), Mcdiarmid’s Lemma
[67] or the bounded difference inequality is applicable, from which we obtain
the following probability bound:

$\displaystyle\mathbb{P}\Big{(}\eqref{eq:37-second}\geq\mathbb{E}[\eqref{eq:37-second}]+\zeta_{0}\Big{)}\leq\exp\Big{(}-\frac{\zeta_{0}^{2}}{2(\ell-k^{*})}\Big{)}.$
(40)

With $\mathbb{E}[\eqref{eq:37-second}]\leq 4\sqrt{\ell-k^{*}}$ and
$\zeta:=\zeta_{0}/\sqrt{(\ell-k^{*})}$, we get

$\displaystyle\mathbb{P}\Big{(}\eqref{eq:37-second}\leq(4+\zeta)\sqrt{\ell-k^{*}}\Big{)}\geq
1-\exp\Big{(}-\frac{\zeta^{2}}{2}\Big{)}.$ (41)

Combining this with (37) finishes proving (i).

Lower Bounding $\eta_{\min}$ (ii). Let $\bm{U}\in\mathbb{R}^{4\times 3}$ have
orthonormal columns and have $\mathcal{S}^{*}$ as its column space, then there
is a unique $\bm{v}\in\mathbb{S}^{2}$ so that $\bm{U}\bm{v}=\bm{w}$ for any
$\bm{w}\in\mathbb{S}^{3}$. Also, since for any $i\in\mathcal{I}^{*}$ every
column of $\bm{Z}_{i}$ is in $\mathcal{S}^{*}$, there is a unique
$\bm{A}_{i}\in\mathbb{R}^{3\times 2}$ with orthonormal columns satisfying
$\bm{Z}_{i}=\bm{U}\bm{A}_{i}$. Moreover, by rotation invariance we know that
each column of $\bm{A}_{i}$ is uniformly distributed on $\mathbb{S}^{2}$. As a
result, we get $\bm{Z}_{i}^{\top}\bm{w}=\bm{A}_{i}^{\top}\bm{v},\forall
i\in\mathcal{I}^{*}$, and $\eta_{\min}$ is equal to

$\displaystyle\eta_{\min}$
$\displaystyle=\frac{1}{k^{*}}\min_{\bm{v}\in\mathbb{S}^{2}}\sum_{i\in\mathcal{I}^{*}}\big{\|}\bm{A}_{i}^{\top}\bm{v}\big{\|}_{2}.$
(42)

Now, lower bounding $\eta_{\min}$ can be done in a similar way to upper
bounding $\eta_{\max}$; thus we only give a proof sketch next. Similarly to
(34) and (35), to bound $\eta_{\min}$ we will find lower bounds respectively
for the two terms

$\displaystyle
k^{*}\min_{\bm{v}\in\mathbb{S}^{2}}\mathbb{E}\Big{[}\big{\|}\bm{A}^{\top}\bm{v}\big{\|}_{2}\Big{]}\text{
and}$ (43)
$\displaystyle\min_{\bm{v}\in\mathbb{S}^{2}}\sum_{i\in\mathcal{I}^{*}}\bigg{(}\big{\|}\bm{A}_{i}^{\top}\bm{v}\big{\|}_{2}-\mathbb{E}\Big{[}\big{\|}\bm{A}^{\top}\bm{v}\big{\|}_{2}\Big{]}\bigg{)},$
(44)

where $\bm{A}$ is an i.i.d. copy of $\bm{A}_{i}$. Similarly to (37), the first
term here is bounded using (33) and Jensen’s inequality:

$\displaystyle\min_{\bm{v}\in\mathbb{S}^{2}}\mathbb{E}\Big{[}\big{\|}\bm{A}^{\top}\bm{v}\big{\|}_{2}\Big{]}\leq\sqrt{\frac{2}{3}}$
(45)

Using the symmetric argument, the vector contraction inequality, and Jensen’s
inequality, the expectation $\mathbb{E}[\eqref{eq:43-second}]$ of the second
term (44) is bounded below by $-4\sqrt{k^{*}}$. Similarly, invoking
Mcdiarmid’s Lemma gives that

$\displaystyle\
\mathbb{P}\Big{(}\eqref{eq:43-second}\leq\mathbb{E}[\eqref{eq:43-second}]-\zeta_{0}\Big{)}\leq\exp\Big{(}-\frac{\zeta_{0}^{2}}{2k^{*}}\Big{)}$
(46) $\displaystyle\Rightarrow$ $\displaystyle\
\mathbb{P}\Big{(}\eqref{eq:43-second}\geq-(4+\zeta)\sqrt{k^{*}}\Big{)}\geq
1-\exp\Big{(}-\frac{\zeta^{2}}{2}\Big{)}$ (47)

where $\zeta_{0}$ is any positive constant and we set
$\zeta:=\zeta_{0}/\sqrt{k^{*}}$. Combining (45) with the above bound finishes
the proof.

##  Appendix C Proof of Proposition 1

Since $\bm{b}^{*}$ is the rotation axis of $\bm{R}^{*}$, we have
$(\bm{R}^{*})^{\top}\bm{b}^{*}=\bm{b}^{*}$. Recall
$\bm{v}_{i}=\bm{y}_{i}-\bm{x}_{i}$ for every $i\in\mathcal{I}$. If
$i\in\mathcal{I}^{*}$ then

$\displaystyle\bm{v}_{i}^{\top}\bm{b}^{*}=(\bm{y}_{i}-\bm{x}_{i})^{\top}\bm{b}^{*}=(\bm{y}_{i}-\bm{R}^{*}\bm{x}_{i})^{\top}\bm{b}^{*},$
(48)

and further more if $(\bm{y}_{i},\bm{x}_{i})$ is an inlier pair we get that

$\displaystyle|\bm{v}_{i}^{\top}\bm{b}^{*}|=\bm{\epsilon}_{i}^{\top}\bm{b}.$
(49)

Clearly $\bm{\epsilon}_{i}^{\top}\bm{b}$ is a Gaussian random variable with
zero mean and variance $\sigma^{2}$. The rest of the proof follows from a
standard probability calculation.

##  Appendix D Interval Stabbing

Here we provide proofs for Propositions 2 and 3. Along the way we will need
multiple temporary variables to illustrate the idea; we use $a_{i,j}$’s to
denote those variables. Here, $i$ denotes the $i$-th point pair, and $j$
denotes the order in which $a_{i,j}$ appears for the first time. In §4.2.1 we
reviewed interval stabbing for closed intervals $\mathcal{J}_{i}$ of the form
$[a,b]$. One should note and verify that this can be easily extended to the
case where $\mathcal{J}_{i}$ is a finite (disjoint) union of closed intervals.

###  D.1 Proof of Proposition 2

Recall $\bm{b}=[\sin(\theta)\cos(\phi),\ \sin(\theta)\sin(\phi),\
\cos(\theta)]^{\top}$ with $\theta\in[0,\pi]$, $\phi\in[0,\pi]$. Denote by
$a_{i,1}:=[\bm{v}_{i}]_{1}\cos(\phi)+[\bm{v}_{i}]_{2}\sin(\phi)$, then
$|\bm{v}_{i}^{\top}\bm{b}|\leq c$ is equivalent to

$\displaystyle|a_{i,1}\sin(\theta)+[\bm{v}_{i}]_{3}\cos(\theta)|\leq c.$ (50)

Without loss of generality we can assume that $a_{i,1}\geq 0$. So there is a
unique $a_{i,2}\in[0,\pi]$ which satisfies

$\displaystyle\cos(a_{i,2})=\frac{[\bm{v}_{i}]_{3}}{\sqrt{[\bm{v}_{i}]_{3}^{2}+a_{i,1}^{2}}},\
\ \ \sin(a_{i,2})=\frac{a_{i,1}}{\sqrt{[\bm{v}_{i}]_{3}^{2}+a_{i,1}^{2}}}.$

Hence (50) is equivalent to

$\displaystyle|\cos(\theta-a_{i,2})|\leq c_{i},\ \ \
c_{i}:=\min\Bigg{\\{}1,\frac{c}{\sqrt{[\bm{v}_{i}]_{3}^{2}+a_{i,1}^{2}}}\Bigg{\\}}$

Since the trigonometric function $\arccos:[0,\pi]\to[-1,1]$ is decreasing and
$|\theta-a_{i,2}|\leq\pi$, the above is equivalent to

$\displaystyle a_{i,3}:=\arccos(-c_{i})\geq|\theta-
a_{i,2}|\geq\arccos(c_{i})=:a_{i,4}.$

Define $a_{i,5}=a_{i,2}-a_{i,3},\ a_{i,6}=a_{i,2}-a_{i,4},\
a_{i,7}=a_{i,2}+a_{i,4},$ and $a_{i,8}=a_{i,2}+a_{i,3}$. Then
$|\bm{v}_{i}^{\top}\bm{b}|\leq c$ is the same as

$\displaystyle\theta\in\big{(}[a_{i,5},a_{i,6}]\cup[a_{i,7},a_{i,8}]\big{)}\cap[0,\pi]$
(51)

To summarize, given $\phi\in[0,\pi]$, the $i$-th constraint of (5) requires
$\theta$ to lie in the union of some disjoint intervals defined in (51). So
maximizing (6) amounts to finding a maximal set of intervals of the form (51)
that overlap a point $\theta$, and can be solved by interval stabbing in
$O(\ell\log\ell)$ time.

###  D.2 Proof of Proposition 3

Assume that the rotation axis $\bm{b}$ of the $3$D rotation

$\displaystyle\bm{R}=\bm{b}\bm{b}^{\top}+[\bm{b}]_{\times}\sin(\omega)+(\bm{I}_{3}-\bm{b}\bm{b}^{\top})\cos(\omega)$
(52)

is given, and we now solve (7). Let
$a_{i,9}=\bm{y}_{i}^{\top}\bm{b}\bm{b}^{\top}\bm{x}_{i}$,
$a_{i,10}=\bm{y}_{i}^{\top}[\bm{b}]_{\times}\bm{x}_{i}$ and
$a_{i,11}=\bm{y}_{i}^{\top}(\bm{I}_{3}-\bm{b}\bm{b}^{\top})\bm{x}_{i}$. Then

$\displaystyle\bm{y}_{i}^{\top}\bm{R}\bm{x}_{i}=a_{i,9}+a_{i,10}\sin(\omega)+a_{i,11}\cos(\omega).$
(53)

Hence the constraint of (LABEL:eq:R-3D) can be written as

$\displaystyle\big{\|}\bm{y}_{i}\big{\|}_{2}^{2}+\big{\|}\bm{x}_{i}\big{\|}_{2}^{2}-c^{2}\leq
2\bm{y}_{i}^{\top}\bm{R}\bm{x}_{i}$ (54) $\displaystyle\Leftrightarrow$
$\displaystyle\ \ a_{i,10}\sin(\omega)+a_{i,11}\cos(\omega)\geq a_{i,12}$ (55)

where we defined
$a_{i,12}=(\big{\|}\bm{y}_{i}\big{\|}_{2}^{2}+\big{\|}\bm{x}_{i}\big{\|}_{2}^{2}-c^{2})/2-a_{i,9}$.
There is a unique angle $a_{i,13}\in[0,2\pi)$ satisfying

$\displaystyle\cos(a_{i,13})=\frac{a_{i,11}}{\sqrt{a_{i,10}^{2}+a_{i,11}^{2}}},\
\sin(a_{i,13})=\frac{a_{i,10}}{\sqrt{a_{i,10}^{2}+a_{i,11}^{2}}}$

Thus, the constraint $\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}\leq c$
of (LABEL:eq:R-3D) is the same as

$\displaystyle\cos(\omega-
a_{i,13})\geq\max\Bigg{\\{}\frac{a_{i,12}}{\sqrt{a_{i,10}^{2}+a_{i,11}^{2}}},-1\Bigg{\\}}=:a_{i,14}.$

Without loss of generality assume $a_{i,14}\leq 1$, for otherwise we could
simply ignore this constraint. Define $a_{i,15}=\arccos(a_{i,14})\in[0,\pi]$.
Since $|\omega-a_{i,13}|\in[0,2\pi]$, we consider two cases, namely $|\omega-
a_{i,13}|\leq\pi$ and $|\omega-a_{i,13}|>\pi$. In the former case, since
$\arccos$ is a decreasing function, the above constraint is equivalent to
$|\omega-a_{i,13}|\leq a_{i,15}$. In the later case the above constraint is
equivalent to

$\displaystyle\cos(2\pi-|\omega-a_{i,13}|)\geq a_{i,14}\Leftrightarrow
2\pi-|\omega-a_{i,13}|\leq a_{i,15}.$

Thus, the constraint $\big{\|}\bm{y}_{i}-\bm{R}\bm{x}_{i}\big{\|}_{2}\leq c$
of (LABEL:eq:R-3D) requires $\omega\in[0,2\pi]$ to lie in the union of the
following intervals.

$\displaystyle[a_{i,13}-a_{i,15},\ a_{i,13}+a_{i,15}]\cap[0,2\pi]$ (56)
$\displaystyle[a_{i,13}-a_{i,15}+2\pi,\ 2\pi]$ (57) $\displaystyle[0,\
a_{i,13}+a_{i,15}-2\pi].$ (58)

In the above, the invalid interval where the right endpoint is smaller than
its left endpoint, if any, should be discarded. To conclude, (7) can be solved
via interval stabbing.

##  Appendix E More Experiments

In this section we present more experiments. Besides rotation errors, we will
also use another metric for evaluation, that is success rate. Given two point
clouds as input, an algorithm succeeds if it outputs a rotation that has error
smaller than a certain threshold; by default the threshold is set to $10$
degree (as in [91]) but we will also vary it when appropriate. The success
rate is the number of success divided by the total number of experiments that
were run. This metric was referred to as recall in other related papers (cf.
[25]).

Note that, like $\operatorname{\texttt{GORE}}$ [19] and
$\operatorname{\texttt{QUASAR}}$ [88],
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ can be applied to image
stitching, because sometimes the translation is negligible and thus the scene
can be justified by a homography $\bm{H}\in\mathbb{R}^{3\times 3}$ that
involves a pure $3$D rotation $\bm{R}$, i.e., $\bm{H}=\bm{K}\bm{R}\bm{K}^{-1}$
(cf. [83]); here $\bm{K}\in\mathbb{R}^{3\times 3}$ is a matrix of intrinsic
camera parameters given by the dataset. However, we noticed that the recent
approaches MAGSAC++ [7, 8, 9] and VSAC [45] achieved surprising performance
and run in fewer than $10$ milliseconds for image stitching, hence we would
recommend them for this task.

###  E.1 Robustness on Gaussian Point Sets

In previous synthetic experiments on robust rotation search (Table 4), we
generated data by ensuring that each point pair $\bm{y}_{i}$ and $\bm{x}_{i}$
has nearly the same norm. This is for fair comparison of the methods, and it
might not be true in practice. Here we show that, without this norm
constraint, $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ can tolerate even
more outliers. In the experiment here we generated point sets
$\\{\bm{y}_{i},\bm{x}_{i}\\}_{i=1}^{\ell}$ as in Table 4 except without the
norm constraint. Then, we first perform a simple step, that removes all point
pairs $(\bm{y}_{i},\bm{x}_{i})$ which satisfy
$\big{|}\big{\|}\bm{y}_{i}\big{\|}_{2}-\big{\|}\bm{x}_{i}\big{\|}_{2}\big{|}>c$,
and then feed the remaining points to
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$. We reported the results in
Figure 3, where we observed that $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$
worked well until there are fewer than $500/10^{7}=0.005\%$ inliers.

050100150$900$$700$$500$$300$$100$# Inlier ($k^{*}$)Rotation
Error$\texttt{ARCS+}_{\texttt{OR}}$ Figure 3: Performance of
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ on Gaussian point sets with
$\ell=10^{7}$ point pairs, $\sigma=0.01$, $100$ trials.
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ works well until there are fewer
than $500/10^{7}=0.005\%$ inliers.

###  E.2 Sensitivity to The Ground-Truth Rotation

In Figure 4 we presented the sensitivity of
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ to the ground-truth rotation
$\bm{R}^{*}$. Figure 4(a) depicted that, with the ground-truth rotation angle
$\omega^{*}$ changing, the mean estimation error of
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ varied from $0.5$ to $1$, while
the standard derivation ranged from $0$ to $0.5$. One the other hand,
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$ refined the estimate from
$\operatorname{\texttt{ARCS+}}_{\texttt{O}}$, so that their combination
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ had much smaller mean error and
standard derivation, nearly imperceivable from Figure 4(a). In Figure 4(b) we
kept $\omega^{*}$ fixed and presented how the errors of
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ vary with $\theta^{*}$ and
$\phi^{*}$, the two angles for the ground-truth rotation axis $\bm{b}^{*}$; we
fixed one of them when varying the other. We observed that
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ is immune to the change of
$\phi^{*}$, as it consistently gave about $0.02$ errors and $0.01$ standard
derivation. This is expected as $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$
selects from multiple $\phi_{j}$’s a best one based on consensus maximization.
On the other hand, varying $\phi^{*}$ does make an impact on the performance
of $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$; the standard deviation
reached its peak, around $0.04$, when $\theta^{*}=\pi/4$. Theoretically
justifying the phenomenon presented here can be an interesting future work.

0.00.51.01.5$0$$90$$180$$270$$350$$\omega^{*}$$\texttt{ARCS+}_{\texttt{O}}$$\texttt{ARCS+}_{\texttt{OR}}$
(a)

0.020.040.060.08$0$$40$$90$$130$$170$$\theta^{*}$/$\phi^{*}$$\texttt{ARCS+}_{\texttt{OR}}$
($\theta^{*}$)$\texttt{ARCS+}_{\texttt{OR}}$ ($\phi^{*}$) (b)

Figure 4: Average rotation errors (in degrees) and standard deviations with
respect to the ground-truth rotation angle $\omega^{*}$ and axis
$\bm{b}^{*}=[\sin(\theta^{*})\cos(\phi^{*}),\ \sin(\theta^{*})\sin(\phi^{*}),\
\cos(\theta^{*})]^{\top}$. Experiments run with $100$ trials, $\ell=10^{5}$,
$k^{*}=1000$, $\sigma=0.01$.

(a)

(b)

(c)

$9e4$$7e4$$5e4$$3e4$$1e4$$1\%$$3\%$$5\%$$7\%$$9\%$Inlier Ratio
$k^{*}/\ell$$\ell$$0.1$$0.5$$1$ (d)

$1\%$$3\%$$5\%$$7\%$$9\%$Inlier Ratio $k^{*}/\ell$ (e)

Figure 5: Average rotation errors in degrees of different robust rotation
search approaches on medium-scale synthetic $3$D point sets of sizes varying
from $10^{4}$ to $9\times 10^{4}$ with inlier ratios ranging from $1\%$ to
$9\%$. Experiments run with $50$ trials, $\sigma=0.01$ fixed.

###  E.3 Phase Transition

In Figure 5, we showed the performances of algorithms for different inlier
ratios $k^{*}/\ell$ and different number $\ell$ of points; whiter means
smaller errors and errors larger than $1$ were truncated to $1$. The major
point we would like to clarify here is that, whether or not an algorithm can
tolerate say $99\%$ outliers might depend on the total number of points (cf.
Figures 5(a) and 5(d)), so sentences such as “our algorithm can tolerate
$99\%$ outliers” might be inaccurate, even though such description has been
widely used in recent papers. Indeed, no algorithm can tolerate $99\%$
outliers if $\ell=100$. Also, as mentioned in §4.3, one theorem of [16] has
shed light on this phenomenon. The other important observation here is that,
$\operatorname{\texttt{GNC-TLS}}$ achieved higher accuracy than
$\operatorname{\texttt{ARCS+}}_{\texttt{R}}$, although they exhibited nearly
the same breaking down points. One reason is that $\operatorname{\texttt{GNC-
TLS}}$ takes advantage of the inlier threshold $c$ as extra information. This
empirically suggests that combing $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$
and $\operatorname{\texttt{GNC-TLS}}$ might further boost the performance for
robust rotation search.

###  E.4 Robustness to Noise

Figure 6 showed that $\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ is sensitive
to noise: In particular, for $c=5.54\sigma$ fixed, the number $\ell$ of output
point pairs grows proportionally as a linear function of $\sigma$. A similar
phenomenon can be found in [91] and its follow-up works: Higher noise leads to
denser graphs, and thus to intractable maximal clique problems (recall Section
1).

$10^{5}$$2\times 10^{5}$$3\times
10^{5}$$10\%$$20\%$$30\%$$40\%$$1\%$$4\%$$7\%$$10\%$Noise Level
($\sigma$)$\ell$$\frac{\ell}{mn}$ Figure 6: Sensitivity of
$\operatorname{\texttt{ARCS+}}_{\texttt{N}}$ to noise: $\ell$ increases
linearly as $\sigma$ grows. $100$ trials, $m=1000,n=800,k^{*}=200$,
$c=5.54\sigma$.

However, $\operatorname{\texttt{ARCS+}}_{\texttt{O}}$ and
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ behave reasonably well as noise
varies. This was shown in Figure 7, where we observed that, for
$k^{*}/\ell=100/1000$, $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ are
competitive to $\operatorname{\texttt{TEASER++}}$ in terms of accuracy (Figure
7(a)) and to $\operatorname{\texttt{RANSAC}}$ in terms of speed (Figure 7(b));
$10\%$ inliers are enough for $\operatorname{\texttt{RANSAC}}$ to be fast.
Also note that the running time of $\operatorname{\texttt{TEASER++}}$
increases exponentially as noise grows, and that
$\operatorname{\texttt{GORE}}$ would achieve higher accuracy if some local
refinement methods were applied.

051015$1\%$$4\%$$7\%$$10\%$Noise Level ($\sigma$)Rotation Error (a)

$0.01$$0.1$$1$$10$$1\%$$4\%$$7\%$$10\%$Noise Level ($\sigma$)Time
(seconds)TEASER++RANSACGORE$\texttt{ARCS+}_{\texttt{OR}}$ (b)

Figure 7: Robustness of various methods to noise. $20$ trials.

###  E.5 Procrustes’s Experiments on Stanford Bunny

Here we use $\operatorname{\texttt{ARCS+}}$ for simultaneous search of
rotation & correspondences on a popular benchmark, the Stanford Bunny dataset
[27].888In view of our opening quote, Bunny here is a victim of Procrustes.
Bunny has $35947$ points with every coordinate of the points located in
$[-1,1]$ (Figure 8(a)). We randomly cut it into two parts, $\mathcal{Q}$ and
$\mathcal{P}$, of sizes $m$ and $n$ respectively and of different overlapping
ratios $k^{*}/\max\\{m,n\\}=k^{*}/m$ (Figures 8(b)-8(c) or 8(e)-8(f)). For
simplicity we set $n=\lceil 35947/2\rceil=17974$ and, $m=\lfloor
35947/2\rfloor+k^{*}$, so the exact values of $m$ and $k^{*}$ can be
calculated as per a given overlapping ratio $k^{*}/m$. We then randomly
rotated $\mathcal{P}$ and added $1\%$ random Gaussian noise to it. The goal is
to align $\mathcal{P}$ and $\mathcal{Q}$. $\operatorname{\texttt{ARCS+}}$ can
be applied directly to this task (Figures 8(d) or 8(g)). For comparison, we
gave $\operatorname{\texttt{GORE}}$ and $\operatorname{\texttt{TEASER++}}$ the
correspondences established by $\operatorname{\texttt{FPFH}}$. For all methods
we set $c=5\times 10^{-5}$. Figure 9 showed the results for different
overlapping ratios, from which we made a few observations:
$\operatorname{\texttt{ARCS+}}$ achieved higher success rates in all
experiments, while the performance of $\operatorname{\texttt{FPFH}}$, and thus
of $\operatorname{\texttt{TEASER++}}$ and $\operatorname{\texttt{GORE}}$,
improved as the overlapping ratios increased. We did not put
$\operatorname{\texttt{RANSAC}}$ into comparison here, because
$\operatorname{\texttt{FPFH}}$ often gave few to none inlier pairs for small
$k^{*}/m$ and so $\operatorname{\texttt{RANSAC}}$ used much longer time to
reach a confidence of $0.99$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 8: Bunny (8(a)) was cut through its body into two parts, $\mathcal{Q}$
(8(b)) and $\mathcal{P}$ (8(c)), with $k^{*}/m=1997/19970=10\%$ overlapping
points in blue. $\mathcal{P}$ was randomly rotated and corrupted by $1\%$
random noise. $\operatorname{\texttt{ARCS+}}$ successfully aligned
$\mathcal{Q}$ and $\mathcal{P}$ (8(d)). For a different cut through the ear of
Bunny (8(e)-8(f)), $\operatorname{\texttt{ARCS+}}$ failed (8(g)).

$0\%$$10\%$$60\%$$0$$2$$4$$6$$8$$10$Threshold (Rot. Deg.)Success
RateTEASER++GOREARCS+ (a)

$0\%$$60\%$$75\%$$0$$2$$4$$6$$8$$10$Threshold (Rot. Deg.)Success
RateTEASER++GOREARCS+ (b)

$0\%$$70\%$$80\%$$0$$2$$4$$6$$8$$10$Threshold (Rot. Deg.)Success
RateTEASER++GOREARCS+ (c)

$0\%$$70\%$$80\%$$0$$2$$4$$6$$8$$10$Threshold (Rot. Deg.)Success
RateTEASER++GOREARCS+ (d)

Figure 9: Success Rates of the methods on Stanford Bunny with different
overlapping ratios $k^{*}/m$. $1000$ trials. In each trial, Bunny was randomly
cut into two parts $\mathcal{Q}$ and $\mathcal{P}$, and $\mathcal{P}$ was then
rotated randomly and corrupted by $1\%$ random Gaussian noise.

##  Appendix F Handling The Translation Case

Utilizing ideas that have been known in prior works, it is easy to extend our
algorithms to the situation where there is an extra unknown translation. As we
did for Problems 1 and 2, we first define the two problems that we will
discuss:

######  Problem 3 (simultaneous pose and correspondences).

Let the two point sets $\mathcal{Q}$ and $\mathcal{P}$ of Problem 1 instead
satisfy

$\displaystyle\bm{q}_{i}=\bm{R}^{*}\bm{p}_{j}+\bm{t}^{*}+\bm{o}_{i,j}+\bm{\epsilon}_{i,j},$
(59)

where $\bm{t}^{*}\in\mathbb{R}^{3}$ is an extra unknown translation vector.
The task is to simultaneously estimate the rotation $\bm{R}^{*}$, translation
$\bm{t}^{*}$, and correspondences $\mathcal{C}^{*}$ from $\mathcal{Q}$ and
$\mathcal{P}$.

######  Problem 4 (robust registration).

Let the $\ell$ pairs of $3$D points
$\\{(\bm{y}_{i},\bm{x}_{i})\\}_{i=1}^{\ell}$ of Problem 2 instead satisfy

$\displaystyle\bm{y}_{i}=\bm{R}^{*}\bm{x}_{i}+\bm{t}^{*}+\bm{o}_{i}+\bm{\epsilon}_{i}.$
(60)

The task is to find $\bm{R}^{*}$, $\bm{t}^{*}$, and correspondences
$\mathcal{I}^{*}$.

We will discuss more about Problem 3 in our future work; here we focus on the
its special case, Problem 4. Specifically, we next extend our
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ algorithm to handle Problem 4
(Appendix F.1), and present its performance on the $3$DMatch dataset [95]
(Appendix F.2).

###  F.1 Extension for Robust Registration

Here we present an extension of our
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ algorithm for solving Problem 4.
In this extension, we essentially combine
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ with known techniques. Thus, the
presentation here serves more as an useful demonstration of concepts, and less
as an entirely novel insight into, or the most efficient method for, solving
Problem 4. Nevertheless, we will show in Appendix F.2 that our extension does
enjoy state-of-the-art performance on the $3$DMatch dataset [95].

We first review three crucial ingredients that are useful for solving Problem
4: translation elimination (TE), rotation elimination (RE), and outlier
removal.

Translation Elimination (TE). For each $i,j\in[\ell]$, $i>j$, define
$\mathrm{\mathbf{y}}_{ij}=\bm{y}_{i}-\bm{y}_{j}$ and
$\mathrm{\mathbf{x}}_{ij}=\bm{x}_{i}-\bm{x}_{j}$, then

$\displaystyle\mathrm{\mathbf{y}}_{ij}=\bm{R}^{*}\mathrm{\mathbf{x}}_{ij}+(\bm{o}_{i}-\bm{o}_{j})+(\bm{\epsilon}_{i}-\bm{\epsilon}_{j}).$
(61)

Here $(\mathrm{\mathbf{y}}_{ij},\mathrm{\mathbf{x}}_{ij})$ is referred to in
the literature as translation invariant measurements, as (61) no longer
involves translation. As a consequence, robust rotation search might be
performed over
$\\{(\mathrm{\mathbf{y}}_{ij},\mathrm{\mathbf{x}}_{ij})\\}_{i>j}$, yielding an
estimate of rotation and correspondences. After this, the translation can be
easily computed. A disadvantage here is that computing all
$(\mathrm{\mathbf{y}}_{ij},\mathrm{\mathbf{x}}_{ij})$’s needs $O(\ell^{2})$
time; also note though that this computation can be implemented in parallel
and thus can be efficient for medium-size datasets (_e.g_., $\ell\leq 3\times
10^{4}$).

Rotation Elimination (RE). Every inlier $(\bm{y}_{i},\bm{x}_{i})$ satisfying
(60) with $\bm{o}_{i}=0$ also necessarily satisfies

$\displaystyle\big{\|}\bm{x}_{i}+\bm{\epsilon}_{i}\big{\|}_{2}=\big{\|}\bm{y}_{i}-\bm{t}^{*}\big{\|}_{2}\Leftrightarrow\big{\|}\bm{x}_{i}\big{\|}_{2}\approx\big{\|}\bm{y}_{i}-\bm{t}^{*}\big{\|}_{2}.$
(62)

If there were no outliers, estimating $t^{*}$ from relation (62) is the
problem of source localization that appears in signal processing applications
[13]. Estimating translation from (62) in the presence of outliers is more
challenging. A possible algorithm is combining the least-squares solvers of
[13] with an iterative reweighting strategy, but this does not have global
optimality guarantee. The other approach, which we employ, is to estimate
$t^{*}$ via branch & bound, solving the following optimization problem:

$\displaystyle\max_{\mathcal{I}\subset[\ell],\bm{t}\in\mathbb{R}^{3}}\ \ \ \ \
\ \ |\mathcal{I}|$ (63) s.t. $\displaystyle\ \
\big{|}\|\bm{y}_{i}-\bm{t}\|_{2}-\|\bm{x}_{i}\|_{2}\big{|}\leq c,\ \ \forall
i\in\mathcal{I}$

If directly applying branch & bound to (63), one would branch over
$\mathbb{R}^{3}$ (cf. [62]). On the other hand, our development in §4.2
implies that branching over $\mathbb{R}^{2}$, where the first two coordinates
of $\bm{t}$ lie, suffices, as the third coordinate can be determined by
interval stabbing. In short, we solve (63) via branching over the two-
dimensional space $\mathbb{R}^{2}$ if needed. As a matter of fact, branch &
bound runs much faster even if the parameter space has smaller dimension.

######  Remark 4 (TE versus RE).

Translation elimination (TE) yields $O(\ell^{2})$ measurements, leads to the
problem of robust rotation search, and is also used in the $2$D-$3$D
perspective-three-point problem (see, _e.g_., [74]); many recent papers on
$3$D-$3$D registration used TE (see, _e.g_., [91] and its follow-up works). RE
yields $O(\ell)$ measurements, leads to a less familiar problem, and receives
fewer attention; [62] is the only paper, which we know, that uses RE (for
Problem 3).

Outlier Removal. Even though rotation or translation can be estimated
independently of each other (using TE or RE respectively), they might not be
able to handle the case of extreme outlier rates. In particular, if using TE
then the inlier ratio decreases from $k^{*}/\ell$ to
$O\big{(}(k^{*}/\ell)^{2}\big{)}$. This is why an outlier removal procedure is
needed prior to estimation. For this, create (in mind) a graph $\mathcal{G}$
with $\ell$ vertices representing the $\ell$ point pairs
$\\{(\bm{y}_{i},\bm{x}_{i})\\}_{i=1}^{\ell}$. Moreover, create an edge between
two vertices $i$ and $j$, if
$|\mathrm{\mathbf{y}}_{ij}-\mathrm{\mathbf{x}}_{ij}|\leq 2c$, where
$\mathrm{\mathbf{y}}_{ij}$ and $\mathrm{\mathbf{x}}_{ij}$ are defined in (61).
Then, find a maximum clique of $\mathcal{G}$, and remove all point pairs whose
corresponding vertices are not contained in the maximum clique. See [71, 91,
79] for more transparent discussion on this idea.

For implementation, we use the code of [71] to create $\mathcal{G}$ and
compute a maximum clique of it.

Algorithms. Having reviewed the three ingredients, we are ready to extend
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ for Problem 4. We have two
extensions, ($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$
and ($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{RE}}$,
summarized in Table 6. Both of them have the same first step, outlier removal
via finding a maximum clique from the constructed graph. Their next steps
proceed by working with point pairs that survive from outlier removal. Step 2
of ($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$ is to
eliminate the translation (TE), and step 3 is to estimate the rotation via
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ from the point pairs
$\\{(\mathrm{\mathbf{y}}_{ij},\mathrm{\mathbf{x}}_{ij})\\}_{i>j}$ (61). Step 4
of ($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$ would
estimate the translation from the remaining point pairs, with an estimated
rotation given by $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$. But we leave
step 4 unspecified, as translation estimation in this situation is
straightforward. On the other hand, step 2 of
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{RE}}$ is to
eliminate the rotation (RE), step 3 is to compute a translation $\hat{\bm{t}}$
by solving (63), and step 4 is to estimate the rotation via
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$, operating on point pairs
$(\bm{y}_{i}-\hat{\bm{t}},\bm{x}_{i})$’s. Finally, one might use an extra step
5, to refine the solution, _e.g_., by singular value decomposition.

Table 6: Two Extensions of $\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ for Problem 4. | $(\operatorname{\texttt{ARCS++}}_{\texttt{OR}})^{\texttt{TE}}$ | $(\operatorname{\texttt{ARCS++}}_{\texttt{OR}})^{\texttt{RE}}$
---|---|---
Step 1 | Outlier Removal
Step 2 | TE | RE
Step 3 |  $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ (61) | Branch & Bound (63)
Step 4 | — | $\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$
Step 5 | Local Refinement (optional)
Table 7: Success rates of methods run on the scene pairs of the 3DMatch dataset [95] for which the ground-truth transformations are provided (rotation error smaller than $10$ degree means a success [91]; see also the first paragraph of Appendix E).  Scene Type | Kitchen | Home 1 | Home 2 | Hotel 1 | Hotel 2 | Hotel 3 | Study Room | MIT Lab | Overall
---|---|---|---|---|---|---|---|---|---
# Scene Pairs | $506$ | $156$ | $208$ | $226$ | $104$ | $54$ | $292$ | $77$ | $1623$
($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{TE}}$ | $98.4\%$ | 92.9% | $89.9\%$ | $98.2\%$ | $92.3\%$ | $94.4\%$ | $\textbf{93.2}\%$ | $88.3\%$ | $94.82\%$
($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{RE}}$ | 99.0% | 92.3% | $89.4\%$ | $98.7\%$ | $91.3\%$ | $94.4\%$ | $92.5\%$ | $88.3\%$ | $94.76\%$
($\operatorname{\texttt{TEASER++}}$)∗ | 99.0% | 98.1% | $94.7\%$ | $98.7\%$ | 99.0% | $98.1\%$ | $97.0\%$ | $94.8\%$ | $97.72\%$
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$ | 98.6% | $92.3\%$ | 90.4% | 98.7% | 93.3% | $94.4\%$ | $92.5\%$ | $88.3\%$ | 94.89%
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\text{{RE}}}$ | $98.4\%$ | $91.7\%$ | 89.9% | 99.1% | 94.2% | $94.4\%$ | $92.5\%$ | $88.3\%$ | 94.82%
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)∗ | $98.4\%$ | $97.4\%$ | 95.7% | $98.7\%$ | $98.1\%$ | 100% | $97.3\%$ | 96.1% | $97.72\%$

###  F.2 Experiments on 3DMatch

Data. The 3DMatch dataset [95] contains more than $1000$ point clouds for
testing, representing $8$ different scenes (such as kitchen, hotel, etc.),
while the number of point clouds for each scene ranges from $77$ to $506$.
Each point cloud has more than $10^{5}$ points, yet in [95] there are $5000$
keypoints for each cloud. We used the pretrained
model999https://github.com/zgojcic/3DSmoothNet of the 3DSmoothNet [35] to
extract descriptors from these key points, and matched them using the Matlab
function pcmatchfeatures, with its parameter MatchThreshold set to the maximum
$1$. It remains to solve Problem 4 using these hypothetical correspondences.

Metrics. We report success rates of the methods. Success rates were defined in
the beginning of Appendix E. The default threshold $10$ on rotation degrees is
the one that was used in $\operatorname{\texttt{TEASER++}}$ [91]. We do not
report errors in terms of translation for two reasons: i) rotation search is
the main theme of the paper, ii) if the rotation is estimated accurately, then
so will be the translation (see, _e.g_., algorithms of [91]).

Methods. We apply
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$ and
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{RE}}$ to restore
the rotation and translation from these correspondences. We use singular value
decomposition as an extra step 5 for
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{RE}}$ to refine
the solution and account for inaccuracy of translation estimation via branch &
bound (63). For reference, we also apply
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)∗, which uses the ground-truth
translation $\bm{t}^{*}$ and point pairs
$(\bm{y}_{i}-\bm{t}^{*},\bm{x}_{i})$’s to estimate a rotation via
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$.

$0\%$$25\%$$50\%$$75\%$$100\%$$0$$2$$4$$6$$8$$10$Rotation Err.
(Deg.)$\texttt{ARCS++}_{\texttt{OR}}$TEASER++ (a)

$0\%$$25\%$$50\%$$75\%$$100\%$$0$$2$$4$$6$$8$$10$Rotation Err.
(Deg.)$\texttt{ARCS++}_{\texttt{OR}}$TEASER++ (b)

$0\%$$25\%$$50\%$$75\%$$100\%$$0$$2$$4$$6$$8$$10$Rotation Err. (Deg.)Success
Rate$\texttt{ARCS++}_{\texttt{OR}}$TEASER++ (c)

Figure 10: Success rates of $\operatorname{\texttt{TEASER++}}$ and
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ on the 3D Match dataset, using
either estimated translation (Fig. 10(a)) or $\operatorname{\texttt{TIM}}$s
(Fig. 10(b)) or ground-truth translation (Fig. 10(c)).

We compare our algorithms with $\operatorname{\texttt{TEASER++}}$ [91].
Similarly, we use three versions of $\operatorname{\texttt{TEASER++}}$. The
first version is ($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{TE}}$. This
is the standard $\operatorname{\texttt{TEASER++}}$, and the difference between
($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{TE}}$ and
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{TE}}$ is that,
($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{TE}}$ estimates the rotation
by $\operatorname{\texttt{GNC-TLS}}$, not
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$. The second version is
($\operatorname{\texttt{TEASER++}}$)${}^{\texttt{TE}}$, where we treat
$\operatorname{\texttt{TEASER++}}$ as a robust rotation search method and let
it play the role of $\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ in
($\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$)${}^{\texttt{RE}}$. The third
version is ($\operatorname{\texttt{TEASER++}}$)∗, where we assume the ground-
truth translation $\bm{t}^{*}$ is given and run
$\operatorname{\texttt{TEASER++}}$ on $(\bm{y}_{i}-\bm{t}^{*},\bm{x}_{i})$’s.
Finally, we did not compare other methods here, as
$\operatorname{\texttt{TEASER++}}$ currently has the best performance (to the
best of our knowledge) on the 3DMatch dataset, see [91] for comparison with
optimization-based methods, and also read from [25] the success rates (recall)
of other deep learning methods.

Results. Following [91], we set $c=0.05$. We presented results in Table 7 and
Figure 10. In Table 7 we observed that
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ and
$\operatorname{\texttt{TEASER++}}$ have very close performance, although
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ has slight advantage (_e.g_.,
in $12$ cases in bold $\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ has
higher success rates). In terms of running times,
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ is slower than
$\operatorname{\texttt{TEASER++}}$. One reason is that we used an industrial-
strength implementation101010https://github.com/MIT-SPARK/TEASER-plusplus of
$\operatorname{\texttt{TEASER++}}$, while
$\operatorname{\texttt{ARCS++}}_{\texttt{OR}}$ was implemented in plain
Matlab. This suggests our current idea of extending
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ into the translation case might
be sub-optimal, and will motivate us to design even faster algorithms for that
purpose, which though will require serious innovations. Finally, in Figure 10,
we reported the success rates averaged over all testing scenes of $3$DMatch
and with the threshold (rotation degree) varying from $0$ to $100$. This
delivers the same message that our direct extension of
$\operatorname{\texttt{ARCS+}}_{\texttt{OR}}$ maintains a state-of-the-art
performance for solving Problem 4.