Improved Pose Graph Optimization for Planar Motions Using Riemannian Geometry on the Manifold of Dual Quaternions

TL;DR

This paper introduces a Riemannian geometry-based method for planar pose graph optimization that inherently considers the nonlinear SE(2) group structure, leading to improved convergence robustness over existing methods.

Contribution

It proposes a novel Riemannian pose graph optimizer (RPG-Opt) using dual quaternions and exponential retraction, enhancing robustness in large uncertainty scenarios.

Findings

Achieves equivalent accuracy to state-of-the-art methods.

Demonstrates better convergence robustness under large odometry uncertainties.

Validated on public planar pose graph datasets.

Abstract

We present a novel Riemannian approach for planar pose graph optimization problems. By formulating the cost function based on the Riemannian metric on the manifold of dual quaternions representing planar motions, the nonlinear structure of the SE(2) group is inherently considered. To solve the on-manifold least squares problem, a Riemannian Gauss-Newton method using the exponential retraction is applied. The proposed Riemannian pose graph optimizer (RPG-Opt) is further evaluated based on public planar pose graph data sets. Compared with state-of-the-art frameworks, the proposed method gives equivalent accuracy and better convergence robustness under large uncertainties of odometry measurements.

Figures

Tables

Table 1. Table 1 : Evaluation results under ordinary noise.

Dataset		RPG-Opt10	g2o10	GTSAM	iSAM
CSAIL	I	$1.07\cdot {10}^{-1}$	$1.07\cdot {10}^{-1}$	$1.07\cdot {10}^{-1}$	$1.07\cdot {10}^{-1}$
	$𝛀$	$3.90\cdot {10}^1$	$3.90\cdot {10}^1$	$3.90\cdot {10}^1$	$9.40\cdot {10}^2$
FR079	I	$7.19\cdot {10}^{-2}$	$7.19\cdot {10}^{-2}$	$7.19\cdot {10}^{-2}$	$7.19\cdot {10}^{-2}$
	$𝛀$	$3.76\cdot {10}^1$	$3.76\cdot {10}^1$	$3.76\cdot {10}^1$	$3.50\cdot {10}^2$
FRH	I	$1.39\cdot {10}^{-6}$	$3.19\cdot {10}^{-4}$	$3.19\cdot {10}^{-4}$	$3.46\cdot {10}^{-4}$
	$𝛀$	$1.93\cdot {10}^{-4}$	$4.18\cdot {10}^{-2}$	$4.18\cdot {10}^{-2}$	$4.61\cdot {10}^{-2}$
M3500	I	$3.02\cdot {10}^0$	$3.02\cdot {10}^0$	$3.02\cdot {10}^0$	$3.02\cdot {10}^0$
	$𝛀$	$1.38\cdot {10}^2$	$1.38\cdot {10}^2$	$1.38\cdot {10}^2$	$1.39\cdot {10}^2$
MITb	I	$6.60\cdot {10}^0$	$2.70\cdot {10}^1$	$2.83\cdot {10}^6$	$2.78\cdot {10}^0$
	$𝛀$	$3.02\cdot {10}^3$	$7.71\cdot {10}^2$	$4.49\cdot {10}^9$	$2.26\times {10}^2$
City10K	I	$8.72\cdot {10}^0$	$8.72\cdot {10}^0$	$6.77\cdot {10}^6$	$8.72\cdot {10}^0$
	$𝛀$	$5.12\cdot {10}^2$	$5.12\cdot {10}^2$	$2.48\cdot {10}^8$	$5.18\cdot {10}^2$
M10K	I	$3.03\cdot {10}^2$	$3.03\cdot {10}^2$	$3.17\cdot {10}^7$	$3.03\cdot {10}^2$
	$𝛀$	$1.99\cdot {10}^5$	$1.98\cdot {10}^5$	$2.91\cdot {10}^{12}$	$2.93\cdot {10}^8$

Table 2. Table 2 : Summarized evaluation results for synthetic data sets with chordal relaxation-based initialization.

	RPG-Opt		chord + RPG-Opt		chord + g2o		chord + GTSAM
Dataset	$e_{\text{t}}$	$e_{\text{r}}$	$e_{\text{t}}$	$e_{\text{r}}$	$e_{\text{t}}$	$e_{\text{r}}$	$e_{\text{t}}$	$e_{\text{r}}$
M3500a+	$2.60\cdot {10}^{-1}$	$8.42\cdot {10}^0$	$2.60\cdot {10}^{-1}$	$8.42\cdot {10}^0$	$2.60\cdot {10}^{-1}$	$8.42\cdot {10}^0$	$3.310\cdot {10}^0$	$1.55\cdot {10}^1$
M3500b+	$3.79\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$3.79\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$3.80\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$4.37\cdot {10}^0$	$2.03\cdot {10}^1$
M3500c+	$3.94\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$3.94\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$3.91\cdot {10}^{-1}$	$1.31\cdot {10}^1$	$1.14\cdot {10}^0$	$1.36\cdot {10}^1$
M3500d+	$8.28\cdot {10}^{-1}$	$1.51\cdot {10}^1$	$4.71\cdot {10}^{-1}$	$1.51\cdot {10}^1$	$4.78\cdot {10}^{-1}$	$1.51\cdot {10}^1$	$7.01\cdot {10}^{-1}$	$1.57\cdot {10}^1$
City10000a	$4.34\cdot {10}^{-2}$	$1.55\cdot {10}^0$	$4.34\cdot {10}^{-2}$	$1.55\cdot {10}^0$	$4.34\cdot {10}^{-2}$	$1.55\cdot {10}^0$	$4.33\cdot {10}^{-2}$	$1.55\cdot {10}^0$
City10000b	$3.20\cdot {10}^{-1}$	$6.95\cdot {10}^0$	$5.32\cdot {10}^{-2}$	$2.66\cdot {10}^0$	$5.32\cdot {10}^{-2}$	$2.66\cdot {10}^0$	$5.29\cdot {10}^{-2}$	$2.66\cdot {10}^0$
City10000c	$2.86\cdot {10}^{-2}$	$5.80\cdot {10}^0$	$5.77\cdot {10}^{-2}$	$2.67\cdot {10}^0$	$5.77\cdot {10}^{-2}$	$2.67\cdot {10}^0$	$5.77\cdot {10}^{-2}$	$2.67\cdot {10}^0$
City10000d	$3.88\cdot {10}^{-1}$	$1.39\cdot {10}^1$	$2.46\cdot {10}^{-1}$	$1.22\cdot {10}^1$	$2.46\cdot {10}^{-1}$	$1.23\cdot {10}^{-1}$	$9.28\cdot {10}^{-1}$	$1.39\cdot {10}^1$

Equations

$$\mathbf{r}={\left[\begin{array}{c}\cos (\theta /2),\sin (\theta /2)\end{array}\right]}^{\top }\in {\mathbb{S}}^1\subset {\mathbb{R}}^2,$$

$${\mathcal{Q}}_{\text{r}}^{└}=\left[\begin{array}{cc}{r}_0& -r_1\\ r_1& r_0\end{array}\right],{\mathcal{Q}}_{\text{s}}^{┘}=\left[\begin{array}{cc}{s}_0& s_1\\ -s_1& s_0\end{array}\right].$$

$$\mathbf{x}=[{\mathbf{x}}_{\text{r}}^{\top },{\mathbf{x}}_{\text{s}}^{\top }{]}^{\top }\in {\mathbb{S}}^1\times {\mathbb{R}}^2:=\mathbb{M},$$

$${\mathbf{x}}_{\text{s}}=\frac{1}{2}\mathbf{t}\otimes {\mathbf{x}}_{\text{r}}=\frac{1}{2}{\mathcal{Q}}_{\text{r}}^{┘}\mathbf{t}.$$

$$\mathcal{Q}_{\text{x}}^{\ulcorner}=\Bigg{[}\begin{smallmatrix}x_{0}&-x_{1}&~{}~{}0&~{}~{}0\\ x_{1}&~{}~{}x_{0}&~{}~{}0&~{}~{}0\\ x_{2}&~{}~{}x_{3}&~{}~{}x_{0}&-x_{1}\\ x_{3}&-x_{2}&~{}~{}x_{1}&~{}~{}x_{0}\end{smallmatrix}\Bigg{]}\,,\quad\mathcal{Q}_{\text{y}}^{\urcorner}=\Bigg{[}\begin{smallmatrix}y_{0}&-y_{1}&~{}~{}0&~{}~{}0\\ y_{1}&~{}~{}y_{0}&~{}~{}0&~{}~{}0\\ y_{2}&-y_{3}&~{}~{}y_{0}&~{}~{}y_{1}\\ y_{3}&~{}~{}y_{2}&-y_{1}&~{}~{}y_{0}\end{smallmatrix}\Bigg{]}\,.$$

$${\mathbf{v}}^{\prime }=\mathbf{x}⊞\mathbf{v}⊞{\mathbf{x}}^{\ast },$$

$$g(\mathbf{x})=\left[\begin{array}{c}{\mathbf{x}}_{\text{r}}^{\top }{\mathbf{x}}_{\text{r}}-1\\ 0\end{array}\right],\mathbf{x}=[{\mathbf{x}}_{\text{r}},{\mathbf{x}}_{\text{s}}{]}^{\top }\in {\mathbb{R}}^4,$$

$$\text{mathds}{T}_{\mathbf{x}}\mathbb{M}=\text{ker}(\nabla g)=\text{ker}\left(\left[\begin{array}{cc}2{\mathbf{x}}_{\text{r}}^{\top }& 0^{\top }\\ 0^{\top }& 0^{\top }\end{array}\right]\right).$$

$${\mathcal{P}}_{\mathbf{x}}=\left[\begin{array}{cc}\mathbf{I}-{\mathbf{x}}_{\text{r}}{\mathbf{x}}_{\text{r}}^{\top }& 0\\ 0& \mathbf{I}\end{array}\right]\in {\mathbb{R}}^{4\times 4},$$

$$\theta \mathbf{m}$$

$$\displaystyle=\frac{1}{\text{sinc}(\theta/2)}\big{(}{\bf t}\times{\bf k}\sin\frac{\theta}{2}+{\bf t}\cos\frac{\theta}{2}\big{)}$$

$$=\frac{1}{\text{sinc}(\theta /2)}\left[\begin{array}{cc}\cos (\theta /2)& \sin (\theta /2)\\ -\sin (\theta /2)& \cos (\theta /2)\end{array}\right]\mathbf{t}.$$

$$0.5\,\theta\,{\bf m}=\frac{\big{[}\begin{smallmatrix}~{}~{}x_{0}&x_{1}\\ -x_{1}&x_{0}\\ \end{smallmatrix}\big{]}{\bf t}}{2\,\text{sinc}(\theta/2)}=\frac{\mathcal{Q}_{\text{r}}^{\lrcorner}{\bf t}}{2\,\text{sinc}(\theta/2)}=\frac{[\,x_{2},\,x_{3}\,]^{\top}}{\text{sinc}(\theta/2)}\,.$$

$${\mathbf{x}}_{\text{t}}={\text{Log}}_{\text{mathds}1}(\mathbf{x})=\frac{1}{\gamma }[x_1,x_2,x_3{]}^{\top }\in {\mathbb{R}}^3,$$

$${\text{Log}}_{\mathbf{x}}(\mathbf{y})=\mathbf{x}⊞[0,{\text{Log}}_{\text{mathds}1}({\mathbf{x}}^{-1}⊞\mathbf{y}{)}^{\top }{]}^{\top }\in \text{mathds}{T}_{\mathbf{x}}\mathbb{M}.$$

$$\mathbf{x}={\text{Exp}}_{\text{mathds}1}({\mathbf{x}}_{\text{t}})=[\cos (x_{\text{t},1}),\gamma {\mathbf{x}}_{\text{t}}^{\top }{]}^{\top }.$$

$$\text{Exp}_{\bf x}({\bf y}_{\text{t}})={\bf x}\boxplus\text{Exp}_{\mathds{1}}\big{(}({\bf x}^{-1}\boxplus{\bf y}_{\text{t}})_{1:3}\big{)}\,.$$

$${\tilde{\mathbf{x}}}^{\ast }=\underset{\tilde{\mathbf{x}}\in {\mathbb{M}}^n}{arg\; min}{\sum }_{⟨i,j⟩\in \text{mathds}C}{\mathbf{e}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathbf{e}}_{ij},$$

$$\mathcal{F}(\tilde{\mathbf{x}})={\sum }_{⟨i,j⟩\in \text{mathds}C}{\mathbf{e}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathbf{e}}_{ij}$$

$${\mathbf{e}}_{ij}={\text{Log}}_{\text{mathds}1}({\mathbf{z}}_{ij}^{-1}⊞{\mathbf{x}}_i^{-1}⊞{\mathbf{x}}_j)\in {\mathbb{R}}^3,$$

$$f_{ij}({\mathbf{x}}_i,{\mathbf{x}}_j)=\frac{1}{2}{\mathbf{e}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathbf{e}}_{ij}=\frac{1}{2}\Vert {\text{Log}}_{\text{mathds}1}({\mathbf{z}}_{ij}^{-1}⊞{\mathbf{x}}_i^{-1}⊞{\mathbf{x}}_j){\Vert }_{{\mathbf{\Omega }}_{ij}}^2,$$

$${\tilde{\mathbf{x}}}^{k+1}={\mathcal{R}}_{{\tilde{\mathbf{x}}}^k}({\tilde{\mathbf{\alpha }}}^k),{\tilde{\mathbf{\alpha }}}^k\in \text{mathds}{T}_{{\tilde{\mathbf{x}}}^k}{\mathbb{M}}^n.$$

$$\text{Hess}F({\tilde{\mathbf{x}}}^k){\tilde{\mathbf{\alpha }}}^k=-\text{grad}F({\tilde{\mathbf{x}}}^k),$$

$$\text{grad}{f}_{ij}({\mathbf{x}}_i)={\mathcal{P}}_{{\mathbf{x}}_i}\nabla f_{ij}({\mathbf{x}}_i),$$

$$\text{grad}{f}_{ij}({\mathbf{x}}_i)={\mathcal{P}}_{{\mathbf{x}}_i}{\mathcal{A}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathbf{e}}_{ij}.$$

$$\text{grad}{f}_{ij}({\mathbf{x}}_j)={\mathcal{P}}_{{\mathbf{x}}_j}{\mathcal{B}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathbf{e}}_{ij},$$

$$\text{grad}{\mathcal{F}}_{[i]}$$

$$\text{grad}{F}_{[j]}$$

$${\mathcal{H}}_{[ii]}={\mathcal{P}}_{{\mathbf{x}}_i}{\mathcal{A}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathcal{A}}_{ij}{\mathcal{P}}_{{\mathbf{x}}_i}^{\top },$$

$${\mathcal{H}}_{[ij]}={\mathcal{P}}_{{\mathbf{x}}_i}{\mathcal{A}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathcal{B}}_{ij}{\mathcal{P}}_{{\mathbf{x}}_j}^{\top },$$

$${\mathcal{H}}_{[ji]}={\mathcal{P}}_{{\mathbf{x}}_j}{\mathcal{B}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathcal{A}}_{ij}{\mathcal{P}}_{{\mathbf{x}}_i}^{\top },$$

$${\mathcal{H}}_{[jj]}={\mathcal{P}}_{{\mathbf{x}}_j}{\mathcal{B}}_{ij}^{\top }{\mathbf{\Omega }}_{ij}{\mathcal{B}}_{ij}{\mathcal{P}}_{{\mathbf{x}}_j}^{\top },$$

$$\mathcal{H}\tilde{\mathbf{\alpha }}=-\text{grad}\mathcal{F},$$

Full text

††thanks: This work is supported by the German Research Foundation (DFG) under grant HA 3789/16-1. Johannes Cox is now with the Technical University of Berlin.

1 Introduction

Pose graph optimization (Grisetti et al. (2010)) plays a fundamental role in robotic and computer vision tasks such as simultaneous localization and mapping (SLAM), structure from motion (SfM), etc. A pose graph is a graph with nodes representing robot poses and edges linking the nodes, between which odometry measurements are available. Pose graph optimization problems are typically formulated to maximize the likelihood of the observed odometry measurements of all the edges with constraints imposed by the group structure. As pointed out in Carlone and Censi (2014), they are non-convex and contain multiple local minima. Mathematically speaking, the planar rigid motions belong to the special Euclidean group $\text{SE(}2\text{)}$. The frequently used representation approach in the context of planar pose graph optimization is the two-dimensional rotation matrix ${\bf R}\in\text{SO(}2\text{)}$ plus the translation ${\bf t}$ . Consequently, the three-DoF planar rigid motions are overparameterized by six elements, which leads to three degrees of redundancy. The cost function for performing the maximum likelihood estimation (MLE) is normally formulated as the sum of individual rotational and translational error, which lacks consideration of the underlying group structure. This can be problematic when the rotation and translation errors are at different scales of uncertainty and a proper scaling ratio of the two error terms is not available. Moreover, for highly nonlinear and uncertain rigid motions, the structure-unaware formulations and conventional solvers for constrained optimization problems may have the risk of non-convergence and are prone to local minima.

The planar dual quaternions are an alternative approach for representing planar rigid motions and can be written as four-dimensional vectors with only one degree of redundancy. Recently, there have been works dedicated to geometry-aware stochastic filters for pose estimation using the (dual) quaternion representation. In Gilitschenski et al. (2016), the Bingham distribution was deployed to stochastically model the uncertainty of unit quaternions and corresponding stochastic filters were proposed for orientation estimation. In Srivatsan et al. (2018) and Li et al. (2018), dual quaternion-based filtering approaches were proposed for rigid body motion estimation. In Cheng et al. (2016), the pose graph optimization problem was re-parameterized by unit dual quaternions and better computational efficiency was shown. However, the deployed objective function heuristically removes entries that might cause singularities without theoretically sound adaptation to the manifold structure of unit dual quaternions. Moreover, the error term in the objective function was merely reformulated from the conventional homogeneous matrix-based scheme. Therefore, the resulting graph optimization scheme was not observed with improved accuracy and robustness, especially in the case of large odometry noise.

From the perspective of differential geometry, the set of planar dual quaternions forms a Riemannian manifold (do Carmo (1992); Busam et al. (2017)) equipped with the so-called Riemannian metric, which allows measuring on-manifold uncertainty in accordance with the manifold structure. Furthermore, optimization approaches proposed on Riemannian manifolds, or equivalently Riemannian optimization approaches, have been gaining in popularity over the years (Absil et al. (2009)). Compared to conventional solvers for constrained optimization problems that use, e.g., the Lagrange multipliers, the family of Riemannian optimizers have shown better convergence efficiency and robustness as they exploit the geometric structure of the manifold.

In this paper, we propose the so-called Riemannian pose graph optimizer (RPG-Opt) for planar motions represented by dual quaternions. Here, odometry errors of edges are measured based on the Riemannian metric on the manifold of planar dual quaternions. Unlike conventional solvers for constrained optimization problems as given in Carlone et al. (2014, 2015a), we apply the Riemannian Gauss–Newton approach on the manifold with exponential retraction for updates. Both the cost function formulation and optimization approach are geometry-aware and inherently consider the underlying nonlinear structure of the SE(2) group. For evaluation, we compare the proposed approach with state-of-the-art pose graph optimization frameworks (e.g., GTSAM Dellaert (2012), g2o Kümmerle et al. (2011) and iSAM Kaess et al. (2008)) based on public planar pose graph data sets under ordinary odometry noise level. We further synthesize new data sets with additionally larger odometry noise for evaluating the robustness. Comparisons show that our approach gives equally accurate results as the state-of-the-art frameworks and gives superior accuracy under large odometry uncertainties.

The remainder of the paper is as follows. In Sec. 2, we introduce the dual quaternion parameterization for planar rigid motions and the geometric structure of the planar dual quaternion manifold. Sec. 3 introduces the proposed cost function and the on-manifold optimizer, the RPG-opt. We further evaluate the proposed approach regarding accuracy and convergence robustness in Sec. 4. The work is finally concluded in Sec. 5.

2 Planar Dual Quaternion Representation and Manifold Structure

2.1 Planar Rigid Motion Represented by Dual Quaternions

By convention, unit quaternions representing planar rotations are written as ${\bf r}=\cos(\theta/2)+{\bf k}\sin(\theta/2)$ ,with the unit vector ${\bf k}$ indicating the $z$-axis, around which a rotation of angle $\theta$ is performed (Brookshire and Teller (2013)). Any ${\bf v}\in\mathbb{R}^{2}$ can be rotated to ${\bf v}^{\prime}$ according to ${\bf r}$ via ${\bf v}^{\prime}={\bf r}\otimes{\bf v}\otimes{\bf r}^{*}$ , with $\otimes$ denoting the Hamilton product (Hamilton (1844)) and ${\bf r}^{*}=\cos(\theta/2)-{\bf k}\sin(\theta/2)$ the conjugate of ${\bf r}$ . Moreover, the planar rotation quaternions can be reformulated into the following vector form

[TABLE]

which are located on the unit circle on the $x,y$-plane. Given two planar quaternions ${\bf r}=[\,r_{0},r_{1}\,]^{\top}$ and ${\bf s}=[\,s_{0},s_{1}\,]^{\top}$ , their Hamilton product can also be represented as ordinary matrix-vector multiplication, namely ${\bf r}\otimes{\bf s}=\mathcal{Q}_{\text{r}}^{\llcorner}\,{\bf s}=\mathcal{Q}_{\text{s}}^{\lrcorner}\,{\bf r}$ , with

[TABLE]

It is then trivial to confirm that both the left and right matrix representation belong to the two-dimensional rotation group $\text{SO(}2\text{)}$, i.e., $\mathcal{Q}\,\mathcal{Q}^{\top}=\mathcal{Q}^{\top}\mathcal{Q}={\bf I}\in\mathbb{R}^{2\times 2}$ and its determinant $\det(\mathcal{Q})=1$ . The dual quaternion representation for planar rigid motions is defined as ${\bf x}={\bf r}+\frac{\epsilon}{2}{\bf t}\otimes{\bf r}\,,$ with $\epsilon$ denoting the dual number which satisfies $\epsilon^{2}=0$ and ${\bf t}$ the translation. The corresponding vector form of planar dual quaternions can then be written as

[TABLE]

with the real part ${\bf x}_{\text{r}}$ defined as in (1) and the dual part

[TABLE]

Therefore, the manifold of planar dual quaternions (denoted as $\mathbb{M}$) is the Cartesian product of the unit circle $\mathbb{S}^{1}$ and the two-dimensional Euclidean space $\mathbb{R}^{2}$. The arithmetics of dual quaternions are the combination of the Hamilton product and the dual number theory. For instance, given two planar dual quaternions ${\bf x}=[\,x_{0},\,x_{1},\,x_{2},\,x_{3}\,]^{\top}$ and ${\bf y}=[\,y_{0},\,y_{1},\,y_{2},\,y_{3}\,]^{\top}$ , their product in the form of matrix-vector multiplication can be written as ${\bf x}\boxplus{\bf y}=\mathcal{Q}_{\text{x}}^{\ulcorner}\,{\bf y}=\mathcal{Q}_{\text{y}}^{\urcorner}\,{\bf x}$ , with

[TABLE]

Similar to the rotation rule of unit quaternions, any ${\bf v}\in\mathbb{R}^{2}$ can be transformed to ${\bf v}^{\prime}$ via a planar rotation of $\theta$ followed by a translation ${\bf t}$ according to

[TABLE]

with ${\bf x}^{*}={\rm diag}\,(1,-1,-1,-1)\cdot{\bf x}$ being the conjugate of ${\bf x}$ . Here, ${\rm diag}(\cdot)$ is a diagonal matrix with the arguments placed at the diagonal entries. It should be noted that the vector form of dual quaternions for planar motions is expressed w.r.t. the coordinate system $\{1,{\bf k},{\bf i},{\bf j}\}$ , which is a reordering of the general quaternion basis $\{1,{\bf i},{\bf j},{\bf k}\}$ .

2.2 Geometry of Planar Dual Quaternion Manifold

The manifold of planar dual quaternions $\mathbb{M}$ can be derived as the zero-level set of the following vector function defined in $\mathbb{R}^{4}$

[TABLE]

with ${\bf 0}\in\mathbb{R}^{2}$ being a zero vector. This can be further used to derive the tangent plane $\mathds{T}_{\bf x}\mathbb{M}$ at ${\bf x}\in\mathbb{M}$ by calculating the null space of the Jacobian of $g({\bf x})$ evaluated at ${\bf x}$, i.e.,

[TABLE]

For any ${\bf y}\in\mathbb{R}^{4}$, its orthogonal projection to $\mathds{T}_{\bf x}\mathbb{M}$ can be obtained via $\mathcal{P}_{\bf x}{\bf y}$, with $\mathcal{P}_{\bf x}$ being the projection matrix evaluated at ${\bf x}\in\mathbb{M}$ . As there only exists the unit norm constraint for the real part, the projection matrix can be derived as the following form as given in Absil and Malick (2012)

[TABLE]

with ${\bf I},{\bf 0}\in\mathbb{R}^{2\times 2}$ . Note that the projection matrix is symmetric and idempotent, namely $\mathcal{P}_{\bf x}\mathcal{P}_{\bf x}=\mathcal{P}_{\bf x}$ .

The logarithm map of unit dual quaternions representing the $\text{SE(}3\text{)}$ states can be obtained via the reparameterization into screw motions of $(\theta,d,{\bf l},{\bf m})$ according to the Lie algebra as introduced in Brookshire and Teller (2013) as well as Busam et al. (2017). Here, $\theta$ denotes the screw angle (identical to the rotation angle), $d$ the translation along the screw axis ${\bf l}$ (identical to rotation axis) and ${\bf m}$ the screw moment. The logarithm map at the identity dual quaternion $\mathds{1}$ (a quaternion representing zero rotation and zero translation) is given as $\text{Log}_{\mathds{1}}({\bf x})=\frac{\theta}{2}{\bf l}+\frac{\epsilon}{2}(\theta\,{\bf m}+d\,{\bf l})$ , with $d={\bf t}^{\top}{\bf l}$ . The screw moment ${\bf m}$ is computed w.r.t. the projected origin ${\bf c}$ on the screw axis, namely ${\bf m}={\bf c}\times{\bf l}$ . Consequently, the logarithm map of dual quaternions representing planar rigid motions can be derived as a degenerate case of the spatial motions. Therefore, the projected point ${\bf c}$ is the intersection of the screw axis ${\bf k}$ with the $x,y$-plane and can be computed according to Busam et al. (2017) as ${\bf c}=\frac{1}{2}({\bf t}+{\bf k}\times{\bf t}\cot\frac{\theta}{2})$ .

Considering that the translation $d$ along the screw axis ${\bf k}$ is zero, we have the logarithm map of planar dual quaternion ${\bf x}=[\,x_{0},x_{1},x_{2},x_{3}\,]^{\top}$ derived as $\text{Log}_{\mathds{1}}({\bf x})=\frac{1}{2}\,[\,\theta,\,\theta\,{\bf m}^{\top}\,]^{\top}$. As a result, it can be further derived that

[TABLE]

According to the definition of the dual part in (3), this can be further derived into a more concise form as

[TABLE]

Therefore, the logarithm of planar dual quaternion can be written as

[TABLE]

with $\gamma=\text{sinc}(0.5\,\theta)$ and $\mathds{1}=[\,1,0,0,0\,]^{\top}$. Moreover, the logarithm map to the tangent plane at ${\bf x}\in\mathbb{M}$ can be computed according to the parallel transport (Busam et al. (2017)) on the manifold of planar dual quaternions as

[TABLE]

Conversely, the exponential map at identity $\mathds{1}$ maps a point ${\bf x}_{\text{t}}=[\,x_{\text{t},1},x_{\text{t},2},x_{\text{t},3}\,]^{\top}$ from the tangent plane at $\mathds{1}$ back to the manifold via

[TABLE]

Here, $\gamma=\text{sinc}(0.5\,\theta)=\text{sinc}(x_{\text{t},1})$ as given in (5). Similar to the logarithm map, the exponential map for arbitrary tangent plane locations can be derived according to the parallel transport

[TABLE]

Here, ${\bf y}_{\text{t}}\in\mathds{T}_{\bf x}\mathbb{M}$ and we take out the last three indices of ${\bf x}^{-1}\boxplus{\bf y}_{\text{t}}$ as its first element is zero according to (6).

3 RPG-Opt:

A Riemannian Pose Graph Optimizer

3.1 Geometry-Aware Cost Function

The pose graph optimization problem is formulated as the maximum likelihood estimation given the observed odometry. The optimized poses on the graph $\mathds{C}$ can then be obtained by minimizing the sum of the distance metrics through each edge of the graph, namely

[TABLE]

with $\tilde{{\bf x}}\in\mathbb{R}^{4n}$ being the poses of $n$ graph nodes concatenated into one vector, where each pose ${\bf x}_{i}\in\mathbb{M}\subset\mathbb{R}^{4}$ is a planar dual quaternion. Here, $\bm{\Omega}_{ij}\in\mathbb{R}^{3\times 3}$ denotes the uncertainty of the odometry measurement in the form of information matrix. The manifold for optimization is thus the Cartesian product of the planar dual quaternion manifold $\mathbb{M}$, which is also a Riemannian manifold according to Birdal et al. (2018). The cost function

[TABLE]

is thus a scalar function proposed on the manifold $\mathbb{M}^{n}$, namely $\mathcal{F}:\mathbb{M}^{n}\rightarrow\mathbb{R}$ . Unlike the existing works in Kaess and Dellaert (2009), Dellaert and Kaess (2006) and Kümmerle et al. (2011), we propose a cost function based on the Riemannian metric (Absil et al. (2009); Arvanitidis et al. (2016)) on the manifold of planar dual quaternions, which inherently considers the geometric structure of the nonlinear manifold. Given the odometry measurement ${\bf z}_{ij}$ of the edge connecting node ${\bf x}_{i}$ and ${\bf x}_{j}$, the error term of the edge is defined as

[TABLE]

such that the cost function for a single edge can be derived according to the Riemannian metric as

[TABLE]

which is the Mahalanobis distance measured on the tangent plane at the identity planar dual quaternion $\mathds{1}$ .

3.2 Riemannian Gauss–Newton Method for On-Manifold Planar Pose Graph Optimization

We apply the Riemannian Gauss–Newton approach in Absil et al. (2009) for solving the nonlinear least square problem formulated in (8). Here, a new iteration $\tilde{{\bf x}}^{k+1}$ is obtained by retracting the on-tangent-plane Newton step back to the manifold, namely

[TABLE]

The iterative step $\bm{\tilde{\bm{\alpha}}}^{k}$ results from the Riemanian gradient $\texttt{grad}\,F(\tilde{{\bf x}}^{k})$ and the Riemannian Hessian $\texttt{Hess}\,F(\tilde{{\bf x}}^{k})$ via

[TABLE]

with $\bm{\tilde{\bm{\alpha}}}^{k}$ being the concatenated on-tangent-plane steps calculated for each node. Without loss of generality, we perform the following derivations based on the cost function $f_{ij}$ of one single edge joining ${\bf x}_{i}$ and ${\bf x}_{j}$. For better readability, we ignore the iteration index $k$ . The Riemannian gradient w.r.t. ${\bf x}_{i}$ can be computed by orthogonally projecting the ordinary gradient of $f_{ij}$ onto the tangent plane at ${\bf x}_{i}$ via

[TABLE]

with $\mathcal{P}_{{\bf x}_{i}}$ being the projection matrix to $\mathds{T}_{{\bf x}_{i}}\mathbb{M}$ given in (4) and $\nabla f_{ij}({\bf x}_{i})=\mathcal{A}_{ij}^{\top}\,\bm{\Omega}_{ij}\,{\bf e}_{ij}$ the classical gradient in the Euclidean space. Here, $\mathcal{A}_{ij}$ is the classical Jacobian of the error metric ${\bf e}_{ij}$ w.r.t. ${\bf x}_{i}$ . Therefore, we have

[TABLE]

Similarly, the Riemannian gradient at node ${\bf x}_{j}$ can be calculated as

[TABLE]

with $\mathcal{B}_{ij}$ denoting the classical Jacobian of ${\bf e}_{ij}$ w.r.t. ${\bf x}_{j}$ . The Riemannian gradient of $\mathcal{F}$ in (9) is then computed by traversing all the nodes through the graph with each entry computed as

[TABLE]

The Riemannian Hessian can be computed by approximation (Absil et al. (2009)) with entries as follows

[TABLE]

Here, $i$ and $j$ indicate the locations of block matrices in the Hessian matrix $\mathcal{H}$ corresponding to the nodes ${\bf x}_{i}$ and ${\bf x}_{j}$ . Finally, the iterative step can be obtained by solving the linear system

[TABLE]

with $\tilde{\bm{\alpha}}$ denoting the on-tangent-plane iteration concatenated through each node.

The retraction maps a point from the tangent plane back to the manifold, namely $\mathcal{R}_{\tilde{{\bf x}}}:\mathds{T}_{\tilde{{\bf x}}}\mathbb{M}^{n}\rightarrow\mathbb{M}^{n}$ , which is not unique. For instance, there have been different projection-like retractions proposed in Absil and Malick (2012). We use the exponential retraction for $\tilde{{\bf x}}$ with each node undergoing the exponential map defined in (7) for update, namely

[TABLE]

for $i=1,\,...\,,n$ . After the retraction back to the manifold, the cost function can then be linearized at the tangent plane at $\tilde{{\bf x}}^{k+1}$. The optimization is stopped until the norm of the Riemannian gradient becomes sufficiently small.

4 Evaluation

We evaluate the proposed RPG-Opt based on publicly available planar pose graph data sets given in Carlone and Censi (2014) and Carlone et al. (2014), and compare it with three state-of-the-art pose graph optimization frameworks: g2o in Grisetti et al. (2010), GTSAM in Dellaert (2012) and iSAM in Kaess et al. (2008). We choose to use Levenberg–Marquardt algorithm for iSAM and the Gauss–Newton method for GTSAM. In g2o, the optimization is performed using the Gauss–Newton algorithm. For most of the publicly available pose graph data sets, the uncertainty of the odometry measurements is described by the information matrix for the state $[\,\theta,\,{\bf t}_{\text{x}},\,{\bf t}_{\text{y}}\,]^{\top}$ , with $\theta$ being the rotation angle and $[\,{\bf t}_{\text{x}},\,{\bf t}_{\text{y}}\,]^{\top}$ the translation. In order to incorporate the estimated uncertainty, the logarithm map in (5) can be reformulated into the following form

[TABLE]

with $\alpha=\delta\theta/2$ and $\beta=\cos(\delta\theta/2)/\text{sinc}(\delta\theta/2)$ . Here, $\delta\theta$ and $[\delta{\bf t}_{\text{x}}\,,\delta{\bf t}_{\text{y}}]^{\top}$ is the rotation angle and translation term of the planar dual quaternion ${\bf z}_{ij}^{-1}\boxplus{\bf x}_{i}^{-1}\boxplus{\bf x}_{j}$ , respectively. The equation in (11) is essentially a nonlinear transformation of the state $[\delta\theta,\delta{\bf t}_{\text{x}},\delta{\bf t}_{\text{y}}]^{\top}$ provided by the data set. For each iteration step, we can assume that the odometry error is small, meaning $\alpha\rightarrow 0$ and $\beta\rightarrow 1$ . Therefore, the information matrix from the raw data set can be directly deployed as the on-tangent-plane information matrix $\bm{\Omega}_{ij}$ in the metric of (11) . When the information matrices for odometry are not available, a typical way to formulate the cost function is to set $\bm{\Omega}_{ij}={\bf I}\in\mathbb{R}^{3\times 3}$ .

In the following evaluations, we first test the proposed approach based on data sets under ordinary odometry noise level using the identity and real information matrix, respectively. Second, we synthesize data sets with additional noise and perform the evaluation w.r.t. the ground truth.

4.1 Evaluation Under Ordinary Odometry Noise Level

We first evaluate the proposed approach using data sets with ordinary level of odometry noise. As no ground truth is available, we use the cost function of g2o for comparison (Kümmerle et al. (2011)), namely

[TABLE]

with ${\bf R}\in\text{SO(}2\text{)}$ being the ordinary rotation matrices and ${\bf t}$ the translations.

Here, the logarithm map of rotation matrices can be derived from Lie algebra as given in Carlone et al. (2015b). Table 1 shows the results for both identity and real information matrices (denoted as ${\bf I}$ and $\bm{\Omega}$, respectively). Here, g2o and RPG-Opt iteration steps are fixed to be $10$ , whereas GTSAM and iSAM use their default stopping criterion. Our approach, the RPG-Opt, reaches comparable accuracy as the state of the art.

4.2 Evaluation Based on Synthetic Data sets

We further synthesize two groups of data sets with known information matrices for the evaluation under large uncertainty of odometry measurements. For that, we add additional odometry noise to the data sets M3500s and City10000 in Carlone and Censi (2014). Fig. 1 shows the results based on the relative pose error (RPE) for both translations and rotations, which are denoted as $e_{\text{t}}$ and $e_{\text{r}}$, respectively. The optimizations are initialized directly from the odometry measurements. Results from GTSAM using the Gauss–Newton approach are not listed because of nonconvergence. For all the sequences, the g2o and iSAM framework are prone to the local minima. However, the proposed approach shows the best accuracy and robustness against local minima, though no special initialization is performed.

For pose graph optimization under large uncertainty of odometry measurements, it is typical to equip the optimizer with an additional initialization block for better convergence. Therefore, we incorporate an initialization method based on chordal relaxation in Carlone et al. (2015b). Fig. 2 visualizes evaluation results with chordal relaxation-based initialization using another two synthetic data sets based on the M3500 and City10000. In this case, the g2o shows significant improvement, achieving the same converged accuracy as the proposed approach. As a summary, Table 2 shows the optimization results with chordal initialization for all the synthetic data sets w.r.t. the ground truth. For data sets M3500a+, M3500b+, M3500c+ and City10000a, the proposed approach initialized with odometry measurements achieves the same errors as the other approaches using the chordal initialization. Results from iSAM are not available, because external initializations are not supported by the framework. We additionally list the results with direct odometry initialization for the RPG-Opt to show the convergence robustness under large uncertainties of the proposed method .

4.3 Convergence

We further compare the convergence behavior of the proposed approach with g2o as shown in Fig. 3. Here, the sequence M3500c+ is used and the RPE for both translation and rotations are plotted in Fig. 3-(a). Fig. 3-(b) shows the convergence of the proposed Riemannian metric in (9). The optimizations are initialized with odometry and the real information matrices are used. In both plots, the proposed approach shows faster convergence and better robustness against local minimum compared to g2o. This mainly results from the geometry-aware cost function formulation and the Riemannian Gauss–Newton of the proposed approach.

5 Conclusion

In this work, we have proposed a Riemannian approach for planar pose graph optimization problems on the manifold of dual quaternions. Here, the cost function is built upon the Riemannian metric of the planar dual quaternion manifold, based on which a Riemannian Gauss–Newton method is applied using the exponential retraction. Both the on-manifold MLE formulation and optimization are geometry-aware, which inherently consider the underlying nonlinear structure of the $\text{SE(}2\text{)}$ group. We further performed evaluations based on real-world and synthetic data sets with extra noise. Compared with the state-of-the-art frameworks, the proposed approach gives equally accurate results under ordinary odometry noise and shows better accuracy and robustness under large uncertainty of odometry measurements.

Based on the presented work, there is still much potential to exploit. The proposed Riemannian approach can be extended to the general dual quaternion-based spatial pose graph optimization. Also, numerous optimizers from the family of Riemannian optimization can be further applied for better convergence robustness and accuracy.