Localization in Internet of Things Network: Matrix Completion Approach Luong Nguyen, Sangtae Kim
Byonghyo Shim
Department of Electrical and Computer Engineering Seoul National University Email: {ltnguyen,stkim}@islab.snu.ac.kr
Department of Electrical and Computer Engineering Institute of New Media and Communications Seoul National University Email:
[email protected] Abstract—In this paper, we propose a matrix completion algorithm for Internet of Things (IoT) localization. In the proposed algorithm, we recast Euclidean distance matrix completion problem as an unconstrained optimization in smooth Riemannian manifold and then propose a nonlinear conjugate gradient method on this manifold to reconstruct Euclidean distance matrix. The empirical results show that the proposed algorithm is effective and also outperforms state-of-the-art matrix completion algorithms both in noise and noiseless scenarios.
I.
I NTRODUCTION
Recently, Internet of Things (IoT) has attracted much attention due to its wide variety of applications, such as healthcare, surveillance, and environmental monitoring. One of essential components in the IoT is wireless sensor network, in which environmental data (e.g. temperature, humidity, and object movements) is collected and processed using hundreds of sensor nodes [1], [2]. In order to respond and react to the environmental data, location information of sensor nodes should be available at the basestation (a.k.a., data center, sensor fusion, access point) [3]. Since the action in IoT networks, such as fire alarm, energy transfer, emergency request, is made primarily on the data center, an approach to identify the location information of whole nodes at the data center is of importance. In this approach, dubbed as localization at data center, each sensor node measures the distance information of adjacent nodes and then send it to the data center. Using the obtained distance information, data center constructs a map of sensor nodes [2]. In order to perform the localization at the data center, pairwise distance information between each sensor pair should be provided. It has been shown that if absolute locations of a few sensor nodes (often called anchor nodes) are provided, one can find out the location information of sensor nodes accurately [3]. Main problem of the localization at the data center is that the data center might not have enough distance information to identify the locations of sensor nodes. In general, one cannot recover the original Euclidean distance matrix D from a subset of its entries since there are infinitely many completion options for the unknown entries. However, it is now well-known that if D is a rank-k matrix, then D can be recovered from Dobs with reasonably number of measured entries which is not far away from information theoretic limit. [4]. The problem to recover a low-rank matrix D from the small
number of known entries is described as e min rank(D), e ∈ Rn×n D
s.t.
(1)
e = Dobs . PE (D)
where PE is the sampling operator defined by Aij if (i, j) ∈ E [PE (A)]ij = 0 otherwise. In our case, E = {(i, j) : kxi − xj k2 ≤ r} is the set of the observed indices where r is the radio communication range. Due to the non-convexity of the rank function, it is computationally infeasible to solve this problem. Recently, many efforts have been made to transform the problem into more tractable form. One popular way is to use the nuclear e ∗ instead of rank(D). e The nuclear norm is the sum norm kDk e ∗ = P σi ). Note that the nuclear of singular values (i.e., kXk i norm minimization is the tightest convex relaxation of the (unrealizable) rank minimization problem in (1). Although this approach can be solved by the semidefinite programming (SDP) [4], computational overhead is still burdensome for the practical systems. Another approach is to use the least-squared minimization: min e ∈ Rn×n D
s.t.
1 e 2 kPE (D)
− Dobs k2F ,
(2)
e ≤ η, rank(D)
where η is the rank of the original matrix. Since the rank constraint in (1) is replaced by the Frobenius norm based cost function, this approach is effective in the noisy scenario. In recent years, various approaches to find a solution of (2) have been suggested. In [5], alternating least squares technique has been proposed. In [6], an approach to solve the problem over the smooth Riemannian manifold of rank-η matrices has been proposed. In this paper, we propose a matrix completion algorithm referred to as conjugate gradient based localization in smooth Riemannian manifold (CGL-SRM) for pursuing efficiency in the reconstruction of the location information of sensor nodes. In CGL-SRM, we formulate the matrix completion problem in (2) into the unconstrained optimization in smooth Riemannian manifold. On the designed manifold, we employ a nonlinear conjugate gradient algorithm. From the numerical results, we show that the proposed CGL-SRM algorithm recovers the distance matrix accurately both in noiseless and noisy conditions.
x2 x3
r
Fig. 1.
x1
The illustration of the sensor networks in IoT where each sensor has the same wireless communication range r.
We summarize notations used in this paper. For a given matrix A, diag(A) is the vector formed by its main diagonal and Sym(A) is defined by Sym(A) = 12 (A + AT ). For an orthogonal matrix Q ∈ Rn×k with n>k, we define its orthogonal complement Q⊥ ∈ Rn×(n−k) such that [ Q Q⊥ ] forms an orthonormal matrix. Given a function f : Y → f (Y), ∇Y f (Y) is the Euclidean gradient of f (Y) with respect to Y, i.e., [∇Y f (Y)]ij = ∂f∂y(Y) . ij II.
r
M ATRIX C OMPLETION OVER S MOOTH R IEMANNIAN M ANIFOLD
In this section, we present the proposed CGL-SRM algorithm to complete the sparse Euclidean distance matrix. In order to transform the matrix completion problem in (2) into an unconstrained optimization problem on the manifold, we exploit the smooth manifold structure of the low-rank symmetric PSD matrices. A. Problem Model From the definition of pairwise distance xj k22 = xTi xi + xTj xj − 2xTi xj , we have D = κ(XXT ),
d2ij
= kxi − (3)
where κ(XXT ) = 1diag(XXT )T + diag(XXT )1T − 2XXT . The next lemma follows immediately from (3). Lemma II.1. When n sensor nodes are distributed in kdimensional Euclidean space and n ≥ k, then rank(D) ≤ k + 2. From this lemma, (2) can be rewritten as min e ∈ Rn×n D
s.t.
1 e 2 kPE (D)
− Dobs k2F ,
e ≤ k + 2. rank(D)
Since the rank constraint in (4) is true, we have min e ∈ Rn×k X
1 e eT 2 kPE (κ(XX ))
− Dobs k2F .
eX e T , we Recalling that Y = X 1 2 min 2 kPE (κ(Y)) − Dobs kF where Y Y∈Y
(5)
further have eX eT : = {X
e ∈ Rn×k }. X B. Optimization over Riemannian Manifold Since we assume that sensor nodes are randomly distributed in k-dimensional Euclidean space (i.e., entries of X are i.i.d.), rank(X) = k is ensured almost surely. Thus, we eX eT : X e ∈ can strengthen the constraint set Y to Ye = {X e = k} and the modified problem is Rn×k , rank(X) min e Y∈Y
1 2 kPE (κ(Y))
− Dobs k2F ,
(6)
In the sequel, we denote f (Y) = 21 kPE (κ(Y)) − Dobs k2F . Now, if we define S = {U ∈ Rn×k : UT U = Ik }1 and L = {eye([λ1 ... λk ]T ) : λ1 ≥ λ2 ≥ ... ≥ λk >0}, then an element Y ∈ Ye can be expressed as Y = QΛQT where Q ∈ S and Λ ∈ L. That is, Ye = {QΛQT : Q ∈ S, Λ ∈ L}.
(7)
It can be shown that Ye is a smooth Riemannian manifold [8, Proposition 1.1]. Since the tools we use in this work is based on the differential geometry, we briefly introduce key ingredients to describe the proposed CGL-SRM algorithm. First, the tangent space of Ye is characterized in the following lemma.
(4) 1S
is an orthogonal Stiefel manifold embedded in Rn×k [7].
Lemma II.2. For the manifold Ye defined in (7), the tangent space TY Ye at Y is B CT QT TY Ye = [Q Q⊥ ] C 0 QT⊥ n = QBQT + QQTp + Qp QT : B ∈ Rk×k , o BT = B, Qp = Q⊥ C, and C ∈ R(n−k)×k . A metric on the tangent space TY Ye is defined as the matrix e that inner product between two tangent vectors β 1 , β 2 ∈ TY Y, is, < β 1 , β 2 >= tr(β T1 β 2 ). Lemma II.3. For a given matrix A, orthogonal projection PTY Ye (A) of A on the tangent space TY Ye is PTY Ye (A) = PQ Sym(A) + Sym(A)PQ − PQ Sym(A)PQ , where PQ = QQT . To express the concept of moving in the direction of a tangent space yet staying on the manifold, we employ an operation called retraction. Note that the retraction operation is a mapping from TY Ye to Ye that preserves the gradient at Y [7]. Definition II.4. The retraction RY (β) of a vector β ∈ TY Ye onto Ye is defined as RY (β) = arg min kY + β − ZkF .
(8)
e Z∈Y
e we focus only on the symmetric part. Since RY (β) is on Y, For a given matrix A, if we take the eigenvalue decomposition of its symmetric part Sym(A) = PΣPT , then Wk (A) is defined as T Wk (A) = PΣ+ (9) kP , T where Σ+ . After σ1+ ... σk+ 0 ... 0 k = eye some computations, one can show that RY (β) = Wk (Y + β). C. Conjugate Gradient based Localization
(10)
In the proposed CGL-SRM method, we apply the CG method to solve (6). Note that the update formula of the conventional CG algorithm is Yi+1 = Yi + αi Pi ,
(11)
where αi is the stepsize and Pi is the conjugate direction. In our problem, since the conjugate direction Pi should be on the tangent space, we use a retraction operator to map this into the element in the manifold. Using (10), we have Yi+1 = Wk (Yi + αi Pi ).
(12)
By this operation, we can ensure that the updated point Yi+1 e is on Y. Now what remains is to compute Pi and αi . First, in order to obtain Pi , we need to obtain the Riemannian gradient gradf (Y) of f (Y). gradf (Y) is expressed as gradf (Y) = PTY Ye (∇Y f (Y)). In the traditional CG algorithm, the conjugate direction Pi has the update formula as Pi = −gradf (Yi ) + βi Pi−1 ,
Algorithm 1: The proposed CGL-SRM algorithm 1 Input: Dobs , PE , τ : tolerance, σ ∈ (0 1) 2 T : number of iterations e i = 1, tangent vector P0 3 Initialize: Y1 ∈ Y, 4 while i ≤ T do 5 Ri = PE (κ(Yi )) − Dobs 6 ∇Y f (Yi ) = 2eye(Sym(Ri )1) − 2Ri 7 gradf (Yi ) = PTY Ye (∇Y f (Yi )) i 8 Hi = gradf (Yi ) − PTY Ye (gradf (Yi−1 )) i 9 h =< Pi , Hi > 10 βi = h12 < hHi − 2Pi kHi k2F , gradf (Yi ) > 11 Pi = −gradf (Yi ) + βi PTY Ye (Pi−1 ) i 12 Find a stepsize αi >0 such that 13 f (Yi )−f (RYi (αi Pi )) ≥ −σαi < gradf (Yi ), Pi > 14 Yi+1 = Wk (Yi + αi Pi ) 15 Di+1 = κ(Yi+1 ) 16 if kPE (Di+1 ) − Dobs kF ≤ τ then 17 break 18 end 19 i=i+1 20 end b = QΛQT 21 Obtain the eigenvalue decomposition Y 1/2 22 Z = QΛ 23 Output: Z
(13)
where βi is the conjugate update parameter. However, gradf (Yi ) and Pi−1 lie in two different vector spaces TYi Ye e so in order to add two vectors together, we need and TYi−1 Y, e Thus, the update formula of CGLto project Pi−1 onto TYi Y. SRM is Pi = −gradf (Yi ) + βi PTY
i
e (Pi−1 ). Y
(14)
For the stepsize αi , we apply Armijo’s rule. In this rule, a finite number of trials are performed to search a stepsize αi that loosely satisfies min f (RYi (αPi )). To be specific, we α>0 find the stepsize αi satisfying f (RYi (αi Pi )) ≤ f (Yi ) + σαi < gradf (Yi ), Pi >, where σ is a given constant 0