1

arXiv:1702.04054v2 [] 3 Aug 2017

Notice: This work has been submitted to the IEEE for possible
publication. Copyright may be transferred without notice, after which
this version may no longer be accessible.

August 4, 2017

DRAFT

2

Localization in the Internet of Things Network:
A Low-Rank Matrix Completion Approach
Luong Trung Nguyen, Sangtae Kim, Byonghyo Shim
Information System Laboratory
Department of Electrical and Computer Engineering, Seoul National University
Email: {ltnguyen,stkim,bshim}@islab.snu.ac.kr

Abstract
Location awareness, providing ability to identify the location of sensor, machine, vehicle, and
wearable device, is a rapidly growing trend of hyper-connected society and one of key ingredients for
internet of things (IoT). In order to make a proper reaction to the collected information from devices,
location information of things should be available at the data center. One challenge for the massive
IoT networks is to identify the location map of whole sensor nodes from partially observed distance
information. This is especially important for massive sensor networks, relay-based and hierarchical
networks, and vehicular to everything (V2X) networks. The primary goal of this paper is to propose
an algorithm to reconstruct the Euclidean distance matrix (and eventually the location map) from
partially observed distance information. By casting the low-rank matrix completion problem into the
unconstrained minimization problem in Riemannian manifold in which a notion of differentiability can
be defined, we are able to solve the low-rank matrix completion problem efficiently using a modified
conjugate gradient algorithm. From the analysis and numerical experiments, we show that the proposed
method, termed localization in Riemannian manifold using conjugate gradient (LRM-CG), is effective
in recovering the Euclidean distance matrix for both noiseless and noisy environments.

This work was sponsored by the National Research Foundation of Korea (NRF) grant funded by the Korean government
(MSIP) (2016R1A2B3015576) and Electronics and Telecommunications Research Institute (ETRI) grant funded by the Korean
government [16ZI1100, Wireless Transmission Technology in Multi-point to Multi-point Communications].
A part of this paper was presented at the Information Theory and Applications (ITA) workshop, 2016 [1].

August 4, 2017

DRAFT

3

I.

I NTRODUCTION

Recently, Internet of Things (IoT) has received much attention for its plethora of applications,
such as healthcare, surveillance, automatic metering, and environmental monitoring. In sensing
the environmental data (e.g. temperature, humidity, pollution density, and object movements),
wireless sensor network consisting of hundreds to thousands sensor nodes is popularly used
[2], [3], [4]. In order to make a proper reaction to the collected environmental data, location
information of sensor nodes should be available at the data center (basestation) [5], [6]. Since
actions in IoT networks, such as fire alarm, energy transfer, and emergency request, are made
primarily on the data center, an approach to identify the location information of whole nodes
at the data center has been received much attention. In this approach, henceforth referred to as
centralized localization, each sensor node measures the distance information of adjacent nodes
and then sends it to the data center. Then the data center constructs a map of sensor nodes using
the collected distance information [7]. In measuring the distance, various modalities, such as
received signal strength indication (RSSI) [8], time difference of arrival (TDoA) [9], and angle
of arrival (AoA) [4], have been popularly used. These approaches are simple and effective in
measuring the short-range distance and can also be used for indoor environments.
When we consider the centralized localization in IoT, there are two major issues to be
addressed. First, location information of the sensor node obtained by this approach is local,
meaning that the location information is true in the relative sense. Thus, proper adjustment
of the location information is needed to identify the absolute (true) location information. In
fact, since the local location of a sensor might be different from the absolute location by some
combinations of translations, rotations, and reflections, absolute locations of a few sensor nodes
(anchor nodes) are needed to transform the local locations into the absolute locations [5], [6].
It has been shown that when the number of anchor nodes is enough (e.g., four or more anchor
nodes in R2 ), one can identify the absolute location of sensor nodes [10]. For the sake of
completeness, we provide a brief summary of the standard procedure to obtain the absolute
location from the pairwise distance information. Let xi ∈ Rk×1 (i = 1, ..., n) be the absolute
locations of n sensor nodes randomly distributed in k-dimensional Euclidean space (typically
k = 2 or 3) and dij = kxi − xj k2 be the pairwise distance between the sensor nodes i and j,

August 4, 2017

DRAFT

4

then the Euclidean distance matrix D is defined as

0 d212 d213

 d2
2
 21 0 d23

2
2
D=
 d31 d32 0
 .
..
..
 ..
.
.

d2n1 d2n2 d2n3

· · · d21n




· · · d22n 


· · · d23n 
.

.
..
.. 
.

··· 0

Without loss of generality, we set the first sensor node as a reference in the local coordinate
ei of sensor node i is x
ei = xi − x1 and the corresponding local
system. Then a local location x

location matrix is

e =
X
where 1 =

h

1 1 ··· 1

iT

h

e1 x
e2 x
e3 · · · x
en
x

= X − 1xT1 ,

∈ Rn×1 and X =

h

iT

x1 x2 · · · xn

(1)
iT

∈ Rn×k is the absolute

e one can easily find a
location matrix of sensor nodes. By forming a Gramian matrix of X,
e To be specific, by denoting Y
e =X
eX
e T , we have
relationship between the entries of D and X.
eTi x
ej
yeij = x

= (xi − x1 )T (xj − x1 )

1
=
kxj − x1 k22 + kxi − x1 k22 − kxi − xj k22
2
1 2
(d + d21j − d2ij ),
=
2 i1
e One can easily show that
where yeij is the (i, j)-th entry of Y.

e = 1 (De1 1T + 1eT D − D),
Y
1
2
h
iT
e = P zi yeij zj = k P zi x
e
ei k22 ≥ 0 for all vector z, Y
where e1 = 1 0 · · · 0 . Since zT Yz
i,j

i

e = QΛQT ,
is positive semi-definite (PSD) matrix. If we express its eigendecomposition as Y
b = QΛ1/2 becomes a local location matrix of sensor nodes. Finally,
then the matrix square root X
when the absolute locations of a few anchor nodes are provided, absolute locations of all sensor
b Readers are referred to [5], [6] for more details.
nodes can be recovered from X.

Second and perhaps more serious issue is that the data center does not have enough distance

information to identify the locations of sensor nodes. For various reasons, such as the power
August 4, 2017

DRAFT

5

x1





7
 
= 
 
9

x3





 11 


=



7

1

x5





 15 


=



6

3

2

5




4

2
x2 =  
7

x4

Fig. 1.





 12 


=



4

Sensor nodes deployed to measure not only environment information but also their pairwise distances. The observed

distances are represented by two-sided arrows. The shadow areas represents the radio communication range of the sensor nodes.

outage of a sensor node or the limitation of radio communication range, only partial distance
information is available at the data center. This situation can also happen in the hierarchical or
relay-based IoT networks where an intermediate node sends partial distance information to the
data center. Also, in the vehicle networks it might not be possible to measure the distance of
all adjacent vehicles when a vehicle is located at the dead zone. To illustrate this scenario, we
consider a sensor network consisting of five sensor nodes in Fig. 1. We see that only a small
number of pairwise distances is measured, and hence there are many unknown entries in the

August 4, 2017

DRAFT

6

observation matrix Dobs :



Dobs

0


 d2
 21

2
=
 d31

 ?

?

d212 d213

?

0

?

?

?

0

d234

?

d243

0

?

d253 d254

?




? 


d235 
,

d245 

0

where the question mark ? indicates unknown entries of D. In general, one cannot recover the
original matrix D from the knowledge of a subset of its entries since there are infinitely many
completion options for the unknown entries. However, it is now well-known that even though
we have only partial measurements on D, as long as D is a low-rank matrix, then D might be
recovered from Dobs [11], [12]. As will be discussed in the next section, rank of D is at most
k + 2 and k is 2 or 3 for most cases. Thus, the Euclidean distance matrix D can be readily
considered as a low-rank matrix. The problem to recover a low-rank matrix D from the small
number of known entries is expressed as
min

e ∈ Rn×n
D

s.t.

e
rank(D),

e = PE (Dobs ).
PE (D)

(2)

where PE is the sampling operator given by

A if (i, j) ∈ E
ij
[PE (A)]ij =
 0
otherwise.

Let r be the radio communication range, then E = {(i, j) : kxi − xj k2 ≤ r} would be the set
of observed indices. Since solving (2) is numerically infeasible due to the non-convexity of the
rank function, many efforts have been made over the years to relax this problem into a tractable
form to solve. Candes and Recht showed that a low-rank matrix can be recovered from partial
measurements by using the nuclear norm minimization (NNM) [12]. The nice feature of this
approach is that the NNM can be cast as a semidefinite programming (SDP) and thus can be
solved via a polynomial time complexity algorithm [12], [13]. However, computational overhead
is still burdensome and hardly scale to the problem size.
An alternative approach is to use the Frobenius-norm minimization given by
min

e ∈ Rn×n
D

s.t.
August 4, 2017

1
e
kPE (D)
2

− PE (Dobs )k2F ,

e ≤ η,
rank(D)

(3)

DRAFT

7

where η is the rank of the original matrix. Since the equality constraint in (2) is replaced by the
Frobenius norm-based cost function, this approach is suited for the noisy scenario. In fact, since
an addition of noise is unavoidable in the distance measurement process of the IoT networks,
an approach robust to the measurement noise is desired. Also, this approach is good fit for the
situation where the rank constraint is known a priori. In recent years, various approaches to find
a solution of (3) have been suggested. In [14], [15], approaches inspired by the greedy recovery
strategy of the compressed sensing (CS) have been proposed. In these approaches, the set of
rank-one matrices that best represents the original matrix is found by the iterative process. In
[16], [17], alternating least squares (ALS) technique has been suggested. In this technique, a lowrank matrix is factorized into the product of lower dimensional matrices and then the problem
is solved via an alternative minimization. In [18], matrix completion problem is modeled as
an unconstrained optimization problem over the smooth Riemannian manifold and a nonlinear
conjugate gradient method has been used to solve the problem.
The main goal of this paper is to propose a Euclidean distance matrix completion technique
optimized for IoT localization. Instead of solving the Frobenius norm minimization problem in
(3) as it is, we express the Euclidean distance matrix D as a function of the fixed rank positive
semidefinite matrix. Since the set of these matrices forms a Riemannian manifold in which the
notation of differentiability can be defined, we can recycle, after a proper modifications, an
optimization algorithm in the Euclidean space. In fact, in order to solve the Frobenius-norm
minimization problem, we propose a modified conjugate gradient algorithm, referred to as localization in Riemannian manifold using conjugate gradient (LRM-CG), whose solution lies on the
Riemannian manifold. We show from the recovery condition analysis that the sequence generated
by LRM-CG converges to the original Euclidean distance matrix under suitable conditions. We
also show from numerical experiments that LRM-CG is effective in recovering the Euclidean
distance matrix from partial measurements and also scalable to the dimension of the matrix.
We briefly summarize notations used in this paper. Pn×n denotes the set of n × n symmetric
positive semidefinite (PSD) matrices. < β 1 , β 2 > is the inner product between two vectors (or
matrices) β 1 and β 2 , i.e., < β 1 , β 2 >= tr(β T1 β 2 ). diag(A) is the vector formed by the main
diagonal of a matrix A. Sym(A) and Skew(A) are the matrices formed by Sym(A) = 21 (A+AT )
and Skew(A) = 12 (A − AT ) for any square matrix A, respectively. Note that A = Sym(A) +
Skew(A). eye(a) is a diagonal matrix whose diagonal entries are elements of a. For an orthogonal
August 4, 2017

DRAFT

8

matrix Q ∈ Rn×k with n>k, we define its orthogonal complement Q⊥ ∈ Rn×(n−k) such that
h
i
Q Q⊥ forms an orthonormal matrix. Given a function f : Y ∈ Rn×n → f (Y) ∈ R,

∇Y f (Y) is the Euclidean gradient of f (Y) with respect to Y, i.e., [∇Y f (Y)]ij = ∂f∂y(Y)
. For a
ij
h
iT
given matrix A = a1 a2 · · · an
∈ Rn×n , the vectorization of A, denoted by vec(A),
iT
h
is defined as vec(A) = aT1 aT2 · · · aTn
. Let f : A → f (A) and g : B → g(B) be

mappings, then the composite g ◦ f of f and g is defined by g ◦ f : A → (g ◦ f )(A) = g(f (A)).

E(i,j) is the standard element
of Rn×n whose the (i, j)-th element is one and the rest are zeros.

0 1
 for R2×2 .
For example, E(1,2) = 
0 0
II.

T HE LRM-CG A LGORITHM

In this section, we present the proposed LRM-CG algorithm. By exploiting the smooth
Riemannian manifold structure for the set of the low-rank symmetric PSD matrices, we formulate
the matrix completion problem (3) as an unconstrained optimization problem on the smooth
Riemannian manifold. Roughly speaking, smooth manifold is a generalization of the Euclidean
space on which a notion of differentiability exists. For more rigorous definition, see, e.g., [19],
[20]. A smooth manifold together with an inner product, often called a Riemannian metric,
forms a smooth Riemannian manifold. Since the smooth Riemannian manifold is a differentiable
structure equipped with an inner product, we can use all necessary ingredients for solving
optimization problems with quadratic cost function, such as Riemannian gradient, Hessian matrix,
exponential map, and parallel translation [19]. Therefore, optimization techniques in Euclidean
vector space (e.g., steepest descent, Newton method, conjugate gradient method) can be readily
applied for solving a problem in the smooth Riemannian manifold.
A. Problem Model
From the definition of pairwise distance d2ij = kxi − xj k22 = xTi xi + xTj xj − 2xTi xj , we have
D = g(XXT ),

(4)

where g(XXT ) = 2Sym(diag(XXT )1T − XXT ). In the example illustrated in Fig. 1, X is

T
h
iT
7 2 11 12 15
 ,
X = x1 x2 x3 x4 x5
=
(5)
9 7 7 4 6
August 4, 2017

DRAFT

9

and D is



0

29

20

50

73





 29 0 81 109 170 




T
.
D = g(XX ) = 
20
81
0
10
17




 50 109 10 0
13 


73 170 17 13
0

(6)

The next lemma follows immediately from (4).

Lemma 2.1: If n sensor nodes are distributed in k-dimensional Euclidean space and n ≥ k,
then rank(D) ≤ k + 2.

Proof: From (4), we have rank(D) = rank(g(XXT )), which gives
1diag(XXT )T + diag(XXT )1T
) + rank(2XXT )
2
≤ 2rank(1diag(XXT )T ) + rank(XXT )

rank(D) ≤ rank(2

≤ 2 + k,
where the last inequality is because rank(XXT ) = rank(X) ≤ k and rank(abT ) ≤ 1 for any
two vectors a and b.
From this lemma, (3) can be rewritten as
min

e ∈ Rn×n
D

s.t.

1
e
kPE (D)
2

− PE (Dobs )k2F ,

e ≤ k + 2.
rank(D)

(7)

eX
e T )) ≤ k + 2 for any X,
e we can further simplify the problem as
Since rank(D) = rank(g(X
min

e ∈ Rn×k
X

1
eX
e T ))
kPE (g(X
2

eX
e T , we have
Recalling that Y = X
min

Y∈Y

1
kPE (g(Y))
2

− PE (Dobs )k2F .

− PE (Dobs )k2F ,

(8)

(9)

eX
eT : X
e ∈ Rn×k }.
where Y = {X

When the sensor nodes are randomly distributed in k-dimensional Euclidean space, rank of

the location matrix X is k almost surely1 . Thus, we can strengthen the constraint set from Y to
1

Consider the case that sensor nodes are randomly distributed in 2D Euclidean space, then rank(X) = 1 if and only if all of

nodes are collinear. This event happens if there exists a constant ρ such that xi1 = ρxi2 for any i-th row. The probability of
n
Q
this event
P (xi1 = ρxi2 ) = [P (x11 = ρx12 )]n becomes negligible when the number of sensor nodes are sufficiently large.
i=1

August 4, 2017

DRAFT

10

t∈R

B = γ ′ (t)
B

t=0

0
RY (B)
Y

Y

γ(t)

γ(t)

TY Ye

Ye

(a)

Fig. 2.
e
Y.

TY Ye
(b)

Ye

e and (b) the retraction operator RY at a point Y in the embedded manifold
Illustration of (a) the tangent space TY Y

eX
eT : X
e ∈ Rn×k , rank(X)
e = k}, and as a result, we have
Ye = {X
min

e
Y∈Y

1
kPE (g(Y))
2

− PE (Dobs )k2F ,

(10)

In the sequel, we denote f (Y) = 21 kPE (g(Y)) − Dobs k2F for notational simplicity.
B. Optimization over Riemannian Manifold
Let S = {Q ∈ Rn×k : QT Q = Ik }2 and L = {eye([λ1 · · · λk ]T ) : λ1 ≥ λ2 ≥ · · · ≥ λk >0}.
e one can express Y = QΛQT , and thus an alternative representation of
Then, for given Y ∈ Y,

Ye is

Ye = {QΛQT : Q ∈ S, Λ ∈ L}.

(11)

It has been shown that Ye is a smooth Riemannian manifold [21, Ch.5] [22]. Our approach to

solve the problem in a smooth Riemannian manifold is beneficial in two major respects: First,
one can easily compute the gradient of the cost function in (10) using the matrix calculus.
Second, we can use an algorithm in the Euclidean space to solve the problem (10).
2

S is an orthogonal Stiefel manifold embedded in Rn×k [20, Theorem 5.12].

August 4, 2017

DRAFT

11

Since our work relies to a large extent on properties and operators of differential geometry,
we briefly introduce tools and ingredients to describe the proposed algorithm. Since Ye is an
embedded manifold in the Euclidean space Rn×n , its tangent spaces3 are determined by the
e Put it formally, for
derivative of its curves, where the curve γ of Ye is a mapping from R to Y.
e the tangent space of Ye at Y, denoted TY Y,
e is defined as TY Ye = {γ ′ (0) :
a given point Y ∈ Y,

e γ(0) = Y} (see Fig. 2). In the following lemma, we characterize the tangent
γ is a curve in Y,
e
space of Y.
Lemma 2.2: For the manifold Ye defined by (11), the tangent space TY Ye at Y is



 
h

i B CT
QT
T
k×k
(n−k)×k
e




TY Y =
.
: B = B ∈ R ,C ∈ R
Q Q⊥


C 0
QT⊥

(12)

Proof: See Appendix A.

A metric on the tangent space TY Ye is defined as the matrix inner product < B1 , B2 >=
e We next define the orthogonal projection
tr(BT B2 ) between two tangent vectors B1 , B2 ∈ TY Y.
1

e which will be used to find the closed-form expression
of a matrix A onto the tangent space TY Y,

of Riemannian gradient in Subsection II-C.

Definition 2.3: The orthogonal projection onto TY Ye is a mapping PTY Ye : Rn×n → TY Ye such
e
that for a given matrix A ∈ Rn×n , < A − P e (A), B >= 0 for all B ∈ TY Y.
TY Y

The following proposition provides a closed form expression of the orthogonal projection

operator.
Proposition 2.4: For a given matrix A, orthogonal projection PTY Ye (A) of A onto the tangent
space TY Ye is
PTY Ye (A) = PQSym(A) + Sym(A)PQ − PQ Sym(A)PQ ,

(13)

where PQ = QQT .
Proof: See Appendix B.
In order to express the concept of moving in the direction of a tangent space while staying
on the manifold, an operation called retraction is used. As illustrated in Fig. 2(b), the retraction
operation is a mapping from TY Ye to Ye that preserves the gradient at Y [23, Definition 4.1.1].
3

Tangent vector on the smooth Riemannian manifold is a generalization of the notion of tangent vector to curves and surfaces

in Euclidean space. A collection of all tangent vectors at each point of smooth manifold forms a tangent space at that point,
which is indeed a vector space.
August 4, 2017

DRAFT

12

Definition 2.5: The retraction RY (B) of a vector B ∈ TY Ye onto Ye is defined as
RY (B) = arg min kY + B − ZkF .
e
Z∈Y

(14)

In obtaining the closed form expression of RY (B), an operator Wk keeping k largest positive
eigenvalues of a matrix, referred to as eigenvalue selection operator, is needed. Since the
e RY (B) should be a symmetric PSD matrix with rank k.
projection RY (B) is an element of Y,

Thus, for given square matrix A, we are only interested in the symmetric part Sym(A). If we
denote the eigenvalue decomposition (EVD) of this as Sym(A) = PΣPT and the k top most

eigenvalues of this as σ1 ≥ σ2 ≥ · · · ≥ σk >0 , then Wk (A) is defined as

where Σk = eye

h

σ1 ... σk

value selection operator Wk .

Wk (A) = PΣk PT ,
(15)

iT
. RY (B) is concisely expressed using the eigen0 ... 0

Lemma 2.6:
RY (B) = Wk (Y + B).

(16)

Proof: See Appendix C.
e we
Finally, to develop the conjugate gradient algorithm over the Riemannian manifold Y,

need the Euclidean gradient of the cost function f (Y).

Lemma 2.7: Euclidean gradient ∇Y f (Y) of f (Y) with respect to Y is
∇Y f (Y) = 2eye(Sym(R)1) − 2Sym(R),

(17)

where R = PE (g(Y)) − PE (Dobs ).
Proof: See Appendix D.
C. Localization in Riemannian Manifold Using Conjugate Gradient (LRM-CG)
In order to solve the problem in (10), we use the conjugate gradient (CG) method. CG method
is widely used to solve sparse symmetric positive definite linear systems [24].

August 4, 2017

DRAFT

13

∇Y f (Y)
gradf (Y)

Y

TP Ye

Ye

Fig. 3. Riemannian gradient gradf (Y) is defined as the projection of the Euclidean gradient ∇Y f (Y) onto the tangent space
TY Ye while the Euclidean gradient is a direction for which the cost function is reduced most in Rn×n , Riemannian gradient is
e
the direction for which the cost function is reduced most in the tangent space TY Y.

First, noting that PE and g are linear mappings, one can easily show that
f (Y) =
=
=
(a)
=

(b)
=

1
kPE (g(Y)) − Dobs k2F
2
!!
X
1
kPE g
− Dobs k2F
yij E(i,j)
2
i,j
1 X
k
yij PE (g(E(i,j))) − Dobs k2F
2 i,j


1 X
k
yij vec PE (g(E(i,j))) − vec(Dobs )k22
2 i,j

1
kAvec(Y) − bk22 ,
2

(18)

h
iT
where (a) is because kMkF = kvec(M)k2 and (b) follows from vec(Y) = y11 y21 · · · ynn ,
h

 i
b = vec(Dobs ), and A = vec PE (g(E(1,1) )) · · · vec PE (g(E(n,n))) .
In (18), we see that the cost function f (Y) has the quadratic form of a sparse symmetric

positive definite system, and thus the CG algorithm would be an efficient means to solve the

August 4, 2017

DRAFT

14

problem. The update equation of the conventional CG algorithm in the Euclidean space is
Yi+1 = Yi + αi Pi ,

(19)

where αi is the stepsize and Pi is the conjugate direction. Note that the stepsize αi is chosen
by the line minimization technique (e.g., Armijo’s rule [25], [26]) and the search direction Pi
of the CG algorithm is chosen as a linear combination of the gradient and the previous search
direction to generate a direction which are conjugate to the previous ones. In doing so, one can
avoid unnecessary searching of directions that have been searched over and achieve the speedup
in the algorithm [27], [24].
e the conjugate
Since we consider the optimization problem over the Riemannian manifold Y,

direction Pi lies on the tangent space. Thus, to make sure that the update point Yi+1 lies on

the manifold, we need a retraction operator. The update equation after applying the retraction
operation is
Yi+1 = RYi (αi Pi )
= Wk (Yi + αi Pi ).

(20)

As studied in Lemma 2.6, the eigenvalue selection operator Wk guarantees that the updated point
Yi+1 lies on the manifold.
We next consider the conjugate direction Pi of LRM-CG. In the conventional nonlinear CG
algorithm, conjugate direction Pi is updated as
Pi = −∇Y f (Yi ) + βi Pi−1 ,

(21)

where βi is the conjugate update parameter4 . Since we optimize over the Riemannian manifold
e (21) needs to be modified. First, we need to use the Riemannian gradient of f (Y) instead of
Y,
e
the Euclidean gradient ∇Y f (Y) since we find the search direction on the tangent space of Y.
Riemannian gradient, denoted gradf (Y), is distinct from ∇Y f (Y) in the sense that it is defined
on the tangent space TY Ye (see Fig. 3).

Definition 2.8: Let f be the function differentiable everywhere in the Riemannian manifold
e The Riemannian gradient gradf (Y) of f at Y is defined as the unique element in TY Ye
Y.
4

There are a number of ways to choose βi . See, e.g., [24], [28], [29], [30].

August 4, 2017

DRAFT

15

satisfying
< B, gradf (Y) >=< B, ∇Y f (Y) >,

(22)

e
where B is any element in TY Y.

As shown in Fig. 3, gradf (Y) ∈ TY Ye is a component of the Euclidean gradient ∇Y f (Y) in
e In other words, gradf (Y) is the projection of ∇Y f (Y) onto TY Y.
e Indeed, from Definition
TY Y.
e Hence,
2.3, < B, ∇Y f (Y) − P e (∇Y f (Y)) >= 0 for any matrix B ∈ TY Y.
TY Y

< B, PTY Ye (∇Y f (Y)) >=< B, ∇Y f (Y) > .

(23)

From (22) and (23), it is clear that
gradf (Y) = PTY Ye (∇Y f (Y)).

(24)

Second, since the Riemannian gradient gradf (Yi ) and previous conjugate direction Pi−1 lie
e we need to project Pi−1 onto the tangent space
on two different vector spaces TY Ye and TY Y,
i

i−1

TYi Ye5 before performing a linear combination between of two. In view of this, the conjugate

direction update equation of LRM-CG is

Pi = −gradf (Yi ) + βi PTY Ye (Pi−1 ).

(25)

i

Finally, in choosing the stepsize αi in (20), we use the Armijo’s rule, a widely used line search
strategy. Note that the Armijo’s rule is a simple yet effective way to find a stepsize αi minimizing
the cost function f , that is, αi ≈ min f (Wk (Yi + αi Pi ) [25], [26].
α>0

The proposed LRM-CG algorithm is summarized in Algorithm 16 .

D. Numerical Experiments
In this subsection, we investigate numerical performance of the proposed LRM-CG for both
noiseless and noisy scenarios. In our simulation, we compare LRM-CG with following matrix
completion algorithms:
5

In transforming a vector from one tangent space to another, an operator called vector transport is used (see Definition 8.1.1

in [19]). For an embedded manifold of Rn×n , vector transport is the orthogonal projection operator [19]. Hence, the vector
transport of Pi−1 is the orthogonal projection of Pi−1 onto TY Ye
i

6

One can also find the matlab source code at http://islab.snu.ac.kr/publication.html.

August 4, 2017

DRAFT

16

Algorithm 1: LRM-CG algorithm
1

Input: Dobs : the observed matrix,

PE : the sampling operator,
ǫ: tolerance,
µ ∈ (0 1): given constant,
T : number of iterations.
2

Initialize: i = 1,

Y1 ∈ Ye : initial point,

P1 : initial conjugate direction.
3

While i ≤ T do

4

Ri = PE (g(Yi )) − PE (Dobs )

// Generate residual matrix

5

∇Y f (Yi ) = 2eye(Sym(Ri )1) − 2Ri

// Compute Euclidean gradient

6

gradf (Yi ) = PTY

// Compute Riemannian gradient

i

7

e (∇Y f (Yi ))
Y

Hi = gradf (Yi ) − PTY

i

8

h =< Pi , Hi >

9

βi =

10

1
h2

e (gradf (Yi−1 ))
Y

< hHi − 2Pi kHi k2F , gradf (Yi ) >

Pi = −gradf (Yi ) + βi PTY

i

11

e (Pi−1 )
Y

Find a stepsize αi >0 such that

// Compute CG coefficient
// Compute conjugate direction
// Perform Armijo’s line search

f (Yi ) − f (RYi (αi Pi )) ≥ −µαi < gradf (Yi ), Pi >
12

Yi+1 = RYi (αi Pi )

// Perform retraction

13

Di+1 = g(Yi+1 )

// Compute updated Euclidean distance matrix

14

If kPE (Di+1 ) − PE (Dobs )kF <ǫ then

15

Exit from while loop

16

End If

17

Obtain Q and Λ using the eigendecomposition

Yi+1 = QΛQT
18
19

b = QΛ1/2
X

// Find updated locations of sensor nodes

i=i+1

20

End While

21

b
Output: X

August 4, 2017

DRAFT

17

• ADMiRA [14]: this algorithm uses greedy projection to identify a set of rank-one matrices
that best represents the original matrix.
• LMaFit [16]: this is a nonlinear successive over-relaxation matrix completion algorithm
based on nonlinear Gauss-Seidel method.
• LRGeomCG [18]: this algorithm is essentially CG method defined over the Riemannian
manifold of low rank matrices (but not necessarily positive definite).
• SDP [12]: this algorithm solves the NNM problem via semidefinite programming.
• TNN-ADMM [31]: as an improved version of the NNM problem, this algorithm solves
the truncated NNM problem via an alternating direction method of multipliers.
In the noiseless scenario, we generate an n × k location matrix X whose entries are sampled
independently and identically from a uniform distribution in unit interval, and then map X into the
Euclidean distance matrix D = g(XXT ). As mentioned, entries dij of D are known (observed)
if dij ≤ r (r is the radio communication range). In the noisy scenario, a noise matrix N ∈ Rn×n
is added to the Euclidean distance matrix D. Elements nij of N are sampled independently and
identically from the Gaussian distribution with zero mean and variance σ 2 =

SNR
kDk2F
10− 10 .
n(n−1)

In

obtaining the performance for each algorithm, we perform at least 1000 independent trials.
1) Convergence Efficiency: As a performance measures, we use two types of the mean square
errors (MSE): MSE at sampled entries (MSEs ) and MSE at all entries (MSEa ), which are
defined respectively as
1
b − PE (D)k2 ,
kPE (D)
F
|E|
1 b
kD − Dk2F ,
=
n2

MSEs =
MSEa

where |E| is the cardinality of the sampling set E (see (2)). In Fig. 4, we plot the log-scale MSE
as a function of the number of iterations for the 2-dimensional location vectors. Note that the
performance results are obtained for the scenario where 200 sensor nodes randomly distributed.
We observe that the log-scale MSE decreases linearly with the number of iteration, which in
turn implies that the MSE decreases exponentially with the number of iterations. For example, if
r = 0.5, it takes less than 25 iterations to achieve 10−2 and 53 iterations to achieve 10−4 . Also,
as expected, required number of iterations to achieve a given performance level gets smaller as
r increases. To further examine the convergence speed of LRM-CG, we measured the running

August 4, 2017

DRAFT

18

101
r = 0.5
r = 0.7
r = 0.8

0

10

10-1

10-2

10-2

MSE

a

r = 0.5
r = 0.7
r = 0.8

0

10-1

MSE

s

10

101

10-3

10-3

10-4

10-4

10-5

10-5

10-6

Fig. 4.

0

10

20

30

40

50

60

70

80

90

10-6

0

10

20

30

40

50

The number of iteration

The number of iteration

(a)

(b)

60

70

80

90

The MSE performance of LRM-CG: (a) M SEs and (b) M SEa for k = 2 (2-dimensional location vectors).

time of LRM-CG for three distinct matrix dimensions (500×500, 1000×1000, and 5000×5000).
The running time is measured using the MATLAB program on a personal computer (Intel Core
i5 CPU with 3.4 GHz). As indicated in Table I, LRM-CG recovers 500 × 500 Euclidean distance
matrix in a few seconds. Even for 1000 × 1000 matrix, it takes less than one minute to recover
the matrix accurately.
2) Recovery Performance: We next investigate the recovery performance of LRM-CG for
both noiseless and noisy scenarios. In Fig. 5, we plot the recovery performance of the noiseless
scenario for k = 2 and 3 as a function of the sampling ratio, where the sampling ratio is defined
as a ratio between the number of observed pairwise distances and total number of pairwise
distances. We observe from the figure that LRM-CG performs well in all cases, achieving target
performance with smaller number of measurements than conventional techniques requires. In
Fig. 6, we plot the performance of LRM-CG as a function of the signal-to-noise ratio (SNR) for
the noisy scenario. We observe that LRM-CG outperforms conventional approaches by a large
margin and the gain increases sharply with SNR.

August 4, 2017

DRAFT

19

TABLE I.

E XPERIMENTAL RESULTS FOR E UCLIDEAN DISTANCE MATRIX COMPLETION .

Size (n × n)

Euclidean dimension

r

time (s)

Number of iterations

M SEa

2

0.5

9

90

1.0 × 10−6

3

0.7

10

92

3.4 × 10−6

2

0.5

38

79

1.3 × 10−6

3

0.7

45

98

1.4 × 10−6

2

0.5

1897

109

4.6 × 10−6

3

0.7

2542

126

7.4 × 10−6

500 × 500

1, 000 × 1, 000

101

100

100

10-1

10-1

10-2

10-2

MSE

a

101

MSE

a

5, 000 × 5, 000

10-3

10-3

LRM-CG
ADMiRA
LMaFit
LRGeomCG
SDP
TNN-ADMM

10-4

10-5

10-6
0.1

0.2

0.3

10-4

10-5

0.4

0.5

Sampling ratio

(a)

Fig. 5.

0.6

0.7

0.8

0.9

10-6
0.1

LRM-CG
ADMiRA
LMaFit
LRGeomCG
SDP
TNN-ADMM

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Sampling ratio

(b)

The MSE performance of the matrix completion algorithms for noiseless scenario for (a) 2-dimensional and (b)

3-dimensional location vectors.

August 4, 2017

DRAFT

20

10-1

10-1

MSE

MSE

a

100

a

100

10-2

10-2
LRM-CG
ADMiRA
LRGeomCG
SDP
TNN-ADMM

LRM-CG
ADMiRA
LRGeomCG
SDP
TNN-ADMM

10-3

5

10

15

20

25

30

35

SNR

(a)

40

45

10-3

5

10

15

20

25

30

35

40

45

SNR

(b)

Fig. 6. The MSE performance of the algorithms for noisy scenario for (a) 2-dimensional and (b) 3-dimensional location vectors.

E. Computational Complexity
In this subsection, we examine the computational complexity of LRM-CG in terms of the
number of floating point operations (flops). As discussed in Section II-B, LRM-CG computes
Euclidean gradient, Riemannian gradient, and the retraction in each iteration.
In order to compute the Euclidean gradient ∇Y f (Yi ) in (17), we first need to consider the

computation of Yi from the (i − 1)-th iteration. Since Yi = QΛQT (Q is a n × k matrix and

Λ is a k × k diagonal matrix), it requires 2k multiplications and (k − 1) additions to compute
k
P
yij =
λt qit qjt so that the associated computational complexity is (3k − 1) flops. Further, from
t=1

(4), we need to compute [g(Y)]ij = yii + yjj − yij , which requires (9k − 1) flops. The residual
matrix Ri = PE (g(Yi )) − Dobs requires 9k|E| flops (|E| is the number of the observed entries

of Dobs ). In addition, since it takes (9k + 2)|E| flops to compute Sym(Ri ) = 12 (Ri + RTi ), it
requires at most (9k + 4)|E| + n − 1 flops to compute 2eye(Sym(Ri )1). Since the cardinality
of Sym(Ri ) is |E|, computational complexity of ∇Y f (Yi ) = 2eye(Sym(Ri )1) − 2Sym(Ri ) in
(17) is at most (9k + 5)|E| + n − 1.
August 4, 2017

DRAFT

21

Second, recalling that the Riemannian gradient gradf (Yi ) is an orthogonal projection of
e we need to estimate the computational complexity of the
∇Y f (Yi ) onto the tangent space TY Y,
i

orthogonal projection operator PTY Ye . In computing PTY Ye (A) for an n × n matrix A, we need
i

i

T

Sym(A), B = Sym(A)Q, and C = Q Sym(A)Q, which require (2k − 1)n2 , 2n2 + (2n − 1)kn,

and (2n − 1)kn + (2n − 1)k 2 flops, respectively. Then, from (13), we have

PTY Ye (A) = QQT Sym(A) + Sym(A)QQT − QQT Sym(A)QQT
i

= QBT + BQT − QCQT ,
which requires O(kn2 + k 2 n + k 3 ) flops.
Finally, in applying Armijo’s rule to find the stepsize αi , we need to compute the retraction
operation RYi (αi Pi ). From (16), the retraction operation is obtained via the eigenvalue selection
operator Wk and this requires the EVD of (Yi + Pi ). In general, computational complexity of
e
the EVD for a n×n matrix is expressed as O(n3 ). However, using Yi = QΛQT and Pi ∈ TY Y,
i

e we have
one can simplify the EVD operation. First, since Pi ∈ TYi Y,

 
T
h
i B i CT
i  Q 

Pi = Q Q⊥
.
Ci 0
QT⊥

Thus,

Yi + Pi


 
T
h
i B i + Λ CT
i  Q 

= Q Q⊥
Ci
0
QT⊥



T
h
i Bi + Λ I
Q


= Q (Q⊥ Ci ) 
T
I
0 (Q⊥ Ci )


T
h
i
Q

= Q (Q⊥ Ci ) KΛ′ KT 
T
(Q⊥ Ci )
= Q′ Λ′ Q′T .

(26)

(27)

From (27), we

 see that the EVD of (Yi + Pi ) is simplified to the EVD of the 2k × 2k matrix
h
i
B +Λ I
 i
, which requires only 8k 3 flops. Also, computing Q⊥ Ci and Q′ = Q (Q⊥ Ci ) K
I
0
needs (2n − 2k − 1)kn and (4k − 1)kn flops, respectively. Thus, computational complexity of
the retraction operation is 8k 3 + (2n − 2k − 1)kn + (4k − 1)kn, which is O(kn2 ) for k ≪ n.
August 4, 2017

DRAFT

22

In summary, computational complexity of the proposed algorithm per iteration is O(k|E| +

kn2 + k 2 n + k 3 ) = O(kn2 ). Since k = 2 or 3 in our problem, complexity per iteration can be

expressed as O(n2 ) flops.

III.

R ECOVERY C ONDITION A NALYSIS

In this section, we analyze a recovery condition under which the LRM-CG algorithm recovers
the Euclidean distance matrices accurately. Overall, our analysis is divided into two parts. In the
first part, we analyze a condition ensuring the successful recovery of the sampled (observed)
b
entries, i.e., kPE (D)−P
E (D)kF = 0. In the second part, we investigate a condition guaranteeing

b − DkF = 0. By exact recovery,
the exact recovery of the Euclidean distance matrices, i.e., kD

we mean that the output Di of LRM-CG converges to the original Euclidean distance matrix D.
∞
Definition 3.1: For a sequence of matrices {Di }∞
i=1 , if lim kDi − DkF = 0, we say {Di }i=1
i→∞

converges to D. Further, we say

{Di }∞
i=1

converges linearly to D with convergent rate λ if there

exists λ (1>λ ≥ 0) satisfying
lim

i→∞

kDi+1 − DkF
= λ.
kDi − DkF

A. Convergence of LRM-CG at Sampled Entries
In this subsection, we show that {PΩ (Di )}∞
generated by LRM1 , sequence of matrices

 0 29 20 


CG at sampled points, converges to PΩ (D). For example, if D =  29 0 81  and E =


20 81 0


 0 29 20 


{(1, 2), (1, 3)}, then PE (D) =  29 0 0 . Thus, we will show that


20 0 0
 


2
2
0
d12 (∞) d13 (∞)   0 29 20 

 

 2
lim PE (Di ) =  d21 (∞)
 =  29 0 0  = PE (D).
0
0
i→∞
 


0
0
20 0 0
d231 (∞)
The minimal set of assumptions used for the analytical tractability are as follows:

A1 :

f (Yi ) − f (RYi (αi Pi )) ≥ −τ αi < gradf (Yi ), Pi > for τ satisfying 0<τ <1/2,

August 4, 2017

DRAFT

23

A2 :

| < gradf (RYi (αi Pi )), Pi > | ≤ −µ < gradf (Yi ), Pi > for µ satisfying τ <µ<1/2,

A3 :

ckgradf (Yi )kF ≥ k∇Y f (Yi )kF for c satisfying c>1.

In essence, A1 and A2 can be considered as extensions of the strong Wolfe’s conditions7
[32]. The assumption A1 says that the cost function f (Yi ) decreases monotonically as long
as Pi is chosen in an opposite direction of gradf (Yi ) on the tangent space TYi Ye (i.e., <
gradf (Yi ), Pi >≤ 0) (see Lemma 3.7). Note that A1 is reasonable assumption since there
always exists a stepsize satisfying this assumption.
Lemma 3.2: There exists αi >0 satisfying A1.
Proof: See Appendix E
Note that if the stepsize αi is chosen to be very small, then Yi+1 = RYi (αi Pi ) ≈ RYi (0) =
Yi , and thus
f (Yi ) − f (RYi (αi Pi )) ≈ 0.
and −τ αi < gradf (Yi ), Pi >≈ 0. In this case, A1 holds true approximately. However, there
would be almost no update of Yi so that the algorithm will converge extremely slowly. To
circumvent this pathological scenario, we use A2, which is in essence an extension of the strong
Wolfe’s condition for the Riemannian manifold. Under this assumption, αi cannot be chosen to
be very small since otherwise we have RYi (αi Pi ) ≈ Yi , and thus
| < gradf (RYi (αi Pi )), Pi > | ≈ | < gradf (Yi ), Pi > | ≥ −µ < gradf (Yi ), Pi >,
which contradicts the assumption A2.
The assumption A3 is needed to guarantee the global convergence of LRM-CG. We will
discuss more on this in Remark 3.6.
7

Consider an unconstrained minimization in Rn with a differentiable cost function f (x) (i.e., minn f (x)). The update equation
x∈R

is given by xi+1 = xi + αi pi for a stepsize αi and a descent direction pi . The well-known strong Wolfe’s conditions is given
by
f (xi ) − f (xi+1 )

≥

−τ αi < ∇x f (xi ), pi >,

| < ∇x f (xi+1 ), pi > |

≤

−µ < ∇x f (xi ), pi >,

for some constants 0<τ <µ<1.

August 4, 2017

DRAFT

24

Our first main result, stating successful recovery condition at sampled entries, is formally
described in the following theorem.
Theorem 3.3 (strong convergence of LRM-CG): Let {Di = g(Yi )}∞
i=1 be the sequence of the
matrices generated by LRM-CG and D be the original Euclidean distance matrix. Under A1,
A2, and A3, {PE (Di )}∞
i=0 converges linearly to PE (D).
Remark 3.4 (strongly convergent condition in Rn ): Note that lim kPE (Di ) − PE (D)kF = 0
i→∞

is equivalent to
lim k∇Y f (Yi )kF = 0.

(28)

i→∞

This condition is often referred to as the strongly convergent condition of the nonlinear CG
algorithms in the vector space. The equivalence can be established by the following sandwich
lemma.
Lemma 3.5:
√
2kPE (Di ) − PE (D)kF ≤ k∇Y f (Yi )kF ≤ (2 n + 2)kPE (Di ) − PE (D)kF .
Proof: See Appendix F
Remark 3.6: Recently, an attempt has been made to extend the convergent analysis of the
conventional CG algorithms (over the Euclidean space Rn ) to the Riemannian manifolds. In [33,
Theorem 4.3], it has been shown that under certain assumption,
lim inf kgradf (Yi )kF = 0.

(29)

i→∞

One can observe that the Euclidean gradient ∇Y f (Yi ) is replaced by the Riemannian gradient
gradf (Yi ). Unfortunately, the convergence of the Riemannian gradient in (29) does not imply
the convergence of Euclidean gradient in (28) because
k∇Y f (Yi )k2F = kPTY Ye (∇Y f (Yi ))k2F + kPT⊥ Ye (∇Y f (Yi ))k2F
Y

= kgradf (Yi )k2F + kPT⊥ Ye (∇Y f (Yi ))k2F ,
Y

(30)

where gradf (Yi ) = PTY Ye (∇Y f (Yi )) (see (24)). One can observe from this that the condition in
(29) is not sufficient to guarantee (28), that is, one cannot guarantee lim kPE (Di )−PE (D)kF =
i→∞

0 just from (29). However, by the introduction of A3, equivalence between (28) and (29) can

August 4, 2017

DRAFT

25

be established. We will show that the assumption A3 holds true with overwhelming probability
in Section III-C.
We are now ready to prove Theorem 3.3.
Proof of Theorem 3.3: First, we show that under A1 and A2, kPE (Di ) − PE (D)kF is nonincreasing. That is, if χ is defined by

kP ⊥ e (∇Y f (Y))kF
TY Y

 sup
k∇Y f (Y)kF
χ = Y∈{Yi }∞
i=1


1

if k∇Y f (Y)kF 6= 0

,

(31)

otherwise

then there exists γ>0 such that γ(1 − χ2 ) ≤ 1 and


kPE (Di+1 ) − PE (D)k2F ≤ 1 − γ(1 − χ2 ) kPE (Di ) − PE (D)k2F .

(32)

We need the following lemma to prove this.

Lemma 3.7: If βi is chosen based on Fletcher-Reeves’ rule8 , that is,
βi =
then

< gradf (Yi ), gradf (Yi ) >
,
< gradf (Yi−1 ), gradf (Yi−1) >

(34)

µi+1
1 − 2µ
< gradf (Yi+1 ), Pi+1 >
−
.
≤
−
kgradf (Yi+1 )k2F
1−µ
1−µ
Proof: See Appendix G.

Lemma 3.8: kgradf (Yi )k2F ≥ 8(1 − χ2 )f (Yi ).
Proof: See Appendix H.
8

In our simulation, we employ Hager-Zhang’s rule in the choice of βi to improve the empirical performance of the CG method

[30]:
βi =
where Hi = gradf (Yi ) − PTY

i

1
< hHi − 2Pi kHi k2F , gradf (Yi ) >
h2

e (gradf (Yi−1 ))
Y

(33)

and h =< Pi , Hi >. In our analysis, however, we use Fletcher-Reeves’ rule

for mathematical tractability.

August 4, 2017

DRAFT

26

We are now ready to prove (32). First, from A1, we have
f (Yi+1) ≤ f (Yi ) + τ αi < gradf (Yi ), Pi >


1 − 2µ
µi
(a)
f (Yi ) − τ αi
kgradf (Yi )k2F
+
≤
1−µ
1−µ


1 − 2µ
kgradf (Yi )k2F ,
≤ f (Yi ) − τ αi
1−µ


1 − 2µ
(b)
f (Yi ) − 8τ αi
(1 − χ2 )f (Yi ),
≤
1−µ
where (a) and (b) follow from Lemma 3.7 and Lemma 3.8, respectively. Let


1 − 2µ
γi = 8τ αi
,
1−µ
then γi >0 (since αi >0) and hence
f (Yi+1 ) ≤ (1 − γi (1 − χ2 ))f (Yi ).
Recalling that f (Yi ) = 21 kPE (Di ) − PE (D)k2F , we have


kPE (Di+1 ) − PE (D)k2F ≤ 1 − γi (1 − χ2 ) kPE (Di ) − PE (D)k2F .

By choosing γ = min γi , we get the desired result.
i

Now, what remains is to show that lim kPE (Di ) − PE (D)kF = 0 under (32). From (30) and
i→∞

(31), we have 0 ≤ χ ≤ 1, and thus we need to consider two cases:
1) χ<1 case: In this case, one can easily show that 1>(1 − γ(1 − χ2 ))1/2 . Using this together
with (32), we have
kPE (Di+1 ) − PE (D)kF
= (1 − h(1 − χ2 ))1/2 <1
kPE (Di ) − PE (D)kF
and hence
lim

i→∞

lim kPE (Di ) − PE (D)kF = 0.

i→∞

Thus, the sequence {PE (Di )}∞
i=1 converges linearly to PE (D).
2) χ = 1 case: In this case, we show that there exists j satisfying k∇Y f (Yj )kF = 0. As
discussed in Remark 3.4 and Lemma 3.5, this is a sufficient condition to guarantee the
strong convergence of LRM-CG. In this case, no further update can be made after j-th
iteration (and thus linear convergence is naturally guaranteed). To show this, we use the

August 4, 2017

DRAFT

27

contradiction argument. Suppose that k∇Y f (Yi )kF 6= 0 for all i. Then, from (31) we
should have

kPT⊥ Ye (∇Y f (Y))kF

sup

Y

k∇Y f (Y)kF

Y∈{Yi }∞
i=1

= χ = 1.

Further, from (30), we have
kPT⊥

2
e (∇Y f (Yi ))kF

= k∇Y f (Yi )k2F − kgradf (Yi )k2F

Yi Y

≤ k∇Y f (Yi )k2F −

1
k∇Y f (Yi )k2F ,
2
c

where the inequality is from A3 (c>1). Thus,
1 =

kPT⊥ Ye (∇Y f (Y))k2F

sup

Y

Y∈{Yi }∞
i=1

≤

sup
Y∈{Yi }∞
i=1

= 1−

1
,
c2

k∇Y f (Y)k2F

k∇Y f (Y)k2F − c12 k∇Y f (Y)k2F
k∇Y f (Y)k2F

which is contradiction. Thus, k∇Y f (Yi )kF = 0 for some j.
B. Exact Recovery of Euclidean Distance Matrices
So far, we have shown that the output of LRM-CG converges to the original Euclidean distance
matrix D at sampled entries (i.e., PE (D∞ ) = PE (D)). In this subsection, we show that all
entries of Di converge to that of the original Euclidean distance matrix D with overwhelming
probability.
Before we proceed, we briefly discuss the probability model of the sampling operator PE .
Let δij be a Bernoulli random variable that takes value 1 if dij ≤ r (recall that r is the radio
communication range) and 0 otherwise. Since the distance is symmetric (i.e., dij = dji ), we have
δij = δji . Also, since the diagonal entries of D are all zeros, we define δii = 0 for all i. Then,
for a matrix A, PE (A) can be expressed as
X
δij < A, ei eTj > ei eTj ,
PE (A) =
i6=j

=

X

δij aij ,

(35)

i6=j

August 4, 2017

DRAFT

28



 0

where ei is the standard basis of Rn . For example, if A =  10

17
then


 0 10 0 


PE (A) =  10 0 3 


0 3 0





 0 1 0 
 0 0 0 
 0





= 10  0 0 0  + 3  0 0 1  + 10  1





0 0 0
0 0 0
0


10 17 

0 3  and E = {(1, 2), (2, 3)},

3 0





0 0 
 0 0 0


0 0 +3 0 0 0


0 0
0 1 0





.


We now characterize the random variables δij using P (dij ≤ r). Since dij = kxi − xj k2 ,

it follows P (dij ≤ r) = P (kxi − xj k2 ≤ r). In this work, we assume that elements of xi
(locations of sensor nodes) are i.i.d. random and uniformly distributed over unit interval. By
denoting p = P (dij ≤ r), the probability mass function (PMF) of δij can be expressed as
f (δij ; p) = pδij (1 − p)1−δij .

(36)

The following lemma provides an explicit expression of p in terms of the radio communication
range r.
Lemma 3.9: If an element of the location vectors xi is i.i.d. and uniform on unit interval,
then
a) If k = 2 (2-dimensional Euclidean space),


1 4
8 3
2

if 0 ≤ r ≤ 1

πr − 3 r + 2 r
√
p=
p1 (r)
if 1 ≤ r ≤ 2 ,




1
else

b) If k = 3 (3-dimensional Euclidean space),

4

πr 3 − 3π
r 4 + 85 r 5 − 61 r 6
if 0 ≤ r ≤ 1

3
2


√


p2 (r)
if 1 ≤ r ≤ 2
p=
√
√ ,


p3 (r)
if 2 ≤ r ≤ 3




1
else

August 4, 2017

DRAFT

29

where
√
2
1
1
p1 (r) = − − 2r 2 − r 4 + (8r 2 + 1) r 2 − 1 + 2r 2 sin−1
3
2
3
2
,

+
sin−1 ( 22 −1)
r
1 + tan
2



2
−1
r2


(37)

r
15π + 37 6π + 1 2 8π 3 3π + 3 4 1 6
1
4
−1
p2 (r) = −
+
r −
r +
r + r + 2r sin
1− 2
30
2
3
2
3
r




 
2
2
44 2 16 4 √ 2
1
4
−1
4
−1
+
− r sin
− r − r
−1
r − 1 − 2r sin
15 15
5
r
r2
16
4
−
+
(38)

2 ,
3
1 + tan 21 sin−1 r22 − 1
3 1 + tan 21 sin−1 r22 − 1
r


26 44 2 8 4 √ 2
16 3
2
π 97 7 2 3 4 1 6
−1
− r − r − r +
+ r + r
r − 2 + r tan
1− 2
p3 (r) = − +
2 30 2
2
6
15 15
5
3
r
r


√
r2 − 2
3 − r2
4
2
−1
+(2 − 8r 2 ) tan−1 r 2 − 2 − (2r 4 − 4r 2 ) sin−1
+
(r
+
2r
)
sin
r2 − 1
r2 − 1
r
r
r
1
1
1
16
2
−1
4
2
−1
+
8r
tan
+
(2r
−
4r
)
sin
− r 3 tan−1
3
r 4 − 2r 2
r2 − 2
r2 − 1
16
1
1
+
−4
.
(39)

3
2
2 2
3 1 + tan 1 sin−1 3−r
1 + tan 1 sin−1 3−r
2
2
2

r −1

2

r −1

Proof: See Appendix I.
It is worth noting that p1 (r), p2 (r), and p3 (r) increase monotonically with r (see Fig. 7).
We now state our main result.
Theorem 3.10: Under the assumption A1, A2, and A3, the output sequence {Di = g(Yi )}∞
1
of LRM-CG converges globally to the Euclidean distance matrix D ( lim kDi − DkF = 0) with
i→∞

the probability at least

 
 

c
1−c
+ c log
1 − exp − (1 − c) log
1−p
p

(40)

for some constant c satisfying 0<c<1 and c<p.
Remark 3.11: From Lemma 3.9, we see that p gets close to 1 as the radio communication
range r increases. Thus, as shown in Fig. 8, the chance of recovering D increases with r.
Remark 3.12: Theoretical guarantee on the recovery of a matrix has been provided by Candes
and Recht in [12], and later improved in [34], [35]. In short, if entries of a matrix are chosen
August 4, 2017

DRAFT

30

1
k=2
k=3

0.9

Sampling parameter p

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

0.5

1

1.5

2

2.5

3

Radio communication range r

Fig. 7.

The sampling parameter p gets close to 1 as r increases. Here, elements of xi are i.i.d. random variables according to

the uniform distribution over unit interval.

at random, then the n × n matrix with rank k can be recovered with overwhelming probability

as long as the number of measurements m follows m = O(kn1.2 log(n)). The analysis in these

works is based on the assumption that observed entries are sampled i.i.d. ( and follow Bernoulli
or uniform distribution). Whereas, our analysis does not require independence assumption among
the sampled entries of D since the elements of D are related. In other words, random variables
δij do not need to be independent. For example, consider the scenario illustrated in Fig. 9. Since
the sensor node 4 is located inside the triangle formed by three sensor nodes (nodes 1, 2, and
3), one can see that d14 ≤ max(d12 , d13 ). Thus, if d12 and d13 are already known (i.e., d12 ≤ r,
d13 ≤ r), then so is d14 . In other words, P (δ14 = 1|δ12 = δ13 = 1) = 1, while P (δ14 = 1) is not
necessarily one. In our work, we do not put any assumption on the independence of the entries
of D yet show that D can be recovered exactly with overwhelming probability when r is large.
Following lemma is useful to prove Theorem 3.10.
Lemma 3.13: For a given matrix A, if the diagonal entries are zeros (i.e., aii = 0 for all i)
and kAkF <∞, then there exists a constant t (0<t<1) satisfying
tkAk2F ≤ kPE (A)k2F ,
August 4, 2017

(41)
DRAFT

31

1

1
c = 0.1
c = 0.2
c = 0.3

0.9
0.8

0.8
0.7
Probability (40)

0.7
Probability (40)

c = 0.1
c = 0.2
c = 0.3

0.9

0.6
0.5
0.4

0.6
0.5
0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Radio communication range r

1.8

2

0

0.2

0.4

0.6

(a)

Fig. 8.

0.8

1

1.2

1.4

1.6

1.8

2

Radio communication range r

(b)

The Euclidean distance matrix D can be recovered with overwhelming probability in (a) 2D and (b) 3D Euclidean

space when r is large.

 


 
mt
with the probability at least 1 − exp − (1 − mt) log 1−mt
+
mt
log
for some con1−p
p
stant m ≥ 1, provided that 0<mt<p<1.
Proof: See Appendix J.
Proof of Theorem 3.10: Let A = Di − D. Then from Lemma 3.13, we have
1
(42)
kDi − DkF ≤ √ kPE (Di ) − PE (D)kF
t
 


 
mt
with the probability at least 1 − exp − (1 − mt) log 1−mt
+
mt
log
for some con1−p
p
stant m satisfying m ≥ 1 and 0<m< pt . Combining this with lim kPE (Di ) − PE (D)kF = 0
i→∞

(Theorem 3.3), we can conclude that
lim kDi − DkF = 0,
 
 
 
1−c
with the probability at least 1 − exp − (1 − c) log 1−p + c log pc
, where c = mt.
i→∞

August 4, 2017

DRAFT

32

2

4
3

1
Fig. 9.

Suppose that the sensor node 4 is inside the triangle formed by three sensor nodes 1, 2, and 3. Then for a given r,

it can be shown that d14 ≤ max(d12 , d13 ), and thus P (d14 ≤ r|d12 ≤ r, d13 ≤ r) = 1 which is not necessarily equivalent to
P (d14 ≤ r).

C. Discussion on A3
In this section, we show that the assumption A3 (ckgradf (Yi )k2F +ǫ>k∇Y f (Yi )k2F 9 for some
c>1 and ǫ>0) holds true with overwhelming probability when r is large. In order to show this,
we first need to define the coherence, a measure of concentration in a matrix [12].
Definition 3.14 (Coherence [12]): Let Q be a subspace of Rn of dimension k and PQ be the
orthogonal projection onto Q. Then the coherence of Q is defined by
µ(Q) =

n
max kPQ ei k22 .
k 1≤i≤n

Consider a matrix A of rank k whose singular value decomposition is given by
A = UΣVT =
where U =

h

u1 · · · uk

i

and V =

h

k
X
i=1

v1 · · · vk

σui viT ,
i

(43)

are the matrices constructed by the left

and right singular vectors, respectively, and Σ is the diagonal matrix whose diagonal entries are
σi . From (43), we see that the concentration on the vertical direction (concentration in the row) is
9

When ǫ is sufficiently small, one can simply put ckgradf (Yi )k2F ≥ k∇Y f (Yi )k2F , which is the strict form of A3. Intuitively,

if the generated sequence of Riemannian gradient goes to zero ( lim kgradf (Yi )kF = 0), so does the corresponding sequence
i→∞

of the Euclidean gradient ( lim k∇Y f (Yi )kF = 0).
i→∞

August 4, 2017

DRAFT

33

determined by ui and that on the horizontal direction (concentration in the column) is determined
h
iT
by vi . For example, if one of the standard basis vector ei , say e1 = 1 0 · · · 0 , lies on
the space spanned by u1 , · · · , uk while others (e2 , e3 , · · · ) are orthogonal to this space, then it
is clear that nonzero entries of the matrix are only on the first row. Since we need to check the
concentration on both vertical and horizontal directions, we need to investigate both µ(U) and
µ(V). In this regard, the coherence of a matrix A is defined by [12]
µ(A) = max(µ(U), µ(V )).

(44)

In particular, if A is a positive semidefinite matrix, then it is clear that U = V and thus
µ(A) = µ(U).
e Then, for any c>1 and ǫ>0,
Theorem 3.15: Suppose µ(Yi ) ≤ µ0 for a given matrix Yi ∈ Y.
c2 kgradf (Yi )k2F + ǫ>k∇Y f (Yi )k2F ,

(45)




 

1 − mǫ
mǫ
+ (1 − mǫ) log
1 − exp − mǫ log
1−p
p

(46)

with probability at least

for some constant m satisfying m>0 and 0<m< 1−p
, provided that n ≥ 2cµ0 k.
ǫ
Remark 3.16: From Lemma 3.9, we see that p gets close to 1 as r increases. Thus, as shown
in Fig. 10, when r is large, (45) holds true with overwhelming probability.
Following lemmas are needed to prove the theorem.
Lemma 3.17:
k∇Y f (Y)k2F − c2 kgradf (Y)k2F
XX
δij | < B, ei eTj >< B, eu eTv >< (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |,
≤

(47)

i6=j u6=v

where B = g(Y) − D and l(A) = 2eye(Sym(A)1) − 2Sym(A).
Proof: See Appendix K.
Lemma 3.18: If n2 ≥ 4cµ(Y)2k 2 and i 6= j, then
| < (I −

cPTY Ye )l(ei eTj ), l(ei eTj )



4cµ(Y)2 k 2
,
>|≥4 1−
n2

(48)

where c>0 and l(A) = 2eye(Sym(A)1) − 2Sym(A).
August 4, 2017

DRAFT

34

1

1
mǫ = 0.9
mǫ = 0.8
mǫ = 0.7

0.9
0.8

0.8
0.7
Probability (46)

0.7
Probability (46)

mǫ = 0.9
mǫ = 0.8
mǫ = 0.7

0.9

0.6
0.5
0.4

0.6
0.5
0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

0.2

Radio communication range r

0.4

0.6

0.8

(a)

Fig. 10.

1

1.2

1.4

1.6

1.8

2

Radio communication range r

(b)

The condition (45) holds true with overwhelming probability in (a) 2D and (b) 3D Euclidean space when the radio

communication range r is large.

Proof: See Appendix L.
Lemma 3.19: Let δ1 , δ2 , · · · , δN be identically (not necessarily independently) distributed Bernoulli
random variables with P (δi = 1) = p and P (δi = 0) = 1 −p. Also, let a1 , a2 , · · · , aN be positive

values. Let q be the largest integer obeying 2q ≤ N. Then, for any ǫ>0,



 

N
X
1 − mǫ
mǫ
+ (1 − mǫ) log
,
P(
δi ai ≥ ǫ) ≤ exp − mǫ log
1−p
p
i=1
where mǫ =

PN

i=1 ai −ǫ
2q amin

(49)

with amin = min ai , provided that 0<mǫ<1 − p.
i

Proof: See Appendix M.
Now, we are ready to prove Theorem 3.15.
Proof of Theorem 3.15: Let
I = k∇Y f (Yi )k2F − c2 kgradf (Yi )k2F ,
sij = < g(Yi ) − D, ei eTj >,
August 4, 2017

DRAFT

35

and
gij =

X
u6=v

|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |.
i

In this proof, we will show that P (I ≤ ǫ) is lower bounded by the quantity in (46). First, since
P
P
δij gij ≤ ǫ) = 1, we have
δij gij from Lemma 3.17 and hence P (I ≤ ǫ |
I≤
i6=j

P(

X
i6=j

δij gij ≤ ǫ) = P (

X
i6=j

i6=j

δij gij ≤ ǫ)P (I ≤ ǫ |

= P (I ≤ ǫ,
= P(

X
i6=j

X
i6=j

X
i6=j

δij gij ≤ ǫ)

δij gij ≤ ǫ)

δij gij ≤ ǫ | I ≤ ǫ)P (I ≤ ǫ)

≤ P (I ≤ ǫ).

(50)

P
What remains is to find out a lower bound of P ( δij gij ≤ ǫ). Equivalently, we find out an
i6=j
P
upper bound of P ( δij gij ≥ ǫ). First, in order to use Lemma 3.19, we need to find a lower
i6=j

bound of gij for (i, j) ∈ Ω (Ω = {(i, j) : sij 6= 0}). Let s = min |sij |, then
gij =

X
u6=v

(i,j)∈Ω

|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |

≥ |s2ij < (I − c2 PTY Ye )l(ei eTj ), l(ei eTj ) > |

(a)
4s2 1 − 4c2 µ(Yi )2 k 2 /n2
≥

≥ 4s2 1 − 4c2 µ20 k 2 /n2 ,

where (a) follows from Lemma 3.18. Now using Lemma 3.19, we have
 




X
mǫ
1 − mǫ
P(
δij gij ≥ ǫ) ≤ exp − mǫ log
+ (1 − mǫ) log
, (51)
1−p
p
i6=j
P
and
where m = ( (i,j)∈Ω gij − ǫ)/(2q ǫc1 ) (c1 = 4s2 (1 − 4c2 µ20 k 2 /n2 )), provided that 0<m< 1−p
ǫ
n ≥ 2cµ0 k. Here, q is the largest integer obeying 2q ≤ |Ω| (|Ω| being the cardinality of Ω).
P
P
From (50) and (51), and noting that P ( δij gij ≤ ǫ) = 1 − P ( δij gij ≥ ǫ), we get the desired
i6=j

i6=j

result.

IV.

C ONCLUSION

In this paper, we have proposed an algorithm to recover the Euclidean distance matrix (and
therefore the location map) from partially observed distance information. In solving the Frobenius
August 4, 2017

DRAFT

36

norm minimization problem with a rank constraint, we expressed the Euclidean distance matrix
as a function of the fixed rank positive semidefinite matrix. By capitalizing on the Riemannian
manifold structure for this set of matrices, we could solve the low-rank matrix completion
problem using the modified conjugate gradient algorithm. We have shown from the recovery
condition analysis that the proposed LRM-CG algorithm recovers the original Euclidean distance
matrix with overwhelming probability under some suitable conditions. We have also demonstrated
from numerical experiments that the LRM-CG algorithm is effective in recovering the original
Euclidean distance matrix while exhibiting computational cost scalable to the matrix dimension.
Given the importance of the location-aware applications and services in the IoT era, we believe
that the proposed LRM-CG algorithm will be a useful tool for various localization problems.
While our work focused primarily on the centralized localization scenario, extension to the
distributed and/or cooperative scenarios would also be interesting direction worth pursuing. We
also would like to point out that our recovery guarantee analysis is not so quantitative due to the
probabilistic distance measurement model (Section III.B). By introducing the non-probabilistic
model (e.g., fixed ratio communication range model), we might come up with more quantitative
and insightful analytic results. We leave these interesting explorations for our future work.

August 4, 2017

DRAFT

37

A PPENDIX A
P ROOF

OF

L EMMA 2.2

e i.e., Ẏ ∈ TY Y.
e Also, let
Proof: Let Ẏ be a tangent vector at Y ∈ Y,



 
h

T
i B CT
Q
   : BT = B ∈ Rk×k , C ∈ R(n−k)×k .
S =
Q Q⊥ 


C 0
QT⊥

(52)

e there
Then, what we need to show is that Ẏ is an element in S. By the definition of TY Y,
exists a curve γ(t) in Ye such that Y = γ(0) and Ẏ = dtd γ(t)
. For convenience, we denote
t=0

γ(t) = Z(t). Using the eigenvalue decomposition Z(t) = Q(t)Λ(t)Q(t)T , we have
Ẏ =

d
Z(t)
dt

t=0

T

= Q̇ΛQ + QΛ̇QT + QΛQ̇T

(53)

Since QT Q = I, Q is an element of the Stiefel manifold Q = {A : AT A = I, A ∈ Rn×k }. The
tangent vector of Q at the point Q is given by [19, Example 3.5.2]
Q̇ = QΩ + Q⊥ K,

(54)

where Ω is the k × k skew-symmetric matrix (i.e., Ω = −ΩT ) and K is the (n − k) × (n − k)
matrix. From (53) and (54), we have
Ẏ = (QΩ + Q⊥ K)ΛQT + QΛ̇QT


+QΛ(QΩ + Q⊥ K)

 
h
i ΩΛ + Λ̇ + ΛΩT ΛKT
QT
 .
= Q Q⊥ 
KΛ
0
QT⊥

If we denote B = ΩΛ + Λ̇ + ΛΩT and C = KΛ, then one can easily see that Ẏ ∈ TY Ye ⊆ S.
e This implies that two vector spaces S
To complete the proof, we need to show that S = TY Y.

and TY Ye have the same dimension. Indeed, from (52), we can easily check that the dimension10
of S is 1 k(2n − k + 1), which is the dimension of Ye [21, Proposition 1.1].
2

10

The dimension of S is obtained by counting the number of independent entries of an element in S. Since B is a k × k

symmetric matrix, the number of independent entries of B is

k(k+1)
.
2

In addition, since C is an arbitrary (n − k) × k matrix,

the number of independent entries of C is (n − k)k. Thus, the dimension of S is

August 4, 2017

k(k+1)
2

+ (n − k)k = 12 k(2n − k + 1).

DRAFT

38

A PPENDIX B
P ROOF

OF

P ROPOSITION 2.4

Proof: First, we partition a matrix A into two parts: A = A1 + A2 where A1 ∈ TY Ye
e ⊥ . Then, it is clear that P e (A) = A1 and thus the goal is to find out the
and A2 ∈ (TY Y)
TY Y

closed form expression of A1 . From Lemma 2.2, there exist a symmetric matrix B ∈ Rk×k and

a matrix C ∈ R(n−k)×k such that


 
h
i B CT
QT
 .
A1 = Q Q⊥ 
C 0
QT⊥

(55)

Since < A1 , A2 >= 0, we have
0 = < A1 , A2 >

Let

= < A1 , A − A1 >

 
 

h
i B CT
h
i B CT
QT
QT



>




, A − Q Q⊥
= < Q Q⊥
C 0
QT⊥
QT⊥
C 0

  

  
h
i B CT
h
i
B CT
QT
QT










= <
A − Q Q⊥
,
Q Q⊥ >
C 0
C 0
QT⊥
QT⊥

  


h
i
B CT
B CT
QT
 ,   A Q Q⊥ − 
>.
= <
T
C 0
C 0
Q⊥

then we have




A21









h
i
QT AQ QT AQ⊥
 =   A Q Q⊥ = 
,
A22
QT⊥
QT⊥ AQ QT⊥ AQ⊥

A11 A12



QT

(56)

(57)

 
 

T
B C
A
A12
B C
 ,  11
−
>
0 = <
C 0
A21 A22
C 0
T

1
= < B, A11 − B > + < C, (A21 + AT12 ) − C >,
2

where the equality is because

 

α1 α2
β β2
, 1
 >=< α1 , β 1 > + < α2 , β 2 > + < α3 , β 3 > + < α4 , β 4 > .
<
α3 α4
β3 β4
August 4, 2017

DRAFT

39

Since B and C are chosen arbitrarily, we should have
< B, A11 − B >= 0,
1
< C, (A21 + AT12 ) − C >= 0.
2

(58)
(59)

First, it is clear from (59) that
1
C = (A21 + AT12 ).
2

(60)

Next, noting that A11 = Sym(A11 ) + Skew(A11 ), (58) becomes
0 = < B, A11 − B >
= < B, Sym(A11 ) − B > + < B, Skew(A11 ) >
(a)
=

< B, Sym(A11 ) − B > + < Sym(B), Skew(A11 ) >

(b)
=

< B, Sym(A11 ) − B >,

where (a) is because B is the symmetric matrix (i.e., B = Sym(B)), and (b) is because <
Sym(C), Skew(D) >= 0 for any matrices C and D. Since B is any symmetric matrix, we
should have
1
B = Sym(A11 ) = (A11 + AT11 ).
2
Substituting (60) and (61) into (55), we have

h
i 1 (A11 + AT )
11
A1 = Q Q⊥  2
1
(A21 + AT12 )
2

1
(AT21
2

(61)

 
+ A12 )
QT
 ,
0
QT⊥

(62)

where A11 , A12 , and A21 are the components of A (see (57)).

Now, what remains is to find a closed form expression for A1 in terms of A. First, we can
rewrite (62) as
h

i



1
(A
 2 11

+

AT11 )

1
(AT21
2

+ A12 )



T



Q
 
1
(A21 + AT12 ) 21 (A22 + AT22 ) QT⊥
2

 
h
i 0
0
QT
 
− Q Q⊥ 
0 12 (A22 + AT22 )
QT⊥


 
T
T
h
i
h
i
A11 A12
A
A12
Q
1
   + 1 Q Q⊥  11

=
Q Q⊥ 
T
2
2
A21 A22
A21 A22
Q⊥

A1 =

Q Q⊥

1
− Q⊥ (A22 + AT22 )QT⊥ .
2

August 4, 2017



T



Q
 
QT⊥

(63)
DRAFT

40

Substituting (57) into (63), we have
 
 
 
 
T
T
T
h
i
h
i
h
i
h
i
Q
QT
Q
Q
1
1
T
A1 =
Q Q⊥   A Q Q⊥   +
Q Q⊥   A Q Q⊥  
2
2
QT
QT
QT
QT
⊥

⊥

⊥

⊥

1
− Q⊥ (QT⊥ AQ⊥ + QT⊥ AT Q⊥ )QT⊥
2
1
1
(A + AT ) − Q⊥ QT⊥ (A + AT )Q⊥ QT⊥
=
2
2
= Sym(A) − Q⊥ QT⊥ Sym(A)Q⊥ QT⊥ .

Since QQT + Q⊥ QT⊥ = I and PQ = QQT , we have
A1 = Sym(A) − (I − QQT )Sym(A)(I − QQT )
= PQ Sym(A) + Sym(A)PQ − PQ Sym(A)PQ,
which is the desired result.
A PPENDIX C
P ROOF

OF

L EMMA 2.6

Proof: Our goal is to find a simple expression of the retraction operator RY (B). First, since
e we have
Z = Sym(Z) for Z ∈ Y,
kY + B − Zk2F = kY + B − Sym(Z)k2F

= kSkew(Y + B) + Sym(Y + B) − Sym(Z)k2F

(64)

= kSkew(Y + B)k2F + kSym(Y + B) − Sym(Z)k2F

(65)

= kSkew(Y + B)k2F + kSym(Y + B) − Zk2F ,

(66)

where (64) is because Sym(A) + Skew(A) = A and (65) is because < Skew(C), Sym(D) >= 0
for any C and D. Since the first term in (66) is unrelated to Z, it is clear that
RY (B) = arg min kSym(Y + B) − ZkF .
e
Z∈Y

Using the eigenvalue decomposition Sym(Y + B) = KΣKT , we have
RY (B) = arg min kKΣKT − ZkF
e
Z∈Y


= arg min kK Σ − KT ZK KT kF
e
Z∈Y

(a)
=

August 4, 2017

arg min kΣ − KT ZKkF ,
e
Z∈Y

DRAFT

41

where (a) is because kKUKT k2F = tr(KUT KT KUKT ) = kUk2F for any matrix U. Now let
RY (B) = Z∗ , Σ∗ = KZ∗ KT , and Q = KT ZK, then

Σ∗ = KZ∗ KT = arg min kΣ − QkF .
Q

(67)

Since Σ is a diagonal matrix, Σ∗ should also be a diagonal matrix. Also, Σ∗  0 and rank(Σ∗ ) =
k.11 Thus, Σ∗ is a diagonal matrix with

σ1

 0

 .
 ..


∗
Σ =
 0

 0

 .
 ..

0

only k positive entries and the rest being zero. That is,

0 ··· 0 0 ··· 0

σ2 · · · 0 0 · · · 0 


.. .. . .
.. . .
. 0 
. . .
.


(68)
0 · · · σk 0 · · · 0 
,

0 ··· 0 0 ··· 0 


.. .. . .
.. . .
. 0 
. . .
.

0 ··· 0 0 ··· 0

where σ1 ≥ σ2 ≥ · · · ≥ σk >0. Recalling that Sym(Y + B) = KΣKT , we finally have
RY (B) = KΣ∗ KT = Wk (Y + B),
where the last equality is from (15).
A PPENDIX D
P ROOF

OF

L EMMA 2.7

Proof: In general, Euclidean gradient ∇Y f (Y) can be obtained by taking partial derivatives
with respect to each coordinate of the Euclidean space. Since ∇Y f (Y) is interpreted as a matrix
whose inner product with an arbitrary matrix H becomes the Frechet differential Df (Y)[H] of
f at Y [36], that is,
Df (Y)[H] =

X
ij

hij

∂
f (Y),
∂yij

it is convenient to compute ∇Y f (Y) as a unique element of Rn×n that satisfies
< ∇Y f (Y), H >= Df (Y)[H],
11

(69)

e Z∗  0 and also Σ∗ = KT Z∗ K  0 and rank(Σ∗ ) = rank(KT Z∗ K) = rank(Z∗ ) = k.
Since Z∗ ∈ Y,

August 4, 2017

DRAFT

42

for all H. We first compute Df (Y)[H] and then use (69) to obtain the expression of ∇Y f (Y).
Note that the cost function f (Y) =

1
kPE (g(Y)
2

− Dobs k2F can be expressed as f (Y) =

h(k(Y)) = (h ◦ k)(Y) where
1
kRk2F ,
2
k(Y) = PE (g(Y)) − PE (D)
h(R) =

= (PE ◦ g)(Y) − PE (D).

(70)

(71)

Thus,
Df (Y)[H] = D(h ◦ k)(Y)[H]
= Dh(k(Y))[Dk(Y)[H]].

(72)

For any two square matrices R and A, we have
X
∂
h(R)
aij
Dh(R)[A] =
∂rij
i,j


X
∂
1
2
=
aij
kRkF
∂rij 2
i,j
=

X
i,j

=

X

aij

∂ 1X 2
(
r )
∂rij 2 p,q pq

aij rij

i,j

= < R, A > .

(73)

By choosing R = k(Y) and A = Dk(Y)[H] in (73), we can rewrite (72) as
Df (Y)[H] = < k(Y), Dk(Y)[H] >
(a)
=

< k(Y), D ((PE (g(Y)) − PE (D)) [H] >

= < k(Y), D ((PE ◦ g)(Y) − PE (D)) [H] >
(b)
=

< k(Y), D(PE ◦ g)(Y)[H] >

(c)
=

< k(Y), DPE (g(Y))[Dg(Y)[H]] >,

where (a) follows (70), (b) is because PE (D) is not a function of Y and thus the Frechet
differential of this is zero, and (c) is due to the chain rule.
August 4, 2017

DRAFT

43

Before we proceed, we remark that if S is a linear operator (i.e., S(α1 A1 +α2 A2 ) = α1 S(A1)+
α2 S(A2 )), then
DS(A)[B] = S(B)

(74)

for all matrices A and B (see Example 4.4.2 [37]).
Since PE is a linear operator, DPE (g(Y))[Dg(Y)[H]] = PE ([Dg(Y)[H]) and hence
Df (Y)[H] = < k(Y), PE (Dg(Y)[H]) >
(a)
=

< PE (k(Y)), Dg(Y)[H] >

(b)
=

< k(Y), Dg(Y)[H] >

(c)
=

< k(Y), g(H) >

(d)
=

< k(Y), 2Sym(1diag(H)T ) − 2Sym(H) >

= 2 < k(Y), Sym(1diag(H)T ) > −2 < k(Y), Sym(H) >,

(75)

where (a) is because PE is a self-adjoint operator12 , (b) is because PE (k(Y)) = PE ((PE ◦
g)(Y) − PE (Y)) = (PE ◦ g)(Y) − PE (Y) = k(Y), (c) is because g is also a linear function
and thus Dg(Y)[H] = g(H), and (d) is due to (4).
Now, the first term in (75) is
2 < k(Y), Sym(1diag(H)T ) >

(a)
=

2 < Sym(k(Y)), 1diag(H)T >

(b)
=

2 < Sym(k(Y)), diag(H)1T >

(c)
=

2 < Sym(k(Y))1, diag(H) >

(d)
=

2 < eye(Sym(k(Y))1), H >,

(76)

where (a) is because Sym() is a self-adjoint operator, (b) is because < U, W >=< UT , WT >,
(c) is because < A, b1T >= tr(AT b1T ) = tr((A1)T b) =< A1, b >, and (d) is because eye()
is the adjoint operator of diag(). Next, the second term in (75) is
− 2 < k(Y), Sym(H) >= −2 < Sym(k(Y)), H > .
12

(77)

Let A and B be two linear operators in Rn×n . If < A(A), B >=< A, B(B) >, we say A and B are adjoint to each other

in Rn×n . Further, if A ≡ B, then we say it is a self-adjoint operator.

August 4, 2017

DRAFT

44

From (75), (76), and (77), we have
Df (Y)[H] = 2 < eye(Sym(k(Y))1), H > −2 < Sym(k(Y)), H >
= < 2eye(Sym(k(Y))1) − 2Sym(k(Y)), H >

(78)

From (69) and (78), we have
∇Y f (Y) = 2eye(Sym(k(Y))1) − 2Sym(k(Y)),
which is the desired result.

A PPENDIX E
P ROOF

OF

L EMMA 3.2

Proof: If Yi is the optimal point (i.e, Yi = arg min f (Y)), then gradf (Yi ) = 0 and
Y

Yi+1 = Yi . For all αi ≥ 0, we have
f (RYi (αi Pi )) = f (Yi+1 ) = f (Yi ) + τ αi < gradf (Yi ), Pi >,
satisfying A1. Next, we consider the case where Yi 6= arg min f (Y). First, we let
Y

g(α) = f (RYi (αPi )),
h(α) = f (Yi ) + τ α < gradf (Yi ), Pi > .

(79)

Note that < gradf (Yi ), Pi >≤ 0 (see Lemma 3.7) and g(α) ≥ 0. Since g(0) = f (RYi (0)) =
f (Yi ) = h(0), g(α) and h(α) intersect at α = 0. Also, when τ varies from 0 to 1, the slope
of h(α) varies from 0 to | < gradf (Yi ), Pi > |. Since

dg(α)
dα

=< gradf (Yi ), Pi >, h(α) is
α=0

the tangential curve of g(α) at α = 0 when τ = 1. Thus, there exits 0<τ <1/2 such that h(α)
intersects g(α) at some point α>0, which means that there exist αi >0 satisfying
f (RYi (αi Pi )) = g(αi ) ≤ h(αi ) = f (Yi ) + τ αi < gradf (Yi ), Pi >,
which completes the proof.

August 4, 2017

DRAFT

45

A PPENDIX F
P ROOF

OF

L EMMA 3.5

Proof: First, a lower bound of k∇Y f (Yi )kF is given by
k∇Y f (Yi )k2F

(a)
=

k2eye(Ri 1) − 2Ri k2F

(b)
=

k2eye(Ri 1)k2F + k2Ri k2F

≥ k2Ri k2F
= 4kPE (Di ) − PE (D)k2F ,

(80)

where (a) is from (17) and (b) is from the fact that diagonal entries of Rj are all zeros and
eye(Rj 1) is a diagonal matrix. That is, positions of nonzero elements in eye(Ri 1) and Ri are
disjoint. An upper bound is obtained as follows.
k∇Y f (Yi )kF

≤ k2eye(Ri 1)kF + k2Ri kF
(a)
≤

k2Ri 1k2 + k2Ri kF

(b)
≤

2kRi kF k1k2 + 2kRi kF
√
≤ (2 n + 2)kRi kF
√
≤ (2 n + 2)kPE (Di ) − PE (D)kF ,

(81)

where (a) is because keye(b)kF = kbk2 for any vector b, and (b) is because kAbk2 ≤ kAkF kbk2
for any matrix A and any vector b. By combining (80) and (81), we obtain the desired result.

A PPENDIX G
P ROOF

OF

L EMMA 3.7

Proof: Recall from (25) that we have
Pi+1 = −gradf (Yi+1 ) + βi+1 PTY

i+1

e (Pi ).
Y

Thus,
< −gradf (Yi+1), Pi+1 > = k − gradf (Yi+1)k2F + βi+1 < −gradf (Yi+1 ), PTY

e (Pi )
Y

>

e (gradf (Yi+1 )), Pi
Y

>

i+1

= k − gradf (Yi+1)k2F + βi+1 < −PTY

i+1

(a)
=
August 4, 2017

k − gradf (Yi+1)k2F + βi+1 < −gradf (Yi+1 ), Pi >,
DRAFT

46

e Then we have
where (a) is because gradf (Yi+1 ) ∈ TYi+1 Y.
< gradf (Yi+1), Pi+1 > +kgradf (Yi+1 )k2F

= βi+1 |< gradf (Yi+1 ), Pi >|

(a)
≤

βi+1 µ < −gradf (Yi ), Pi >,

where (a) is from the assumption A2.
(Yi ),Pi >
If we denote ζi = − <gradf
, then
kgradf (Yi )k2
F

−ζi+1 kgradf (Yi+1)k2F + kgradf (Yi+1)k2F ≤ βi+1 µζi kgradf (Yi )k2F ,
and also

kgradf (Yi )k2F
| − ζi+1 + 1| ≤ µβi+1
ζi .
kgradf (Yi+1 )k2F

(82)

kgradf (Y )k2

i F
From Fletcher-Reeves rule in (34), we have βi+1 kgradf (Yi+1
= 1 and thus
)k2
F

| − ζi+1 + 1| ≤ µζi .
In other words,
ζi+1 ≥ 1 − µζi ,

(83)

and
ζi+1 ≤ 1 + µζi
ζi ≤ 1 + µζi−1
..
.

ζ2 ≤ 1 + µζ1 ,
where we set ζ1 = 1. Thus,
ζi ≤

From (83) and (84), we finally have

i−1
X

µj

j=0

1 − µi
=
.
1−µ

(84)

1 − µi
1−µ
1 − 2µ + µi+1
,
=
1−µ

ζi+1 ≥ 1 − µ

which is the desired result.
August 4, 2017

DRAFT

47

A PPENDIX H
P ROOF

OF

L EMMA 3.8

Proof: From (24), we have
gradf (Yi ) = PTY Ye (∇Y f (Yi )),
i

where ∇Y f (Yi ) is the Euclidean gradient. Let PT⊥ Ye be the orthogonal operator on the comYi
e
plement space of TY Y, then we obtain
i

k∇Y f (Yi )k2F = kPTY Ye (∇Y f (Yi )) + PT⊥

2
e (∇Y f (Yi ))kF

Yi Y

i

= kPTY Ye (∇Y f (Yi ))k2F + kPT⊥

2
e (∇Y f (Yi ))kF ,

Yi Y

i

(85)

and hence
kgradf (Yi )k2F = kPTY Ye (∇Y f (Yi ))k2F
i

= k∇Y f (Yi )k2F − kPT⊥

2
e (∇Y f (Yi ))kF .

Yi Y

Now, we define
χ=





kP ⊥

e
TY Y

sup

Y∈{Yi }∞
i=1




(∇Y f (Y))kF

if k∇Y f (Y)kF 6= 0

k∇Y f (Y)kF

1

.

(86)

(87)

otherwise

Note that 1 ≥ χ ≥ 0 because k∇Y f (Y)kF ≥ kPT⊥ Ye (∇Y f (Y))kF (see (85)). From (86) and
(87), we have

Y



kgradf (Yi )k2F = 1 −

2
e (∇Y f (Yi ))kF
Yi Y

kPT⊥

k∇Y f (Yi )k2F

≥ (1 − χ2 )k∇Y f (Yi )k2F .



 k∇Y f (Yi )k2F

(88)
(89)

Now, what remains is to show that k∇Y f (Yi )k2F ≥ 8f (Yi ). Indeed, from Lemma 2.7, we have
k∇Y f (Yi )k2F = keye((R + RT )1) − 2Rk2F ,

August 4, 2017

DRAFT

48

where R = PE (g(Yi )) − PE (D). Noting that R is symmetric with zero diagonal entries rii = 0,
we have
1
k∇Y f (Yi )k2F = keye(R1)k2F + kRk2F − 2 < eye(R1), R >
4
!
X X
= kR1k22 + kRk2F − 2
rij rii
i

j

= kR1k22 + kRk2F

≥ kRk2F
= 2f (Yi ).

(90)

By substituting (90) into (89), we obtain the desired result.
A PPENDIX I
P ROOF
Proof:
By denoting xi =

h

xi1 xi2 · · · xik

OF

iT

L EMMA 3.9

, we have

p = P (dij ≤ r)
= P (kxi − xj k22 ≤ r 2 )

k
X
= P ( (xit − xjt )2 ≤ r 2 ).
t=1

Before finding the general form of p, we compute the distribution of Y = (X1 − X2 )2 where
X1 and X2 are i.i.d. uniformly distributed random variables at unit interval. Let Z = X1 − X2 ,
then the cdf of Z is given by
FZ (z) = P (Z ≤ z)
= P (X1 − X2 ≤ z)
Z 1
=
P (X1 ≤ z + x2 X2 = x2 )fX2 (x2 )dx2
=
=

August 4, 2017

Z

Z

0

1

P (X1 ≤ z + x2 )fX2 (x2 )dx2

0
1

FX1 (z + x2 )fX2 (x2 )dx2 .
0

DRAFT

49

Thus,
d
FZ (z)
dz

1 − |z| if |z| ≤ 1
=
 0
otherwise

fZ (z) =

Now, we can easily obtain the cdf of Y = Z 2



1
if 1 ≤ y


√
FY (y) = 2 y − y if 0 ≤ y ≤ 1 .




0
if y ≤ 0

and also the pdf of Y


 √1 − 1 if 0 ≤ y ≤ 1
y
.
fY (y) =
 0
otherwise

Using this, we can compute p as

k
X

p = Pr

t=1

Z

=

···

yt ≤ r 2

Z

(91)

!

fY1 ,...,Yk (α1 , ..., αk )dα1 ...dαk

α1 +α2 +...+αk ≤r 2

Z

=

···

Z

fY1 (α1 )...fYk (αk )dα1 ...dαk ,

(92)

α1 +α2 +...+αk ≤r 2

where yt = (xit − xjt )2 . When the sensor nodes are located in two dimensional Euclidean space
(k = 2), we have

ZZ

p=

fY1 (α1 )fY2 (α2 )dα1 dα2 .

α1 +α2 ≤r 2

Let t = α1 + α2 , then we have

ZZ

p =

fY1 (α1 )fY2 (α2 )dα1 dα2

α1 +α2 ≤r 2

=

Z

r2

0

=

Z

0

August 4, 2017

Z

0

r2

1



fY1 (α1 )fY2 (t − α1 )dα1 dt

fY1 (t) ∗ fY2 (t)dt.

(93)

DRAFT

50

After some manipulations, we have

√


π
−
4
t+t
if 0 ≤ t ≤ 1


√
fY1 (t) ∗ fY2 (t) = 2 sin−1 ( 2t − 1) + 4 t − 1 − t − 2 if 1 ≤ t ≤ 2 .




0
otherwise

(94)

Ru
Ru
√
√
Let h1 (u) = (π − 4 t + t)dt (0 ≤ u ≤ 1) and h2 (u) = (2 sin−1 ( 2t − 1) + 4 t − 1 − t − 2)dt
1

0

(1 ≤ u ≤ 2), then we have

p=











h1 (r 2 )

if 0 ≤ r ≤ 1
√
h1 (1) + h2 (r 2 ) if 1 ≤ r ≤ 2 ,
1

otherwise

where h1 (u) and h2 (u) are given by
u

8
1
= πu − u3/2 + u2,
3
2
0
u


Z
u
2
1 2
8
−1
3/2
(t − 1) − t − 2t + 2 sin
− 1 dt
h2 (u) =
3
2
t
1
1
8
h1 (u) = πt − t3/2 + t2
3
2

1




2
sin
− 1 dt,
t
1


8
2
1
3/2
−1
=
(u − 1) − (u − 1)(u + 5) + 2u sin
−1
3
2
u
√
2
 + 3 u − 1 − π − 1.

+
2
sin−1 ( u −1)
1 + tan
2

1
8
(u − 1)3/2 − (u − 1)(u + 5) + 2
=
3
2

Zu

−1

Denoting p1 (r) = h1 (1) + h2 (r 2 ), we get the desired result for k = 2.

Similarly, when the sensor nodes are located in three dimensional Euclidean space (k = 3),
we have

ZZZ

p =

fY1 (α1 )fY2 (α2 )fY3 (α3 )dα1 dα2 dα3

α1 +α2 +α3 ≤r 2

=

Z

r2

0

=

Z

=

0

August 4, 2017

1

0

r2

0

Z

Z

Z

0

r2

Z

1

fY1 (α1 )fY2 (α2 )fY3 (t − α1 − α2 )dα1 dα2 dt

0

1

Z



1
0





fY1 (α1 )fY2 (u − α1 )dα1 fY3 (t − u)du dt

[fY1 (t) ∗ fY2 (t)] ∗ fY3 (t)dt.
DRAFT

51

After some manipulations, we have

e

h3 (t) if 0 ≤ t ≤ 1




e
h4 (t) if 1 ≤ t ≤ 2
,
[fY1 (t) ∗ fY2 (t)] ∗ fY3 (t) =

e

h
(t)
if
2
≤
t
≤
3
5




0
otherwise

where

(95)

√
√
1
e
h3 (t) = 4t t + 2π t − 3πt − t2 ,
2
√
√
√
1
1
e
h4 (t) = 3π − 4π t + (3π + 4)t + t2 − 8t t − 1 − 3 t − 1 + (t − 1)2
2
2
r
r


2
1
2
1
,

−4t sin−1
+ 4t sin−1 1 − − 2t sin−1
−1 −
t
t
t
sin−1 ( 2t −1)
1 + tan
2
r
√
√
√
1
t−2
−1
2
e
h5 (t) = 3 t − 2 + 4t t − 2 − 4t − (t − 1) − 3 + 8 t tan
2
t
r


√
t−2
3−t
−1
−1
−1
+ 2(t + 1) sin
t − 2 − 4(t − 1) sin
−8 tan
t−1
t−1
s
r
r
√
1
1
1
−8 t tan−1
+ 8 tan−1
+ 4(t − 1) sin−1
t(t − 2)
t−2
t−1
+
1 + tan
Now, we let h3 (u) =
Ru
2



Ru
0

2

sin−1 ( 3−t
t−1 )
2

.

e
h3 (t)dt (0 ≤ u ≤ 1), h4 (u) =

e
h5 (t)dt (2 ≤ u ≤ 3), then we have

August 4, 2017









Ru
1

e
h4 (t)dt (1 ≤ u ≤ 2), and h5 (u) =

h3 (r 2 )

if 0 ≤ r ≤ 1
√
h3 (1) + h4 (r 2 )
if 1 ≤ r ≤ 2
p=
√ ,
√


h3 (1) + h4 (2) + h5 (r 2 ) if 2 ≤ r ≤ 3




1
otherwise

(96)

DRAFT

52

where h3 (u), h4 (u), and h5 (u) are given by
4π √
8 √
3π 2 1 3
u u + u2 u −
u − u,
3
5
2
6
√
π 5
3π + 4 2 1 3 1
8π √
u u+
u + u + (u − 1)3 − − − 6 u − 1
h4 (u) = 3πu −
3
2
6
6
3 2
r
√
√
16
1
28
− (u − 1) u − 1 − (u − 1)2 u − 1 − 2u2 sin−1
3
5
u
r


2
16
1
2u2 sin−1 1 − − u2 sin−1
−1 −
3
u
u
3 1 + tan 21 sin−1 u2 − 1
4
+
2 ,
1 + tan 21 sin−1 u2 − 1
√
√
1
29 8
+ (u − 2)2 u − 2
h5 (u) = 2(u − 2) u − 2 − 2u2 − (u − 1)3 − 3u +
6
2 !5
√
√
√
√
22
16 2
+ (u − 2) u − 2 + 14 u − 2 +
− 12 π + 2 tan−1 u − 2
3
3
r
r
√
√
16
u
−
2
u−2
− 2u(u − 2) sin−1
−8u tan−1 u − 2 + u u tan−1
3
u
u−1
s
r


3−u
1
1
16 √
−1
−1
−1
+u(u + 2) sin
− u u tan
+ 8u tan
u−1
3
u(u − 2)
u−2
r
1
16
1
+2u(u − 2) sin−1
+

1
u−1
3 1 + tan sin−1 3−u 3

h3 (u) =

2

−4

1

1 + tan

1
2

sin−1

3−u
u−1

u−1

2 .

By denoting p2 (r) = h3 (1) + h4 (r 2 ) and p3 (r) = h3 (1) + h4 (2) + h5 (r 2 ), we get the desired
result for k = 3.

A PPENDIX J
P ROOF

OF

L EMMA 3.13

Proof: Since kAk2F = kPE (A)k2F + kPE⊥ (A)k2F , we can rewrite tkAk2F ≤ kPE (A)k2F as
tkAk2F ≤ kAk2F − kPE⊥ (A)k2F ,
and also
kPE⊥ (A)k2F ≤ (1 − t)kAk2F .
August 4, 2017

(97)
DRAFT

53

To show that (97) holds true with overwhelming probability, we first have
P (kPE⊥ (A)k2F ≥ (1 − t)kAk2F ) = P (exp(ǫkPE⊥ (A)k2F ) ≥ exp(ǫ(1 − t)kAk2F ))


(a)
2
⊥
2
exp(−ǫ(1
−
t)kAk
)E
exp(ǫkP
(A)k
)
F
E
F
≤
"
#
Y
(b) exp(−ǫ(1 − t)kAk2 )E
exp(ǫ(1 − δij )a2ij ) , (98)
F
=
i6=j

for any ǫ>0, where (a) follows from the Markov inequality and (b) is from kPE⊥ (A)k2F =
P
(1 − δij )a2ij (see (35)). Let Ω = {(i, j) : aij 6= 0} (i.e., Ω is the index set of nonzero entries

i6=j

e be a subset of Ω such that
of A), and N = 2⌊log2 |Ω|⌋ (|Ω| is the cardinality of Ω). Also, let Ω
e = N, then
|Ω|


"
#
Y
Y
Y
E
exp(ǫa2ij )E 
exp(ǫ(1 − δij )a2ij ) ≤
exp(ǫ(1 − δij )a2ij )
i6=j

e
(i,j)∈(Ω\Ω)

(a)
=

Y

e
(i,j)∈(Ω\Ω)

=

Y



exp(ǫa2ij ) 

=
=

e
(i,j)∈Ω

exp(ǫa2ij )

= exp(ǫkAk2F )

Y

e
(i,j)∈Ω

Y

e
(i,j)∈Ω

exp(ǫkAk2F )

Y

e
(i,j)∈Ω

exp(ǫa2ij )

(i,j)∈Ω

(b)
≤

Y

exp(ǫa2ij )

e
(i,j)∈(Ω\Ω)

Y

Y

e
(i,j)∈Ω

e
(i,j)∈(Ω\Ω)

Y

e
(i,j)∈Ω

Y

e
(i,j)∈Ω

1/N


E exp(Nǫ(1 − δij )a2ij ) 

(1 − p) exp(Nǫa2ij ) + p

1/N

exp(ǫa2ij ) 1 − p + p exp(−Nǫa2ij )

1 − p + p exp(−Nǫa2ij )

1 − p + p exp(−Nǫa2ij )

1/N

1/N

1/N

1 − p + p exp(−Nǫa2min )

1/N


(99)
exp(ǫkAk2F ) 1 − p + p exp(−Nǫa2min ) ,
M  M

Q
Q  M  1/M
where (a) is because E
Ai ≤
E Ai
for positive random variable Ai and M = 2q
(c)
=

i=1

i=1

e = N.
(q ≥ 1), (b) is because amin = min |aij |, and (c) is because |Ω|
(i,j)∈Ω

In summary, we have

P (kPE⊥ (A)k2F ≥ (1 − t)kAk2F ) ≤ g(ǫ),

August 4, 2017

DRAFT

54

where g(ǫ) = exp(mtNǫa2min ) (1 − p + p exp(−Nǫa2min )) (m = kAk2F /(a2minN)). If 0<mt<1,
we obtain the minimum value of g(ǫ) at ǫ∗ = 1/(Na2min ) log((1 − mt)p/((1 − p)mt)). Thus,

P (kPE⊥ (A)k2F ≥ (1 − t)kAk2F ) ≤ g(ǫ∗ )

1−mt  
1−p
p mt
=
1 − mt
mt

 
 

mt
1 − mt
+ mt log
,
= exp − (1 − mt) log
1−p
p
which is the desired result.
A PPENDIX K
P ROOF

OF

L EMMA 3.17

Proof: We first denote B = g(Y) − D, and then express B =
which gives
PE (B) =

X

P

ij

< B, ei eTj > ei eTj ,

δij < B, ei eTj > ei eTj .

i6=j

Recalling that l(A) = 2eye(Sym(A)1) − 2Sym(A), we have
∇Y f (Y) = l(PE (B))
X
= l(
δij < B, ei eTj > ei eTj )
i6=j

(a)
=

X

δij < B, ei eTj > l(ei eTj )

i6=j

=

X

δij sij l(ei eTj ),

i6=j

where sij =< B, ei eTj > and (a) is because l(αC + βD) = αl(C) + βl(D). Now, if we let
I = k∇Y f (Y)k2F − c2 kgradf (Y)k2F , then we have
XX
I =
δij δuv sij suv < l(ei eTj ), l(eu eTv ) >
i6=j u6=v

−c2

(a)
=

XX
i6=j u6=v

XX

δij δuv sij suv < PTY Ye (l(ei eTj )), PTY Ye (l(eu eTv )) >

δij δuv sij suv < l(ei eTj ), l(eu eTv ) >

i6=j u6=v

−c2
August 4, 2017

XX
i6=j u6=v

δij δuv sij suv < PTY Ye (l(ei eTj )), l(eu eTv ) >
DRAFT

55

=

XX
i6=j u6=v

≤
(b)
≤

=
=

XX
i6=j u6=v

δij δuv |sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |

1 XX 2
2
(δ + δuv
)|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |
2 i6=j u6=v ij

1 XX
(δij + δuv )|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |
2 i6=j u6=v
1X X
δij
|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |
2 i6=j
u6=v
+

(c)
=

δij δuv sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) >

1X X
|sij suv < (I − c2 PTY Ye )l(eu eTv ), l(ei eTj ) > |
δij
2 i6=j
u6=v

1X X
δij
|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |
2 i6=j
u6=v

1X X
δij
|sij suv < l(eu eTv ), (I − c2 PTY Ye )l(ei eTj ) > |
2 i6=j
u6=v
X X
|sij suv < (I − c2 PTY Ye )l(ei eTj ), l(eu eTv ) > |,
δij
=
+

u6=v

i6=j

where (a) is because < PTY Ye (E), F >=< PTY Ye (E), PTY Ye (F) > + < PTY Ye (E), PT⊥ Ye (F) >=<
Y

2

2

PTY Ye (E), PTY Ye (F) >, (b) is because x + y ≥ 2xy (x, y ≥ 0), and (c) is because (I − c2 PTY Ye )
is a self-adjoint operator.
A PPENDIX L
P ROOF

OF

L EMMA 3.18

Proof: Recalling that l(A) = 2eye(Sym(A)1) − 2Sym(A), we have
l(ei eTj ) = 2eye(Sym(ei eTj )1) − 2Sym(ei eTj )
= eye(ei eTj 1 + ej eTi 1) − ei eTj − ej eTi
(a)
=

eye(ei + ej ) − ei eTj − ej eTi

= ei eTi + ej eTj − ei eTj − ej eTi
= (ei − ej )(ei − ej )T ,

August 4, 2017

DRAFT

56

where (a) is because eye(ei eTj 1) = eye(ei ). Let δ = ei − ej , then l(ei eTj ) = δδ T . Also, if

n2 ≥ 4cµ(Y)2k 2 , then

| < (I − cPTY Ye )l(ei eTj ), l(ei eTj ) > |
= | < (I − cPTY Ye )δδ T , δδ T > |
= | < δδ T − cPTY Ye (δδ T ), δδ T > |
= | < δδ T , δδ T > −c < PTY Ye (δδ T ), δδ T > |
(a)
=

| < δδ T , δδ T > −c < PQ δδ T + δδ T PQ − PQ δδ T PQ , δδ T > |

= | < δδ T , δδ T > −c(< PQ δδ T , δδ T > + < δδ T PQ , δδ T > − < PQ δδ T PQ , δδ T >)|
(b)
=

|δ T δδ T δ − c(δ T PQ δδ T δ + δ T δδ T PQ δ − δ T PQ δδ T PQ δ)|

(c)
=

|kδk42 − c(2kPQδk22 kδk22 − kPQ δk42 )|

(d)
=

|4 − 4ckPQδk22 + ckPQ δk42 |.

4µ(Y)2r 2
16µ(Y)4r 4
4 − 4c
+
c
n2
n4


4cµ(Y)2r 2
≥ 4 1−
,
n2

(e)
≥

where (a) follows from Proposition 2.4, (b) is because < X, zzT >= tr(XzzT ) = zT Xz, (c) is
because PTQ PQ = PQ and thus δ T PQ δ = δ T PTQ PQ δ = kPQ δk22 , (d) is because kδk22 = 2 for
i 6= j, and (e) is because kPQ δk2 = kPQ ei − PQ ej k2 ≤ kPQ ei k2 + kPQ ej k2 ≤ 2µ(Y)r/n (see
Definition 3.14).
A PPENDIX M
P ROOF

OF

L EMMA 3.19

Proof: For any t>0, we have
P(

N
X
i=1

δi ai ≥ ǫ) = P (exp(t
(a)
≤
(b)
≤

August 4, 2017

N
X
i=1

δi ai ) ≥ exp(tǫ))

"

exp(−tǫ)E exp(t

exp(−tǫ + t

N
X
i=1

N
X

i=2n +1

δi ai )
"

#

ai )E exp(t

n

2
X
i=1

δi ai )

#
DRAFT

57

= exp(−tǫ + t

N
X

ai )E

i=2n +1
(c)
≤

exp(−tǫ + t

N
X

= exp(−tǫ + t

n

2
Y

ai )

(d)
=

exp(−tǫ + t

i=2n +1

= exp(−tǫ + t

N
X

ai )

= exp(−tǫ + t

= exp(−tǫ + t

i=2n +1

= exp(−tǫ + t

N
X

i=1

exp(−tǫ + t

n

2
Y

ai )

i=1

exp(t2n ai )((1 − p) exp(−t2n ai ) + p)

n

2
Y

ai )

exp(t2n ai )

i=1

2n

ai )

Y
i=1

! 21n

n

i=1

2n

Y
exp(tai )
((1 − p) exp(−t2n ai ) + p)
i=1

n

ai ) exp(t

2
X
i=1

n

i=1

n

i=1
n

i=1

! 21n

! 21n

2
Y
ai )
((1 − p) exp(−t2n amin ) + p)

i=1

! 21n

2
Y
ai )
((1 − p) exp(−t2n ai ) + p)

N
X

N
X

! 21n

2
Y
((1 − p) exp(−t2n ai ) + p)

2
Y
ai )
((1 − p) exp(−t2n ai ) + p)

i=1

= exp(−tǫ + t

! 21n

N
X
i=1

(e)
≤

! 21n

n

i=2n +1

= exp(−tǫ + t

! 21n

2
Y
ai )
(1 − p + p exp(t2n ai ))

i=2n +1
N
X


2n

E [exp(tδi ai 2n )]

i=1

i=2n +1
N
X

E (exp(tδi ai ))

n

2
Y

#



i=1

i=2n +1
N
X

exp(tδi ai )

i=1

i=2n +1
N
X

" 2n
Y

! 21n

ai )((1 − p) exp(−t2n amin ) + p),

where (a) follows from the Markov’s inequality, (b) is because δi ≤ 1 for all i, (c) is because
M
M
Q
Q
E[ Xi ] ≤ ( E[XiM ])1/M for positive random variables Xi and M = 2q (q ≥ 1), (d) is
i=1

i=1

because δi is Bernoulli random variable with P (δi = 1) = p, and (e) is because amin = min ai .
i
N
P
Let g(t) = exp(−tǫ+t ai )((1−p) exp(−t2n amin )+p), then the minimum of g(t) is obtained
i=1

August 4, 2017

DRAFT

! 21n

58

at t∗ = 1/(2n amin ) ln((1 − mǫ)(1 − p)/(mǫp)) (m = (
P(

N
X
i=1

δi ai ≥ ǫ) ≤ g(t∗ )

PN

i=1

ai − ǫ)/(2n ǫamin )). Thus, we have

1−mǫ 
mǫ
p
1−p
=
1 − mǫ
mǫ
 




mǫ
1 − mǫ
= exp − mǫ log
+ (1 − mǫ) log
1−p
p


when 0<mǫ<1 − p, which establishes the lemma.
R EFERENCES
[1]

L. Nguyen, S. Kim, and B. Shim, “Localization in internet of things network: Matrix completion approach,” Proc. Inform.
Theory Applicat. Workshop, 2016.

[2]

M. Delamo, S. Felici-Castell, J. J. Perez-Solano, and A. Foster, “Designing an open source maintenance-free environmental
monitoring application for wireless sensor networks,” J. Syst. Softw., vol. 103, pp. 238–247, May 2015.

[3]

G. Hackmann, W. Guo, G. Yan, Z. Sun, C. Lu, and S. Dyke, “Cyber-physical codesign of distributed structural health
monitoring with wireless sensor networks,” IEEE Trans. Parallel Distrib. Syst., vol. 25, pp. 63–72, Jan. 2014.

[4]

A. Pal, “Localization algorithms in wireless sensor networks: Current approaches and future challenges,” Netw. Protocols
Algorithms, vol. 2, no. 1, pp. 45–74, 2010.

[5]

V. J. Hodge, S. O’Keefe, M. Weeks, and A. Moulds, “Wireless sensor networks for condition monitoring in the railway
industry: A survey,” IEEE Trans. Intell. Transp. Syst., vol. 16, pp. 1088–1106, Jun. 2015.

[6]

P. Rawat, K. D. Singh, H. Chaouchi, and J. M. Bonnin, “Wireless sensor networks: a survey on recent developments and
potential synergies,” J. Supercomput., vol. 68, no. 1, pp. 1–48, Apr. 2014.

[7]

W. S. Torgerson, “Multidimensional scaling: I. theory and method,” Psychometrika, vol. 17, no. 4, pp. 401–419, 1952.

[8]

R. Parker and S. Valaee, “Vehicular node localization using received-signal-strength indicator,” IEEE Trans. Veh. Technol.,
vol. 56, pp. 3371–3380, Nov. 2007.

[9]

H. C. So, Y. T. Chan, and F. K. Chan, “Closed-form formulae for time-difference-of-arrival estimation,” IEEE Trans.
Signal Process., vol. 56, pp. 2614–2620, Jun. 2008.

[10]

Y. Shang, W. Ruml, Y. Zhang, and M. Fromherz, “Localization from mere connectivity,” in Proc. ACM Symp. Mobile Ad
Hoc Netw. Comput., Annapolis, Maryland, USA, Jun. 2003, pp. 201–212.

[11]

M. Fazel, “Matrix rank minimization with applications,” Ph.D. dissertation, Standford Univ., Standford, CA, 2002.

[12]

E. J. Candes and B. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math., vol. 6, pp. 717–772,
2009.

[13]

J. F. Cai, E. J. Candes, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM J. Optimiz.,
vol. 20, no. 4, pp. 1956–1982, 2010.

[14]

K. Lee and Y. Bresler, “Admira: Atomic decomposition for minimum rank approximation,” IEEE Trans. Inf. Theory,
vol. 56, no. 9, pp. 4402–4416, 2010.

August 4, 2017

DRAFT

59

[15]

Z. Wang, M. J. Lai, Z. Lu, W. Fan, H. Davulcu, and J. Ye, “Orthogonal rank-one matrix pursuit for low rank matrix
completion,” SIAM J. Sci. Comput., no. 1, pp. A488–A514, 2015.

[16]

Z. Wen, W. Yin, and Y. Zhang, “Solving a low-rank factorization model for matrix completion by a nonlinear successive
over-relaxation algorithm,” Math. Program. Comput., no. 4, pp. 333–361, 2012.

[17]

T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh, “Matrix completion and low-rank svd via fast alternating least squares,”
J. Mach. Learning Research, no. 16, pp. 3367–3402, 2015.

[18]

B. Vandereycken, “Low-rank matrix completion by riemannian optimization,” SIAM J. Optimiz., vol. 23, no. 2, pp. 1214–
1236, 2013.

[19]

P. A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds. Princeton Univ. Press, 2008.

[20]

J. Lee, Introduction to Smooth Manifolds, 2nd ed. New York: NY:Springer, 2013, vol. 218.

[21]

U. Helmke and J. B. Moore, Optimization and Dynamical Systems. London: Springer-Verlag, 1994.

[22]

B. Mishra, G. Meyer, S. Bonnabel, and R. Sepulchre, “Fixed-rank matrix factorizations and riemannian low-rank
optimization,” Comput. Stat., no. 3-4, pp. 591–621, 2014.

[23]

P. A. Absil and J. Malick, “Projection-like retractions on matrix manifolds,” SIAM J. Optimiz., vol. 22, pp. 135–158,
2012.

[24]

Y. H. Dai, “Nonlinear conjugate gradient methods,” Wiley Encyclopedia of Operations Research and Manage. Sci., 2011.

[25]

L. Armijo, “Minimization of functions having lipschitz continuous first partial derivatives,” Pacific J. Math., no. 1, 1966.

[26]

C. T. Kelley, Iterative methods for optimization. Philadelphia: PA: Frontiers in Applied Mathematics, 1999.

[27]

M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” NBS, 1952.

[28]

R. Fletcher and C. M. Reeves, “Function minimization by conjugate gradients,” J. Comput., no. 2, pp. 149–154, 1964.

[29]

Y. H. Dai and Y. Yuan, “A nonlinear conjugate gradient method with a strong global convergence property,” SIAM J.
Optimiz., no. 1, pp. 177–182, 1999.

[30]

W. W. Hager and H. Zhang, “A new conjugate gradient method with guaranteed descent and an efficient line search,”
SIAM J. Optimiz., no. 1, pp. 170–192, 2005.

[31]

Y. Hu, D. Zhang, J. Ye, X. Li, and X. He, “Fast and accurate matrix completion via truncated nuclear norm regularization,”
IEEE Trans. Pattern Anal. Mach. Intell., no. 9, pp. 2117–2130, Sep. 2013.

[32]

P. Wolfe, “Convergence conditions for ascent methods,” SIAM Rev., no. 2, pp. 226–235, 1969.

[33]

H. Sato and T. Iwai, “A new, globally convergent riemannian conjugate gradient method,” Optim. J. Math. Program. Oper.
Res., no. 4, pp. 1011–1031, 2015.

[34]

E. J. Candes and T. Tao, “The power of convex relaxation: Near-optimal matrix completion,” IEEE Trans. Inf. Theory,
no. 5, pp. 2053–2080, May 2010.

[35]

B. Recht, “A simpler approach to matrix completion,” J. Mach. Learning Research, vol. 12, pp. 3413–3430, Dec. 2011.

[36]

W. Cheney, Analysis for applied mathematics. New York: Springer, 2013.

[37]

V. Hutson, J. Pym, and M. Cloud, Applications of functional analysis and operator theory. Elsevier, 2005.

August 4, 2017

DRAFT

