An Efficient Solution to Non-Minimal Case Essential Matrix Estimation

TL;DR

This paper introduces a globally optimal, efficient solver for the N-point problem in essential matrix estimation, improving speed and robustness over existing methods for relative pose estimation in computer vision.

Contribution

The paper proposes a certifiably globally optimal solver for the N-point problem formulated as a QCQP, with theoretical guarantees and a robust extension using M-estimators.

Findings

Solver is 2-3 orders of magnitude faster than state-of-the-art methods.

Robust N-point method outperforms existing methods in robustness and accuracy.

Theoretical guarantees include tightness and local stability of the relaxation.

Abstract

Finding relative pose between two calibrated images is a fundamental task in computer vision. Given five point correspondences, the classical five-point methods can be used to calculate the essential matrix efficiently. For the case of $N$ ($N > 5$) inlier point correspondences, which is called $N$-point problem, existing methods are either inefficient or prone to local minima. In this paper, we propose a certifiably globally optimal and efficient solver for the $N$-point problem. First we formulate the problem as a quadratically constrained quadratic program (QCQP). Then a certifiably globally optimal solution to this problem is obtained by semidefinite relaxation. This allows us to obtain certifiably globally optimal solutions to the original non-convex QCQPs in polynomial time. The theoretical guarantees of the semidefinite relaxation are also provided, including tightness and local stability. To deal with outliers, we propose a robust $N$-point method using M-estimators. Though global optimality cannot be guaranteed for the overall robust framework, the proposed robust $N$-point method can achieve good performance when the outlier ratio is not high. Extensive experiments on synthetic and real-world datasets demonstrated that our $N$-point method is $2\sim3$ orders of magnitude faster than state-of-the-art methods. Moreover, our robust $N$-point method outperforms state-of-the-art methods in terms of robustness and accuracy.

Figures

Tables

Table 1. TABLE I: Efficiency comparison with other globally optimal methods. The last column is the normalized runtime by setting ours as 1 1 1 .

method	optimization	runtime	norm. runtime
Hartley & Kahl [18]	BnB	$>7$ s	$>$ $1000$
Chesi [8]	LMI	$1.15$ s	$190$
SDP-Briales [10]	SDP	$1$ s	$160$
ours	SDP	$6$ ms	$1$

Table 2. TABLE II: SDP formulations comparison. For domain 𝒮 n superscript 𝒮 𝑛 \mathcal{S}^{n} , the number of variable is n ​ ( n + 1 ) / 2 𝑛 𝑛 1 2 n(n+1)/2 .

method	domain	#variable	#constraint
SDP-Briales [10]	${𝒮}^{40}$	$820$	$536$
ours	${𝒮}^{12}$	$78$	$7$

Equations

$${\mathcal{S}}^2\triangleq \{\mathbf{t}\in {\mathbb{R}}^3|{\mathbf{t}}^{\top }\mathbf{t}=1\}.$$

$$\text{SO}(3)\triangleq \{\mathbf{R}\in {\mathbb{R}}^{3\times 3}|{\mathbf{R}}^{\top }\mathbf{R}=\mathbf{I},\text{det}(\mathbf{R})=1\},$$

$$\mathbf{E}=[\mathbf{t}{]}_{\times }\mathbf{R},$$

$$[\mathbf{t}{]}_{\times }={\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]}_{\times }=\left[\begin{array}{ccc}0& -t_3& t_2\\ t_3& 0& -t_1\\ -t_2& t_1& 0\end{array}\right].$$

$$\mathbf{E}=\left[\begin{array}{c}{\mathbf{e}}_1^{\top }\\ {\mathbf{e}}_2^{\top }\\ {\mathbf{e}}_3^{\top }\end{array}\right]=\left[\begin{array}{ccc}{e}_{11}& e_{12}& e_{13}\\ e_{21}& e_{22}& e_{23}\\ e_{31}& e_{32}& e_{33}\end{array}\right].$$

$$\mathbf{e}$$

$$=[e_{11};e_{12};e_{13};e_{21};e_{22};e_{23};e_{31};e_{32};e_{33}].$$

$${\mathcal{M}}_{\mathbf{E}}\triangleq \{\mathbf{E}|\mathbf{E}=[\mathbf{t}{]}_{\times }\mathbf{R},\exists \mathbf{R}\in \text{SO}(3),\mathbf{t}\in {\mathcal{S}}^2\}.$$

$$\text{(i)}\mathbf{E}{\mathbf{E}}^{\top }=[\mathbf{t}{]}_{\times }[\mathbf{t}{]}_{\times }^{\top }\text{and}\text{(ii)}{\mathbf{t}}^{\top }\mathbf{t}=1.$$

$${\mathbf{f}}_i^{\top }\mathbf{E}{\mathbf{f}}_i^{\prime }=0.$$

$$\underset{\mathbf{E}\in {\mathcal{M}}_{\mathbf{E}}}{\min }$$

$$\sum _{i=1}^N({\mathbf{f}}_i^{\top }\mathbf{E}{\mathbf{f}}_i^{\prime }{)}^2={\mathbf{e}}^{\top }\mathbf{Ce},$$

$$\mathbf{C}=\sum _{i=1}^N\left({\mathbf{f}}_i\otimes {\mathbf{f}}_i^{\prime }\right){\left({\mathbf{f}}_i\otimes {\mathbf{f}}_i^{\prime }\right)}^{\top },$$

$$\underset{\mathbf{E},\mathbf{R},\mathbf{t}}{\min }$$

$$\mathbf{E}=[\mathbf{t}{]}_{\times }\mathbf{R},\mathbf{R}\in \text{SO}(3),\mathbf{t}\in {\mathcal{S}}^2,$$

$$\underset{\mathbf{E},\mathbf{t}}{\min }$$

$$\mathbf{E}{\mathbf{E}}^{\top }=[\mathbf{t}{]}_{\times }[\mathbf{t}{]}_{\times }^{\top },{\mathbf{t}}^{\top }\mathbf{t}=1.$$

$$h_1={\mathbf{e}}_1^{\top }{\mathbf{e}}_1-(t_2^2+t_3^2)=0,$$

$$h_2={\mathbf{e}}_2^{\top }{\mathbf{e}}_2-(t_1^2+t_3^2)=0,$$

$$h_3={\mathbf{e}}_3^{\top }{\mathbf{e}}_3-(t_1^2+t_2^2)=0,$$

$$h_4={\mathbf{e}}_1^{\top }{\mathbf{e}}_2+t_1{t}_2=0,$$

$$h_5={\mathbf{e}}_1^{\top }{\mathbf{e}}_3+t_1{t}_3=0,$$

$$h_6={\mathbf{e}}_2^{\top }{\mathbf{e}}_3+t_2{t}_3=0,$$

$$h_7={\mathbf{t}}^{\top }\mathbf{t}-1=0.$$

$$\mathbf{C}=\sum _{i=1}^N{w}_i\left({\mathbf{f}}_i\otimes {\mathbf{f}}_i^{\prime }\right){\left({\mathbf{f}}_i\otimes {\mathbf{f}}_i^{\prime }\right)}^{\top }.$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }$$

$${\mathbf{x}}^{\top }{\mathbf{A}}_i\mathbf{x}=b_i,i=1,\cdots ,m.$$

$$\mathbf{x}\triangleq [\mathbf{e};\mathbf{t}]$$

$$\{\begin{array}{l}{\mathbf{x}}^{\top }{\mathbf{C}}_0\mathbf{x}=\mathrm{trace}({\mathbf{x}}^{\top }{\mathbf{C}}_0\mathbf{x})=\mathrm{trace}({\mathbf{C}}_0\mathbf{x}{\mathbf{x}}^{\top }),\\ {\mathbf{x}}^{\top }{\mathbf{A}}_i\mathbf{x}=\mathrm{trace}({\mathbf{x}}^{\top }{\mathbf{A}}_i\mathbf{x})=\mathrm{trace}({\mathbf{A}}_i\mathbf{x}{\mathbf{x}}^{\top }).\end{array}$$

$$\underset{\mathbf{X}\in {\mathcal{S}}^n}{\min }$$

$$\mathrm{trace}({\mathbf{A}}_i\mathbf{X})=b_i,i=1,\cdots ,m,$$

$$\mathbf{X}⪰0,\mathrm{rank}(\mathbf{X})=1.$$

$$\underset{\mathbf{X}\in {\mathcal{S}}^n}{\min }$$

$$\mathrm{trace}({\mathbf{A}}_i\mathbf{X})=b_i,i=1,\cdots ,m,$$

$$\mathbf{X}⪰0.$$

$$\underset{\mathbf{\lambda }}{\max }$$

$$\mathbf{Q}(\mathbf{\lambda })={\mathbf{C}}_0-\sum _{i=1}^m{\lambda }_i{\mathbf{A}}_i⪰0,$$

$${\mathbf{E}}^{\star }=\mathrm{mat}({\mathbf{e}}^{\star },[3,3]),$$

Full text

1 Introduction

FINDNG relative pose between two images using 2D-2D point correspondences is a cornerstone in geometric vision. It makes a basic building block in many structure-from-motion (SfM), visual odometry, and simultaneous localization and mapping (SLAM) systems [1]. Relative pose estimation is a difficult problem since it is by nature non-convex and known to be plagued by local minima and ambiguous solutions. Relative pose for uncalibrated and calibrated cameras are usually characterized by fundamental matrix and essential matrix, respectively [1]. A matrix is a fundamental matrix if and only if it has two non-zero singular values. An essential matrix has an additional property that the two nonzero singular values are equal. Due to these strict constraints, essential matrix estimation is arguably thought to be more challenging than fundamental matrix estimation. This paper focuses on optimal essential matrix estimation.

Due to the scale ambiguity of translation, relative pose for calibrated cameras has $5$ degrees-of-freedom (DoFs), including $3$ for rotation and $2$ for translation. Except for degenerate configurations, $5$ point correspondences are hence enough to determine the relative pose. Given five point correspondences, the five-point methods using essential matrix [2, 3] or rotation matrix parametrization [4] can be used to calculate the relative pose efficiently. The aforementioned solvers are the so-called minimal solvers. When point correspondences contain outliers, minimal solvers are usually integrated into a hypothesize-and-test framework, such as RANSAC [5], to find the solution corresponding to the maximal consensus set. Hence this framework can provide high robustness.

Once the maximal consensus set has been found, the standard RANSAC optionally re-estimates a model by using all inliers to reduce the influence of noise [1, 6]. Thus a non-minimal solver is needed in this procedure. This RANSAC framework is called the gold standard algorithm [1]. Figure 1 illustrates this framework by taking line fitting as an example. In addition to the usage as post-processing, the non-minimal solvers can also be integrated more tightly into RANSAC variants. For example, LO-RANSAC [7] attempts to enlarge the consensus set of an initial RANSAC estimate by generating hypotheses from larger-than-minimal subsets of the consensus set. The rationale is that hypotheses fitted on a larger number of inliers typically lead to better estimates with higher support.

In this paper, the non-minimal case relative pose estimation is called $N$-point problem, and its solver is called $N$-point method. As it has been investigated in [8, 9, 10], $N$-point methods usually lead to more accurate results than five-point methods. Thus the $N$-point method is useful for scenarios that require accurate pose estimation, such as visual odometry and image-based visual servoing. The well-known direct linear transformation (DLT) technique [1] with proper normalization [11] can be used to estimate the essential matrix using $8$ or more point correspondences. However, DLT ignores the inherent nonlinear constraints of the essential matrix. To deal with this problem, an essential matrix is recovered after an approximated essential matrix is obtained from the DLT solution [1]. An eigenvalue-based formulation and its variant were proposed to solve $N$-point problem [9, 10]. However, all the aforementioned methods fail to guarantee global optimality and efficiency simultaneously.

Since the $N$-point problem is challenging, the progress of its solvers is far behind the absolute pose estimation (perspective-$n$-point (P$n$P)) and point cloud registration. For example, the EP$n$P algorithm  [12] in P$n$P area is very efficient and has linear complexity with the number of observations. It has been successfully used in RANSAC framework given an arbitrary number of 2D-3D correspondences. Though relative pose estimation is more difficult than absolute pose estimation, it is desirable to find a practical solver to $N$-point problem, whose efficiency and global optimality are both satisfactory. This is one motivation of this paper.

Another motivation of this paper is developing an M-estimator based method for relative pose estimation. In practice, it arises frequently that the data have been contaminated by large noise and outliers. Existing robust estimation methods are mainly classified into two main categories, i.e., inlier set maximization [13] and M-estimator based method [14]. Inlier set maximization can be achieved by randomized sampling methods [5, 7] or deterministic optimization methods [15, 16]. For M-estimator based methods, the associated optimization problems are non-convex and difficult to solve. In this paper, the proposed robust $N$-point method uses M-estimators. To solve the associated optimization problem effectively, the line process is adopted [17]. The line process uses a continuous iterative optimization strategy which takes a weighted version of non-minimal solver as a vital requirement. Due to the lack of an efficient and globally optimal $N$-point method, there did not exist a practical M-estimator based relative estimation method before.

Based on the aforementioned motivations, in this paper we propose a novel $N$-point method and integrate this method into M-estimators. The contributions of this paper are three-fold.

• •

Efficient and globally optimal $N$-point method. A simple parameterization is proposed to characterize the essential manifold. Based on this parameterization, a certifiably globally optimal $N$-point method is proposed, which is $2\sim 3$ orders of magnitude faster than state-of-the-art methods.

• •

Robust $N$-point method. We propose a robust essential matrix estimation method by integrating $N$-point method into M-estimators. Considering that the robust components of the overall framework are not certifiably and provably optimal, we can only demonstrated empirical performance assurances.

• •

Theoretical aspects. We provide theoretical proofs of the semidefinite relaxation (SDR) in the proposed $N$-point method, including SDR tightness and local stability with small observation noise.

The paper is organized as follows. Section 2 introduces the related work. In Section 3, we propose a simple parameterization of essential manifold and provide novel formulations of $N$-point methods. Based on these formulations, Section 4 derives a convex optimization approach by SDR. Section 5 proves tightness and local stability of SDR. A robust $N$-point method based on robust loss function is proposed in Section 6. Section 7 presents the performance of our method in comparison to other approaches, followed by a concluding discussion in Section 8.

2 Related Work

Estimating an essential matrix for a calibrated camera from point correspondences is an active research area in computer vision. Finding the optimal essential matrix by $L_{\infty}$ norm cost and branch-and-bound (BnB) was proposed in [18]. It achieves global optima but is inefficient. There are several works for $N$-point fundamental/essential matrix estimation using local optimization [19, 20] or manifold optimization [21, 22, 23]. A method for minimizing an algebraic error was investigated in [8]. Its global optimality was obtained by a square matrix representation of homogeneous forms and relaxation. An eigenvalue-based formulation was proposed to estimate rotation matrix [9], in which the problem is optimized by local gradient descent or BnB search. It was later improved by a certifiably globally optimal solution by relaxation and semidefinite programming (SDP) [10]. However, none of the aforementioned methods can find the globally optimal solution efficiently. Though BnB search methods [18, 9] can obtain global optima in theory, they have the exponential time complexity in the worst case. In [8, 10], a theoretical guarantee of the convexification procedures is not provided. The most related paper to this work is [10], which converts the $N$-point problem of an eigenvalue-based formulation to a QCQP. However, its efficiency is not satisfactory and tightness of its SDR has not been proved.

There are also several works on fundamental matrix estimation for uncalibrated cameras. The eight-point method [11] uses a linear solution, then recovers a valid fundamental matrix by SVD decomposition. This method ignores the rank constraint in the fundamental matrix, thus the solution is not optimal. In [24], a method for minimizing an algebraic error was proposed which ensures the rank constraint. However, it does not guarantee global minima. In [25, 26], the constraint for a fundamental matrix is imposed by setting its determinant as [math], leading to a cubic polynomial constraint. In [27], the fundamental matrix estimation problem is reduced to one or several constrained polynomial optimization problems. Unfortunately, the aforementioned methods deal with uncalibrated cameras only, where the underlying Euclidean constraints of an essential matrix are not exploited. Thus they cannot be applied to essential matrix estimation.

For both the essential matrix and fundamental matrix, optimal pose estimation can be formulated as a polynomial optimization problem [28]. A polynomial optimization problem can be converted to a QCQP. In multiple view geometry, SDR for polynomial optimization problems was first studied in [29]. Later, a large number of methods using QCQP formulations were developed in computer vision and robotics. For example, SDR or Lagrangian duality of QCQPs has been used in point set registration [30], triangulation [31], relative pose estimation [10], rotation averaging [32], pose synchronization [33], and the Wahba Problem [34]. However, SDR does not guarantee a priori that it generates an optimal solution. If a QCQP satisfies certain conditions and data noise lies within a critical threshold, a recent study proved that the solution to the SDP optimization algorithm is guaranteed to be globally optimal [35, 36]. A noise threshold that guarantees tightness of SDR is given for the rotation averaging problem [32]. For general QCQPs, the global optimality still remains an open problem.

Existing robust estimation methods in geometric vision are mainly based on inlier set maximization [13] or M-estimators [14]. Inlier set maximization was proven to be NP-hard [37]. BnB search can be used to find the globally optimal solution [38, 39, 40, 41, 42], but its efficiency is not satisfactory. There are a variety of methods to approximately and efficiently solve the inlier set maximization problem. The most popular algorithms belong to a class of randomized sampling techniques, i.e., RANSAC [5] and its variants [7, 6][15][16]. A hybrid method of BnB and mixed integer programming (MIP) [43] was proposed to solve the inlier set maximization in relative pose estimation. Another alternative robust framework, which is based on M-estimators, has been successfully applied to many fields such as bundle adjustment [44], registration [45, 46], and data clustering [47]. An important technique to optimize M-estimators is the line process [17], which is also a building block of the robust version of the proposed method. In the line process, progress is hindered by a lack of an efficient and globally optimal non-minimal solver. The proposed $N$-point method in this paper can be integrated into the line process to make a robust $N$-point method.

3 Formulations of $N$-Point Method

Denote $(\mathbf{p}_{i},\mathbf{p}^{\prime}_{i})$ as the $i$-th point correspondence of the same 3D world point from two distinct viewpoints. Point observations $\mathbf{p}_{i}$ and $\mathbf{p}^{\prime}_{i}$ are represented as homogeneous coordinates in normalized image plane111Bold capital letters denote matrices (e.g., $\mathbf{E}$ and $\mathbf{R}$); bold lower-case letters denote column vectors (e.g., $\mathbf{e},\mathbf{t}$); non-bold lower-case letters represent scalars (e.g., $\lambda$). By Matlab syntax, we use semicolon/comma in matrix concatenation to arrange entries vertically/horizontally. For example, $[[a],[b]]=[a,b]$ and $[[a];[b]]=\begin{bmatrix}a\\ b\end{bmatrix}$.. Each point in the normalized image plane can be translated into a unique unit bearing vector originating from the camera center. Let ($\mathbf{f}_{i}$, $\mathbf{f}^{\prime}_{i}$) denote a correspondence of bearing vectors pointing at the same 3D world point from two distinct viewpoints, where $\mathbf{f}_{i}$ represents the observation from the first viewpoint, and $\mathbf{f}^{\prime}_{i}$ the observation from the second viewpoint. The bearing vectors are determined by $\mathbf{f}_{i}=\frac{\mathbf{p}_{i}}{\|\mathbf{p}_{i}\|}$ and $\mathbf{f}^{\prime}_{i}=\frac{\mathbf{p}^{\prime}_{i}}{\|\mathbf{p}^{\prime}_{i}\|}$.

The relative pose is composed of rotation $\mathbf{R}$ and translation $\mathbf{t}$. Rotation $\mathbf{R}$ transforms vectors from the second into the first frame. Translation $\mathbf{t}$ is expressed in the first frame and denotes the position of the second frame with respect to the first one. The normalized translation $\mathbf{t}=[t_{1},t_{2},t_{3}]^{\top}$ will be identified with points in the $2$-sphere $\mathcal{S}^{2}$, i.e.,

[TABLE]

The 3D rotation will be featured as $3\times 3$ orthogonal matrix with positive determinant belonging to the special orthogonal group $\text{SO}(3)$, i.e.,

[TABLE]

where $\mathbf{I}$ is a $3\times 3$ identity matrix.

3.1 Parametrization for Essential Manifold

The essential matrix $\mathbf{E}$ is defined as [1]

[TABLE]

where $[\cdot]_{\times}$ defines the corresponding skew-symmetric matrix for a $3$-dimensional vector, i.e.,

[TABLE]

Denote the essential matrix $\mathbf{E}$ as

[TABLE]

where $\mathbf{e}_{i}^{\top}$ is the $i$-th row of $\mathbf{E}$. Denote its corresponding vector as

[TABLE]

where $\operatorname*{vec}(\cdot)$ means stacking all the entries of a matrix by row-first order. There is not any essential difference between different orders of entry stacking in this paper. We adopt row-first order to make the notations in Section 5 more convenient.

In this paper, an essential matrix set is defined as

[TABLE]

This essential matrix set is called normalized essential manifold [21, 22, 23]. Theorem 1 provides equivalent conditions to define $\mathcal{M}_{\mathbf{E}}$, which will greatly simplify the optimization in the proposed methods.

Theorem 1**.**

A real $3\times 3$ matrix, $\mathbf{E}$, is an element in $\mathcal{M}_{\mathbf{E}}$ if and only if there exists a vector $\mathbf{t}\in\mathbb{R}^{3}$ satisfying the following two conditions:

[TABLE]

Proof.

For if direction, first it can be verified that $\text{det}\left([\mathbf{t}]_{\times}[\mathbf{t}]_{\times}^{\top}-\sigma\mathbf{I}\right)=-\sigma[\sigma-(t_{1}^{2}+t_{2}^{2}+t_{3}^{2})]^{2}=-\sigma(\sigma-1)^{2}$. By combining this result with condition (i), we can see that $\mathbf{E}\mathbf{E}^{\top}$ has an eigenvalue $1$ with multiplicity $2$ and an eigenvalue [math]. According to the definition of singular value, the nonzero singular values of $\mathbf{E}$ are the square roots of the nonzero eigenvalues of $\mathbf{E}\mathbf{E}^{\top}$. Thus the two nonzero singular values of $\mathbf{E}$ are equal to $1$. According to Theorem 1 in [48], $\mathbf{E}$ is an essential matrix. By combining condition (ii), $\mathbf{E}$ is an element in $\mathcal{M}_{\mathbf{E}}$.

For only if direction, $\mathbf{E}$ is supposed to be an essential matrix $\mathcal{M}_{\mathbf{E}}$. According to the definition of $\mathcal{M}_{\mathbf{E}}$, there exists a vector $\mathbf{t}$ satisfying condition (ii). In addition, there exists a rotation matrix $\mathbf{R}$ such that $\mathbf{E}=[\mathbf{t}]_{\times}\mathbf{R}$. It can be verified that $\mathbf{E}\mathbf{E}^{\top}=([\mathbf{t}]_{\times}\mathbf{R})([\mathbf{t}]_{\times}\mathbf{R})^{\top}=[\mathbf{t}]_{\times}[\mathbf{t}]_{\times}^{\top}$, thus condition (i) is also satisfied. ∎

It is worth mentioning that a necessary condition for general essential matrix, which is similar to the only if direction in Theorem 1, was presented in [48, Proposition 2] and [49, Lemma 7.2]. In Theorem 1, we further prove that this condition is also sufficient and propose a novel parameterization for the normalized essential manifold.

3.2 Optimizing Essential Matrix by Minimizing an Algebraic Error

For noise-free cases, the epipolar constraint [1] implies that

[TABLE]

Under the presence of noise, this constraint does not strictly hold. We pursue the optimal pose by minimizing an algebraic error

[TABLE]

The algebraic error in objective has been widely used in previous literature [50, 24, 27, 8].

The objective in problem (7) can be reformulated as a standard quadratic form

[TABLE]

where

[TABLE]

and “$\otimes$” means Kronecker product. Note that $\mathbf{C}$ is a Gram matrix, so it is positive semidefinite and symmetric.

3.3 QCQP Formulations

By explicitly writing the constraints for the essential manifold $\mathcal{M}_{\mathbf{E}}$, we reformulate problem (7) as

[TABLE]

This problem is a QCQP: The objective is positive semidefinite quadratic polynomials; the constraint on the translation vector, $\mathbf{t}^{\top}\mathbf{t}=1$, is also quadratic; a rotation matrix $\mathbf{R}$ can be fully defined by $20$ quadratic constraints [4, 10]; and the relationship between $\mathbf{E}$, $\mathbf{R}$ and $\mathbf{t}$, $\mathbf{E}=[\mathbf{t}]_{\times}\mathbf{R}$, is also quadratic. This formulation has $21$ variables and $30$ constraints.

According to Theorem 1, an equivalent QCQP form of minimizing the algebraic error is

[TABLE]

There are $12$ variables and $7$ constraints in this problem. The constraints can be written explicitly as below

[TABLE]

Problems (10) and (11) are equivalent since their objectives are the same and their feasible regions are equivalent. Both of the two problems are nonconvex. Appendix A provides another equivalent optimization problem. In the following, we will only consider problem (11) due to its fewer variables and constraints. Moreover, it is homogeneous without the need of homogenization, which makes it simpler than the alternative formulations.

Remark: Both problems (10) and (11) are equivalent to an eigenvalue-based formulation [9, 10]. A proof of the equivalence is available in [10], see its supplementary material. Since all the mentioned formulations essentially utilize the normalized essential manifold as feasible regions and have the equivalent objectives, all these formulations are equivalent. Our formulations (10) and (11) have the following two advantages:

(1) Our formulations have fewer variables and constraints. In contrast, the eigenvalue based formulation in [10] involves $40$ variables and $536$ constraints. As shown in the following sections, the simplicity of our formulations will result in much more efficient solvers and enable the proof of tightness and local stability.

(2) Our formulations are easy to associate priors for each point correspondence by simply introducing weights in the objective. For example, we may introduce a weight for each sample by slightly changing the objective tp $\sum_{i=1}^{N}w_{i}(\mathbf{f}_{i}^{\top}\mathbf{E}\mathbf{f}^{\prime}_{i})^{2}$, where $w_{i}\geq 0$ is the weight for $i$-th observation. For general cases in which $w_{i}\geq 0$, we can keep current formulation by changing the construction of $\mathbf{C}$ to

[TABLE]

4 Optimization of $N$-Point Method

QCQP is a long-standing problem in optimization literature with many applications. Solving its general case is an NP-hard problem. Global optimization methods for QCQP are typically based on convex relaxations of the problem. There are two main relaxations for QCQP: SDR and the reformulation-linearization technique. In this paper, we use SDR since it usually has better performance [51] and it is convenient for tightness analysis.

Let us consider a QCQP in a general form as

[TABLE]

where $\mathbf{C}_{0}$, $\mathbf{A}_{1},\cdots,\mathbf{A}_{m}\in\mathcal{S}^{n}$ and $\mathcal{S}^{n}$ denotes the set of all real symmetric $n\times n$ matrices. In our problem,

[TABLE]

is a vector stacking all entries in essential matrix $\mathbf{E}$ and translation vector $\mathbf{t}$; $n=12$; $m=7$; $\mathbf{C}_{0}=\begin{bmatrix}\mathbf{C}&\mathbf{0}_{9\times 3}\\ \mathbf{0}_{3\times 9}&\mathbf{0}_{3\times 3}\end{bmatrix}$; $\mathbf{A}_{1},\cdots,\mathbf{A}_{7}$ correspond to the canonical form $\mathbf{x}^{\top}\mathbf{A}_{i}\mathbf{x}$ of Eqs. (12a)$\sim$(12g), respectively.

A crucial first step in deriving an SDR of problem (14) is to observe that

[TABLE]

It can be seen that both the objective and constraints in problem (14) are linear in matrix $\mathbf{x}\mathbf{x}^{\top}$. Thus, by introducing a new variable $\mathbf{X}=\mathbf{x}\mathbf{x}^{\top}$ and noting that $\mathbf{X}=\mathbf{x}\mathbf{x}^{\top}$ is equivalent to $\mathbf{X}$ being a rank one symmetric positive semidefinite (PSD) matrix, we obtain the following equivalent formulation of problem (14)

[TABLE]

Here, $\mathbf{X}\succeq\mathbf{0}$ means that $\mathbf{X}$ is PSD. Solving rank constrained semidefinite programs (SDPs) is NP-hard [52]. SDR drops the rank constraint to obtain the following relaxed version of problem (17)

[TABLE]

Problem (18) turns out to be an instance of SDP [52], which belongs to convex optimization and can be readily solved using primal-dual interior point methods [53]. Its dual problem is

[TABLE]

where $\mathbf{b}=[b_{1},\cdots,b_{m}]^{\top}$, $\boldsymbol{\lambda}=[\lambda_{1},\cdots,\lambda_{m}]^{\top}\in\mathbb{R}^{m}$. Problem (19) is called the Lagrangian dual problem of problem (14), and $\mathbf{Q}(\boldsymbol{\lambda})$ is the Hessian of the Lagrangian. In our problem, $b_{i}=0$ for $i=1,\cdots,6$; $b_{7}=1$; and $\mathbf{b}^{\top}\boldsymbol{\lambda}=\lambda_{7}$.

In summary, the relations between different formulations are demonstrated by Fig. 2.

4.1 Essential Matrix and Relative Pose Recovery

Once the optimal $\mathbf{X}^{\star}$ of the SDP primal problem (18) has been obtained by an SDP solver, we need to recover the optimal essential matrix $\mathbf{E}^{\star}$. Denote $\mathbf{X}_{e}$ as the top-left $9\times 9$ submatrix of $\mathbf{X}$; and denote $\mathbf{X}_{t}$ as the bottom-right $3\times 3$ submatrix of $\mathbf{X}$, i.e., $\mathbf{X}_{e}\triangleq\mathbf{X}_{[1:9,1:9]}$ and $\mathbf{X}_{t}\triangleq\mathbf{X}_{[10:12,10:12]}$. Empirically, we found that the largest singular value of $\mathbf{X}_{e}^{\star}$ is near $2$ and others are close to zero. It is common to set the rank of a matrix as the number of singular values larger than a threshold. In our method, the threshold depends on the accuracy of SDP solver (usually $10^{-7}\sim 10^{-5}$), leading to $\operatorname*{rank}(\mathbf{X}_{e}^{\star})=1$. Denote the eigenvector that corresponding to the nonzero eigenvalue of $\mathbf{X}_{e}^{\star}$ as $\mathbf{e}^{\star}$, then the optimal essential matrix is recovered by

[TABLE]

where $\operatorname*{mat}(\mathbf{e},[r,c])$ means reshape the vector $\mathbf{e}$ to an $r\times c$ matrix by row-first order.

After the essential matrix has been obtained, we can recover the rotation and translation by the standard method in literature [1]. A recent work proved that the rotation matrix can be accurately recovered from the essential matrix for pure rotation scenarios [54]. Moreover, a statistic, the mean of $\left\{\frac{\mathbf{x}_{i}\times\mathbf{R}^{\star}\mathbf{x}^{\prime}_{i}}{\|\mathbf{x}_{i}\|\|\mathbf{x}^{\prime}_{i}\|}\right\}_{i=1}^{N}$, was proposed to identify the pure rotation scenarios.

In Section 4.2, the theoretical guarantee of such a pose recovery method will be provided. In Section 5, the proof of tightness and local stability that guarantees the global optimality will be provided. The outline of $N$-point method is shown in Algorithm 1.

4.2 Necessary and Sufficient Conditions for Global Optimality

The following Theorem 2 provides a theoretical guarantee for the proposed pose recovery method.

Theorem 2**.**

For QCQP (11), its SDR is tight if and only if: the optimal solution $\mathbf{X}^{\star}$ to its primal SDP problem (18) satisfies $\operatorname*{rank}(\mathbf{X}_{e}^{\star})=\operatorname*{rank}(\mathbf{X}_{t}^{\star})=1$.

Proof.

First, we prove the if part. Note that $\mathbf{X}_{e}^{\star}$ and $\mathbf{X}_{t}^{\star}$ are real symmetric matrices because they are in the feasible region of the primal SDP. In addition, it is given that $\operatorname*{rank}(\mathbf{X}_{e}^{\star})=\operatorname*{rank}(\mathbf{X}_{t}^{\star})=1$, thus there exist two vectors $\mathbf{e}^{\star}$ and $\mathbf{t}^{\star}$ satisfying $\mathbf{e}^{\star}({\mathbf{e}^{\star}})^{\top}=\mathbf{X}_{e}^{\star}$ and $\mathbf{t}^{\star}(\mathbf{t}^{\star})^{\top}=\mathbf{X}_{t}^{\star}$. Since the constraints in problem (11) do not include any cross term between $\mathbf{E}$ and $\mathbf{t}$, the intersection part of $\mathbf{E}$ and $\mathbf{t}$ in matrix $\mathbf{A}_{i}$ is zero in SDP problem. By substituting $\mathbf{e}^{\star}$ and $\mathbf{t}^{\star}$ into Eq. (16), it can be verified that $\mathbf{e}^{\star}$ and $\mathbf{t}^{\star}$ satisfy the constraints in primal problem (11). Now we can see that $\mathbf{X}^{\star}$ and its uniquely determined derivatives ($\mathbf{e}^{\star}$ and $\mathbf{t}^{\star}$) are feasible solutions for the semidefinite relaxation problem and the primal problem, respectively. Thus the relaxation is tight.

Then we prove the only if part. Since the semidefinite relaxation is tight, we have $\operatorname*{rank}(\mathbf{X}^{\star})=1$. Then $\operatorname*{rank}(\mathbf{X}_{e}^{\star})\leq 1$ and $\operatorname*{rank}(\mathbf{X}_{t}^{\star})\leq 1$. Since $\mathbf{X}_{e}^{\star}$ and $\mathbf{X}_{t}^{\star}$ cannot be zero matrices (otherwise $\mathbf{X}^{\star}$ is not in the feasible region), the equalities should hold, i.e., $\operatorname*{rank}(\mathbf{X}_{e}^{\star})=\operatorname*{rank}(\mathbf{X}_{t}^{\star})=1$. ∎

Theorem 2 provides a necessary and sufficient condition to recover the globally optimal solution for the primal problem. Empirically, the optimal $\mathbf{X}^{\star}$ by the SDP problem always satisfies this condition. Specifically, the optimal $\mathbf{X}^{\star}$ has the following structure

[TABLE]

The $\Diamond$ parts could be arbitrary matrices making $\mathbf{X}^{\star}$ symmetric.

Remark: The block diagonal structure of $\mathbf{X}^{\star}$ in Eq. (21) is caused by the sparsity pattern of the problem. The aggregate sparsity pattern in our SDP problem, which is the union of individual sparsity patterns of the data matrices, $\{\mathbf{C}_{0},\mathbf{A}_{1},\cdots,\mathbf{A}_{m}\}$, includes two cliques: one includes the $1\sim 9$-th entries of $\mathbf{x}$, and the other includes the $10\sim 12$-th entries of $\mathbf{x}$, see Fig. 3(a). There is no common node in these two cliques, see Fig. 3(b). The chordal decomposition theory of the sparse SDPs can explain the structure of $\mathbf{X}^{\star}$ well. The interested reader may refer to [55, 56] for more details.

4.3 Time Complexity

First, we consider the time complexity of problem construction. The construction of the optimization problem (i.e., calculating $\mathbf{C}$ in Eq. (9) or Eq. (13)) is linear with the number of point correspondences, so its time complexity is $\mathcal{O}(N)$. It is worth mentioning that optimization is independent of the number of point correspondences once the problem has been constructed.

Second, we discuss the time complexity of SDP optimization. Most SDP solvers use an interior-point algorithm. The SDP problem (18) can be solved with a worst case complexity of

[TABLE]

flops given a solution accuracy $\epsilon>0$ [53]. It can be seen that the time complexity can be largely reduced given a smaller variable number $n$ and constraint number $m$. Since our formulations have much fewer variables and constraints than previous work [10], they own much lower time complexity.

Finally, we discuss the time complexity of pose recovery. In our method, the essential matrix is recovered by finding the eigenvector corresponding to the largest eigenvalue of a $9\times 9$ matrix $\mathbf{X}_{e}^{\star}$. Thus the time complexity of pose recovery is $\mathcal{O}(1)$.

5 Tightness and Local Stability of $N$-Point Method

In this section, we prove tightness and local stability of the SDR for our problem. The readers who are not interested in theory can safely skip this section. To understand the proofs in this section, preliminary knowledge about convex optimization, manifold, and algebraic geometry is necessary. We recommend the readers to refer to [35, Chapter 6][36] for more details.

In the following, the bar on a symbol stands for a value under the noise-free case. In other words, it represents ground truth. For example, $\bar{\mathbf{C}}$ denotes the matrix in objective constructed by noise-free observations; $\bar{\mathbf{x}}$ is the optimal state estimation from noise-free observations; $\bar{\mathbf{t}}=[\bar{t}_{1},\bar{t}_{2},\bar{t}_{3}]^{\top}$ is the optimal translation estimation from noise-free observations.

5.1 Tightness of the Semidefinite Relaxation

In this subsection, we prove that the SDR is tight given noise-free observations. The proof is mainly based on Lemma 2.4 in [36].

Lemma 1**.**

If point correspondences $\{(\mathbf{f}_{i},\mathbf{f}^{\prime}_{i})\}_{i=1}^{N}$ are noise-free, the matrix $\bar{\mathbf{C}}$ in Eq. (9) satisfies that $\operatorname*{rank}(\bar{\mathbf{C}})\leq\min(N,8)$, where $N$ is the number of point correspondences. The equality holds except for degenerate configurations including points on a ruled quadratic, points on a plane, and no translation (explanation of these degeneracies can be found in [57]).

Proof.

(i) When $N\leq 8$, from the construction of $\mathbf{C}$ in Eq. (9), each point correspondence adds a rank-$1$ matrix to the Gram matrix $\mathbf{C}$. The rank-$1$ matrices are linear independent except for the degenerate cases [57]. Thus the rank of $\mathbf{C}$ is $N$ for non-degenerate cases. When a degeneracy occurs, there will be linear dependence between these rank-1 matrices, and the rank will be below $N$. (ii) When $N>8$, there exists the stacked vector $\bar{\mathbf{e}}$ of an essential matrix satisfying that $\bar{\mathbf{e}}^{\top}\bar{\mathbf{C}}\bar{\mathbf{e}}=0$, so the upper limit of $\operatorname*{rank}(\bar{\mathbf{C}})$ is $8$. Then we complete the proof by combining these two properties. ∎

Given a set of point correspondences, if $N\geq 8$ after excluding points that belong to planar degeneracy, the rank of $\bar{\mathbf{C}}$ is $8$. This principle is the basis of the eight-point method [11]. In our methods, we do not need to distinguish the points that belong to the degenerate configurations, so the rank-$8$ assumption in the following text can be easily satisfied in $N$-point problem in which $N\gg 8$.

Lemma 2**.**

Let $\mathbf{C}\in\mathbb{R}^{n\times n}$ be positive semidefinite. If $\mathbf{x}^{\top}\mathbf{C}\mathbf{x}=0$ for a given vector $\mathbf{x}$, then $\mathbf{C}\mathbf{x}=\mathbf{0}$.

Proof.

Since $\mathbf{C}$ is positive semidefinite, its eigenvalues are non-negative. Suppose the rank of $\mathbf{C}$ is $r$. The eigenvalues can be listed as $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{r}>0=\sigma_{r+1}=\cdots=\sigma_{n}$. Denote the $i$-th eigenvector as $\mathbf{x}_{i}$. The eigenvectors are orthogonal to each other, i.e., $\mathbf{x}_{i}^{\top}\mathbf{x}_{j}=0$ when $i\neq j$. Vector $\mathbf{x}$ can be expressed as $\mathbf{x}=\sum_{i=1}^{n}a_{i}\mathbf{x}_{i}$. Thus, $\mathbf{C}\mathbf{x}=\mathbf{C}\sum_{i=1}^{n}a_{i}\mathbf{x}_{i}=\sum_{i=1}^{n}a_{i}\sigma_{i}\mathbf{x}_{i}$, and $\mathbf{x}^{\top}\mathbf{C}\mathbf{x}=\sum_{i=1}^{n}a_{i}\mathbf{x}_{i}^{\top}\cdot\sum_{i=1}^{n}a_{i}\sigma_{i}\mathbf{x}_{i}=\sum_{i=1}^{r}\sigma_{i}a_{i}^{2}\|\mathbf{x}_{i}\|^{2}$. Given $\mathbf{x}^{\top}\mathbf{C}\mathbf{x}=0$, we have $a_{i}=0$ for $i=1,\cdots,r$. Thus $\mathbf{C}\mathbf{x}=\sum_{i=1}^{n}a_{i}\sigma_{i}\mathbf{x}_{i}=\mathbf{0}$ is obtained. ∎

Lemma 3**.**

If point correspondences $\{(\mathbf{f}_{i},\mathbf{f}^{\prime}_{i})\}_{i=1}^{N}$ are noise-free, there is zero-duality-gap between problem (11) and its Lagrangian dual problem (19).

Proof.

Our proof is an application of the Lemma 2.4 in [36]. Let $\bar{\mathbf{x}}=[\bar{\mathbf{e}};\bar{\mathbf{t}}]$ be a feasible point in the primal problem, where $\bar{\mathbf{e}}$ and $\bar{\mathbf{t}}$ are ground truth. And let $\boldsymbol{\lambda}=\mathbf{0}$ be a feasible point in its Lagrangian dual problem (19). The three conditions needed in Lemma 2.4 in [36] are satisfied: (i) Primal feasibility. By substituting $\bar{\mathbf{x}}$ in the primal problem, the constraints are satisfied since $\bar{\mathbf{x}}$ is ground truth and the point correspondences are noise-free. (ii) Dual feasibility. $\mathbf{Q}(\boldsymbol{\lambda})=\mathbf{C}_{0}-\sum_{i=1}^{m}\lambda_{i}\mathbf{A}_{i}=\begin{bmatrix}\bar{\mathbf{C}}&\mathbf{0}_{9\times 3}\\ \mathbf{0}_{3\times 9}&\mathbf{0}_{3\times 3}\end{bmatrix}\succeq 0$. (iii) Lagrangian multiplier. It satisfies that $(\mathbf{C}_{0}-\sum_{i=1}^{m}\lambda_{i}\mathbf{A}_{i})\bar{\mathbf{x}}=\begin{bmatrix}\bar{\mathbf{C}}\bar{\mathbf{e}}\\ \mathbf{0}_{3\times 1}\end{bmatrix}$. Since $\bar{\mathbf{e}}$ is the ground truth, $\bar{\mathbf{e}}^{\top}\mathbf{C}\bar{\mathbf{e}}=0$ according to Eq. (8). Recall that $\bar{\mathbf{C}}$ is a Gram matrix and thus it is positive semidefinite. According to Lemma 2, $\bar{\mathbf{C}}\bar{\mathbf{e}}=\mathbf{0}$ is obtained. ∎

The zero-duality gap still holds for the case of noisy observations. A proof is provided in Appendix B.

5.2 Local Stability of the Semidefinite Relaxation

In this subsection, we prove that the SDR has local stability near noise-free observations. In other words, our QCQP formulation has a zero-duality-gap regime when its observations are perturbed (e.g., with noise in the case of sensor measurements). The proof is based on Theorem 5.1 in [36]. Following [36], we will use the following notations in the remains of this section to make notation simplicity.

• •

$\bar{\theta}\in\Theta$ is a zero-duality-gap parameter. In our problem, $\bar{\theta}=\{\bar{\mathbf{C}}\}$.

• •

Given noise-free observations, let $\bar{\mathbf{x}}\in\mathbb{R}^{n}$ be optimal for the primal problem, and $\bar{\boldsymbol{\lambda}}\in\mathbb{R}^{m}$ be optimal for the dual problem. In our problem, $n=12$ and $m=7$. According to the proof procedure of Lemma 3, we have $\bar{\mathbf{x}}=[\bar{\mathbf{e}};\bar{\mathbf{t}}]$ and $\bar{\boldsymbol{\lambda}}=\mathbf{0}$.

• •

Denote $\bar{\mathbf{Q}}\triangleq\mathbf{Q}_{\bar{\theta}}(\bar{\boldsymbol{\lambda}})\in\mathcal{S}^{n}$ as the Hessian of the Lagrangian at $\bar{\theta}$. In our problem, $\bar{\mathbf{Q}}=\begin{bmatrix}\bar{\mathbf{C}}&\mathbf{0}_{9\times 3}\\ \mathbf{0}_{3\times 9}&\mathbf{0}_{3\times 3}\end{bmatrix}\succeq 0$.

• •

Denote $X_{\theta}\triangleq\{\mathbf{x}\in\mathbb{R}^{n}|h_{i|\theta(\mathbf{x})}=0,i=1,\cdots,m\}$ as the primal feasible set given $\theta$, and denote $\bar{X}\triangleq X_{\bar{\theta}}$.

• •

Denote $\mathbf{h}(\mathbf{x})=[h_{1}(\mathbf{x}),\cdots,h_{m}(\mathbf{x})]$.

In QCQP (11), the objective $\mathbf{e}^{\top}\mathbf{C}\mathbf{e}$ is convex with $\mathbf{e}$. However, the presence of the auxiliary variables $\mathbf{t}$ makes the objective is not strictly convex. Theorem 5.1 in [36] provides a framework to prove the local stability for such kinds of problems.

Theorem 3** (Theorem 5.1 in [36]).**

Assume that the following $4$ conditions are satisfied:

RS (restricted Slater): There exists $\boldsymbol{\mu}\in\mathbb{R}^{m}$ such that $\boldsymbol{\mu}^{\top}\nabla\mathbf{h}_{\bar{\theta}}(\bar{\mathbf{x}})=0$ and $(\sum_{i=1}^{m}\mu_{i}\mathbf{A}_{i|\bar{\theta}})|_{V}\succ 0$, where $V\triangleq\{\mathbf{v}\in\mathbb{R}^{n}|\bar{\mathbf{Q}}\mathbf{v}=\mathbf{0},\bar{\mathbf{x}}^{\top}\mathbf{v}=0\}$.

R1 (constraint qualification): Abadie constraint qualification $\operatorname*{ACQ}_{\bar{\mathbf{X}}}(\bar{\mathbf{x}})$ holds.

R2 (smoothness): $\mathcal{W}\triangleq\{(\theta,\mathbf{x})|\mathbf{h}_{\theta}(\mathbf{x})=\mathbf{0}\}$ is a smooth manifold nearby $\bar{w}\triangleq(\bar{\theta},\bar{\mathbf{x}})$, and $\dim_{\bar{w}}\mathcal{W}=\dim\Theta+\dim_{\bar{\mathbf{x}}}\bar{X}$.

R3 (not a branch point): $\bar{\mathbf{x}}$ is not a branch point of $\bar{X}$ with respect to $\mathbf{v}\mapsto\bar{\mathbf{Q}}\mathbf{v}$.

Then the SDR relaxation is tight when $\theta$ is close enough to $\bar{\theta}$. Moreover, the QCQP has a unique optimal solution $\mathbf{x}_{\theta}$, and the SDR problem has a a unique optimal solution $\mathbf{x}_{\theta}\mathbf{x}_{\theta}^{\top}$.

Among the four conditions in Theorem 3, the RS (restricted Slater) is the main assumption, which is related to the convexity of the Lagrangian function. It corresponds to the strict feasibility of an SDP. R1$\sim$R3 are regularity assumptions which are related to the continuity of the Lagrange multipliers. In the following text, we prove that the restricted Slater and R3 are satisfied in problem (11), which builds an important foundation to prove the local stability. In our problem, $\mathbf{A}_{i}$ and $h_{i}$ are independent of $\theta$, thus $\mathbf{A}_{i|\bar{\theta}}=\mathbf{A}_{i}$ and $h_{i|\bar{\theta}}=h_{i}$.

Lemma 4** (Restricted Slater).**

For QCQP (11), suppose $\operatorname*{rank}(\bar{\mathbf{C}})=8$. Then there exists $\boldsymbol{\mu}\in\mathbb{R}^{m}$ such that $\boldsymbol{\mu}^{\top}\nabla\mathbf{h}_{\bar{\theta}}(\bar{\mathbf{x}})=0$ and $(\sum_{i=1}^{m}\mu_{i}\mathbf{A}_{i|\bar{\theta}})|_{V}\succ 0$, where $V\triangleq\{\mathbf{v}\in\mathbb{R}^{n}|\bar{\mathbf{Q}}\mathbf{v}=\mathbf{0},\bar{\mathbf{x}}^{\top}\mathbf{v}=0\}$.

Proof.

For constraint of Eqs. (12a)$\sim$(12g), its gradient is

[TABLE]

Let

[TABLE]

Note that for noise-free cases we have

[TABLE]

or equivalently

[TABLE]

By combining Eqs. (5.2)$\sim$(25), it can be verified that $\boldsymbol{\mu}^{\top}\nabla\mathbf{h}_{\bar{\theta}}(\bar{\mathbf{x}})=0$.

It remains to check the positivity condition. According to the definition of $V$, $\begin{bmatrix}\bar{\mathbf{Q}}\\ \bar{\mathbf{x}}^{\top}\end{bmatrix}\mathbf{v}=\mathbf{0}\Leftrightarrow\begin{bmatrix}\bar{\mathbf{C}}&\mathbf{0}\\ \bar{\mathbf{e}}^{\top}&\bar{\mathbf{t}}^{\top}\end{bmatrix}\mathbf{v}=\mathbf{0}$. Since $\bar{\mathbf{C}}$ is constructed by noise-free observations, $\bar{\mathbf{C}}\bar{\mathbf{e}}=0$. In other words, $\bar{\mathbf{e}}$ is orthogonal to the space spanned by $\bar{\mathbf{C}}$. It is given that $\operatorname*{rank}(\bar{\mathbf{C}})=8$, thus $\operatorname*{rank}\left(\begin{bmatrix}\bar{\mathbf{C}}\\ \bar{\mathbf{e}}^{\top}\end{bmatrix}\right)=9$. Considering $\mathbf{v}$ as a non-trivial solution of a homogeneous linear equation system, $\mathbf{v}$ can be expressed by a coordinate system $\mathbf{v}=\begin{bmatrix}\mathbf{0}_{9\times 1}\\ \mathbf{t}\end{bmatrix}$.

It can be seen that only $10\sim 12$-th entries in coordinate system $\mathbf{v}$, which correspond to $\mathbf{t}$, are nonzero. Take Hessian for variable $\mathbf{t}$ and calculate the linear combination with coefficient $\boldsymbol{\mu}$, then we have

[TABLE]

Recall that in the proof of Theorem 1 we have proved that the eigenvalues of $\bar{\mathbf{E}}\bar{\mathbf{E}}^{\top}$ are $1$, $1$, and [math]. So the eigenvalues of $\mathcal{A}(\boldsymbol{\mu})$ are $1$, $1$, and [math]. In addition, it can be verified that $\bar{\mathbf{t}}=[\bar{t}_{1},\bar{t}_{2},\bar{t}_{3}]^{\top}$ is the normalized eigenvector corresponding to eigenvalue [math] of $\mathcal{A}(\boldsymbol{\mu})$. Hence $V$ is the orthogonal complement of $\bar{\mathbf{t}}$.

For any vector $\mathbf{v}\in V\backslash\{\mathbf{0}\}$, its $10\sim 12$-th entries are orthogonal to $\bar{\mathbf{t}}$, so $\mathbf{v}_{[10:12]}^{\top}\mathcal{A}(\boldsymbol{\mu})\mathbf{v}_{[10:12]}$ is strictly positive. It follows that $(\sum_{i=1}^{m}\mu_{i}\mathbf{A}_{i|\bar{\theta}})|_{V}=\mathcal{A}(\boldsymbol{\mu})|_{(\bar{\mathbf{t}})^{\perp}}\succ 0$. ∎

Definition 1** (Branch Point [36]).**

Let $\pi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{k}$ be a linear map: $\mathbf{v}\mapsto\bar{\mathbf{Q}}\mathbf{v}$. Let $\bar{X}\subseteq\mathbb{R}^{n}$ be the zero set of the equation system $\mathbf{h}(\mathbf{x})=[h_{1}(\mathbf{x}),\cdots,h_{m}(\mathbf{x})]$, and let $T_{\mathbf{x}}\bar{X}\triangleq\ker(\nabla\mathbf{h}(\mathbf{x}))$ denote the tangent space of $\bar{X}$ at $\mathbf{x}$. We say that $\mathbf{x}$ is a branch point of $X$ with respect to $\pi$ if there is a nonzero vector $\mathbf{v}\in T_{\mathbf{x}}\ \bar{X}$ with $\pi(\mathbf{v})=\mathbf{0}$.

Lemma 5**.**

For QCQP (11), suppose $\operatorname*{rank}(\bar{\mathbf{C}})=8$. Then $\bar{\mathbf{x}}=[\bar{\mathbf{x}};\bar{\mathbf{t}}]$ is not a branch point of $\bar{X}$ with respect to the mapping $\pi:\mathbf{v}\mapsto\bar{\mathbf{Q}}\mathbf{v}$.

Proof.

It can be verified that $\pi(\mathbf{v})=\mathbf{0}\Leftrightarrow\bar{\mathbf{Q}}\mathbf{v}=\mathbf{0}\Leftrightarrow\begin{bmatrix}\bar{\mathbf{C}}&\mathbf{0}_{9\times 3}\\ \mathbf{0}_{3\times 9}&\mathbf{0}_{3\times 3}\end{bmatrix}\mathbf{v}=\mathbf{0}\Leftrightarrow\mathbf{v}=\begin{bmatrix}c\bar{\mathbf{e}}\\ \mathbf{t}\end{bmatrix}$, where $c$ and $\mathbf{t}$ are free parameters. The last equivalence takes advantage of $\operatorname*{rank}(\bar{\mathbf{C}})=8$. If $\bar{\mathbf{x}}$ is a branch point, there should exist a nonzero vector $\mathbf{v}\in\ker(\nabla\mathbf{h}(\bar{\mathbf{x}}))$, i.e., $\nabla\mathbf{h}(\bar{\mathbf{x}})\mathbf{v}=\mathbf{0}$. We will prove that such nonzero $\mathbf{v}$ does not exist. By substituting $\mathbf{v}=\begin{bmatrix}c\bar{\mathbf{e}}\\ \mathbf{t}\end{bmatrix}$ into equation $\nabla\mathbf{h}_{\bar{\theta}}(\bar{\mathbf{x}})\mathbf{v}=\mathbf{0}$ (see Eq. (5.2)), we obtain a homogeneous linear system with unknowns $\mathbf{t}=[t_{1},t_{2},t_{3}]^{\top}$ and $c$

[TABLE]

By eliminating $\mathbf{t}$ from Eqs. (27a)(27b)(27c)(27g), we obtain that $(\bar{\mathbf{e}}_{1}^{\top}\bar{\mathbf{e}}_{1}+\bar{\mathbf{e}}_{2}^{\top}\bar{\mathbf{e}}_{2}+\bar{\mathbf{e}}_{3}^{\top}\bar{\mathbf{e}}_{3})c=0$ and $c=0$. By substituting $c=0$ into this equation system, it can be further verified that this equation system only has zeros as its solution. So $\mathbf{v}$ can only be a zero vector. ∎

Theorem 4**.**

For QCQP (11) and its Lagrangian dual problem (19), let $\bar{\mathbf{C}}$ being constructed by noise-free observations and assume that $\operatorname*{rank}(\bar{\mathbf{C}})=8$. (i) There is zero-duality-gap whenever $\mathbf{C}$ is close enough to $\bar{\mathbf{C}}$. In other words, there exists a hypersphere of nonzero radius $\epsilon$ with center $\bar{\mathbf{C}}$, i.e., $\mathcal{B}(\bar{\mathbf{C}},\epsilon)=\|\mathbf{C}-\bar{\mathbf{C}}\|\leq\epsilon$, such that for any $\mathbf{C}\in\mathcal{B}(\bar{\mathbf{C}},\epsilon)$ there is zero-duality-gap. (ii) Moreover, the semidefinite relaxation problem (19) recovers the optimum of the original QCQP problem (11).

Proof.

From Lemma 3, when the point correspondences are noise-free, the relaxation is tight. From the proof procedure, the optimum is $\bar{\mathbf{x}}=[\bar{\mathbf{e}};\bar{\mathbf{t}}]$, which uniquely determines the zero-duality-gap parameter $\bar{\mathbf{C}}$. Our proof is an application of the Theorem 3. The four conditions needed in Theorem 3 are satisfied: (RS) From Lemma 4, the restricted Slater is satisfied. (R1) The equality constraints in the primal problem forms a variety. Abadie constraint qualification (ACQ) holds everywhere since the variety is smooth and the ideal is radical, see [36, Lemma 6.1]. (R2) smoothness. The normalized essential manifold is smooth [23]. (R3) From Lemma 5, $\bar{\mathbf{x}}$ is not a branch point of $\bar{X}$. ∎

Now we complete the proof that the SDR of the proposed QCQP is tight under low noise observations. In other words, when the observation noise is small, Theorem 4 guarantees that the optimum of the original QCQP can be found by optimizing its SDR. Finding the noise bounds that SDR can tolerate is still an open problem. In Section 7, we will demonstrate that the SDR is tight for large noise levels which are much larger than that in actual occurrence.

6 Robust $N$-Point Method

In previous sections, we propose a solution to non-minimal case essential matrix estimation. However, the presence of outliers in point correspondences is inevitable due to the ambiguities of feature points’ local appearance. When the correspondences are contaminated by outliers, the solution to minimize an algebraic error will be biased. To apply our method for outlier-contamination scenarios, a robust loss instead of least-square loss can be used in the objective function. We take advantage of M-estimators in robust statistics [14], and modify problem (7) as

[TABLE]

where $\rho(\cdot)$ is a robust loss function.

The selection of an appropriate robust loss function is critical to this problem. A natural choice for loss function would be the $\ell_{0}$ norm, i.e., $\rho(x;\tau)=\begin{cases}0,&\text{if }\ |x|\leq\tau,\\ 1,&\text{otherwise}.\end{cases}$. It turns out to be the inlier set maximization problem, which has been proved to be a computationally expensive problem [37]. The presented approach below can accommodate many M-estimators in a computationally efficient framework. First we will take the scaled Welsch (Leclerc) function [58, 59] as an example

[TABLE]

where $\tau\in(0,+\infty)$ is a scale parameter.

Problem (28) can be optimized by local optimization methods directly. However, local optimization is prone to local optima. To solve this problem, our approach is based on Black-Rangarajan duality between M-estimators and line processes [17]. The duality introduces an auxiliary variable $w_{i}$ for $i$-th point correspondence and optimizes a joint objective over the essential matrix $\mathbf{E}$ and the line process variables $\mathbb{W}=\{w_{i}\}_{i=1}^{N}$

[TABLE]

Here $\Psi(w_{i})$ is a loss on ignoring $i$-th correspondence. $\Psi(w_{i})$ tends to zero when the $i$-th correspondence is active and to one when it is inactive. A broad variety of robust estimators $\rho(\cdot)$ have corresponding loss functions $\Psi(\cdot)$ such that problems (28) and (30) are equivalent with respect to $\mathbf{E}$: optimizing either of the two problem yields the same essential matrix. The form of problem (30) enables efficient and scalable optimization by an iterative solution of the weighted $N$-point method. This yields a general approach that can accommodate many robust nonconvex functions $\rho(\cdot)$.

The loss function that makes problems (28) and (30) equivalent with respect to $\mathbf{E}$ is

[TABLE]

Objective (30) is biconvex on $\mathbf{E},\mathbb{W}$. When $\mathbf{E}$ is fixed, the optimal value of each $w_{i}$ has a closed-form solution. When variables $\mathbb{W}$ are fixed, objective (30) turns into a weighted $N$-point problem, which can be efficiently solved by the proposed method. Specifically, we exploit this special structure and optimize the objective by alternatively updating variable sets $\mathbf{E}$ and $\mathbb{W}$. As a block coordinate descent algorithm, this alternating minimization scheme provably converges. By fixing $\mathbf{E}$, the optimal value of each $w_{i}$ is given by

[TABLE]

This can be verified by substituting Eq. (32) into objective (30), which yields a constant $1$ with respect to $\mathbf{E}$. Thus optimizing objective (30) yields a solution $\mathbf{E}$ that is also optimal for the original objective (28). The line process theory [17] provides general formulations to calculate $w_{i}$ for a broad variety of robust loss functions $\rho(x;\tau)$. The alternating minimization of a line process is illustrated in Fig. 4.

It is worth mentioning that the loss function in our robust $N$-point method is not limited to Welsch function. The Black-Rangarajan duality [17] provides a theory to build the relation between M-estimators and line processes. Based on this general-purpose framework, integrating other robust loss functions is essentially the same as the Welsch function. Typical robust loss functions include Cauchy (Lorentzian), Charbonnier (pseudo-Huber, $\ell_{1}$-$\ell_{2}$), Huber, Geman-McClure, smooth truncated quadratic, truncated quadratic, Tukey’s biweight functions, etc, see Fig. 5(a). The Black-Rangarajan duality of these loss functions can be found in [17, 59].

Objective (30) is non-convex and its shape is controlled by the parameter $\tau$ of loss function $\rho(\cdot)$. To alleviate the effect of local minima, we employ graduated non-convexity (GNC) which is dating from 1980s [60, 46]. The idea is to start from a problem that is easier to solve, then progressively deformed to the actual objective while tracking the solution along the way. Large $\tau$ makes the objective function smoother and allows many correspondences to participate in the optimization even when they are not fit well by the essential matrix $\mathbf{E}$, see Fig. 5(b). Our method begins with a very large value $\tau$. Over the iterations, $\tau$ is automatically decreased, gradually introducing non-convexity into the objective. The line process together with the GNC strategy is not guaranteed to obtain the global optimum for the original non-convex problem (28). They can obtain a good solution if not the globally optimal one in most cases when the outlier ratio is below a threshold. The outline of the robust $N$-point method is shown in Algorithm 2.

7 Experimental Results

Setting for the $N$-point method: We compared the proposed $N$-point method with several state-of-the-art methods on synthetic and real data. Specifically, we compared our method with $5$ classical or state-of-the-art methods:

• •

a general five-point method 5pt [2] for relative pose estimation.

• •

two general methods for fundamental matrix estimation, including seven-point method 7pt [1] and eight-point method 8pt [11].

• •

eigenvalue-based method proposed by Kneip and Lynen [9] which is referred to as eigen and a certifiably globally optimal solution by Briales et al. [10] which is referred to as SDP-Briales.

Among these methods, the implementation of SDP-Briales is provided by the authors. The implementations of other comparison methods are provided by OpenGV [61]. The method eigen needs initialization, and it is initialized using 8pt method. Our method and all comparison methods in openGV are implemented in C++. The SDP-Briales method is implemented in Matlab, and the optimization in it relies on an SDP solver with hybrid Matlab/C++ programming.

The eigen method in openGV implementation provides rotation only. Once the relative rotation has been obtained, we can calculate the translation $\mathbf{t}$. Recall that $\mathbf{f}_{i}$ and $\mathbf{f}^{\prime}_{i}$ represent bearing vectors of a point correspondence across two images. From Eqs. (1) and (6), the epipolar constraint can be written as

[TABLE]

Since $\mathbf{R}$ has been calculated, each point correspondence provides a linear constraint on the entries of the translation vector $\mathbf{t}$. Due to the scale-ambiguity, the translation has only two DoFs. After DLT, a normalized version of $\mathbf{t}$ can be recovered by simple linear derivation of the right-hand null-space vector (e.g. via singular value decomposition). Given $N$ ($N\geq 2$) point correspondences, the least squares fit of $\mathbf{t}$ can be determined by considering the singular vector corresponding to the smallest singular value.

Setting for robust $N$-point method: We compared the proposed robust $N$-point method with RANSAC+5pt, RANSAC+7pt and RANSAC+8pt, which stands for integrating 5pt, 7pt and 8pt into the RANSAC framework [5], respectively. All comparison methods are provided by openGV [61], and the default parameters are used. In the experiment on real-world data, we also compare our method with a branch-and-bound method BnB-Yang222The C++ code is available from http://jlyang.org/ [41]. The angular error threshold in this method was set as $0.002$ radians.

To evaluate the performance of the proposed method, we separately compared the relative rotation and translation accuracy. We follow the criteria defined in [62] for quantitative evaluation. Specifically,

• •

the angle difference for rotations is defined as

[TABLE]

• •

and the translation direction error is defined as

[TABLE]

In above criteria, $\mathbf{R}_{\text{true}}$ and $\mathbf{R}^{\star}$ are the ground truth and estimated rotation, respectively; $\mathbf{t}_{\text{true}}$ and $\mathbf{t}^{\star}$ are the ground truth and estimated translation, respectively.

7.1 Efficiency of $N$-Point Method

The SDPA solver [63] was adopted as an SDP solver in our methods. All experiments were performed on an Intel Core i7 CPU running at $2.40$ GHz. The number of point correspondences is fixed to $100$. The SDP optimization takes about $5$ ms. In addition, it takes about $1$ ms for other procedures in our method, including problem construction, optimal essential matrix recovery, and pose decomposition. In summary, the runtime of our method is about $6$ ms.

We compared the proposed method with several state-of-the-art methods [18, 8, 10], which also aim to find the globally optimal relative pose. The efficiency comparison is shown in Table I. It can be seen that our method is $2\sim 3$ orders of magnitude faster than comparison methods. The superior efficiency makes our method the first globally optimal method that can be applied to large scale structure-from-motion and realtime SLAM applications.

Since both the methods of SDP-Briales and ours take advantage of SDP optimization, further comparison between them is reported in Table II. It can be seen that our method has a much simpler formulation in terms of numbers of variables and constraints. It is not surprising that our method has significantly better efficiency.

7.2 Accuracy of $N$-Point Method

7.2.1 Synthetic Data

To thoroughly evaluate our method’s performance, we perform experiments on synthetic scenes in a similar manner to [9]. We generate random scenes by first fixing the position of the first frame to the origin and its orientation to the identity. The translational offset of the second frame is chosen with uniformly distributed random direction and a maximum magnitude of $2$. The orientation of the second frame is generated with random Euler angles bounded to $0.5$ radians in absolute value. This generates random relative poses as they would appear in practical situations. Point correspondences result from uniformly distributed random points around the origin with a distance varying between $4$ and $8$, transforming those points into both frames. Then a virtual camera is defined with a focal length of $800$ pixels. Gaussian noise is added by perturbing each point correspondence. The standard deviation of the Gaussian noise is referred to as the noise level.

First, we test image noise resilience. For each image noise level, we randomly generate synthetic scenes and repeat the experiments $1000$ times. The number of correspondences is fixed to $10$, and the step size of the noise level is $0.1$ pixels. The results for all methods with varying image noise levels are shown in Fig. 6. It can be seen that our method consistently has smaller rotation error $\varepsilon_{\text{rot}}$ and translation error $\varepsilon_{\text{tran}}$ than other methods. Moreover, our method and eigen significantly outperform 5pt, 7pt and 8pt. This result demonstrates the advantage of non-minimal solvers in terms of accuracy.

Second, we set the noise level as $5$ pixels and vary the number $N$ of point correspondences. The step size of the correspondence number is $5$. The methods 5pt, 7pt, and 8pt can only take a small subset of the point correspondences. In contrast, eigen and ours utilize all the point correspondences. To make a fair comparison, we randomly sample minimal number of point correspondences for 5pt, 7pt, and 8pt, and repeat $20$ times for each method. Then we find the optimal relative pose among them. Since all point correspondences are inliers, we cannot use the maximal inlier criterion to find the optimal rotation as that in the traditional RANSAC framework. Instead, we use the algebraic error to find the optimal relative rotation for these methods.

The pose estimation accuracy with respect to the number of point correspondences is shown in Fig. 7. We have the following observations: (1) The errors of eigen and ours decrease when increasing the numbers of point correspondences. It further verifies the effectiveness of non-minimal solvers. (2) When $N>20$, eigen and ours have significantly smaller errors than other methods, and our method has the smallest rotation and translation error among all methods. (3) The error curve of Eigen has oscillation due to local minima when the noise level is large. In contrast, the error curve of our method is smooth for any noise level.

7.2.2 Real-World Data

We further provide an experiment on real-world images from the EPFL Castle-P19 dataset [64]. It contains $19$ images in this dataset. We generate $18$ wide-baseline image pairs by grouping adjacent images. For each image pair, putative point correspondences are determined by SIFT feature [65]. Then we use RANSAC with iteration number $2000$ and Sampson distance threshold $1.0\times 10^{-3}$ for outlier removal. Since the iteration number is sufficiently large and the distance threshold is small, the preserved correspondences can be treated as inliers.

Given correct point correspondences, we compare the relative pose accuracy of different methods. For 5pt, 7pt and 8pt methods, they only use a small portion of the correspondences. For a fair comparison, we repeat them $10$ times using randomly sampled subsets. The rotation and translation errors of different methods are shown in Fig. 9. It can be seen that the mean and median error of our $N$-point method is significantly smaller than those errors produced by the comparison methods. Specifically, our method achieves a median rotation error of $0.15^{\circ}$ and a median translation error of $0.56^{\circ}$. In contrast, 5pt, 7pt and 8pt achieve a median rotation error of $0.33^{\circ}$, $0.99^{\circ}$ and $1.02^{\circ}$, respectively; and they achieve a median translation error of $1.07^{\circ}$, $3.07^{\circ}$ and $3.01^{\circ}$, respectively. In this experiment, eigen and our method have a negligible difference. It means that when the noise level is small, local optimization methods might work as well as global optimization methods.

7.2.3 Performance for Pure Rotation

We use synthetic data to validate the performance of our method for pure rotation. The synthetic scenes are generated as that in Section 7.2.1. The field-of-view of the camera is $180^{\circ}$. We set the noise level as $0.5$ pixels and the number of point correspondences as $100$. The translation length is varied. For each translation length, we repeat the experiments $1000$ times using randomly generated scenes. The pose estimation errors with respect to translation length is shown in Fig. 10. It can be seen that the rotation error is not affected by translation length. In contrast, the translation error is greatly affected by translation length333When the ground truth of translation is zero, the translation error $\varepsilon_{\text{rot}}$ is ill-defined. We define $\varepsilon_{\text{rot}}$ as $90^{\circ}$, because the expectation of $\varepsilon_{\text{rot}}$ is $90^{\circ}$ when the estimated translation is uniformly random. . The larger translation length tends to result in more accurate translation estimation.

In Figure 10(c), we plot a pure rotation statistic [54] with respect to translation length. The statistic is defined as the mean of $\frac{\mathbf{x}_{i}\times\mathbf{R}^{\star}\mathbf{x}^{\prime}_{i}}{\|\mathbf{x}_{i}\|\|\mathbf{x}^{\prime}_{i}\|}$, $i=1,\cdots,N$. It can be seen that this statistic is discriminative to identify pure rotation scenarios. When (near) pure rotation occurs, this statistic is (near) zero. From the above experiments, it certifies that the proposed $N$-point method can be applied to pure rotation cases.

7.3 Global Optimality of $N$-Point Method

Recall that $\operatorname*{rank}(\mathbf{X}_{e}^{\star})$ is used to verify the optimality of the proposed $N$-point method. Thus the second largest singular value is the key to ensure rank-$1$ condition and global optimality. We use the synthetic scenes defined in Section 7.2.1 to demonstrate the second largest singular value of $\operatorname*{rank}(\mathbf{X}_{e}^{\star})$. Figure 11 shows the cumulative distribution functions (CDF) of the second largest singular values under different settings. Given a threshold to determine the rank of $\operatorname*{rank}(\mathbf{X}_{e}^{\star})$, the proportion of instances with global optimality can be obtained from corresponding CDF. The success rate of the global optimality depends on the threshold of singular value and the accuracy of the SDP solver.

In Figure 11(a)(b), the number of point correspondences is fixed to $100$, and the noise level is varied for forward and sideways motion modes. In Figure 11(c)(c), the noise level is fixed to $0.5$ pixels, and the number of point correspondences is varied for forward and sideways motion modes. We have the following observations. (1) When the threshold of singular value is $1.0\times 10^{-6}$, the global optimality can be obtained in most cases. (2) Smaller noise levels or more point correspondences will result in higher success rates of global optimality. (3) For usual cases whose noise level is below $1$ pixel and the number of point correspondences is above $20$, the global optimality can be obtained in most cases.

7.4 Performance of Robust $N$-Point Method

We test the robust $N$-point method on both synthetic data and real-world data. The parameter $\tau^{2}_{\text{min}}$ is set as $6.0\times 10^{-7}$.

7.4.1 Synthetic Data

The synthetic scene is generated as that in Section 7.2.1. The noise level is fixed to $0.5$ pixels, and the correspondence number is fixed to $100$. The outlier ratio is varied from $0\%$ to $100\%$ with a step size of $5\%$. For each outlier ratio, we repeat the experiments $100$ times using randomly generated data.

We evaluate $8$ robust loss functions in Fig. 5(a), and report their success rates. The success rate is the ratio of successful cases to the overall trials. A trial is treated as a success if $\varepsilon_{\text{rot}}\leq 0.15^{\circ}$ and $\varepsilon_{\text{tran}}\leq 0.5^{\circ}$. The results are shown in Fig. 12. It can be seen that truncated quadratic, Welsch, smooth truncated quadratic, and Turkey’s biweight functions can tolerate up to $45\%$ outliers. Cauchy and Geman-McClure functions can tolerate up to $40\%$ outliers; Charbonnier and Huber functions can tolerate up to $15\%$ outliers. In the following experiments, we will use Welsch function in the robust $N$-point method since it has superior performance.

The performance of the robust $N$-point method and RANSAC-based methods is shown in Fig. 13. Our method can consistently obtain smaller errors than RANSAC-based methods in terms of both mean and median errors. From logarithmic scale plots in Fig. 13(e)$\sim$(h), our method has significantly smaller rotation and translation error than other methods when the outlier ratio is below a threshold.

7.4.2 Real-World Data

We use the EPFL Castle-19 dataset [64] for real-world data. The image pairs and their putative correspondences are generated in the same way as that in Section 7.2.2. Each image pair contains $4408$ putative correspondences on average. Though it is theoretically sound, the BnB-Yang method [41] cannot handle a large number of point correspondences. It fails to return a solution for an image pair in a day, therefore we randomly select $100$ putative correspondences as input for it.

The relative pose accuracy of different methods is shown in Fig. 14. It can be seen that the proposed robust $N$-point method has higher overall accuracy than RANSAC-based methods and BnB-Yang method. The efficiency comparison between our method and RANSAC-based methods is summarized as below: (1) The robust $N$-point method takes a roughly constant number of iterations. According to the GNC strategy and its parameter setting, its iteration number is at most $79$. In contrast, the iteration in RANSAC-based methods increases drastically with a high outlier ratio or a high confidence level. For example, given $45\%$ outliers and $99\%$ confidence level, 5pt, 7pt and 8pt methods need to iterate at least $90$, $301$ and $548$ times, respectively. In practice, the RANSAC-based methods need more iterations than the least required to achieve high accuracy. (2) The non-minimal solver in the robust $N$-point method takes much more time than minimal solvers in RANSAC-based methods ($6$ milliseconds vs. $30\sim 180$ microseconds). As a result, our method takes about $480$ ms, and RANSAC-based methods take $52$ ms on average. (3) The BnB-Yang method takes $118\sim 9007$ seconds for randomly selected $100$ point correspondences.

8 Conclusions

This paper introduces a novel non-minimal solver for $N$-point problem in essential matrix estimation. First, we reformulate $N$-point problem as a simple QCQP by proposing an equivalent form of the essential manifold. Second, semidefinite relaxation is exploited to convert this problem to an SDP problem, and pose recovery from an optimal solution of SDP is proposed. Finally, a theoretical analysis of tightness and local stability is provided. Our method is stable, globally optimal, and relatively easy to implement. In addition, we propose a robust $N$-point method by integrating the non-minimal solver into M-estimators. Extensive experiments demonstrate that the proposed $N$-point method can find and certify the global optimum of the optimization problem, and it is $2\sim 3$ orders of magnitude faster than state-of-the-art non-minimal solvers. Moreover, the robust $N$-point method outperforms state-of-the-art methods in terms of robustness and accuracy.

Appendix A Another Formulation of $N$-Point Problem

The proposition below provides another equivalent form to define essential manifold $\mathcal{M}_{\mathbf{E}}$, which will derivate another simple optimization problem.

Lemma 6**.**

For an essential matrix $\mathbf{E}$ which can be decomposed by $\mathbf{E}=[\mathbf{t}]_{\times}\mathbf{R}$, it satisfies that $\emph{trace}(\mathbf{E}\mathbf{E}^{\top})=\sum_{i=1}^{3}\sum_{j=1}^{3}\mathbf{E}_{ij}^{2}=2\|\mathbf{t}\|^{2}$.

Proof.

Note that the norm of each row of $\mathbf{R}$ is $1$, and the rows of $\mathbf{R}$ are orthogonal to each other. Taking advantage of $\mathbf{E}=[\mathbf{t}]_{\times}\mathbf{R}$, it can be verified the $\operatorname*{trace}(\mathbf{E}\mathbf{E}^{\top})=\sum_{i=1}^{3}\sum_{j=1}^{3}\mathbf{E}_{ij}^{2}=2(t_{1}^{2}+t_{2}^{2}+t_{3}^{2})=2\|\mathbf{t}\|^{2}$. ∎

Proposition 1** (Proposition 7.3 in [49]).**

A real $3\times 3$ matrix, $\mathbf{E}$, is an essential matrix if and only if it satisfies the equation:

[TABLE]

Proposition 2**.**

A real $3\times 3$ matrix, $\mathbf{E}$, is an essential matrix in $\mathcal{M}_{\mathbf{E}}$ if and only if it satisfies the following two conditions:

[TABLE]

Proof.

For the if direction, by combining conditions (i) and (ii) we obtain Eq. (34). According to Theorem 1, $\mathbf{E}$ is a valid essential matrix. So there exist (at least) a pair of $\mathbf{t}\in\mathbb{R}^{3}$ and $\mathbf{R}\in\text{SO}(3)$ such that $\mathbf{E}=[\mathbf{t}]_{\times}\mathbf{R}$. According to condition (i) and Lemma 6, we have $\operatorname*{trace}(\mathbf{E}\mathbf{E}^{\top})=2\|\mathbf{t}\|^{2}=2$, which means $\|\mathbf{t}\|=1$. Thus we prove $\mathbf{E}\in\mathcal{M}_{\mathbf{E}}$.

For only if direction, since $\mathbf{E}\in\mathcal{M}_{\mathbf{E}}$, it is straightforward that condition (i) is satisfied according to Lemma 6. In addition, Eq. (34) is satisfied since $\mathbf{E}$ is an essential matrix. By substituting condition (i) in Eq. (34), we obtain condition (ii). ∎

According to Proposition 2, an equivalent form of minimizing the algebraic error is

[TABLE]

By introducing an auxiliary matrix $\mathbf{G}$, this problem can be reformulated as a QCQP

[TABLE]

Note that $\mathbf{G}$ is a symmetric matrix which introduces $6$ variables. Thus there are $15$ variables and $16$ constraints in this QCQP.

Appendix B Strong Duality Between Primal SDP and Its Duality

Lemma 7**.**

For QCQP (11), there is no duality gap between the primal SDP problem (18) and its dual problem (19).

Proof.

Denote the optimal value for problem (18) and its dual problem (19) as $f_{\text{primal}}$ and $f_{\text{dual}}$. The inequality $f_{\text{primal}}\geq f_{\text{dual}}$ follows from weak duality. Equality, and the existence of $\mathbf{X}^{\star}$ and $\boldsymbol{\lambda}^{\star}$ which attain the optimal values follow if we can show that the feasible regions of both the primal and dual problems have nonempty interiors, see [52, Theorem 3.1] (also known as Slater’s constraint qualification [66].)

For the primal problem, let $\mathbf{E}_{0}$ be an arbitrary point on the essential manifold $\mathcal{M}_{\mathbf{E}}$: $\mathbf{E}_{0}=[\mathbf{t}_{0}]_{\times}\mathbf{R}_{0}$, where $\|\mathbf{t}_{0}\|=1$. Denote $\mathbf{x}_{0}=[\text{vec}(\mathbf{E}_{0});\mathbf{t_{0}}]$. It can be verified that $\mathbf{X}_{0}=\mathbf{x}_{0}\mathbf{x}_{0}^{\top}$ is an interior in the feasible domain of the primal problem. For the dual problem, let $\boldsymbol{\lambda}_{0}=[-1,-1,-1,0,0,0,-3]^{\top}$. Recall that $\mathbf{C}\succeq 0$ and it can be verified that $\mathbf{Q}(\boldsymbol{\lambda}_{0})=\begin{bmatrix}\mathbf{C}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\end{bmatrix}+\mathbf{I}\succ 0$. That means $\boldsymbol{\lambda}_{0}$ is an interior in the feasible domain of the dual problem. ∎

Acknowledgments

The author would like to thank Prof. Laurent Kneip at ShanghaiTech, Prof. Qian Zhao at XJTU, and Haoang Li at CUHK for fruitful discussions. The author also thanks Dr. Jesus Briales at the University of Malaga for providing the code of [10] and Dr. Yijia He at CASIA for his help in the experiments.