A Factor Graph Approach to Clock Offset Estimation in Wireless Sensor Networks

0.5 MB·English

by  Aitzaz Ahmad

1
A Factor Graph Approach to Clock Offset
Estimation in Wireless Sensor Networks
Aitzaz Ahmad, Davide Zennaro, Erchin Serpedin, and Lorenzo Vangelista
Abstract—The problem of clock offset estimation in a two for timing synchronization is presented in [2] and [3]. The
way timing message exchange regime is considered when the one-way message exchange mechanism involves a reference
likelihood function of the observation time stamps is Gaussian,
node broadcasting its timing information to other nodes in
exponentialorlog-normallydistributed.Aparametrizedsolution
a network. The receiver nodes record the arrival of these
tothemaximumlikelihood(ML)estimationofclockoffset,based
2 on convex optimization, is presented, which differs from the messages with respect to their own clock. After several such
1 earlier approaches where the likelihood function is maximized time stamps have been exchanged, the nodes estimate their
0 graphically. In order to capture the imperfections in node offsets based on these observations. A particular case of this
2 oscillators,whichmayrenderatime-varyingnaturetotheclock approachisthefloodingtimesynchronizationprotocol(FTSP)
offset,anovelBayesianapproachtotheclockoffsetestimationis
n [4] which uses regression to estimate the clock offset. On
proposedbyusingafactorgraphrepresentationoftheposterior
a density.Messagepassingusingthemax-productalgorithmyields the other hand, through a two-way timing exchange process,
J
a closed form expression for the Bayesian inference problem. adjacent nodes aim to achieve pairwise synchronization by
1 Several lower bounds on the variance of an estimator are communicatingtheirtiminginformationwitheachother.After
3 derived for arbitrary exponential family distributed likelihood
a round of N messages, each nodes tries to estimate its
functions which, while serving as stepping stones to benchmark
] the performance of the proposed clock offset estimators, can be own clock parameters. The timing-sync protocol for sensor
T useful in their own right in classical as well Bayesian param- networks (TPSNs) [5] uses this strategy in two phases to
I eter estimation theory. To corroborate the theoretical findings, synchronize clocks in a network. The level discovery phase
.
s extensive simulation results are discussed for classical as well involves a spanning tree based representation of a WSN
c
as Bayesian estimators in various scenarios. It is observed that
[ while nodes attempt to synchronize with their immediate
the performance of the proposed estimators is fairly close to the
1 fundamental limits established by the lower bounds. parents using a two-way message exchange process in the
synchronization phase. In receiver-receiver synchronization,
v
Index Terms—Clock synchronization, factor graphs, message
5 nodes collect time stamps sent from a common broadcasting
passing, estimation bounds, wireless sensor networks.
8 node and utilize them to adjust their clocks. The reference
6 broadcast synchronization (RBS) protocol [6] uses reference
6 I. INTRODUCTION beacons sent from a master node to establish a common
.
1 WIRELESS sensor networks (WSNs) typically consist notion of time across a network. An alternative framework
0 of a large number of geographically distributed sensor for network-wide distributed clock synchronization consists
2
nodes,deployedtoobservesomephenomenonofinterest.The of recasting the problem of agreement on oscillation phases
1
: nodes constituting such a network are low cost sensors that and/or frequencies as a consensus based recursive model in
v have limited abilities of data processing and communication. which only local message passing is required among nodes.
i
X WSNs envisage tremendous applications in such diverse areas By assuming a connected network, it is possible to design ef-
r as industrial process control, battlefield surveillance, health ficient distributedalgorithms bycarefully choosing theupdate
a monitoring, target localization and tracking, etc., [1]. With function. Under this framework, [9] proposed a Laplacian-
the recent advances in digital circuit technology, WSNs are based algorithm for establishing agreement on oscillation
expected to play a pivotal role in future wireless communica- frequenciesalloverthenetworkbasedonstandardconsensus.
tions. Acombinedagreementoverbothclockphasesandfrequencies
Clock synchronization in sensor networks is a critical hasbeenstudiedin[10],bymakinguseofstate-of-the-artfast
component in data fusion and duty cycling operations, and consensustechniques.Scalablesynchronizationalgorithmsfor
has gained widespread interest in the past few years. Most large sensor networks are developed in [11] and [12] inspired
of the current methods consider sensor networks exchanging by mathematical biology models justifying synchrony in the
time stamps based on the time at their respective clocks. biological agents.
A survey of the popular approaches employed in practice The clock synchronization problem in a WSN offers a
natural statistical signal processing framework whereby, the
A. Ahmad and E. Serpedin are with the Department of Electrical and clockparametersaretobeestimatedusingtiminginformation
ComputerEngineering,TexasA&MUniversity,Texas,TX,77843USA.
from various sensors [22]. A model based synchronization
D. Zennaro and L. Vangelista are with Department of Information Engi-
neering,UniversityofPadova,Padova,Italy. approach to arrest the clock drifts is explored in [8]. The
A.AhmadandE.SerpedinaresupportedbyQtel. impairments in message transmission arise from the various
TheworkofD.Zennaroispartiallysupportedbyan“A.Gini”fellowship
delays experienced by the messages as they travel through
and has been performed while on leave at Texas A&M University, College
Station,TX(USA). the transmission medium. Therefore, a crucial component of
2
efficient clock parameter estimation is accurate modeling of Gauss-Markovprocess.Bayesianinferenceisperformed
thenetworkdelaydistributions.Severaldistributionshavebeen using factor graphs and the max-product algorithm. The
proposed that aim to capture the random queuing delays in message passing strategy yields a closed form solution
a network [7]. Some of these candidate distributions include for Gaussian, exponential and log-normally distributed
exponential, Weibull, Gamma and log-normal distributions. likelihood functions. This extends the current literature
Assuminganexponentialdelaydistribution,severalestimators tocaseswheretheclockoffsetmaynotbedeterministic,
wereproposedin[13].Itwasarguedthatwhenthepropagation but is in fact a random process.
delay d is unknown, the maximum likelihood (ML) estimator 3) In order to evaluate the performance of the proposed
for the clock offset θ is not unique. However, it was later estimators, classical as well as Bayesian bounds are
shown in [14] that the ML estimator of θ does exist uniquely derived for arbitrary exponential family distributed like-
for the case of unknown d. The performance of these estima- lihood functions, which is a wide class and contains
torswascomparedwithbenchmarkestimationboundsin[15]. almost all distributions of interest. While these results
Considering an offset and skew model, Chaudhari et.al. pre- aid in comparing various estimators in this work, they
sented algorithms for the joint ML estimation of clock offset canbeusefulintheirownrightinclassicalandBayesian
and skew in [16] when the network delays are exponentially estimation theory.
distributed.Clockoffsetandskewestimatorsweredetermined
This paper is organized as follows. The system model is
in [17] based on the assumption that the network delays
outlined in Section II. The ML estimation of clock offset
arise from the contribution of several independent processes
based on convex optimization is proposed in Section III. The
and as such, were modeled as Gaussian. The convergence
factorgraphbasedinferencealgorithmforthesynchronization
of distributed consensus time synchronization algorithms is
problem in a Bayesian paradigm is detailed in Section IV and
investigated in [18] and [19], assuming a Gaussian delay
a closed form solution is obtained. Section V presents several
between sensor nodes. More recently, the minimum variance
theoretical lower bounds on the variance of an estimator eval-
unbiased estimator (MVUE) for the clock offset under an
uated in the classical as well as Bayesian regime. Simulation
exponential delay model was proposed in [20]. The timing
studies are discussed in Section VI which corroborate the
synchronization problem for the offset-only case was also
earlier results. Finally, the paper is concluded in Section VII
recastasaninstanceofconvexoptimizationin[21]forWeibull
along with some directions for future research.
distributed network delays. A recent contribution [23] has
investigatedthefeasibilityofdeterminingtheclockparameters
by studying the fundamental limits on clock synchronization II. SYSTEMMODEL
for wireline and wireless networks.
Theprocessofpairwisesynchronizationbetweentwonodes
In this work, considering the classic two-way message
S and R is illustrated in Fig. 1. At the jth message exchange,
exchangemechanism,aunifiedframeworkfortheclockoffset
node S sends the information about its current time through
estimation problem is presented when the likelihood function
a message including time stamp T1. Upon receipt of this
of the observations is Gaussian, exponential or log-normally j
message, Node R records the reception time T2 according to
distributed. A parameterized solution is proposed for the j
itsowntimescale.Thetwo-waytimingmessageexchangepro-
ML estimation of clock offset by recasting the likelihood
cessiscompletedwhennodeRreplieswithasynchronization
maximization as an instance of convex optimization. In order
packet containing time stamps T2 and T3 which is received
to incorporate the effect of time variations in the clock j j
at time T4 by node S with respect to its own clock. After N
offset between sensor nodes, a Bayesian inference approach j
such messages have been exchanged between nodes S and R,
is also studied based on a factor graph representation of the
node S is equipped with time stamps T1,T2,T3,T4 N .
posterior density. The major contributions of this work can be { j j j j}j=1
The impairments in the signaling mechanism occur due to a
summarized as follows.
fixed propagation delay, which accounts for the time required
1) A unified framework for ML estimation of clock offset, by the message to travel through the transmission medium,
based on convex optimization, is presented when the and a variable network delay, that arises due to queuing
likelihood function of the observations is Gaussian, delay experienced by the messages during transmission and
exponential and log-normally distributed. The proposed reception. By assuming that the respective clocks of nodes S
frameworkrecoversthealreadyknownresultsforGaus- and R are related by CR(t)=θ+CS(t), the two-way timing
sian and exponentially distributed likelihood functions message exchange model at the jth instant can be represented
and determines the ML estimate in case of log-normal as
distribution. Hence, the proposed convex optimization
based approach represents a simpler alternative, and a Tj2 =Tj1+d+θ+Xj
more general derivation of ML estimator, which by- T4 =T3+d θ+Y (1)
j j − j
passes the graphical analysis used in [14] to maximize
the likelihood function. where d represents the propagation delay, assumed symmetric
2) In order to capture the time variations in clock offsets in both directions, and θ is offset of the clock at node R
due to imperfect oscillators, a Bayesian framework is relative to the clock at node S. X and Y are the indepen-
j j
presented by considering the clock offset as a random dent and identically distributed variable network delays. By
3
T2 T3 V =∆ [V1,...,VN]T is Gaussian or log-normally distributed
R j j
is given below.
Unconstrained Likelihood:
 
N
(cid:88)
f(U;ξ) expξ ηξ(Uj) Nφξ(ξ) (5)
∝ −
j=1
S Tj1 Tj4  (cid:88)N 
f(V;ψ) expψ ηψ(Vj) Nφψ(ψ) (6)
d+θ+X d θ+Y ∝ −
j j
− j=1
whereη (U )andη (V )aresufficientstatisticsforestimating
ξ j ψ j
Fig.1. Atwo-waytimingmessageexchangemechanism ξ and ψ, respectively. The log-partition functions φ (.) and
ξ
φ (.) serve as normalization factors so that f(U;ξ) and
ψ
f(V;ψ) are valid probability distributions. The likelihood
defining [13] function is called ‘unconstrained’ since its domain is inde-
pendent of the parameters ξ and ψ.
Uj =∆ Tj2−Tj1 Similarly, the general notation used for an exponentially
V =∆ T4 T3 , distributed likelihood function is given below.
j j − j
the system in (1) can be equivalently expressed as Constrained Likelihood:
 
N N
U =d+θ+X (cid:88) (cid:89)
j j f(U;ξ) expξ ηξ(Uj) Nφξ(ξ) I(Uj ξ)
V =d θ+Y . (2) ∝ − −
j − j j=1 j=1
(7)
By further defining
 
N N
ξ =∆ d+θ ψ =∆ d−θ , (3) f(V;ψ)∝expψ(cid:88)ηψ(Vj)−Nφψ(ψ)(cid:89)I(Vj −ψ)
j=1 j=1
the model in (2) can be written as
(8)
Uj =ξ+Xj where the indicator function I(.) is defined as
V =ψ+Y (cid:26)
j j 1 x 0
I(x)= ≥ .
0 x<0
for j = 1,...,N. The goal is to determine precise estimates
of ξ and ψ using observations {Uj,Vj}Nj=1. An estimate of θ andtherolesofηξ(Uj),ηψ(Vj),φξ(.)andφψ(.)aresimilarto
can, in turn, be obtained using (3) as follows (5)and(6).Thelikelihoodfunctioniscalledconstrainedsince
itsdomaindependsontheparametersξandψ.Itmustbenoted
ξ ψ
θ = − . (4) that the likelihood functions (5)-(8) are expressed in terms of
2
general exponential family distributions. This approach helps
Accurate modeling of the variable delays, X and Y , has
j j to keep the exposition sufficiently general and also allows us
been a topic of interest in recent years. Several distributions
to recover the known results for the ML estimation of clock
have been proposed that aim to capture the random effects
offset for Gaussian and exponentially distributed likelihood
caused by the queuing delays [7]. These distributions include
functions [14] [15], and determine the ML estimator of the
exponential, gamma, log-normal and Weibull. In addition, the
clock offset in case of log-normally distributed likelihood
authorsin[17]arguedthatX andY resultfromcontributions
j j function, as shown in Section III. The proposed approach will
ofnumerousindependentrandomprocessesandcan,therefore,
alsoproveusefulininvestigatingaunifiednovelframeworkfor
be assumed to be Gaussian. The ML estimate of d and θ for
clock offset estimation in the Bayesian setting for Gaussian,
the case of exponential distribution was determined in [14].
exponential or log-normally distributed likelihood functions,
Recently, the minimum variance unbiased estimate (MVUE)
as will be shown in Section IV.
of the clock offset under an exponentially distributed network
Some key ingredients of the proposed solution for the
delay was proposed in [20]. In this work, instead of working
clock offset estimation problem, based on the properties of
with a specific distribution, a general framework of the clock
exponential family, can be summarized as follows [25].
synchronization problem is proposed that yields a parameter-
1) The mean and variance of the sufficient statistic η (U )
ξ j
ized solution of the clock offset estimation problem in the
are expressed as
classicalaswellBayesianregimewhenthelikelihoodfunction
∂φ (ξ)
of the observations, Uj and Vj, is Gaussian, exponential or E[ηξ(Uj)]= ξ (9)
log-normally distributed. ∂ξ
In particular, the general notation used when the likelihood σ2 =∆ Var[η (U )]= ∂2φξ(ξ) . (10)
function of the observations U =∆ [U ,...,U ]T and ηξ ξ j ∂ξ2
1 N
4
2) The moment generating function (MGF) of the statistic
η (U ) is given by Theorem 1: The likelihood maximization problems (15)
ξ j
and (16) are strictly concave and the ML estimates are given
M (h)=exp(φ (ξ+h) φ (ξ)) . (11)
ηξ ξ − ξ by
3) Tthhaet tnhoenl-onge-gpaatirvtiittiyonoffutnhcetivoanriφanξc(.e)σisη2ξcoinnv(e1x0.) implies ξˆML = (cid:80)Nj=N1ση2ξ(Uj) (17a)
4) For Gaussian, exponential and log-normally distributed ηξ
likelihoodfunctions,thelog-partitionfunctionφξ(ξ)can (cid:80)N η (V )
be expressed as a second degree polynomial given by ψˆ = j=1 ψ j . (17b)
ML Nσ2
φ (ξ)=a ξ2 . (12) ηψ
ξ ξ
Hence, the ML estimator θˆ for the clock offset can be
Thecoefficienta inthisapproximationcanbeobtained ML
ξ
written as
using the variance of the statistic η (U ), which is
ξ j ξˆ ψˆ
assumed known. Using (10), aξ is given by θˆML = ML− ML . (18)
2
σ2
a = ηξ . Proof: The ML estimate of ξ in (15) can be equivalently
ξ
2
determined by maximizing the exponent in the likelihood
If the statistical moment in (10) is not available, the
function. It can be easily verified that the exponent in (15),
empirical moment can be substituted since it readily
which is a quadratic function of ξ, is strictly concave [29].
follows from the weak law of large numbers that
A similar explanation applies to the ML estimate of ψ. The
(cid:80)N η (U ) ML estimates in (17) can be obtained by setting the first
j=1 ξ j p E[ηξ(U)], N derivative of the exponent with respect to ξ (resp. ψ) to zero.
N → →∞
(cid:80)N η2(U ) The estimate θˆML in (18) can be obtained by invoking the
j=1 ξ j p E(cid:2)η2(U)(cid:3), N . invariance principle [28]. Hence, the proof readily follows.
N → ξ →∞ 1) Gaussian Distributed Likelihood Function: A particular
Similar expressions can also be written for ηψ(Vj), Mηψ(h) application of Theorem 1 is the case when the likelihood
and φψ(ψ), respectively. Using the aforementioned properties functions f(Uj;ξ) and f(Vj;ψ) have a Gaussian distribution
of the exponential family, the ML as well as Bayesian esti- i.e., f(U ;ξ) (ξ,σ2) and f(V ;ψ) (ψ,σ2) [15].
mates of ξ and ψ are to be determined utilizing the data set Hence, itjfollo∼wsNthat ξ j ∼ N ψ
U ,V N , based on (12).
{ j j}j=1 f(U;ξ)= 1 exp(cid:32) (cid:80)Nj=1(Uj −ξ)2(cid:33) (19)
III. MAXIMUMLIKELIHOODESTIMATION (2πσξ2)N2 − 2σξ2
In this section, the ML estimates of θ are obtained by
which can be rearranged as
recastingthelikelihoodmaximizationasaninstanceofconvex
ompetnimtsizuasteiodnt.oTmhiasxaimppizroeatchhedliifkfeelrishoforodminth[e14g]r.aTphhiecaslpeacrgifiuc- f(U;ξ) exp(cid:32)ξ(cid:80)Nj=1Uj N ξ2(cid:33) .
∝ σ2 − 2σ2
cases of unconstrained and constrained likelihood functions ξ ξ
are considered separately.
By comparing with (13), we have
A. Unconstrained Likelihood η (U )= Uj, σ2 = 1 (20)
ξ j σ2 ηξ σ2
Using (5), (6) and (12), the unconstrained likelihood func- ξ ξ
tions are given by
and the ML estimate using (17a) is given by
 
f(U;ξ)∝expξ(cid:88)N ηξ(Uj)−Nσ2η2ξξ2 (13) ξˆML = (cid:80)Nj=1Uj . (21)
N
j=1
 
f(V;ψ) expψ(cid:88)N ηψ(Vj) Nση2ψψ2 . (14) Bbyy a similar reasoning, the ML estimate for ψ (17b) is given
∝ − 2
j=1 (cid:80)N V
ψˆ = j=1 j .
The ML estimates of ξ and ψ can now be expressed as ML N
 
ξˆML =arg maxexpξ(cid:88)N ηξ(Uj) Nση2ξξ2 (15) Uassing (18), the ML estimate for the offset θ can be expressed
ξ  j=1 − 2  θˆML = (cid:80)Nj=1(Uj −Vj) . (22)
ψˆML =arg maxexpψ(cid:88)N ηψ(Vj) Nση2ψψ2 . (16) The above estimate coincides ex2aNctly with the one reported in
ψ − 2
j=1 [15].
5
2) Log-NormallyDistributedLikelihoodFunction: Another where U and V denote the first order statistics of the
(1) (1)
application of Theorem 1 is when the samples Uj and Vj are samples Uj and Vj, respectively. The ML estimator θˆML for
log-normally distributed. In this case, we have the clock offset is given by
f(U;ξ)= (cid:113)21πσξ2 j(cid:89)N=1Uj−1exp(cid:32)−(cid:80)Nj=1(l2oσgξ2Uj −ξ)(cid:33) Proof: By usingθˆMaLrg=umξˆeMnLts−2sψiˆmMiLla.r to Theorem 1(a3n0d)
exp(cid:32)ξ(cid:80)Nj=1logUj N ξ2(cid:33) . (23) noting that the N constrains Uj ≥ξ (resp. Vj ≥ψ) are linear
∝ σ2 − 2σ2 functionsofξ (resp.ψ),theproofofconcavityreadilyfollows.
ξ ξ
Also,thelikelihoodmaximizationproblems(27)and(28)can
A comparison with (13) yields be equivalently expressed as
logU 1  
ηξ(Uj)= σξ2 j, ση2ξ = σξ2 . ξˆML =arg mξaxexpξ(cid:88)N ηξ(Uj)−Nσ2η2ξξ2
j=1
Following a similar line of reasoning,
such that U ξ
(1)
ηψ(Vj)= loσgψ2Vj, ση2ψ = σ1ψ2 . ψˆML =arg maxexp≥ψ(cid:88)N ηψ(Vj) Nση2ψψ2
ψ − 2
The ML estimator for θ can obtained from (18) using (17a) j=1
and (17b), and is given by such that V ψ .
(1)
≥
(cid:80)N (logU logV ) The unconstrained maximizer ξ¯ of the objective function in
θˆML = j=1 j − j . (24) (27)isgivenby(17a).Ifξ¯ U ,thenξˆ =ξ¯.Ontheother
N (1) ML
≤
hand, if U ξ¯, then the ML estimate is given by ξˆ =
(1) ML
≤
B. Constrained Likelihood U(1) using concavity of the objective function. Combining the
two cases, the ML estimate in (29a) is obtained. A similar
Using(7),(8)and(12),theconstrainedlikelihoodfunctions explanation applies to the ML estimate ψˆ in (29b) and θˆ
ML ML
are given by
follows from the invariance principle.
  1) Exponentially Distributed Likelihood Function: For the
f(U;ξ) expξ(cid:88)N ηξ(Uj) Nση2ξξ2(cid:89)N I(Uj ξ) case when the likelihood functions are exponentially dis-
∝ − 2 − tributed,thedensityfunctionofthesamplesU canbewritten
j=1 j=1 j
(25) as [14]
   
f(V;ψ)∝expψ(cid:88)N ηψ(Vj)−Nσ2η2ψψ2(cid:89)N I(Vj −ψ). f(U;ξ)=λNξ exp−λξ(cid:88)N (Uj −ξ)I(U(1)−ξ) (31)
j=1 j=1 j=1
(26)
whereλ−ξ1 isthemeanofthedelaysXj.Thedensityfunction
The resulting ML estimates of ξ and ψ can be obtained as can be rearranged as
 
ξˆML =arg maxexpξ(cid:88)N ηξ(Uj) Nση2ξξ2 f(U;ξ)∝exp(Nλξξ) .
ξ − 2 Comparing the above formulation with (25),
j=1
such that U ξ (27) η (U )=λ , σ2 =0. (32)
j ≥ ξ j ξ η
 
ψˆML =arg maxexpψ(cid:88)N ηψ(Vj) Nση2ψψ2 Using (29a), the ML estimate is given by
ψ − 2 ξˆ =U . (33)
j=1 ML (1)
such that Vj ψ . (28) Employing a similar reasoning,
≥
ψˆ =V .
ML (1)
Theorem 2: The likelihood maximization problems (27)
Using (30), the ML estimate of θ is given by
and (28) are strictly concave and the ML estimates can be
expressed as U V
θˆ = (1)− (1) , (34)
ML
(cid:32)(cid:80)N η (U ) (cid:33) 2
ξˆ =min j=1 ξ j ,U (29a) which coincides exactly with the one reported in [14], where
ML Nσ2 (1)
ηξ it is derived using graphical arguments.
(cid:32)(cid:80)N η (V ) (cid:33)
ψˆ =min j=1 ψ j ,V (29b) Remark 1: The ML estimation method outlined above dif-
ML Nσ2 (1) fers from the previous work [14] in that it is based on convex
ηψ
6
optimization,while[14]maximizedthelikelihoodgraphically. where w and v are i.i.d such that w ,v (0,σ2). The
k k k k
∼N
Hence, this approach presents an alternative view of the ML posterior pdf can be expressed as
estimation of the clock offset. It also allows us the determine
f(ξ,ψ U,V) f(ξ,ψ)f(U,V ξ,ψ)
the ML estimator of θ when the likelihood function is log- | ∝ |
N N
normally distributed. In addition, Theorem 2 will also be (cid:89) (cid:89)
=f(ξ ) f(ξ ξ )f(ψ ) f(ψ ψ )
useful in Section IV where estimation of θ in the Bayesian 0 k| k−1 0 k| k−1
regime is discussed . k=1 k=1
N
(cid:89)
f(U ξ )f(V ψ ) (37)
k k k k
· | |
k=1
IV. AFACTORGRAPHAPPROACH
where uniform priors f(ξ ) and f(ψ ) are assumed. Define
0 0
The imperfections introduced by environmental conditions δk =∆ f(ξ ξ ) (ξ ,σ2), νk =∆ f(ψ ψ )
inthequartzoscillatorinsensornodesresultsinatime-varying k−(ψ1 ,σ2k),| fk−1=∆∼f(NU kξ−)1, h =∆k−f(1V ψ ),kw|hke−re1 th∼e
k 1 k k k k k k
clock offset between nodes in a WSN. To cater for such a lNikelih−ood functions are giv|en by |
temporal variation, in this section, a Bayesian approach to
(cid:32) (cid:33)
the clock synchronization problem is adopted by representing f(U ξ ) exp ξ η (U ) ση2kξ2
the a-posteriori density as a factor graph. The inference is k| k ∝ k ξ k − 2 k
performed on the resulting factor graph by message passing (cid:32) (cid:33)
σ2
using max-product algorithm. To ensure completeness, a brief f(V ψ ) exp ψ η (V ) ηkψ2 , (38)
description of factor graphs and the max-product algorithm is k| k ∝ k ψ k − 2 k
provided below. basedon(12).Theresultingfactorgraphrepresentationofthe
A factor graph is a bipartite graph that represents a factor- posterior pdf is shown in Fig. 2.
ization of a global function as a product of local functions
called factors, each factor being dependent on a subset of Remark 2: A few important observations of this represen-
variables. Factor graphs are often used to produce a graphical tation can be summarized below.
model depicting the various inter-dependencies between a Noticethatthesubstitutionin(3)rendersacycle-freena-
collection of interacting variables. Each factor is represented • ture to the factor graph. Therefore, inference by message
byafactornodeandeachvariablehasanedgeorahalf-edge. passing on such a factor graph is indeed optimal [24].
An edge connects a particular variable to a factor node if and The two factor graphs shown in Fig. 2 have a similar
only if it is an argument of the factor expressed by the factor • structure and hence, message computations will only be
node [24]. shown for the estimate ξˆ . Clearly, similar expressions
N
Inference can be performed by passing messages will apply to ψˆ .
N
(sometimes called beliefs) along the edges of a factor In addition, only the case of constrained likelihood will be
graph. In particular, max-product algorithm is used to considered, since the case of an unconstrained likelihood is
compute the messages exchanged between variables and subsumed,aswillbeshownshortly.Theclockoffsetestimator
factor nodes. These messages can be summarized as follows θˆ can be obtained from ξˆ and ψˆ using (4).
N N N
By defining α =∆ ση2ξ,k and β =∆ η (U ), the
variable to factor node: ξ,k − 2 ξ,k ξ k
constrained likelihood function for the samples U can be
k
(cid:89) written as
m (x)= m (x) (35)
x f h x
→ h∈n(x)\f → fk ∝exp(cid:0)αξ,kξk2+βξ,kξk(cid:1)I(Uk−ξk). (39)
factor node to variable: The message passing strategy starts by sending a message
from the factor node f to the variable ξ . The variable
N N
 
ξ relays this message to the factor node δN . The factor
(cid:89) N N 1
mf x(x)=maxf(Z) mz f(z) (36) node computes the product of this message w−ith the factor
→ \{x} z∈n(f)\{x} → aδfNNte−r1‘saunmdmseanridzsintgh’eorveesrultthiengvamrieasbslaegξe t.oInthethevamriaaxb-leproξNdu−c1t
N
where Z =n(f) is the set of arguments of the local function algorithm,a‘max’functionisusedasasummarypropagation
f.Themarginaldistributionsassociatedwitheachvariablecan operator (cf. (36)). These messages are computed as
be obtained by the product of all incoming messages on the
m =f
variable. fN→ξN N
In order to sufficiently capture the temporal variations, the mξN→δNN−1 =fN
parameters ξ and ψ are assumed to evolve through a Gauss- m max δN m
Markov process given by δNN−1→ξN−1 ∝ ξN N−1· ξN→δNN−1
1 (cid:18) (ξ ξ )2(cid:19)
N N 1
=max exp − − −
ξk =ξk 1+wk ξN √2πσ2 2σ2
ψ =ψ− +v for k =1,...,N exp(cid:0)α ξ2 +β ξ (cid:1)I(U ξ )
k k−1 k · ξ,N N ξ,N N N − N
7
f1 fk fN
ξ0 ξ1 ξk ξN
δ0 δ01 ... δkk−1 δkk+1 ... δNN−1
h1 hk hN
ψ0 ψ1 ψk ψN
ν0 ν01 ... νkk−1 νkk+1 ... νNN−1
Fig.2. Factorgraphrepresentationoftheposteriordensity(37)
which can be rearranged as Upon receipt of this message, the factor node δNN−21 delivers
m max exp(cid:0)A ξ2 +B ξ2 + a product of this message and the factor δNN−21 to t−he variable
δNN−1→ξN−1 ∝ξN≤UN C ξ,Nξ Nξ +ξ,NDN−1ξ (cid:1) (40) nboedeexpξrNe−ss2edafatser summarizing over ξN−1. T−his message can
ξ,N N N 1 ξ,N N
−
where Aξ,N =∆ −2σ12 +αξ,N, Bξ,N =∆ −2σ12 m=δmNNa−−x21→ξN1−2 ∝eξxNp−(cid:18)m1≤aUx(Nξ−N1−δ1NN−−−21ξN·m−2ξ)N2−(cid:19)1→δNN−−21
∆ 1 ∆ ξN−1 √2πσ2 − 2σ2
Cξ,N = σ2, Dξ,N =βξ,N . (41) (cid:40)(cid:32) C2 (cid:33) C D (cid:41)
exp B ξ,N ξ2 ξ,N ξ,Nξ
Letξ¯N betheunconstrainedmaximizeroftheexponentinthe · ξ,N − 4Aξ,N N−1− 2Aξ,N N−1
objective function above, i.e., exp(cid:0)α ξ2 +β ξ (cid:1)I(U ξ ).
ξ¯ =argmax (cid:0)A ξ2 +B ξ2 +C ξ ξ + · ξ,N−1 N−1 ξ,N−1 N−1 N−1− N−1
N ξN ξ,N N ξ,N N−1 ξ,N N N−1 After some algebraic steps, the message above can be com-
D ξ (cid:1). pactly represented as
ξ,N N
This implies that m max exp(A ξ2 +
C ξ +D δNN−−21→ξN−2 ∝ξN−1≤UN−1 ξ,N−1 N−1
ξ¯N = ξ,N N−1 ξ,N . (42) Bξ,N 1ξN2 2+Cξ,N 1ξN 1ξN 2+Dξ,N 1ξN 1) (44)
− 2Aξ,N − − − − − − −
where
Following aline of reasoning similar toTheorem 2, itfollows
that 1 C2
ξˆN =min(cid:0)ξ¯N,UN(cid:1) . Aξ,N−1 =∆ −2σ2 +αξ,N−1+Bξ,N − 4Aξξ,N,N,
sHtaogwee.vHere,nξ¯cNe,dweepennedesdotonfξuNrt−h1e,rwtrahvicehrsiestuhnedcehtearinmbinaecdkwatarthdiss. Bξ,N−1 =∆ −2σ12, Cξ,N−1 =∆ σ12
C D
Assuming that ξ¯ U , ξ¯ from (42) can be plugged back D =∆ β ξ,N ξ,N .
in (40) which afNter≤somNe simNplification yields ξ,N−1 ξ,N−1− 2Aξ,N
(cid:40)(cid:32) C2 (cid:33) Proceeding as before, the unconstrained maximizer ξ¯N 1 of
m exp B ξ,N ξ2 the objective function above is given by −
δNN−1→ξN−1 ∝ ξ,N − 4Aξ,N N−1−
Cξ2,ANξD,Nξ,NξN−1(cid:41). (43) ξ¯N−1 =−Cξ,N−12ξNA−ξ,2N+−1Dξ,N−1
andthesolutiontothemaximizationproblem(44)isexpressed
Themessagepassedfromthevariableξ tothefactornode as
N 1
δNN−21 is the product of the message (4−3) and the message ξˆN 1 =min(cid:0)ξ¯N 1,UN 1(cid:1) .
rec−eived from the factor node f , i.e., − − −
N 1
− Again, ξ¯ depends on ξ and therefore, the solution
N 1 N 2
mξN−1→δNN−−21 =mδNN−1→ξN−1 ·mfN−1→ξN−1 . demands a−nother traversal ba−ckwards on the factor graph
8
representation in Fig. 2. By plugging ξ¯ back in (44), it Proof: The constants A , C and D are defined in
N 1 ξ,k ξ,k ξ,k
follows that − (41) and (46). The proof follows by noting that −Cξ,k > 0
2Aξ,k
which implies that g (.) is a monotonically increasing func-
m ξ,k
δNN−−(cid:40)21(cid:32)→ξN−2 ∝ (cid:33) (cid:41) tion.
exp Bξ,N−1− 4CAξ2ξ,N,N−11 ξN2−2− Cξ,N2A−ξ1,DNξ,1N−1ξN−2 Wbeitchonthveennioetnattliyonwgrξit,tke(n.)a,sthefollowingchainofequalitiescan
− − (45) ξ¯ =g (cid:16)ξˆ(cid:17)
1 ξ,1 0
which has a form similar to (43). It is clear that one can keep (cid:16) (cid:16) (cid:17)(cid:17)
ξˆ =min U ,g ξˆ
traversing back in the graph yielding messages similar to (43) 1 1 ξ,1 0
and (45). In general, for i=1,...,N −1 we can write ξ¯2 =gξ,2(cid:16)ξˆ1(cid:17)
1 C2 (cid:16) (cid:16) (cid:17)(cid:17)
Aξ,N−i =∆ −2σ2 +αξ,N−i+Bξ,N−i+1− 4Aξξ,N,N−ii++11 ξˆ2 =min U2,gξ,2 ξˆ1
−
B =∆ 1 , C =∆ 1 (46) where
ξ,N−i −2σ2 ξ,N−i σ2 (cid:16) (cid:17) (cid:16) (cid:16) (cid:16) (cid:17)(cid:17)(cid:17)
∆ Cξ,N i+1Dξ,N i+1 gξ,2 ξˆ1 =gξ,2 min U1,gξ,1 ξˆ0
Dξ,N−i =βξ,N−i− 2−Aξ,N−i+1− =min(cid:16)gξ,2(U1),gξ,2(cid:16)gξ,1(cid:16)ξˆ0(cid:17)(cid:17)(cid:17) (54)
and
ξ¯N i = Cξ,N−iξN−i−1+Dξ,N−i (47) wexhperreess(e5d4)asfollows from Lemma 1. The estimate ξˆ2 can be
− − 2Aξ,N i
ξˆN i = min(cid:0)ξ¯N i,UN i−(cid:1) . (48) ξˆ =min(cid:16)U ,min(cid:16)g (U ),g (cid:16)g (cid:16)ξˆ(cid:17)(cid:17)(cid:17)(cid:17)
− − − 2 2 ξ,2 1 ξ,2 ξ,1 0
Using (47) and (48) with i=N −1, it follows that =min(cid:16)U2,gξ,2(U1),gξ,2(cid:16)gξ,1(cid:16)ξˆ0(cid:17)(cid:17)(cid:17) .
C ξ +D
ξ¯1 = ξ,1 0 ξ,1 (49) By following the same procedure, one can write
− 2A
ξ,1
ξˆ1 =min(cid:0)ξ¯1,U1(cid:1) . (50) ξˆ3 =min(cid:16)U3,gξ,3(U2),gξ,3(gξ,2(U1)),
(cid:16) (cid:16) (cid:16) (cid:17)(cid:17)(cid:17)(cid:17)
Similarly, by observing the form of (43) and (45), it follows g g g ξˆ .
ξ,3 ξ,2 ξ,1 0
that
m exp(cid:40)(cid:32)B Cξ2,1 (cid:33)ξ2 Cξ,1Dξ,1ξ (cid:41) . For m≥j, define
δ01→ξ0 ∝ ξ,1− 4Aξ,1 0 − 2Aξ,1 0 Gmξ,j(.)=∆ gξ,m(gξ,m 1...gξ,j(.)) . (55)
(51) −
The estimate ξˆ can be obtained by maximizing the received The estimate ξˆ can, therefore, be compactly represented as
0 3
message in (51). It can be noticed from the structure of the (cid:16) (cid:16) (cid:17)(cid:17)
factorgraphthatthismaximizationisinherentlyunconstrained ξˆ3 =min U3,G3ξ,3(U2),G3ξ,2(U1),G3ξ,1 ξˆ0 .
i.e.,
Hence, one can keep estimating ξˆ at each stage using this
k
ξˆ0 =ξ¯0 =mξa0x mδ01→ξ0 strategy. Note that the estimator only depends on functions of
data and can be readily evaluated.
C D
⇒ξˆ0 = 4A Bξ,1 ξ,1C2 . (52) Inordertoderiveanalogousexpressionsforψ,asimilarline
ξ,1 ξ,1− ξ,1 of reasoning should be followed. In particular, the constants
The estimate in (52) can now be substituted in (49) to yield A , B , C and D for i = 0,...,N 1,
ξ,N i ξ,N i ξ,N i ξ,N i
ξ¯, which can then be used to solve for ξˆ in (50). Clearly, can b−e obtain−ed straig−htforwardly−from (41) and (46)−by
1 1
this chain of calculations can be continued using recursions substituting α and β with α and β ,
ξ,N i ξ,N i ψ,N i ψ,N i
(47) and (48). respectively. Usin−g these cons−tants, ψˆ , g − and Gm c−an
0 ψ,k ψ,j
Define be defined analogously to (52), (53) and (55).
C x+D
g (x)=∆ ξ,k ξ,k . (53) Generalizingthisframework,theclosedformexpressionfor
ξ,k
− 2Aξ,k theclockoffsetestimateθˆN isgivenbythefollowingtheorem.
A key property of the function g (.), which proves useful
in the quest for a closed form soluξt,ikon, can be summarized in Theorem 3: ThestateestimatesξˆN andψˆN fortheposterior
the following lemma. pdf in (37) can be expressed as
Lemma 1: For real numbers a and b, the function g (.) (cid:16) (cid:16) (cid:17)(cid:17)
ξ,k ξˆ =min U ,GN (U ),...,GN (U ),GN ξˆ
defined in (53) satisfies N N ξ,N N 1 ξ,2 1 ξ,1 0
−
(cid:16) (cid:16) (cid:17)(cid:17)
ψˆ =min V ,GN (V ),...,GN (V ),GN ψˆ
gξ,k(min(a,b))=min(gξ,k(a),gξ,k(b)) . N N ψ,N N−1 ξ,2 1 ξ,1 0
9
and the factor graph based clock offset estimate (FGE) θˆ is Using these values for α and β , (41) and (46) can be
N ξ,k ξ,k
given by written as
1 1 1
θˆN = ξˆN −ψˆN . (56) Aξ,N =−2σ2 − 2σξ2,N, Bξ,N =−2σ2
2
1 U
C = , D = N (59)
ξ,N σ2 ξ,N σ2
ξ,N
Proof: The proof follows from the discussion above and
1 1 C2
using (4). Aξ,N−i =−2σ2 − 2σξ2,N i +Bξ,N−i+1− 4Aξξ,N,N−ii++11
− −
Remark 3: The closed form expressions in Theorem 3 1 1
B = , C =
enabletheestimationoftheclockoffset,whenitmaybetime- ξ,N−i −2σ2 ξ,N−i σ2
U C D
hvaarvyeinagGaanudssthiaenl,ikeexlpihoonoendtifaulnoctriolongs,-nfo(rUmka|lξkd)isatrnidbufti(oVnk.|ψk), Dξ,N−i = σξ2N,N−ii − ξ,N2−Ai+ξ,1N ξi,+N1−i+1
− −
fori=1,...,N 1.Usingsimilararguments,itcanbeshown
that the estimate−ψˆ is given by
N
(cid:16) (cid:17)
A. Gaussian Distributed Likelihood Function ψˆN =GNψ,1 ψˆ0 .
It follows from (4) that the FGE, θˆ , can be expressed as
A particular case of the Bayesian framework described N
fab(o(Vξvke|,ψσko2)cc)huaravsnedwfah(eVGnauψlsisk)iealnihodois(dtψribf,uuσtni2cotnio,),nisi..ee..,f, (fU(kU|kξk|ξ)k)an∼d θˆN = GNξ,1(cid:16)ξˆ0(cid:17)−2GNψ,1(cid:16)ψˆ0(cid:17) . (60)
N k ξ,k k| k ∼N k ψ,k
The behavior of θˆ can be further investigated for the case
N
(cid:40) (cid:41) when the noise variance σ2 in the Gauss-Markov model goes
1 (U ξ )2
f(Uk|ξk)= (cid:113)2πσξ2,k exp − k2σ−ξ2,kk to zero. Considerg (ξ)= Cξ,Nξ+Dξ,N
(cid:32) (cid:33) ξ,N − 2A
ξ U ξ2 ξ,N
∝exp 2kσ2k − 2σk2 . (57) wheretheconstantsAξ,N,Bξ,N,Cξ,N andDξ,N aregivenby
ξ,k ξ,k (59). After some algebraic steps, we have
σ2ξ+σ2U
The aforementioned Gaussian distribution constitutes an un- g (ξ)= ξ N .
constrained likelihood function, i.e., the domain of the pdf is ξ,N σ2+σ2
ξ
independentoftheunknownparameterξ .Itisclearfromthe
k As σ2 0, g (ξ) ξ. Similarly, it can be shown that
ξ,N
message passing approach that at each stage k of the factor → →
g (ξ) ξ as σ2 0. Hence, it follows that in the low
graph, the unconstrained maximizer ξ¯k is the actual solution syξ,sNte−m1nois→e regime, a→s σ2 0
to the likelihood maximization problem →
C D
ξˆ ξˆ = ξ,1 ξ,1 .
N → 0 4A B C2
max exp(cid:0)A ξ2+B ξ2 +C ξ ξ +D ξ (cid:1) ξ,1 ξ,1− ξ,1
ξ,k k ξ,k k 1 ξ,k k k 1 ξ,k k
ξk − − Similarly, it can be shown that
C D
i.e., ξˆk = ξ¯k ∀k = 1,...,N. Hence, the unconstrained ψˆN →ψˆ0 = 4Aψ,1Bψ,ψ1,1ψ−,1Cψ2,1 .
likelihood maximization problem is subsumed in the message
Therefore
passing framework for constrained likelihood maximization. ξˆ ψˆ
It follows from Theorem 3 that ξˆ for Gaussian distributed θˆ 0− 0 ,
N N → 2
observations U in (57) is given by
k
which can be proven to be equal to the ML estimator (22).
(cid:16) (cid:17)
ξˆN =GNξ,1 ξˆ0 B. Log-Normally Distributed Likelihood Function
The log-normally distributed likelihood function in the
whereξˆ andGN (.)aredefinedin(52)and(55),respectively. Bayesian regime can be expressed as
0 ξ,1 (cid:32) (cid:33)
Evaluating ξˆN requires to determine the constants in (41) and 1 (logUk ξk)2
f(U ξ )= exp −
(46). By comparing (57) with (39), we have k| k U σ √2π − 2σ2
k ξ,k ξ,k
(cid:32) (cid:33)
ξ log(U ) ξ2
αξ,k =−2σ12 , βξ,k = σU2k . (58) ∝exp k2σξ2,kk − 2σkξ2,k . (61)
ξ,k ξ,k
10
By comparing (61) and (39), we have and so on. The estimator ξˆ at the last step can be evaluated
0
as
1 logU
k C D
αξ,k =−2σξ2,k βξ,k = σξ2,k . ξˆ0 = 4Aξ,1Bξ,1ξ,1ξ−,1Cξ2,1 =+∞.
Clearly,theonlydifferenceherewiththeGaussiandistribution
This implies that
is a redefinition of βξ,k. The expression of ξˆN in this case is GN (ξˆ)=+ . (64)
again ξ,1 0 ∞
(cid:16) (cid:17)
ξˆ =GN ξˆ Using (64) and Theorem 3, it readily follows that
N ξ,1 0
ξˆ =min(U ,U +λ σ2,U +2λ σ2,
where GNξ,1(.) and ξˆ0 are given by (55) and (52), respectively. N N N−1...,ξU +(NN−2 1)λξσ2). (65)
The recursively evaluated constants in (41) and (46) can be 1 ξ
−
written as Similarly, for the pdf [14]
1 1 1
Aξ,N =−2σ2 − 2σ2 , Bξ,N =−2σ2 f(Vk|ψk)=λψexp(−λψ(Vk−ψk))I(Vk−ψk)
ξ,N exp(λ ψ )I(V ψ ),
1 logUN ∝ ψ k k− k
C = , D =
ξ,N σ2 ξ,N σ2 the estimate ψˆ is given by
ξ,N N
Aξ,N−i =−2σ12 − 2σξ21,N i +Bξ,N−i+1− 4CAξ2ξ,N,N−ii++11 ψˆN =min(VN,VN−1.+..λ,ψVσ2+,V(NN−2+1)2λλψσσ22), (66)
1 − 1 − 1 − ψ
Bξ,N−i =−2σ2, Cξ,N−i = σ2 and the estimate θˆN can be obtained using (56), (65) and (66)
logU C D as
N i ξ,N i+1 ξ,N i+1
Dξ,N−i = σξ2,N −i − 2−Aξ,N i+1− θˆ = 1min(U ,U +λ σ2,U +2λ σ2,
− − N 2 N N−1 ξ N−2 ξ
for i = 1,...,N 1. Similar arguments apply to ψ. Hence,
the FGE θˆN can b−e expressed as ...,U1+(N −1)λξσ2)−
1
(cid:16) (cid:17) (cid:16) (cid:17) min(V ,V +λ σ2,V +2λ σ2,
θˆN = GNξ,1 ξˆ0 −2GNψ,1 ψˆ0 . (62) 2 N N−1...,ψV1+(NN−−2 1)λψψσ2). (67)
Again, as the Gauss-Markov system noise σ2 0, the above As the Gauss-Markov system noise σ2 0, (67) yields
→ →
estimator approaches its ML counterpart (24).
min(U ,...,U ) min(V ,...,V )
θˆ θˆ = N 1 − N 1 ,
N ML
→ 2
C. Exponential Distribution
which is the ML estimator given by (34).
Theorem 3 can also be used to derive a Bayesian estimator
ξˆN fortheexponentiallydistributedlikelihoodcaseconsidered V. CLASSICALANDBAYESIANBOUNDS
in [14]. In this case, we have
Toevaluatetheperformanceoftheestimatorsderivedinthe
f(U ξ )=λ exp( λ (U ξ ))I(U ξ ) precedingsections,classicalaswellasBayesianlowerbounds
k k ξ ξ k k k k
| − − −
exp(λ ξ )I(U ξ ) (63) onthevarianceoftheestimatorsarediscussed.Theplacement
ξ k k k
∝ − ofalowerboundallowsonetocompareestimatorsbyplotting
where λ−ξ1 is the mean network delay of Xk. A comparison their performance against the bound. It must be emphasized
of (63) with (39) reveals that herethat theresultsinthis sectionassumeno specificformof
the log-partition function and are therefore, valid for arbitrary
α =0, β =λ .
ξ,k ξ,k ξ distributionsfromtheexponentialfamily,whichisawideclass
For these values of α and β , the constants A , B , and contains almost all distributions of interest. Hence, these
ξ,k ξ,k ξ,k ξ,k
C and D are given by results are fairly general and can be useful in their own right
ξ,k ξ,k
in classical as well as Bayesian parameter estimation theory,
1 1
Aξ,k =−2σ2, Bξ,k =−2σ2 and at the same time will be used as a stepping stone towards
1 comparing the estimators developed thus far.
Cξ,k = σ2, Dξ,k =λξ The likelihood function of the data is considered an
arbitrary member of the exponential family of distributions.
for all k =1,...,N. Using Theorem 3, we have
In addition, depending on whether the domain of the
GNξ,N(UN−1)=−Cξ,NU2NA−ξ1,N+Dξ,N lciakseelsihooofdudnecpoennsdtrsaionnedthaespwarealml eatesr ctoonsbteraiensetidmalitkedel,ihboootdh
=U +λ σ2 . functions are discussed to maintain full generality. The
N 1 ξ
− general expressions for the unconstrained and constrained
Similarly it can be shown that
likelihood functions for observations Z =∆ [Z ,...,Z ]T are
1 N
GN (U )=U +2λ σ2 given by
ξ,N 1 N 2 N 2 ξ
− − −