# A Comparison of Parametric and Sample-Based Message Representation in Cooperative Localization

Authors: Jaime Lien; Ulric J. Ferner; Warakorn Srichavengsup; Henk Wymeersch; Moe Z. Win
Journal: International Journal of Navigation and Observation
 (2012)

---

## Abstract

Location awareness is a key enabling feature and fundamental challenge in present and future wireless networks.
Most existing localization methods rely on existing infrastructure and thus lack the flexibility and robustness necessary for large ad hoc networks. In this paper, we build upon SPAWN (sum-product algorithm over a wireless network), which determines node locations through iterative message passing, but does so at a high computational cost. We compare different message representations for SPAWN in terms of performance and complexity and investigate several types of cooperation based on censoring. Our results, based on experimental data with ultra-wideband (UWB) nodes, indicate that parametric message representation combined with simple censoring can give excellent performance at relatively low complexity.


---

## Body

## 1. Introduction

Location awareness has the potential to revolutionize a diverse array of present and future technologies. Accurate knowledge of a user's location is essential for a wide variety of commercial, military, and social applications, including next-generation cellular services [1, 2], sensor networks [3, 4], search-and-rescue [5, 6], military target tracking [7, 8], health care monitoring [9, 10], robotics [11, 12], data routing [13, 14], and logistics [15, 16]. Typically, only a small fraction of the nodes in the network, known as anchors, have prior knowledge about their location. The remaining nodes, known as agents, must determine their locations through a process of localization or positioning. The ad hoc and often dynamic nature of wireless networks requires distributed and autonomous localization methods. Moreover, location-aware wireless networks are frequently deployed in unknown environments and hence can rely only on minimal (if any) infrastructure, human maintenance, and a priori location information.Cooperation is an emerging paradigm for localization in which agents take advantage of network connections and interagent measurements to improve their location estimates. Non-Bayesian cooperative localization in wireless sensor networks is discussed in [17]. Different variations of Bayesian cooperation have been considered, including Monte-Carlo sequential estimation [18] and nonparametric belief propagation in static networks [19]. For a comprehensive overview of Bayesian and non-Bayesian cooperative localization in wireless networks, we refer the reader to [20], which also introduces a distributed cooperative algorithm for large-scale mobile networks called SPAWN (sum-product algorithm over a wireless network). This message-passing algorithm achieves improved localization accuracy and coverage compared to other methods and will serve as the basic algorithm in this paper.The complexity and cost associated with the SPAWN algorithm depend largely on how messages are represented for computation and transmission. As wireless networks typically operate under tight power and resource constraints, the choice of message representation heavily impacts the feasibility and ease of implementation of the algorithm. The method of message representation and ensuing tradeoff between communication cost and localization performance are thus of great practical importance in the deployment of realistic localization systems. Particle methods do not necessarily lend themselves well in practice to be exchanged wirelessly between devices, due to their high computational complexity and communication overhead [21]. Other message-passing methods have been developed that rely on parametric message representation, thus alleviating these drawbacks but limiting representational flexibility. In particular, in [22], the expectation propagation algorithm is considered with Gaussian messages, while in [23], variational message passing with parametric messages is shown to exhibit low complexity. A variation of SPAWN combining GPS and UWB was evaluated in [24], using a collection of parametric distributions with ellipsoidal, conic, and cylindrical shapes.This paper addresses the need for accurate, resource-efficient localization with an in-depth comparison of various message representations for SPAWN. We describe and evaluate different parametric and nonparametric message representations in terms of complexity and accuracy. Additionally, we analyze the performance of various cooperative schemes and message representations in a simulated large-scale ultra-wide bandwidth (UWB) network using experimental UWB ranging data. UWB is an attractive choice for ranging and communication due to its ability to resolve multipath [25, 26], penetrate obstacles [27], and provide high resolution distance measurements [28, 29]. Recent research advances in UWB signal acquisition [30, 31], multiuser interference [32, 33], multipath channels [34, 35], non-line-of-sight (NLOS) propagation [28, 29], and time-of-arrival estimation [36] increase the potential for highly accurate UWB-based localization systems in harsh environments. Consequently, significant attention has been paid to both algorithm design [37–42] and fundamental limits of accuracy [43–46] for UWB localization. It is expected that UWB will be exploited in future location-aware systems that utilize coexisting networks of sensors, controllers, and peripheral devices [47, 48].


## 2. Problem Formulation

We consider a wireless network ofN nodes in an environment ℰ. Time is slotted with nodes moving independently from time slot to time slot. The position of node i at time t is described by the random variable xi(t); the vector of all positions is denoted by x(t). At each time t, node i may collect internal position-related measurements zi,self(t), for example, from an inertial measurement unit. The set of all internal measurements is denoted by zself(t). Within the network, nodes communicate with each other via wireless transmissions. We denote the set of nodes from which node i can receive transmissions at time t by 𝒮→i(t). Note that the communication link may not be bidirectional; that is, j∈𝒮→i(t) does not imply i∈𝒮→j(t). Using packets received from j∈𝒮→i(t), node i may collect a set of relative measurements, represented by the vector zj→i(t), which we will limit to distance measurements. We denote the set of all relative measurements made in the network at time t by zrel(t). The full set of relative and internal measurements is denoted z(t).The objective of the localization problem is for each nodei to determine the a posteriori distribution p(xi(t)∣z(1:t)) of its position xi(t) at each time t, given information up to and including t.


## 3. A Brief Introduction to SPAWN

In [20], we proposed a cooperative localization algorithm by factorizing the joint distribution p(x(0:T)∣z(1:T)), formulating the problem as a factor graph with temporal and spatial constraints, and applying the sum-product algorithm. This leads to a distributed algorithm, known as SPAWN, presented in Algorithm 1. The aim of SPAWN is to compute a beliefbi(t)(xi(t)) available to node i at the end of any time slot t, which serves as an approximation of the marginal a posteriori distribution p(xi(t)∣z(1:t)). Note that each operation of SPAWN requires only information local to an individual node. Information is shared between nodes via physical transmissions. Each node can therefore perform the computations in Algorithm 1 using its local information and transmissions received from neighboring nodes.Algorithm 1:SPAWN.
(
1)  Initialize belief b(0)(xi(0))=p(xi(0)), ∀i
(
1
)
for
t
=
1 to Tdo{time index}
(
2
)  forallido{mobility update}
(
5
)   Mobility update:
b~i(t)(xi(t))∝
∫p(xi(t)|xi(t-1))p(zi,self(t)|xi(t-1),xi(t))
×bi(t-1)(xi(t-1))dxi(t-1)      (A-1)
(
6
)  end for
(
7
)  Initialize bi(t)(xi(t))=b~i(t)(xi(t)), ∀i
(
8
)  forl=1 to Nitdo{iteration index; begin cooperative
update}
(
9
)   for allido
(
10
)   for allj∈S→i(t)do
(
11
)   Receive and convert bj(t)(xj(t)) to a distribution
cj→i(l)(xi(t)):
cj→i(l)(xi(t))∝
∫p(zj→i(t)|xi(t),xj(t))bj(t)(xj(t))dxj(t)    (A-2)
(
12
)  Update and broadcast bi(t)(xi(t)):
bi(t)(xi(t))∝b~i(t)(xi(t))∏k∈S→i(t)ck→i(l)(xi(t))  (A-3)
(
13
)   end for
(
14
)   end for
(
15
)  end for  {end cooperative update}
(
16
)
end for  {end current time step}Observe that Algorithm1 contains a number of key steps.(i)
Mobility update (line 4), requiring knowledge of mobility models p(xi(t)∣xi(t-1)) and self-measurement likelihood functions p(zi,self(t)∣xi(t-1),xi(t)).(ii)
Message conversion (line 10) of position information from neighboring devices to account for relative measurements, requiring knowledge of the neighbors and of relative measurement likelihood functions p(zj→i(t)∣xi,xj).(iii)
Belief update (line 11), to fuse information from the mobility update with information from the current neighbors.The first two operations can be interpreted as message filtering, while the latter operation is a message multiplication. How these operations can be implemented in practice will be the topic of Section 4.


## 4. Message Representation

### 4.1. Key Operations

In SPAWN, probabilistic information is exchanged and computed throughmessages. The manner in which these messages are represented for transmission between nodes and internal computation is closely related to the complexity and performance of the localization algorithm. In traditional communications problems, such as decoding, messages can be represented efficiently and exactly through, for instance, log-likelihood ratios [49]. In SPAWN, exact representation is impossible, so we must resort to different types of approximate message representations. Any representation must be able to capture the salient properties of the true message and must enable efficient computation of the key steps in SPAWN, namely, message filtering(A-1)-(A-2) and message multiplication(A-3). We consider three types of message representation: discretized, sample-based, and parametric.For convenience, we will introduce a set of new notations. For the filtering operation, the incoming message is denoted bypX(x), the filtering operation by h(x,y), and the outgoing message by pY(y), with
(1)pY(y)∝∫h(x,y)pX(x)dx.
For the multiplication operation, we assume M incoming messages pX(i)(x) (i=1,…,M) over a single variable X, and an outgoing message
(2)ϕX(x)∝∏i=1M‍pX(i)(x).
Note that (1) maps to (A-1) through the following association:
(3)x→xi(t-1),y→xi(t),h(x,y)→p(xi(t)∣xi(t-1))p(zi,self(t)∣xi(t-1),xi(t)),pX(x)→bi(t-1)(xi(t-1)).
Similarly, (1) maps to (A-2) through the following association:
(4)x→xj(t),y→xi(t),h(x,y)→p(zj→i(t)∣xi(t),xj(t)),pX(x)→bi(t)(xi(t)).


### 4.2. Discretized Message Representation

A naive but simple approach to represent a continuous distributionpX(x) is to uniformly discretize the domain of X, yielding a set of quantization points Q={x1,…,xR}. The distribution is then approximated as a finite list of values, {pX(xk)}k=1R. The filtering operation then becomes
(5)pY(yk)∝∑l=1Rh(xl,yk)pX(xl),
requiring 𝒪(R2) operations. The multiplication becomes
(6)ϕX(xk)∝∏i=1M‍pX(i)(xk),
requiring 𝒪(RM) operations. Because R scales exponentially with the dimensionality of X and a large number of points are required in every dimension to capture fine features of the messages, discretization is impractical for SPAWN in UWB localization.


### 4.3. Sample-Based Message Representation

A sample-based message representation, as used in [19, 50], overcomes the drawback of discretization by representing messages as samples, concentrated where the messages have significant mass. Before describing the detailed implementation of the filtering and multiplication operations, we give a brief overview of generic sampling techniques (see also [51, 52]) and kernel density estimation (KDE).


#### 4.3.1. Background: Sampling and Kernel Density Estimation

We say that a list of samples with associated weights{xk,wk}k=1R is a representation for a distribution pX(x) if, for any integrable function g(x), we have the following approximation:
(7)I=∫g(x)pX(x)dx≈∑k=1Rwkg(xk).
Popular methods for obtaining the list of weighted samples include (i) direct sampling, where we draw R i.i.d. samples from pX(x), each with weight 1/R; and (ii) importance sampling, where we draw R i.i.d. samples from a distribution qX(x), with a support that includes the support of pX(x), and set the weight corresponding to sample xk as wk=pX(xk)/qX(xk). In both cases, it can easily be verified that the approximation is unbiased with mean I and variance that reduces with R (and that depends on qX(x), for importance sampling). Most importantly, the variance does not depend on the dimensionality of x.A variation of importance sampling that is not unbiased but that often has smaller variance is obtained by setting the weights as follows:wk∝pX(xk)/qX(xk), ∑k‍wk=1. This approach has the additional benefit that it does not require knowledge of the normalization constants of pX(x) or qX(x). A list of R equally weighted samples can be obtained from {xk,wk}k=1R through resampling, that is, by drawing (with repetition) R samples from the probability mass function defined by {xk,wk}k=1R.For numerical stability reasons, weights are often computed and stored in the logarithmic domain, that is,λk=logpX(xk)-logqX(xk). When the distributions involved contain exponentials or products, the log-domain representation is also computationally efficient. Operations such as additions can be evaluated efficiently in the log-domain as well, using the Jacobian logarithm [49, pages 90–94]. Once all R log-domain weights are computed, they are translated, exponentiated, and normalized: wk∝exp(λk-maxlλl).Given a sample representation{xk,wk}k=1R of a distribution pX(x), we obtain a kernel density estimate of pX(x) as
(8)p^X(x)=∑k=1RwkKσ(x-xk),
where Kσ(x) is the so-called kernel with bandwidth σ. The kernel is a symmetric distribution with a width parameter that is tuned through σ. For instance, a two-dimensional Gaussian kernel is given by
(9)Kσ(x)=12πσ2exp(-∥x∥22σ2).
While the choice of kernel affects the performance of the estimate to some limited extent (e.g., in an MMSE sense, where the error is ∫|pX(x)-p^X(x)|2pX(x)dx), the crucial parameter is the bandwidth σ, which needs to be estimated from the samples {xk,wk}k=1R. A large choice of σ makes p^X(x) smooth, but it may no longer capture the interesting features of pX(x). When σ is too small, p^X(x) may exhibit artificial structure not present in pX(x) [53].With this background in sampling techniques and KDE, we return to the problem at hand: filtering and multiplication of messages.


#### 4.3.2. Message Filtering

We assume a message representation ofpX(x) as {xk,wk}k=1R and wish to obtain a message representation of pY(y)∝∫‍h(x,y)pX(x)dx. Let us interpret h(x,y) as a conditional distribution pY∣X(y∣x), up to some arbitrary constant. Suppose we can draw samples {[xk,yk],wk}k=1R~pY∣X(y∣x)pX(x); then {yk,wk}k=1R will form a sample representation of pY(y). Now the problem reverts to drawing samples from pY∣X(y∣x)pX(x). This can be accomplished as follows: first, for every sample xk, draw yk~qY∣X(y∣xk) from some distribution qY∣X(y∣xk). Second, set the weight of sample [xk,yk] as
(10)vk=wkpY∣X(yk∣xk)qY∣X(yk∣xk).
Finally, renormalize the weights vk to vk/∑lvl. The complexity of the filtering operation scales as 𝒪(R), a significant improvement from 𝒪(R2) for discretization. In addition, R can generally be much smaller in a particle-based representation.Let us consider some examples of the filtering operation in SPAWN.(i)
Mobility update(A-1): let pX(x) be the belief before movement (represented by {xk,wk}k=1R) and pY(y) the belief after movement. Assume that we are able to measure perfectly the distance traveled (given by zself), but have no information regarding the direction, and furthermore that the direction is chosen uniformly in (0,2π]. In that case,
(11)h(x,y)∝δ(zself-∥x-y∥),
where δ is a Dirac delta function, so that qY∣X(y∣xk)=pY∣X(y∣xk)∝δ(zself-∥xk-y∥) is a reasonable choice. For every xk, we can now draw values for y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and setting r=zself, leading to yk=xk+zself×[cosθksinθk]T, with vk∝wk.(ii)
Ranging update(A-2): let pX(x) be a message (represented by {xk,wk}k=1R) from a node with which we have performed ranging, resulting in a range estimate z. Let h(x,y)=pZ∣D(z∣d), where d=∥x-y∥. Note that pZ∣D(z∣d) is a likelihood function, since the measurement z is known. Assume that we have a model for the ranging performance in the form of distributions pZ∣D(z∣d) for any value of d. We then sample pY(y) as follows: for every xk, draw y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and r~qR∣Z(r∣z), for some well-chosen qR∣Z(r∣z)(e.g., a Gaussian distribution with mean equal to the distance estimate, z, and a standard deviation that is sufficiently large with respect to the standard deviation of p(z|d) for any d). The weights are set as
(12)vk=wkpZ∣D(z∣rk)qR∣Z(rkz).


#### 4.3.3. Message Multiplication

Here we assume message representations{xk(i),wk(i)}k=1R for pX(i)(xk), i=1,…,M. In contrast to the discretization approach, we cannot directly compute ∏i=1M‍pX(i)(x) for arbitrary values of x. Rather, for every message pX(i)(xk), we create a KDE p^X(i)(x)=∑k=1R‍wkKσ(i)(x-xk(i)) with a Gaussian kernel and a bandwidth estimated using the methods from [53]. Suppose we now draw R samples from a distribution qX(x); then the weights are
(13)vk∝∏i=1M‍p^X(i)(xk)qX(xk),
which can be computed efficiently in the log-domain. A reasonable choice for qX(x) could be one of the incoming messages pX(i)(x) (e.g., the one with the smallest entropy) or a mixture of the incoming messages. The computational complexity of the message multiplication operation scales as 𝒪(MR2). This appears worse than the discretized case (complexity 𝒪(MR)), but note that R is much smaller for sample-based representations than for discretization (e.g., R=103 or R=104 for the sample-based representation compared to R=108 in the discretization).


### 4.4. Parametric Message Representation

#### 4.4.1. Choosing a Suitable Parameterization

From the previous section, it is clear that the bottleneck of the sample-based message representation lies in the message multiplication, which scales quadratically with the number of samples. An alternative approach is to represent each message as a set of parameters (e.g., a Gaussian distribution characterized by a mean and covariance matrix). In contrast to the sample-based message representation, which can represent messages of any shape, parametric representations must be specially tailored to the problem at hand. For example, single two-dimensional Gaussian parametric messages are utilized in [22] for localization with both range and angle measurements. Our choice of parametric message is based on the following observations.(i)
For the filtering operation with a two-dimensional Gaussian inputpX(x), the output pY(y) can be approximated by a circular distribution with the same mean for both the mobility update (A-1) and the ranging update (A-2).(ii)
Multiplying Gaussian distributions yield a Gaussian distribution.(iii)
The multiplication of multiple circular distributions can be approximated by a Gaussian distribution or a mixture of Gaussian distributions.We will use as a basic building block the following distribution in two dimensions:(14)𝒟(x;m1,m2,σ2,ρ)=1C(σ2,ρ)exp{-[(x1-m1)2+(x2-m2)2-ρ]22σ2},
where [m1,m2] is the midpoint of the distribution, ρ is the radius, σ2 is the variance, and C(σ2,ρ) is a normalization constant equal to
(15)C(σ2,ρ)=2πσ2[exp(-ρ22σ2)+122πρ2σ2(1+erfρ22σ2)].
As a special case, we note that, when ρ=0, (14) reverts to a two-dimensional Gaussian. Moreover, we will represent all messages as a mixture of two distributions of the type (14), so that
(16)pX(x)=12𝒟(x;m1(a),m2(a),σ2,ρ)+12𝒟(x;m1(b),m2(b),σ2,ρ),
which can be represented by the six-dimensional vector [m1(a),m2(a),m1(b),m2(b),ρ,σ2]. We will denote the family of distributions of the form (16) by 𝒟2. Note that it is trivial to extend this distribution, which is designed for two-dimensional localization systems, for use in three-dimensional systems. Before we describe the message filtering and message multiplication operations, let us first show how the parameters of (14) can be estimated from a list of samples.


#### 4.4.2. ML Estimation of the Parametersm1,m2,σ2,ρ.

Given a list of samples{xk}k=1R, we can estimate the parameters [m1,m2,σ2,ρ] as follows. The midpoint m=[m1,m2] is estimated by
(17)m^=1R∑k=1Rxk.
To find the radius ρ and variance σ2 of 𝒟-distribution, we use maximum likelihood (ML) estimation, assuming the R samples are independent. Introducing α=[ρσ2]T, we find that
(18)α^ML=argmaxα∑k=1RlogpX(xk;m^,α)︸Λ(α).
Treating the log-likelihood function (LLF) Λ(α) as an objective function, we find its maximum through the gradient ascent algorithm
(19)α^(n+1)=α^(n)-ε∇αΛ(α)|α=α^(n),
where ε is a suitably small step size, and the gradient vector can be approximated using finite differences. To initialize (19), we consider two initial estimates for σ2 and ρ: one assuming ρ≪σ and a second assuming ρ≫σ. The LLF is evaluated for both preliminary solutions, and the one with the largest log-likelihood is used as the initial estimate in (19).


#### 4.4.3. Message Filtering

To perform message filtering, we use the fact that sample-based message filtering is a low-complexity operation. We decomposepX(x), represented in parametric form, into its two mixture components. From each component, we draw R/2 samples and perform sample-based message filtering, as outlined in Section 4.3.2. We can then estimate the new 𝒟-parameters for each mixture component using the ML method described above. We thus have pY(y) in parametric form. The complexity of this operation scales as 𝒪(R).


#### 4.4.4. Message Multiplication

The motivation for using the parametric message representation is to avoid the complexity associated with sample-based message multiplication. GivenM distributions px(i)(x)∈𝒟2, our goal is to compute
(20)ϕX(x)=∏i=1M‍pX(i)(x).
Typically, ϕX(x)∉𝒟2, so we will approximateϕX(x) by qX*(x)∈𝒟2 by projecting ϕX(x) onto the family 𝒟2:
(21)qX*(·)=argminqX∈𝒟2DKL(qx∥ϕx),
where DKL(·||·) denotes the Kullback Leibler (KL) divergence, defined as
(22)DKL(qX∥ϕX)=∫qX(x)logqX(x)ϕX(x)dx.Observe that all elements of𝒟2 are characterized by the parameters p≜[m1(a),m2(a),m1(b),m2(b),ρ,σ2]T and that the optimization (21) is therefore a six-dimensional problem over all possible p. The divergence DKL(q∥ϕ) for an arbitrary p∈ℝ4×ℝ+2 can be determined using Monte-Carlo integration as follows. We rewrite (22) as
(23)DKL(qx∥ϕx)=∫qX(x)f(x)dx,
where f(x)=logqX(x)-∑i=1M‍logpX(i)(x). By drawing R weighted samples {wk,xk}k=1R from qX(x) (e.g., through importance sampling), we can approximate (23) by
(24)DKL(qx∥ϕx)≈∑k=1Rwkf(xk).
Using this approximation, the six-dimensional optimization problem (21) is solved through gradient descent, similar to (19). The complexity of this operation scales as 𝒪(RM). The initial estimate of p is obtained through a set of heuristics: we first decide whether ϕX(x) can reasonably be represented by a distribution in 𝒟2. If not, the outgoing message is not computed. Otherwise, we use a geometric argument to find at most two midpoints. The initial estimates for ρ and σ2 are set to a small constant value.


### 4.5. Comparison of Message Representations

The complexities of the discretized, sample-based, and parametric message representations are compared in Table1.Table 1
Comparison of complexity of message representations, whereR is the number samples taken from each distribution and M is the number of distributions or messages.
Approach
Operation
Complexity
Value ofR
Discretized
Filtering
𝒪
(
R
2
)
Large
Discretized
Multiplication
𝒪
(
R
M
)
Large
Sample-based
Filtering
𝒪
(
R
)
Small
Sample-based
Multiplication
𝒪
(
R
2
M
)
Small
Parametric
Filtering
𝒪
(
R
)
Small
Parametric
Multiplication
𝒪
(
R
M
)
Small


## 4.1. Key Operations

In SPAWN, probabilistic information is exchanged and computed throughmessages. The manner in which these messages are represented for transmission between nodes and internal computation is closely related to the complexity and performance of the localization algorithm. In traditional communications problems, such as decoding, messages can be represented efficiently and exactly through, for instance, log-likelihood ratios [49]. In SPAWN, exact representation is impossible, so we must resort to different types of approximate message representations. Any representation must be able to capture the salient properties of the true message and must enable efficient computation of the key steps in SPAWN, namely, message filtering(A-1)-(A-2) and message multiplication(A-3). We consider three types of message representation: discretized, sample-based, and parametric.For convenience, we will introduce a set of new notations. For the filtering operation, the incoming message is denoted bypX(x), the filtering operation by h(x,y), and the outgoing message by pY(y), with
(1)pY(y)∝∫h(x,y)pX(x)dx.
For the multiplication operation, we assume M incoming messages pX(i)(x) (i=1,…,M) over a single variable X, and an outgoing message
(2)ϕX(x)∝∏i=1M‍pX(i)(x).
Note that (1) maps to (A-1) through the following association:
(3)x→xi(t-1),y→xi(t),h(x,y)→p(xi(t)∣xi(t-1))p(zi,self(t)∣xi(t-1),xi(t)),pX(x)→bi(t-1)(xi(t-1)).
Similarly, (1) maps to (A-2) through the following association:
(4)x→xj(t),y→xi(t),h(x,y)→p(zj→i(t)∣xi(t),xj(t)),pX(x)→bi(t)(xi(t)).


## 4.2. Discretized Message Representation

A naive but simple approach to represent a continuous distributionpX(x) is to uniformly discretize the domain of X, yielding a set of quantization points Q={x1,…,xR}. The distribution is then approximated as a finite list of values, {pX(xk)}k=1R. The filtering operation then becomes
(5)pY(yk)∝∑l=1Rh(xl,yk)pX(xl),
requiring 𝒪(R2) operations. The multiplication becomes
(6)ϕX(xk)∝∏i=1M‍pX(i)(xk),
requiring 𝒪(RM) operations. Because R scales exponentially with the dimensionality of X and a large number of points are required in every dimension to capture fine features of the messages, discretization is impractical for SPAWN in UWB localization.


## 4.3. Sample-Based Message Representation

A sample-based message representation, as used in [19, 50], overcomes the drawback of discretization by representing messages as samples, concentrated where the messages have significant mass. Before describing the detailed implementation of the filtering and multiplication operations, we give a brief overview of generic sampling techniques (see also [51, 52]) and kernel density estimation (KDE).


### 4.3.1. Background: Sampling and Kernel Density Estimation

We say that a list of samples with associated weights{xk,wk}k=1R is a representation for a distribution pX(x) if, for any integrable function g(x), we have the following approximation:
(7)I=∫g(x)pX(x)dx≈∑k=1Rwkg(xk).
Popular methods for obtaining the list of weighted samples include (i) direct sampling, where we draw R i.i.d. samples from pX(x), each with weight 1/R; and (ii) importance sampling, where we draw R i.i.d. samples from a distribution qX(x), with a support that includes the support of pX(x), and set the weight corresponding to sample xk as wk=pX(xk)/qX(xk). In both cases, it can easily be verified that the approximation is unbiased with mean I and variance that reduces with R (and that depends on qX(x), for importance sampling). Most importantly, the variance does not depend on the dimensionality of x.A variation of importance sampling that is not unbiased but that often has smaller variance is obtained by setting the weights as follows:wk∝pX(xk)/qX(xk), ∑k‍wk=1. This approach has the additional benefit that it does not require knowledge of the normalization constants of pX(x) or qX(x). A list of R equally weighted samples can be obtained from {xk,wk}k=1R through resampling, that is, by drawing (with repetition) R samples from the probability mass function defined by {xk,wk}k=1R.For numerical stability reasons, weights are often computed and stored in the logarithmic domain, that is,λk=logpX(xk)-logqX(xk). When the distributions involved contain exponentials or products, the log-domain representation is also computationally efficient. Operations such as additions can be evaluated efficiently in the log-domain as well, using the Jacobian logarithm [49, pages 90–94]. Once all R log-domain weights are computed, they are translated, exponentiated, and normalized: wk∝exp(λk-maxlλl).Given a sample representation{xk,wk}k=1R of a distribution pX(x), we obtain a kernel density estimate of pX(x) as
(8)p^X(x)=∑k=1RwkKσ(x-xk),
where Kσ(x) is the so-called kernel with bandwidth σ. The kernel is a symmetric distribution with a width parameter that is tuned through σ. For instance, a two-dimensional Gaussian kernel is given by
(9)Kσ(x)=12πσ2exp(-∥x∥22σ2).
While the choice of kernel affects the performance of the estimate to some limited extent (e.g., in an MMSE sense, where the error is ∫|pX(x)-p^X(x)|2pX(x)dx), the crucial parameter is the bandwidth σ, which needs to be estimated from the samples {xk,wk}k=1R. A large choice of σ makes p^X(x) smooth, but it may no longer capture the interesting features of pX(x). When σ is too small, p^X(x) may exhibit artificial structure not present in pX(x) [53].With this background in sampling techniques and KDE, we return to the problem at hand: filtering and multiplication of messages.


### 4.3.2. Message Filtering

We assume a message representation ofpX(x) as {xk,wk}k=1R and wish to obtain a message representation of pY(y)∝∫‍h(x,y)pX(x)dx. Let us interpret h(x,y) as a conditional distribution pY∣X(y∣x), up to some arbitrary constant. Suppose we can draw samples {[xk,yk],wk}k=1R~pY∣X(y∣x)pX(x); then {yk,wk}k=1R will form a sample representation of pY(y). Now the problem reverts to drawing samples from pY∣X(y∣x)pX(x). This can be accomplished as follows: first, for every sample xk, draw yk~qY∣X(y∣xk) from some distribution qY∣X(y∣xk). Second, set the weight of sample [xk,yk] as
(10)vk=wkpY∣X(yk∣xk)qY∣X(yk∣xk).
Finally, renormalize the weights vk to vk/∑lvl. The complexity of the filtering operation scales as 𝒪(R), a significant improvement from 𝒪(R2) for discretization. In addition, R can generally be much smaller in a particle-based representation.Let us consider some examples of the filtering operation in SPAWN.(i)
Mobility update(A-1): let pX(x) be the belief before movement (represented by {xk,wk}k=1R) and pY(y) the belief after movement. Assume that we are able to measure perfectly the distance traveled (given by zself), but have no information regarding the direction, and furthermore that the direction is chosen uniformly in (0,2π]. In that case,
(11)h(x,y)∝δ(zself-∥x-y∥),
where δ is a Dirac delta function, so that qY∣X(y∣xk)=pY∣X(y∣xk)∝δ(zself-∥xk-y∥) is a reasonable choice. For every xk, we can now draw values for y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and setting r=zself, leading to yk=xk+zself×[cosθksinθk]T, with vk∝wk.(ii)
Ranging update(A-2): let pX(x) be a message (represented by {xk,wk}k=1R) from a node with which we have performed ranging, resulting in a range estimate z. Let h(x,y)=pZ∣D(z∣d), where d=∥x-y∥. Note that pZ∣D(z∣d) is a likelihood function, since the measurement z is known. Assume that we have a model for the ranging performance in the form of distributions pZ∣D(z∣d) for any value of d. We then sample pY(y) as follows: for every xk, draw y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and r~qR∣Z(r∣z), for some well-chosen qR∣Z(r∣z)(e.g., a Gaussian distribution with mean equal to the distance estimate, z, and a standard deviation that is sufficiently large with respect to the standard deviation of p(z|d) for any d). The weights are set as
(12)vk=wkpZ∣D(z∣rk)qR∣Z(rkz).


### 4.3.3. Message Multiplication

Here we assume message representations{xk(i),wk(i)}k=1R for pX(i)(xk), i=1,…,M. In contrast to the discretization approach, we cannot directly compute ∏i=1M‍pX(i)(x) for arbitrary values of x. Rather, for every message pX(i)(xk), we create a KDE p^X(i)(x)=∑k=1R‍wkKσ(i)(x-xk(i)) with a Gaussian kernel and a bandwidth estimated using the methods from [53]. Suppose we now draw R samples from a distribution qX(x); then the weights are
(13)vk∝∏i=1M‍p^X(i)(xk)qX(xk),
which can be computed efficiently in the log-domain. A reasonable choice for qX(x) could be one of the incoming messages pX(i)(x) (e.g., the one with the smallest entropy) or a mixture of the incoming messages. The computational complexity of the message multiplication operation scales as 𝒪(MR2). This appears worse than the discretized case (complexity 𝒪(MR)), but note that R is much smaller for sample-based representations than for discretization (e.g., R=103 or R=104 for the sample-based representation compared to R=108 in the discretization).


## 4.3.1. Background: Sampling and Kernel Density Estimation

We say that a list of samples with associated weights{xk,wk}k=1R is a representation for a distribution pX(x) if, for any integrable function g(x), we have the following approximation:
(7)I=∫g(x)pX(x)dx≈∑k=1Rwkg(xk).
Popular methods for obtaining the list of weighted samples include (i) direct sampling, where we draw R i.i.d. samples from pX(x), each with weight 1/R; and (ii) importance sampling, where we draw R i.i.d. samples from a distribution qX(x), with a support that includes the support of pX(x), and set the weight corresponding to sample xk as wk=pX(xk)/qX(xk). In both cases, it can easily be verified that the approximation is unbiased with mean I and variance that reduces with R (and that depends on qX(x), for importance sampling). Most importantly, the variance does not depend on the dimensionality of x.A variation of importance sampling that is not unbiased but that often has smaller variance is obtained by setting the weights as follows:wk∝pX(xk)/qX(xk), ∑k‍wk=1. This approach has the additional benefit that it does not require knowledge of the normalization constants of pX(x) or qX(x). A list of R equally weighted samples can be obtained from {xk,wk}k=1R through resampling, that is, by drawing (with repetition) R samples from the probability mass function defined by {xk,wk}k=1R.For numerical stability reasons, weights are often computed and stored in the logarithmic domain, that is,λk=logpX(xk)-logqX(xk). When the distributions involved contain exponentials or products, the log-domain representation is also computationally efficient. Operations such as additions can be evaluated efficiently in the log-domain as well, using the Jacobian logarithm [49, pages 90–94]. Once all R log-domain weights are computed, they are translated, exponentiated, and normalized: wk∝exp(λk-maxlλl).Given a sample representation{xk,wk}k=1R of a distribution pX(x), we obtain a kernel density estimate of pX(x) as
(8)p^X(x)=∑k=1RwkKσ(x-xk),
where Kσ(x) is the so-called kernel with bandwidth σ. The kernel is a symmetric distribution with a width parameter that is tuned through σ. For instance, a two-dimensional Gaussian kernel is given by
(9)Kσ(x)=12πσ2exp(-∥x∥22σ2).
While the choice of kernel affects the performance of the estimate to some limited extent (e.g., in an MMSE sense, where the error is ∫|pX(x)-p^X(x)|2pX(x)dx), the crucial parameter is the bandwidth σ, which needs to be estimated from the samples {xk,wk}k=1R. A large choice of σ makes p^X(x) smooth, but it may no longer capture the interesting features of pX(x). When σ is too small, p^X(x) may exhibit artificial structure not present in pX(x) [53].With this background in sampling techniques and KDE, we return to the problem at hand: filtering and multiplication of messages.


## 4.3.2. Message Filtering

We assume a message representation ofpX(x) as {xk,wk}k=1R and wish to obtain a message representation of pY(y)∝∫‍h(x,y)pX(x)dx. Let us interpret h(x,y) as a conditional distribution pY∣X(y∣x), up to some arbitrary constant. Suppose we can draw samples {[xk,yk],wk}k=1R~pY∣X(y∣x)pX(x); then {yk,wk}k=1R will form a sample representation of pY(y). Now the problem reverts to drawing samples from pY∣X(y∣x)pX(x). This can be accomplished as follows: first, for every sample xk, draw yk~qY∣X(y∣xk) from some distribution qY∣X(y∣xk). Second, set the weight of sample [xk,yk] as
(10)vk=wkpY∣X(yk∣xk)qY∣X(yk∣xk).
Finally, renormalize the weights vk to vk/∑lvl. The complexity of the filtering operation scales as 𝒪(R), a significant improvement from 𝒪(R2) for discretization. In addition, R can generally be much smaller in a particle-based representation.Let us consider some examples of the filtering operation in SPAWN.(i)
Mobility update(A-1): let pX(x) be the belief before movement (represented by {xk,wk}k=1R) and pY(y) the belief after movement. Assume that we are able to measure perfectly the distance traveled (given by zself), but have no information regarding the direction, and furthermore that the direction is chosen uniformly in (0,2π]. In that case,
(11)h(x,y)∝δ(zself-∥x-y∥),
where δ is a Dirac delta function, so that qY∣X(y∣xk)=pY∣X(y∣xk)∝δ(zself-∥xk-y∥) is a reasonable choice. For every xk, we can now draw values for y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and setting r=zself, leading to yk=xk+zself×[cosθksinθk]T, with vk∝wk.(ii)
Ranging update(A-2): let pX(x) be a message (represented by {xk,wk}k=1R) from a node with which we have performed ranging, resulting in a range estimate z. Let h(x,y)=pZ∣D(z∣d), where d=∥x-y∥. Note that pZ∣D(z∣d) is a likelihood function, since the measurement z is known. Assume that we have a model for the ranging performance in the form of distributions pZ∣D(z∣d) for any value of d. We then sample pY(y) as follows: for every xk, draw y=xk+r×[cosθsinθ]T by drawing θ~𝒰(0,2π) and r~qR∣Z(r∣z), for some well-chosen qR∣Z(r∣z)(e.g., a Gaussian distribution with mean equal to the distance estimate, z, and a standard deviation that is sufficiently large with respect to the standard deviation of p(z|d) for any d). The weights are set as
(12)vk=wkpZ∣D(z∣rk)qR∣Z(rkz).


## 4.3.3. Message Multiplication

Here we assume message representations{xk(i),wk(i)}k=1R for pX(i)(xk), i=1,…,M. In contrast to the discretization approach, we cannot directly compute ∏i=1M‍pX(i)(x) for arbitrary values of x. Rather, for every message pX(i)(xk), we create a KDE p^X(i)(x)=∑k=1R‍wkKσ(i)(x-xk(i)) with a Gaussian kernel and a bandwidth estimated using the methods from [53]. Suppose we now draw R samples from a distribution qX(x); then the weights are
(13)vk∝∏i=1M‍p^X(i)(xk)qX(xk),
which can be computed efficiently in the log-domain. A reasonable choice for qX(x) could be one of the incoming messages pX(i)(x) (e.g., the one with the smallest entropy) or a mixture of the incoming messages. The computational complexity of the message multiplication operation scales as 𝒪(MR2). This appears worse than the discretized case (complexity 𝒪(MR)), but note that R is much smaller for sample-based representations than for discretization (e.g., R=103 or R=104 for the sample-based representation compared to R=108 in the discretization).


## 4.4. Parametric Message Representation

### 4.4.1. Choosing a Suitable Parameterization

From the previous section, it is clear that the bottleneck of the sample-based message representation lies in the message multiplication, which scales quadratically with the number of samples. An alternative approach is to represent each message as a set of parameters (e.g., a Gaussian distribution characterized by a mean and covariance matrix). In contrast to the sample-based message representation, which can represent messages of any shape, parametric representations must be specially tailored to the problem at hand. For example, single two-dimensional Gaussian parametric messages are utilized in [22] for localization with both range and angle measurements. Our choice of parametric message is based on the following observations.(i)
For the filtering operation with a two-dimensional Gaussian inputpX(x), the output pY(y) can be approximated by a circular distribution with the same mean for both the mobility update (A-1) and the ranging update (A-2).(ii)
Multiplying Gaussian distributions yield a Gaussian distribution.(iii)
The multiplication of multiple circular distributions can be approximated by a Gaussian distribution or a mixture of Gaussian distributions.We will use as a basic building block the following distribution in two dimensions:(14)𝒟(x;m1,m2,σ2,ρ)=1C(σ2,ρ)exp{-[(x1-m1)2+(x2-m2)2-ρ]22σ2},
where [m1,m2] is the midpoint of the distribution, ρ is the radius, σ2 is the variance, and C(σ2,ρ) is a normalization constant equal to
(15)C(σ2,ρ)=2πσ2[exp(-ρ22σ2)+122πρ2σ2(1+erfρ22σ2)].
As a special case, we note that, when ρ=0, (14) reverts to a two-dimensional Gaussian. Moreover, we will represent all messages as a mixture of two distributions of the type (14), so that
(16)pX(x)=12𝒟(x;m1(a),m2(a),σ2,ρ)+12𝒟(x;m1(b),m2(b),σ2,ρ),
which can be represented by the six-dimensional vector [m1(a),m2(a),m1(b),m2(b),ρ,σ2]. We will denote the family of distributions of the form (16) by 𝒟2. Note that it is trivial to extend this distribution, which is designed for two-dimensional localization systems, for use in three-dimensional systems. Before we describe the message filtering and message multiplication operations, let us first show how the parameters of (14) can be estimated from a list of samples.


### 4.4.2. ML Estimation of the Parametersm1,m2,σ2,ρ.

Given a list of samples{xk}k=1R, we can estimate the parameters [m1,m2,σ2,ρ] as follows. The midpoint m=[m1,m2] is estimated by
(17)m^=1R∑k=1Rxk.
To find the radius ρ and variance σ2 of 𝒟-distribution, we use maximum likelihood (ML) estimation, assuming the R samples are independent. Introducing α=[ρσ2]T, we find that
(18)α^ML=argmaxα∑k=1RlogpX(xk;m^,α)︸Λ(α).
Treating the log-likelihood function (LLF) Λ(α) as an objective function, we find its maximum through the gradient ascent algorithm
(19)α^(n+1)=α^(n)-ε∇αΛ(α)|α=α^(n),
where ε is a suitably small step size, and the gradient vector can be approximated using finite differences. To initialize (19), we consider two initial estimates for σ2 and ρ: one assuming ρ≪σ and a second assuming ρ≫σ. The LLF is evaluated for both preliminary solutions, and the one with the largest log-likelihood is used as the initial estimate in (19).


### 4.4.3. Message Filtering

To perform message filtering, we use the fact that sample-based message filtering is a low-complexity operation. We decomposepX(x), represented in parametric form, into its two mixture components. From each component, we draw R/2 samples and perform sample-based message filtering, as outlined in Section 4.3.2. We can then estimate the new 𝒟-parameters for each mixture component using the ML method described above. We thus have pY(y) in parametric form. The complexity of this operation scales as 𝒪(R).


### 4.4.4. Message Multiplication

The motivation for using the parametric message representation is to avoid the complexity associated with sample-based message multiplication. GivenM distributions px(i)(x)∈𝒟2, our goal is to compute
(20)ϕX(x)=∏i=1M‍pX(i)(x).
Typically, ϕX(x)∉𝒟2, so we will approximateϕX(x) by qX*(x)∈𝒟2 by projecting ϕX(x) onto the family 𝒟2:
(21)qX*(·)=argminqX∈𝒟2DKL(qx∥ϕx),
where DKL(·||·) denotes the Kullback Leibler (KL) divergence, defined as
(22)DKL(qX∥ϕX)=∫qX(x)logqX(x)ϕX(x)dx.Observe that all elements of𝒟2 are characterized by the parameters p≜[m1(a),m2(a),m1(b),m2(b),ρ,σ2]T and that the optimization (21) is therefore a six-dimensional problem over all possible p. The divergence DKL(q∥ϕ) for an arbitrary p∈ℝ4×ℝ+2 can be determined using Monte-Carlo integration as follows. We rewrite (22) as
(23)DKL(qx∥ϕx)=∫qX(x)f(x)dx,
where f(x)=logqX(x)-∑i=1M‍logpX(i)(x). By drawing R weighted samples {wk,xk}k=1R from qX(x) (e.g., through importance sampling), we can approximate (23) by
(24)DKL(qx∥ϕx)≈∑k=1Rwkf(xk).
Using this approximation, the six-dimensional optimization problem (21) is solved through gradient descent, similar to (19). The complexity of this operation scales as 𝒪(RM). The initial estimate of p is obtained through a set of heuristics: we first decide whether ϕX(x) can reasonably be represented by a distribution in 𝒟2. If not, the outgoing message is not computed. Otherwise, we use a geometric argument to find at most two midpoints. The initial estimates for ρ and σ2 are set to a small constant value.


## 4.4.1. Choosing a Suitable Parameterization

From the previous section, it is clear that the bottleneck of the sample-based message representation lies in the message multiplication, which scales quadratically with the number of samples. An alternative approach is to represent each message as a set of parameters (e.g., a Gaussian distribution characterized by a mean and covariance matrix). In contrast to the sample-based message representation, which can represent messages of any shape, parametric representations must be specially tailored to the problem at hand. For example, single two-dimensional Gaussian parametric messages are utilized in [22] for localization with both range and angle measurements. Our choice of parametric message is based on the following observations.(i)
For the filtering operation with a two-dimensional Gaussian inputpX(x), the output pY(y) can be approximated by a circular distribution with the same mean for both the mobility update (A-1) and the ranging update (A-2).(ii)
Multiplying Gaussian distributions yield a Gaussian distribution.(iii)
The multiplication of multiple circular distributions can be approximated by a Gaussian distribution or a mixture of Gaussian distributions.We will use as a basic building block the following distribution in two dimensions:(14)𝒟(x;m1,m2,σ2,ρ)=1C(σ2,ρ)exp{-[(x1-m1)2+(x2-m2)2-ρ]22σ2},
where [m1,m2] is the midpoint of the distribution, ρ is the radius, σ2 is the variance, and C(σ2,ρ) is a normalization constant equal to
(15)C(σ2,ρ)=2πσ2[exp(-ρ22σ2)+122πρ2σ2(1+erfρ22σ2)].
As a special case, we note that, when ρ=0, (14) reverts to a two-dimensional Gaussian. Moreover, we will represent all messages as a mixture of two distributions of the type (14), so that
(16)pX(x)=12𝒟(x;m1(a),m2(a),σ2,ρ)+12𝒟(x;m1(b),m2(b),σ2,ρ),
which can be represented by the six-dimensional vector [m1(a),m2(a),m1(b),m2(b),ρ,σ2]. We will denote the family of distributions of the form (16) by 𝒟2. Note that it is trivial to extend this distribution, which is designed for two-dimensional localization systems, for use in three-dimensional systems. Before we describe the message filtering and message multiplication operations, let us first show how the parameters of (14) can be estimated from a list of samples.


## 4.4.2. ML Estimation of the Parametersm1,m2,σ2,ρ.

Given a list of samples{xk}k=1R, we can estimate the parameters [m1,m2,σ2,ρ] as follows. The midpoint m=[m1,m2] is estimated by
(17)m^=1R∑k=1Rxk.
To find the radius ρ and variance σ2 of 𝒟-distribution, we use maximum likelihood (ML) estimation, assuming the R samples are independent. Introducing α=[ρσ2]T, we find that
(18)α^ML=argmaxα∑k=1RlogpX(xk;m^,α)︸Λ(α).
Treating the log-likelihood function (LLF) Λ(α) as an objective function, we find its maximum through the gradient ascent algorithm
(19)α^(n+1)=α^(n)-ε∇αΛ(α)|α=α^(n),
where ε is a suitably small step size, and the gradient vector can be approximated using finite differences. To initialize (19), we consider two initial estimates for σ2 and ρ: one assuming ρ≪σ and a second assuming ρ≫σ. The LLF is evaluated for both preliminary solutions, and the one with the largest log-likelihood is used as the initial estimate in (19).


## 4.4.3. Message Filtering

To perform message filtering, we use the fact that sample-based message filtering is a low-complexity operation. We decomposepX(x), represented in parametric form, into its two mixture components. From each component, we draw R/2 samples and perform sample-based message filtering, as outlined in Section 4.3.2. We can then estimate the new 𝒟-parameters for each mixture component using the ML method described above. We thus have pY(y) in parametric form. The complexity of this operation scales as 𝒪(R).


## 4.4.4. Message Multiplication

The motivation for using the parametric message representation is to avoid the complexity associated with sample-based message multiplication. GivenM distributions px(i)(x)∈𝒟2, our goal is to compute
(20)ϕX(x)=∏i=1M‍pX(i)(x).
Typically, ϕX(x)∉𝒟2, so we will approximateϕX(x) by qX*(x)∈𝒟2 by projecting ϕX(x) onto the family 𝒟2:
(21)qX*(·)=argminqX∈𝒟2DKL(qx∥ϕx),
where DKL(·||·) denotes the Kullback Leibler (KL) divergence, defined as
(22)DKL(qX∥ϕX)=∫qX(x)logqX(x)ϕX(x)dx.Observe that all elements of𝒟2 are characterized by the parameters p≜[m1(a),m2(a),m1(b),m2(b),ρ,σ2]T and that the optimization (21) is therefore a six-dimensional problem over all possible p. The divergence DKL(q∥ϕ) for an arbitrary p∈ℝ4×ℝ+2 can be determined using Monte-Carlo integration as follows. We rewrite (22) as
(23)DKL(qx∥ϕx)=∫qX(x)f(x)dx,
where f(x)=logqX(x)-∑i=1M‍logpX(i)(x). By drawing R weighted samples {wk,xk}k=1R from qX(x) (e.g., through importance sampling), we can approximate (23) by
(24)DKL(qx∥ϕx)≈∑k=1Rwkf(xk).
Using this approximation, the six-dimensional optimization problem (21) is solved through gradient descent, similar to (19). The complexity of this operation scales as 𝒪(RM). The initial estimate of p is obtained through a set of heuristics: we first decide whether ϕX(x) can reasonably be represented by a distribution in 𝒟2. If not, the outgoing message is not computed. Otherwise, we use a geometric argument to find at most two midpoints. The initial estimates for ρ and σ2 are set to a small constant value.


## 4.5. Comparison of Message Representations

The complexities of the discretized, sample-based, and parametric message representations are compared in Table1.Table 1
Comparison of complexity of message representations, whereR is the number samples taken from each distribution and M is the number of distributions or messages.
Approach
Operation
Complexity
Value ofR
Discretized
Filtering
𝒪
(
R
2
)
Large
Discretized
Multiplication
𝒪
(
R
M
)
Large
Sample-based
Filtering
𝒪
(
R
)
Small
Sample-based
Multiplication
𝒪
(
R
2
M
)
Small
Parametric
Filtering
𝒪
(
R
)
Small
Parametric
Multiplication
𝒪
(
R
M
)
Small


## 5. Performance Analysis

In this section, we compare the performance of the SPAWN algorithm with sample-based versus parametric message representation in a simulated wireless network. We also analyze the use of different subsets of information in the algorithm and its effect on localization performance.


### 5.1. Simulation Setup and Performance Measures

We simulate a large-scale ultra-wide bandwidth (UWB) network in a 100 m × 100 m homogeneous environment, with 100 uniformly distributed agents and 13 fixed anchors in a grid configuration. Each node is able to measure its range to other nodes within 20 meters. The simulated ranging measurements are independently drawn from the UWB ranging model developed in [54]. The model, based on data collected in a variety of indoor scenarios, consists of three component Gaussian densities, where the mean and variance of each component are experimentally determined functions of the true distance between the ranging nodes. To decouple the effect of mobility with the message representation, we consider a single time slot, where every agent has a uniform a priori distribution over the environment ℰ. SPAWN was run for Nit=20 iterations, though convergence was generally achieved well before 10 iterations. For the sample-based representation, the number of samples is set to R=2048 unless otherwise stated.We quantify localization performance using the complementary cumulative distribution function (CCDF) of the localization errore=∥xi-x^i∥, where x^i is the estimated location of node i, taken as the mean of the belief, similar to [18]. To estimate the CCDF, we consider 50 random network topologies and collect position estimates at every iteration for every agent. Note that a CCDF of 0.01 at an error of, say, e=1 m means that 99% of the nodes have an error less than 1 meter.


### 5.2. Cooperation with Censoring

In Section3, we considered processing messages between all neighboring pairs of nodes. However, information from neighbors may not always be useful: (i) when the receiving node's belief is already very informative (e.g., concentrated around the mean); or (ii) when the transmitting node's belief is very uninformative. To better understand how much cooperative information is beneficial to localization, we will consider varying the subset of nodes that broadcast and update their location beliefs at each iteration. We distinguish between these subsets by the level of cooperative information they induce in the algorithm. The level of cooperative information indicates how each node utilizes information from its neighbors at each iteration.We introduce the following terminology: a distribution is said to be “sufficiently informative” when 95% of the probability mass is located within 2 m of the mean; a node becomes avirtual anchor when its belief is sufficiently informative; a virtual bianchor is a node with a bimodal belief, with each mode being sufficiently informative; a node that is neither a virtual anchor nor a virtual bianchor will be called a blind agent. We are now ready to introduce four levels of cooperative information at each iteration.(i)
Level 1 (L1): virtual anchors broadcast their beliefs, while all other nodes censor their belief broadcast. Virtual anchors do not update their beliefs.(ii)
Level 2 (L2): virtual anchors and virtual bianchors broadcast their beliefs, while blind nodes censor their belief broadcast. Virtual anchors do not update their beliefs.(iii)
Level 3 (L3): all nodes broadcast their beliefs. Virtual anchors do not update their beliefs.(iv)
Level 4 (L4): all nodes broadcast their beliefs. All nodes update their beliefs. In terms of cooperation, note that L4 utilizes more cooperative information than L3, L3 utilizes more cooperative information than L2, and L2 utilizes more cooperative information than L1. In this sense, the levels of cooperative information are strict subsets.From previous sections, we know that the algorithm complexity scales linearly inM, the number of incoming messages in the multiplication operation. Hence, the level of cooperative information directly affects the algorithm's computational cost, with lower levels requiring less computation.


### 5.3. Numerical Results

We now examine how localization performance varies with the algorithm parameters. In particular, numerical results show the effect of message representation (sample-based or parametric) and level of cooperative information (L1, L2, L3, or L4) on the CCDF of the localization error.We first consider the localization performance as a function of the number of samplesR and level of cooperative information. Figure 1 displays the CCDF at e=1 m after 10 iterations. As expected, for any level of cooperative information, the CCDF decreases as the number of samples is increased. However, the decrease in CCDF comes with a cost in computation time; as R is increased, the per-node complexity increases quadratically. Figure 1 also shows that levels L1, L2, and L4 are not as sensitive to R as L3 and that each generally outperforms L3. This effect is particularly pronounced when R is small. L3 broadcasts more complex distributions than L2 and L1, and these elaborate distributions are not accurately represented with a small number of samples.Figure 1
The effect ofR and the level of cooperative information on the CCDF at e=1 meter. All results are after Nit=10 iterations.Secondly, we investigate level of cooperative information and its effect on localization performance, with numerical results represented in Figures2 and 3, after Nit=20 iterations. Note that each curve exhibits a “floor” because there is always some subset of nodes that have insufficient information to localize without ambiguity. This may be due to lack of connectivity or large flip ambiguities. Let us focus on the sample-based representation in Figure 2 and consider the effect of the level of cooperative information on localization performance. In general, L4 has the best performance in terms of accuracy and floor. Intuitively, one might expect L3 to have the next best performance, followed by L2, and then L1. However, Figure 2 demonstrates that in some cases L3 has poorer accuracy than L2 and a similar floor. This effect can be explained as follows. Agents that do not become virtual anchors within Nit=20 tend to have large localization errors, creating a floor. Such agents comprise 1.7% of the total nodes for L1 and 0.3% for both L2 and L3. Since L2 and L3 have a similar fraction of agents that do not become virtual anchors, they have similar floors. In addition, the accuracy of beliefs belonging to agents that have become virtual anchors turns out to be highest for L1, followed by L2, and then L3. This is because L3 uses less reliable information than L2, which in turn is less reliable than L1. The final CCDF depends both on the fraction of virtual anchors (lowest for L1) and the accuracy of those virtual anchors (highest for L1). Note that we cannot compare L4 in this context, since there is no concept of a virtual anchor in L4.Figure 2
Effect of the level of cooperative information on localization accuracy for the distributed sample-based and parametric algorithms, whenNit=20, and averaged over 50 networks.Figure 3
Effect of the level of cooperative information on the convergence speed fore=1 m, for the distributed sample-based and parametric algorithms. Results are averaged over 50 networks.We now move on to the parametric representation, still in Figure2. We observe that L4 has the lowest overall CCDF for any e, for both types of message representation. For the parametric messages, the differences among different levels of cooperative information are smaller, and we generally obtain better performance (for e<1 m) compared to sample-based messages.Finally, in Figure3, we evaluate the convergence speed of the different message representations and levels of cooperation, for a fixed error of 1 meter. We see that the parametric messages generally lead to faster convergence and lower CCDF than their sample-based counterparts. Levels L2, L3, and L4 all converge in around 5 iterations with a final CCDF at e=1 m of around 0.01 for the parametric representation. Our results show that more cooperative information leads to faster improvement in terms of accuracy. The lowest level of cooperative information, L1, is consistently slower to converge and less accurate. However, higher levels of cooperative information also require the computation and representation of more complicated distributions. As a possible consequence, convergence issues may occur for levels L3 and L4. We also see that the parametric message representation performs approximately equal to or better than the sample-based messages in terms of both convergence and accuracy, while requiring much less execution time. Overall, parametric message representations yield a better performance/complexity tradeoff. This is due to the fact that the parametric distributions are well tailored to the localization problem and the homogeneous simulation environment.


## 5.1. Simulation Setup and Performance Measures

We simulate a large-scale ultra-wide bandwidth (UWB) network in a 100 m × 100 m homogeneous environment, with 100 uniformly distributed agents and 13 fixed anchors in a grid configuration. Each node is able to measure its range to other nodes within 20 meters. The simulated ranging measurements are independently drawn from the UWB ranging model developed in [54]. The model, based on data collected in a variety of indoor scenarios, consists of three component Gaussian densities, where the mean and variance of each component are experimentally determined functions of the true distance between the ranging nodes. To decouple the effect of mobility with the message representation, we consider a single time slot, where every agent has a uniform a priori distribution over the environment ℰ. SPAWN was run for Nit=20 iterations, though convergence was generally achieved well before 10 iterations. For the sample-based representation, the number of samples is set to R=2048 unless otherwise stated.We quantify localization performance using the complementary cumulative distribution function (CCDF) of the localization errore=∥xi-x^i∥, where x^i is the estimated location of node i, taken as the mean of the belief, similar to [18]. To estimate the CCDF, we consider 50 random network topologies and collect position estimates at every iteration for every agent. Note that a CCDF of 0.01 at an error of, say, e=1 m means that 99% of the nodes have an error less than 1 meter.


## 5.2. Cooperation with Censoring

In Section3, we considered processing messages between all neighboring pairs of nodes. However, information from neighbors may not always be useful: (i) when the receiving node's belief is already very informative (e.g., concentrated around the mean); or (ii) when the transmitting node's belief is very uninformative. To better understand how much cooperative information is beneficial to localization, we will consider varying the subset of nodes that broadcast and update their location beliefs at each iteration. We distinguish between these subsets by the level of cooperative information they induce in the algorithm. The level of cooperative information indicates how each node utilizes information from its neighbors at each iteration.We introduce the following terminology: a distribution is said to be “sufficiently informative” when 95% of the probability mass is located within 2 m of the mean; a node becomes avirtual anchor when its belief is sufficiently informative; a virtual bianchor is a node with a bimodal belief, with each mode being sufficiently informative; a node that is neither a virtual anchor nor a virtual bianchor will be called a blind agent. We are now ready to introduce four levels of cooperative information at each iteration.(i)
Level 1 (L1): virtual anchors broadcast their beliefs, while all other nodes censor their belief broadcast. Virtual anchors do not update their beliefs.(ii)
Level 2 (L2): virtual anchors and virtual bianchors broadcast their beliefs, while blind nodes censor their belief broadcast. Virtual anchors do not update their beliefs.(iii)
Level 3 (L3): all nodes broadcast their beliefs. Virtual anchors do not update their beliefs.(iv)
Level 4 (L4): all nodes broadcast their beliefs. All nodes update their beliefs. In terms of cooperation, note that L4 utilizes more cooperative information than L3, L3 utilizes more cooperative information than L2, and L2 utilizes more cooperative information than L1. In this sense, the levels of cooperative information are strict subsets.From previous sections, we know that the algorithm complexity scales linearly inM, the number of incoming messages in the multiplication operation. Hence, the level of cooperative information directly affects the algorithm's computational cost, with lower levels requiring less computation.


## 5.3. Numerical Results

We now examine how localization performance varies with the algorithm parameters. In particular, numerical results show the effect of message representation (sample-based or parametric) and level of cooperative information (L1, L2, L3, or L4) on the CCDF of the localization error.We first consider the localization performance as a function of the number of samplesR and level of cooperative information. Figure 1 displays the CCDF at e=1 m after 10 iterations. As expected, for any level of cooperative information, the CCDF decreases as the number of samples is increased. However, the decrease in CCDF comes with a cost in computation time; as R is increased, the per-node complexity increases quadratically. Figure 1 also shows that levels L1, L2, and L4 are not as sensitive to R as L3 and that each generally outperforms L3. This effect is particularly pronounced when R is small. L3 broadcasts more complex distributions than L2 and L1, and these elaborate distributions are not accurately represented with a small number of samples.Figure 1
The effect ofR and the level of cooperative information on the CCDF at e=1 meter. All results are after Nit=10 iterations.Secondly, we investigate level of cooperative information and its effect on localization performance, with numerical results represented in Figures2 and 3, after Nit=20 iterations. Note that each curve exhibits a “floor” because there is always some subset of nodes that have insufficient information to localize without ambiguity. This may be due to lack of connectivity or large flip ambiguities. Let us focus on the sample-based representation in Figure 2 and consider the effect of the level of cooperative information on localization performance. In general, L4 has the best performance in terms of accuracy and floor. Intuitively, one might expect L3 to have the next best performance, followed by L2, and then L1. However, Figure 2 demonstrates that in some cases L3 has poorer accuracy than L2 and a similar floor. This effect can be explained as follows. Agents that do not become virtual anchors within Nit=20 tend to have large localization errors, creating a floor. Such agents comprise 1.7% of the total nodes for L1 and 0.3% for both L2 and L3. Since L2 and L3 have a similar fraction of agents that do not become virtual anchors, they have similar floors. In addition, the accuracy of beliefs belonging to agents that have become virtual anchors turns out to be highest for L1, followed by L2, and then L3. This is because L3 uses less reliable information than L2, which in turn is less reliable than L1. The final CCDF depends both on the fraction of virtual anchors (lowest for L1) and the accuracy of those virtual anchors (highest for L1). Note that we cannot compare L4 in this context, since there is no concept of a virtual anchor in L4.Figure 2
Effect of the level of cooperative information on localization accuracy for the distributed sample-based and parametric algorithms, whenNit=20, and averaged over 50 networks.Figure 3
Effect of the level of cooperative information on the convergence speed fore=1 m, for the distributed sample-based and parametric algorithms. Results are averaged over 50 networks.We now move on to the parametric representation, still in Figure2. We observe that L4 has the lowest overall CCDF for any e, for both types of message representation. For the parametric messages, the differences among different levels of cooperative information are smaller, and we generally obtain better performance (for e<1 m) compared to sample-based messages.Finally, in Figure3, we evaluate the convergence speed of the different message representations and levels of cooperation, for a fixed error of 1 meter. We see that the parametric messages generally lead to faster convergence and lower CCDF than their sample-based counterparts. Levels L2, L3, and L4 all converge in around 5 iterations with a final CCDF at e=1 m of around 0.01 for the parametric representation. Our results show that more cooperative information leads to faster improvement in terms of accuracy. The lowest level of cooperative information, L1, is consistently slower to converge and less accurate. However, higher levels of cooperative information also require the computation and representation of more complicated distributions. As a possible consequence, convergence issues may occur for levels L3 and L4. We also see that the parametric message representation performs approximately equal to or better than the sample-based messages in terms of both convergence and accuracy, while requiring much less execution time. Overall, parametric message representations yield a better performance/complexity tradeoff. This is due to the fact that the parametric distributions are well tailored to the localization problem and the homogeneous simulation environment.


## 6. Conclusions and Extensions

In this paper, we considered different message representations for Bayesian cooperative localization in wireless networks: a generic sample-based representation and a tailored parametric representation. We used experimentally derived UWB ranging models to evaluate the performance of SPAWN as a function of message representation and level of cooperative information. Our results show that the tradeoffs between message representation, cooperative information, localization accuracy, and algorithm convergence are not straightforward and should be tailored to the scenario.Through large-scale network simulations, we demonstrated that more cooperative information may improve localization accuracy but also increase the complexity of messages. Higher levels of cooperative information do not always correspond to an improvement in localization accuracy or convergence rate. As complicated distributions associated with location-uncertain nodes are computed and transmitted, the resulting increases in computational complexity and signal interference can actually reduce localization performance. It may therefore be advantageous to broadcast only confident information in cooperative localization networks, especially considering the resources saved by a node censorship policy.We also demonstrated that though parametric messages have less representational flexibility, they can outperform nonparametric message representation at a much lower computational cost. In our simulations, the parametric representation achieved a lower probability of outage for errors under 1 meter while converging in equal or fewer iterations than the sample-based representation. Clearly, a parametric representation well tailored to the localization scenario is desirable in terms of both resource efficiency and localization accuracy.The use of parametric distributions for localization can be extended to (i) different ranging models; (ii) different types of measurements; (iii) more general scenarios. In terms of ranging models, the proposed distributions can be applied as long as typical distributions in SPAWN roughly resemble a distribution in the𝒟2-family. Note that a Gaussian ranging error satisfies this criteria, as would many other, more realistic, models. Other models, such as those derived from received signal strength, will require different types of parametric distributions. The same comment applies to the use of different types of measurements. For instance, with angle-of-arrival measurements, the parametric distributions should include a collection of linear distributions. Finally, more general scenarios may require tailor-made distributions. With NLOS measurements that can be modeled as biased Gaussians [20], for example, mixtures of 𝒟2 distributions would easily accommodate LOS/NLOS propagation, without relying on explicit NLOS identification.



---

*Source: 281592-2012-09-12.xml*