Global registration of multiple point clouds using semidefinite programming

Introduction

The problem of point-cloud registration comes up in computer vision and graphics , and in distributed approaches to molecular conformation and sensor network localization . The registration problem in question is one of determining the coordinates of a point cloud $P$ from the knowledge of (possibly noisy) coordinates of smaller point cloud subsets (called patches) $P_{1},\ldots,P_{M}$ that are derived from $P$ through some general transformation. In certain applications , one is often interested in finding the optimal transforms (one for each patch) that consistently align $P_{1},\ldots,P_{M}$. This can be seen as a sub-problem in the determination of the coordinates of $P$ .

In this paper, we consider the problem of rigid registration in which the points within a given $P_{i}$ are (ideally) obtained from $P$ through an unknown rigid transform. Moreover, we assume that the correspondence between the local patches and the original point cloud is known, that is, we know beforehand as to which points from $P$ are contained in a given $P_{i}$. In fact, one has a control on the correspondence in distributed approaches to molecular conformation and sensor network localization . While this correspondence is not directly available for certain graphics and vision problems, such as multiview registration , it is in principle possible to estimate the correspondence by aligning pairs of patches, e.g., using the ICP (Iterative Closest Point) algorithm .

The problem is to infer $O$ and $t$ from the above equations. To uniquely determine $O$ and $t$, one must have at least $N\geq d+1$ non-degenerate pointsBy non-degenerate, we mean that the affine span of the points is $d$ dimensional.. In this case, $O$ can be determined simply by fixing the first equation in (1) and subtracting (to eliminate $t$) any of the remaining $d$ equations from it. Say, we subtract the next $d$ equations:

By the non-degeneracy assumption, the matrix on the right of $O$ is invertible, and this gives us $O$. Plugging $O$ into any of the equations in (1), we get $t$.

In practical settings, (1) would hold only approximately, say, due to noise or model imperfections. A particular approach then would be to determine the optimal $O$ and $t$ by considering the following least-squares program:

in which $x_{c}=(x_{1}+\cdots+x_{N})/N$ and $y_{c}=(y_{1}+\cdots+y_{N})/N$ are the centroids of the respective point clouds. The optimal translation is $t^{\star}=y_{c}-O^{\star}x_{c}$.

The fact that two-patch registration has a closed-form solution is used in the so-called incremental (sequential) approaches for registering multiple patches . The most well-known method is the ICP algorithm (note that ICP uses other heuristics and refinements besides registering corresponding points). Roughly, the idea in sequential registration is to register two overlapping patches at a time, and then integrate the estimated pairwise transforms using some means. The integration can be achieved either locally (on a patch-by-patch basis), or using global cycle-based methods such as synchronization . More recently, it was demonstrated that, by locally registering overlapping patches and then integrating the pairwise transforms using synchronization, one can design efficient and robust methods for distributed sensor network localization and molecular conformation . Note that, while the registration phase is local, the synchronization method integrates the local transforms in a globally consistent manner. This makes it robust to error propagation that often plague local integration methods .

2. Multi-patch Registration

In this paper, we assume that the local coordinates of a given patch can (ideally) be related to the global coordinates through a single rigid transform, that is, through some rotation, reflection, and translation. More precisely, with every patch $P_{i}$ we associate some (unknown) orthogonal transform $O_{i}$ and translation $t_{i}$. If point $x_{k}$ belongs to patch $P_{i}$, then its representation in $P_{i}$ is given by (cf. (1) and Figure 1)

Alternatively, if we fix a particular patch $P_{i}$, then for every point belonging to that patch,

In particular, a given point can belong to multiple patches, and will have a different representation in the coordinate system of each patch.

The premise of this paper is that we are given the membership graph and the local coordinates (referred to as measurements), namely

and the goal is to recover the coordinates $x_{1},\ldots,x_{N}$, and in the process the unknown rigid transforms $(O_{1},t_{1}),\ldots,(O_{M},t_{M})$, from (5). Note that the global coordinates are determined up to a global rotation, reflection, and translation. We say that two points clouds (also referred to as configurations) are congruent if one is obtained through a rigid transformation of the other. We will always identify two congruent configurations as being a single configuration.

Under appropriate non-degeneracy assumptions on the measurements, one task would be to specify appropriate conditions on $\Gamma$ under which the global coordinates can be uniquely determined. Intuitively, it is clear that the patches must have enough points in common for the registration problem to have an unique solution. For example, it is clear that the global coordinates cannot be uniquely recovered if $\Gamma$ is disconnected.

In practical applications, we are confronted with noisy settings where (4) holds only approximately. In such cases, we would like to determine the global coordinates and the rigid transforms such that the discrepancy in (4) is minimal. In particular, we consider the following quadratic loss:

The input to the problem are the measurements in (5). Note that our ultimate goal is to determine $x_{1},x_{2},\ldots,x_{N}$; the rigid transforms can be seen as latent variables.

The problem of multipatch registration is intrinsically non-convex since one is required to optimize over the non-convex domain of orthogonal transforms. Different ideas from the optimization literature have been deployed to attack this problem, including Lagrangian optimization and projection methods. In the Lagrangian setup, the orthogonality constraints are incorporated into the objective; in the projection method, the constraints are forced after every step of the optimization . Following the observation that the registration problem can be viewed as an optimization on the Grassmanian and Stiefel manifolds, researchers have proposed algorithms using ideas from the theory and practice of manifold optimization . A detailed review of these methods is beyond the scope of this paper, and instead we refer the interested reader to these excellent reviews . Manifold-based methods are, however, local in nature, and are not guaranteed to find the global minimizer. Moreover, it is rather difficult to certify the noise stability of such methods.

3. Contributions

The main contributions of the paper can be organized into the following categories.

Algorithm: We demonstrate how the translations can be factored out of (6), whereby the least-squares problem can be reduced to the following optimization:

Exact Recovery: We present conditions on the coefficient matrix $C$ in (7) for exact recovery using Algorithm 2. In particular, we show that the exact recovery questions about Algorithm 2 can be mapped into rigidity theoretic questions that have already been investigated earlierThe authors thank the anonymous referees for pointing this out. in . The contribution of this section is the connection made between the $C$ matrix in (7) and various notions of rigidity considered in these papers. We also present an efficient randomized rank test for $C$ than can be used to certify exact recovery (motivated by the work in ).

Stability Analysis: We study the stability of Algorithms 1 and 2 for the noise model in which the patch coordinates are perturbed using noise of bounded size (note that the stability of the spectral relaxation was not investigated in ). Our main result here is Theorem 13 which states that, if $C$ satisfies a particular rank condition, then the registration error for Algorithm 2 is within a constant factor of the noise level. To the best of our knowledge, there is no existing algorithm for multipatch registration that comes with a similar stability guarantee.

Empirical Results: We present numerical results on simulated data to numerically verify the exact recovery and noise stability properties of Algorithms 1 and 2. Our main empirical findings are the following: (1) The semidefinite relaxation performs significantly better than spectral and manifold-based optimization (say, with the spectral solution as initialization) in terms of the reconstruction quality (cf. first plot in Figure 7). (2) Up to a certain noise level, we are actually able to solve the original non-convex problem using the semidefinite relaxation (cf. second plot in Figure 7).

4. Broader Context and Related Work

The objective (6) is a straightforward extension of the objective for two-patches . In fact, this objective was earlier considered by Zhang et al. for distributed sensor localization . The present work is also closely tied to the work of Cucuringu et al. on distributed localization , where a similar objective is implicitly optimized. The common theme in these works is that some form of optimization is used to globally register the patches, once their local coordinates have been determined by some means. There is, however, some fundamental differences between the various algorithms used to actually perform the optimization. Zhang et al. use alternating least-squares to iteratively optimize over the global coordinates and the transforms, which to the best of our knowledge has no convergence guarantee. On the other hand, Cucuringu et al. first optimize over the orthogonal transforms (using synchronization ), and then solve for the translations (in effect, the global coordinates) using least-squares fitting. In this work, we combine these different ideas into a single framework. While our objective is similar to the one used in , we jointly optimize the rigid transforms and positions. In particular, the algorithms considered in Section 2 avoid the convergence issues associated with alternating least-squares in , and are able to register patch systems that cannot be registered using the approach in .

Another closely related work is the paper by Krishnan et al. on global registration , where the optimal transforms (rotations to be specific) are computed by extending the objective in (2) to the multipatch case. The subsequent mathematical formulation has strong resemblance with our formulation, and, in fact, leads to a subproblem similar to (7). Krishnan et al. propose the use of manifold optimization to solve (7), where the manifold is the product manifold of rotations. However, as mentioned earlier, manifold methods generally do not offer guarantees on convergence (to the global minimum) and stability. Moreover, the manifold in (7) is not connected. Therefore, any local method will fail to attain the global optimum of (7) if the initial guess is on the wrong component of the manifold.

It is exactly at this point that we depart from , namely, we propose to relax (7) into a tractable semidefinite program (SDP). This was motivated by a long line of work on the use of SDP relaxations for non-convex (particularly NP-hard) problems. See, for example, , and these reviews . Note that for $d=1$, (7) is a quadratic Boolean optimization, similar to the MAX-CUT problem. An SDP-based algorithm with randomized rounding for solving MAX-CUT was proposed in the seminal work of Goemans and Williamson . The semidefinite relaxation that we consider in Section 2 is motivated by this work. In connection with the present work, we note that provably stable SDP algorithms have been considered for low rank matrix completion , phase retrieval , and graph localization .

In other words, the goal is to achieve the best possible alignment of the $M$ patches through orthogonal transforms. This can be seen as an instance of the global registration problem without the translations ($t_{1}=\cdots=t_{M}=0$), and in which $\Gamma$ is complete. It is not difficult to see that (8) can be reduced to (7). On the other hand, using the analysis in Section 2, it can be shown that (6) is equivalent to (8) in this case. While the Procrustes problem is known to be NP-hard, several polynomial-time approximations with guarantees have been proposed. In particular, SDP relaxations of (8) have been considered in , and more recently in . We use the relaxation of (7) considered in for reasons to be made precise in Section 2.

5. Notations

The Kronecker product between matrices $A$ and $B$ is denoted by $A\otimes B$ . The all-ones vector is denoted by $e$ (the dimension will be obvious from the context), and $e^{N}_{i}$ denotes the all-zero vector of length $N$ with $1$ at the $i$-th position.

6. Organization

In the next section, we present the semidefinite relaxation of the least-squares registration problem described in the introduction. For reference, we also present the closely related spectral relaxation that was already considered in . Exact recovery questions are addressed in section 3, followed by a randomized test in section 4. Stability analysis for the spectral and semidefinite relaxations are presented in section 5. Numerical simulations can be found in section 6, and a discussion of certain open questions in section 7.

Spectral and Semidefinite Relaxations

The minimization of (6) involves unconstrained variables (global coordinates and patch translations) and constrained variables (the orthogonal transformations). We first solve for the unconstrained variables in terms of the unknown orthogonal transformations, representing the former as linear combinations of the latter. This reduces (6) to a quadratic optimization problem over the orthogonal transforms of the form (7).

In particular, we combine the global coordinates and the translations into a single matrix:

Similarly, we combine the orthogonal transforms into a single matrix,

To express (6) in terms of $Z$ and $O$, we write $x_{k}-t_{i}=Ze_{ki}$, where

Similarly, we write $O_{i}=O(e^{M}_{i}\otimes I_{d})$. This gives us

The matrix $L$ is the combinatorial graph Laplacian of $\Gamma$ , and is of size $(N+M)\times(N+M)$. The matrix $B$ is of size $Md\times(N+M)$, and the size of the block diagonal matrix $D$ is $Md\times Md$.

The Hessian of $\psi(Z)$ equals $2L$, and it is clear from (12) that $L\succeq 0$. Therefore, $Z^{\star}$ is a minimizer of $\psi(Z)$.

If $\Gamma$ is connected, then $e$ is the only vector in the null space of $L$ . Let $L^{\dagger}$ be the Moore-Penrose pseudo-inverse of $L$, which is again positive semidefinite. It can be verified that

If we right multiply (13) by $L^{\dagger}$, we get

Note that (16) has the global translation $t$ taken out. This is not a surprise since $\phi$ is invariant to global translations. Moreover, note that we have not forced the orthogonal constraints on $O$ as yet. Since $\phi(Z,O)\geq 0$ for any $Z$ and $O$, it necessarily follows from (16) that $C\succeq 0$. We will see in the sequel how the spectrum of $C$ dictates the performance of the convex relaxation of (16).

In analogy with the notion of stress in rigidity theory , we can consider (6) as a sum of the “stress” between pairs of patches when we try to register them using rigid transforms. In particular, the $(i,j)$-th term in (16) can be regarded as the stress between the (centered) $i$-th and $j$-th patches generated by the orthogonal transforms. Keeping this analogy in mind, we will henceforth refer to $C$ as the patch-stress matrix.

2. Optimization over orthogonal transforms

3. Spectral relaxation and rounding

This is precisely a spectral problem in that the global minimizers are determined from the spectral decomposition of $C$. More precisely, let $\mu_{1}\leq\ldots\leq\mu_{Md}$ be eigenvalues of $C$, and let $r_{1},\ldots,r_{Md}$ be the corresponding eigenvectors. Define

As noted earlier, this has a closed-form solution, namely $O_{i}^{\star}=UV^{T},$ where $U\Sigma V^{T}$ is the SVD of $W^{\star}_{i}$. We now put the rounded blocks back into place and define

In the final step, following (15), we define

The first $N$ columns of $Z^{\star}$ are taken to be the reconstructed global coordinates.

We will refer to this spectral method as the “Global Registration over Euclidean Transforms using Spectral Relaxation” (GRET-SPEC). The main steps of GRET-SPEC are summarized in Algorithm 1. We note that a similar spectral algorithm was proposed for angular synchronization by Bandeira et al. , and by Krishnan et al. for initializing the manifold optimization.

4. Semidefinite relaxation and rounding

Introducing the variable $G=O^{T}O$, (22) is equivalent to

This is a standard semidefinite program which can be solved using software packages such as SDPT3 and CVX . We provide details about SDP solvers and their computational complexity later in Section 2.5.

We now proceed as in the GRET-SPEC, namely, we define $O^{\star}$ and $Z^{\star}$ from $W^{\star}$ as in (20) and (21). We refer to the complete algorithm as “Global Registration over Euclidean Transforms using SDP” (GRET-SDP). The main steps of GRET-SDP are summarized in Algorithm 2.

Similar to Observation 1, we note the following for GRET-SDP.

where $Q=BL^{\dagger}B^{T}\succeq 0$. Bandeira et al. show that the orthogonal transforms (which we continue to denote by $O^{\star}$) obtained by a certain random rounding of $G^{\star}$ satisfy

5. Computational complexity

The main computations in GRET-SPEC are the Laplacian inversion, the eigenvector computation, and the orthogonal rounding. The cost of inverting $L$ when $\Gamma$ is dense is $O((N+M)^{3})$. However, for most practical applications, we expect $\Gamma$ to be sparse since every point would typically be contained in a small number of patches. In this case, it is known that the linear system $Lx=b$ can be solved in time almost linear in the number of edges in $\Gamma$ . Applied to (14), this means that we can compute $L^{\dagger}$ in $O(|E(\Gamma)|)$ time (up to logarithmic factors). Note that, even if $L$ is dense, it is still possible to speed up the inversion (say, compared to a direct Gaussian elimination) using the formula :

The speed up in this case is however in terms of the absolute run time. The overall complexity is still $O((N+M)^{3})$, but with smaller constants. We note that it is also possible to speed up the inversion by exploiting the bipartite nature of $\Gamma$ , although we have not used this in our implementation.

The complexity of the eigenvector computation is $O(M^{3}d^{3})$, while that of the orthogonal rounding is $O(Md^{3})$. The total complexity of GRET-SPEC, say, using a linear-time Laplacian inversion, is (up to logarithmic factors)

Exact Recovery

To express exact recovery in the matrix notation introduced earlier, define

Henceforth, we will always assume that $\Gamma$ is connected (clearly one cannot have exact recovery otherwise).

At this point, we note that if a patch has less than $d+1$ points, then even when $\bar{x}_{1},\ldots,\bar{x}_{N}$ are the unique set of coordinates that satisfy 25, we cannot guarantee $\bar{O}_{1},\ldots,\bar{O}_{M}$ and $\bar{t}_{1},\ldots,\bar{t}_{M}$ to be unique. Therefore, we will work under the mild assumption that each patch has at least $d+1$ non-degenerate points, so that the patch transforms are uniquely determined from the global coordinates.

The patch framework $\Theta=(\Gamma,(x_{k,i}))$ is affinely rigid if $y_{1},\ldots,y_{N}$ is identical to $\bar{x}_{1},\ldots\bar{x}_{N}$ up to a global affine transform.

Since each patch has $d+1$ points, we now give a characterization of affine rigidity that will be useful later on.

A patch framework $\Theta=(\Gamma,(x_{k,i}))$ is affinely rigid if and only if the rank of $C$ is $(M-1)d$.

The corollary gives an easy way to check for affine rigidity. However, it is not clear what construction of $\Gamma$ will ensure such property. In , the notion of graph lateration was introduced that guarantees affine rigidity. Namely, $\Gamma$ is said to be a graph lateration (or simply laterated) if there exists an reordering of the patch indices such that, for every $i\geq 2$, $P_{i}$ and $P_{1}\cup\cdots\cup P_{i-1}$ have at least $d+1$ non-degenerate nodes in common. An example of a graph lateration is shown in Figure 2.

If $\Gamma$ is laterated and the local coordinates are non-degenerate then the framework $\Theta$ is affinely rigid.

Following exactly the same arguments used to establish Proposition 5, one can derive the following.

The question then is under what conditions is the patch framework universally rigid? This was also addressed in using a graph construction derived from $\Gamma$ called the body graph. This is given by $\Gamma_{B}=(V_{B},E_{B})$, where $V_{B}=\{1,2,\ldots,N\}$ and $(k,l)\in E_{B}$ if and only if $x_{k}$ and $x_{l}$ belong to the same patch (cf. Figure 3). Next, the following distances are associated with $\Gamma_{B}$:

Finally, we note that universal rigidity is a weaker condition on $\Gamma$ than affine rigidity.

If a patch framework is affinely rigid, then it is universally rigid.

In , it was also shown that the reverse implication is not true using an counter-example for which the patch framework fails to be affinely rigid, but for which the body graph (a Cauchy polygon) has an unique realization in any dimension . This means that GRET-SDP can solve a bigger class of problems than GRET-SPEC, which is perhaps not surprising.

Randomized Rank Test

Corollary 6 tells us by checking the rank of the patch stress matrix $C$, we can tell whether a patch framework is affinely rigid. In this regard, the patch-stress matrix serves the same purpose as the so-called alignment matrix in and the affinity matrix in . The only difference is that the kernel of $C$ represents the degree of freedom of the affine transform, whereas kernel of alignment or affinity matrix directly tell us the degree of freedom of the point coordinates. As suggested in , an efficient randomized test for affine rigidity using the concept of affinity matrix can be easily derived. In this section, we describe a randomized test based on patch stress matrix, which parallels the proposal in . This procedure is also similar in spirit to the randomized tests for generic local rigidity by Hendrickson , for generic global rigidity by Gortler et al. , and for matrix completion by Singer and Cucuringu .

Let us continue to denote the patch-stress matrix obtained from $\Gamma$ and the measurements (25) by $C$. We will use $C_{0}$ to denote the patch-stress matrix obtained from the same graph $\Gamma$, but using the (unknown) original coordinates as measurements, namely,

The advantage of working with $C_{0}$ over $C$ is that the former can be computed using just the global coordinates, while the latter requires the knowledge of the global coordinates as well as the clean transforms. In particular, this only requires us to simulate the global coordinates. Since the coordinates of points in a given patch are determined up to a rigid transform, we claim the following (cf. Section 8.1 for a proof).

For a fixed $\Gamma$, $C$ and $C_{0}$ have the same rank.

In other words, the rank of $C_{0}$ can be used to certify exact recovery. The proposed test is based on Proposition 8.1, and the fact that if two different generic configurations are used as input in (32) (for the same $\Gamma$), then the patch-stress matrices they produce would have the same rank. By generic, we mean that the coordinates of the configuration do not satisfy any non-trivial algebraic equation with rational coefficients .

The complete test called “GRET-Randomized Rank Test” (GRET-RRT) is described in Algorithm 3. Note that the main computations in GRET-RRT are the Laplacian inversion (which is also required for the registration algorithm) and the rank computation.

Stability Analysis

We have so far studied the problem of exact recovery from noiseless measurements. In practice, however, the measurements are invariably noisy. This brings us to the question of stability, namely how stable are GRET-SPEC and GRET-SDP to perturbations in the measurements? Numerical results (to be presented in the next Section) show that both the spectral and semidefinite relaxations are quite stable to perturbations. In particular, the reconstruction error degrades quite gracefully with the increase in noise (reconstruction error is the gap between the outputs with clean and noisy measurements). In this Section, we try to quantify these empirical observations. In particular, we prove that, for a specific noise model, the reconstruction error grows at most linearly with the level of noise for the semidefinite relaxation.

The noise model we consider is the “bounded” noise model. Namely, we assume that the measurements are obtained through bounded perturbations of the clean measurements in (25). More precisely, we suppose that we have a membership graph $\Gamma$, and that the observed local coordinates are of the form

In other words, every coordinate measurement is offset within a ball of radius $\varepsilon$ around the clean measurements. Here, $\varepsilon$ is a measure of the noise level per measurement. In particular, $\varepsilon=0$ corresponds to the case where we have the clean measurements (25).

Since the coordinates of points in a given patch are determined up to a rigid transform, it is clear that the above problem is equivalent to the one where the measurements are

By equivalent, we mean that the reconstruction errors obtained using either (33) or (34) are equal. The reason we use the latter measurements is that the analysis in this case is much more simple.

The reconstruction error is defined as follows. Generally, let $Z^{\star}$ be the output of Algorithms 1 and 2 using (34) as input, and let

where we assume that the centroid of $\{\bar{x}_{1},\cdots,\bar{x}_{N}\}$ is at the origin.

Ideally, we would require that $Z^{\star}=Z_{0}$ (up to a rigid transformation) when there is no noise, that is, when $\varepsilon=0$. This is the exact recovery phenomena that we considered earlier. In general, the gap between $Z_{0}$ and $Z^{\star}$ is a measure of the reconstruction quality. Therefore, we define the reconstruction error to be

Note that we are not required to factor out the translation since $Z_{0}$ is centered by construction.

Here $\lambda_{2}(L)$ is the second smallest eigenvalue of $L$.

Under the conditions of Theorem 12, we have the following for GRET-SDP$:$

The rest of this Section is devoted to the proofs of Theorem 12 and 13. First, we introduce some notations.

Note that every $d\times d$ block of $G_{0}$ is $I_{d}$, and that we can write

We first present an estimate that applies generally to both algorithms. The proof is provided in Section 8.2.

Let $R$ be the radius of the smallest Euclidean ball that encloses the clean configuration. Then, for any arbitrary $\Theta$,

In other words, the reconstruction error in either case is controlled by the rounding error:

The rest of this Section is devoted to obtaining a bound on $\delta$ for GRET-SPEC and GRET-SDP. In particular, we will show that $\delta$ is of the order of $\varepsilon$ in either case. Note that the key difference between the two algorithms arises from the eigenvector rounding, namely the assignment of the “unrounded” orthogonal transform $W^{\star}$ (respectively from the patch-stress matrix and the optimal Gram matrix). The analysis in going from $W^{\star}$ to the rounded orthogonal transform $O^{\star}$, and subsequently to $Z^{\star}$, is however common to both algorithms.

We now bound the error in (39) for both algorithms. Note that we can generally write

where $u_{1},\ldots,u_{d}$ are orthonormal. In GRET-SPEC, each $\alpha_{i}=M$, while in GRET-SDP we set $\alpha_{i}$ using the eigenvalues of $G^{\star}$.

Our first result gives a control on the quantities obtained using eigenvector rounding in terms of their Gram matrices.

Next, we use a result by Li to get a bound on the error after orthogonal rounding.

The proofs of Lemma 15 and 16 are provided in Appendices 8.3 and 8.4. At this point, we record a result from which is repeatedly used in the proof of these lemmas and elsewhere.

In particular, the above result holds for the Frobenius and spectral norms.

By combining Lemma 15 and 16, we have the following bound for (39):

We now bound the quantity on the right in (40) for GRET-SPEC and GRET-SDP.

For the spectral relaxation, this can be done using the Davis-Kahan theorem . Note that from (18), we can write

Following [7, Ch. 7], let $A$ be some symmetric matrix and $S$ be some subset of the real line. Denote $P_{A}(S)$ to be the orthogonal projection onto the subspace spanned by the eigenvectors of $A$ whose eigenvalues are in $S$. A particular implication of the Davis-Kahan theorem is that

where $S_{1}^{c}$ is the complement of $S_{1}$, and $\rho(S_{1},S_{2})=\min\{|u-v|:u\in S_{1},v\in S_{2}\}$.

Now, it is not difficult to verify that for the noise model (34),

Combining Proposition 14 with (40),(43), and (44), we arrive at Theorem 12.

2. Bound for GRET-SDP

We record a result about this decomposition from Wang and Singer .

Suppose $G_{0}+\Delta\succeq 0$ and $\Delta_{ii}=0\ (1\leq i\leq M)$. Let $\Delta=P+Q+T$ as in (45). Then

and $O^{\star\star}_{i}$ to be the $i$-th $Md\times d$ block of $O^{\star\star}$, that is, $O^{\star\star}\stackrel{{\scriptstyle\text{def}}}{{=}}[O^{\star\star}_{1}\ \cdots\ O^{\star\star}_{M}]$.

By construction, $G^{\star}={O^{\star\star}}^{T}O^{\star\star}$. Moreover, by feasibility,

In particular, we will use the fact that $(Z^{\star\star},O^{\star\star})$ are the minimizers of the unconstrained program

We are done if we can bound the term on the right. To do so, we note from (47) that

and use $\lVert x+y\rVert^{2}\leq 2(\lVert x\rVert^{2}+\lVert y\rVert^{2})$ to get

Finally, using the optimality of $(Z^{\star\star},O^{\star\star})$ for (47), we have

The desired result follows from (48), (49), and (50). ∎

We will heavily use decomposition (45) and its properties. Let $G^{\star}=G_{0}+\Delta$. By triangle inequality,

Moreover, since the bottom eigenvalues of $G_{0}$ are zero, it follows from Lemma 17 that the norm of the diagonal matrix is bounded by $\|\Delta\|_{F}$. Therefore,

That is, $(A(p,q))$ are the coordinates of $A$ in the basis $\{s_{1},...,s_{d}\}\cup\ \{s_{d+1},\ldots,s_{Md}\}$.

Decompose $\Delta=P+Q+T$ as in (45). Note that $P,Q,$ and $T$ are represented in the above basis as follows: $P$ is supported on the upper $d\times d$ diagonal block, $T$ is supported on the lower $(M-1)d\times(M-1)d$ diagonal block, and $Q$ on the off-diagonal blocks. The matrix $G_{0}$ is diagonal in this representation.

where we have used the properties $T\succeq 0$ and $T_{ll}\succeq 0\ (1\leq l\leq M)$. In particular,

On the other hand, since $G_{0}+\Delta\succeq 0$, we have $(G_{0}+\Delta)(p,q)^{2}\leq(G_{0}+\Delta)(p,p)(G_{0}+\Delta)(q,q)$. Therefore,

Combining (51), (52), (54), and (53), we get the desired bound. ∎

Putting together (40) with Propositions (14),(19), and (20), we arrive at Theorem (13).

Numerical Experiments

We now present some numerical results on multipatch registration using GRET-SPEC and GRET-SDP. In particular, we study the exact recovery and stability properties of the algorithm. We define the reconstruction error in terms of the root-mean-square deviation (RMSD) given by

In other words, the RMSD is calculated after registering (aligning) the original and the reconstructed configurations. We use the SVD-based algorithm for this purpose.

In the left plot in Figure 4, we consider a patch system with $N=10$ points. The points that belong to two or more patches are marked red, while the rest are marked black. The patches taken in the order $P_{1},P_{2},P_{3}$ form a lateration in this case. As predicted by Corollary 6 and Theorem 7, the rank of the patch-stress matrix $C_{0}$ for this system must be $2(3-1)=4$. This is indeed confirmed by our experiment. We expect GRET-SPEC and GRET-SDP to recover the exact configuration. Indeed, we get a very small RMSD of the order of 1e-7 in this case. As shown in the figure, the reconstructed coordinates obtained using GRET-SDP perfectly match the original ones after alignment.

We next consider the example shown in the center plot in Figure 4. The patch system is not laterated in this case, but the rank of $C_{0}$ is $4$. Again we obtain a very small RMSD of the order 1e-7 for this example. This example demonstrates that lateration is not necessary for exact recovery.

Finally, we plot the rank of the SDP solution $G^{\star}$ and notice an interesting phenomenon. Up to a certain noise level, $G^{\star}$ has the desired rank and rounding is not required. This means that the relaxation gap is zero for the semidefinite relaxation, and that we can solve the original non-convex problem using GRET-SDP up to a certain noise threshold. It is therefore not surprising that the RMSD shows no improvement after we refine the SDP solution using manifold optimization. We have noticed that the rank of the SDP solution is stable with respect to noise for other numerical experiments as well (not reported here).

Discussion

There are several directions along which the present work could be extended and refined. We summarize some of these below.

Conditions on $\Gamma$. We have seen that the universal rigidity of the body graph (derived from $\Gamma$) is both necessary and sufficient for exact recovery using GRET-SDP. However, to test unique rigidity, we need to run a semidefinite program . Unfortunately, the complexity of this program is much more than GRET-SDP itself. This led us to consider the rank criteria that could be tested efficiently. The rank test is nonetheless not necessary for exact recovery, and weaker conditions can be found. In particular, an interesting question is whether we could find an efficiently-testable condition that would hold true for the extreme example in Figure 4, in which $\Gamma$ fails the rank test?

Tighter Bounds. The stability in Theorem 13 was for the bounded noise model, which made the subsequent analysis quite straightforward. The goal was to establish that the reconstruction error is within $C\varepsilon$ for some constant $C$ independent of the noise. In particular, the bounds in Theorem 13 are quite loose. One possible direction would be to consider a stochastic noise model with statistically independent perturbations to tighten the bound.

Anchor Points. In sensor network localization, one has to infer the coordinates of sensors from the knowledge of distances between sensors and its geometric neighbors. In distributed approaches to sensor localization , one is faced exactly with the multipatch registration problem described in this paper. Besides the distance information, one often has the added knowledge of the precise positions of selected sensors known as anchors . This is often by design and is used to improve the localization accuracy. The question is can we incorporate the anchor constraints into the present registration algorithm? One possible way of leveraging the existing framework is to introduce an additional patch (called anchor patch) for the anchor points. The anchor coordinates are assigned to the points in the anchor patch (treating them as local coordinates). This gives us an augmented bipartite graph $\Gamma_{a}$ which has one more patch vertex than $\Gamma$, and extra edges connecting the anchor patch to the anchor vertices. We then proceed exactly as before, that is, we solve for the global coordinates of both the anchor and non-anchor points given the measurements on $\Gamma_{a}$.

Acknowledgements

K. N. Chaudhury was partially supported by the Swiss National Science Foundation under Grant PBELP2-135867 and by Award Number R01GM090200 from the National Institute of General Medical Sciences while he was at Princeton University, and more recently by a Startup Grant from the Indian Institute of Science.

Y. Khoo was partially supported by Award Number R01 GM090200-01 from the National Institute for Health.

A. Singer was partially supported by Award Numbers FA9550-13-1-0076 and FA9550-12-1-0317 from AFOSR, Award Number R01GM090200 from the National Institute of General Medical Sciences, and Award Number LTR DTD 06-05-2012 from the Simons Foundation.

The authors are grateful to the anonymous referees for their comments and suggestions. They particularly thank the referees for pointing out references and their relation to the present work. The authors also thank Mihai Cucuringu, Lanhui Wang, and Afonso Bandeira for useful discussions, and Nicolas Boumal for advice on the usage of the Manopt toolbox.

References

Technical proofs

In this Section, we give the proof of Propositions 11 and 14, and Lemmas 15 and 16,.

We are done if we can show that there exists a bijection between the nullspace of $C$ and that of $C_{0}$. To do so, we note that the associated quadratic forms can be expressed as

Now, it follows from (25) that there is a one-to-one correspondence between $u$ and $v$, namely

such that $u^{T}Cu=v^{T}C_{0}v$. In other words, the null space of $C$ is related to the null space of $C_{0}$ through an orthogonal transform, as was required to be shown.

2. Proof of Proposition 14

Without loss of generality, we assume that the smallest Euclidean ball that encloses the clean configuration $\{\bar{x}_{1},\ldots,\bar{x}_{N}\}$ is centered at the origin, that is,

Let $B_{0}$ be the matrix $B$ in (12) computed from the clean measurements, i.e., from (34) with $\varepsilon=0$. Let $B_{0}+H$ be the same matrix obtained from (34) for some $\varepsilon>0$.

Recall that $Z_{0}=O_{0}B_{0}L^{\dagger}$ (by the centering assumption in (35)). Therefore,

where $\lambda_{2}(L)$ is the smallest non-zero eigenvalue of $L$. On the other hand,

As for the other term in (57), we can write

Therefore, using Cauchy-Schwarz, the orthonormality of the columns of $O^{\star}_{i}$’s, and the noise model (34), we get

Combining (57),(58), and (59), we get the desired estimate.

3. Proof of Lemma 15

The proof is mainly based on the observation that if $u$ and $v$ are unit vectors and $0\leq u^{T}v\leq 1$, then

$(\Omega_{1}s_{i})^{T}(\Omega_{2}u_{j})=0$ for $i\neq j$, and $0\leq(\Omega_{1}s_{i})^{T}(\Omega_{2}u_{i})\leq 1$ for $1\leq i\leq d$.

Here $\Theta_{1}$ and $\Theta_{2}$ are the orthogonal transforms that map $\{u_{1},\ldots,u_{d}\}$ and $\{s_{1},\ldots,s_{d}\}$ into the corresponding principal vectors.

Now, using (60) and the principal angle property $3$, we get

Moreover, using triangle inequality and properties $1$ and $2$, we have

Combining the above relations, and setting $\Theta=\Theta_{1}^{T}\Theta_{2}$, we arrive at Lemma 15.

4. Proof of Lemma 16

This is done by adapting the following result by Li : If $A,B$ are square and non-singular, and if $\mathcal{R}(A)$ and $\mathcal{R}(B)$ are their orthogonal rounding (obtained from their polar decompositions ), then

We recall that if $A=U\Sigma V^{T}$ is the SVD of $A$, then $\mathcal{R}(A)=UV^{T}$.

Note that it is possible that some of the blocks of $W^{\star}$ are singular, for which the above result does not hold. However, the number of such blocks can be controlled by the global error. More precisely, let $\mathcal{B}\subset\{1,2,\ldots,M\}$ be the index set such that, for $i\in\mathcal{B}$, $\|W^{\star}_{i}-\Theta\|_{F}\geq\beta$. Then

This gives a bound on the size of $\mathcal{B}$. In particular, the rounding error for this set can trivially be bounded as

On the other hand, we known that, for $i\in\mathcal{B}^{c}$, $\|W^{\star}_{i}-\Theta\|_{F}<\beta$. From Lemma 17, it follows that

Fix $\beta\leq 1$. Then $\sigma_{\min}(W^{\star}_{i})>1-\beta$, and we have from (62),

Fixing $\beta=1/\sqrt{2}$ and combining (63) and (64), we get the desired bound.