Localization
Paul Tseng
Mathematics, University of Washington Seattle
7th Int. Conf. Num. Optim. Num. Lin. Algeb., Lijiang August 9, 2009
(joint work with Ting Kei Pong)
Talk Outline
• Sensor network localization
• SDP, ESDP relaxations: properties and soln accuracy certificate
Talk Outline
• Sensor network localization
• SDP, ESDP relaxations: properties and soln accuracy certificate
• A robust version of ESDP to handle noises
• Log-barrier penalty CGD method
Talk Outline
• Sensor network localization
• SDP, ESDP relaxations: properties and soln accuracy certificate
• A robust version of ESDP to handle noises
• Log-barrier penalty CGD method
• Numerical simulations
• Conclusion & Ongoing work
Sensor Network Localization
Basic Problem
:• n pts in <2.
• Know last n − m pts (‘anchors’) xm+1, ..., xn and Eucl. dist. estimate for pairs of ‘neighboring’ pts
dij ≥ 0 ∀(i, j) ∈ A with A ⊆ {(i, j) : 1 ≤ i, j ≤ n}.
• Estimate first m pts (‘sensors’).
Sensor Network Localization
Basic Problem
:• n pts in <2.
• Know last n − m pts (‘anchors’) xm+1, ..., xn and Eucl. dist. estimate for pairs of ‘neighboring’ pts
dij ≥ 0 ∀(i, j) ∈ A with A ⊆ {(i, j) : 1 ≤ i, j ≤ n}.
• Estimate first m pts (‘sensors’).
History? Graph realization/rigidty, Euclidean matrix completion, position estimation in wireless sensor network, ...
Optimization Problem Formulation
υopt := min
x1,...,xm
X
(i,j)∈A
kxi − xjk2 − d2ij
Optimization Problem Formulation
υopt := min
x1,...,xm
X
(i,j)∈A
kxi − xjk2 − d2ij
• Objective function is nonconvex. m can be large (m ≥ 1000).
6. .
_
• Problem is NP-hard (reduction from PARTITION).
6. .
_
• Local improvement heuristics can fail badly.
6. .
_
Optimization Problem Formulation
υopt := min
x1,...,xm
X
(i,j)∈A
kxi − xjk2 − d2ij
• Objective function is nonconvex. m can be large (m ≥ 1000).
6. .
_
• Problem is NP-hard (reduction from PARTITION).
6. .
_
• Local improvement heuristics can fail badly.
6. .
_
• Use a convex (SDP, SOCP) relaxation (& local improvement).
Low soln accuracy OK. Distributed computation preferred.
SDP Relaxation
Let X := [x1 · · · xm]. Y = XTX ⇐⇒ Z = Y XT
X I
0, rankZ = 2
SDP Relaxation
Let X := [x1 · · · xm]. Y = XTX ⇐⇒ Z = Y XT
X I
0, rankZ = 2 SDP relaxation (Biswas,Ye ’03):
υsdp := min
Z
X
(i,j)∈A,i≤m<j
yii − 2xTj xi + kxjk2 − d2ij
+ X
(i,j)∈A,i<j≤m
yii − 2yij + yjj − d2ij s.t. Z = Y XT
X I
0
Adding the nonconvex constraint rankZ = 2 yields original problem.
But SDP relaxation is still expensive to solve for m large..
ESDP Relaxation
ESDP relaxation (Wang, Zheng, Boyd, Ye ’06):
υesdp := min
Z
X
(i,j)∈A,i≤m<j
yii − 2xTj xi + kxjk2 − d2ij
+ X
(i,j)∈A,i<j≤m
yii − 2yij + yjj − d2ij s.t. Z = Y XT
X I
yii yij xTi yij yjj xTj xi xj I
0 ∀(i, j) ∈ A, i < j ≤ m
0 ≤ υesdp ≤ υsdp ≤ υopt. In simulation, ESDP is nearly as strong as SDP, and solvable much faster by IP method.
Example 1
n = 3, m = 1, d12 = d13 = 2
Problem
:0 = min
x1∈<2
|kx1 − 1
0 k2 − 4| + |kx1 − −1
0 k2 − 4|
SDP/ESDP Relaxation
:0 = min
x1=[αβ]∈<2
y11∈<
|y11 − 2α − 3| + |y11 + 2α − 3|
s.t.
y11 α β
α 1 0
β 0 1
0
If solve SDP/ESDP by IP method, then likely get analy. center y11 = 3, x1 = h
0 0
i
Example 2
n = 4, m = 1, d12 = d13 = 2, d14 = 1
Problem
:0 = min
x1∈<2
|kx1 − 1
0 k2 − 4| + |kx1 − −1
0 k2 − 4| + |kx1 − h
√1 3
ik2 − 1|
SDP/ESDP Relaxation
:0 = min
x1=[αβ]∈<2
y11∈<
|y11 − 2α − 3| + |y11 + 2α − 3| + |y11 − 2α − 2√
3β + 3|
s.t.
y11 α β
α 1 0
β 0 1
0
SDP/ESDP has unique soln y11 = 3, x1 = h
√0 3
i
Properties of SDP & ESDP Relaxations
Assume each i ≤ m is conn. to some j > m in the graph ({1, ..., n}, A).
Fact 0
:• Sol(SDP) and Sol(ESDP) are nonempty, closed, convex.
• If
dij = kxtruei − xtruej k ∀ (i, j) ∈ A “noiseless case”
(xtruei = xi ∀ i > m), then
υopt = υsdp = υesdp = 0 and
Ztrue := Xtrue I T
Xtrue I
is a soln of SDP and ESDP (i.e., Ztrue ∈ Sol(SDP) ⊆ Sol(ESDP)).
Let tri[Z] := yii − kxik2, i = 1, ..., m. “ith trace”
Fact 1
(Biswas,Ye ’03, T ’07, Wang et al ’06): For each i,tri[Z] = 0 ∃Z ∈ ri(Sol(ESDP)) =⇒ xi is invariant over Sol(ESDP) (so xi = xtruei in noiseless case) Still true with “ESDP” changed to “SDP”.
Let tri[Z] := yii − kxik2, i = 1, ..., m. “ith trace”
Fact 1
(Biswas,Ye ’03, T ’07, Wang et al ’06): For each i,tri[Z] = 0 ∃Z ∈ ri(Sol(ESDP)) =⇒ xi is invariant over Sol(ESDP) (so xi = xtruei in noiseless case) Still true with “ESDP” changed to “SDP”.
Fact 2
(Pong, T ’09): Suppose υopt = 0. For each i,tri[Z] = 0 ∀Z ∈ Sol(ESDP) ⇐= xi is invariant over Sol(ESDP).
Proof is by induction, starting from sensors that neighbor anchors.
(Q: True for SDP?)
Proof idea
:• If (i, j) ∈ A and xi, xj are invar. over Sol(ESDP), then tri[Z] = trj[Z]
∀Z ∈ Sol(ESDP).
• Suppose ∃i ≤ m such that xi is invar. over Sol(ESDP) but tri[ ¯Z] > 0 for some ¯Z ∈ Sol(ESDP). Consider maximal ¯I ⊂ {1, . . . , m} such that xi is invar. over Sol(ESDP) and tri[ ¯Z] > 0 ∀i ∈ ¯I.
• Then xi is not invar. over Sol(ESDP) ∀i ∈ N (¯I).
So ∃Z ∈ ri(Sol(ESDP)) with xi 6= ¯xi ∀i ∈ N (¯I).
• Let Zα = α ¯Z + (1 − α)Z with α > 0 suff. small.
Can rotate xαi ∀i ∈ ¯I and Zα still remains in Sol(ESDP). ⇒⇐
In practice, there are measurement noises:
d2ij = kxtruei − xtruej k2 + δij ∀(i, j) ∈ A.
When δ := (δij)(i,j)∈A ≈ 0, does tri[Z] = 0 (with Z ∈ ri(Sol(ESDP))) imply xi ≈ xtruei ?
In practice, there are measurement noises:
d2ij = kxtruei − xtruej k2 + δij ∀(i, j) ∈ A.
When δ := (δij)(i,j)∈A ≈ 0, does tri[Z] = 0 (with Z ∈ ri(Sol(ESDP))) imply xi ≈ xtruei ? No!
6. .
_
Fact 3
(Pong, T ’09): For δ ≈ 0 and for each i,tri[Z] = 0 ∃Z ∈ ri(Sol(ESDP)) 6=⇒ xi ≈ xtruei . Still true with “ESDP” changed to “SDP”.
Proof is by counter-example.
An example of sensitivity of ESDP solns to measurement noise:
Problem data: m = 2, n = 6;
d12 = p4 + (1 − )2, d13 = 1 + , d14 = 1 − , d25 = d26 = √
2 ( > 0)
Thus, even when Z ∈ Sol(ESDP) is unique, tri[Z] = 0 fails to certify accuracy of xi in the noisy case!
Robust ESDP
Fix any ρij > |δij| ∀(i, j) ∈ A (ρ > |δ|).
Let Sol(ρESDP) denote the set of Z = Y XT
X I
satisfying
|yii − 2xTj xi + kxjk2 − d2ij| ≤ ρij ∀(i, j) ∈ A, i ≤ m < j
|yii − 2yij + yjj − d2ij| ≤ ρij ∀(i, j) ∈ A, i < j ≤ m
yii yij xTi yij yjj xTj xi xj I
0 ∀(i, j) ∈ A, i < j ≤ m
Note: Ztrue = Xtrue I T
Xtrue I ∈ Sol(ρESDP).
Let
Zρ,δ := arg min
Z∈Sol(ρESDP)
X
(i,j)∈A,i<j≤m
− ln det
yii yij xTi yij yjj xTj xi xj I
Let
Zρ,δ := arg min
Z∈Sol(ρESDP)
X
(i,j)∈A,i<j≤m
− ln det
yii yij xTi yij yjj xTj xi xj I
Fact 4
(Pong, T ’09): ∃ η > 0 and ¯ρ > 0 such that for each i, tri[Zρ,δ] < η ∃|δ| < ρ ≤ ¯ρe =⇒ lim|δ|<ρ→0xρ,δi = xtruei
tri[Zρ,δ] > 10η ∃|δ| < ρ ≤ ¯ρe =⇒ xi not invar. over Sol(ESDP) when δ = 0 Moreover,
kxρ,δi − xtruei k ≤ p
2|A| + m q
tri[Zρ,δ] ∀ |δ| < ρ.
Log-barrier Penalty CGD Method
Efficiently compute Zρ,δ? Let ha(t) := 1
2(t − a)2+ + 1
2(−t − a)2+ (|t| ≤ a ⇐⇒ ha(t) = 0) and
fµ(Z) := X
(i,j)∈A,i≤m<j
hρij(yii − 2xTj xi + kxjk2 − d2ij)
+ X
(i,j)∈A,i<j≤m
hρij(yii − 2yij + yjj − d2ij)
+µ X
(i,j)∈A,i<j≤m
− ln det
yii yij xTi yij yjj xTj xi xj I
• fµ is partially separable, strictly convex & diff. on its domain.
• For each fixed ρ > |δ|, argminfµ → Zρ,δ as µ → 0.
• fµ is partially separable, strictly convex & diff. on its domain.
• For each fixed ρ > |δ|, argminfµ → Zρ,δ as µ → 0.
Idea
: Minimize fµ approx. by block-coordinate gradient descent (BCGD). (T,Yun ’06)
Log-barrier Penalty CGD Method
:Given Z in domfµ, compute gradient ∇Zifµ of fµ w.r.t.
Zi := {xi, yii, yij : (i, j) ∈ A} for each i.
• If k∇Zifµk ≥ max{µ, 10−7} for some i, update Zi by moving along the Newton direction −
∂Z2
iZifµ−1
∇Zifµ with Armijo stepsize rule.
• Decrease µ when k∇Zifµk < max{µ, 10−6} ∀ i.
µinitial = 10, µfinal = 10−14. Decrease µ by a factor of 10 each time.
Coded in Fortran. Compute Newton direc. by sparse Cholesky.
Computation easily distributes.
Simulation Results
• Compare ρESDP as solved by LPCGD method with ESDP as solved by Sedumi 1.05 Sturm (with the interface to Sedumi coded by Wang et al).
Simulation Results
• Compare ρESDP as solved by LPCGD method with ESDP as solved by Sedumi 1.05 Sturm (with the interface to Sedumi coded by Wang et al).
• Anchors and sensors xtrue1 , ..., xtruen uniformly distributed in [−.5, .5]2, m = .9n. (i, j) ∈ A whenever kxtruei − xtruej k ≤ rr. Set
dij = kxtruei − xtruej k · |1 + σ · ij|, where ij ∼ N (0, 1).
Simulation Results
• Compare ρESDP as solved by LPCGD method with ESDP as solved by Sedumi 1.05 Sturm (with the interface to Sedumi coded by Wang et al).
• Anchors and sensors xtrue1 , ..., xtruen uniformly distributed in [−.5, .5]2, m = .9n. (i, j) ∈ A whenever kxtruei − xtruej k ≤ rr. Set
dij = kxtruei − xtruej k · |1 + σ · ij|, where ij ∼ N (0, 1).
• Sensor i is judged as “accurately positioned” if
tri[Zfound] < (.01 + 30σ)davgij .
ρESDPLPCGD ESDPSedumi
n m σ rr cpu/map/errap cpu(cpus)/map/errap 1000 900 0 .06 7/662/1.7e-3 182(104)/669/2.1e-3 1000 900 .01 .06 5/660/2.2e-2 119(42)/720/3.1e-2 2000 1800 0 .06 26/1762/3.1e-4 1157(397)/1742/3.9e-4 2000 1800 .01 .06 20/1699/1.4e-2 966(233)/1746/2.4e-2 10000 9000 0 .02 77/7844/2.3e-3 16411(1297)/6481/2.5e-3 10000 9000 .01 .02 63/8336/1.0e-2 16368(1264)/8593/8.7e-3
• cpu(sec) times are on a HP DL360 workstation, running Linux 3.5. ESDP is solved by Sedumi; cpus:= run time for Sedumi.
• Set ρij = d2ij · ((1 − 2σ)−2 − 1).
• map := # accurately positioned sensors.
errap := maxiaccurate. pos. kxi − xtruei k.
900 sensors, 100 anchors, rr = 0.06, σ = 0.01, solve ρESDP by LPCGD method. xtruei (shown as ∗) and xρ,δi (shown as •) are joined by blue line segment; anchors are shown as ◦.
60 sensors, 4 anchors at corners, rr = 0.3, σ = 0.1. xi (shown as ∗) and xρ,δi (shown as
•) are joined by blue line segment; anchors are shown as ◦. Left: Soln of ρESDP found by LPCGD method. Right: After local gradient improvement.
−0.5 0 0.5
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4 0.5
−0.5 0 0.5
−0.5
−0.4
−0.3
−0.2
−0.1 0 0.1 0.2 0.3 0.4 0.5
Conclusion & Ongoing work
• SDP and ESDP solns are sensitive to measurement noise. Has soln accuracy certificate under no noise only (though it works well enough in simulation).
• ρESDP solns are more stable. Has soln accuracy certificate under low noise (which works well enough in simulation). Needs to estimate the noise level δ to set ρ. Can ρ > |δ| be
relaxed?
• SDP, ESDP, ρESDP solns can be further refined by local improvement. This improves the rmsd when noise level is high (e.g., σ = 0.1).
• Approximation bounds? Extensions to handle lower bounds on distances (e.g., (i, j) 6∈ A imply kxtruei − xtruej k > rr)?
Thanks for coming!
6. .
^