Matrix Completion Optimization for Localization in Wireless Sensor Networks for Intelligent IoT

Localization in wireless sensor networks (WSNs) is one of the primary functions of the intelligent Internet of Things (IoT) that offers automatically discoverable services, while the localization accuracy is a key issue to evaluate the quality of those services. In this paper, we develop a framework to solve the Euclidean distance matrix completion problem, which is an important technical problem for distance-based localization in WSNs. The sensor network localization problem is described as a low-rank dimensional Euclidean distance completion problem with known nodes. The task is to find the sensor locations through recovery of missing entries of a squared distance matrix when the dimension of the data is small compared to the number of data points. We solve a relaxation optimization problem using a modification of Newton’s method, where the cost function depends on the squared distance matrix. The solution obtained in our scheme achieves a lower complexity and can perform better if we use it as an initial guess for an interactive local search of other higher precision localization scheme. Simulation results show the effectiveness of our approach.


Localization in Wireless Sensor Networks
Localization of sensor nodes is a challenging issue in wireless sensor networks (WSNs) for intelligent Internet of Things (IoT). Localization systems are not only for location identification but also provide information for routing, density control, tracking and a number of other communication network applications which integrate in many technologies of IoT. In general, the sensor positioning process has two steps. First, we choose the signal parameters to describe the location information between sensors. Second, we use some parametric methods for estimating the sensor positions based on the signal parameters in the first step. For the first step, GPS-based localization systems have a high degree of accuracy and offer global location information. However, alternative solutions for GPS are required, which are cost effective, rapidly deployable and can operate in diverse conditions, especially for indoor or non-line-of-sight environments. For these reasons, more suitable localization algorithms for WSNs need to be investigated. Some existing localization approaches such as time-of-arrival (TOA), time-difference-of-arrival, and angle-of-arrival achieve accurate localization results, but require high cost, complicated timing and synchronization. On the other hand, received signal strength (RSS) measurements are quite simple to obtain [1][2][3].
In this paper, we study sensor network localization problems in embedding dimension, given anchors and the RSS information between sensors. The anchors are located at fixed known sensor positions, and distances between unknown sensors and anchors are estimated from RSS measurements. Our goal is to approximate all sensor locations by using only a partial Euclidean distance matrix for the second step. Recently, many solution techniques for this problem have been proposed, such as semidefinite programming relaxation (SDP) and solvers [4][5][6], multidimensional scaling (MDS) and its improvements [7][8][9], heuristics [10], Euclidean distance matrix completion (EDMC) [11], to name a few. However, most of previous approaches are not scalable and require high computational complexity, so that we try to reduce computational complexity by using a modified iterative Newton's algorithm for a cost function.
Notations: The following notations are used throughout our paper. (·) T and (·) −1 denote the transpose and the inverse operations. • denotes the Hadamard product. || · || 2 and || · || F denote the 2 -norm and the Frobenius matrix norm, respectively. The notation w ∼ N (µ, σ 2 ) means that w is distributed according to the normal Gaussian distribution with the mean µ and the covariance σ 2 . The operator diag(A) returns a column vector of the main diagonal elements of the matrix A. For vectors, arg min(x) returns the indices of the smallest elements in x. For two arbitrary symmetric matrices A and B, A B means A − B is positive semidefinite. grad f and Hess f denote the gradient and the Hessian vectors representing the first and the second partial derivatives.

Euclidean Distance Matrix Completion
Consider a set of n points (or sensors) with locations x 1 , · · · , x n ∈ R r (in practice, r = 2 or 3). Denote D ∈ R n×n as an Euclidean distance matrix whose entries are the squared pairwise distances between n points by setting It is easy to see that the rank of D is upper bound by r + 2 which is very small compared to the number of data points n, especially when n becomes large. The set of all Euclidean distance matrices in R n×n is denoted as EDM(n). We associate a weighted undirected graph G = (V, E, W) with D, where the vertex set V = {1, · · · , n}, the edge set E = {(i, j) : i = j, and D is specified}, and the positive edge weight W = (w ij ) with w ij = D ij for all (i, j) ∈ E. Let H be the matrix whose entries are Given the observation data D, the low-rank distance matrix completion problem states as min D∈EDM(n) To find out the relationship between positive semidefinite matrices and Euclidean distance matrices, let Then, we can observe that Hence, D = κ(Y), where κ(·) is a linear map defined as Here, S n = {Y ∈ R n×n : Y = Y T } and e ∈ R n is the vector of all ones. The problem Equation (3) is equivalent to min Since rank(D) ≤ r + 2, we can solve a sequence of non-convex problems as the following rank constrained semidefinite optimization problem By screening the value from ρ = 1 to ρ = r + 2, the solution presented in [12] guarantees a monotonic convergence to the original solution of Equation (7).

Contribution of the Paper
In order to solve Equation (7), several approaches have been studied in [4][5][6][7][8][9]. The conventional MDS method transforms the pairwise distance information into the relative coordinates of sensor nodes [4][5][6]. Thus, the global solution only obtained when we use all pairwise distance measurements of sensors, which is impractical. To overcome this problem, the authors in [7,8] proposed the distributed weighted MDS (WMDS/dwMDS) method, while the authors in [9] proposed the MDS-MAP. Both methods are basically based on either the linearized least square estimator for TOA information or the biased RSS-based estimators, so that the estimated solutions can be obtained explicit and much easier than the conventional MDS. However, their accuracy is not good when the variance of the measurement noise is large.
In this paper, our goal is to design a numerical estimator especially for the RSS-based localization systems. This estimator is based on the RSS measurements rather than pairwise distance information. We first formulate the sensor localization problem as a relaxation of SDP, then the solution can be numerically obtained via a modification of the Newton's method as long as the RSS measurements are valid. This result is mainly used for coarse localization strategy, i.e., reducing the region of interest and computation time for the fine localization stage, where the estimated solutions can be used as a starting point for other fine localization algorithms such as dwMDS and MDS-MAP, etc. This localization scheme provides a low-cost, low-complexity, and easy-implementation solution, thus there also exists tradeoff between the location accuracy and the computation complexity compared with other techniques.
The organization of this paper is as follows. In Section 2, we formulate the sensor network localization problem as a reduction of SDP. In Section 3, we describe our approach for node localization in WSNs by using a modification of the Newton's method. Then, numerical evaluation results are presented in Section 4, followed by concluding remarks in Section 5.

Problem Formulation of Sensor Network Localization
Consider a sensor network consisting of n wireless sensors at the locations p 1 , · · · , p n ∈ R r , where r is the embedding dimension. As defined in the previous section, D represents the Euclidean distances between sensors within given radio range R. The embedding dimension is the smallest integer r satisfying We denote the locations of m known anchor nodes as a 1 , · · · , a m ∈ R r . Let A T = [a 1 , · · · , a m ], and X T = [p 1 , · · · , p n−m ] represents the unknown sensor nodes. In the sensor network localization problems, we ignore the distinction between the anchors and the other sensors, hence we identify a i with p n−m+i (i = 1, · · · , m) and set We also assume that there are a sufficient number of anchors, which makes our problem not be realized in a smaller embedding dimension.
Denoting the distance between the unknown node i and the anchor j as τ ij , the corresponding RSS measurement can be expressed according to the following radio propagation path loss model in dB's [13].
where L ij = P T − P ij is the path loss, P T is the transmission power, L 0 denotes the path loss value at the reference distance τ 0 , γ is the path loss exponent, and ν ij is a Gaussian random variable representing the log-normal shadowing effect, ν ij ∼ N (0, σ 2 ij ). The distance τ ij can be computed through the RSS measurement by the maximum likelihood (ML) estimation. Then, the corresponding ML estimator τ ij is given by [5] Note that if we define ξ ij = exp ν ij 10 log 10 e , where e is the Euler's number, then 10 log 10 ξ ∼ N (0, σ 2 ij ). The Equation (12) can be rewritten as By using the Gaussian moment generating function, we can derive which shows that the estimator is biased. In order to obtain the unbiased estimator, we let , which yields the following variance of Var(τ ij ) = After obtaining the estimated distanceτ ij , we can apply the SDP estimator [6] or its relaxation [5,14] to reconstruct the distance matrix D from the RSS measurements. By partitioning, the n × n Euclidean squared distance matrix can be expressed as follows.
where D 11 is the (n − m) × (n − m) distance sub-matrix between the unknown locations, D 21 = D T 12 is the distance sub-matrix between the anchor locations and the unknown locations, and D 22 is the distance sub-matrix between the anchor locations. However, since at any time an unknown node i is only in the communication range of small subset of the anchor nodes, the matrix of RSS measurements is partially known. Thus, the matrix D in Equation (8) is incomplete and is affected by noises, i.e., D = D + N. Our goal is to reconstruct the complete distance matrix D from the RSS measurements. From this matrix, we can recover the unknown node locations. The problem formulation is briefly described as follows. Given A ∈ R m×r and n × n distance matrix D with corresponding adjacency matrix H, we need to solve min D∈EDM(n)

Main Results
The problem in Equation (16) cannot be solved easily due to its non-convex nature and high computational complexity. Then, we try to reformulate it into another form without leaving the minimum. We suppose that P ∈ R n×r satisfies Pe = 0 (centroid) and H • (PP) T = H • D. Based on the fact that any rank-r positive semidefinite matrix admits a factorization Y = PP T [11], we rewrite Equation (16) as For simplicity, let , and Y = [y ij ]. Therefore, we denote the objective function as Following Equation (5), we have where i, j = 1, · · · , n. This problem is now reformulated as an constrained optimization for the cost function Equation (18), which can be solved efficiently by using the Newton's method [15][16][17]. One of the motivations behind this method for optimization is to describe it as a sequence of second-order Taylor expansions and minimization.
where ξ = P − P k . The iteration that computes an approximate solution to the system of equations grad f (P) = 0, is updated as In practice, one does not compute P k+1 by explicitly computing [Hess f (P k )] −1 and then multiplying by grad f (P k ), since it is computationally inefficient. Instead, it is more practical to solve the system of linear equations Hess f (P k )ξ k = − grad(P k ) for unknown ξ k . There are also many approaches such as matrix decompositions or other algebra techniques [18] in calculating [Hess f (P k )] −1 to reduce its complexity when computing matrix operators. In general, the Newton's method will find the minimum in one iteration and has extremely fast convergence, if Hess f (·) is positive definite at the minimum and the initial guess for P is close enough to this point. However, in some cases it may return any critical points of f , that is, minima, maxima, or saddle points. To avoid these situations, we make some modifications to the Newton's method, so that it finds the critical points which are actually the minima. For search direction, we perform a line search, i.e., for each iteration we solve Equation (22) and update where t is small enough to get fast convergence. It can be 1, 1 2 or 2 −j , j = 0, 1, 2, · · · for some special cases [19].
It is also necessary to find an effective method to solve Equation (22). The Cholesky factorization is one of standard effective methods for solving linear systems compared to the Gaussian elimination or the LU decomposition. Note that any symmetric and positive definite matrix A can be expressed as A = LΣL T for some unit lower triangular matrix L and diagonal matrix Σ with positive entries on the diagonal. For each step, we decompose Hess f (P k ) as Hess f (P k ) = L k Σ k L T k , then solve the following three sub-equations.
However, when Hess f (P k ) is not positive definite, a Cholesky decomposition cannot be performed since its eigenvalues are sometimes very small negative numbers caused by rounding in computing and noise in the data. To overcome this kind of problems, we define a positive integer number p to avoid the computation error by replacing Σ k by p when Σ k is non-positive definite, and continue the factorization. If the absolute value of grad f (P k ) is small or P k+1 is close to P k , then we terminate P k+1 as the minimum, otherwise we continue the iteration. We always wish to construct P k with the convergence to the minimum at a rapid rate, so that few iterations are needed until the stopping criterion is satisfied. However, this has to be counterbalanced with the computation cost per iteration because of the trade-off between fast convergence and higher computational cost per iteration.
In Equation (17), we can discard the constraint Pe = 0 during the computational process. The reason is explained as follows. Let us assume that P ∈ R n×r is partitioned as where P 1 ∈ R (n−m)×r and P 2 ∈ R m×r . Since H • (PP) T = H • D, we have κ(P 2 P T 2 ) = κ(AA T ). If A T e = P 2 e = 0, it imples that P 2 P T 2 = AA T , i.e., there exists an orthogonal transformation Q satisfying P 2 Q = A. It was well-known as the ProCrustes problem [4]. Hence, the problem Equation (17) becomes an unconstrained optimization and it can be solve efficiently by the modified Newton's method. However, the equality constraint Pe = 0 reduces the flexibility of P and can help the convergence to the global minimum in some cases.
In summary, we use the following procedures to solve the localization problem Equation (17) as illustrated in Figure 1.

Remark 1.
The modification of the Newton's method is illustrated by the solution path on the function in [19], which is able to follow the shape of the U-valley and converges to the minimum using only finite difference gradients. For the deployment of other adaptive methods for the mesh refinement or the discrete Newton-like method, we will consider these methods as well as their convergence issues in the future work.

Remark 2.
It has been shown that the conventional MDS has the time complexity of O(n 3 ), while the WMDS [7] performed well for sparse networks but it is about two orders of magnitude slower than the MDS for larger networks (more than 100 node-networks). This is because the refinement strategy in [7] leads to a tradeoff between solution quality and computation cost. For the dwMDS algorithm [8], the authors proved that it takes O(nL + n 2 Ld thr e dwMDS ) time, where L is the number of iterations required until the stopping rule is satisfied, d thr = O(n 1/r ) is the threshold distance, and e dwMDS depends on d thr in a nonlinear way. A drawback of the dwMDS is that complexity, convergence time and initial estimate requirements for each transmission time depend sensitively on heavy weights and range measurement accuracy. For our approach, it can be verified that the computing for matrix decomposition takes O(n κ ), where 1 ≤ κ < 3, thus the total cost can be computed as O(Kn κ ), where K is the number of iterations. Thus, our approach compares favorably to the conventional methods in terms of the complexity.

Numerical Evaluation
In this section, the effectiveness of the proposed localization approach based on the modified Newton's method is numerically investigated. Moreover, its location estimator is employed and compared. We consider a WSN consisting of 20 nodes placed in an area [−20, 20] × [−20, 20] m 2 with 8 of them are anchors and its geometry is illustrated in Figure 2 In the first experiment, we investigate the performance of the proposed algorithm as shown in Figure 3. It shows the effectiveness of the location estimation of our proposed scheme for a single trial with σ = 3.98 dB [13] in two cases: (a) presence of anchor locations and (b) absence of anchor locations. By ''presence of anchor locations", we mean that these locations are used in the initial guess for P, and inverse for ''absence of anchor locations" case. A maximum RMSEs of 1.4854 meter and 1.8833 meter are observed over 100 simulations in Figure 3a,b, respectively. It indicates that better initial points leads to better final results for our proposed scheme.  For the purpose of comparison, the RMSEs of our scheme and several schemes versus σ are shown in Figure 4. Each result is collected after 1000 independent runs. After obtaining the distance information in Equation (12), we get the standard SDP estimator (SDP is used in the context of linear matrix inequalities and can be efficiently solved by interior point methods, which has the complexity of O(N 4.5 ln(1/ )), here, is a given solution precision) by using the CVX package [18,19]. From the figure, we observe that our solution achieve the solution that better than the conventional MDS, WMDS, and dwMDS methods when the noise is large, and quite comparable estimation accuracy with the standard SDP for a small value of σ, since our solution depends on how good initial points are. It can be improved if we can choose a good parameter to set up. Figure 4b give us such an example. However, note that in this particular simulation with MATLAB implementation, the proposed method is faster than the others in running time.
Standard deviation (dB) 1 1.5 2 2.5   The results indicate that our scheme is more suitable for the case of limited anchor node information. This is a good property for some scenarios where anchor node positions are not available. Thus, the obtained solution can be provided as a good initial point for other fine localization algorithms.

Conclusions
In contrast to many existing algorithms developed for the Euclidean distance matrix completion problem or its relaxation, we cast the problem under unconstrained formulation in terms of the location vectors directly and propose a new scheme using a modified Newton's method instead. We manage the complexity by organizing the gradient and Hessian matrices. The theory is simple, but the results are effective and robust for localization in WSNs for intelligent IoT.