Convex Global 3D Registration with Lagrangian Duality

Convex Global 3D Registration with Lagrangian Duality Jesus Briales Javier Gonzalez-Jimenez MAPIR-UMA Group University of Malaga, Spain {jesusbriales,javiergonzalez}@uma.es Abstract The registration of 3D models by a Euclidean transformation is a fundamental task at the core of many application in computer vision. This problem is non-convex due to the presence of rotational constraints, making traditional local optimization methods prone to getting stuck in local minima. This paper addresses finding the globally optimal transformation in various 3D registration problems by a unified formulation that integrates common geometric registration modalities (namely point-to-point,point-to-line and point-to-plane). This formulation renders the optimization problem independent of both the number and nature of the correspondences. The main novelty of our proposal is the introduction of a strengthened Lagrangian dual relaxation for this problem, which surpasses previous similar approaches [32] in effectiveness. In fact, even though with no theoretical guarantees, exhaustive empirical evaluation in both synthetic and real experiments always resulted on a tight relaxation that allowed to recover a guaranteed globally optimal solution by exploiting duality theory. Thus, our approach allows for effectively solving the 3D registration with global optimality guarantees while running at a fraction of the time for the state-of-the-art alternative [34], based on a more computationally intensive Branch and Bound method. 1. Introduction The problem of registering 3D geometric data is a classical problem in numerous fields, including computer vision, robotics, photogrammetry or medical imaging [23,12,37]. It seeks the transformation that brings closest together different surfaces according to some meaningful distance function. Consider the ubiquituous scenario in which a system (usually a sensor) returns 3D points {xi}m i=1 of an object and a model of the same is available consisting of 3D primitives {Pi}(typically points, lines and/or planes) [3,14,47,28]. Assuming the correspondences between the Figure 1. This paper presents a unified formulation for the 3D registration with point-to-plane,point-to-line and point-to-point correspondences, and then provides a certifiable globally optimal solution using Lagrangian duality. sets are given, xi↔Pi, the general problem reduces to finding the optimal roto-translation T= (R,t)∈SE(3) as T?= arg min T∈SE(3) m X i=1 dPi(T⊕xi)2.(1) Here SE(3) stands for the usual Special Euclidean group in 3D, T⊕xidenotes the Euclidean transformation of the point xiand dPi(·)is the distance to the primitive Pi[34]. Even with known correspondences, the registration problem (1) is tough to solve in a global fashion due to the nonconvexity of the constraints in the rotation R∈SO(3). A closed-form solution exists only if all the correspondences are point-to-point, as given by Horn et al. [24]. For the other cases, most pipelines resort to local approximations, with the inherent risk of getting stuck in local minima [33,34]. In this paper we address the global optimization of the general 3D registration problem (1). After providing a thorough overview of potential alternatives for this task in Section 2, we present in Section 3a unified formulation integrating point-to-point,point-to-line and point-to-plane correspondences into a single quadratic objective that is a function of the rotation Ronly. Our main contribution is the development of a novel convex relaxation for this formulation under the usual framework of the Lagrangian dual problem, IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii which results on a small constant size semidefinite program (SDP). This novel relaxation, fully characterized in Section 4, turned tight for all the evaluations performed on synthetic and real data, under an extensive range of different problem conditions (Section 5). This allowed us to solve the original non-convex problem in a global fashion using its connection to the convex relaxation through duality theory [4]. In summary, the proposed algorithm Alg. 1provides an iterative optimization framework for 3D registration problems that, unlike current local iterative approaches, is able to provide a certifiable global solution and, in fact, our empirical observation is that it always does so. We remark that this global performance comes without any theoretical guarantees so far. Whereas we acknowledge the importance of the missing formal proof that justifies this behavior, we hope the empirical performance of the proposed relaxation is motivation enough to encourage further exploitation of this approach. 2. Related work As a result of its widespread interest, geometric registration has been object of extensive research for decades now. The number of works aiming at guaranteing global optimality in registration problems is much smaller, though, and focuses mainly on the case of point cloud registration with unknown correspondences [27,35,50,52]. Multimodal registration Registration across diverse 3D geometric primitives arises in a variety of problems and applications. In the Iterative Closest Point (ICP) framework, since the original proposal using point-to-point correspondences by Besl and McKay [3], different types of correspondences have been introduced for improved performance: e.g. point-toplane [14] or plane-to-plane [45,46]. Multimodal correspondences play also an important role in the extrinsic calibration of sensors of different nature, such as camera and lidar [51,49,16,18]. Numerous state-of-the-art SLAM (Simultaneous Localization and Mapping) solutions work with models that include plane and line primitives too [2,47,17,43,28,19]. In the context of multimodal registration a significant effort has been devoted for decades to the resolution of minimal problems [26], e.g. for line-to-plane correspondences [14,49,6,5] and point-to-plane correspondences [21,39,38]. These solutions find application mainly in random sample and consensus (RANSAC) frameworks, but they do not face the more general least squares problem (1). Instead, this optimization problem is traditionally addressed through local linear approximations on the rotation, assuming an initialization close to the solution is given [14], which can easily lead to suboptimal local solutions [33]. Global optimization Next we review the main approaches available in the literature for global optimization [22] and connect them to the registration problem (1) at hand. Analytical solutions A classical approach involves the computation of all the stationary points (among which there is the global minimum). This approach is used by Censi [13] for the global resolution of the 2D registration problem, reducing the problem to solving a 4-th order polynomial equation. However, this approach does not generalize well to the 3D case due to the higher complexity of the rotation space, which produces an explosion in the complexity of the resulting polynomial system [29,30]. Branch and Bound A commonly used tool for NP-hard optimization problems is Branch and Bound (BnB). This is used by Olsson et al. [33,34] for solving the same problem addressed in this paper, yielding a provably global solution. However, the resolution time is notably high (in the order of seconds) due to the exploratory nature and exponential worst-case performance of BnB. Convex relaxations Convex relaxation techniques consider approximate, simpler versions of the problem whose global optimum is much easier to reach. If the approximation is good, the solution of the relaxed problem may then provide valuable information about the original problem. Thus, the main task is to find an appropriate relaxation. For problems such as ours (1), which are affected by the nonconvexity of rotation constraints, a possible relaxation is the search into the convex hull of SO(3). This has proved to give good approximations in rotation synchronization [44], SLAM [41] and 2D/3D registration [25]. Another generic (and successful) tool for providing relaxations of difficult constrained problems is the Lagrangian dual relaxation [4]. This provides particularly good approximations for many problems that can be reformulated as aQuadratically Constrained Quadratic Program (QCQP), where the relaxed problem becomes a Semidefinite Program (SDP) [4,15]. Some problems involving rotations can be characterized as QCQPs, such as Pose Graph Optimization, for which recent literature applying the Lagrangian dual relaxation has shown impressive results finding globally optimal solutions based solely on convex relaxations [11,10,7,40,8]. The Lagrangian dual relaxation has been applied before to the QCQP formulation of 3D registration [32]. In this case the approximation can be very good and even provides the global solution in a certain range of problems, but it deteriorates when the noise level increases or the number IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii of correspondences approaches the minimal cases, that is, when the problem becomes inherently more difficult [33]. Still in the context of the Lagrangian dual relaxation, the relaxation can be strengthened (improved) by introducing additional valid constraints [31, Chap. 13]. This has found applicability in the optimization literature [36,42], and in QCQP problems involving orthonormal constraints, improving the obtained relaxation considerably [1]. To the best of our knowledge, this “trick” has not been considered before for relaxations involving rotations. In this paper, we extend the QCQP formulation of [32] with a whole set of valid quadratic rotational constraints to achieve an improved relaxation that (empirically) returned a globally optimal solution in all problem instances, regardless of the noise level or the number of correspondences. 3. Formulation of multimodal registration In this section we revisit the formulation of the 3D registration problem (1), rewriting it in a suitable form for later applying the Lagrangian dual relaxation. 3.1. Generalized distance function First, we provide a unified formulation of the distances for the point-to-point,point-to-line and point-to-plane correspondences on the basis of the results presented in [34]. In the registration problem (1), the square distance from a 3D point xto a 3D primitive Pminimized in the registration problem (1) is typically that to the closest point in P, defined by d2 P(x) = min y0∈Pkx−y0k2 2.(2) For all the primitives considered here (see Fig. 1) the closest distance problem (2) has a simple closed-form solution that fits a generalized distance function of the form dP(x)2=kx−yk2 C= (x−y)>C(x−y),(3) namely min y0∈Pkx−y0k2 2=(4) =kx−yk2 2=kx−yk2 I3,(point) =k(I−vv>)(x−y)k2 2=kx−yk2 (I−vv>),(line) =(n>(x−y))2=kx−yk2 nn>.(plane) Here y∈R3is any point lying in the primitive1,vis the unit direction vector for a line, niis the unit normal vector for a plane, and C∈S3is a symmetric matrix whose expression depends on the primitive as reflected by (4). These results stem from applying elementary algebra to each primitive [34], but we provide the full proof in the supplementary material for completeness. 1For the point-to-point case yiis the primitive. 3.2. Quadratic formulation and marginalization The distances minimized in the registration problem (1) depend on the transformation T, d2 Pi(T⊕xi)=(T⊕xi−yi)>Ci(T⊕xi−yi).(5) If a matrix representation is chosen for the rotation R, the expression of the transformed point T⊕xiis linear in the elements of Rand t: T⊕xi=Rxi+t= (˜ x>⊗I3) | {z } Xi vec(T),(6) where ˜ x=x>,1>refers to the homogeneized version of x,⊗is the Kronecker product and vec(T)is the vectorization (applied column-wise) of the transformation matrix, vec(T) = vec(R) t.(7) Proof. Supplementary material With this linear parameterization of the transformed point, it is easy to see that the generalized distance to minimize is a quadratic function of τ= vec(T), writable as d2 Pi(T⊕xi) = ˜ τ>N> iCiNi | {z } ˜ Mi ˜ τ,(8) with Ni=˜ x> i⊗I3| − yiand ˜ τ=vec(T) 1> . Because of the quadratic nature of the cost (8) it is possible to accumulate the observations, compressing all the data into a single 13 ×13 matrix term ˜ M: f(T) = m X i=1 d2 Pi(T⊕xi) = ˜ τ> m X i=1 ˜ Mi! | {z } ˜ M ˜ τ.(9) Thanks to this compression step, the size of the following reformulated problem is independent of m. When minimizing the quadratic objective f(T)in (9), the problem can be further reduced if we apply marginalization on the unconstrained part of the unknown T, that is, in the translation t. It is well known from previous work that tcan be derived in terms of R[24,33]. In the quadratic formulation this is straightforward: Lemma 1. The optimal translation for a fixed Ris t?(R) = −˜ M−1 t,t˜ Mt,!t˜ r,˜ r=vec(R) 1.(10) IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Input: List of correspondences {xi↔Pi}m i=1 Output: Global optimum T? /*build compressed quadratic form */ initialize ˜ M←zeros(13); for (xi, Pi)∈ {xi↔Pi}m i=1 do ˜ Mi←buildTerm(xi, Pi);// depends on Pi ˜ M←˜ M+˜ Mi; end /*minimize quadratic objective */ marginalize ˜ Q=˜ M/˜ Mt,t;// Schur complement solve optimal rotation R?using Alg. 2; get optimal translation t?=t?(R?)from (10); return T?= (R?,t?) Algorithm 1: Globally optimal 3D registration Here the subindex tstands for the set of indexes corresponding to translation variables, whereas !tis its complement. The marginalized optimization problem is then f?= min R∈SO(3) ˜ r>˜ Q˜ r |{z} q(˜ r) ,˜ r=vec(R) 1,(11) where the marginalized quadratic form ˜ Q=˜ M/˜ Mt,tis the Schur complement of the block ˜ Mt,tin the matrix ˜ M. Proof. Supplementary material The registration problem is addressed then following the pipeline depicted in Algorithm 1, where the main complexity remains in solving the marginalized problem (11). This is still a non-convex optimization problem in R, for which a convex relaxation is provided in the next session. This relaxation empirically proves to be tight in all the evaluated cases, allowing us to recover a globally optimal solution for the marginalized problem (11). 4. Tight dual relaxation The key ingredient to solve the non-convex problem (11) lies upon an adequate application of Lagrangian duality. We will apply fundamental results from duality theory [4,9], so some basic properties and notions regarding the Lagrangian dual problem are provided in the supplementary material for completeness. Then we present a specific formulation (˜ P) of the constrained problem (11) in Section 4.1 and its corresponding dual problem (D) in Section 4.2. Finally we show how to recover the globally optimal solution in Section 4.3, provided that strong duality holds. The experiments of Section 5show that, empirically, this relaxation is always tight (strong duality holds), even in extreme conditions. Once the whole dual framework has been developed, its implementation is straightforward and the resolution of the marginalized problem (11) is done following the relatively simple pipeline depicted in Algorithm 2. 4.1. Primal problem We will address now the task of formulating the optimization problem (11) in such a manner that the approach described above produces successful results. The constraint Ri∈SO(3) in this problem states that the 3×3rotation matrix Rfulfills the orthonormality and determinant constraints, that is, SO(3) ≡ {R∈R3×3:R> R=I3,det(R) = +1}.(12) In order to apply the Lagrangian dual relaxation, it is particularly appealing to formulate the primal problem as a Quadratically Constrained Quadratic Program (QCPQ). In the usual characterization of SO(3), the orthonormality constraints are all quadratic but the determinant constraint is cubic. Because of this, it has been customary in other problems involving rotations to relax the constraints by dropping the determinant constraint det(R) = +1 and keeping only the orthonormality constraints R>R=I3, which amounts to performing the optimization in O(3) rather than in SO(3). This approach has provided tight relaxations for other problems [11,10,7,40]. For the registration problem however it works well only in a certain range of problems [32]. Duality strengthening Let us now make a quick annotation about an important fact concerning the Lagrangian duality that will be key for the success of our proposal: By construction, every time a new scalar constraint ck+1(·)is introduced into the Lagrangian a new dual variable λk+1 appears and the domain of the dual problem increases its dimension in one. As a result, the bound d? k+1 provided by the new dual problem is at least as good as that of the previous one, d? k≤d? k+1 ≤f?. As a consequece the dual problem is not intrinsic [36,4]: it depends on the particular formulation of the primal problem. In particular, it depends on the specific characterization of the optimization domain or feasible region: Adding appropiate redundant valid constraints has actually shown to be remarkably effective for improving the quality of the dual relaxations in other problems [31, Chap. 13]. Following this idea, our approach is to characterize the feasible set SO(3) with the largest possible amount of quadratic constraints. The chosen constraints need to be linearly independent to introduce any potential improvement [36]. It is important that we keep the complexity of the constraints quadratic in order to maintain the Lagrangian dual problem simple. For the set of orthogonal matrices O(3) it has been shown [1] that a complete set of quadratic constraints is given by the combination of both columnbased and row-based orthogonality constraints: {R>R= I3,RR>=I3}. However, Tron et al. show in [48] that the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Table 1. Table of constraints, Lagrange multipliers and penalizations for Problem ( ˜ P) Constraint type Constraint equation Dual variable Penalization term Orthonormal rows y2I3−RR>=0 Λr=  λ1λ6λ5 λ6λ2λ4 λ5λ4λ3 ∈S3˜ r>˜ Pr(Λr)˜ r Orthonormal columns y2I3−R>R=0 Λc=  λ7λ12 λ11 λ12 λ8λ10 λ11 λ10 λ9 ∈S3˜ r>˜ Pc(Λc)˜ r Handedness R(1) ×R(2) −yR(3) =0λd123 =λ13 λ14 λ15>∈R3 ˜ r>˜ Pdijk (λdijk )˜ r R(2) ×R(3) −yR(1) =0λd231 =λ16 λ17 λ18>∈R3 R(3) ×R(1) −yR(2) =0λd312 =λ19 λ20 λ21>∈R3 Homogeneization 1−y2= 0 γ≡λ22 ∈Rγ+˜ r>˜ Ph(γ)˜ r (a) Orthonormal rows: ˜ Pr(Λr)(b) Orthonormal columns: ˜ Pc(Λc)(c) Handedness: ˜ Pdijk (λdijk )(d) Homogeneization: ˜ Ph(γ) Figure 2. Pattern of the penalization matrices (top row) and dual variables (bottom row) for the different sets of constraints. A coloured cell indicates its value depends (linearly) only on the corresponding dual variable λi. Black cells in (a) and (b) stand for values involving a linear combination of several dual variables. rotation space SO(3) has additional quadratic constraints due to the handedness property that forces the positive unit determinant: Since R∈O(3) ⇒det(R) = ±1, provided that R∈O(3) the positive sign is guaranteed if the matrix columns fulfill the well-known right-hand rule, R(1) ×R(2) =R(3), where R(k)is the k-th column of R. Taking any of the three possible cyclic permutations of the column indexes for the right-hand rule provides exactly three independent quadratic constraints. Altogether we end up with 2·6 + 3 ·3 = 21 scalar rotational constraints: 6for each symmetric matrix constraint from orthonormality, and 3for each vector constraint from the handedness constraint. Problem homogeneization The optimization objective in (11) as well as the gathered constraints are all quadratic functions but in general not homogeneous. It is very convenient for simplifying the derivation of the dual problem to homogeneize the problem by introducing an auxiliary variable ywith the constraint y2= 1 [11,48]. We define the equivalent,homogeneous,strengthened primal problem: min R˜ r>˜ Q˜ r,˜ r=vec(R) y s.t. R>R=y2I3, RR>=y2I3, R(i)×R(j)=yR(k), i, j, k =cyclic(1,2,3) y2= 1. (˜ P) 4.2. Dual problem Once the primal problem has been clearly defined, the derivation of the dual problem is a mechanical work, basically reduced to the derivation of the penalization term corresponding to each constraint. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii The primal problem ( ˜ P) is a QCQP, so the Lagrangian is L(˜ r,˜ λ) = γ+˜ r>(˜ Q+˜ P(˜ λ)) | {z } ˜ Z ˜ r,(13) where the “homogeneous” dual vector ˜ λ=λ>, γ>∈ R22 gathers the dual variables λcorresponding to all the rotation constraints (altogether 21) and the dual variable γfrom the homogeneization constraint y2= 1, shown in Tab. 4. The penalized matrix ˜ Zis the sum of two terms: ˜ Qthat contains all the data from the original problem, and ˜ P(˜ λ)that accumulates all the penalization terms corresponding to the different kinds of constraints: ˜ P(˜ λ) = ˜ Pr(Λr) + ˜ Pc(Λc) + ˜ Pd({λdijk }) + ˜ Ph(γ). This matrix is (by definition) a linear function of the dual variables, and the pattern of the different matrix components can be seen in Fig. 2. A detailed overview of the construction and formulae for ˜ P(˜ λ)is available in the accompanying supplementary material. With this particularly simple expression for the Lagrangian function, the Lagrangian relaxation is an unconstrained problem which can be solved in closed-form as d(˜ λ) = min ˜ rL(˜ r,˜ λ) = min ˜ rγ+˜ r>˜ Z˜ r(14) =(γif ˜ Z<0, −∞ otherwise.(15) The Lagrangian relaxation is unbounded below unless the penalized matrix ˜ Zis positive semidefinite (PSD). As a result, the maximization of the dual objective d(˜ λ)can be safely restricted to those vectors ˜ λpreserving the positive semidefiniteness of ˜ Z. Thus, the dual problem corresponding to the homogeneous primal problem ( ˜ P) is a Semidefinite Program (SDP): d?= max ˜ λ γ, s.t. ˜ Z(˜ λ) = ˜ Q+˜ P(˜ λ)<0.(D) This problem is convex and off-the-shelf specialized solvers exist for it [20]. 4.3. Primal-via-dual resolution In this section we begin by assuming that the duality gap for our primal-dual pair is zero (we will see in the experiments that this assumption always holds in practice). By duality theory [4], ˜ r?must be a minimizer of the Lagrangian (13) evaluated at ˜ λ?, x?= arg min x L(x,λ?)⇒(˜ r?)>˜ Z?˜ r?= 0.(16) Since Z?<0, this means that the primal optimum ˜ r?must lie in the nullspace of ˜ Z?: x?= arg min x L(x,λ?)⇒˜ r?∈null( ˜ Z?).(17) Input: Marginalized quadratic form q(˜ r)(11) Output: Global optimum R?(if no error occurs) build symbolic penalized matrix ˜ Z(˜ λ)←˜ Q+˜ P(˜ λ); solve the SDP problem (D) to get ˜ λ?and ˜ Z?;// CVX get nullspace V= null( ˜ Z?); ASSERT(rank(V) == 1 ); dehomogeneize ˜ r?←V/V(end); ASSERT(q(˜ r?)−d?≈0) ; // strong duality R?˜ r?;// reshape return R Algorithm 2: Convex global optimizer If the nullspace has rank 1, the solution ˜ r?is recovered up to a scale factor. Then, since the solution must also fulfill the original constraints in the primal problem ( ˜ P), we fully determine the solution by setting y= 1, which in practice reduces to dehomogeneizing the solution ˜ r?. Then with the obtained primal solution r?we can check that our initial assumption holds and the duality gap is effectively zero, d?=f(R?). As we will see in the experiments, both conditions rank(null( ˜ Z?)) = 1 and d?=f? were fulfilled for absolutely all the experiments considered, even under the most extreme situations (in terms of noise and number of correspondences). 5. Experiments In this section we show that in practice the strong duality assumption holds in our primal-dual formulation for any considered problem. As a result, using Algorithm 1it is always possible to recover the primal optimal solution, and we do this at a fraction of the time necessary for the more complex exploratory techniques. The SDP (D) is solved using CVX [20]. We assess the performance of our method, Ours, in both synthetic and real data and compare it to that of two different state-of-the-art approaches for solving the 3D registration problem globally: a provably optimal exploratory approach based on Branch and Bound, BnB [34], and a dualbased approach which provides a relaxation, referred to as Olsson [32]. But first of all, we define the main parameters that will characterize our general registration problem (1) as well as the metrics used in the assessment of performance. Effective number of correspondences The number of measurements has a notable impact in the complexity of the registration problem. Namely, Olsson and Eriksson [32] justify and illustrate that the problem becomes easier to solve as the number of measurements increase. Following this intuition, the most difficult registration problems should be those close to minimal cases. In orIEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii der to measure how close a problem is to being minimal we consider the unifying framework presented by Ramalingam and Taguchi [38]: Point-to-point and point-to-line correspondences can be transformed into equivalent sets of 3and 2point-to-plane correspondences, respectively. We define then the effective number of correspondences,ˆm, as the equivalent number of independent point-to-plane correspondences ˆm= 3mpoint + 2mL+ 1mΠ.(18) The minimum value of ˆmfor which a general 3D registration problem may have a unique global minimum is 7[38]. Geometry Even if the effective number of correspondences is higher than 6, degeneracies and symmetries with multiple global minima may still occur depending on the geometric distribution of the correspondences. These cases are identified as well in [38]. We took care during the evaluation on both synthetic and real data to discard these degenerate configurations where the true global minimum cannot be found from the data only. Measurement noise This models the quality of the measurements (we consider no oultliers). If there is no noise, a simple linear relaxation would provide the global solution. Then, as the noise level σincreases the problem becomes harder to solve [32]. Metrics Several metrics can be used to measure the effectiveness of the relaxation approaches, Ours and Olsson. Due to space issues we choose to show here the most significant metric, which is in our view the optimality ratio, that is, the percentage of cases in which a globally optimal solution was attained. A solution is considered globally optimal if the suboptimality gap ∆ = f−f?is zero up to numerical precision. The global minimum f?can be found from a provably global algorithm such as BnB or, as we will see next, also from our tight relaxation. In order to measure the computational performance of the different methods we also plot the resolution times. In particular, we use shaded error bars to display the median values plus the 1st and 3rd quartiles reflecting the distribution of the values. The statistics shown in the figures were generated from a population of 100 registration problems in each case. Other interesting metrics, as well as other parameters ranges beyond those displayed in this document are available in the supplementary material. These additional results support the same conclusions about the methods reached in this document and previous works [32,34]. 0 0.1 0.2 0.3 0.4 0.5 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (a) (b) Figure 3. Results in the synthetic scenario for increasing level of noise σ, with ˆm= 10. (a) The effectivenes of the reference relaxation Olsson [32] drops steadily as the level of noise increases. Ours remained optimal in all cases. (b) Both convex relaxations Ours and Olsson take roughly the same time, whereas BnB gets two orders of magnitude slower. 1 10 100 1000 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (a) (b) Figure 4. Results in the synthetic scenario for extremely severe σ noise conditions and reduced amount of effective number measurements ˆm= 7. These conditions are of reduced practical interest, but are used to show the general robustness of the tight condition for our relaxation compared to the reference relaxation Olsson [32]. In these regimes BnB suffers of numerical issues. 5.1. Synthetic data For obtaining synthetic problems we generated a set of random model primitives that added up to ˆmeffective correspondences. Similarly to [32], each primitive was determined by randomly taking a point inside a sphere of radius 10 m (plus a random unit direction for the case of lines and planes). Then a “measured” point was randomly picked from the set defined by the primitive and we corrupted it with a Gaussian noise of standard deviation σ. We show the behaviour w.r.t. the noise level in Fig. 3in a challenging case with ˆm= 10. The results for a varying value of ˆmwere similar to those shown for the real data in Fig. 6(a). These and other evaluations for different ranges of parameters are shown in the supplementary material. To sum up, our method attained the globally optimal solution in all the considered cases, without exception, even in the most severe cases where the number of features ˆmremained almost minimal and the measurement noise σwas raised way beyond any expectable value in real scenarios (see Fig. 4). The reference relaxation Olsson in contrast was rarely tight in the challenging scenarios and returned a suboptimal solution. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Figure 5. The experimental setup for obtaining real data, taken from [32], and optimal registration of the measurements to the model. 5.2. Real data Thanks to the corresponding authors, it was possible to exploit the same real data employed in the references [32,34]. Their experimental setup consisted of using a MicroScribe-3DLX 3D scanner to measure the 3D coordinates of some points on the real object, as shown in Fig. 5. For the Space Station model 49 points were measured on different primitives of the object, namely 27 on planes, 12 in lines and 10 on corner points. The registration of the complete set to the computer model is shown in 5. In [32] it was shown that for this complete problem Olsson works fine, attaining the global optimum as BnB but in much less time. We use the same data to generate a more extensive set of challenging real problems. In this case, we can produce significantly more difficult problems by sampling a smaller set of measurements from the data: We choose different combinations of point-to-point,point-to-line and point-to-plane correspondences that result in a particular effective number of measurements ˆm. The precision σof the measurements in this case is fixed by the sensor, with errors of about 0.5 millimeters according to the authors of the dataset. The obtained results are consistent with those observed from the evaluation on synthetic data. The optimality ratio is displayed in Fig. 6(a). The behaviour for Olsson w.r.t. the parameter ˆmwas consistent with that predicted in the original work [32]: It hardly attained global optimality in near-minimal cases, and the performance improved steadily with the increase on ˆm. Meanwhile, our approach succeeded again in all the cases, always returning the globally optimal solution. Again, both convex relaxations Ours and Olsson took roughly the same time, whereas BnB is two orders of magnitude slower (see Fig. 6(b)). In conclusion, our approach attained the same optimality ratio guaranteed so far only for the provably optimal BnB method, whereas Ours was two orders of magnitudes faster than BnB. 7 8 9 10 11 12 13 14 15 ˆm 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (a) (b) Figure 6. Results in the real scenario for a low effective number of measuremens ˆm. (a) The effectivenes of the reference relaxation Olsson [32] improves with the increase in ˆm, yet it hardly attains optimality in most cases whereas ours remained optimal in all cases. (b) Both convex relaxations Ours and Olsson take roughly the same time, whereas BnB is two orders of magnitude slower. 6. Conclusions We have presented a unified formulation for the 3D registration problem involving point-to-point, point-to-line and point-to-plane correspondences that compresses the objective into a single quadratic function of the rotation. Thanks to its generality and flexibility, this formulation should have the potential to introduce further types of correspondences beyond those explored in this work. The remaining optimization problem has then been characterized as a Quadratically Constrained Quadratic Program. Exploiting a full set of quadratic rotational constraints we obtain a Lagrangian dual relaxation from which a globally optimal solution could be recovered in 100% of the tested cases, although it remains open the theoretical question of why strong duality holds for this relaxation in virtually all cases. Finally, even though the current approach is already 2orders of magnitude faster than the competing BnB approach, we are just taking an off-the-shelf generic SDP solver so this performance could be improved further by using specialized solvers that exploit the low-rank structure of the underlying SDP problem. Acknowledgements The authors would like to thank Carl Olsson for kindly providing code and data for the experiments in [32]. This work has been supported by a University of Malaga travel grant, the Spanish grant program FPU14/06098 and the project PROMOVE (DPI2014-55826-R), funded by the Spanish Government and the ”European Regional Development Fund”. References [1] K. Anstreicher and H. Wolkowicz. On Lagrangian relaxation of quadratic matrix constraints. SIAM J. Matrix Anal. Appl., 22(1):41–55, 2000. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii [2] A. Bartoli and P. Sturm. Structure-from-motion using lines: Representation, triangulation, and bundle adjustment. Comput. Vis. Image Underst., 100(3):416–441, 2005. [3] P. J. Besl and N. D. McKay. Method for registration of 3D shapes. In Robot. Tentat., pages 586–606. International Society for Optics and Photonics, 1992. [4] S. Boyd and L. Vandenberghe. Convex optimization. 2004. [5] J. Briales and J. Gonzalez-Jimenez. A Minimal Closedform Solution for the Perspective Three Orthogonal Angles (P3oA) Problem: Application To Visual Odometry. J. Math. Imaging Vis., pages 1–18, 2015. [6] J. Briales and J. Gonzalez-Jimenez. A minimal solution for the calibration of a 2D laser-rangefinder and a camera based on scene corners. In Intell. Robot. Syst. (IROS), 2015 IEEE/RSJ Int. Conf., pages 1891–1896, sep 2015. [7] J. Briales and J. Gonz´ alez-Jim´ enez. Fast Global Optimality Verification in 3D SLAM. In Int. Conf. Intell. Robot. Syst. IEEE/RSJ, 2016. [8] J. Briales and J. Gonz´ alez-Jim´ enez. Initialization of 3D Pose Graph Optimization using Lagrangian duality. In Int. Conf. on Robotics and Automation (ICRA). IEEE, 2017. [9] G. C. Calafiore and L. El Ghaoui. Optimization Models. Cambridge university press, 2014. [10] L. Carlone, G. C. Calafiore, C. Tommolillo, and F. Dellaert. Planar Pose Graph Optimization: Duality, Optimal Solutions, and Verification. IEEE Trans. Robot., 32(3):545–565, 2016. [11] L. Carlone, D. M. Rosen, G. Calafiore, J. J. Leonard, and F. Dellaert. Lagrangian duality in 3D SLAM: Verification techniques and optimal solutions. In Intell. Robot. Syst. (IROS), 2015 IEEE/RSJ Int. Conf., pages 125–132, sep 2015. [12] U. Castellani and A. Bartoli. 3d shape registration. In 3D Imaging, Anal. Appl., pages 221–264. Springer, 2012. [13] A. Censi. An ICP variant using a point-to-line metric. Robot. Autom. 2008. ICRA 2008. IEEE, 2008. [14] Y. Chen and G. Medioni. Object modelling by registration of multiple range images. Image Vis. Comput., 10(3):145–155, 1992. [15] Y. Ding. On efficient semidefinite relaxations for quadratically constrained quadratic programming. 2007. [16] E. Fern´ andez-Moral, J. Gonz´ alez-Jim´ enez, and V. Ar´ evalo. Extrinsic calibration of 2D laser rangefinders from perpendicular plane observations. Int. J. Rob. Res., 34(11):1401– 1417, 2015. [17] E. Fern´ andez-Moral, W. Mayol-Cuevas, V. Ar´ evalo, and J. Gonz´ alez-Jim´ enez. Fast place recognition with planebased maps. In Robot. Autom. (ICRA), 2013 IEEE Int. Conf., pages 2719–2724. IEEE, 2013. [18] R. Gomez-Ojeda, J. Briales, E. Fern´ andez-Moral, and J. Gonz´ alez-Jim´ enez. Extrinsic Calibration of a 2D LaserRangefinder and a Camera based on Scene Corners. In IEEE Int. Conf. Robot. Autom., pages 3611–3616. IEEE, 2015. [19] R. Gomez-ojeda, J. Briales, and J. Gonzalez-jimenez. PLSVO: Semi-Direct Monocular Visual Odometry by Combining Points and Line Segments. [20] M. Grant and S. Boyd. {CVX}: Matlab Software for Disciplined Convex Programming, version 2.1. \url{http://cvxr.com/cvx}, 2014. [21] W. Grimson and T. Lozano-Perez. Model-based recognition and localization from sparse range or tactile data. Int. J., 1984. [22] R. Hartley and F. Kahl. Optimal Algorithms in Multiview Geometry. In Comput. Vis. ACCV 2007, pages 13–34. Springer Berlin Heidelberg, Berlin, Heidelberg. [23] R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision. Cambridge University Press, New York, NY, USA, 2 edition, 2003. [24] B. K. P. Horn. Closed-form solution of absolute orientation using unit quaternions. JOSA A, 4(4):629–642, 1987. [25] Y. Khoo and A. Kapoor. Non-iterative rigid 2D/3D point-set registration using semidefinite programming. IEEE Trans. Image Process., 25(7):2956–2970, 2016. [26] Z. Kukelova, M. Bujnak, and T. Pajdla. Automatic generator of minimal problem solvers. Comput. VisionECCV 2008, 2008. [27] H. Li and R. Hartley. The 3D-3D registration problem revisited. In 2007 IEEE 11th Int. Conf. Comput. Vis., pages 1–8. IEEE, 2007. [28] L. Ma, C. Kerl, J. St¨ uckler, and D. Cremers. CPASLAM: Consistent Plane-Model Alignment for Direct RGBD SLAM. [29] F. M. Mirzaei and S. I. Roumeliotis. Globally optimal pose estimation from line correspondences. In IEEE Int. Conf. Robot. Autom. (ICRA), 2011, pages 5581–5588. Ieee, may 2011. [30] F. M. Mirzaei and S. I. Roumeliotis. Optimal estimation of vanishing points in a Manhattan world. In IEEE Int. Conf. Comput. Vis. (ICCV), 2011, pages 2454–2461. IEEE, 2011. [31] Y. Nesterov, H. Wolkowicz, and Y. Ye. Semidefinite programming relaxations of nonconvex quadratic optimization. In Handb. semidefinite Program., pages 361–419. Springer, 2000. [32] C. Olsson and A. Eriksson. Solving quadratically constrained geometrical problems using lagrangian duality. In Pattern Recognition, 2008. ICPR 2008. 19th Int. Conf., pages 1–5. IEEE, 2008. [33] C. Olsson, F. Kahl, and M. Oskarsson. The registration problem revisited: Optimal solutions from points, lines and planes. In Comput. Vis. Pattern Recognition, 2006 IEEE Comput. Soc. Conf., volume 1, pages 1206–1213. IEEE, 2006. [34] C. Olsson, F. Kahl, and M. Oskarsson. Branch-and-Bound Methods for Euclidean Registration Problems. IEEE Trans. Pattern Anal. Mach. Intell., 31(5):783–794, 2009. [35] C. Papazov and D. Burschka. Stochastic global optimization for robust point set registration. Comput. Vis. Image Underst., 115(12):1598–1609, 2011. [36] P. A. Parrilo and S. Lall. Semidefinite programming relaxations and algebraic optimization in control. Eur. J. Control, 9(2):307–321, 2003. [37] F. Pomerleau, F. Colas, and R. Siegwart. A review of point cloud registration algorithms for mobile robotics. Found. Trends Robot., 4(1):1–104, 2015. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii 5. Lagrangian derivation This section shows the construction of the Lagrangian function for the equivalent,homogeneous,strengthened primal problem: min R,y ˜ r>˜ Q˜ r | {z } q(˜ r) ,˜ r=vec(R) y(˜ P) s.t. R>R=y2I3,(43) RR>=y2I3,(44) R(i)×R(j)=yR(k),(i, j, k)∈(123),(45) y2= 1.(46) We use the notation (123) to refer to the set of cyclic permutations: (123) = {(1,2,3),(2,3,1),(3,1,2)}. The Lagrangian of the primal problem ( ˜ P) is then built by introducing the equality constraints as penalization terms L(˜ r,˜ λ) = q(˜ r)[objective] (47) +hΛr, y2I3−RR>i[orthonormal rows] (48) +hΛc, y2I3−R>Ri[orthonormal columns] (49) +X (i,j,k)∈(123) hλdijk ,R(i)×R(j)−yR(k)i[handedness] (50) +hγ, 1−y2i.[homogeneization] (51) Every penalization term is written as the inner product of a Lagrange multiplier with a set of constraints. If the constraint is scalar (51) the inner product is a simple scalar product ha, bi=ab. For vector constraints (50) it becomes the dot product for vectors, ha,bi=a>b. Finally if the constraints form a matrix (48)(49), we use the inner product for matrices hA,Bi= tr(A>B). Since the constraints in the primal problem ( ˜ P) are (homogeneous) quadratic functions of R, a linear combination of them will still be quadratic w.r.t. Rand the Lagrangian can be refactored into L(˜ r,˜ λ) = γ+˜ r>(˜ Q+˜ P(˜ λ))˜ r.(52) The offset γoriginates from the homogeneization constraint (51), the only constraint with a non-homogeneous term. The augmented vector of dual variables ˜ λ∈R22 is formed from the concatenation of all the dual variables, and the penalization matrix ˜ Pis the sum of the different components coming from the constraints in ( ˜ P), namely orthonormality of rows (r) (48), orthonormality of columns (c) (49), positive determinant or handedness (d) (50) and homogeneization (h) (51): ˜ P(˜ λ) = ˜ Pr(Λr) + ˜ Pc(Λc) + ˜ Pd({λdijk } |{z } λd ) + ˜ Ph(γ).(53) Next we will fully characterize these penalization terms. A complete summary of the constraints, associated dual variables and penalizations can be seen in Tab. 5. A visualization of the pattern for the resulting matrices ˜ P(·)is available in Fig. 2. 5.1. Orthonormality of rotation rows (48) The matrix constraint for the orthonormality of rows, RR>=I3, provides the penalization term hΛr, y2I3−RR>i= tr(Λ> r(y2I3−RR>)) (54) [using Λr=Λ> rand cyclic property of tr(·)]=y2tr(Λr)−tr(R>ΛrR)(55) [using tr(IX>BY )(81)] =y2tr(Λr)−r>(I3⊗Λr)r(56) [reordering into block matrix] =˜ r>˜ Pr(Λr)˜ r.(57) IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Table 3. Table of constraints and penalizations for ( ˜ P) Constraint type Constraint equation Dual variable Penalization term Orthonormal rows y2I3−RR>=0 Λr=  λ1λ6λ5 λ6λ2λ4 λ5λ4λ3  ∈S3˜ r>˜ Pr(Λr)˜ r, see (58) Orthonormal columns y2I3−R>R=0 Λc=  λ7λ12 λ11 λ12 λ8λ10 λ11 λ10 λ9  ∈S3˜ r>˜ Pc(Λc)˜ r, see (63) Handedness R(1) ×R(2) −yR(3) =0λd123 =λ13 λ14 λ15>∈R3 ˜ r>˜ Pdijk (λdijk )˜ r, see (74) R(2) ×R(3) −yR(1) =0λd231 =λ16 λ17 λ18>∈R3 R(3) ×R(1) −yR(2) =0λd312 =λ19 λ20 λ21>∈R3 Homogeneization 1−y2= 0 γ≡λ22 ∈Rγ+˜ r>˜ Ph(γ)˜ r, see (77) (a) Orthonormal rows, see (58) (b) Orthonormal columns, see (63) (c) Handedness, see (74) (d) Homogeneization, see (77) Figure 2. Sparsity pattern of the penalization matrices (top row) and dual variables (bottom row) due to the different sets of constraints. A coloured cell indicates its value depends (linearly) only on the corresponding dual variable λi. Black cells (in (10,10)) stand for values involving a linear combination of several dual variables. The matrix Lagrange multiplier Λris symmetric because only 6effective different scalar constraints stand here. The quadratic penalization matrix ˜ Pris shown in Fig. 2(a) and takes the expression ˜ Pr(Λr) =  −I3⊗Λr09×1 01×9tr(Λr) .(58) 5.2. Orthonormality of rotation columns (49) Similar to the previous orthonormality constraint but with exchanged transposes, R>R=I3, the orthonormality of rotation columns provides the penalization hΛc, y2I3−R>Ri= tr(Λ> c(y2I3−R>R)) (59) =y2tr(Λc)−tr(Λ> cR>R)(60) [using tr(A>X>IY )(81)] =y2tr(Λc)−r>(Λc⊗I3)r(61) [reordering into block matrix] =˜ r>˜ Pr(Λc)˜ r.(62) IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Again, there are only 6different constraints, so Λc∈S3. The quadratic penalization matrix ˜ Pcis shown in Fig. 2(b) and takes the expression ˜ Pc(Λc) =  −Λc⊗I309×1 01×9tr(Λc) .(63) 5.3. Determinant constraints (50) The right-hand rule applied to the columns of the rotation matrix Rtakes the form R(i)×R(j)=yR(k),(64) where R(u)refers to the u-th column in the matrix R. Let us introduce the canonical vector eu∈R3, defined as the vector that is zero everywhere except for the u-th entry which is 1. It will be useful as well to define the canonical matrix euv as the matrix which has entry 1at the position (u, v)and is zero everywhere else. Note euv =eue> v. Then the handedness constraint can be written in a more convenient form as hλdijk ,R(i)×R(j)−yR(k)i=λ> dijk ((Rei)×(Rej)−y(Rek)) = (65) [using a·(b×c) = −b·(a×c)]=−(Rei)>(λdijk ×(Rej)) −yλ> dijk (Rek)(66) [using a×b= [a]×b]=−e> iR>λdijk ×Rej−yλ> dijk Rek(67) [using tr(a) = a, a ∈R]=−tr(e> iR>λdijk ×Rej)−ytr(λ> dijk Rek)(68) [using cyclic property of tr(·)]=−tr(eje> iR>λdijk ×R)−ytr(ekλ> dijk R)(69) [reordering transposes] =−tr(e> ijR>λdijk ×R)−ytr((λdijk e> k)>R)(70) [using tr(A>X>IY )(81) and tr(A>B)(80)] =−r>(eij ⊗λdijk ×)r−yvec(λdijk e> k)>r(71) [using vec(ab>)(79)] =−r>(eij ⊗λdijk ×)r−y(ek⊗λdijk )>r(72) [reordering into block matrix] =˜ r>˜ Pdijk (λdijk )˜ r.(73) The (symmetric) quadratic penalization matrix ˜ Pdijk (λdijk )takes the expression ˜ Pdijk = −eij ⊗λdijk ×−(ek⊗λdijk ) 01×90 .(74) Note this expression is generic for (i, j, k)∈(123) = {(1,2,3),(2,3,1),(3,1,2)}. The complete penalization matrix ˜ Pd(λd) = ˜ Pd123 (λd123 ) + ˜ Pd231 (λd231 ) + ˜ Pd312 (λd312 )from all the determinant constraints, after taking an equivalent symmetric matrix for the representation, is shown in Fig. 2(c). 5.4. Homogeneization constraint (51) Finally, the simple scalar constraint y2= 1 yields hγ, 1−y2i=γ−γy2(75) =γ+˜ r>˜ Ph(γ)˜ r(76) For compatibility with the rest of quadratic expressions this has been written in matrix form too, where ˜ Ph(γ) = −γe10,10 is the matrix with all zeros but the lower right corner, where it contains the value −γ: ˜ Ph(γ) =   09×909×1 01×9−γ .(77) 6. Additional experimental results Due to the lack of space, only the most relevant experimental results and metrics were shown in the paper. In this section we include all the obtained results comparing the BnB solution [4], and the dual-based relaxations, Ours and Olsson [3]. The results support the claims in the main document. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii Evaluation metrics We used several metrics to evaluate the performance of the methods. Next we describe them and give some notes on the representation used for the results. All the statistics were obtained from 200 samples. Suboptimality gap The suboptimality gap is defined as ∆ = f−f?, where f?is the minimum objective value attained in the global solution. Since our solution always provided a globally optimal solution (in a certifiable way, since f=d?⇒f= f?), it is possible to plot this value for all the methods giving a clear reference for the performance: A solution is optimal if ∆=0, suboptimal otherwise. We displayed the suboptimality gap in the oncoming figures using boxplots. To show as much information as possible, we also superimposed a sorted set of points displaying the original values underlying each boxplot. Note that for the globally optimal methods, Ours and BnB, the boxplots and points degenerate to a single line in ∆=0. In order to cope with the high variability of the suboptimality values for the different methods, we applied a custom scale in the Y axis, where the log(1 + x)of the values is displayed in order to better visualize different orders of magnitude while keeping the lower threshold ∆=0in the origin. Global optimality ratio The global optimality ratio provides a refined view of suboptimality metric, where we directly show a bar plot with the percentage of cases in which a globally optimal solution (∆ = 0) was attained. Time The computational performance of the different algorithms is evaluated by measuring the CPU time necessary to solve each problem. Note these times are just for orientation, as the implementations were not highly optimized. Nevertheless, the conclusions when comparing the different methods remain valid. The obtained values are represented with a shaded error bar plot displaying the median values plus a band defined by the 1st and 3rd quartiles to reflect the distribution of the values. Note that we always used a logarithmic scale for the time values, in order to better compare all the considered methods. Problem parameters As discussed in the paper, the main problem parameters defining the overall difficulty of solving the optimization problem globally are the effective number of correspondences ˆmand the level of the noise σaffecting the measurements. A more detailed explanation of these is available in the main document. 6.1. Evaluation on synthetic data 6.1.1 Challenging range of ˆm The evaluation in a challenging range of ˆmvalues, going from the almost minimal case ˆm= 7 up to ˆm= 15, shows that as the number of correspondences increases the level of difficulty of the problem decreases (see Fig. 3). This is consistent with the claim in [3]. This explains why the looser relaxation Olsson has better performance as ˆmgrows. This trend holds for different levels of noise σ, although as expected the general difficulty keeps higher when the noise is larger, making the performance of Olsson worse. In all cases, as the provably optimal method BnB, our method provided a globally optimal solution. However, the computational burden for the BnB method was notably higher, and in the easier scenario (σ= 0.1) it increased with the decrease in the number of correspondences. This is consistent with the remarks in [4], where the nearminimal cases were tagged as harder to explore in the Branch and Bound framework. The time for the dual-based algorithms (Ours and Olsson) remained practically constant in all cases and was notably lower than for BnB (around 2orders of magnitude). IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii (a) Suboptimality gap, σ= 0.1m 7 8 9 10 11 12 13 14 15 ˆm 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (b) % of optimality, σ= 0.1m (c) Times, σ= 0.1m (d) Suboptimality gap, σ= 0.5m 7 8 9 10 11 12 13 14 15 ˆm 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (e) % of optimality, σ= 0.5m (f) Times, σ= 0.5m (g) Suboptimality gap, σ= 1 m 7 8 9 10 11 12 13 14 15 ˆm 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (h) % of optimality, σ= 1 m (i) Times, σ= 1 m Figure 3. Evaluation in a challenging range of effective number of correspondences, from the almost minimal case ˆm= 7 to ˆm= 15. Different levels of noise were considered: A low value of σ= 0.1m (top row), a middel value of σ= 0.5m (middle row), and a more challenging scenario with σ= 1 m (bottom row). 6.1.2 Evaluation w.r.t. noise level As in [3], we provide a continuous analysis on how the performance of the methods are affected for different levels of noise σwithin usual application values, going from the ideal noise-free case with σ= 0 to the significantly high value of σ= 0.5 m. This evaluation is repeated several times under a fixed effective number of correspondences, namely ˆm={7,10,14,21}. The whole set of results can be seen in Fig. 4. As expected, the general trend is that the difficulty of the problem increases with the noise level σ, making it harder for Olsson to reach the global solution, and also rendering the resolution with BnB slower. In all cases our solution provided the globally optimal solution, at a constant low resolution time. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii (a) Suboptimality gap, ˆm= 7 0 0.1 0.2 0.3 0.4 0.5 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (b) % of optimality, ˆm= 7 (c) Times, ˆm= 7 (d) Suboptimality gap, ˆm= 10 0 0.1 0.2 0.3 0.4 0.5 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (e) % of optimality, ˆm= 10 (f) Times, ˆm= 10 (g) Suboptimality gap, ˆm= 14 0 0.1 0.2 0.3 0.4 0.5 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (h) % of optimality, ˆm= 14 (i) Times, ˆm= 14 (j) Suboptimality gap, ˆm= 21 0 0.1 0.2 0.3 0.4 0.5 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (k) % of optimality, ˆm= 21 (l) Times, ˆm= 21 Figure 4. Evaluation w.r.t. noise level from σ= 0 to σ= 0.5m, for different effective numbers of measurements: ˆm={7,10,14,21}. Note the decreasing scale in the suboptimality gap (first column of figures) as the effective number of correspondences ˆmincreases. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii 6.1.3 Evaluation in extreme conditions Finally, we used the synthetic framework to test our method in a substantially difficult scenario. For that purpose, we took the most challenging conditions for the registration problem, that is almost minimal number of correspondences ˆm= 7 and very high noise level (far beyond any expectable value found in a real problem). Specifically, we increased the noise level σ exponentially from σ= 1 m to σ= 1000 m. With such a high level of noise the results of the problem would of little use in a real application, but again we do this just with the main aim of testing the behaviour and robustness of all the algorithms from a purely mathematical point of view. The results, displayed in Fig. 5, show that under these severe conditions the relaxation used by Olsson never worked in the higher regimes of noise, and even the BnB method incurred in some suboptimalities probably due to numerical issues and exit conditions in the algorithm. Meanwhile, Ours worked in 100% of the cases without exception. (a) Suboptimality gap, ˆm= 7 1 10 100 1000 σ(m) 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (b) % of optimality, ˆm= 7 (c) Times, ˆm= 7 Figure 5. Evaluation in the most severe regime of parameters: Lowest possible number of correspondences ˆm= 7 and significantly high measurement noise. Note the different scales in the axes. 6.2. Evaluation on real data We use the same real dataset used in [3]. Following the procedure explained in the main document, we were able to generate a large set of different registration problems by taking different effective numbers of correspondences ˆm. The noise level was intrinsic to the sensor used for the measurements, so it is out of our control. As a result, as in Section 6.1.1, the evaluation on real data considers an increasing value for ˆm. Note however that, even though the sensor noise should be well below the synthetic noise considered in the top row of Fig. 3,σ= 0.1m, the results are closer to those of a much higher synthetic noise (compare to σ= 0.5m in Fig. 3). This may be a consequence of the non-Gaussian noise in the real measurements, which might be turning the problem even harder. In any case, our method still reached a globally optimal solution in all cases, regardless of the effective number of measurements ˆm, and at a fraction of the time for BnB. (a) Suboptimality gap 7 8 9 10 11 12 13 14 15 ˆm 0 10 20 30 40 50 60 70 80 90 100 % optimal BnB Ours Olsson (b) % of optimality (c) Times Figure 6. Evaluation in real data taken from the Space Station dataset [3]. Different subsets of the complete dataset are sampled to get more challenging registration problems with a variable effective number of correspondences ˆm. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii A. Some matrix calculus Next we provide some properties related to the vectorization of matrix expressions. These relations were taken from [2]: vec(AXB)=(B>⊗A) vec(X)(78) vec(ab>) = b⊗a(79) tr(A>B) = vec(A)>vec(B)(80) tr(A>X>BY ) = vec(X)>(A⊗B) vec(Y)(81) References [1] S. Boyd and L. Vandenberghe. Convex optimization. 2004. 5 [2] P. L. Fackler. Notes on Matrix Calculus. pages 1–14, 2005. 3,13 [3] C. Olsson and A. Eriksson. Solving quadratically constrained geometrical problems using lagrangian duality. In Pattern Recognition, 2008. ICPR 2008. 19th Int. Conf., pages 1–5. IEEE, 2008. 8,9,10,12 [4] C. Olsson, F. Kahl, and M. Oskarsson. Branch-and-Bound Methods for Euclidean Registration Problems. IEEE Trans. Pattern Anal. Mach. Intell., 31(5):783–794, 2009. 8,9 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017, Honolulu, Hawaii