Outlier-Robust Spatial Perception: Hardness, General-Purpose Algorithms, and Guarantees

TL;DR

This paper investigates the computational hardness of outlier rejection in spatial perception, introduces a novel adaptive trimming algorithm with theoretical guarantees, and demonstrates its effectiveness across various robotics applications.

Contribution

It proves the inapproximability of simple outlier rejection problems, provides per-instance bounds for solution quality, and proposes a versatile algorithm that outperforms existing methods.

Findings

Adaptive trimming outperforms state-of-the-art methods in experiments.

The problem of outlier rejection is computationally hard to approximate.

The proposed method is applicable to multiple spatial perception tasks.

Abstract

Spatial perception is the backbone of many robotics applications, and spans a broad range of research problems, including localization and mapping, point cloud alignment, and relative pose estimation from camera images. Robust spatial perception is jeopardized by the presence of incorrect data association, and in general, outliers. Although techniques to handle outliers do exist, they can fail in unpredictable manners (e.g., RANSAC, robust estimators), or can have exponential runtime (e.g., branch-and-bound). In this paper, we advance the state of the art in outlier rejection by making three contributions. First, we show that even a simple linear instance of outlier rejection is inapproximable: in the worst-case one cannot design a quasi-polynomial time algorithm that computes an approximate solution efficiently. Our second contribution is to provide the first per-instance sub-optimality bounds to assess the approximation quality of a given outlier rejection outcome. Our third contribution is to propose a simple general-purpose algorithm, named adaptive trimming, to remove outliers. Our algorithm leverages recently-proposed global solvers that are able to solve outlier-free problems, and iteratively removes measurements with large errors. We demonstrate the proposed algorithm on three spatial perception problems: 3D registration, two-view geometry, and SLAM. The results show that our algorithm outperforms several state-of-the-art methods across applications while being a general-purpose method.

Tables

Table 1. TABLE I: Parameters for ADAPT .

Experiment	$T$	$\delta$	$g$	$\gamma$
Bunny	2	$1\text{×}{10}^{-4}$	10	0.99
ETH	2	$1\text{×}{10}^{-2}$	10	0.99
MIT	2	10	20	0.99
Intel	2	60	20	0.5
CSAIL	2	60	20	0.5
Grid	2	20	20	0.99
Sphere 2500	2	10	200	0.9
Two-view (synthetic)	2	$1.5\text{×}{10}^{-4}$	20	0.99
EuRoC	2	$1.5\text{×}{10}^{-4}$	20	0.99

Equations

$$\underset{x\in \mathbb{X}}{\min }\sum _{i\in \mathcal{M}}\Vert h_i(y_i,x){\Vert }^2,$$

$$\underset{\mathcal{O}\subseteq \mathcal{M}}{\min }\underset{x\in \mathbb{X}}{\min }|\mathcal{O}|,\text{s.t.}\sum _{i\in \mathcal{M}\setminus \mathcal{O}}\Vert h_i(y_i,x){\Vert }^2\le {ϵ}_{\mathcal{M}\setminus \mathcal{O}},$$

$$\underset{x\in {\mathbb{R}}^n}{\min }\underset{\mathcal{O}\subseteq \mathcal{M}}{\min }|\mathcal{O}|,\text{s.t.}\sum _{i\in \mathcal{M}\setminus \mathcal{O}}\Vert y_i-a_i^{\top }x{\Vert }^2\le {ϵ}_{\mathcal{M}\setminus \mathcal{O}}.$$

$$\underset{\begin{array}{c}R\in SO(3)\\ t\in {\mathbb{R}}^3\end{array}}{\min }\sum _{(i,j)\in \mathcal{M}}\Vert R{p}_i+t-p_j^{\prime }{\Vert }^2.$$

$$\underset{\begin{array}{c}R\in SO(3)\\ t\in {\mathbb{R}}^3\end{array}}{\min }\underset{\mathcal{O}\subseteq \mathcal{M}}{\min }|\mathcal{O}|,\text{s.t.}\sum _{(i,j)\in \mathcal{M}\setminus \mathcal{O}}\Vert R{p}_i+t-p_j^{\prime }{\Vert }^2\le {ϵ}_{\mathcal{M}\setminus \mathcal{O}}.$$

$$\underset{\begin{array}{c}R\in SO(3)\\ t\in {\mathbb{S}}^2\end{array}}{\min }\sum _{(i,j)\in \mathcal{M}}{\left[f_i^{\top }\left(t\times (R{f}_j^{\prime })\right)\right]}^2,$$

$$\underset{\begin{array}{c}R\in SO(3)\\ t\in {\mathbb{S}}^2\end{array}}{\min }\underset{\mathcal{O}\subseteq \mathcal{M}}{\min }|\mathcal{O}|,\text{s.t.}\sum _{(i,j)\in \mathcal{M}}{\left[t^{\top }\left(f_i\times (R{f}_j^{\prime })\right)\right]}^2\le {ϵ}_{\mathcal{M}\setminus \mathcal{O}}.$$

$$\underset{\begin{array}{c}{T}_i\in \mathrm{SE}(3)\\ i=1,\dots ,n\end{array}}{\min }\sum _{(i,j)\in \mathcal{M}}\Vert T_j-T_i{\bar{T}}_{ij}{\Vert }_F^2,$$

$$\underset{\begin{array}{c}{T}_i\in \mathrm{SE}(3)\\ i=1,\dots ,n\end{array}}{\min }\underset{\mathcal{O}\subseteq {\mathcal{E}}_{lc}}{\min }|\mathcal{O}|,\text{s.t.}$$

$$\sum _{(i,j)\in {\mathcal{E}}_o}\Vert T_j-T_i{\bar{T}}_{ij}{\Vert }_F^2+\sum _{(i,j)\in {\mathcal{E}}_{lc}\setminus \mathcal{O}}\Vert T_j-T_i{\bar{T}}_{ij}{\Vert }_F^2\le {ϵ}_{\mathcal{M}\setminus \mathcal{O}},$$

$$\frac{r(\mathcal{O})-r_{|\mathcal{O}|}^{\star }}{r(\varnothing )-r_{|\mathcal{O}|}^{\star }}\le {\chi }_{\mathcal{O}},$$

$${\chi }_{\mathcal{O}}\doteq \frac{r(\mathcal{O})}{r(\varnothing )-r(\mathcal{O})}.$$

$$\frac{r(\mathcal{O})-r^{\star }}{r(\varnothing )-r^{\star }}\le {\chi }_{\mathcal{O}}.$$

$$\underset{d\in {\mathbb{R}}^l}{\min }\Vert d{\Vert }_0,\text{s.t.}\Vert Ud-z{\Vert }_2\le \Delta .$$

$$\underset{d\in {\mathbb{R}}^l}{\min }\Vert d{\Vert }_0,\text{s.t.}Ud=1_{m(l)}.$$

$$\underset{d\in {\mathbb{R}}^l,x\in {\mathbb{R}}^n}{\min }\Vert d{\Vert }_0,\text{s.t.}y=Ax+d.$$

$$\min_{{\cal O}\in{\cal M}\!,~{}x\in\mathbb{R}^{n}}\;\;|{\cal O}|,\;\text{s.t.}\;\;\|y_{{{\cal M}}\setminus{{\cal O}}}-A_{{{\cal M}}\setminus{{\cal O}}}x\|^{2}\leq\epsilon.\$$

$$f^{\star }$$

$$\le r(\varnothing )$$

$$=r(\varnothing )+r(\mathcal{O})-r(\mathcal{O})$$

$$=f(\mathcal{O})+r(\mathcal{O}).$$

$$\frac{f(\mathcal{O})}{f^{\star }}$$

$$\ge 1-\frac{r(\mathcal{O})}{f(\mathcal{O})},$$

$$\frac{f(\mathcal{O})}{f^{\star }}$$

$$r(\varnothing )-r(\mathcal{O})$$

$${\chi }_{\mathcal{O}}(r(\varnothing )-r_{|\mathcal{O}|}^{\star })$$

$$\frac{r(\mathcal{O})-r^{\star }}{r(\varnothing )-r^{\star }}\le \frac{r(\mathcal{O})-r_{|\mathcal{O}|}^{\star }}{r(\varnothing )-r_{|\mathcal{O}|}^{\star }},$$

Full text

I Introduction

Spatial perception is concerned with the estimation of a geometric model that describes the state of the robot, and/or the environment the robot is deployed in. As such, spatial perception includes a broad set of robotics problems, including motion estimation [1], object detection, localization and tracking [2], multi-robot localization [3], dense reconstruction [4], and Simultaneous Localization and Mapping (SLAM) [5]. Spatial perception algorithms find applications beyond robotics, including virtual and augmented reality, and medical imaging [2], to mention a few.

Safety-critical applications, including self-driving cars, demand robust spatial perception algorithms that can estimate correct models (and assess their performance) in the presence of measurement noise and outliers. While we currently have several approaches that can tolerate large measurement noise (e.g., [6, 7, 8]), these algorithm tend to catastrophically fail in the presence of outliers resulting from incorrect data association, sensor malfunction, or even adversarial attacks.

In this paper, we focus on the analysis and design of outlier-robust general-purpose algorithms for robust estimation applied to spatial perception. Our proposal is motivated by three observations. First, recent years have seen a convergence of the robotics community towards optimization-based approaches for spatial perception. Therefore, despite the apparent heterogeneity of the perception landscape, it is possible to develop general-purpose methods to reject outliers (e.g., M-estimators [9] and consensus maximization [10] can be thought as general estimation tools). Second, the research community has developed global solutions to many perception problems without outliers, from well-established techniques for point cloud registration [8], to very recent solvers for SLAM [6] and two-view geometry [7]. These global solvers offer unprecedented opportunities to tackle robust estimation with outliers. Third, the literature still lacks a satisfactory answer to provably-robust spatial perception.

The literature on outlier-robust spatial perception is currently divided between fast approaches (that mainly work in the low-outlier regime, without performance guarantees) and provably-robust approaches (that can tolerate many outliers, but have exponential run-time). While we postpone a comprehensive literature review to Section VII, it is instructive to briefly review this dichotomy. Fast approaches include RANSAC [11], M-estimators [9], and measurement-consistency checking [12, 13]. These methods fall short of providing performance guarantees. In particular, RANSAC is known to become slow and brittle with high outlier rates ($>50\%$) [10], and does not scale to high-dimensional problems, while M-estimators have a breakdown point of zero, meaning that a single “bad” outlier can compromise the results. On the other hand, provably-robust methods, typically based on branch-and-bound [14, 15, 16, 17, 10], can tolerate more than $50\%$ of outliers [18], but do not scale to large problems and are relatively slow for robotics applications. Overall, the first goal of this paper is to understand whether we can resolve this divide, and design algorithms that are both efficient and provably robust.

Contributions. We propose a Minimal Trimmed Squares (MTS) formulation for outlier-robust estimation. MTS encapsulate a wide spectrum of commonly-used outlier-robust formulations in the literature, such as the popular maximum consensus [19], Linear Trimmed Squares [20], and truncated least-squares [21]. In particular, MTS aims to compute a “good” estimate by rejecting a minimal set of measurements.

Our first contribution (Section III) is a negative result: we show that outlier rejection is inapproximable. In the worst-case, there exist no quasi-polynomial algorithm that can compute (even an approximate) solution to the outlier rejection problem. We prove that this remains true, surprisingly, even if the algorithm knows the true number of outliers and even if we allow the algorithm to reject more measurements than necessary. Our conclusions largely extend previously-known negative results [19], which already ruled-out the possibility of designing polynomial-time approximation methods.

Our second contribution (Section IV) is to derive the first per-instance sub-optimality bounds to assess the quality of a given outlier rejection solution. While in the worst case we expect efficient algorithms to perform poorly, we can still hope that in typical problem instances a polynomial-time algorithm can compute good solutions, and we can use the proposed sub-optimality bounds to assess the performance of such an algorithm. Our bounds are algorithm-agnostic (e.g., they also apply to RANSAC) and can be computed efficiently.

Our third contribution (Section V) is a general-purpose algorithm for outlier rejection, named Adaptive Trimming (ADAPT). ADAPT leverages recently-proposed global solvers that solve outlier-free problems and adaptively removes measurements with large residual errors. Despite its simplicity, our experiments show that it outperforms RANSAC and even specialized state-of-the-art methods for robust estimation.

We conclude the paper by providing an experimental evaluation across multiple spatial perception problems (Section VI). The experiments show that ADAPT can tolerate up to 90% outliers in 3D registration (with a run-time similar to existing methods), and up to 50% outliers in two-view geometry and most SLAM datasets. The experiments also show that the proposed sub-optimality bounds are effective in assessing the outlier rejection outcomes. We report extra results and proofs in the Appendix.

II Outlier Rejection: a Minimally Trimmed Squares Formulation

Many estimation problems in robotics and computer vision can be formulated as non-linear least squares problems:

[TABLE]

where we are given measurements $y_{i}$ of an unknown variable $x$, with $i\in{\cal M}$ (${\cal M}$ is the measurement set), and we want to estimate $x$, potentially restricted to a given domain $\mathbb{X}$ (e.g., $x$ is a pose, and $\mathbb{X}$ is the set of 3D poses). The least squares problem in eq. (1) looks for the $x$ that minimizes the (squares of) the residual errors $h_{i}(y_{i},x)$, where the $i$-th residual error captures how well $x$ explains the measurement $y_{i}$. The problem in eq. (1) typically results from maximum likelihood and maximum a posteriori estimation [5, 22], under the assumption that the measurement noise is Gaussian.

Both researchers and practitioners are well-aware that least squares formulations are sensitive to outliers, and that the estimator in eq. (1) fails to produce a meaningful estimate of $x$ in the presence of gross outliers $y_{i}$. Therefore, in this paper we address the following question:

Can we compute an accurate estimate of $x$ that is insensitive to the presence of outlying measurements?

We formulate the resulting robust estimation problem as the problem of selecting a small number of outliers, such that the remaining measurements (the inliers) can be explained with small error. In other words, a good estimate (in the presence of outliers) is one that explains as many measurements as possible while disregarding outliers. This intuition leads to the following formulation.

Problem 1** (Minimally Trimmed Squares (MTS))**

Let ${\cal M}$ denote a set of measurements of an unknown variable $x$, and let $y_{i}$ denote the $i$-th measurement. Also denote with $h_{i}(y_{i},x)$ the residual error that quantifies how well $x$ fits the measurement $y_{i}$. Then, the minimally trimmed squares problem consists in estimating the unknown variable $x$ by solving the following optimization problem:

[TABLE]

where one searches for the smallest set of outliers ${\cal O}$ ($|\hskip 2.27621pt\cdot\hskip 2.27621pt|$ is the cardinality of a set) among the given measurements ${\cal M}$, such that the remaining measurements ${\cal M}\setminus{\cal O}$ (i.e., the inliers) can be explained with small error, i.e., $\sum_{i\in{\cal M}\setminus{\cal O}}\|h_{i}(y_{i},x)\|^{2}\leq\epsilon_{{\cal M}\setminus{\cal O}}$ for some $x\in\mathbb{X}$, and where $\epsilon_{{\cal M}\setminus{\cal O}}$ is a given outlier-free bound. $\lrcorner$

Example 1** (Robust linear estimation and bound $\epsilon_{{\cal M}\setminus{\cal O}}$)**

In linear estimation one wishes to recover a parameter $x\in\mathbb{R}^{n}$ from a set of noisy measurements $y_{i}=a_{i}^{\top}x+d_{i}$, $i\in{\cal M}$, where $a_{i}$ is a known vector, and $d_{i}\in{{\mathbb{R}}}$ models the unknown measurement noise. Some of the measurements (the inliers) are such that the corresponding noise $d_{i}$ can be assumed to follow a Gaussian distribution, while others (the outliers) may be affected by large noise. Therefore, our MTS estimator can be written as:

[TABLE]

Evidently, $\epsilon_{{\cal M}\setminus{\cal O}}$ must increase with the number of inliers, since each inlier adds a positive summand $\|y_{i}-a_{i}^{\top}x\|^{2}$ due to the presence of noise. Moreover, since the sum is restricted to the inliers, for which the noise is assumed to be Gaussian, we can compute the desired outlier-free bound explicitly: if $d_{i}$ follows a Gaussian distribution, then each $\|y_{i}-a_{i}^{\top}x\|^{2}_{2}$ follows a $\chi^{2}$ distribution with $1$ degree of freedom. Thus, with desired probability $p_{\epsilon}$ (e.g., $0.99$), $\|y_{i}-a_{i}^{\top}x\|^{2}_{2}\leq\epsilon$ where $\epsilon$ is the $p_{\epsilon}$-quantile of the $\chi^{2}$ distribution, and the outlier-free bound is $\epsilon_{{\cal M}\setminus{\cal O}}=|{\cal M}\setminus{\cal O}|\epsilon$. $\lrcorner$

Remark 2** (Generality ​and ​applicability)**

In this paper we address robustness in non-linear and non-convex estimation problems as the ones arising in robotics and computer vision. Therefore, while the linear estimation Example 1 is instructive (and indeed we will prove in Section III that even in such a simple case, it is not possible to even approximate the MTS estimator in polynomial time), the algorithms and bounds presented in this paper hold for any function $h_{i}(y_{i},x)$ and any domain $\mathbb{X}$. In contrast with related work [23, 24], we do not assume the number of outliers to be known in advance (an unrealistic assumption in perception problems). Indeed, our MTS formulation looks for the smallest set of outliers. Finally, while the formulation (2) requires to set an outlier-free threshold, we will propose an algorithm (Section V) that will automatically compute a suitable threshold without any prior knowledge about the measurement noise. $\lrcorner$

In summary, MTS is a general non-linear and non-convex outlier rejection framework. We exemplify its generality by discussing its application to three core perception problems: 3D registration, two-view geometry, and SLAM.

II-A *Outlier rejection for robust spatial perception:

3D registration, two-view geometry, and SLAM*

Here we review three core problems in spatial perception, and show how to tailor the framework of Section II to these examples. The expert reader can safely skip this section.

Outlier rejection for 3D registration. Point cloud registration consists in finding the rigid transformation that aligns two point clouds. Formally, we are given two sets of points ${\cal P}\doteq\{p_{1},\ldots,p_{n}\}$ and ${\cal P}^{\prime}\doteq\{p^{\prime}_{1},\ldots,p^{\prime}_{n}\}$ (with $p_{i},p^{\prime}_{i}\in{{\mathbb{R}}^{3}}$, for $i=1,\ldots,n$), as well as a set ${\cal M}$ of putative correspondences $(i,j)$, such that the point $p_{i}\in{\cal P}$ and the point $p^{\prime}_{j}\in{\cal P}^{\prime}$ are (putatively) related by a rigid transformation, for all $(i,j)\in{\cal M}$. Point correspondences are typically obtained by descriptor matching [18].

Given the points and the putative correspondences, 3D registration looks for a rotation $R\in SO(3)$, and a translation $t\in\mathbb{R}^{3}$, that align (i.e., minimize the sum of the squared distances between) corresponding points:

[TABLE]

The problem in eq. (4) can be solved in closed form [8]. However, eq. (4) fails to produce a reasonable pose (rotation and translation) estimate when some of the correspondences are outliers [18, 25], and related work resorts to robust estimators (reviewed in Section VII). Here we rephrase robust registration as an MTS problem:

[TABLE]

Outlier rejection for two-view geometry. Two-view geometry estimation consists in finding the relative pose (up to scale) between two camera images picturing a static scene, and it is crucial for motion estimation [26], object localization [26], and reconstruction [27, Chapter 1]. We consider a feature-based calibrated setup where the camera calibration is known and one extracts features (keypoints) ${\cal F}=\{f_{1},\ldots,f_{n}\}$ and ${\cal F}^{\prime}=\{f^{\prime}_{1},\ldots,f^{\prime}_{n}\}$ from the first and second image, respectively. We are also given a set of putative correspondences ${\cal M}$ between pairs of features $(i,j)$, such that features $f_{i}$ and $f^{\prime}_{j}$ (putatively) picture the same 3D point observed in both images.

Given the features and the putative correspondences, two-view geometry looks for the rotation $R\in SO(3)$ and the translation $t\in{{\mathbb{R}}^{3}}$ (up to scale) that minimizes the violation of the epipolar constraint:

[TABLE]

where $t$ is restricted to the unit sphere $\mathbb{S}^{2}$ to remove the scale ambiguity. In the absence of outliers, problem (6) can be solved globally using convex relaxations [7].

In the presence of outliers, the non-robust formulation in eq. (6) fails to compute accurate pose estimates, hence we rephrase two-view estimation as an MTS problem:

[TABLE]

Outlier rejection for SLAM. Here we consider one of the most popular SLAM formulations: Pose graph optimization (PGO). PGO estimates a set of robot poses $T_{i}\in\mathrm{SE}(3)$ ($i=1,\ldots,n$) from pairwise relative pose measurements $\bar{T}_{ij}\in\mathrm{SE}(3)$ between pairs of poses $(i,j)\in{\cal M}$. The measurement set ${\cal M}$ includes odometry (ego-motion) measurements as well as loop closures. In the absence of outliers, one can compute the pose estimates as:

[TABLE]

where $\|\hskip 2.27621pt\cdot\hskip 2.27621pt\|_{F}^{2}$ denotes the Frobenius norm.

In practice, many loop closure measurements are outliers (e.g., due to failures in place recognition). Therefore, we rephrase PGO as an MTS problem over the loop closures:

[TABLE]

where we split the measurement set ${\cal M}$ into odometric edges ${\cal E}_{o}$ (these can be typically trusted), and loop closures ${\cal E}_{lc}$ (typically containing outliers).

III Outlier Rejection is Inapproximable

We show that MTS is inapproximable even by quasi-polynomial-time algorithms. To this end, we find worst-case instances for which there is no algorithm that can reject a few measurements to achieve a prescribed residual error $\epsilon$ (subject to a widely believed conjecture in complexity theory, similar to $NP\neq P$). We start with some definitions and present our key result in Theorem 5.

Definition 3** (Approximability)**

Consider the MTS Problem 1. Let ${\cal O}^{\star}$ be an optimal solution, let $k^{\star}\doteq|{\cal O}^{\star}|$, and $\epsilon\doteq\epsilon_{{\cal M}\setminus{\cal O}^{\star}}$, that is, $\epsilon$ is the outlier-free bound when the measurements ${\cal O}^{\star}$ are the rejected outliers. Also, consider a number $\lambda>1$. We say that an algorithm makes MTS $(\lambda,\epsilon)$-approximable if it returns a set ${\cal O}$, and a parameter $x$, such that $|{\cal O}|\leq\lambda k^{\star}$ and $\sum_{i\in{\cal M}\setminus{\cal O}}\|h_{i}(y_{i},x)\|^{2}\leq\epsilon$. $\lrcorner$

The definition of $(\lambda,\epsilon)$-approximability allows some slack in the quality of the MTS’s solution: rather than solving Problem (2) exactly ($\lambda=1$), Definition 3 only requires, for MTS to be approximable, to find an algorithm that computes an estimate close to the optimal solution. Indeed, Definition 3 includes algorithms that can reject more outliers than necessary (since $\lambda k^{\star}>k^{\star}$).

Definition 4** (Quasi-polynomial algorithm)**

An algorithm is said to be quasi-polynomial if it runs in $2^{O\left[(\log m)^{c}\right]}$ time, where $m$ is the size of the input and $c$ is constant. $\lrcorner$

Any polynomial algorithm is also quasi-polynomial, since $m^{k}=2^{k\log m}$. Yet, a quasi-polynomial algorithm is asymptotically faster than an exponential-time algorithm, since exponential algorithms run in $O(2^{m^{c}})$ time, for some $c>0$.

Theorem 5** (Inapproximability)**

Consider the linear MTS problem (3). Let $x^{\star}$ be the optimal value of the variable to be estimated, $m$ be the number of measurements ($m\doteq|{\cal M}|$), ${\cal O}^{\star}$ be the optimal solution, and set $k^{\star}\!\doteq\!|{\cal O}^{\star}|$. Then, for any $\delta\!\in\!(0,1)$, there exist a polynomial $\lambda_{1}(m)$ and a $\lambda_{2}(m)=2^{\Omega(\log^{1-\delta}m)}$ and instances of MTS (i.e., measurements $y_{i}$, vectors $a_{i}$, and outlier-free bound $\epsilon$) where $\epsilon=\lambda_{2}(m)$, such that unless ${\rm NP}{\in}{\rm BPTIME(\textit{m}^{\rm poly\log\textit{m}})}$,111The complexity hypothesis ${\rm NP}{\notin}{\rm BPTIME(\textit{m}^{\rm poly\log\textit{m}})}$ means there is no randomized algorithm which outputs solutions to problems in ${\rm NP}$ with probability $2/3$, after running for $O(m^{(\log m)^{c}})$ time, for a constant $c$ [28]. there is no quasi-polynomial algorithm making MTS $(\lambda_{1}(m),\lambda_{2}(m))$-approximable. This holds true even if the algorithm knows $k^{\star}$​, and that $x^{\star}$ ​exist.  $\lrcorner$

Theorem 5 stresses the extreme hardness of MTS. Even if we inform the algorithms with the true number of outliers, it is impossible in the worst-case for even quasi-polynomial algorithms to find a good set of inliers. Surprisingly, this remains true even if we allow the algorithms to cheat by rejecting more measurements than $k^{\star}$ (i.e., $\lambda_{1}k^{\star}$).

Thinking beyond the worst-case, our inapproximability result suggests that to obtain a good solution efficiently, our only hope is that nature (which picks the outliers) is not adversarial, thus fast algorithms can compute good solutions in practice. Hence, it becomes important to derive per-instance bounds that, for a given MTS problem (i.e., given $y_{i},h_{i}$, and $\epsilon_{{\cal M}\setminus{\cal O}}$ in (1)), can evaluate how far an algorithm is from the optimal MTS solution. In order words, since we cannot guarantee than any efficient algorithm will do well in the worst-case, we are happy with evaluating (a posteriori) if an algorithm computed a good solution for a given problem instance. For this reason, in the next section we develop the first per-instance sub-optimality bound for Problem 1.

IV Performance Guarantees

We present the first per-instance (i.e., a posteriori) sub-optimality bound for the MTS Problem 1. The bound is algorithm-agnostic (does not take assumption on the way ${{\cal O}}$ is generated), and is computable in $O(1)$ time. Also, we demonstrate its informativeness via simulations.

Theorem 6** (A posteriori sub-optimality bound)**

Consider the MTS problem (2) and let ${\cal O}^{\star}$ be an optimal solution to (2). Also, for any candidate solution ${\cal O}$, let:

• •

$r({\cal O})\doteq\min_{x\in\mathbb{X}}\;\;\sum_{i\in{\cal M}\setminus{\cal O}}\|h_{i}(y_{i},x)\|^{2}$; i.e., $r({\cal O})$ is the minimum residual error given the rejection ${\cal O}$;

• •

$r^{\star}_{k}\doteq\min_{{\cal O}\subseteq{\cal M},|{\cal O}|\leq k}r({\cal O})$; i.e., $r^{\star}_{k}$ is the optimal residual error when at most $k$ measurements are rejected;

• •

$r^{\star}\doteq r({\cal O}^{\star})$; i.e., $r^{\star}$ is the residual error for the optimal outlier rejection ${\cal O}^{\star}$.

Then, given a candidate solution ${\cal O}$, the following bound relates the residual error $r({\cal O})$ of the candidate solution with the residual error of an optimal solution rejecting the same number of outliers:

[TABLE]

where

[TABLE]

Moreover, if it is also known that $|{\cal O}|\;\geq|{\cal O}^{\star}|$, then it holds:

[TABLE]

$\lrcorner$**

Eq. (10) quantifies the distance between the residual of the candidate solution and the residual of an optimal solution rejecting the same number of outliers $|{\cal O}|$. Intuitively, if we incorrectly pick outliers and obtain a residual error $r({\cal O})$, there might exist a more clever selection that instead obtains $r^{\star}_{|{\cal O}|}\ll r({\cal O})$; on the other hand, we would like $r({\cal O})$ and $r^{\star}_{|{\cal O}|}$ to be as close as possible. For this reason, the smaller $\chi_{\cal O}$, the closer the candidate selection is to the optimal selection. For example, when $\chi_{\cal O}=0$, then $r({\cal O})=r^{\star}_{|{\cal O}|}$, i.e., we conclude that the algorithm returned a globally optimal solution (restricted to the ones rejecting $|{\cal O}|$ measurements).

Eq. (12) completes the picture by stating that if the algorithm rejects at least as many measurements as the optimal solution ($|{\cal O}|\;\geq|{\cal O}^{\star}|$), then the bound in eq. (12) compares the quality of ${\cal O}$ directly with the optimal residual error of ${\cal O}^{\star}$, the optimal solution of the MTS Problem 1.

Remark 7** (Quality of the bound)**

We showcase the quality of the bound (11) by considering the linear estimation Example 1. We generate small instances for which we can compute the optimal solution and evaluate the corresponding residual error $r^{\star}$. In particular, we compute the optimal solution using CPLEX [30], a popular library for mixed-integer linear programming, and we compare the optimal solution against a candidate solution ${\cal O}$. We generate candidate solutions using a greedy algorithm [29]. The greedy algorithm, at each iteration, rejects the measurement that induces the largest decrease in the residual error. Fig. 2 shows in blue the actual (true) approximation performance of the greedy algorithm ($\frac{r({\cal O})-r^{\star}}{r(\emptyset)-r^{\star}}$) and in red the bound $\chi_{{\cal O}}$ in eq. (10). The results are averaged over 10 Monte Carlo simulations. The figure shows that the bound predicts well the actual sub-optimality ratio (note the scale of the $y$-axis) and its quality improves for increasing number of outliers. $\lrcorner$

We remark that the bound (11) can be also used to quantify the performance of existing algorithms, including RANSAC. Having introduced our per-instance sub-optimality bounds, we step forward to a novel general-purpose algorithm for outlier rejection that empirically returns accurate solutions (and for which our bound $\chi_{\cal O}$ is typically close to zero).

V A General-Purpose Algorithm: ADAPT

We introduce a novel algorithm for outlier rejection that we name Adaptive Trimming (ADAPT). The algorithm starts by processing all measurements and at each iteration it trims measurements with residuals larger than a threshold. It is adaptive in that it dynamically decides the threshold at each iteration (hence relaxing the need for a threshold $\epsilon_{{\cal M}\setminus{\cal O}}$). Moreover, it is not greedy in that it can reject multiple measurements at each iteration while it keeps revisiting the quality of previously rejected outliers.222In our tests we found that a greedy algorithm similar to [29] tends to converge to poor outlier rejection decisions and is typically slow for practical applications, since it has quadratic runtime in the number of measurements.

Assumption 8** (Global solver)**

ADAPT* assumes the availability of a black-box solver that can (even approximately) solve the outlier-free problem (1) to optimality.*

Luckily, for all problems in Section II-A, there exist (outlier-free) global solvers, including [6, 7, 8].

Description of ADAPT. The preudo-code of ADAPT is given in Algorithm 1. Here, we use the additional notation:

• •

Let $x^{\star}({\cal O})\in\arg\min_{x\in\mathbb{X}}\;\;\sum_{i\in{\cal M}\setminus{\cal O}}\|h_{i}(y_{i},x)\|^{2}$; i.e., $x^{\star}({\cal O})$ is an estimator of $x$ given an outlier selection ${\cal O}$.

• •

Let $r_{i}({\cal O})\doteq\|h_{i}(y_{i},x^{\star}({\cal O}))\|^{2}$; i.e., $r_{i}({\cal O})$ is the residual of the measurement $i$, given an outlier selection ${\cal O}$.

Per Algorithm 1, ADAPT executes five distinctive operations:

Initialization (line 1)

ADAPT initializes the iteration counter $t=0$ and the current candidate outlier set ${\cal O}_{t}$ with the empty set. It also initializes $g$, the outlier group size, which constrains the maximum number of measurements that can be deemed as outliers in a single iteration of the algorithm. Moreover, ADAPT initializes the counter $c=0$: this is used to decide whether convergence has been reached. Finally, it initializes the outlier threshold $\tau$ with the value of the largest residual across all measurements. Note that computing $r_{i}({\cal O})$ (for any ${\cal O}\subseteq{\cal M}$) requires calling the global solver on the measurements ${\cal M}\setminus{\cal O}$, and then evaluating the residual errors for all measurements in ${\cal M}$.

Outlier set update (lines 5-6)

Given the current threshold $\tau$ and group size $g$, the algorithm updates the outlier set in two steps: first (line 5), it finds the set of measurements ${\cal I}$ with the $g$ largest residuals among all measurements in ${\cal M}$;333Note that the selection is performed over all measurements ${\cal M}$, potentially revisiting measurements that were previously deemed to be outliers. second (line 6), the algorithm updates the outlier set as the collection of all the measurements in ${\cal I}$ whose residual exceeds the outlier threshold $\tau$.

Outlier threshold update (lines 7-8)

If the updated outlier set ${\cal O}_{t}$ remains the same as that in the previous iteration ${\cal O}_{t-1}$ (line 7), then the outlier threshold $\tau$ is not tight enough. As a result, the algorithm updates $\tau$ with a discounted value $\gamma<1$ (line 8). This process is repeated as long as necessary, as indicated by the “while” loop in line 4.

Outlier group size update (line 9)

After each iteration $t$, ADAPT increases the outlier group size $g$ by $\bar{g}$. This has the effect of increasing the maximum number of measurements that can be deemed as outliers in future iterations: intuitively, ADAPT is conservative in rejecting measurements at the beginning (small initial $g=\bar{g}$) , while it gets more and more aggressive by gradually increasing $g$.

Termination

ADAPT terminates when one of the following two conditions is satisfied. First (lines 10-11), it may terminate when all but a number $v$ of measurements have been rejected, where $v$ is the minimum number of measurement that the global solver needs to solve the problem (for example, in 3D registration, $v=3$). Second (lines 14-15), ADAPT may terminate if convergence has been achieved. In particular, if the algorithm observes for $T$ consecutive times that the absolute value of the residuals function changes by less than $\delta$, then it terminates ($c$ counts the number of consecutive times a decrease smaller than $\delta$ is observed).

Remark 9** (Complexity and practicality)**

The termination condition in line 10 guarantees the termination of the algorithm with at most $|{\cal M}|-v$ calls of the global solver. ADAPT terminates faster as one increases the outlier group size $\bar{g}$, the convergence thresholds $\delta$, and/or as one decreases the discount factor $\gamma$ and the number $T$ of iterations to decide convergence. Overall, the linear runtime (in the number of measurements) of ADAPT makes the algorithm practical in real-time applications where fast global solvers are available.

Remark 10** (vs. RANSAC)**

While RANSAC builds an inlier set by sampling small (minimal) sets of measurements, ADAPT iteratively prunes the overall set of measurements. Arguably, this gives ADAPT a “global vision” of the measurement set as we showcase in the experimental section. RANSAC assumes the availability of fast minimal solvers, while ADAPT assumes the availability of fast global (non-minimal) solvers. Finally, RANSAC is not suitable for high-dimensional problems where is becomes exponentially more difficult to sample an outlier-free set [18]. On the other hand, ADAPT is deterministic and guaranteed to terminate in a number of iterations bounded by the number of measurements.

VI Experiments and Applications

We evaluate ADAPT against the state of the art in three spatial perception problems: 3D registration (Section VI-A), two-view geometry (Section VI-B), and SLAM (Section VI-C). The results show that ADAPT outperforms RANSAC in terms of accuracy and scalability, and often outperforms specialized outlier rejection methods (in particular for SLAM) while being a general-purpose algorithm. Finally, the tests show that the performance bounds of Section IV are informative and can be used to assess the outlier rejection outcomes. All results are averaged over 10 Monte Carlo runs.

VI-A Robust Registration

Experimental setup. We test ADAPT on two standard datasets for 3D registration: the Stanford Bunny and the ETH Hauptgebaude [31]. In both cases we downsample the point clouds obtaining 453 points for Bunny and 3617 points for ETH. For each point cloud ${\cal P}$ we generate a second point cloud ${\cal P}^{\prime}$ by applying a random rigid transformation and adding noise and outliers. The (inlier) noise standard deviation is set to 0.025% and 0.05% of the point cloud diameter respectively. Outliers are generated by replacing a subset of the points in ${\cal P}^{\prime}$ with random points uniformly sampled in the bounding box containing ${\cal P}$. In each iteration, ADAPT uses Horn’s method [8] as global solver. We benchmark ADAPT against Fast Global Registration (FGR) [25] and the three-point RANSAC. We set the maximum number of iterations in RANSAC to 1000 and use default parameters for FGR. All methods are implemented in MATLAB.

Results. Fig. 3 shows the (average) translation and rotation errors for the estimates computed by ADAPT, FGR, and RANSAC on the Bunny dataset for increasing outlier percentages. ADAPT performs on-pair with FGR which is a specialized robust solver for 3D registration and they both achieve practically zero error for up to 90% of outliers, after which they both break. RANSAC starts performing distinctively worse early on and is dominated by the other techniques (after 90% all techniques fail to provide a satisfactory estimate). We obtain similar results on the ETH dataset hence for space reasons we report them in Appendix C.

For both the Bunny and ETH datasets, we compute the sub-optimality bound for the result of ADAPT, using Theorem 6. The plot of the bound is given in Appendix C; the bound remains around $10^{-5}$, confirming that ADAPT remains close to the optimal outlier selection. The runtime of ADAPT is comparable to FGR and is reported in Appendix C.

VI-B Robust Two-view Geometry

Experimental setup. We tested ADAPT on both synthetic data and on the MH_01 sequence of the EuRoC dataset [32]. To generate the synthetic data we place the first camera at the origin (identity pose) and place the second camera randomly within a bounded region. Then, we generate a random point cloud within the field-of-view of both camera. The points projected on the camera frame are finally corrupted with Gaussian noise, and outliers are added. For the EuRoC dataset, we extract Harris corners in each frame and track keypoints in consecutive frames using Lucas-Kanade feature tracking as in [33]; the results are then averaged across all pairs of consecutive frames in the sequence. In each iteration, ADAPT uses Briales’ QCQP relaxation [7] as global solver. We benchmarked ADAPT against Nister’s five-point [26] and the eight-points algorithm [34] within RANSAC.

Results. Fig. 4 shows the box-plot of translation and rotation errors for the estimates of ADAPT, the five- and eight-point RANSAC on the synthetic dataset. ADAPT outperforms the other techniques across all the spectrum. ADAPT and five-point perform on-pair till 40% of outliers. Beyond that point, the five-point method attains considerably higher errors than ADAPT (50% to 100% more in rotation; and more than 300% more in translation). The eight-point method results in higher errors than the five-point across the spectrum.

Fig. 5 shows the results on the EuRoC dataset focusing on the comparison between ADAPT and the five-point RANSAC. ADAPT achieves a mean rotation error of $2.5\cdot 10^{-3}$rad versus $2.8\cdot 10^{-3}$rad of the five-point. Similarly for the translation error: 0.075m for ADAPT versus 0.09m for the five-point. For visualization purposes we cut the translation box-plot in Fig. 5 above 0.25m error: in reference to the rest of the plot, we report that ADAPT exhibits translation errors larger than $1$ in 10% of the frames, whereas only $1\%$ of the five-point estimates have translation errors larger than $1$.

For the synthetic dataset, the typical value for the sub-optimality bound achieved by ADAPT is $0.2$. That is, ADAPT makes a rejection that achieves an error that is at most 20% away from the optimal, even in the presence of 90% of outliers. The runtime of ADAPT is reported in Appendix C: our approach is one order of magnitude slower than the five-point method, mainly due to the relatively high runtime of the global solver [7], which is called in each iteration.

VI-C Robust SLAM

Experimental setup. We test ADAPT on standard 2D and 3D SLAM benchmarking datasets and report extra results on synthetic datasets in Appendix C. We spoil existing datasets with spurious loop closures: we sample random pairs of nodes and we add an outlier relative pose measurement between them, where the relative translation is sampled in the ball of radius 5m and the rotation is sampled uniformly at random in $\mathrm{SO}(2)$ or $\mathrm{SO}(3)$. The ground truth trajectory is generated by optimizing the problem with SE-Sync [6] before adding outliers. We also use SE-Sync as the global solver for ADAPT. We test the following datasets, described in [35, 6]: MIT (2D), Intel (2D), CSAIL (2D), and Sphere2500 (3D). We also test a simulated $5\times 5\times 5$ 3D grid dataset (results in Appendix C). We benchmark ADAPT against DCS [36]; we report DCS results for three choices of the robust kernel size: {1,10,100} (the default value is 1, see [36]).

Results. ADAPT outperforms DCS (independently on the choice of the kernel size) across all outlier percentages.

2D SLAM: In the MIT dataset (Fig. 6a), a particularly challenging dataset, ADAPT is insensitive to up to 20% of outliers. All variants of DCS fail to produce an error smaller than 10 meters even in the absence of outliers: DCS is an iterative solver, hence is may converge to local minima when bootstrapped from a bad initial guess (the odometric estimate is particularly noisy in MIT). ADAPT leverages SE-sync, which is a global solver, hence is able to converge to the correct solution. In the CSAIL dataset (Fig. 6b), ADAPT also dominates DCS while DCS performance is acceptable when the kernel size is equal to 100. Unfortunately, the choice of a “good” kernel size is dataset-dependent and it can make the difference between an adequate performance in DCS and a catastrophic failure. ADAPT also dominates DCS in the INTEL dataset, whose statistics are reported in Appendix C.

3D SLAM: In the Sphere dataset (Fig. 7), ADAPT achieves and error below 0.25m across all outlier percentages. On the other hand, DCS starts with an error of 0.25 meters (0% outliers) and ends up with at least 2.5 meters error (50% outliers). Again, ADAPT, a general-purpose approach for outlier rejection, outperforms specialized techniques for robust SLAM. Appendix C also reports similar conclusions and results for the 3D grid dataset.

For 2D SLAM, the typical value for the sub-optimality bound achieved by ADAPT, per Theorem 6, is 0.1 (MIT) and 0.01 (Intel and CSAIL) across all spectrum of outlier percentages. For 3D SLAM, the typical value of the bound achieved by ADAPT is 0.1 (3D grid) and 0.01 (sphere).

ADAPT is one to two orders slower than DCS. This is due to the repeated calls to SE-sync and is further aggravated in the 3D case by the fact that SE-sync tends to be slow in the presence of outliers (the Riemannian staircase method [6] requires multiple staircase iterations since the rank of the relaxation increases in the presence of outliers).

VII Related Work

We extend the literature review in Section I, across robotics and computer vision (Section VII-A), and statistics and control (Section VII-B), as well as, across sub-optimality guarantees in outlier-robust estimation (Section VII-C).

VII-A Outlier-robust estimation in robotics and computer vision

Outlier-robust estimation is an active research area in robotics and computer vision for decades [37, 38, 9]. Traditionally, low-dimensional problems are solved using RANSAC [11] (e.g., 3D registration, and two-view geometry). RANSAC is efficient and effective for low outlier rates [10]. However, the need to cope with higher outlier rates pushed research towards global optimization methods, such as branch-and-bound (BnB) [14, 15, 16, 17, 10], and mixed-integer programming (MIP) [39]. Nevertheless, for high-dimensional problems (e.g., SLAM), RANSAC, BnB, and MIP are typically slow to be practical: in the worst-case, they require exponential run-time. Hence, research has also focused on M-estimators, in conjunction with either non-convex optimization [40, 36, 41], or convex relaxations [42], possibly also including decision variables for the outlier rejection [43, 42, 40, 36]. In more detail, representative outlier-robust methods for 3D registration, two-view geometry, and SLAM are the following.

Robust 3D registration. The goal is to find the transformation (translation and rotation) that aligns two point clouds. In the presence of outliers, created by incorrect point correspondences, one typically resorts to RANSAC [44, 11], along with a 3-point minimal solver —the outlier-free problem admits well-known closed form solutions [45, 8]. But when the outliers’ number is more than 50%, RANSAC tends to be slow, and brittle [10, 18]. Thereby, recent approaches adopt either robust cost functions [46, 9, 25], or BnB [24, 47, 18]: Zhou et al. [25] propose fast global registration (FGR), which is based on the Geman-McClure robust cost function; Yang et al. [24] propose a approach; Campbell et al. [47] employ BnB to search the space of camera poses, guaranteeing global optimality without requiring a pose prior; and Bustos et al. [18] add a pre-processing step, that removes gross outliers before RANSAC or BnB. Other approaches, that iteratively compute point correspondences, include iterative closest point (ICP) [48, 49], and trimmed iterative closest point algorithm [50]; all require an accurate initial guess [25].

Robust two-view geometry. The problem consists in estimating the relative pose (up to scale) between two images, given pixel correspondences. In robotics, RANSAC is again the go-to approach [1], typically in conjunction with Nister’s 5-point method [51], a minimal solver —other minimal solvers exist when one is given a reference direction [52], the relative rotation [53], or motion constraints [54]. In computer vision, recent approaches are investigating the use of provably-robust techniques, typically based on BnB: Hartley et al. [15] propose a BnB approach for $\ell_{\infty}$ optimization in one-view and two-view geometry; Li [17] uses BnB and mixed-integer programming for two-view geometry; Bazin et al. [14] use BnB for rotation-only estimation; Chin et al. [39] propose a method to remove outliers in conjunction with mixed-integer linear programming; Zheng et al. [16] use BnB to estimate the fundamental matrix; and Speciale et al. [10] improve BnB approaches by including linear matrix inequalities. BnB is typically slow [17], but it is able to tolerate high outlier rates.

Robust SLAM. Outlier-robust SLAM has traditionally relied on M-estimators; e.g., [9]. Olson and Agarwal [41] use a max-mixture distribution to approximate multi-modal measurement noise. Sünderhauf and Protzel [40, 55] augment the problem with latent binary variables responsible for deactivating outliers. Tong and Barfoot [56, 57] propose algorithms to classify outliers via chi-square statistical tests that account for the effect of noise in the estimate. Latif et al. [12] propose realizing, reversing, and recovering (RRR), which performs loop-closure outlier rejection, by clustering measurements together and checking for consistency using the chi-square inverse test as an outlier-free bound. A pair-wise consistency check for multi-robot SLAM, pair-wise consistency maximization (PCM), was proposed by Mangelson et al. [13]. Agarwal et al. [36] propose dynamic covariance scaling (DCS), which adjusts the measurement covariances to reduce the influence of outliers. Lee et al. [58] use expectation maximization. The papers above rely either on the availability of an initial guess for optimization, or on parameter tuning. Recent work also includes convex relaxations for outlier-robust SLAM [59, 60, 61, 42]. Currently, only [42] provides so far sub-optimality guarantees, which however degrade with the quality of the relaxation. Additionally, [42] requires parameter tuning.

VII-B Outlier-robust estimation in statistics and control

Outlier-robust estimation has received long-time attention in statistics and control [21, 62]. It has a fundamental applications, such as prediction and learning [63], linear decoding [64], and secure state estimation for control [65].

Outlier-robust estimation in statistics, in its simplest form, aims to learn the mean and covariance of an unknown distribution, given both a portion of noiseless i.i.d. samples, and a portion of arbitrarily corrupted samples (outliers); particularly, the number of outliers is assumed known.444In contrast, eq. (1)’s framework (and, correspondingly, MTS’s) considers the estimation of an unknown parameter given some measurements, which framework is equivalent to the estimation of the parameter given non-i.i.d. samples (the samples depend on the unknown parameter). Then, researchers provide polynomial time near-optimal algorithms [63, 66].

In scenarios where one aims to estimate an unknown parameter given noisy, and possibly outlying, measurements, then researchers focus on outlier-robust reformulations of eq. (1). For example, Rousseeuw [23] proposed a celebrated algorithm to solve a dual reformulation of MTS, that aims to minimize the residual errors of the remaining measurements given a maximum number of measurement rejections. The algorithm assumes the outliers’ number known, and as such, requires parameter tuning. Similar celebrated algorithms, that also assume the outliers’ number known, are the forward greedy by Nemhauser et al. [29], and forward-backward greedy by Zhang [67] (notably, both algorithms have quadratic run-time, a typically prohibitive run-time for robotics and computer vision applications, such as 3D registration, two-view, and SLAM). In contrast to [23, 29, 67], the greedy-like algorithm proposed in [68] considers the outliers’ number unknown. However, it still requires parameter tuning, this time for an outlier-free bound parameter.

Outlier-robust estimation in control takes typically the form of secure state estimation in the presence of outliers (caused by sensor malfunctions, or measurement attacks) [65, 69, 70, 71]. In [65, 69, 70, 71], the researchers propose optimal algorithms, that achieve exact state estimation when the non-outlying measurements are noiseless. However, the algorithms have exponential run-time.

VII-C Sub-optimality guarantees in outlier-robust estimation

Additionally to the aforementioned algorithms that offer sub-optimality gurantees (or certificates of optimality), researchers have also provided conditions for exact estimation when some of the measurements are outliers, while the rest are noiseless [64, 67]. However, the conditions’ evaluation is NP-hard. Additionally, the conditions are restricted on a linear and convex framework, where the measurements are linear in the unknown parameter $x$, and $x\in\mathbb{R}^{n}$.555In contrast, Section IV provides a posteriori sub-optimality bounds, computable in $O(1)$ run-time, and applicable to even non-convex and non-linear frameworks, such as for 3D registration, two-view geometry, and SLAM.

VIII Conclusion

We proposed a minimally trimmed squares (MTS) formulation to estimate an unknown variable from measurements plagued with outliers. We proved that the resulting outlier rejection problem is inapproximable: one cannot compute even an approximate solution in quasi-polynomial time. We derived theoretical performance bounds: while polynomial-time algorithms may perform poorly in the worst-case, the bounds allow assessing the algorithms’ post-run performance on any given problem instances (which are typically more favorable than the worst-case). Finally, we proposed a linear-time, general-purpose algorithm for outlier rejection, and showed that it outperforms several specialized methods across three spatial perception problems (3D registration, two-view geometry, SLAM). This work paves the way for several research avenues. While we focused on a non-linear least squares cost function, many of our conclusions extend to other norms, and robust costs. We also plan to explore applications of the proposed bounds to other algorithms, including RANSAC, and to other perception problems.

Appendix A Proof of Theorem 5

Here, we show the inapproximability of MTS by reducing it to the variable selection problem, which we define next.

Problem 2** (Variable Selection)**

Let $U\in\mathbb{R}^{m\times l}$, $z\in\mathbb{R}^{m}$, and let $\Delta$ be a non-negative number. The variable selection problem asks to pick $d\in\mathbb{R}^{l}$ that is an optimal solution to the following optimization problem:

[TABLE]

$\lrcorner$**

Variable selection is inapproximable in quasi-polynomial time. We summarise the result in Lemm 13 below. To this end, we first review basic definitions from complexity theory.

Definition 11** (Big O notation)**

Let $\mathbb{N}_{+}$ be the set of non-negative natural numbers, and consider two functions $h:\mathbb{N}_{+}\mapsto\mathbb{R}$ and $g:\mathbb{N}_{+}\mapsto\mathbb{R}$. The big O notation in the equality $h(n)=O(g(n))$ means there exists some constant $c>0$ such that for all large enough $n$, $h(n)\leq cg(n)$. $\lrcorner$

That is, $O(g(n))$ denotes the collection of functions $h$ that are bounded asymptotically by $g$, up to a constant factor.

Definition 12** (Big $\Omega$ notation)**

Consider two functions $h:\mathbb{N}_{+}\mapsto\mathbb{R}$ and $g:\mathbb{N}_{+}\mapsto\mathbb{R}$. The big $\Omega$ notation in the equality $h(n)=\Omega(g(n))$ means there exists some constant $c>0$ such that for all large enough $n$, $h(n)\geq cg(n)$. $\lrcorner$

That is, $\Omega(g(n))$ denotes the collection of functions $h$ that are lower bounded asymptotically by $g$, up to a constant.

Lemma 13** ([72, Proposition 6])**

For each $\delta\in(0,1)$, unless it is ${\rm NP}{\in}{\rm BPTIME(\textit{l}^{\rm poly\log\textit{l}})}$, there exist:

• •

a function $q_{1}(l)$ which is in $2^{\Omega(\log^{1-\delta}l)}$;

• •

a polynomial $p_{1}(l)$ which is in $O(l)$;666In this context, a function with a fractional exponent is considered to be a polynomial, e.g., $l^{1/5}$ is considered to be a polynomial in $l$.

• •

a polynomial $\Delta(l)$;

• •

a polynomial $m(l)$,

and a zero-one $m(l)\times l$ matrix $U$ such that even if it is known that a solution to $Ud=1_{m(l)}$ exists, no quasi-polynomial algorithm can distinguish between the next cases for large $l$:

1. $S_{1}$)

There exists a vector $d\in{\mathbb{R}^{l}}$ such that $Ud=\mathbf{1}_{m(l)}$ and $||d||_{0}\leq p_{1}(l)$. 2. $S_{2}$)

For any vector $d\in\mathbb{R}^{l}$ such that $||Ud-\mathbf{1}_{m(l)}||_{2}^{2}\leq\Delta(l)$, we have $||d||_{0}\geq p_{1}(l)q_{1}(l)$.

Unless ${\rm NP}{\in}{\rm BPTIME(\textit{l}^{\rm poly\log\textit{l}})}$, Theorem 2 says that variable selection is inapproximable even in quasi-polynomial time. This is in the sense that for large $l$ there is no quasi-polynomial algorithm that can distinguish between the two mutually exclusive statements $S_{1}$ and $S_{2}$. These statements are indeed mutually exclusive for large $l$, since then $q_{1}(l)>1$, since it is $q_{1}(l)=2^{\Omega\left(\log^{1-\delta}l\right)}$.

From the inapproximability of variable selection, we can now infer the inapproximability for the following problem:

[TABLE]

Proof that problem in eq. (13) is inapproximable Indeed, it suffices to set $\Delta=0$ in the definition of the variable selection problem, and then consider Lemma 13. $\blacksquare$

Next, from the inapproximability of the problem in eq. (13), we next infer the inapproximability for the problem below:

[TABLE]

To this end, consider the instance of Lemma 13, and let: $\Delta^{\prime}(l)=m^{2}(l)\Delta(l)$; $y$ be any solution to $Uy=1_{m(l)}$ (per Lemma 13 we know there exists a solution to this equation); and $A$ be a matrix in $\mathbb{R}^{l\times n}$, where $n=l-\text{rank}(U)$,777By this construction, t is $l>n$. That is, $A$ is a tall matrix with more rows than columns. such that the columns of $A$ span the null space of $U$ (hence, $A$ is such that $UA=0$). This instance of the problem in eq. (14) is constructed in polynomial time in $l$, since solving a system of equations (as well as finding eigenvectors that span the null space of a matrix) happens in polynomial time.

Given the above instance of the problem in eq. (14), we next prove that the following two statements are indistinguishable to prove that the problem is inapproximable:

1. $S_{1}^{\prime}$)

There exist vectors $d\in{\mathbb{R}^{l}}$ and $x\in{\mathbb{R}^{n}}$ such that $y=Ax+d$ and $||d||_{0}\leq p_{1}(l)$. 2. $S_{2}^{\prime}$)

For any vectors $d\in{\mathbb{R}^{l}}$ and $x\in{\mathbb{R}^{n}}$ such that $||y-Ax-d||^{2}\leq\Delta^{\prime}(l)$,888The norm in $S_{2}^{\prime}$ (namely, the $||y-Ax-d||^{2}$) can be any norm that is polynomially close (in $l$) to $\ell_{2}$-norm, such as the $\ell_{1}$-norm. we have $||d||_{0}\geq p_{1}(l)q_{1}(l)$.

Proof that $S_{1}^{\prime}$ and $S_{2}^{\prime}$ are indistinguishable: We prove that whenever statements $S_{1}$ and $S_{2}$ in Theorem 13 are true, then also statements $S_{1}^{\prime}$ and $S_{2}^{\prime}$ above are true, respectively. That is, all true instances of $S_{1}$ and $S_{2}$ are mapped to true instances of $S_{1}^{\prime}$ and $S_{2}^{\prime}$. Then, since also the mapping is done in polynomial time, this implies that no algorithm can solve the problem in eq. (14) in quasi-polynomial time and distinguish the cases $S_{1}^{\prime}$ and $S_{2}^{\prime}$, because that would contradict that $S_{1}$ and $S_{2}$ are indistinguishable.

Proof that when $S_{1}$ is true then also $S_{1}^{\prime}$ also is

Since $Uy=UAx+Ud$ implies $1_{m(l)}=Ud$, when $S_{1}$ is true, then $S_{1}^{\prime}$ is also true with $x$ being the unique solution of $Ax=y-d$ (it is unique since $A$ is full column rank).

Proof that when $S_{2}$ is true then also $S_{2}^{\prime}$ also is

By contradiction: Assume that there are vectors $d\in{\mathbb{R}^{l}}$ and $x\in{\mathbb{R}^{n}}$ such that $||y-Ax-d||^{2}\leq\Delta^{\prime}(l)$ and $||d||_{0}<p_{1}(l)q_{1}(l)$. Without loss of generality, assume $\ell_{1}$ norm for $||y-Ax-d||^{2}$. Then, $||y-Ax-d||^{2}_{1}\leq\Delta^{\prime}(l)$ implies $||U||^{2}_{1}||y-Ax-d||^{2}_{1}\leq||U||^{2}_{1}\Delta^{\prime}(l)$, which implies $||U(y-Ax-d)||^{2}_{1}\leq||U||^{2}_{1}\Delta^{\prime}(l)$, or $||1_{m(l)}-Ud||^{2}\leq||U||^{2}\Delta^{\prime}(l)$, or $||1_{m(l)}-Ud||^{2}\leq m^{2}(l)\Delta^{\prime}(l)$, where the last holds true because $U$ is a zero-one matrix. Finally, because by definition of $\Delta^{\prime}(l)$ it is $m^{2}(l)\Delta^{\prime}(l)=\Delta(l)$, we have $||1_{m(l)}-Ud||^{2}\leq\Delta(l)$. Overall, we just proved that there exist $d$ such that $||1_{m(l)}-Ud||^{2}\leq\Delta(l)$ and $||d||_{0}<p_{1}(l)q_{1}(l)$, which contradicts $S_{2}$. $\blacksquare$

The final step for the proof of Theorem 5 is to reduce MTS to the problem in eq. (14). To this end, per the statement of Theorem 5 we focus on the linear framework of Example 1. Also, we use the following notation:

• •

Let $y_{{\cal M}\setminus{\cal O}}$ be the vector-stack of measurements $y_{i}$ such that $i\in{\cal M}\setminus{\cal O}$, given a measurement rejection set ${\cal O}$.

• •

Similarly, let $d_{{\cal M}\setminus{\cal O}}$ be the vector-stack of noises $d_{i}$ such that $i\in{\cal M}\setminus{\cal O}$.

• •

And similarly, let $A_{{\cal O}}$ be the matrix-stack of measurement row-vectors $a_{i}^{\top}$.

Given the above notation, we can now write MTS per the framework of Example 1 as follows:

[TABLE]

To prove now the inapproximability of the MTS problem in eq. (15), consider an instance of Lemma 13, along with a corresponding instance considered above for the inapproximability of the problem in (14). Then, set $\epsilon=\Delta^{\prime}(l)$ in eq. (15). We prove that the following two statements are indistinguishable:

1. $S_{1}^{\prime\prime}$)

There exist ${{\cal O}}\in\{1,\ldots,l\}$ and $x\in{\mathbb{R}^{n}}$ such that $y_{{{\cal M}}\setminus{{\cal O}}}=A_{{{\cal M}}\setminus{{\cal O}}}x$ and $|{{\cal O}}|\leq p_{1}(l)$. 2. $S_{2}^{\prime\prime}$)

For any ${{\cal O}}\in\{1,\ldots,l\}$ and $x\in{\mathbb{R}^{n}}$ such that $||y_{{{\cal M}}\setminus{{\cal O}}}-A_{{{\cal M}}\setminus{{\cal O}}}x||^{2}\leq\Delta^{\prime}(l)$,999The norm in $S_{2}^{\prime\prime}$ can be any norm that is polynomially close (in $l$) to $\ell_{2}$-norm, such as the $\ell_{1}$-norm. we have $|{{\cal O}}|\geq p_{1}(l)q_{1}(l)$.

Proof that $S_{1}^{\prime\prime}$ and $S_{2}^{\prime\prime}$ are indistinguishable We prove that whenever statements $S_{1}^{\prime}$ and $S_{2}^{\prime}$ above are true, then also statements $S_{1}^{\prime\prime}$ and $S_{2}^{\prime\prime}$ above are true, respectively.

Proof that when $S_{1}^{\prime}$ is true then also $S_{1}^{\prime\prime}$ also is

Assume $S_{1}$ is true. Let ${{\cal O}}=\{\text{all }i\text{ such that }d_{i}\neq 0,i\in\{1,\ldots,l\}\}$. Then, $y_{{{\cal M}}\setminus{{\cal O}}}=A_{{{\cal M}}\setminus{{\cal O}}}x$ since $d_{{{\cal M}}\setminus{{\cal O}}}=0$, and $|{{\cal O}}|=||d||_{0}\leq p_{1}(l)$.

Proof that when $S_{2}^{\prime}$ is true then also $S_{2}^{\prime\prime}$ also is

By contradiction: Assume there are ${{\cal O}}\in\{1,\ldots,l\}$ and $x\in{\mathbb{R}^{n}}$ such that $||y_{{{\cal M}}\setminus{{\cal O}}}-A_{{{\cal M}}\setminus{{\cal O}}}x||^{2}\leq\Delta^{\prime}(l)$ and $|{{\cal O}}|<p_{1}(l)q_{1}(l)$. Let $d_{{{\cal M}}\setminus{{\cal O}}}=0$, and $d_{{\cal O}}=y_{{\cal O}}-A_{{\cal O}}x$. Then, $||d||_{0}=|S|<p_{1}(l)q_{1}(l)$ and $||y-Ax-d||^{2}=||y_{{{\cal M}}\setminus{{\cal O}}}-A_{{{\cal M}}\setminus{{\cal O}}}x||^{2}\leq\Delta^{\prime}(l)$, which contradicts $S_{2}^{\prime}$. $\blacksquare$

Overall, we found instances for MTS, per Example 1, and for a number $l$ of measurements, where the statements $S_{1}^{\prime\prime}$ and $S_{2}^{\prime\prime}$ are indistinguishable even in quasi-polynomial time, and as a result, the proof of Theorem 5 is now complete.

Appendix B Proof of Theorem 6

First observe that $r^{\star}_{|{\cal O}|}\leq r({\cal O})$. This holds true due to the definition of $r^{\star}_{|{\cal O}|}$ as the smallest value of $r(\cdot)$ among all sets with cardinality $|{\cal O}|$. Next, define the quantities:

• •

$f({\cal O})\doteq r(\emptyset)-r({\cal O})$;

• •

$f^{\star}\doteq r(\emptyset)-r^{\star}_{|{\cal O}|}$.

We observe that:

[TABLE]

The above now imples:

[TABLE]

where the latter holds since $f({\cal O})\leq f^{\star}$, which in turns holds because $r({\cal O})\geq r^{\star}$. Finally, the above is:

[TABLE]

which gives:

[TABLE]

which gives eq. (10).

To prove eq. (12), we first observe that $r^{\star}\geq r^{\star}_{|{\cal O}|}$, since $|{\cal O}|\geq|{\cal O}^{\star}|$. Now, it can be verified that:

[TABLE]

since it also is $r(\emptyset)\geq r({\cal O})$. Substituting eq. (16) in eq. (10), the proof of the theorem is now complete.

Appendix C Supplemental for Experiments and Applications

We first provide ADAPT’s selected input parameters for each experiment, along with the methodology we applied to select them. Then, we list all missing plots from the experimental section of the paper (Experiments and Applications).

For our experiments, we selected the input parameters for ADAPT by testing the algorithm on each dataset on scenarios where we added 20% of outliers to the dataset. We did not add outliers to the EuRoc dataset. instead, we used the parameters selected for the synthetic two-view dataset. In particular, the parameters $g$ and $\gamma$ were chosen to make the algorithm terminate with a few iterations (e.g., around 20 iterations for SLAM) —we recall that $g$ (outlier group size) and $\gamma$ (outlier-threshold discount factor) control the maximum number of rejected measurements at each iteration. The parameter $\delta$ (convergence threshold) was chosen so the algorithm rejects an accurate number of measurements upon termination. Specifically, 20% of measurements, since 20% was the percentage of the added outliers by us to the datatset. We kept $T$ (number of iterations to decide convergence) fixed to its default value 2 for all the experiments. Overall, the selected parameters for each experiment are shown in Table I. Finally, ADAPT requires the minimum number of measurements required by global solver. This number for registration is 3; for two-view is 5; and for SLAM is 0.

We list all missing plots from the experimental section of the paper in the following pages.

C-A Supplemental for Robust Registration

C-B Supplemental for Robust Two-view Geometry

C-C Supplemental for Robust SLAM