An as-rigid-as-possible approach to sensor network localization

We present a novel approach to localization of sensors in a network given a subset of noisy inter-sensor distances. The algorithm is based on “stitching” together local structures by solving an optimization problem requiring the structures to fit together in an “As-Rigid-As-Possible” manner, hence the name ARAP. The local structures consist of reference “patches” and reference triangles, both obtained from inter-sensor distances. We elaborate on the relationship between the ARAP algorithm and other state-of-the-art algorithms, and provide experimental results demonstrating that ARAP is significantly less sensitive to sparse connectivity and measurement noise. We also show how ARAP may be distributed.


INTRODUCTION
Sensor networks are a distributed collection of small devices, capable of a limited amount of local processing and wireless communication [19].They are drawing more and more attention for use in a variety of applications, such as environment monitoring, military surveillance, smart places etc [12].In many typical applications, sensors are deployed over a physical area, sometimes at random, without any prior knowledge of their locations.Lacking GPS components, a sensor cannot know its final location in the arena, which may be necessary for it to be useful.Thus autonomous computation of the spatial coordinates of the sensors, based on just the limited information available to the sensors, is an important practical problem, known as the localization problem.
One useful source of information for the localization problem is inter-sensor distances.Two sensors with a communication link between them may be able to measure the distance between them and use it for localization.This type of information is the typical input to a localization algorithm.The communication links between the sensors may be modeled as a sensor graph, and the distances measured along the edges of this graph are called the measurement graph.The coordinates of the sensors computed during localization are sometimes called an embedding, realization or layout of the measurement graph.
More formally, this paper treats the following problem: Given a set of sensors with unknown spatial distribution, and a mechanism to measure the distances between a sensor and its nearby (neighbor) sensors, determine the coordinates of all the sensors through local sensor-to-sensor communication.
The localization problem is NP-Hard [16,20], i.e. there is no known efficient algorithm that is guaranteed to find the correct embedding (or even an approximate solution) for all measurement graphs.This is true even when the measurement graph is generically globally rigid, i.e., when the coordinates are completely determined, up to a rigid transformation, from the measured distance data.In this case the embedding is called globally rigid.Note that generic global rigidity (assuming that the embedding is also "generic") is a property of the graph structure, not of the actual distance measurements [8].Only if the measurement graph is known to belong to a special subclass of globally rigid graphs, such as trilateration graphs [4], efficient localization may be provably possible.
Despite the inherent difficulty in the general case, the importance of this problem has led to many heuristic and numerical algorithms for its approximate solution.The realistic, and most challenging, scenario is when the measurement graph is sparse and noisy.This means that there is only small amount of unreliable information available.Many existing localization algorithms are extremely sensitive to this and degrade very quickly when presented with these types of inputs.This paper presents an algorithm which is relatively easy to implement, can be distributed among the sensors, and, most important, is extremely robust, even when the measurement graph is sparse and noisy.

II. RELATED WORK
Assume the sensor network is modeled by a graph G = (V={1,..,n}, E), such that for every edge (i,j)∈E we know the (noisy) distance d ij between sensors i and j.We assume, without loss of generality, that the distances are always symmetric, so d ij =d ji .This can be achieved by requiring neighboring sensors to reach an agreement about their distance (e.g. by averaging d ij and d ji ).We would like to generate a two-dimensional embedding P = (p 1 ,..,p n ) that realizes all known distances.This embedding can be formally characterized as the minimum of the stress function: Lei Zhang Ligang Liu Zhejiang University Hangzhou, China {ethan_zhang, ligangliu}@zju.edu.cn The stress function and its variants are widely used in problems in distance geometry and in statistics, where distances measure some sort of affinity between subjects.When all pairwise distances are available, it is possible to easily minimize the strain function -a function very similar to stress.Strain is a convex function, thus may be minimized using Classical Multi-Dimensional Scaling (MDS) [2], which, in practice, is an eigenvector computation.The situation changes dramatically when only a subset of the pairwise distances is available.Stress is a non-convex function, hence, in general, we cannot find its global minimum efficiently, and slow iterative methods must be employed.The main phenomenon observed in the (suboptimal) solutions is that of foldovers, where entire pieces of the embedding fold over on top of others, while still realizing well most of the measurement graph.
In the general case, stress may be minimized using iterative stress majorization, also called the SMACOF algorithm [2].Majorization is a modern optimization technique achieving better results than applying traditional optimization methods, such as gradient descent, to the stress function.Although all techniques will eventually get stuck in local stress minima, and require a good foldover-free initial embedding to start from [9], stress majorization is faster and seems to "jump over" many local minima to a generally better solution.
In an attempt to depart from the cumbersome (and slow) global optimization strategies, Moore et al. [13] describe a localization algorithm which takes advantage of the special nature of the connectivity of sensor network measurement graphs.Since each sensor typically can communicate only with its local neighborhood, the connectivity tends to be local, with only short edges.These graphs can even be mathematically modeled as "disk graphs", where two sensors are connected if (and only if) the distance between them is less than some parameter.This characteristic of sensor network connectivity implies that the connectivity within small neighborhoods is typically dense, even frequently making those neighborhoods generically globally rigid.Thus Moore et al. [13] propose an incremental algorithm which attempts to first localize only small "patches" of the network.Each such patch is a set of four sensors forming a rigid quadrilateral, which is localized using a procedure called "trilateration" and then improved by running stress minimization on that patch.Given one quad which has been localized, another quad is found which has some sensors in common with the first.This is then laid out relative to the first by applying the best possible rigid transformation between the two.In such a manner, a sequence of quads is laid out in breadth-first order until no more sensors can be localized.
While a very promising approach, the algorithm of Moore et al. will localize only those sensors contained in a trilateralizable component of the network.The algorithm will not produce anything useful for the other sensors.More unfortunately, there exist generically globally rigid graphs that are not trilateralizable.Goldenberg et al. [6] describe a method to find generically globally rigid subgraphs (on which they run an MDS approach), and thus avoid contaminating their answer with underdetermined portions of the network.
A similar "patching" algorithm was described by Shang and Ruml [17].However, they localize an individual patch by first computing all pairwise shortest paths (in the measurement graph) between sensors in the patch.Classical MDS is then applied to these distances to yield an initial embedding, which is subsequently improved by stress minimization.Similarly to Moore et al. [13], the patches are "stitched" together incrementally in a greedy order by finding the best affine transformation between a new patch and the global embedding.A postprocessing stage of the algorithm further improves the embedding by minimizing the stress energy of the complete graph.These incremental algorithms seem to be sound, but, unfortunately, like any other incremental algorithm, may accumulate error indefinitely, especially when the measurement graph is very noisy.
The PATCHWORK algorithm of Koren et al. [11] solves the error accumulation error problem which may affect the methods of Moore et al. [13] and Shang and Ruml [17].It also first localizes small patches of sensors.However, instead of then gluing them together incrementally, and possibly accumulating error in the process, it applies a global process which attempts to map the patches to a global coordinate system via per-patch affine transformations.Since the patches overlap, this involves solving a linear least-squares problem in order that the transformations agree well on the locations of sensors they have in common.Koren et al. [11] observe that using affine transformations between each patch and the global coordinate system, may be interpreted as affine relationships between patches.Some algebraic manipulation shows that PATCHWORK reduces to a method very similar to dimension reduction methods [3,15], used in machine learning.
A more recent approach of Singer [18], called Locally-Rigid Embedding (LRE), makes an explicit connection between localization and the machine learning technique of Locally-Linear Embedding (LLE).LRE tries to preserve the local affine relationships, present within patches, in the global coordinate system.Because of the overlap between patches, this imposes affine relationships between patches.A linear system is set up, each sensor contributing one equation relating its location to those of its neighbors.Solving this system results in an embedding of all sensors, from which a global affine transformation must be removed.
The method we propose here, called As-Rigid-As-Possible (ARAP) belongs to the same family as PATCHWORK and LRE.We also start off by localizing small patches in a very similar manner.However, instead of then embedding them in a global coordinate system using affine mappings, we embed them using rigid mappings.Here too, the overlap between patches constrains the mappings.Using rigid mappings has the advantage of preserving better the local relationships between patches, but has the disadvantage of resulting in a (sparse) nonlinear system of equations.Fortunately, this non-linear system may be solved rather simply and efficiently using a two-phase alternating least-squares method.Along the way, we show how to improve LRE and PATCHWORK at just a slightly more expensive price of solving a larger (but still sparse) linear system, which we call As-Affine-As-Possible (AAAP).In fact, we use AAAP as the starting point for the iterations of ARAP.Experimental results show that ARAP is more robust to sparsity and noise in the measurement graph than all its predecessors.So, while ARAP might seem to be a straightforward extension of LRE, it provides a very large increase in perform-ance in return for a very small increase in computational complexity.

III. THE ARAP ALGORITHM OVERVIEW
A sensor network can be modeled as a measurement graph G = (V={1,..,N}, E), where the sensors are treated as the nodes of the graph, and edge (i,j) ∈E is associated with the distance d ij measured between sensors i and j.We denote by N(i) = {j:(i,j)∈E }∪{i} the set of sensors connected to sensor i, and their number by N i .Sensor i has no prior knowledge about its location, and knows only the distances to its neighbors in N(i).We would like to find a localization of these sensors in the plane, i.e., P = {p 1 ,.., p n }, where p i ∈ℜ 2 , such that ||p i -p j || is as close as possible to the given d ij , essentially a global minimum of the stress energy (1).
Since a triangle is generically globally rigid, any sensor i, along with its two neighbors j and k, whose accurate intersensor distances d ij , d jk , d ki are known, may be uniquely positioned (or localized) in an arbitrary coordinate system, up to a rigid transformation.Call a triplet of sensors for which this data is given an available triangle, and its localization a reference triangle.Number the available triangles by t=1,..,T, and denote their localizations in the local coordinate systems by Q t = {q t i , q t j , q t k }.More complicated graphs, such as the subgraph of G induced by N(i), are not necessarily generically globally rigid, even if G is.And even if they are generically globally rigid, they still may not be localizable using known methods.But they may be, if the subgraph contains sufficient edges.Thus we may attempt to localize the sensors in N(i) in an arbitrary coordinate system, again up to a rigid transformation.The localization of N(i) is called a reference patch, and denoted by Q i = {q i ,q i1 , q iNi }.See Figure 1 for an illustration of reference triangles and patches.
The spirit of our algorithm is to construct as many of these reference triangles and reference patches as possible, and "stitch" them together in a global coordinate system to localize the entire measurement graph.The latter is done by mapping each patch to the global coordinate system using a rigid transformation.Both phases of the algorithm try to respect the measurement graph as much as possible, thus minimizing the stress energy.

IV. REFERENCE TRIANGLE AND PATCH LOCALIZATOIN
Reference triangles can be obtained easily.If the distances between the three nodes i, j and k are given, they have a stable localization in the plane up to a rigid transformation: a translation, rotation and reflection (sometimes called "flip").This follows from elementary geometry.Localizing a patch N(i) is more complicated.The distances between i and the other sensors in N(i) are all known.Only some of the distances between the neighboring sensors are known.We use the three-stage method of PATCHWORK [11], which proceeds as follows.
Estimating missing distances.We estimate the distance d jk for all (j,k)∉E, j,k∈N(i).From the triangle inequality, we obtain an upper bound B jk on d jk , and from the disk graph assumption a lower bound b jk on d jk : Then, the estimate for d jk is ( ) Classical MDS.After estimating all the missing distances between sensors in N(i), the reference patch can be obtained by Classical MDS [2].
Stress majorization.This step is to improve the embedding of the reference patch by using only the known distances between the sensors in N(i).It involves iteratively applying the following update rule, for each i∈N(i): ( ) V. RIGID ALIGNMENT Once the reference triangles and patches are obtained, it remains to align ("stitch") all the triangles and patches together into one coherent layout.If we assume that we have properly localized the triangles and patches up to a rigid transformation, ideally the alignment should allow only rigid transformations between the patches and the global coordinate system.Reference triangles and patches undergo the same "stitch" operation, so we give a detailed derivation of the rigid alignment for reference patches, which is applicable also to the reference triangles.
Denote by the 2×N i matrix Q i = (q i , q i1 ,.., q iNi ), where q i = (q i x , q i y ) T are the sensor locations in the reference patch N(i) and by P = (p 1 ,.., p N ), a 2×N matrix, where p i = (p i x , p i y ) T is the global localization of the i'th sensor.Note that we assume that the location vectors are 2D column vectors.We would like to find a P such that all the reference patch localizations relate to it via rigid transformations, i.e. find also R i and T i for Q i , such that: P i is the localization of N(i) induced by P, R i is a 2×2 rigid transformation matrix, and T i is a translation vector.
To eliminate the translational component T i immediately, we work only with vector differences between the coordinates of the sensors in N i and the coordinates of the i'th sensor.Thus, to perform localization of the entire sensor network, we have to solve the following optimization problem simultaneously for P (from which all P i are induced), and all the R i (one for each patch): ( , ,.., ) arg min : where || ⋅|| F is the Frobenius matrix norm, and the P i and Q i are the appropriately translated vectors.Once this cost function is minimized, we may discard the R i .
We now make two key observations: 1. Given P i and Q i , (3) may be solved separately (and locally) for each R i .
{ } 2 arg min : This is the least-squares version of (2), since it will not always be possible to find a R i satisfying (2) exactly, as the system may be overconstrained.Sadly, ( 4) is a non-linear least squares problem, because of the orthogonality constraint on R i .Fortunately, it still has a relatively simple solution.Solving (4) for the best rigid matrix R i for each N(i) is an instance of the Rotation Orthogonal Procrustes Problem, which can be solved by Procrustes analysis [10].A closed form solution is , where U i and V i are the orthogonal components of the Singular Value Decomposition (SVD) [7,14] of Q i P i T : (5) If all R i are given, then the remaining unknown P can be easily obtained by solving a quadratic optimization problem: By setting the gradients to zero, we obtain the following 2N×2N linear system of normal equations in the (x and y) components of P: These two observations mean that we can solve (3) using a two-phase Alternating Least-Squares (ALS) method: In the first local phase, P is assumed to be fixed, and we solve (4) for each of the R i .In the second global phase, all the R i are assumed to be fixed, and we solve (6) for P.
We call the use of the ALS technique to minimize (3) the "As-Rigid-As-Possible" (ARAP) localization algorithm.Since the entries of the coefficient matrix of ( 7) depend only on the reference patch, this matrix is a sparse matrix which is fixed throughout ARAP.Thus we may pre-factor the matrix, and reuse the factorization to efficiently solve (7) with a different left-hand side at each iteration.The observant reader may recognize this matrix as the combinatorial Laplacian matrix of the measurement graph [1], and (7) as the celebrated Poisson equation on a graph [14].
To run ARAP, we need an initial localization which the ARAP will refine as it alternates through the two phases.We describe how to do this in the next section.

A. AAAP Localization
At this point we make another key observation.
If the measurement graph is connected and generically globally rigid and all patches have been correctly localized (up to a rigid transformation), then it is possible to correctly localize the entire graph (up to a global rigid transformation) by solving (3), while allowing R i to be a similarity, or even affine, transformation.
This observation is non-trivial.It relies on the concept of stress matrices in rigidity theory, and is implied by the main result of Gortler et al. [8] that generic global rigidity of a graph implies the existence of a symmetric equilibrium stress matrix of co-rank 3 (and no more) for the graph.This observation, surprising at first glance, becomes less surprising when we realize that most of the localization work has already been done in the reference patch localization, and the fact that the graph is generically globally rigid implies that these localizations essentially determine the global localization.Thus relaxing the conditions on R i do not really allow the stitching stage to misbehave.The advantage of relaxing the condition on R i is that now (3) becomes a single global linear system in the parameters of R i and P, so may be solved in "one shot", as opposed to the iterative ARAP.In fact, the resulting per-patch affine transformations will probably be all rigid transformations composed with a common global affine transformation.There is, however, one caveat: The global localization is unique only up to a global transformation taken from the relevant family used: similarities or affine.In particular, this means that (3), without the condition on R i , will be homogeneous, thus admit the trivial solutions P = 0 and R i = 0 for all i.To avoid this, either a small number of "anchored" sensor positions must be known in advance, and these "pinned down" when solving the system, or the linear system can be solved as an eigenvector problem, which effectively factors out the nullspace.In addition, the global transformation can be "factored out" (up to a rigid transformation) in a postprocessing stage.This is done by comparing the inter-sensor distances of the result with the input measurement graph and finding the best transformation that minimizes the distortion of these values.
To compute an initial localization, we solve (3) when R i are required to only be affine transformations.If we parameterize R i as follows: (8) ..., , ( ) which is linear (and homogeneous) in the variables.
We call this algorithm the As-Affine-As-Possible (AAAP) localization procedure.
Should we want to restrict the transformations R i to belong to the class of similarities, this would impose the additional constraints a i = d i and c i = -b i , which are still linear, as opposed to imposing that R i are rigid transformations, which would additionally imply that a i 2 +b i 2 =1, a non-linear equation.
The observant reader might ask why not solve (3) for R i restricted to the family of similarities, since the system would still be linear, more compact, and closer to the desired class of linear transformations -the rigid ones ?The answer is that, surprisingly, affine transformations are more suitable, since they also allow reflections, which similarities do not allow.These reflections are crucial for the stitching process, since we no information on the orientation of patches, and their reference localizations must be allowed to reflect ("flip") if needed.
It is relatively simple to see that our initial AAAP localization procedure is quite similar to the PATCHWORK method of Koren et al. [11], but a little more general.The interested reader is referred to Section 3.2 of [11] for details.We have already stated that in a noiseless scenario this method would suffice.However, as we shall demonstrate later, in a realistic, noisy, scenario, the method degrades rather quickly.Koren et al. [11] compensate for this by applying standard stress majorization in a post-processing step, but this is quite slow, and has limited effect.Thus AAAP suffices only as an initial embedding for ARAP.

B. A Distributed Solution
order for a localization algorithm to be practical in a real scenario, the algorithm must be distributable, meaning that it can be computed on the sensor network through exchange of information between a sensor and its immediate neighbors only.This is the case for our ARAP algorithm.Obviously, reference patch computation is local, as is the local stage (4).The global linear systems, solved during the AAAP initialization (9) and the global phase (7), are sparse.The associated matrices have non-zero values only between nodes and their immediate neighbors, thus can be solved using a fully distributed version of the well-known iterative Gauss-Seidel method [14,21].

VI. EXPERIMENTAL RESULTS AND COMPARISON
We now present the results of our ARAP approach applied to some sample inputs.As done by previous authors, we generated synthetic inputs on 200 points distributed uniformly within 2D regions of different shapes, endowed with the connectivity of a "disk graph", meaning that all pairs of sensors closer than a parameter r are connected.We used two regions: a simple square, and a non-convex C-shaped region, and three different values of r, corresponding to sparse, medium and dense graphs, whose average number of neighbors is 4, 5 and 6, respectively.The more dense the graph, the more localizable it should be.See Figures 2 and 3. ARAP converged in less than 40 iterations in all the experiments.
To measure the localization performance of the algorithms, we simply compute the average normalized error in localization per sensor: where p i is the localization of sensor i, g i is the gound truth location of the sensor (available in our simulation experiments), N is the number of sensors and D is the diameter of the embedding.Before applying this measure, the best rigid transformation between the embedding and the ground truth was computed and factored out.To test the sensitivity of the algorithm to noise, random noise ε ij in the range [-σl ij , σl ij ] is added to the true edge lengths l ij when constructing the measurement graph: We compared our ARAP approach with the traditional SMACOF approach, when initialized using a Laplacian eigenvector embedding, as in [9], the LRE approach [18], and the AAAP approach, which is essentially our initial embedding for the ARAP procedure.We use the same three-stage reference patch localization procedure for LRE, AAAP and ARAP (typically requiring 20 stress minimization iterations per patch), as described in Section IV.The linear systems involved in LRE and AAAP are solved by an eigenvector calculation, and a global affine transformation is factored out after that, as described by Singer [18].
Our results are summarized in Figs.2-4, which provide a consistent, and clear, picture of the comparative performance of the algorithms.LRE, AAAP and ARAP, who are all based on the same local reference patch localization, are capable of "reproducing" the correct embedding (up to a rigid transformation) when no noise is present and the measurement graph is dense enough.The difference between the three becomes clear immediately upon departure from this optimistic scenario.The graphs in Fig. 4 show the performance for noise levels between 0% and 10%, at densities equivalent to 4, 5, and 6 neighbors per sensor on the average.The embeddings in Figs.2-3 were generated for inputs contaminated with 5% noise, at the same densities.These are the most interesting scenarios, as this is more or less where the transition to generic global rigidity typically occurs.The runtime required to generate each embedding on a 1.86GHz PC with 1G RAM, along with the resulting localization error, are provided for each embedding.
Both AAAP and LRE, who allow affine transformations between patches, cannot overcome the inaccurate localizations of the patches, and even amplify the noise.AAAP is typically a little better than LRE.ARAP, on the other hand, compensates for inaccurate patch localization with rigid transformations between the patches, and is able to filter out the noise.It results in localizations which are up to an order of magnitude more accurate than the others, at the price of no more than a factor of two in runtime.Remember, also, that our versions of LRE and AAAP require an eigenvector computation and global affine transformation removal, which are not easily distributable.From our experience, distributable algorithms for achieving a similar effect, such as pinning down a number of anchor vertices and solving a simple linear system (which may be distributed by a iterative Gauss-Seidel procedure) typically produce less accurate results, so the results presented here for LRE and AAAP may be interpreted as an upper bound on the performance of any algorithm within that family.SMACOF, which starts off with a reasonable, but nowhere near accurate, embedding, tries to directly minimize stress relative to the measurement graph, but quickly gets stuck in one of many local minima, typically involving "foldovers" in the embedding.

VII. DISCUSSION
Our new ARAP approach follows the paradigm of "think globally, fit locally".As only information about local distances is available, we are able to localize small reference patches and triangles, and then localize the entire network by fitting the patches together with rigid transformations to a global coordinate system.
It is interesting to explore the relationships between the various methods used in this paper.The iterative SMACOF tries to directly minimize the stress cost function (1).Close inspection of the SMACOF algorithm (see the excellent derivation by Gansner et al. [5]) shows that it is actually quite similar in spirit to our iterative ARAP algorithm, and can be considered an edge-based version of ARAP.By edge-based, we mean that the algorithm tries to stitch together edges rather than triangles or patches.But apart from its length, an edge does not contain much information on the local rigid structure of the graph, so is a very weak version of ARAP.Since the cost function ( 1) is known to be plagued with many local minima, SMACOF will quickly get stuck in a suboptimal solution.
AAAP and LRE are extremely similar.Rather than minimize a global stress function, as SMACOF does, AAAP tries to minimize the cost function (3), without constraints on R i , which is therefore a convex quadratic function, having a unique global minimum, obtained by solving a linear system of equations.The main difference between AAAP and LRE are whether the optimal affine transformation between a patch and the global coordinate system is solved for explicitly along with the global embedding, or solved for separately and hard-wired into the linear system for the global embedding.Koren et al. [11] also distinguish between versions of the PATCHWORK algorithm generating N i linear equations per patch, which is the spirit of AAAP, and just one linear equation per patch, which is the spirit of LRE.
ARAP also minimizes cost function (3), with the non-linear constraint on R i .This cost function may have local minima, which ARAP will get stuck in, but these seem to be much fewer than those of the stress cost function (1), and provide more accurate localizations than the global minimum corresponding to the cost function of AAAP or LRE.

Figure 1 .
Figure 1.Possible reference patch and triangles associated with a given vertex (red dot) in a measurement graph.Solid lines are edges in the measurement graph (i.e. the length of the edge is known).Dashed lines are edges not in the measurement graph (i.e. the length of the edge is unknown).

0 i
In this results in a quadratic optimization problem with unknowns p i = (p i x , p i y ) T , a i , b i , c i and d i , whose normal equations are: