Linearly Converging Quasi Branch and Bound Algorithms for Global Rigid Registration

TL;DR

This paper introduces Quasi BnB, a novel framework that uses quadratic quasi-lower bounds to achieve linear convergence in global rigid registration, outperforming traditional BnB algorithms especially at high accuracy levels.

Contribution

The paper proposes Quasi BnB, a new approach that replaces linear bounds with quadratic quasi-lower bounds, enabling linear convergence and improved efficiency in global rigid registration.

Findings

Quasi BnB achieves $ ilde{O}( ext{log}(1/ extepsilon))$ time complexity.

Quasi BnB outperforms state-of-the-art BnB algorithms in efficiency.

The method is especially effective for high-accuracy registration problems.

Abstract

In recent years, several branch-and-bound (BnB) algorithms have been proposed to globally optimize rigid registration problems. In this paper, we suggest a general framework to improve upon the BnB approach, which we name Quasi BnB. Quasi BnB replaces the linear lower bounds used in BnB algorithms with quadratic quasi-lower bounds which are based on the quadratic behavior of the energy in the vicinity of the global minimum. While quasi-lower bounds are not truly lower bounds, the Quasi-BnB algorithm is globally optimal. In fact we prove that it exhibits linear convergence -- it achieves $\epsilon$-accuracy in $~O(\log(1/\epsilon)) $ time while the time complexity of other rigid registration BnB algorithms is polynomial in $1/\epsilon $. Our experiments verify that Quasi-BnB is significantly more efficient than state-of-the-art BnB algorithms, especially for problems where high accuracy is desired.

Figures

Equations

$$F(g)=\underset{\pi }{\min }E(g,\pi )$$

$$\underset{x\in {\mathbb{R}}^D,y\in \mathcal{Y}}{\min }E(x,y),$$

$$\underset{\begin{array}{c}(R,t)\in G_{\mathrm{CP}}\\ \pi \in {\Pi }_{\mathrm{CP}}\end{array}}{\min }{E}_{\mathrm{CP}}(R,t,\pi )=\frac{1}{n}\sum _{i=1}^n\Vert R{p}_i+t-q_{\pi (i)}{\Vert }^2$$

$$R_r=\exp ([r])\in \mathrm{SO}(d).$$

$$\underset{(r,t)\in {\mathbb{R}}^s\times {\mathbb{R}}^d,\pi \in {\Pi }_{\mathrm{CP}}}{\min }{E}_{\mathrm{CP}}(R_r,t,\pi ).$$

$$\underset{r\in {\mathbb{R}}^s,\pi \in {\Pi }_{\mathrm{bi}}}{\min }{E}_{\mathrm{bi}}(R_r,\pi )=\frac{1}{n}\sum _i\Vert R_r{p}_i-q_{\pi (i)}{\Vert }^2$$

$$F(x)=\underset{y\in \mathcal{Y}}{\min }E(x,y).$$

$$F(x_1)-F(x_2)\le \Delta (\delta )=L\delta +O({\delta }^2).$$

$${\mathrm{lb}}_i\equiv F(x_i)-\Delta (\sqrt{D}{h}_i).$$

$$F(x)-F(x_{\ast })\le C\Vert x-x_{\ast }{\Vert }^2,\forall x\in {\mathcal{C}}^0$$

$$F(x)-F(x_{\ast })$$

$$\le C\Vert x-x_{\ast }{\Vert }^2\text{qed}$$

$$F(x)-F(x_{\ast })\le {\Delta }_{\ast }(\delta )=C{\delta }^2+O({\delta }^3).$$

$${\mathrm{lb}}_i\equiv F(x_i)-{\Delta }_{\ast }(\sqrt{D}{h}_i).$$

$$C_1{ϵ}^{-D/2}\le n_{\mathrm{BnB}}\le C_2{ϵ}^{-D}.$$

$$n_{\mathrm{qBnB}}\le C_3{ϵ}^{-D/2}.$$

$$n_{\mathrm{qBnB}}\le C_4{\log }_2(1/ϵ).$$

$${\psi }_k(x)=e^x-\sum _{j=0}^{k-1}\frac{x^j}{j!}$$

$$F_{\mathrm{bi}}(r)-F_{\mathrm{bi}}(r_{\ast })\le {\Delta }_{\ast }(\delta )\equiv \frac{2}{n}{\sigma }_{\mathcal{P}}{\sigma }_Q{\psi }_2(\delta )$$

$$E_{\mathrm{bi}}(R,\pi )=\frac{{\sigma }_{\mathcal{P}}^2+{\sigma }_{\mathcal{Q}}^2-2{\sum }_i⟨R{p}_i,q_{\pi (i)}⟩}{n}$$

$$F_{\mathrm{bi}}(r)$$

$$\stackrel{\text{(}\ast )}{=}-\frac{2}{n}\sum _{i=1}^n\sum _{k=2}^{\infty }⟨\frac{1}{k!}[r{]}^k{p}_i,q_{\pi (i)}⟩$$

$$\le \frac{2}{n}\sum _{i=1}^n\sum _{k=2}^{\infty }\frac{1}{k!}\Vert [r]{\Vert }_{\mathrm{op}}^k\Vert p_i\Vert \Vert q_{{\pi }_{\ast }(i)}\Vert$$

$$\le \frac{2}{n}{\psi }_2(\Vert [r]{\Vert }_{\mathrm{op}}){\sigma }_{\mathcal{P}}{\sigma }_{\mathcal{Q}}.$$

$$F_{\mathrm{CP}}(r,t)-F_{\mathrm{CP}}(r_{\ast },t_{\ast })\le {\Delta }_{\ast }({\delta }_1,{\delta }_2)$$

$${\Delta }_{\ast }({\delta }_1,{\delta }_2)$$

$$\displaystyle+2\delta_{2}\psi_{1}(\delta_{1})\sum_{i}\|p_{i}\|+n\delta_{2}^{2}\bigg{]}$$

$$F(x_1)-F(x_2)\le L\Vert x_1-x_2\Vert ,\forall x_1,x_2\in {\mathcal{C}}^0$$

$$F(x_1)-F(x_2)\le E(x_1,y_2)-E(x_2,y_2)\le L\Vert x_1-x_2\Vert$$

$${\Delta }_{\alpha }(\delta )=M{\delta }^{\alpha },\alpha =1,2.$$

$$n_1=n_{\mathrm{BnB}},n_2=n_{\mathrm{qBnB}}.$$

$$F(x_{\ast })-F^{\ast }=\mathrm{ub}-F^{\ast }\le \mathrm{ub}-\mathrm{lb}\le ϵ.$$

$${\Delta }_{\alpha }(\sqrt{D}h(g))\le ϵ.$$

Full text

1 Introduction

Rigid registration is a fundamental problem in computer vision and related fields. The input to this problem are two similar shapes and the goal is to find the rigid motion that best aligns the two shapes, as well as a good mapping between the aligned shapes. There are several different formulations of the rigid registration problem. In this paper we focus on two popular formulations, the rigid closest point (rigid-CP) [7] problem which is appropriate for partial matching problems (e.g.where a partial point cloud obtained from a scanning device is mapped to a reconstructed model), and the rigid-bijective problem [28] which is applicable to full matching problems where it naturally defines a distance between shapes [4].

Both of these rigid registration formulations optimize an appropriate energy $E(g,\pi)$ that depends on the chosen rigid motion $g$ and a correspondence $\pi$ between the two given point clouds. This energy is typically non-convex, however it is conditionally tractable in that there exist efficient algorithms for globally optimizing for either of its variables when the other is fixed. This is put to good use by the well-known iterative closest point (ICP) algorithm [7] that optimizes $E$ using a very efficient alternating approach. It converges, however, to a local minimum and thus it strongly depends on a good initialization.

Conditional tractability can also be useful for global optimization algorithms. For a fixed rigid motion $g$ it is possible to compute

[TABLE]

and so to optimize $E$ globally it suffices to optimize $F$. Figure 2(a) shows an example function $F$ derived from a 2D rigid registration problem. The function $F$ is non-convex and non-differentiable, and the complexity of global optimization algorithms for optimizing $F$ is typically exponential in the dimension $D$ of the rigid motion space. However, in many cases of interest $D$ is a small constant and so these algorithms are in fact tractable. Thus, while high dimensional rigid registration problems can be computationally hard [13], they are fixed parameter tractable [12].

Branch and bound (BnB) algorithms [11, 33, 25, 21] are perhaps the only algorithms with a deterministic guarantee to globally optimize the rigid registration problem. The key ingredient in BnB algorithms is the ability to give a lower bound for the value of $F$ in a certain region, based on a single evaluation of $F$. For example, a lower bound for the grey interval in Figure 2(c) would attempt to estimate the minimum of $F$ on the interval (denoted by the purple line (d)), based on the blue point in the center of the interval.

In this paper we introduce the term quasi-lower bound, which is a lower bound for the value of the function, over a given region, under the possibly-false assumption that the region contains a global minimum. While quasi-lower bounds are not true bounds, the global optimality of qBnB algorithms is not compromised if lower bounds are replaced by quasi-lower bounds - leading to algorithms that we name quasi-BnB (qBnB).

The advantage of quasi-lower bounds is that for conditionally smooth energies (as defined in Definition 1) $F$ exhibits quadratic behavior near the global minimum, and this can be leveraged to define quadratic quasi-lower bounds, in contrast with lower bounds that are typically linear. This leads to much improved bounds, as illustrated in Figure 2, where the quasi-lower bound (e) is able to provide a bound that is higher than the true minimal value (d) in the grey interval, and hence tighter than any lower bound.

The tightness of quasi-lower bounds yields substantial benefits in terms of complexity –we show that the complexity of achieving an $\epsilon$-optimal solution using BnB is at best $\sim\epsilon^{-D/2}$. In contrast, under some weak assumptions, the complexity of the qBnB algorithm is $\sim\log(1/\epsilon)$, that is, it achieves a linear convergence rate (in the sense of [8]). Our experiments (e.g.Figure 1) verify these theoretical results, and show that qBnB can be considerably more efficient than state-of-the-art BnB algorithms, especially in the presence of high noise or when high accuracy is required.

2 Related Work

We focus on global optimization methods for rigid registration. For other aspects of the rigid registration problem we refer the reader to surveys such as [30].

A variety of methods from the global optimization toolbox have been used to address the rigid registration problem, including Monte Carlo [19], RANSAC based [3], particle filtering [29], particle swarming [32], convex relaxations [22, 20], and graduated non-convexity [35]. Unlike BnB algorithms, none of these algorithms come with a deterministic guarantee for global convergence.

In [25, 23, 10] BnB algorithms for rigid matching in 2D are suggested. In [16] an algorithm for rigid matching in 3D is proposed, based on a combination of matching rotation-invariant descriptors and BnB. In [21] a BnB algorithm for a modified bijective matching objective in 3D is suggested, based on Lipschitz optimization.

In [33, 34] the Go-ICP algorithm for rigid matching in 3D is suggested. Their algorithm searches the 6D transformation space very efficiently in the low noise regime, although it can be slow in the high noise regime. A branch and bound procedure for the rigid matching energy of [10] in 3D is proposed in [11, 24]. For this energy they suggest various speedups that are not available for the classic maximum likelihood objective considered in ICP.

The quasi-BnB approach we propose here can boost the convergence of any of the methods above which consider conditionally smooth objective, such as [33, 25, 21], by simply replacing their first order lower bound with an appropriate second order quasi-lower bound. However it is not directly applicable for optimizing objectives which are not conditionally smooth such as the geometric matching objective of [11, 10].

BnB algorithms with quadratic lower bounds have been proposed for smooth unconstrained optimization [15], and for camera pose estimation [18]. The latter approach obtains quadratic lower bounds by jointly optimizing over a linear approximation of a rotation, and the remaining variables; This approach does not straightforwardly apply to our scenario. In contrast, our qBnB approach only requires solving optimization problems with a fixed rotation.

[26] show that ICP achieves linear converge rate, and suggest alternative second order methods with super-linear convergence rate. These algorithms only converge to a local minimum, whereas our algorithm is shown to exhibit linear convergence to the global minimum.

3 Method

3.1 Problem statement

We consider the application of the quasi BnB framework as described in Section 1 to the problem of rigid alignment of point clouds. This problem has several variants, and in this paper we will focus on two of them that we will refer to as the rigid-closest-point (rigid-CP) problem and the rigid-bijective problem. We also wish to focus on the common structure of problems for which our quasi BnB framework is applicable, which we will name in this paper $D$-quasi-optimizable. We begin by defining $D$-quasi-optimizable problems.

Consider optimization problems of the form

[TABLE]

where $\mathcal{Y}$ is a finite set. Denote the cube centered at $x$ with half-edge length $h$ by $\mathcal{C}_{h}(x)$.

Definition 1**.**

We say that an optimization problem of the form (1) is $D$-quasi-optimizable in a cube $\mathbf{\mathcal{C}}^{0}=\mathcal{C}_{h_{0}}(x_{0})$ if it satisfies the following conditions:

Existence of a minimizer:* There exists a minimizer $(x_{*},y_{*})$ of $E$ such that $x_{*}\in\mathbf{\mathcal{C}}^{0}$.* 2. 2.

$D$-tractability*: For fixed $x$, minimizing $E(x,\cdot)$ over $\mathcal{Y}$ can be performed in polynomial time.* 3. 3.

Conditional smoothness* For fixed $y\in\mathcal{Y}$, the function $E(\cdot,y)$ is smooth.*

We now turn to defining the rigid-CP and rigid-bijective problems and explaining why they are $D$-quasi-optimizable.

The input to the rigid-CP problem are two point clouds $\mathcal{P}=\{p_{1},\ldots,p_{n}\}$ and $\mathcal{Q}=\{q_{1},\ldots,q_{m}\}$ in $\mathbb{R}^{d}$ where we assume $d=3$ or $d=2$. We assume the point clouds are normalized to be in $[-1,1]^{d}$ and have zero mean. Our goal is to find a rigid motion that aligns the points as-well-as-possible. Namely, denoting the group of (orientation preserving) rigid motions $\mathrm{SO}(d)\times\mathbb{R}^{d}$ by $G_{\mathrm{CP}}$, and denoting by $\Pi_{\mathrm{CP}}$ the collection of all functions $\pi:\mathcal{P}\to\mathcal{Q}$, our goal is to solve the minimization problem

[TABLE]

Following [33, 21] we simplify the domain of (2) by using the exponential map: We use $s$ to denote the intrinsic dimension of $\mathrm{SO}(d)$. For a vector $r\in\mathbb{R}^{s}$, we define $[r]$ to be the unique $d\times d$ skew-symmetric matrix ($[r]^{T}$=$-[r]$) whose entries under the diagonal are given by $r$. By applying the matrix exponential to $[r]$ we get a matrix

[TABLE]

Every rotation can be represented as $R_{r}$ for some $r$ in the closed ball $B_{\pi}(0)$ and so (2) can be reduced to the problem

[TABLE]

Our next step is to identify a cube $\mathbf{\mathcal{C}}^{0}$ in $\mathbb{R}^{s}\times\mathbb{R}^{d}$ in which $E_{\mathrm{CP}}$ must have a global minimum. For the rotation component we can take the cube $\mathcal{C}_{\pi}(0)$ that bounds the ball $B_{\pi}(0)$ . For the translation component, we note that if $(R_{*},\pi_{*},t_{*})$ is a minimizer of $E$, then the optimal translation is the difference between the average of $R_{*}p_{i}$, which is zero by assumption, and the average of $q_{\pi(i)}$, which will be in the unit cube. Thus there exists a minimizer $(r_{*},t_{*},\pi_{*})$ such that $(r_{*},t_{*})$ is in $\mathcal{C}_{\pi}(0)\times\mathcal{C}_{1}(0)$. This shows that (3) satisfies the first condition for being $D$-quasi-optimizable, where ${x=(r,t)}$, $y=\pi$ and $D=s+d$ (i.e., $D=6$ or $3$ for 3D and 2D problems respectively).

$D$-tractability follows from the fact that the optimal mappings $\pi$ minimizing $E_{\mathrm{CP}}(R,t,\cdot)$ for fixed $R,t$ are just the mappings that take each rigidly transformed point $Rp_{i}+t$ to its closest point in $\mathcal{Q}$, hence the term rigid-CP. Conditional smoothness is obvious.

The rigid bijective problem is similar to rigid-CP, but focuses on the case $n=m$, and only allows mappings between $\mathcal{P}$ and $\mathcal{Q}$ that are bijective (i.e., permutations). We denote this set of mappings by $\Pi_{\mathrm{bi}}$. In this scenario optimizing for $\pi$ while holding the rigid transformation component fixed is tractable, but more computationally intensive than rigid-CP since it requires solving a linear assignment problem. On the other hand, in this scenario the optimal translation is always $t_{*}=0$, since the mean of the points $p_{i}$ and $q_{\pi(i)}$ is zero. Therefore we can reduce our problem to lower-dimensional optimization over $G_{\mathrm{bi}}=\mathrm{SO}(d)$. As before, we use the exponential map to reparameterize the problem as

[TABLE]

It follows from our discussion above that this optimization problem is $D$-quasi-optimizable over $\mathcal{C}_{\pi}(0)$ with $x=r,y=\pi$ and $D=s$ (i.e., $D=3$ or $1$ for 3D and 2D problems respectively).

3.2 Optimizing quasi-optimizable functions

BnB algorithms.

$D$-tractability implies that we can reduce (1) to the equivalent problem of minimizing the $D$-dimensional function

[TABLE]

BnB algorithms for rigid registration and related problems are typically based on the ability to show that for $\delta>0$ and $x_{1},x_{2}\in\mathbf{\mathcal{C}}^{0}$ satisfying $\|x_{1}-x_{2}\|<\delta$,

[TABLE]

Using a bound of this form $F$ can be bounded from below in a cube $\mathcal{C}_{h_{i}}(x_{i})$ by evaluating $F(x_{i})$ and noting that for any $x$ in the cube, (6) implies that $F(x)$ is larger than

[TABLE]

In the appendix we show that any $F$ arising from a $D$-quasi-optimizable optimization problem is Lipschitz and thus a bound of the form (6) always exists.

Based on the lower bound (7), a simple breadth-first-search (BFS) BnB algorithm starts from a coarse partitioning of $\mathbf{\mathcal{C}}^{0}$ into sub-cubes $\mathcal{C}_{h_{i}}(x_{i})$, and then evaluates $F$ at each of the $x_{i}$. Each such evaluation gives a global upper bound $\mathrm{ub}_{i}=F(x_{i})$ for the minimum $F^{*}$, and a local lower bound $\mathrm{lb}_{i}$ for the minimum on the cube as defined in (7). At each step the algorithm keeps track of the best global upper bound found so far, which we denote by $\mathrm{ub}$, and on the global lower bound $\mathrm{lb}$ that is defined as the smallest local lower bound found in this partition. Now for every $i$, If $\mathrm{lb}_{i}>\mathrm{ub}$ then $F$ is not minimized in $\mathcal{C}_{h_{i}}(x_{i})$ and this cube can be excluded from the search. The BnB algorithm then refines the partition into smaller cubes for cubes that have not yet been eliminated, and this process is continued until $\mathrm{ub}-\mathrm{lb}<\epsilon$. This is the BnB strategy used in [33, 25, 21, 23, 10] for rigid registration problems, though these papers vary in the search strategy they employ, the rigid energy the consider, the bound (6) they compute, and other aspects.

Quasi-BnB.

The qBnB algorithm we suggest in this paper is based on the observation that the non-differentiable $F$ resembles smooth functions in that its behavior near minimizers is quadratic:

Lemma 1**.**

If $x_{*}\in\mathbf{\mathcal{C}}^{0}$ is a global minimizer of $F$, then there exists $C>0$ such that

[TABLE]

Proof.

Assume $(x_{*},y_{*})$ is a global minimizer of $E$, Let $N_{x,y}$ be the operator norm of the Hessian of the function $E(\cdot,y)$ at a point $x$, and let $C$ be the maximum of $N_{x,y}$ over $\mathbf{\mathcal{C}}^{0}\times\mathcal{Y}$. Since the first order approximation of $E(\cdot,y_{*})$ at $x_{*}$ vanishes, we obtain

[TABLE]

Let us assume that we have an explicit quadratic bound of the general form (9). That is, if $x_{*}$ is a global minimizer, and $x,x_{*}\in\mathbf{\mathcal{C}}^{0}$ satisfy $\|x-x_{*}\|<\delta$, then

[TABLE]

In Subsection 3.4 we give an explicit computation of $\Delta_{*}$ for rigid-bijective and rigid-CP problems. While $\Delta_{*}$ only induces a bound on $F$ near a global minimum, it can be used similarly to the Lipschitz bounds to eliminate cubes $\mathcal{C}_{h_{i}}(x_{i})$, for if the global minimizer $x_{*}$ is contained in $\mathcal{C}_{h_{i}}(x_{i})$ then $F(x_{*})$ cannot be larger than

[TABLE]

Thus, $\mathcal{C}_{h_{i}}(x_{i})$ can be eliminated whenever ${\mathrm{lb}_{i}>\mathrm{ub}}$. We emphasize that $\mathrm{lb}_{i}$ is not a lower bound for the value of $F$ in the cube since if the cube does not contain a global minimizer $F$ may vary linearly inside the cube – for this reason we name it a quasi-lower bound. Nonetheless from the point of view of BnB algorithms (see Figure 2), linear lower bounds can be replaced by quadratic quasi-lower bounds, leading to a significantly more efficient qBnB algorithm, still guaranteed to converge to an $\epsilon$-optimal solution. A simple BFS qBnB algorithm for optimizing $D$-quasi-optimizable functions, for a given quasi-lower bound $\Delta_{*}(\delta)$, is provided in Algorithm 1.

3.3 Complexity

Theorem 2 (below, proof in appendix) provides complexity bounds for the qBnB algorithm with a quadratic $\Delta_{*}(\delta)$, as in Algorithm 1, and for a similar BnB algorithm, where $\Delta_{*}$ is replaced by a Lipschitz bound $\Delta(\delta)$. We denote the number of $F$-evaluations in the algorithms by $n_{\mathrm{qBnB}}$ and $n_{\mathrm{BnB}}$ respectively. We consider the limit where the prescribed accuracy $\epsilon$ tends to zero and all other parameters of the problem are held fixed:

Theorem 2**.**

There exist positive constants $C_{1},\ldots,C_{4}$, such that

[TABLE]

[TABLE]

Furthermore, if $E$ has a finite number of minimizers $(x_{\ell}^{*},y_{\ell}^{*})_{\ell=1}^{N}$, and the Hessian of $E(\cdot,y_{\ell}^{*})$ is strictly positive definite for all $\ell$, then

[TABLE]

3.4 Computing quasi-lower bounds

To conclude our discussion we now need to state explicit quadratic quasi-lower bounds for the rigid-CP and rigid-bijective problems. We denote by $F_{\mathrm{CP}}$ and $F_{\mathrm{bi}}$ the functions obtained by applying (5) to $E_{\mathrm{CP}}$ and $E_{\mathrm{bi}}$.

Rigid bijective quasi-lower bounds.

For $k\in\mathbb{N}$, let $\psi_{k}$ denote the truncated Taylor expansion of $e^{x}$,

[TABLE]

Theorem 3**.**

Let $\delta>0,r\in\mathbb{R}^{D}$ and $r_{*}$ be a global minimizer of $F_{\mathrm{bi}}$, and assume $\|r-r_{*}\|\leq\delta$. Let $\sigma_{\mathcal{P}},\sigma_{\mathcal{Q}}$ denote the Frobenius norm of the matrices whose columns are the points in $\mathcal{P}$ and $\mathcal{Q}$ respectively. Then $\Delta_{*}(\delta)$ is given by

[TABLE]

Proof.

Let us consider the case where $r_{*}=0$, so $(I_{d},\pi_{*})$ minimizes $E_{\mathrm{bi}}$ for some appropriate $\pi_{*}$. Note that for all $R,\pi$,

[TABLE]

it follows that

[TABLE]

To obtain the equality $(*)$ note that $r_{*}=0$ minimizes $E_{\mathrm{bi}}(\cdot,\pi_{*})$ and so the first order approximation of this function vanishes. In the appendix we show that $\|[r]\|_{\mathrm{op}}\leq\|r\|$ which concludes the proof for the case $r_{*}=0$, and then use this case to prove the theorem for general $r_{*}$. ∎

Based on this quasi-lower bound, our algorithm for optimizing the rigid-bijective problem is just Algorithm 1, where as stated in Subsection 3.1 we take $x=r,y=\pi,\mathcal{C}_{h_{0}}(x_{0})=\mathcal{C}_{\pi}(0)$, and we use $\Delta_{*}$ defined in (15).

Rigid CP quasi-lower bounds.

For the rigid-CP problem we obtain the following quasi-lower bound that we prove in the appendix.

Theorem 4**.**

Let $(r_{*},t_{*})$ be a minimizer of $F_{\mathrm{CP}}$, and let $(r,t)\in\mathbb{R}^{s}\times\mathbb{R}^{d}$, and $\delta_{1},\delta_{2}>0$ which satisfy $\|r-r_{*}\|\leq\delta_{1}$and$\|t-t_{*}\|\leq\delta_{2}$. Let $f_{*}$ be some upper bound for the global minimum of $F_{\mathrm{CP}}$. Then

[TABLE]

where

[TABLE]

Using this quasi-lower bound the most straightforward way to construct a qBnB algorithm is by using Algorithm 1 setting $x=(r,t),y=\pi$ and using the initial cube $\mathcal{C}_{\pi}(0)\times\mathcal{C}_{1}(0)$ as described in Subsection 3.1 and the bound $\Delta_{*}(\delta_{1},\delta_{2})$ from (4). However, to enable simple comparison with the Go-ICP algorithm [33] we use their BnB architecture, wherein two nested BnB are used – an outer BnB for the rotation component and a separate inner BnB for the translation component. The quadratic quasi-lower bounds for these BnBs can be computed by setting $\delta_{2}=0$ or $\delta_{1}=0$, respectively. Further details are provided in the appendix.

4 Results

Evaluation of qBnB for rigid-CP.

We evaluate the performance of the proposed qBnB algorithm for the rigid-CP problem. We compare with Go-ICP [33], a state-of-the-art approach for global minimization of this problem. Our implementation is based on a modified version of the Go-ICP code, which utilizes a nested BnB (see section 3.4) and an efficient approximate CP computation (see implementation details below).

We ran both algorithms on synthetic rigid-CP problems that were generated by uniformly sampling $n$ points on the unit cube to form $\mathcal{P}$, and then applying a random rigid transformation and Gaussian noise with std $\sigma$ to form $\mathcal{Q}$. Finally $\mathcal{P}$ and $\mathcal{Q}$ are translated and scaled to have zero mean and reside in the unit cube.

Figure 3 shows comparisons performed with varying (a) prescribed accuracy and (b) noise. Both algorithms are comparable in the low noise/low accuracy regime. This is consistent with [33] that report better performance in the regime where the optimal energy $E^{*}$ is lower than the error tolerance $\epsilon$. In contrast, when accuracy or noise are increased, the complexity of the Go-ICP algorithm grows rapidly, while the complexity of qBnB stabilizes at an almost constant number of $\small{\approx}10^{6}$ function evaluations. This results is consistent with the complexity analysis of Theorem 2.

In Figure 3(a) we took $n=50$, $\sigma=0.05$, and prescribed accuracy $\epsilon$ varying between $10^{-1}$ and $10^{-6}$, and in Figure 3(b) we took $n=100$, $\epsilon=10^{-3}$ and $\sigma$ varying between [math] and $0.5$. Each data point in the graph represents an average over 100 instances with the same parameters. We exclude results for Go-ICP with $\epsilon<10^{-4}$ or $\sigma>0.2$ due to excessive run-time (over six hours). The maximal run-time of our algorithm over all experiments presented was $9$ minutes.

Figure 4 exemplifies the applicability of our algorithm for solving larger rigid-CP problems. We used our algorithm and Go-ICP to register $500$ points sampled from the Stanford bunny [31] to a full 3D model consisting of nearly 36K vertices. Both methods are global, thus they obtained the same results up to small indiscernible differences; therefore, we only show our alignments, as well as the evaluation counts of both algorithms. As in Figure 3(b), qBnB is more efficient than Go-ICP, especially in the high noise regime. For noise level of $\sigma=0.1$ our algorithm required under $7$ minutes, whereas Go-ICP took over $6$ hours to complete.

Evaluation of qBnB for rigid-bijective.

We compare our qBnB algorithm for the rigid bijective problem with the BFS version of the BnB algorithm proposed in [21]. This algorithm modifies the standard $\ell^{2}$ energy we study here, and globally optimizes it over the same feasible set $G_{\mathrm{bi}}\times\Pi_{\mathrm{bi}}$. in Figure 1 we examined the performance of both algorithms on the problem of registering 2D frog silhouettes [2] and 3D chair meshes [17], sampled at $50$ points, with accuracy set at $\epsilon=10^{-6}$. Both algorithms returned identical results, which are visualized on the top of Figure 1. However qBnB took $0.5$ and $6.5$ seconds to solve the 2D and 3D problems respectively, while [18] required $11$ seconds to solve the 2D problem, and did not converge to the required accuracy for the 3D problem after $10^{7}$ evaluations (over $1$ hour). The graphs in Figure 1 show the number of evaluations the algorithms performs in each depth level of the BFS. In both cases qBnB required significantly less evaluations than [21], and furthermore we note that asymptotically qBnB requires a constant number of evaluations at each depth, while [21] requires $\sim 2^{Dg/2}$ iterations at depth $g$, as shown by the log-scale aymptotes in the figure. This illustrates the theoretical results of Theorem 2: As for BnB algorithms $\epsilon\sim 2^{-g}$, they require $\sim\epsilon^{-D/2}$ evaluations. In contrast, the dependence of qBnB on $g\sim\log(1/\epsilon)$ is linear.

Application to biological morphology.

Biological morphology is concerned with the geometry of anatomical shapes. In this context, Auto3dgm [9, 27] is a popular algorithm that seeks for a plausible, consistent alignment of a collection of scanned morphological shapes. It is based on solving the rigid bijective matching problem between pairs of shapes, and on a global step that synchronizes the rotations. The rigid bijective matching problem is solved using an alternating ICP-like algorithm, initialized from one of the $2^{3}$ rotations that align the principal axes of the two shapes. In this application, reflections are typically allowed, and our algorithm can be easily extended to this case as $\mathrm{O}(d)$ is just a disjoint union of two copies of $\mathrm{SO}(d)$.

Figure 5 compares the performance of qBnB to the pairwise matching algorithm of Auto3dgm for aligning a challenging collection of ten almost-symmetric teeth of spider monkeys (Ateles) from [1]. Details on the data for this experiment can be found in Appendix C. Each bone surface is sampled at $200$ points using farthest point sampling. The pairwise energy between all pairs, and the alignment computed between the first tooth and all remaining teeth is shown for qBnB (top) and Auto3dgm (bottom). We see that qBnB obtains lower energy solutions and more plausible alignments which have been verified as semantically correct by biological experts. This comes at the cost of higher complexity ($2-5$ minutes as opposed to $\sim 3$ seconds).

Implementation details.

The qBnB algorithm for rigid-CP was implemented in C++, based on the implementation of [33]. The algorithm for rigid-bijective matching was implemented in Matlab; We solved for the bijective mapping used the excellent MEX implementation of [5] for the auction algorithm [6]. Timings were measured using a Intel 3.10GHz CPU.

For the optimization of large scale rigid-CP with moderate accuracy $\epsilon=10^{-3}$ as in Figure 3(b) and Figure 4, we followed the implementation of GoICP [33] and used the 3D Euclidean distance transform (DT) from [14] with a $300\times 300\times 300$ grid; this provides a very fast approximation of the closest point computation, compared with the standard (accurate) KD-tree based computation. For evaluating the complexity as a function of accuracy in Figure 3(a), we used both GoICP and qBnB without DT transform.

5 Future work and conclusions

We presented the qBnB framework for globally optimizing $D$-quasi-optimizable functions, and demonstrated theoretically and empirically the advantage of this framework over existing BnB algorithm. Future challenges include applying the qBnB framework to other rigid registration problems that can handle outliers better than the standard $\ell^{2}$ energy, as well as using this framework for global optimization of $D$-quasi-optimizable functions in other knowledge domains.

Acknowledgements

The authors would like to thank Prof. Ingrid Daubechies for helpful discussions. Research is partially supported by Simons Math+X Investigators Award 400837.

Appendix A Proofs

We begin by proving that $D$-optimizable functions are Lipschitz on the cube $\mathbf{\mathcal{C}}^{0}$ as we claimed in the beginning of Subsection 3.2. In the remainder of the Section we restate and prove the theorems and lemmas stated in the main text.

Lemma 5**.**

If $F$ is $D$-quasi optimizable in $\mathbf{\mathcal{C}}^{0}$, then there exists $L>0$ such that.

[TABLE]

Proof.

If $F(x_{1})\leq F(x_{2})$ then (A.1) is obvious. Now assume $F(x_{1})>F(x_{2})$ and let $y_{2}$ be the minimizer of $E(x_{2},\cdot)$. Then due to the differentiability of $E(\cdot,y_{2})$,

[TABLE]

where $L$ is the maximum of the norm of the gradient of $E(\cdot,y)$ over all $x\in\mathbf{\mathcal{C}}^{0}$ and $y\in\mathcal{Y}$. ∎

See 2

Proof.

Part 1. We begin with a general discussion of the complexity of Algorithm 1 in both the BnB and quasi-BnB version, and prove the upper bound on $n_{\mathrm{BnB}}$ and $n_{\mathrm{qBnB}}$ simultaneously. To do so we denote

[TABLE]

We denote by $n_{\alpha}$ the number of $F$-evaluations used by the algorithm in each case, so

[TABLE]

We also note that while our bounds (15),(4) include besides $M\delta^{2}$ also higher terms in $\delta$ these terms can be absorbed into $M\delta^{2}$ with a larger value of $M$.

In the proof we will call $g$ form Algorithm 1 the generation of the algorithm. We begin with showing the algorithm terminates and bounding the final value of $g$ which we denote by $g_{f}$. First note if $\mathcal{C}_{h_{i}}(x_{i})\in L_{g}$ contains a global minimum $x_{*}$, then its sub-cubes will always be added to $L_{g+1}$. This is because since $\mathcal{C}_{h_{i}}(x_{i})$ contains a global minimum, $\mathrm{lb}_{i}$ is in fact a true lower bound for the minimal value of $F$ on the cube, which in this case is the global minimum $F_{*}$. Thus the value of $\mathrm{lb}_{i}\leq F_{*}\leq\mathrm{ub}$ and so the cube will be preserved. It also follows that the global lower bound $\mathrm{lb}$ is in fact a true lower bound, since by definition it is smaller than the quasi-lower bound in a cube containing a global minimizer. It follows that if the algorithm does terminate, then the output $x_{*}$ is an $\epsilon$-optimal solution as its difference from the minimum $F^{*}$ satisfies

[TABLE]

Now to bound $g_{f}$ note that cubes in $L_{g}$ have half-edge length $h(g)\equiv h_{0}2^{-g}$. the algorithm must terminate once it visits all cubes of generation $g$ whose edge length $h(g)$ satisfies

[TABLE]

This is because for all cubes $\mathcal{C}_{h_{j}}(x_{j})$ in this generation,

[TABLE]

and by taking the minimum over $j$ we obtain that $\mathrm{ub}-\mathrm{lb}\leq\epsilon$.

Some algebraic manipulation shows that (A.4) occurs when

[TABLE]

where we use $\log=\log_{2}$ throughout this proof. The number of $F$ evaluations $n_{\alpha}$ is bounded by the worst case scenario where all cubes need to be divided in all generations

[TABLE]

this proves the upper bound on $n_{1},n_{2}$.

Part 2. We now show a lower bound on $n_{\mathrm{BnB}}$. Our first step is to show that the sub-cubes of a given cube have better=larger lower bounds. To see this let $\mathcal{C}^{1}=\mathcal{C}_{h_{1}}(x_{1})$ be a cube, and let $\mathcal{C}^{2}=\mathcal{C}_{h_{2}}(x_{2})$ be one of its sub-cubes. Then $h_{2}=h_{1}/2$ and $|x_{1}-x_{2}|=\sqrt{D}h_{1}/2$. It follows that, using the notation of (A.2) with $\alpha=1$,

[TABLE]

Now denoting by $\mathrm{lb}_{1}$ and $\mathrm{lb}_{2}$ the lower bounds computed for the cubes $\mathcal{C}^{1}$ and $\mathcal{C}^{2}$ respectively, we have

[TABLE]

and so we have $\mathrm{lb}_{2}\ \geq\mathrm{lb}_{1}$ as we stated.

Next note that by the quadratic bound (8) we have that for a global minimum $x_{*}$ and $\eta=\sqrt{\epsilon/C}$,

[TABLE]

Now, let $g_{F}$ denote the value of $g_{f}(2\epsilon)$ from (A.5), for the case $\alpha=1$. Recall that $g_{f}(2\epsilon)$ is defined as the first integer for which (A.4) holds, where $\epsilon$ is replaced by $2\epsilon$, and $\alpha=1$. Thus for $g<g_{F}$ we have that

[TABLE]

Let $\mathcal{C}_{h_{i}}(x_{i})$ be a cube of generation $g_{0}$ containing $x_{*}$, where $g_{0}$ is large enough so that the diameter of the cube is smaller than $\eta$ and thus it is contained in $B_{\eta}(x_{*})$. This occurs for

[TABLE]

Every sub-cube $\mathcal{C}_{h_{j}}(x_{j})$ of $\mathcal{C}_{h_{i}}(x_{i})$, from any generation $g_{0}\leq g<g_{F}$, satisfies

[TABLE]

In particular it follows that the cube $\mathcal{C}_{h_{i}}(x_{i})$, and all its sub-cubes, will be visited during the BnB search. This is because we saw that lower bounds improve by refinement, and so by (A.9) any cube from the earlier generations $g<g_{0}$ which contains $\mathcal{C}_{h_{i}}(x_{i})$, also has lower bounds which are lower that the global minimum (by at least $\epsilon$) and so such a cube would not be removed from the search..

We can now bound $n_{\mathrm{BnB}}$ by the number of subcubes of $\mathcal{C}_{h_{i}}(x_{i})$ at the $g_{F}-1$ generation alone:

[TABLE]

Now

[TABLE]

So returning to (A.10) we obtain

[TABLE]

Part 3. We now turn to prove the last part of the theorem. Let $\mathcal{J}(\ell)$ denote the set of indices $k$ for which $(x_{\ell}^{*},y_{k}^{*})$ is a minimizer. Note the we always have that $\ell$ is in $\mathcal{J}_{\ell}$. Let $m$ be half of the minimun over the minimal eigenvalue of the hessian of $E(\cdot,y_{\ell}^{*})$ at $x_{\ell}^{*}$ for all $\ell$. The assumption that $E$ has a finite number $N$ of minimizers $x_{\ell}^{*},y_{\ell}^{*}$, with strictly positive definite hessian, implies that $m>0$, and so for small enough positive $\eta$,

[TABLE]

The minimum of $F$ on $\mathbf{\mathcal{C}}^{0}\setminus\bigcup_{i}B_{\eta/2}(x_{\ell}^{*})$ is strictly larger than $F^{*}$. Therefore there exists some $g_{0}$ independent of $\epsilon$, such that all cubes of generation $g_{0}$ which are not contained in one of the balls $B_{\eta}(x_{\ell}^{*})$ will be removed in the $g_{0}$-th stage.

We now claim that for $g\geq g_{0}$, $g$-th generation cubes $\mathcal{C}_{h_{i}}(x_{i})$ contained in one of the balls $B_{\eta}(x_{\ell}^{*})$ will be removed if

[TABLE]

This is because

[TABLE]

where $(*)$ follows from the fact that if $\mathcal{C}_{h_{i}}(x_{i})$ is the $g$-th generation cube containing $x_{\ell}^{*}$, then

[TABLE]

Now for $g\geq g_{0}$, the condition (A.12) is not fulfilled in at most $\bar{C}=(\sqrt{\frac{2MD}{m}}+2)^{D}$ cubes surrounding each minimizer, and so in total only $N\bar{C}$ cubes can survive each generation $g>g_{0}$. The important point is that this number is independent of $\epsilon$. So the total number of iterations is bounded by the sum of the total number of cubes in all generations $g\leq g_{0}$, which is some constant indepentent of $\epsilon$ which we denote by $b$, and the constant $N\bar{C}$ multiplied by the remaining number of iterations $g_{f}-g_{0}$, that is

[TABLE]

This bound can be replaced with a bound of the form (13) with an appropriate constant. ∎

See 3

Proof.

To conclude the proof of the theorem for the case $r_{*}=0$ we need to show

Lemma 6**.**

For all $r\in\mathbb{R}^{D}$,

[TABLE]

Proof.

The non-zero eigenvalues of a skew-symmetric real matric $[r]$ can be written as

[TABLE]

where $a_{1}\geq a_{2}\ldots>0$. Therefore

[TABLE]

∎

For the general case $r_{*}\neq 0$, we define a change of variable $\tilde{p}_{i}=R_{r_{*}}p_{i}$ and denote by $\tilde{E}_{\mathrm{bi}}$ the energy resulting by replacing $p_{i}$ with $\tilde{p}_{i}$ in the definition of $E_{\mathrm{bi}}$. Then for all $R_{0},\pi$ we have

[TABLE]

In particular $\tilde{r}_{*}=0$ is a minimizer of $\tilde{F}_{\mathrm{bi}}$ which is defined by replacing $E_{\mathrm{bi}}$ with $\tilde{E}_{\mathrm{bi}}$ in the definition of $F_{\mathrm{bi}}$. We claim that there exists $r_{1}$ such that

[TABLE]

In the case $d=2$ we can identify $R_{r}R_{r_{*}}^{T}$ with $e^{i(r-r_{*})}$ and so we can simply choose $r_{1}=r-r_{*}$. For $d=3$ it is proven in Lemma 3.2 in [18] that the angular distance between $R_{r}$ and $R_{r_{*}}$ is smaller or equal to $\|r-r_{*}\|$. This means that the angular distance between $R_{r}R_{r_{*}}^{T}$ and the identity is less than $\|r-r_{*}\|$ and thus implies the existence of $r_{1}$ satisfying (A.14).

Now using the bound from 15 for $\tilde{F}$ which is minimized at zero we obtain

[TABLE]

∎

See 4

Proof.

The proof is very similar to the proof of Theorem 3. Let us first consider the case $(r_{*},t_{*})=(0,0)$, and let $\pi_{*}$ be the corresponding mapping so that $(I_{d},0,\pi_{*})$ minimizes $E_{\mathrm{CP}}$. Then

[TABLE]

Here $(*)$ follows from the fact that $E_{\mathrm{CP}}(\cdot,\cdot,\pi_{*})$ is minimized at the origin and so the first order terms cancel out, and the next inequalities follow from the Cauchy-Schwarz inequality and from Lemma 6.

For general $(r_{*},t_{*})$, we use a change of variables $\tilde{p}_{i}=R_{*}p_{i},\,\tilde{q}_{i}=q_{i}-t_{*}$, and denote by $\tilde{E}_{\mathrm{CP}}$ and $\tilde{F}_{\mathrm{CP}}$ the functions obtained by replacing $p_{i},q_{i}$ by $\tilde{p}_{i},\tilde{q}_{i}$ in the definition of these functions. For given $(r,t,\pi)$ we have

[TABLE]

Let $r_{1}\in\mathbb{R}^{D}$ satisfying (A.14). By applying the theorem to $\tilde{F}_{\mathrm{CP}}$ which is minimized at $(0,0)$ we obtain

[TABLE]

∎

Appendix B BnB for rigid closest point

In the following we explain how we construct a quasi-BnB framework for the rigid CP problem based on the BnB architecture proposed in Go-ICP [33]. Following Go-ICP, we use a nested BnB structure: We perform an “outer” BnB search on the rotation space, wherein the upper and lower bounds are functions of the translation component $t$; in turn, to compute these bounds we perform an “inner” BnB search over the variable $t$. Namely, for given $r_{i}$, we define an upper bound for the outer BnB by

[TABLE]

To compute a quasi-lower bound for the outer BnB we note that if $(r_{*},t_{*})$ minimizes $F_{\mathrm{CP}}$ and $r_{*}\in\mathcal{C}_{h}(r_{i})$, then by using (4) where we set $r=r_{i}$ , take $\delta_{1}$ to be the maximal distance of a point in the cube from the center, $t=t_{*}$ and $\delta_{2}=0$ we obtain

[TABLE]

and since $F_{\mathrm{CP}}(r_{i},t_{*})\geq\bar{E}_{\mathrm{CP}}(r_{i})$ it follows that if $r_{*}\in\mathcal{C}_{h}(r_{i})$ then

[TABLE]

The RHS of the equation above gives us our quasi lower bound for the rotation quasi BnB. To compute $\bar{E}_{\mathrm{CP}}(r_{i})$ we compute a BnB in translation space, where throughout the translation BnB the rotation coordinate $r_{i}$ is fixed. For a given translation cube $\mathcal{C}_{h}(t_{j})$ an upper bound for the value of $\bar{E}_{\mathrm{CP}}(r_{i})$ is given by evaluation of $F_{\mathrm{CP}}(r_{i},t_{j})$. If $t_{*}$ is a minimizer of $F_{\mathrm{CP}}(r_{i},\cdot)$ then a quasi-lower bound in the cube is given by

[TABLE]

We note that this bound is similar to what we would get by setting $\delta_{1}=0$ and $\delta_{2}$ to be the maximal distance in the cube from $t_{j}$ in (4). Although the bound does not follow directly from this equation the derivation is similar, and can be obtained by studying the behavior of a minimizer of $E(r_{i},\cdot,\pi_{*})$, so we do not go into the details. Finally we note that when the quasi-lower bounds in the outer or inner BnB is lower than zero we replace it with zero.

The rest of the architecture of the BnB is also borrowed from Go-ICP. We use best-first-search, where the cube with the lowest lower bound is visited first. Every time the upper bound is improved, an ICP algorithm is run to improve the resolution of the solution. For more details see [33].

Appendix C Morphological data

The morphological data for the experiment shown in Figure 5 comes from the MorphoSource dataset. The figure shows ten different second mandibular molars of spider monkeys (Ateles), which come from three different taxonomical groups. More details on the data are shown in the table below.