Rotation Averaging and Strong Duality

TL;DR

This paper leverages Lagrangian duality and spectral graph theory to solve rotation averaging problems efficiently, providing certifiably global solutions with no duality gap in many cases, and demonstrates superior scalability and performance.

Contribution

It introduces a duality-based approach that guarantees global optimality in rotation averaging, with an efficient algorithm applicable to large-scale problems.

Findings

No duality gap in many cases unless noise is severe

Proposed algorithm outperforms general solvers in scalability

Successfully applied to synthetic and real-world data

Abstract

In this paper we explore the role of duality principles within the problem of rotation averaging, a fundamental task in a wide range of computer vision applications. In its conventional form, rotation averaging is stated as a minimization over multiple rotation constraints. As these constraints are non-convex, this problem is generally considered challenging to solve globally. We show how to circumvent this difficulty through the use of Lagrangian duality. While such an approach is well-known it is normally not guaranteed to provide a tight relaxation. Based on spectral graph theory, we analytically prove that in many cases there is no duality gap unless the noise levels are severe. This allows us to obtain certifiably global solutions to a class of important non-convex problems in polynomial time. We also propose an efficient, scalable algorithm that out-performs general purpose numerical solvers and is able to handle the large problem instances commonly occurring in structure from motion settings. The potential of this proposed method is demonstrated on a number of different problems, consisting of both synthetic and real-world data.

Tables

Table 1. Table 1: Comparison of running times and resulting errors on synthetic data. Here the errors are given with respect to the lowest feasible objective function value found. The fraction of the times the global optima was reached by the LM algorithm is indicated along side the average error.

	LM [29]			Alg. 1			SeDuMi [24]
$n$         $\sigma$ [rad]	$avg.error$ (%)	$time\left[s\right]$		$avg.error$	$time\left[s\right]$		$avg.error$	$time\left[s\right]$
$20$              $0.2$	1.49 (0.48)	0.012		9.34e-10	0.028		4.30e-09	0.501
$0.5$	0.56 (0.73)	0.008		3.94e-08	0.023		3.72e-09	0.553
$50$              $0.2$	0.55 (0.50)	0.026		1.3e-09	0.17		6.85e-09	5.91
$0.5$	0.17 (0.58)	0.017		1.83e-07	0.33		2.00e-09	6.32
$100$              $0.2$	0.15 (0.55)	0.042		1.46e-07	8.89		5.31e-09	47.0
$0.5$	0.15 (0.45)	0.039		6.64e-08	7.97		7.41e-10	49.51
$200$              $0.2$	0.099 (0.40)	0.082		4.02e-08	17.01		4.15e-10	419.04
$0.5$	0.031 (0.33)	0.071		6.79e-08	29.4		6.91e-10	391.23

Table 2. Table 2: The average run time and largest resulting angular residual ( | α i ​ j | subscript 𝛼 𝑖 𝑗 |\alpha_{ij}| ) and bound ( α max subscript 𝛼 \alpha_{\max} ) on five different real-world datasets.

		$time\left[s\right]$
$Dataset$	$n$	Alg. 1	SeDuMi	$\left|{\alpha }_{ij}\right|$	${\alpha }_{\max }$
Gustavus	57	3.25	8.28	${6.33}^{\circ }$	${8.89}^{\circ }$
Sphinx	70	3.87	14.40	${6.14}^{\circ }$	${12.13}^{\circ }$
Alcatraz	133	12.73	117.19	${7.68}^{\circ }$	${43.15}^{\circ }$
Pumpkin	209	9.23	688.65	${8.63}^{\circ }$	${3.59}^{\circ }$
Buddha	322	16.71	1765.72	${7.29}^{\circ }$	${14.01}^{\circ }$

Equations

$$\displaystyle a_{ij}=\left\{\begin{array}[]{ccccccccccc}0&(i,j)\notin E\\ 1&(i,j)\in E\end{array}\mbox{ for }i,j=1,\ldots,n.\right.$$

$$L_G=D-A.$$

$$\mathrm{SO}(3)\in \{R\in {\mathbb{R}}^{3\times 3}|R^{T}R=I,\det (R)=1\}.$$

$$R_i{\tilde{R}}_{ij}=R_j,$$

$$\underset{R_1,...,R_n}{\min }\sum _{(i,j)\in E}d(R_i{\tilde{R}}_{ij},R_j{)}^p,$$

$$d(R,S)=\Vert R-S{\Vert }_F.$$

$$\underset{R_1,...,R_n\in \mathrm{SO}(3)}{arg\; min}\sum _{(i,j)\in E}\Vert R_i{\tilde{R}}_{ij}-R_j{\Vert }_F^2,$$

$$\underset{R_1,...,R_n\in \mathrm{SO}(3)}{arg\; min}-\sum _{(i,j)\in E}\text{tr}\left(R_i{\tilde{R}}_{ij}{R}_j^T\right),$$

$$\tilde{R}=\left[\begin{array}{cccc}0& a_{12}{\tilde{R}}_{12}& \dots & a_{1n}{\tilde{R}}_{1n}\\ a_{21}{\tilde{R}}_{21}& 0& \dots & a_{2n}{\tilde{R}}_{2n}\\ ⋮& & \ddots & ⋮\\ a_{n1}{\tilde{R}}_{n1}& a_{n2}{\tilde{R}}_{n2}& \dots & 0\end{array}\right],$$

$$R=\left[\begin{array}{cccc}{R}_1& R_2& \dots & R_n\end{array}\right].$$

$$\begin{array}[]{lll}(P)&\min&-\text{tr}\left(R\tilde{R}R^{T}\right)\\ &\text{s.t.}&R\in{\rm SO}(3)^{n}.\end{array}$$

$$\begin{array}[]{lll}(P^{\prime})&\min&-\text{tr}\left(R\tilde{R}R^{T}\right)\\ &\text{s.t.}&R\in{\rm O}(3)^{n}.\end{array}$$

$$\Lambda =\left[\begin{array}{cccc}{\Lambda }_1& 0& 0& \dots \\ 0& {\Lambda }_2& 0& \dots \\ 0& 0& {\Lambda }_3& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

$$\begin{array}[]{rl}L(R,\Lambda)&=-\text{tr}\left(R\tilde{R}R^{T}\right)-\text{tr}\left(\Lambda(I-R^{T}R)\right)\\ &=-\text{tr}\left(R(\Lambda-\tilde{R})R^{T}\right)-\text{tr}\left(\Lambda\right).\end{array}$$

$$({\Lambda }^{\ast }-\tilde{R})R^{{\ast }^T}=0$$

$$R^{\ast }\in \mathrm{SO}(3{)}^n.$$

$${\Lambda }_i^{\ast }{R}_i^{\ast T}=\sum _{j\ne i}{a}_{ij}{\tilde{R}}_{ij}{R}_j^{\ast T}⟺{\Lambda }_i^{\ast }=\sum _{j\ne i}{a}_{ij}{\tilde{R}}_{ij}{R}_j^{\ast T}{R}_i^{\ast }$$

$$\underset{\Lambda -\tilde{R}⪰0}{\max }\underset{R}{\min }L(R,\Lambda ).$$

$$(D)\underset{\Lambda -\tilde{R}⪰0}{\max }-\text{tr}\left(\Lambda \right).$$

$$x^T\left({\Lambda }^{\ast }-\tilde{R}\right)x\ge 0,$$

$$\begin{array}[]{rcl}0&\leq&\text{tr}\left(R(\Lambda^{*}-\tilde{R})R^{T}\right)=\text{tr}\left(\Lambda^{*}\right)-\text{tr}\left(R\tilde{R}R^{T}\right)\\ &=&\text{tr}\left(R^{*}\Lambda^{*}R^{*T}\right)-\text{tr}\left(R\tilde{R}R^{T}\right),\end{array}$$

$$|{\alpha }_{ij}|\le {\alpha }_{\max }\forall (i,j)\in E,$$

$${\alpha }_{\max }=2\arcsin \left(\sqrt{\frac{1}{4}+\frac{{\lambda }_2(L_G)}{2d_{\max }}}-\frac{1}{2}\right),$$

$$D_{R^{\ast }}=\left[\begin{array}{cccc}{R}_1^{\ast }& 0& 0& \dots \\ 0& R_2^{\ast }& 0& \dots \\ 0& 0& R_3^{\ast }& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right].$$

$$\left[\begin{array}{cccc}{\sum }_{j\ne 1}{a}_{1j}{\mathcal{E}}_{1j}& -a_{12}{\mathcal{E}}_{12}& -a_{13}{\mathcal{E}}_{13}& \dots \\ -a_{12}{\mathcal{E}}_{12}^T& {\sum }_{j\ne 2}{a}_{2j}{\mathcal{E}}_{2j}& -a_{23}{\mathcal{E}}_{23}& \dots \\ -a_{13}{\mathcal{E}}_{13}^T& -a_{23}{\mathcal{E}}_{23}^T& {\sum }_{j\ne 3}{a}_{3j}{\mathcal{E}}_{3j}& \dots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right].$$

$$D_{R^{\ast }}({\Lambda }^{\ast }-\tilde{R})D_{R^{\ast }}^T=L_G\otimes I_3.$$

$$\Delta =D_{R^{\ast }}({\Lambda }^{\ast }-\tilde{R})D_{R^{\ast }}^T-L_G\otimes I_3.$$

$$|\lambda |\le \sum _{j=1}^n\Vert {\Delta }_{ij}\Vert \text{mbox}forsomei=1,\dots ,n.$$

$$|\lambda |=\Vert \lambda \frac{x_i}{\Vert x_i\Vert }\Vert =\Vert \sum _{j=1}^n{\Delta }_{ij}\frac{x_j}{\Vert x_i\Vert }\Vert \le \sum _{j=1}^n\Vert {\Delta }_{ij}\Vert .$$

$$\Vert {\Delta }_{ii}\Vert \le 2d_i{\sin }^2(\frac{{\alpha }_{\max }}{2})\forall i=1,\dots n,$$

$${\mathcal{E}}_{ij}=V_{ij}\left[\begin{array}{ccc}\cos ({\alpha }_{ij})& -\sin ({\alpha }_{ij})& 0\\ \sin ({\alpha }_{ij})& \cos ({\alpha }_{ij})& 0\\ 0& 0& 1\end{array}\right]V_{ij}^T,$$

Full text

We also propose an efficient, scalable algorithm that out-performs general purpose numerical solvers and is able to handle the large problem instances commonly occurring in structure from motion settings. The potential of this proposed method is demonstrated on a number of different problems, consisting of both synthetic and real-world data.

1 Introduction

Rotation averaging appears as a subproblem in many important applications in computer vision, robotics, sensor networks and related areas. Given a number of relative rotation estimates between pairs of poses, the goal is to compute absolute camera orientations with respect to some common coordinate system. In computer vision, for instance, non-sequential structure from motion systems such as [21, 11, 22] rely on rotation averaging to initialize bundle adjustment. The overall idea is to consider as much data as possible in each step to avoid suboptimal reconstructions. In the context of rotation averaging this amounts to using as many camera pairs as possible.

The problem can be thought of as inference on the camera graph. An edge $(i,j)$ in this undirected graph represents a relative rotation measurement $\tilde{R}_{ij}$ and the objective is to find the absolute orientation $R_{i}$ for each vertex $i$ such that $R_{i}\tilde{R}_{ij}=R_{j}$ holds (approximately in the presence of noise) for all edges. The problem is generally considered difficult due to the need to enforce non-convex rotation constraints. Indeed, both $L_{1}$ and $L_{2}$ formulations of rotation averaging can have local minima, see Fig. 1. Wilson et al. [28] studied local convexity of the problem and showed that instances with large loosely connected graphs are hard to solve with local, iterative optimization methods.

In contrast, our focus is on global optimality. In this paper we show that convex relaxation methods can in fact overcome the difficulties with local minima in rotation averaging. We utilize Lagrangian duality to handle the quadratic non-convex rotation constraints. While such an approach is normally not guaranteed to provide a tight relaxation we give analytical error bounds that guarantee there will be no duality gap. For instance, it is sufficient that each angular residual is less than $42.9^{\circ}$ to ensure optimality for complete camera graphs. Additionally, we develop a scalable and efficient algorithm, based on block coordinate descent, that outperforms standard semidefinite program (SDP) solvers for this problem.

Related work.

Rotation averaging has been under intense study in recent years, see [19, 20, 21, 2, 25, 8]. Despite progress in practical algorithms, they largely come without guarantees. One of the earliest averaging methods was due to Govindu [15], who showed that when representing the rotations with quaternions the problem can be viewed as a linear homogeneous least squares problem. There is however a sign ambiguity in the quaternion representation that has to be resolved before the formulation can be applied. It was observed by Fredriksson and Olsson in [14] that since both the objective and the constraints are quadratic, the Lagrange dual can be computed in closed form. The resulting SDP was experimentally shown to have no duality gap for moderate noise levels.

A more straightforward rotation representation is $3\times 3$ matrices. Martinec and Pajdla [21] approximately solve the problem by ignoring the orthogonality and determinant constraints. A similar relaxation was derived by Arie-Nachimson et al. in [1]. In addition, an SDP formulation was presented which is equivalent to the one we address here, but with no performance guarantees. The tightness of SDP relaxations for 2D rotation averaging is studied in [30].

A number of robust approaches have been developed to handle outlier measurements. A sampling scheme over spanning trees of the camera graph is developed by Govindu in [16]. Enqvist et al. [11] also start from a spanning tree and add relative rotations that are consistent with the solution. In [17] the Weiszfeld algorithm is applied to single rotation averaging with the $L_{1}$ norm. In [18] convexity properties of the single rotation averaging problem are given. To our knowledge these results do not generalize to the case of multiple rotations. In [9] a robust formulation is solved using IRLS and in [3] Cramér-Rao lower bounds are computed for maximum likelihood estimators, but neither with any optimality guarantees.

A closely related problem is that of pose graph estimation, where camera orientations and positions are jointly optimized. In this context Lagrangian duality has been applied [6, 7]. In [26] a consensus algorithm that allows for efficient distributed computations is presented. A fast verification technique for pose graph estimation was given in [5]. In a recent paper [23] an SDP relaxation for pose graph estimation with performance guarantees is analyzed. It is shown that there is a noise level $\beta$ for which the relaxation is guaranteed to provide the optimal solution. However, the result only shows the existence of $\beta$. Its value which is dependent on the problem instance is not computed. In contrast our result for rotation averaging gives explicit noise bounds.

The main contributions of this paper are:

• •

We apply Lagrangian duality to the rotation averaging problem with the chordal error distance and study the properties of the obtained relaxations.

• •

We develop strong theoretical bounds on the noise level that guarantee exact global recovery based on spectral graph theory.

• •

We develop a conceptually simple and scalable algorithm which is able to handle large problem instances occurring in structure from motion problems.

• •

We present experimental results that confirm our theoretical findings.

1.1 Notation and Conventions

Let $G=(V,E)$ denote an undirected graph with vertex set $V$ and edge set $E$ and let $n=|V|$. The adjacency matrix $A$ is by definition the $n\times n$ matrix with elements

[TABLE]

The degree $d_{i}$ is the number of edges that touch vertex $i$, and the degree matrix $D$ is the diagonal matrix $D=\text{diag}\left(d_{1},\ldots,d_{n}\right)$. The Laplacian $L_{G}$ of $G$ is defined by

[TABLE]

It is well-known that $L_{G}$ has a zero eigenvalue with multiplicity $1$. The second smallest eigenvalue $\lambda_{2}$ of $L_{G}$, also known as the Fiedler value, reflects the connectivity of $G$. For a connected graph $G$, which is the only case of interest to us, we always have $\lambda_{2}>0$.

The group of all rotations about the origin in three dimensional Euclidean space is the Special Orthogonal Group, denoted ${\rm SO}(3)$. This group is commonly represented by rotation matrices, orthogonal $3\times 3$ real-valued matrices with positive determinant, i.e.,

[TABLE]

If we omit $\det(R)\!=\!1$, we get the Orthogonal Group, ${\rm O}(3)$.

We will use the convention that $\lambda_{i}(A)$ is the $i$:th smallest eigenvalue of the symmetric matrix $A$. The trace of matrix $A$ is denoted by $\text{tr}\left(A\right)$ and the Kronecker product of matrices $A$ and $B$ by $A\otimes B$. The norm $\|A\|$ is the standard operator 2-norm and $\|A\|_{F}$ the Frobenius norm.

2 Problem Statement

The problem of rotation averaging is defined as the task of determining a set of $n$ absolute rotations $R_{1},...,R_{n}$ given distinct estimated relative rotations $\tilde{R}_{ij}$. Available relative rotations are represented by the edge set $E$ of the camera graph $V$. Under ideal conditions this amounts to finding the $n$ rotations compatible with the linear relations,

[TABLE]

for all $(i,j)\in E$. However, in the presence of noise, a solution to (6) is not guaranteed to exist. Instead, it is typically solved in a least-metric sense,

[TABLE]

where $p\geq 1$ and $d(\cdot,\cdot)$ is a distance function.

A number of distinct choices of metrics on ${\rm SO}(3)$ exist, see Hartley et al. [19] for a comprehensive discussion. In this work we restrict ourselves to the chordal distance, the most commonly used metric when analyzing Lagrangian duality in rotation averaging. It has proven to be a convenient choice as it is quadratic in its entries leading to a particularly simple derivation and form of the associated dual problem.

The chordal distance between two rotations $R$ and $S$ is defined as their Euclidean distance in the embedding space,

[TABLE]

It can be shown [19] that the chordal distance can also be written as $d(R,S)=2\sqrt{2}\sin\frac{|\alpha|}{2}$, where $\alpha$ is the rotation angle of $RS^{-1}$. With the this choice of metric, the rotation averaging problem is defined as

[TABLE]

which, with trace notation, can be simplified to

[TABLE]

which constitutes our primal problem.

It will be convenient with a compact matrix formulation. Let

[TABLE]

where $\tilde{R}_{ij}=\tilde{R}_{ji}^{T}$ and $a_{ij}$ are the elements of the adjacency matrix $A$ of the camera graph $G$ and let

[TABLE]

We may now write the primal problem as

[TABLE]

3 Optimality Conditions

3.1 Necessary Local Optimality Conditions

We now turn to the KKT conditions of our primal problem $(P)$. The constraint set $R\in{\rm SO}(3)^{n}$ consists of two types of constraints; the orthogonality constraints $R_{i}^{T}R_{i}=I$ and the determinant constraints $\det(R_{i})=1$.

Consider relaxing the rotation averaging problem by removing the determinant constraint,

[TABLE]

The constraint $R\in{\rm O}(3)^{n}$ still requires the $R_{i}$’s to be orthogonal. The orthogonal matrices consist of two disjoint, non-connected sets, with determinants $1$ and $-1$ respectively. Hence, any local minimizer to the problem $(P)$ also has to be a local minimizer, and therefore a KKT point, to $(P^{\prime})$. We note that orthogonality can be enforced by restricting the $3\times 3$ diagonal blocks of the symmetric matrix $R^{T}R$ to be identity matrices. If

[TABLE]

is a symmetric matrix then the Lagrangian can be written

[TABLE]

Taking derivatives gives the KKT equations

[TABLE]

Equation (17a) states that the rows of a local minimizer $R^{*}$ will be eigenvectors of the matrix $\Lambda^{*}-\tilde{R}$ with eigenvalue zero. This allows us to compute the optimal Lagrange multiplier $\Lambda^{*}$ from a given minimizer $R^{*}$. By (17a) we see that

[TABLE]

for $i=1,\ldots,n$.

Lemma 3.1**.**

For a stationary point $R^{*}$ to the primal problem ($P$), we can compute the corresponding Lagrangian multiplier $\Lambda^{*}$ in closed form via (18).

3.2 Sufficient Global Optimality Conditions

We begin this section by deriving the Lagrange dual of ($P$) which is a semidefinite program that we will use for optimization in later sections. The dual problem is defined by

[TABLE]

Since the (unrestricted) optimum of $\min_{R}L(R,\Lambda)$ is either $-\text{tr}\left(\Lambda\right)$, when $\Lambda-\tilde{R}\succeq 0$, or $-\infty$ otherwise, we get

[TABLE]

It is clear (through standard duality arguments) that ($D$) gives a lower bound on ($P$). Furthermore, if $R^{*}$ is a stationary point with corresponding Lagrangian multiplier $\Lambda^{*}$ that satisfies $\Lambda^{*}-\tilde{R}\succeq 0$ then $\Lambda^{*}$ is feasible in ($D$) and by (18), $-\text{tr}\left(\Lambda^{*}\right)=-\text{tr}\left(R^{*}\tilde{R}{R^{*}}^{T}\right)$, which shows that there is no duality gap between ($P$) and ($D$). Thus, the convex program ($D$) provides a way of solving the non-convex ($P$) when $\Lambda^{*}-\tilde{R}\succeq 0$.

It also follows that for the stationary point $R^{*}$ we have $\text{tr}\left(R^{*}\Lambda^{*}R^{*T}\right)=\text{tr}\left(R^{*}\tilde{R}R^{*T}\right)$ due to $\eqref{eq:KKT}$. We further note that if $\Lambda^{*}-\tilde{R}\succeq 0$ then by definition it is true that

[TABLE]

for any $3n$-vector $x$. In particular, for any $R\in{\rm O}(3)^{n}$,

[TABLE]

which shows that $-\text{tr}\left(R^{*}\tilde{R}R^{*T}\right)\leq-\text{tr}\left(R\tilde{R}R^{T}\right)$ for all $R\in{\rm O}(3)^{n}$, that is, $R^{*}$ is the global optimum.

Lemma 3.2**.**

If a stationary point $R^{*}$ with corresponding Lagrangian multiplier $\Lambda^{*}$ fulfills $\Lambda^{*}-\tilde{R}\succeq 0$ then:

There is no duality gap between ($P$) and ($D$). 2. 2.

$R^{*}$* is a global minimum for ($P$).*

In the remainder of this paper we will study under which conditions $\Lambda^{*}-\tilde{R}\succeq 0$ holds and derive an efficient implementation for solving ($D$).

4 Main Result

In this section, we will state our main result which gives error bounds that guarentee that that strong duality holds for our primal and dual problems. From a practical point of view, the result means that it is possible to solve a convex semidefinite program and obtain the globally optimal solution to our non-convex problem, which is quite remarkable.

4.1 Strong Duality Theorem

Returning to our initial, primal rotation averaging problem (9). The goal is to find rotations $R_{i}$ and $R_{j}$ such that the sum of the residuals $\|R_{i}\tilde{R}_{ij}-R_{j}\|_{F}^{2}$ is minimized. For strong duality to hold, we need to bound the residual error.

Theorem 4.1** (Strong Duality).**

Let $R_{i}^{*}$, $i=1,\ldots,n$ denote a stationary point to the primal problem ($P$) for a connected camera graph $G$ with Laplacian $L_{G}$. Let $\alpha_{ij}$ denote the angular residuals, i.e., $\alpha_{ij}=\angle(R^{*}_{i}\tilde{R}_{ij},R^{*}_{j})$. Then $R_{i}^{*}$, $i=1,\ldots,n$ will be globally optimal and strong duality will hold for ($P$) if

[TABLE]

where

[TABLE]

and $d_{\max}$ is the maximal vertex degree.

Note that any local minimizer that fulfills this error bound will be global, and conversely there are no non-global minimizers with error residuals fulfilling (23). It is clear that (24) will give a positive bound $\alpha_{\max}$ for any graph. Thus for any given problem instance, $\alpha_{\max}$ gives an explicit bound on the error residuals for which strong duality is guaranteed to hold. The strength of the bound will depend on the particular graph connectivity encapsulated by the Fiedler value $\lambda_{2}(L_{G})$ and the maximal vertex degree $d_{\max}$. We will see that for tightly connected graphs the bound ensures strong duality under surprisingly generous noise levels. In [28] it was observed that local convexity at a point holds under similar circumstances.

Example.

Consider a graph with $n=3$ vertices that are connected, and all degrees are equal, $d_{max}=2$. Now from the Laplacian matrix $L_{G}$, one easily finds that $\lambda_{2}=3$. This gives $\alpha_{\max}=\frac{\pi}{3}\mathrm{rad}=60^{\circ}$. So, any local minimizer which has angular residuals less than $60^{\circ}$ is also a global solution.

Complete graphs.

Let us turn to a more general class of graphs, namely complete graphs with $n$ vertices, see Fig. 2. As every pair of vertices is connected, it follows that $d_{\max}=n-1$. Further, it is well-known (and easy to show) that $\lambda_{2}(L_{G})=n$, see [13]. Again, for $n=3$, we retrieve $\alpha_{\max}=\frac{\pi}{3}\mathrm{rad}$. As $n$ becomes larger, we get a decreasing series of upper bounds which in the limit tends to $2\arcsin(\frac{\sqrt{3}-1}{2})\approx 0.749\mathrm{rad}=42.9^{\circ}$. Hence, as long as the residual angular errors are less than $42.9^{\circ}$ - which is quite generous from a practical point of view - we can compute the optimal solution via a convex program. Also note that this bound holds independently of $n$.

Corollary 4.1**.**

For a complete graph $G$ with $n$ vertices, the residual upper bound $\alpha_{\max}=2\arcsin(\frac{\sqrt{3}-1}{2})\approx 0.749\mathrm{rad}=42.9^{\circ}$ ensures global optimality and strong duality for any $n$.

Cycle graphs.

Now consider the other spectrum in terms of graph connectivity, namely cycle graphs. A cycle graph has a single cycle, or in other words, every vertex in the camera graph has degree two ($d_{\max}=2$) and the vertices form a closed chain (Fig. 2). From the literature, we have that the Fiedler value $\lambda_{2}=2(1-\cos{\frac{2\pi}{n}})$. Inserting into (24) and simplifying, we get $\alpha_{\max}=2\arcsin\left(\sqrt{\frac{1}{4}+\sin^{2}(\frac{\pi}{n})}-\frac{1}{2}\right)$. Again, for $n=3$, we retrieve $\alpha_{\max}=\frac{\pi}{3}\mathrm{rad}$. For larger values of $n$, the upper bound decreases rapidly. In fact, the upper bound is quite conservative and it is possible to show a much stronger upper bound using a different analysis. In the appendix, we prove the following theorem.

Theorem 4.2**.**

Let $R_{i}^{*}$, $i=1,\ldots,n$ denote a stationary point to the primal problem ($P$) for a cycle graph with $n$ vertices. Let $\alpha_{ij}$ denote the angular residuals, i.e., $\alpha_{ij}=\angle(R^{*}_{i}\tilde{R}_{ij},R^{*}_{j})$. Then, $R_{i}^{*}$, $i=1,\ldots,n$ will be globally optimal and strong duality will hold for ($P$) if $|\alpha_{ij}|\leq\frac{\pi}{n}$ for all $(i,j)\in E$.

Requiring that the angular residuals $|\alpha_{ij}|$ must be less than $\pi/n$ for the global solution may seem like a restriction, but it is actually not. To see this, note that a non-optimal solution to the rotation averaging problem can be obtained by choosing $R_{1}$ such that the first residual $\alpha_{12}$ is zero, and then continuing in the same fashion such that all but the last residual $\alpha_{1n}$ in the cycle is zero. In the worst case, $\alpha_{1n}=\pi$. However, this is (obviously) non-optimal. A better solution is obtained if we distribute the angular residual error evenly so that $\alpha_{ij}=\alpha=\frac{\alpha_{1n}}{n}$ (which is always possible, see Theorem 23 in [10]). In conclusion, the angular residuals $|\alpha_{ij}|$ of the globally optimal solution for a cycle graph is always less than or equal to $\frac{\pi}{n}$, and conversely, if the angular residual is larger than $\frac{\pi}{n}$ for a local minimizer, then it does not correspond to the global solution.

In Fig. 1, we have a real example of an orbital camera motion which is close to a cycle. It may seem hard to determine if the camera motion consists of one or more loops around the object - we give three different local minima for this example. Still, applying formula (24) for this instance gives $\alpha_{\max}=8.89^{\circ}$ which is typically sufficient in practice to ensure that the optimal solution can be obtained by solving a convex program. Before developing an actual algorithm, we shall prove our main result on strong duality.

4.2 Proof of Theorem 4.1

Recall that a sufficient condition for strong duality to hold is that $\Lambda^{*}-\tilde{R}\succeq 0$ (Lemma 3.2). To prove Theorem 4.1 we will show that this is true under the conditions of the theorem.

To simplify the presentation we denote the residual rotations ${\mathcal{E}}_{ij}=R^{*}_{i}\tilde{R}_{ij}R_{j}^{*T}$ and define

[TABLE]

Then $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}=$

[TABLE]

Note that $\sum_{j\neq i}a_{ij}{\mathcal{E}}_{ij}=\frac{1}{2}\sum_{j\neq i}a_{ij}({\mathcal{E}}_{ij}+{\mathcal{E}}^{T}_{ij})$ by symmetry of $\Lambda^{*}$. Since $D_{R^{*}}$ is orthogonal, the matrix $\Lambda^{*}-\tilde{R}$ is positive semidefinite if and only if $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}$ is.

In the noise free case we note that the residual rotations will fulfill ${\mathcal{E}}_{ij}=I$ and therefore

[TABLE]

In the general noise case our strategy will therefore be to bound the eigenvalues of $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}$ by those of $L_{G}$ for which well-known estimates exist. Thus, we will analyze the difference and define the matrix

[TABLE]

The following results characterize the eigenvalues of $\Delta$.

Lemma 4.1**.**

Let $\Delta_{ij}$, $i=1,...,n$, $j=1,...,n$ be the $3\times 3$ sub-blocks of $\Delta$. If $\lambda$ is an eigenvalue of $\Delta$ then

[TABLE]

Proof.

The proof is similar to that of Gerschgorin’s theorem [12]. Let $\Delta x=\lambda x$, with $\|x\|=1$. Then $\lambda x_{i}=\sum_{j}\Delta_{ij}x_{j}.$ Now pick $i$ such that $\|x_{i}\|\geq\|x_{j}\|$ for all $j$. Then

[TABLE]

∎

Lemma 4.2**.**

Denote $\alpha_{\max}$ the largest (absolute) residual angle of all ${\mathcal{E}}_{ij}$ and assume $0\leq\alpha_{\max}\leq\frac{\pi}{2}$. Then

[TABLE]

where $d_{i}$ is the degree of vertex $i$.

Proof.

It is easy to see that by applying a change of coordinates ${\mathcal{E}}_{ij}$ can be written

[TABLE]

and therefore

[TABLE]

This gives

[TABLE]

and since $\Delta_{ii}=\sum_{j\neq i}a_{ij}\left(\frac{1}{2}({\mathcal{E}}_{ij}+{\mathcal{E}}_{ij}^{T})-I\right)$ we get

[TABLE]

Thus $\|\Delta_{ii}\|\leq d_{i}(1-\cos(\alpha_{\max}))=2d_{i}\sin^{2}(\frac{\alpha_{\max}}{2})$. ∎

Lemma 4.3**.**

If $0\leq\alpha_{max}\leq\frac{\pi}{2}$ and $i\neq j$ then

[TABLE]

Proof.

To estimate the off-diagonal blocks $\|\Delta_{ij}\|=a_{ij}\|I-{\mathcal{E}}_{ij}\|$ we note that for a unit vector $v$ we have

[TABLE]

where $\angle(v,{\mathcal{E}}_{ij}v)$ is the angle between $v$ and ${\mathcal{E}}_{ij}v$. Furthermore, we will have equality if $v$ is perpendicular to the rotation axis of ${\mathcal{E}}_{ij}$. Therefore

[TABLE]

∎

Summarizing the results in Lemmas 4.1- 4.3 we get that the eigenvalues $\lambda$ of $\Delta$ fulfill

[TABLE]

where $d_{\max}$ is the maximal vertex degree. Note that the same bound holds for all eigenvalues of $\Delta$, in particular, the one with the largest magnitude $\lambda_{\max}(\Delta)$.

Now returning to our goal of showing that $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}\succeq 0$. Let $N=\begin{bmatrix}I&I&\ldots\end{bmatrix}^{T}$. The columns of $N$ will be in the nullspace of $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}$. Therefore $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}$ is positive semidefinite if $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}+\mu NN^{T}$ is, and hence it is enough to show that

[TABLE]

for sufficiently large $\mu$. The Laplacian $L_{G}$ is positive semidefinite with smallest eigenvalue $\lambda_{1}=0$ and corresponding eigenvector $v=\left(\begin{matrix}1&1&\ldots&1\end{matrix}\right)^{T}$. Furthermore, as $N=v\otimes I_{3}$, it is clear that for sufficiently large $\mu$ we have $\lambda_{1}(L_{G}\otimes I_{3}+\mu NN^{T})=\lambda_{1}(L_{G}+\mu vv^{T})=\lambda_{2}(L_{G})$. Since

[TABLE]

we therefore get

[TABLE]

If the right-hand side is positive, then so is the left-hand side. Using (39) for $\lambda_{\max}(\Delta)$ yields the following result.

Lemma 4.4**.**

The matrix $\Lambda^{*}-\tilde{R}$ is positive semidefinite if

[TABLE]

By completing squares, one obtains the equivalent condition

[TABLE]

which shows Theorem 4.1.

5 Solving the Rotation Averaging Problem

The dual problem ($D$) is a convex semidefinite program, and although it is theoretically sound and provably solvable in polynomial time by interior point methods [4], in practice such problems quickly become intractable as the dimension of the entering variables grow.

In this section we present a first-order method for solving semidefinite programs with constant block diagonals. Our approach solves the dual of ($D$) and consists of two simple matrix operations only, matrix multiplication and square roots of $3\times 3$ symmetric matrices, the latter which can be solved in closed form. Consequently, these two operations permit a simple and efficient implementation without the need for dedicated numerical libraries.

The dual of ($D$) is given by

[TABLE]

Let the matrix $Y$ be partitioned as follows,

[TABLE]

where each block $Y_{ij}\in{\mathbb{R}}^{3\times 3}$ for $i,j=1,\ldots,n$. Since $\Lambda$ is block-diagonal (15) it is clear that the inner maximization is unbounded when $Y_{ii}-I_{3\times 3}\neq 0$ and zero otherwise. We therefore get

[TABLE]

Since $Y\succeq 0$ it is clear that

[TABLE]

for all $\Lambda$ of the form (15). Therefore $(DD)\geq(D)$ and assuming strong duality holds $(D)=(P)$. Furthermore if $R^{*}$ is the global optimum of $(P)$ then $Y={R^{*}}^{T}R^{*}$ is feasible in (47) which shows that $(DD)=(P)$.

Thus, when strong duality holds, recovering a primal solution to $(P)$ is then achieved by simply reading off the first three rows of $Y^{*}$ and choosing their signs to ensure positive determinants of the resulting rotation matrices, see supplementary material for further details.

5.1 Block Coordinate Descent

In this section we present a block coordinate descent method for solving semidefinite programs with block diagonal constraints on the form (47). This method is a generalization of the row-by-row algorithms derived in [27].

Consider the following semidefinite program,

[TABLE]

This is a subproblem that arises when attempting to solve ($DD$) in (47) using a block coordinate descent approach, i.e., by fixing all but one row and column of blocks in (46) and reordering as necessary. It turns out that this subproblem has a particularly simple, closed form solution, established by the following lemma.

Lemma 5.1**.**

Let $B$ be a positive semidefinite matrix. Then, the solution to (48) is given by,

[TABLE]

Here † denotes the Moore–Penrose pseudoinverse.

Proof.

See supplementary material. ∎

6 Experimental Results

In this section we present an experimental study aimed at characterizing the performance and computational efficiency of the proposed algorithm compared to existing standard numerical solvers.

Synthetic data.

In our first set of experiments we compared the computational efficiency of the Levenberg-Marquardt (LM) algorithm [29], a standard nonlinear optimization method, Algorithm 1 and that of SeDuMi [24], a publicly available software package for conic optimization.

We constructed a large number of synthetic problem instances of increasing size, perturbed by varying levels of noise. Each absolute rotation was obtained by rotation about the z-axis by $2\pi/n$ rad and by construction, forming a cycle graph. The relative rotations were perturbed by noise in the form of a random rotation about an axis sampled from a uniform distribution on the unit sphere with angles normally distributed with mean [math] and variance $\sigma$. The absolute rotations were initialized (if required) in a similar fashion but with the angles uniformly distributed over $[0,2\pi]$ rad.

The results, averaged over $50$ runs, can be seen in Table 1. As expected, the LM algorithm significantly outperforms our algorithm as well as SeDuMi, but it only manages to obtain the global optima in about $30-70\%$ of the time. As predicted by Theorem 4.2 and the discussion in Section 4.1 on cycle graphs, both Algorithm 1 and SeDuMi produce globally optimal solutions at every single problem instance, independent of the noise level and independent on the number of cameras. From this table we also observe that Algorithm 1 does appear to outperform SeDuMi quite significantly with respect to computational efficiency.

Real-world data.

In our second set of experiments we compared the computational efficiency on a number of publicly available real-world datasets [11]. The results, again averaged over $50$ runs, are presented in Table 2. Here, as in the previous experiment, both methods correctly produce the global optima at each instance. Algorithm 1 again significantly outperforms SeDuMi in computational cost, providing further evidence of the efficiency of the proposed algorithm. It can further be seen that Theorem 4.1 provides bounds sufficiently large to guarantee strong duality, and hence global optimality, in all the real-world instances except for one, the Pumpkin dataset. Although strong duality does indeed hold in this case, the resulting certificate is less than the largest angular residual obtained. The camera graph is comprised both of densely as well as sparsely connected cameras, resulting in a large value of $d_{\max}$ in combination with a small value of $d_{\min}$ (minimum degree). Since $\lambda_{2}\leq d_{\min}$ a limited bound on $\alpha_{\max}$ follows directly from (24). This instance serves as a representative example of when the bounds of Theorem 4.1, although still valid and strictly positive, become too conservative in practice.

7 Conclusions

In this paper we have presented a theoretical analysis of Lagrangian duality in rotation averaging based on spectral graph theory. Our main result states that for this class of problems strong duality will provably hold between the primal and dual formulations if the noise levels are sufficiently restricted. In many cases the noise levels required for strong duality not to hold can be shown to be quite severe. To the best of our knowledge, this is the first time such practically useful sufficient conditions for strong duality have been established for optimization over multiple rotations.

A scalable first-order algorithm, a generalization of coordinate descent methods for semidefinite cone programming, was also presented. Our empirical validation demonstrates the potential of this proposed algorithm, significantly outperforming existing general purpose numerical solvers.

Proof of Theorem 4.2

Theorem 4.2**.**

Let $R_{i}^{*}$, $i=1,\ldots,n$ denote a stationary point to the primal problem ($P$) for a cycle graph with $n$ vertices. Let $\alpha_{ij}$ denote the angular residuals, i.e., $\alpha_{ij}=\angle(R^{*}_{i}\tilde{R}_{ij},R^{*}_{j})$. Then, $R_{i}^{*}$, $i=1,\ldots,n$ will be globally optimal and strong duality will hold for ($P$) if

[TABLE]

Proof.

A sufficient condition for strong duality to hold is that $\Lambda^{*}-\tilde{R}\succeq 0$ (Lemma 3.2), which is equivalent to $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}\succeq 0$ with the same notation and argument as in (25) and (26). For a cycle graph, we get $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}=$

[TABLE]

As this matrix is symmetric, it implies for the first diagonal block that ${\mathcal{E}}_{12}-{\mathcal{E}}^{T}_{12}={\mathcal{E}}^{T}_{1n}-{\mathcal{E}}_{1n}$. As all ${\mathcal{E}}_{ij}\in{\rm SO}(3)$, it follows that ${\mathcal{E}}_{12}={\mathcal{E}}^{T}_{1n}={\mathcal{E}}$ for some rotation ${\mathcal{E}}\in{\rm SO}(3)$. Similarly, for the second diagonal block ${\mathcal{E}}_{12}={\mathcal{E}}^{T}_{23}={\mathcal{E}}$ and by induction, the matrix $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}$ has the following tridiagonal (Laplacian-like) structure

[TABLE]

Note that this means that the total error is equally distributed in an optimal solution among all the residuals, in particular, $\alpha_{ij}=\alpha$ for all $(i,j)\in E$, where $\alpha$ is the residual rotation angle of ${\mathcal{E}}$.

Let $v$ denote the rotation axis of ${\mathcal{E}}$ and let $u$ and $w$ be an orthogonal base which is orthogonal to $v$. Then, define the two vectors $v_{\pm}=(\,v_{\pm,1}\quad v_{\pm,2}\enspace\ldots\enspace v_{\pm,n}\,)^{T}$, where $v_{\pm,i}=\cos(\frac{2\pi i}{n})u\pm\sin(\frac{2\pi i}{n})w$ for $i=1,\ldots,n$. Now it is straight-forward to check that $v_{\pm}$ are eigenvectors to (51) with eigenvalues $4\sin(\frac{\pi}{n}\pm\alpha)\sin(\frac{\pi}{n})$. The sign of the smallest of these two eigenvalues determines the positive definiteness of the matrix in (51). In other words, we have shown that if $|\alpha|\leq\frac{\pi}{n}$ then $D_{R^{*}}(\Lambda^{*}-\tilde{R})D_{R^{*}}^{T}\succeq 0$. ∎

Proof of Lemma 5.1

Lemma 5.1**.**

Let $B$ be a positive semidefinite matrix. Then, the solution to (48) is given by,

[TABLE]

Proof.

From the Schur complement, we have that the $2\times 2$ block matrix in (48) is positive semidefinite if and only if

[TABLE]

Hence the problem (48) is equivalent to

[TABLE]

The KKT conditions for (55), with Lagrangian multipliers $\Gamma$ and $\Upsilon$ , become

[TABLE]

Rewrite (56) and (60) as

[TABLE]

Since the pseudoinverse fulfills $B^{\dagger}BB^{\dagger}=B^{\dagger}$, combining (61) and (62) we obtain

[TABLE]

Here the last equality follows since $B(I-BB^{\dagger})=0$. This gives

[TABLE]

Inserting (66) in (61)

[TABLE]

multiplying with $B$ form the left on both sides and using (58), $BB^{\dagger}S=S$, we arrive at

[TABLE]

and consequently

[TABLE]

Finally, since

[TABLE]

the conditions (57) and (59) are satisfied then (52) must be a feasible and optimal solution to (55) and consequently also to (48). ∎