ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

1

Multiple scan data association by
convex variational inference (extended version)

arXiv:1607.07942v2 [] 23 Jan 2018

Jason L. Williams, Senior Member, IEEE and Roslyn A. Lau, Student Member, IEEE

Abstract—Data association, the reasoning over correspondence
between targets and measurements, is a problem of fundamental
importance in target tracking. Recently, belief propagation (BP)
has emerged as a promising method for estimating the marginal
probabilities of measurement to target association, providing fast,
accurate estimates. The excellent performance of BP in the particular formulation used may be attributed to the convexity of the
underlying free energy which it implicitly optimises. This paper
studies multiple scan data association problems, i.e., problems
that reason over correspondence between targets and several sets
of measurements, which may correspond to different sensors or
different time steps. We find that the multiple scan extension of
the single scan BP formulation is non-convex and demonstrate
the undesirable behaviour that can result. A convex free energy
is constructed using the recently proposed fractional free energy
(FFE). A convergent, BP-like algorithm is provided for the single
scan FFE, and employed in optimising the multiple scan free
energy using primal-dual coordinate ascent. Finally, based on a
variational interpretation of joint probabilistic data association
(JPDA), we develop a sequential variant of the algorithm that is
similar to JPDA, but retains consistency constraints from prior
scans. The performance of the proposed methods is demonstrated
on a bearings only target localisation problem.

I. I NTRODUCTION
Multiple target tracking is complicated by data association,
the unknown correspondence between measurements and
targets. The classical problem arises under the assumption
that measurements are received in scans (i.e., a collection of
measurements made at a single time), and that within each
scan, each target corresponds to at most one measurement, and
each measurement corresponds to at most one target.
Techniques for addressing data association may be classified
as either single scan (considering a single scan of data at a
time) or multiple scan (simultaneously considering multiple
scans), and as either maximum a posteriori (MAP) (finding the
most likely correspondence), or marginal-based (calculating the
full marginal distribution for each target). Common methods
include:
1) Global nearest neighbour (GNN), e.g., [1], is a single
scan MAP method, which finds the MAP correspondence
in the latest scan, and proceeds to the next scan assuming
that correspondence was correct
J. L. Williams (e-mail: Jason.Williams@dst.defence.gov.au) is with the
National Security, Intelligence, Surveillance and Reconnaissance Division,
Defence Science and Technology Group, Australia, and the School of Electrical
Engineering and Computer Science, Queensland University of Technology,
Australia. R. A. Lau (e-mail: Roslyn.Lau@dst.defence.gov.au) is with the
Maritime Division, Defence Science and Technology Group, Australia and
the Research School of Computer Science, Australian National University,
Australia.

2) Multiple hypothesis tracking (MHT) [2]–[4] is a multiple
scan MAP method, which in each scan seeks to find the
MAP correspondence over a recent history of scans
3) Joint probabilistic data association (JPDA) [5] is a
single scan marginal-based method, which calculates
the marginal distribution of each target, and proceeds by
approximating the joint distribution as the product of its
marginals
Classical JPDA additionally approximates the distribution of
each target as a moment-matched Gaussian distribution; in
this paper, we use the term JPDA more generally to refer to
the approach that retains the full marginal distribution of each
target in a manner similar to [6]–[10].
Compared to GNN, MHT and JPDA, multiple scan variants
of JPDA, e.g., [11], have received less attention. One may posit
that this is due to their formidable computational complexity:
While there exist fast approximations to the multiple scan
MAP problem such as Lagrangian relaxation [12], [13], no
such equivalents have existed for either single scan or multiple
scan JPDA.
Variational inference (e.g., [14], [15]) describes the collection
of methods that use optimisation (or calculus of variations)
to approximate difficult inference problems in probabilistic
graphical models (PGMs). Methods within this framework
include belief propagation (BP) [16], [17], mean field (MF)
[18],1 hybrid BP/MF approaches [19], [20], tree-reweighted
sum product (TRSP) [21] and norm-product BP (NPBP) [22].
Excellent performance has been demonstrated in a variety of
problems, typified by the recognition that turbo coding is an
instance of BP [23].
Variational inference was first applied to data association in
[24]–[26], addressing the problem of distributed tracking using
wireless sensor networks. The problem was formulated with
vertices corresponding to targets and sensors, where sensor
nodes represent the joint association of all sensor measurements
in the scan to targets. A related sensor network application
was studied in [27].
In contrast to these methods, which hypothesise the joint
association of all sensor measurements via a single variable,
the approach in [28]–[32] formulates the single scan problem
in terms of a bipartite graph, where vertices hypothesise
the measurement associated with a particular target, or the
target associated with a particular measurement. Empirically, it
was found that BP converges reliably and produces excellent
estimates of both the marginal association probabilities [30],
1 Mean field is also referred to as variational Bayes; following [14], we
use the term variational inference more generally, to refer to the entire family
of optimisation-based methods.

2

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

[32] and the partition function [31]. Convergence of BP in the
II. BACKGROUND AND M ODEL
two related formulations was proven in [33], [34].
We use the abbreviated notation p(x), p(z|x), etc, to
In [35] it was shown that, when correctly parameterised, represent the probability density function (PDF) or point mass
the variational inference problem that underlies the bipartite function (PMF) of the random variable corresponding to the
formulation is convex. This may be understood to be the source value x, z conditioned on x, etc.
of the empirically observed robustness of the approach when
applied to the bipartite model. For example, it was shown in
[36] that the approximations produced by BP tend to be either A. Probabilistic graphical models and variational inference
very good, or very bad. This may be understood through the PGMs [14], [15], [47] aim to represent and manipulate the
intuition that BP converges to a “good” approximation if the joint probability distributions of many variables efficiently
objective function that it implicitly optimises is locally convex by exploiting factorisation. The Kalman filter [48] and the
in the area between the starting point and the optimal solution, hidden Markov model (HMM) [49] are two examples of
and a “bad” approximation if it is non-convex. Accordingly, if algorithms that exploit sparsity of a particular kind (i.e., a
the problem is globally convex, a major source of degenerate Markov chain) to efficiently conduct inference on systems
cases is eliminated.
involving many random variables. Inference methods based on
Various applications merging the BP formulation in [32] with the PGM framework generalise these algorithms to a wider
MF (as proposed in [20]) and expectation maximisation (EM) variety of state spaces and dependency structures.
PGMs have been developed for undirected graphical models
were examined in [37]–[41]. PGM methods provide a path for
extending the bipartite model to multiple scan problems; this (Markov random fields), directed graphical models (Bayes
has been studied in [42], [43], which extends the single scan BP nets) and factor graphs. In this work we consider a subclass
formulation of [32] to multiple scans using a restricted message of pairwise undirected models, involving vertices (i.e., random
passing schedule, rather than optimising the variational problem variables) v ∈ V, and edges (i.e., dependencies) e ∈ E ⊂ V ×V,
to convergence. An alternative PGM formulation of the same and where the joint distribution can be written as:2
Y
Y
problem was utilised in [44] for the purpose of parameter
p(x
)
∝
ψ
(x
)
ψi,j (xi , xj ).
V
v
v
identification. However, these approaches lose the convexity
v∈V
(i,j)∈E
property of the single scan bipartite formulation, which is
understood to be the source of the robustness of this special As an example, a Markov chain involving variables
case.
(x1 , . . . , xT ) may be formulated by setting ψ1 (x1 ) = p(x1 ),
ψt (xt ) = 1 for t > 1, and edges ψt−1,t (xt−1 , xt ) =
p(xt |xt−1 ), t ∈ {2, . . . , T } representing the Markov transition
A. Contributions
kernels, although other formulations are possible.
Exact inference can be conducted on tree-structured graphs
This paper addresses the multiple scan data association
using belief propagation (BP), which operates by passing mesproblem using a convexification of the multiple scan model.
sages between neighbouring vertices. We denote by µi→j (xj )
We consider methods that optimise the variational problem to
the message sent from vertex i ∈ V to vertex j ∈ V where
convergence, rather than using restricted message schedules as
(i, j) ∈ E. The iterative update equations are then:
proposed in [42], [43]. A preliminary look at the convergence
X
Y
and performance of BP in multiple scan problems was also
µi→j (xj ) ∝
ψi,j (xi , xj )ψi (xi )
µj 0 →i (xi ).
included in [30]. In section IV, we perform a thorough
xi
(j 0 ,i)∈E,j 0 6=j
evaluation of these methods on a bearings only localisation
(1)
problem, and find that each of these methods can, in a This is also known as the sum-product algorithm. If the
challenging environment, give vanishingly small likelihood summations are replaced with maximisations, then we arrive
to the true solution in a significant portion of cases.
at max-product BP, which generalises the Viterbi algorithm
While our preliminary study [45] applying an extension of [50], providing the MAP joint state of all variables in the
the fractional free energy (FFE) (introduced in [46]) to tracking tree-structured graph. At convergence of sum-product BP, the
problems showed promise, the absence of a rapidly converging, marginal distribution at a vertex v can be calculated as:
Y
BP-like algorithm for solving it has limited its practical use.
p(x
)
∝
ψ
(x
)
µi→v (xv ).
(2)
v
v
v
In this paper, we provide such an algorithm, and prove its
(v,i)∈E
convergence.
Subsequently, the FFE is used as a building block in a In the case of a Markov chain, if all vertices are jointly
multiple scan formulation, which is shown in section IV to Gaussian, BP is equivalent to a Kalman smoother. Similarly,
address the undesirable behaviour of the previous methods. The if all vertices are discrete, BP is equivalent to inference on an
proposed method results in a convex variational problem for the HMM using the forward-backward algorithm. BP unifies these
multiple scan model, and a convergent algorithm for solving algorithms and extends them from chains to trees.
the problem is developed. Importantly, association consistency
2 In the general setting, the joint distribution is a product of maximal
constraints from previous scans are retained. The sequential
cliques [14, p9]. Since the graph is undirected, we assume that E is symmetric,
version of the algorithm is motivated by a new variational i.e., if (i, j) ∈ E then (j, i) ∈ E. We need only incorporate one of these two
interpretation of JPDA.
factors in the distribution.

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

3

Inference in cyclic graphs (graphs that have cycles, i.e., that For tree-structured graphs, BP can be shown to converge to the
are not tree-structured) is far more challenging. Conceptually, optimal value of this convex, variational optimisation problem.
one can always convert an arbitrary cyclic graph to a tree by It has further been shown that the feasible set described by
merging vertices (e.g., so-called junction tree representations) (8)-(11) is exact, i.e., any feasible solution can be obtained by
[15], [47], but in practical problems, the dimensionality of a valid joint distribution, and any valid joint distribution maps
the agglomerated variables tends to be prohibitive. BP may to a feasible solution.
If a graph contains a leaf vertex xi that is connected only
be applied to cyclic graphs; practically, this simply involves
repeated application of (1) until convergence occurs (i.e., until to vertex xj via factor ψi,j (xi , xj ),3 then inference can be
the maximum change between subsequent messages is less than performed equivalently by eliminating vertex xi , and replacing
a pre-set threshold). Unfortunately, this is neither guaranteed the vertex factor for xj with [15, ch 9]
X
to converge to the right answer, nor to converge at all.
ψ̃j (xj ) =
ψi,j (xi , xj ).
(12)
Nevertheless, it has exhibited excellent empirical performance
xi
in many practical problems [36]. For example, the popular
iterative turbo decoding algorithm has been shown to be an Given the resulting marginal distribution for p(xj ) (or an
approximation thereof), the pairwise joint distribution (or belief)
instance of BP applied to a cyclic graph [23].
The current understanding of BP in cyclic graphs stems of (xi , xj ) can be reconstructed as:
from [17]. It has been shown (e.g., [14, Theorem 3.4]) that one
ψi,j (xi , xj )
.
(13)
p(xi , xj ) = p(xj )
can recover exact marginal probabilities from an optimisation
ψ̃j (xj )
of a convex function known as the Gibbs free energy. In the
single-vertex case (or if all variables are merged into a single Lemma 2 provides a variational viewpoint of vertex elimination,
interpreting it as a partial minimisation of the pairwise joint of
vertex), this can be written as described in lemma 1.
the neighbour and leaf, conditioned on the marginal distribution
Lemma 1. The Gibbs free energy variational problem for a
of the neighbour of the leaf. The theorem uses the conditional
single random variable x can be written as:
entropy, defined as: [51]
X
X
minimise − H(x) − E[log ψ(x)]
(3)
H(xi |xj ) = −
q(xj )
q(xi |xj ) log q(xi |xj )
(14)
q(x)
X
xj
xi
subject to q(x) ≥ 0,
q(x) = 1,
(4)
XX
q(xi , xj )
x
=−
q(xi , xj ) log P
. (15)
0
P
x0i q(xi , xj )
xi xj
where E[log ψ(x)] ,
x q(x) log ψ(x), and H(x) =
−E[log q(x)] is the entropy of the distribution q(x) (all Conditional entropy was shown to be a concave function of
expectations and entropies are under the distribution q). The the joint in [52]. Note that, due to marginalisation constraints,
solution of the optimisation is q(x) = P ψ(x)
∝ ψ(x).
0
x0 ψ(x )
H(xi |xj ) = H(xi , xj ) − H(xj ) = H(xi ) − I(xi ; xj ). (16)
Similar expressions apply for continuous random variables,
Lemma 2. Let J[q(xj )] be the solution of the following
replacing sums with integrals. The objective in (3) can be
optimisation problem:
recognised as the Kullback-Leibler (KL) divergence between
q(x) and the (unnormalised) distribution ψ(x). If the graph is
J[q(xj )] = minimise − H(xi |xj ) − E[log ψi,j (xi , xj )]
q(xi ,xj )≥0
a tree, the entropy can be decomposed as:
(17)
X
X
X
H(x) =
H(xv ) −
I(xi ; xj ),
(5)
subject to
q(xi , xj ) = q(xj ) ∀ xj .
(18)
v∈V

(i,j)∈E

I(xi ; xj ) = H(xi ) + H(xj ) − H(xi , xj ).

(6)

I(xi ; xj ) is the mutual information between xi and xj [51].
Accordingly, the variational problem can be written as:
X
minimise −
{H(xv ) + E[log ψv (xv )]}
q(xv ),q(xi ,xj )

−

X

v∈V

{−I(xi ; xj ) + E[log ψi,j (xi , xj )]}

(7)

subject to q(xi , xj ) ≥ 0 ∀ (i, j) ∈ E, ∀ xi , xj
X
q(xv ) = 1 ∀ v ∈ V

(8)

(i,j)∈E

(9)

xv

X

q(xi , xj ) = q(xi ) ∀ (i, j) ∈ E

(10)

q(xi , xj ) = q(xj ) ∀ (i, j) ∈ E.

(11)

xj

X
xi

xi

Then
J[q(xj )] = −

X

q(xj ) log

xj

X

ψi,j (xi , xj ),

(19)

xi

and the minimum of the optimisation of (17) is attained at
q(xi , xj ) = q(xj )q(xi |xj ), where
ψi,j (xi , xj )
q(xi |xj ) = P
.
0
x0 ψi,j (xi , xj )

(20)

i

Proof of this result can be found in [22, App C]. On treestructured graphs, BP may be viewed as successive applications
of variable elimination, followed by reconstruction.
If the graph has cycles (i.e., is not a tree), then the entropy
does not decompose into the form in (5)-(7). Furthermore, the
3 Without loss of generality, assume that neither vertex has a vertex factor,
as this can be incorporated into the edge factor.

4

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

feasible set in (8)-(11) is an outer bound to the true feasible set, state of all targets is X = (x1 , . . . , xn ). We denote the
i.e., there are feasible combinations of marginal distributions set of measurements received in scan s ∈ S = {1, . . . , S}
s
that do not correspond to a valid joint distribution. Nevertheless, by Zs = {z 1s , . . . , z m
s }. We consider a batch-processing
as an approximation, one may solve the optimisation in (7)- algorithm, where all scans s ∈ S are processed at once, and a
(11). The objective in (7) is referred to as the Bethe free energy sequential method, where scans are introduced incrementally,
(BFE) after [53], a connection identified in [17]. For a cyclic but some reprocessing is performed on each scan in the
graph, it differs from the Gibbs free energy, but is a commonly window after a new scan is revealed. Our goal is to avoid the
utilised approximation.
need to explicitly enumerate or reason over global association
It was shown in [17] that, if it converges, the solution hypotheses, i.e., hypotheses in the joint state space of all
obtained by BP is a local minimum of the BFE. The BP targets. Instead, marginal distributions of each target are stored,
message iterates in (1) can be viewed as a general iterative and dependencies between targets are accounted for using
method for solving a series of fixed point equations derived variational methods.
from the optimality conditions of the BFE variational problem
We assume that the prior information for each target is
(see [14, 4.1.3]). The marginal probability estimates obtained independent, such that the prior distribution of X is
using BP, denoted in this paper by the symbol q, are referred
n
Y
to as beliefs.
p(X) ∝
ψ i (xi ),
(21)
TRSP [14], [21] provides a convex alternative to the BFE, by
i=1
applying weights γi,j ∈ [0, 1] to the mutual information terms
i
i
I(xi ; xj ). If the weights correspond to a convex combination where the factors ψ (x ) collectively represent the joint.
We use the symbol i ∈ {1, . . . , n} to refer to a target index,
of embedded trees (i.e., a convex combination of weighted,
j
∈
{1, . . . , ms } to refer to a measurement index, and s ∈
tree-structured sub-graphs, where in a given graph γi,j = 1 if
S
to
refer to a measurement scan index. Each target i is
the edge is included, and γi,j = 0 otherwise), then the resulting
detected
in each scan s with probability Psd (xi ), target-related
free energy is convex. A rigorous method for minimising energy
measurements follow the model ps (z s |xi ), and false alarms
functions of this form was provided in [22].
fa
Finally, MF [18] approaches the problem by approximating occur according to a PPP with intensity λs (z s ).
The relationship between targets and measurements is
the entropy and expectation in (3) assuming that the joint
described
via a set of latent association variables, comprising:
distribution is in a tractable form, e.g., the product of the
1) For each target i ∈ {1, . . . , n}, an association variable
marginal distributions. Consequently, the expectations (the
ais ∈ {0, 1, . . . , ms }, the value of which is an index to
second term of (3)) are non-convex, and resulting methods
the measurement with which the target is hypothesised to
tend to underestimate the support of the true distribution, e.g.,
be associated in scan s (zero if the target is hypothesised
finding a single mode of a multi-modal distribution. It is also
to have not been detected)
possible to use a hybrid of MF and BP, as proposed in [19],
2) For each measurement j ∈ {1, . . . , ms }, an association
[20].
variable bjs ∈ {0, 1, . . . , n}, the value of which is an index
B. Multiple scan data association
to the target with which the measurement is hypothesised
to be associated (zero if the measurement is hypothesised
The problem we consider is that of data association across
to be a false alarm)
multiple scans, involving many targets, the state of which is
to be estimated through point measurements. We assume that This redundant representation implicitly ensures that each
many targets are present, each target gives rise to at most one measurement corresponds to at most one target, and each target
measurement (excluding so-called extended target problems, corresponds to at most one measurement. It was shown in [32]
where targets may produce multiple measurements), and each that, for the single scan case, this choice of formulation results
measurement is related to at most one target (excluding so- in an approximate algorithm with guaranteed convergence and
called merged measurement problems, e.g., where multiple remarkable accuracy.
s
targets fall within a resolution cell). We assume that false
Denoting as = (a1s , . . . , ans ) and bs = (b1s , . . . , bm
s ), the
alarms occur according to a Poisson point process (PPP).
joint distribution of the measurements and association variables
For clarity of presentation, we assume that all measurements for scan s can be written as:
related to the i-th target are independent conditioned on the


 Y

target state xi . If xi is the target state vector at a given time,
i
as i
d
i
p(Z
,
a
,
b
|X)
∝
P
(x
)p
(z
|x
)
s
s s
s s
s
this effectively restricts the problem such that the multiple
 i

i|as >0
scans correspond to different sensors at the same time instant,

 

or the target state is static. Problems involving multiple time
 Y
  Y

j
steps can be addressed by replacing xi with the joint state
×
[1 − Psd (xi )] ×
λfa
(z
)
s
s
 i
  j

over a time window, e.g., xi = (xi1 , . . . , xit ). An alternative
i|as =0
j|bs =0
approach involving association history hypotheses is discussed
× ψs (as , bs ), (22)
in [54, app A].
We assume that the number of targets n is known, though the where ψs (as , bs ) = 1 if as and bs form a consistent
method can easily be extended to an unknown, time-varying association event (i.e., if ais = j > 0 then bjs = i and vice
number of targets using the ideas in [10], [55]. The joint versa), and ψs (as , bs ) = 0 otherwise.

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

5

Fig. 2. Bipartite formulation of a single scan data association problem.

optimisation, but first we review the single scan formulation
and the underlying variational problem.
Following similar lines to [35], we show in [54, app B-A] that
the Bethe variational problem for the single scan formulation
can be solved by minimising the objective:

Fig. 1. Graphical model formulation of multiple scan problem.

The quantity of interest is the posterior distribution
Y
p(X, aS , bS |ZS ) ∝ p(X)
p(Zs , as , bs |X).

minimise
i,j
qs

(23)

s∈S
fa j
j=1 λs (z s )

Dividing (22) through by
(since the measurement
values are constants in (23)), this can be written as
p(X, aS , bS |ZS ) ∝



ms
n 

Y
Y
Y
ψsi (xi , ais )
ψ i (xi )
ψsi,j (ais , bjs ) , (24)


j=1

s∈S

qsi,j log

i=1 j=0
ms
n X
X

where

qsi,j
wsi,j

+

ms
X

qs0,j log qs0,j

j=1

(1 − qsi,j ) log(1 − qsi,j )

−

Qms

i=1

ms
n X
X

(27)

i=1 j=1

subject to

ms
X
j=0
n
X

qsi,j = 1 ∀ i ∈ {1, . . . , n}

(28)

qsi,j = 1 ∀ j ∈ {1, . . . , ms }

(29)

i=0

0 ≤ qsi,j ≤ 1

(30)

where qsi,j = q(ais = j) = q(bjs = i) is the belief that target i
is associated with measurement j (or, if j = 0, that target i
ψsi (xi , ais ) =
(25)
is missed, or if i = 0, that measurement j is a false alarm4 ),
1 − Psd (xi ),
ais = 0
and wsi,j = ψ i (ais = j) is the node factor, which will be
and the functions
defined subsequently in (32). The constraints in (28) and
(
i
j
j
i
(29)
are referred to as consistency constraints, as they are
0, as = j, bs 6= i or bs = i, as 6= j
(26) a necessary condition for the solution to correspond to a valid
ψsi,j (ais , bjs ) =
1, otherwise
joint association event distribution.
It is not obvious that the optimisation in (27)-(30) is convex,
provide a factored form of ψs (as , bs ), collectively ensuring
that the redundant sets of association variables (a1s , . . . , ans ) but it can be proven using the result in [35], which shows that
s
a closely related objective (excluding terms involving qsi,0 and
and (b1s , . . . , bm
s ) are consistent (i.e., setting the probability
of any event in which the collections are inconsistent to zero). qs0,j ) is convex on the subset in which (30) and either (28) or
A graphical model representation of (24) is illustrated in (29) apply. Details can be found in the proof of lemma 3.
In [46] it was shown that, if the correct fractional coefficient
figure 1.
Over multiple time steps, JPDA operates by calculating the γ ∈ [−1, 1] is incorporated on the final term in the objective
marginal distribution of each target, pi (xi ), fitting a Gaussian (27) (excluding false alarms and missed detections, and setting
to this distribution, and proceeding to the next time step ms = n), the value of the modified objective function at
approximating the joint prior distribution by the product of the optimum is the same as the Gibbs free energy objective.
these approximated marginal distributions. It may be applied In the formulation that incorporates false alarms and missed
similarly to multiple sensor problems, introducing an arbitrary detections, the fractional free energy (FFE) objective is:5
order to the sensors.
ms
ms
n X
X
X
q i,j
C. Single scan BP data association
FBγ ([qsi,j ]) =
qsi,j log si,j + γ
qs0,j log qs0,j
w
s
i=1 j=0
j=1
The single scan version of the model in section II-B was
ms
n X
X
studied in [30], [32], with similar formulations (excluding false
−γ
(1 − qsi,j ) log(1 − qsi,j ). (31)
alarms and missed detections) examined in [31], [35]; the
i=1 j=1
graphical model for the single scan data association problem
is illustrated in figure 2. In [32], [34], simplified BP equations
4 If j = 0 then there is no corresponding q(bj ), and if i = 0 then there
s
were provided, and convergence was proven; these results do
is no corresponding q(ais ).
not apply to the multiple scan problem. In the present work, we
5 Note that we reverse the sign of γ in comparison to [46], so that γ = 1
seek to address the multiple scan problem using ideas in convex yields the regular BFE.
( P d (xi )p
s

j
i
s (z s |x )
j
λfa
(z
)
s
s

,

ais = j > 0

6

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

Despite the fact that the “right” value of γ is not known for any
particular problem, our preliminary investigations in [45], and
the results in section IV, show that fractional values γ ∈ [0, 1]
can yield improved beliefs.6 Practically, the value could be
chosen a priori based on the problem parameters. It is straightforward to show that the inclusion of the fractional coefficient
γ ∈ [−1, 1] retains convexity of (31) (on the appropriate
subset); again, details are in the proof of lemma 3.
We now consider how these results may be applied to solve
the single scan problem incorporating the kinematic states xi ,
as illustrated in figure 3(a) (or, equivalently, the single scan
version of figure 1). Whilst a solution for the belief of xi
can be recovered from a solution of (27)-(30) using lemma
2, a variational problem formulation incorporating xi admits
extension to multiple scan problems. As in (24), this model
involves factors ψ i (xi ) and ψsi (xi , ais ), where in (27),
Z
wsi,j = ψsi (ais = j) = ψ i (xi )ψsi (xi , ais = j)dxi . (32)
It is shown in [54, app B-A] that we can arrive at the problem in
(27)-(30) by performing a partial minimisation of the following
objective over q i (xi ) and the pairwise joint qsi (xi , ais ):
FB ([q i (xi )], [qsi (xi , ais )], [qsi,j ]) = −

n
X


H(xi )

f (q) is closed and convex, then the conjugate dual of f ∗ (λ) is
the original function f (q), thus f (q) and f ∗ (λ) are alternate
representations of the same object.
Dual functions are useful in constrained convex optimisation
since the optimal value of the primal
min f (q)
q

max −f ∗ (AT λ),

i=1 j=1

subject to the constraints:

and if f (q) is strictly convex then, given the optimal solution
λ∗ of the dual, the optimal value q ∗ of the primal can be
recovered as the solution of the unconstrained optimisation
min f (q) − (AT λ∗ )T q.

qsi,j

Z
=
n
X

qs0,j

≥ 0,

= q i (xi ),

(34)
(35)

qsi (xi , j)dxi ,

(36)

qsi,j = 1,

(37)

One additional usefulness of conjugate duality over Lagrangian duality is its ability to tractably address objectives
that decompose additively. In this work, we utilise the primaldual framework developed in [22], which addresses problems
of the form:
n
X
min f (q) +
hi (q),
(43)
i=1

where f (q) and hi (q) are proper, closed, convex functions.
Constraints are addressed by admitting extended real-valued
functions. It is shown in [22], [57] that the dual of (43) is

λ1 ,...,λn

Pn

i=1 λi ) −

n
X

h∗i (λi ).

(44)

i=1

Thus, assuming smoothness of f ∗ (or strict convexity of f ),
the dual optimisation can be performed via block coordinate
ascent, iteratively performing the following steps for each i:
X
µ :=
λj ,
(45)
j6=i

i=0
ms Z
X

(42)

q

max −f ∗ (−

qsi (xi , ais ) ≥ 0,
ms
X
qsi (xi , ais )
i
as =0

(41)

λ

q

+ E[log ψ i (xi )] + H(ais |xi ) + E[log ψsi (xi , ais )]
m
ms
n X
s
X
X
0,j
0,j
+
qs log qs −
(1 − qsi,j ) log(1 − qsi,j ), (33)

(40)

is the same as the optimal value of the dual optimisation

i=1

j=1

subject to Aq = 0

λi := arg max −f ∗ (−λi − µ) − h∗i (λi ).

(46)

λi

qsi (xi , ais )dxi = 1.

(38)

ais =0

D. Conjugate duality and primal-dual coordinate ascent
The Fenchel-Legendre conjugate dual of a function f (q) is
defined as: [56]
∗

T

f (λ) = sup λ q − f (q).

(39)

q

The dual f ∗ (λ) is convex regardless of convexity of f (q),
since it is constructed as the supremum of a family of linear
functions. The key outcome of conjugate duality is that, if
6 e.g., in high SNR cases (very high P d , very low λfa ), the BFE objective
tends to yield solutions that are almost integral, i.e., are closer to MAP
solutions rather than marginal probabilities, as illustrated later in figure 6(a).
The inclusion of the γ coefficient on the term involving qs0,j retains the
property of the BFE that the optimisation provides a near-exact result when
targets are well-spaced.

The method in [22] shows that the block optimisations
required in (46) can be performed via primal minimisations, i.e.,
the updated value λi can be obtained through the optimisation
q ∗ := arg min f (q) + hi (q) + µT q,

(47)

q

λi := −µ − ∇f (q ∗ ).

(48)

In [22], the authors develop the norm-product belief propagation algorithm for convexifications of general PGM inference
problems. While these methods could be applied to the problem
of interest, they do not exploit the unique problem structure and
resulting convexity discussed in section II-C. Consequently,
the necessary convexification procedure would produce an
unnecessarily large change to the BFE. Thus we adopt the
optimisation framework of [22], but the solution does not
exactly fit the NPBP algorithm, and so it is necessary to develop
it from the basic framework.

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

III. VARIATIONAL MULTIPLE SCAN DATA ASSOCIATION

7

It can be shown that the modification of the factors ψ̄1i (xi , ai1 )
effectively incorporates the Lagrange multipliers for the constraints that are being relaxed, and a linearisation of the concave
term (the third line in the objective in figure 3(a)), such that
the change in 3) above attempts to nullify changes 1) and
2), i.e., the solution of the modified variational problem is
the same as the original one. The linearisation of the concave
term is reminiscent of a single iteration of the convex-concave
procedure in [58].
The proposed solution is illustrated in figure 3(d). Compared
to JPDA-BP (figure 3(b)), it makes the following modifications:

The standard approach to tracking using JPDA and related
methods is to calculate the marginal distribution of each target,
and proceed to the next scan approximating the posterior as
the product of the single-target marginal distributions.7 In
many cases, this approach is surprisingly effective. The method
proposed in this work is based on a variational interpretation of
the JPDA approach, which gives rise to a family of formulations
that includes the JPDA-like approach and MSBP. We refer to the
JPDA-like approach using the BP approximation of marginal
association probabilities (beliefs), as JPDA-BP. In comparison,
true JPDA uses exact marginal association probabilities, and
1) As in the FFE objective in (31), weP
applyPfractional
n
ms
approximates the posterior as a Gaussian at each step.
weights γs to the concave terms − i=1 j=1
(1 −
P
i,j
i,j
In particular, we examine the single scan BP approach (e.g.,
qs ) log(1 − qs ), such that if
s γs < 1 then the
[10]) which calculates association beliefs, and approximates
objective function willPbe strictly convex
n
the joint as the product of the beliefs (although the steps taken
2) We retain constraints i=0 qsi,j = 1 for both scans
are not unique to the BP estimate of marginal probabilities).
3) If the coefficient γs differs from the value used for that
Consider the example in figure 3(a), where in scan 1 (S = {1})
scan at the previous time step, we modify the factor
the beliefs are calculated through the optimisation of (33).
ψsi (xi , ais ) such that the solution remains unchanged
Denote the solution as q i (xi , ai1 ) and q1i,j . JPDA-BP moves
(before the following scan is incorporated)
forward to the next scan approximating the joint as:
The reason for the latter step is that it was found to be desirable
n
Y
to use values γs closer to 1 in the current (most recent) scan,
p(X, a1 ) ≈
q i (xi , ai1 ).
(49) and reducing to zero in earlier scans. In our experiments, we use
i=1
the selection γs = 0.55 in the current scan, and γs = 0 in past
This approximated joint distribution can be formulated as a scans. In this case, the algorithm approximates the objective
PGM using the graph in the left-hand side of figure 3(b), function in a similar manner to JPDA-BP, but the consistency
where the factor ψ̄1i (xi , ai1 ) is modified such that the simplified constraints from previous scans are retained. Thus, unlike
formulation results in the same solution as the original problem JPDA, when later scans are processed, the constraints which
in figure 3(a). This does not mean to say that the objective on ensure that the origin of past measurements is consistently
the left-hand side of figure 3(b) is equivalent to that in figure explained remain enforced. In section IV, we will see that this
3(a). Note that we could eliminate the nodes ai1 from the graph can result in improved performance.
as they will subsequently remain leaf nodes, but we choose to
JPDA and JPDA-BP operate sequentially, at each stage
retain them for comparison to the proposed algorithm.
introducing a new set of nodes, as in figure 3(b). Following
When a second scan of measurements is introduced (S = approximation, association variables from past scans are leaf
{1, 2}), JPDA-BP effectively solves the problem in the right- nodes, and thus can be eliminated. Like MSBP, the proposed
hand side of figure 3(b) where, for the newest scan (s = method maintains past association variables, and introduces a
2), ψ̄2i (xi , ai2 ) , ψ2i (xi , ai2 ). Although the data for q1i,j (i.e., new set in each scan. The sequential variant of the algorithm
ψ̄1i (xi , ai1 )) is unchanged from the left-hand side of figure 3(b), operates by performing one or more forward-backward sweeps
introduction of new information in scan 2 will, in general, over recent scans. Complexity could be mitigated by performing
modify the belief values for the first scan, q1i,j . Consistency the operations on past scans intermittently, rather than upon
constraints (specifically, (37)) will no longer hold. If BP was receipt of every scan.
applied directly without the approximations made by JPDA (i.e.,
In what follows (and in figure 3), we assume that the target
(49)), we would arrive at the problem in figure 3(c). Comparing state xi is discrete. This can be achieved by using a particle
the variational problems at scan 2 in figures 3(b) (JPDA-BP) representation, such that the prior distribution ψ i (xi ) contains
and 3(c) (MSBP), we note the following differences:
the
weights of the particles. An expression such as
P prior
Pn Pms
i
i
i,j
i,j
ψ
(x
)ψsi (xi , ais ) should be interpreted as the sum of
1) The concave term − i=1 j=1 (1 − qs ) log(1 − qs )
xi
appears only for scan 2 in figure 3(b), but for both scans the function evaluated at the particle locations. The output
of the inference procedure developed is used to reweight the
in figure 3(c) P
n
2) The constraint i=0 qsi,j = 1 appears only for scan 2 in particles, e.g., recovering an estimate of the PDF of continuous
kinematic state ξ i as:
figure 3(b), but for both scans in figure 3(c)
X
3) The factors ψ̄1i (xi , ai1 ) have been modified in figure 3(b)
q i (ξ i ) ≈
q i (xi )δ(ξ i − xi ).
(50)
but remain at their original values in figure 3(c)
4) The objective function in figure 3(b) is convex, whereas
the objective in figure 3(c) is non-convex
7 It can be shown (e.g., [15, p277]) that the product of the marginals is
the distribution with independent targets that best matches the exact joint
distribution.

xi

The association history formulation of appendix A may be
applied simply by replacing xi with the single target association
history aiS . Some simplifications can be made in this case since
H(ais |aiS ) = 0.

8

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

Fig. 3. Probabilistic graphical model and Bethe free energy (a) single scan formulation, (b) multiple scan formulation using JPDA-BP approximation, (c)
standard multiple scan approach using MSBP, and (d) convex multiple scan approach, where shading depicts the weighting of the corresponding entropy and
mutual information terms. The formulation in the diagrams excludes false alarm events (qs0,j ) for simplicity; these events are modelled in the formulation in
the text. For the newest scan (s = 2) in the right of (b) and (d), ψ̄2i (xi , ai2 ) , ψ2i (xi , ai2 ).

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

A. Bethe free energy function
In [54, app B-B], we show that the Bethe free energy for
the multiple scan formulation in section II-B (and figure 3(c))
can be written as:
FB ([q i (xi )], [qsi (xi , ais )], [qsi,j ]) =
n
X

H(xi ) + E[log ψ i (xi )]
−
i=1

−

n
XX


H(ais |xi ) + E[log ψsi (xi , ais )]

s∈S i=1

+

ms
XX

qs0,j log qs0,j −

s∈S j=1

ms
n X
XX

(1 − qsi,j ) log(1 − qsi,j ),

s∈S i=1 j=1

(51)

9

The use of re-weighting to obtain a convex free energy is
closely related to the TRSP algorithm [21].
In the development that follows, we provide a decomposition
of the convex free energy that permits application of primaldual coordinate ascent (PDCA). First we give the basic form,
and show that the components are convex. Then we state
weights which ensure that the objective is the same as (58).
Then, in section III-C, we provide algorithms to minimise each
block that needs to be solved in PDCA.
In order to optimise (58), we consider decompositions of
the expression of the form:
FBγ,β ([q i (xi )], [qsi (xi , ais )]) = f ([q i (xi )], [qsi (xi , ais )])
Xn
+
hs,1 ([q i (xi )], [qsi (xi , ais )])
s∈S

subject to the constraints:

o
+ hs,2 ([q i (xi )], [qsi (xi , ais )]) , (60)

qsi (xi , ais ) ≥ 0, q i (xi ) ≥ 0
ms
XX
qsi (xi , ais ) = 1

(52)
where
(53)

xi ais =0
ms
X

qsi (xi , ais ) = q i (xi ) ∀ i, xi ∀ s ∈ S

(54)

f ([q i (xi )], [qsi (xi , ais )])
n
X

=−
κf,x H(xi ) + E[log ψ i (xi )]
i=1

ais =0

qsi,j

=

X

qsi (xi , j)

xi
i,j
qs ≥
n
X
qsi,j
i=0

∀ i, j, ∀ s ∈ S

(55)
(56)

= 1 ∀ j, ∀ s ∈ S,

(57)

B. Convexification of energy function
In this work, we consider convex free energies of the form:
FBγ,β ([q i (xi )], [qsi (xi , ais )], [qsi,j ])
n
X

i

=

i=1

−

H(ais |xi ) + E[log ψsi (xi , ais )]

s∈S i=1

+

X
s∈S

−

X
s∈S

γs

βs

ms
X

qs0,j log qs0,j

j=1

ms
n X
X

(1 −

qsi,j ) log(1

−

qsi,j ).

(58)

i=1 j=1

The difference between (58) and (51) is the incorporation of
the coefficients γs ∈ [0, 1) and βs ∈ (0, 1] in the final two
terms. We will show that (58) is convex if η ≥ 0 (strictly
convex if η > 0), where
X
η =1−
γs .
(59)
s∈S

hs,1 ([q i (xi )], [qsi (xi , ais )])
n
X

=−
κs,1,x H(xi ) + κs,1,s H(xi , ais )
i=1
ms
X

+ βs

qs0,j log qs0,j

j=1

− γs

ms
n X
X

(1 − qsi,j ) log(1 − qsi,j ),

i=1 j=1
i
i
hs,2 ([q (x )], [qsi (xi , ais )])
n
X

κs,2,x H(xi )
=−
i=1

+ κs,2,s H(xi , ais ) .

(62)

(63)

Note that we do not consider (60) or (62) to depend on [qsi,j ]
as these are uniquely determined from [qsi (xi , ais )] by the
constraints in (55)–(57) (which will be enforced whenever
we consider the block hs,1 , the only block that includes these
variables).
Immediate statements that can be made regarding convexity
of (61)–(63) include:

H(x ) + E[log ψ i (xi )]

n
XX


κf,s H(xi , ais ) + E[log ψsi (xi , ais )] , (61)

s∈S i=1

0 ∀ i, j, ∀ s ∈ S

where q i (xi ) is the belief (i.e., approximate probability) of the
state of target i, qsi,j is the belief that target i is associated
with measurement z js , qs0,j is the belief that measurement z jS
is not associated with any target, and the set S indexes the
measurement scans under consideration.

−

−

n
XX


1) f is strictly convex if κf,x > 0 and κf,s > 0; this is the
consequence of the convexity of entropy
2) hs,2 is convex if the consistency constraints (54) are
enforced for scan s, κs,2,s ≥ 0, and κs,2,x + κs,2,s ≥ 0;
this is the consequence of the convexity of entropy, and
of conditional entropy
Convexity of hs,1 is proven in the lemma below. Subsequently,
the decomposition in terms of f , hs,1 and hs,2 (for each s ∈ S)
is utilised in the PDCA framework introduced in section II-D.

10

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

Lemma 3. If constraints (52)–(57) are enforced for scan s,
κs,1,s ≥ 0, βs ≥ 0 and κs,1,s + κs,1,x ≥ γs , then hs,1 is
convex.
The proof of lemma 3 is in [54, app C]. The following
lemma provides coefficients which fit (58) into the form (60),
in order to permit solution using PDCA.
Lemma 4. Given the following, (60) is equivalent to (58):
η
κf,x = κf,s =
,
(64)
S+1
η
κs,1,x = −
,
(65)
S+1
η
,
(66)
κs,1,s = γs +
S
+
1


2η
κs,2,x = − 1 − γs −
,
(67)
S+1
2η
κs,2,s = 1 − γs −
,
(68)
S+1
where S = |S| is the number of scans in the problem
Lemma 4 can be shown by substituting the coefficients into
(60), and showing that the coefficients of H(xi ) sum to −(S −
1), and the coefficients of H(xi , ais ) sum to 1. Examining the
values in lemma 4, we find (60) is convex (under the constraints
(52)–(57)). This in turn shows that (58) is convex.

where
φi (xi ) = log ψ i (xi ) − µi (xi ),

= log ψsi (xi , ais ) − µis (xi , ais ),
(72)
#
"


X
φi (xi ) + φis (xi , ais )
. (73)
φ̃is (ais ) = log
exp
κf,s + κs,1,s
i
x

The single scan marginal qsi (ais ) is the solution of the following
sub-problem:
minimise

− γ̃s

minimise f ([q i (xi )], [qsi (xi , ais )])
)], [qsi (xi , ais )])

+ hs,1 ([q (x
n X
X
+
µi (xi )q i (xi )
i=1 xi

+

τ ∈S i=1

xi

qsi,j
wsi,j

+ β̃s

ms
X

qs0,j log qs0,j

j=1

(1 − qsi,j ) log(1 − qsi,j ),

(74)

i=1 j=1

Pms i,j
qs =
subject to (56)–(57) and the additional constraint j=0
βs
1 ∀ i, derived from (53) and (55). In (74), β̃s = κf,s +κ
,
s,1,s
γs
i,j
i
γ̃s = κf,s +κs,1,s , and log ws = φ̃s (j). This sub-problem is
studied in theorem 1. The update to λ is
c

λis,1,x (xi ) = φi (xi ) − κf,x log q i (xi ),
c
λis,1,s (xi , ais ) =

q i (xi ) =

φis (xi , ais ) − κf,s
ms
X
qsi (xi , ais ),
ais =0

log qsi (xi , ais ),

(75)
(76)
(77)

where = denotes equality up to an additive constant. For other
scans τ ∈ S, τ 6= s, λis,1,τ (xi , aiτ ) = 0.
Lemma 6. The solution of block (s, 2):
minimise f ([q i (xi )], [qsi (xi , ais )])
+ hs,2 ([q i (xi )], [qsi (xi , ais )])
n X
X
+
µi (xi )q i (xi )
i=1 xi

+

Lemma 5. Consider the problem to be solved in block (s, 1):

mτ
n X X
XX

qsi,j log

c

The convex free energy in (60) is of the form (43), so we
propose a solution using PDCA. As discussed in section II-D,
this is achieved by iterating (45), (47) and (48), commencing
with λs = 0 ∀ s ∈ S. To begin, we note that if the gradient of
a block hi with respect to a subset of variables is zero,8 then
those updated dual variables in (48) will be zero.
The algorithm proceeds by repeatedly cycling through blocks
hs,1 and hs,2 for each s ∈ S. The following lemmas provide
algorithms for solving each block in turn; the proofs can be
found in [54, app C].

i

ms
n X
X

i=1 j=0
ms
n X
X

C. Solution of convex energy

i

(71)

φis (xi , ais )

µiτ (xi , aiτ )qτi (xi , aiτ )

(69)

aiτ =0

under the constraints (52)–(57), where (54)–(57) is applied
only for scan s, and κf,x + κs,1,x = 0. The solution of this
problem is given by:
n i i i i i o
φ (x )+φs (x ,as )
exp
κf,s +κs,1,s
i
i i
i i
qs (x , as ) = qs (as ) ×
,
(70)
exp φ̃is (ais )
8 As discussed in [22], constraints for each block can be incorporated into
hi , so this condition also implies that no constraints are incorporated for the
subset of variables in block hi .

mτ
n X X
XX
τ ∈S i=1

xi

µiτ (xi , aiτ )qτi (xi , aiτ ),

(78)

aiτ =0

under the constraints (52)–(54) (including (54) only for scan
s) is:
(
)
φi (xi ) + φ̃is (xi )
i
i
q (x ) ∝ exp
, (79)
κf,x + κf,s + κs,2,x + κs,2,s
n i i i o
φs (x ,as )
exp κf,s
+κs,2,s
i
i i
i
i
n
o,
qs (x , as ) = q (x ) ×
(80)
φ̃is (xi )
exp κf,s +κs,2,s
where
φi (xi ) = log ψ i (xi ) − µi (xi ),
φis (xi , ais )
φ̃is (xi )

= log ψsi (xi , ais ) − µis (xi , ais ),
= (κf,s + κs,2,s )


 i i i 
ms
X
φ
(x
,
a
)
s
s
.
exp
× log 
κ
+
κ
f,s
s,2,s
i
as =0

(81)
(82)

(83)

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

1
2
3
4
5
6

7
8
9

10

input : n, m ∈ N; wi,j ∀ i ∈ {0, . . . , n}, j ∈ {0, . . . , m};
γ ∈ [0, 1); α ∈ (0.5, ∞) ∩ [γ, ∞);
β ∈ (0.5, ∞) ∩ [γ, ∞)
output : q i,j ∀ i ∈ {0, . . . , n}, j ∈ {0, . . . , m}
κ := −1 − γ + α + β
Initialise
y i,j := 1 ∀ i ∈ {0, . . . , n}, j ∈ {0, . . . , m}
Iterate to convergence
repeat
−(1−γ)

P
′
′
×
xi,j := w0,j y 0,j + i′ wi ,j y i ,j
−γ

P
′
′
× eκ ∀ i, j
w0,j y 0,j + i′ 6=i wi ,j y i ,j
1

xi,0 := (y i,0 ) α −1 ∀ i
−1

P
′
′
∀j
x0,j := w0,j y 0,j + i′ wi ,j y i ,j


P
′
′ −(1−γ)
y i,j := wi,0 xi,0 + j ′ wi,j xi,j
×


P
′
′ −γ
× eκ ∀ i, j
wi,0 xi,0 + j ′ 6=j wi,j xi,j


P
′
′ −1
∀i
y i,0 := wi,0 xi,0 + j ′ wi,j xi,j

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

1

11
12
13
14
15

y 0,j := (x0,j ) β −1 ∀ j
until sufficiently small change in y i,j , y i,0 , y 0,j
Calculate outputs (q 0,0 , 0)
q i,j := wi,j y i,j x0,j ∀ i, j
1
1
q i,0 := wi,0 (y i,0 ) α ∀ i, q 0,j := w0,j (y 0,j ) 1−β ∀ j

Fig. 4. Fractional BP algorithm for optimising single scan fractional free energy
q i,0
i,0
(74), where α is the coefficient
P of q log wi,0 (i.e., α = 1 in this instance).
Calculations marked ∀i or i are over
P the range i ∈ {1, . . . , n} (excluding
i = 0), while those marked ∀j or j are over the range j ∈ {1, . . . , m}
(excluding j = 0).

18

11

input : S; ψ i (xi ), ψsi (xi , ais ); γs , βs , ∆γs , ∆βs
output : q i (xi ), qsi (xi , ais ), ψ̄si (xi , ais )
Initialise
Set κ values according to lemma 4
λis,1,x (xi ) := 0, λis,1,s (xi , ais ) := 0
λis,2,x (xi ) := 0, λis,2,s (xi , ais ) := 0
Perform primal-dual coordinate ascent until convergence
repeat
for s ∈ S do
Solve blockP(s, 1)
P
µs (xi ) := τ ∈S\{s} λis,1,x (xi ) + τ ∈S λis,2,x (xi )
µs (xi , ais ) := λis,2,s (xi , ais )
Calculate λis,1,x (xi ), λis,1,s (xi , ais ) using lemma 5
Solve blockP(s, 2)
P
µs (xi ) := τ ∈S λis,1,x (xi ) + τ ∈S\{s} λis,2,x (xi )
µs (xi , ais ) := λis,1,s (xi , ais )
Calculate λis,2,x (xi ), λis,2,s (xi , ais ), q i (xi ) and
qsi (xi , ais ) using lemma 6
end
until sufficiently small change in λis,1,x (xi ), λis,2,x (xi ),
λis,1,s (xi , ais ), λis,2,s (xi , ais )
Calculate modified ψ̄si (xi , ais ) using theorem 3 (if
required)

Fig. 5. PDCA algorithm for minimising convex free energy based on
decomposition described in lemma 4.

D. Modification of factors between scans

As discussed at the beginning of this section (and illustrated
in figure 3), we propose solving a problem at scan S involving
a recent history of scans of measurements s ∈ S with fractional
weights configured to give high accuracy in the newest scan, and
using lower values in earlier scans. The main goal of retaining
historical scans is to ensure that consistency constraints from
The update to λ is
past scans remain enforced.
Suppose that we solve the multiple scan problem at scan S.
c
λis,2,x (xi ) = φi (xi ) − κf,x log q i (xi ),
(84) When we move to scan S 0 = S + 1, we will set γS = 0, this
c
λis,2,s (xi , ais ) = φis (xi , ais ) − κf,s log qsi (xi , ais ).
(85) time using the larger value for γS 0 . Thus we seek to modify
the problem parameters at time S to counteract the change of
reducing γs to zero. This is analogous to the approximation
For other scans τ ∈ S, τ 6= s, λis,2,τ (xi , aiτ ) = 0.
that JPDA makes, approximating the posterior as the product
of the marginals.
Theorem 1. The iterative procedure in figure 4 converges to
More generally, suppose we have been using weights γs ,
the minimum of the problem in (74), provided that β̃s > 0.5, s ∈ S, and at the next time, we will change these to γ̄s .
γ̃s ∈ [0, 1), and a feasible interior solution exists.
Similarly, suppose that the coefficients of the terms involving
qs0,j were βs , and will be changed to β̄s . The following theorem
This theorem is proven in [54, app D]. Implementation of gives
the modification to the problem parameters necessary to
the algorithm can be challenging due to numerical underflow. ensure that the solution of the problem remains unchanged.
This can be mitigated by implementing the updates in the log
domain using well-known numerical optimisations for log-sum- Theorem 3. Let [q i (xi )], [qsi (xi , ais )] and [qsi,j ] be the solution
exp [59, p844]. An alternative method for solving this form of of the problem in (58) using fractional weights γs and βs , s ∈ S.
problem based on Newton’s method was provided in [45]; the Suppose that the weights are changed to γ̄s = γs + ∆γs and
iterative BP-like method in figure 4 is significantly faster in β̄s = βs + ∆βs , and the problem parameters are changed as
follows:
most cases.
Theorem 2. The iterative procedure in figure 5 converges to
the minimum of the overall convex free energy,P
provided that
weights are as given in lemma 4, γs ≥ 0, and s∈S γs < 1.
This theorem is a corollary of claim 8 in [22], recognising
that the algorithms are an instance of this framework.

log ψ̄si (xi , ais = j) = log ψsi (xi , j)
+ ∆γs [1 + log(1 − qsi,j )] − ∆βs [1 + log qs0,j ] (86)
for j > 0, and log ψ̄si (xi , 0) = log ψsi (xi , 0) remains unchanged. Then the solution of the modified problem, denoted
[q̄ i (xi )], [q̄si (xi , ais )] and [q̄si,j ], is unchanged.

12

The proof of the theorem is in [54, app E]. Figure 5 includes
a step to incorporate these modifications.
E. Uses and limitations of proposed method

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

alarms are uniform over the bearing range covering the region
of interest.
We compare to two variants of JPDA, both of which maintain
a particle representation of each target location, and utilise BP
to approximate data association probabilities. The first variant
(which we refer to as JPDA-PBP, with ‘P’ denoting parallel)
processes measurements from different sensors in parallel,
solving each single-sensor problem once, as in [43] (this is
better suited to maintaining track rather than initial localisation).
The second variant (JPDA-BP) approaches multi-sensor data
by sequentially processing individual sensors, similar to the
IC-TOMB/P approach in [43].
For the proposed method, we compare:
1
1) γs = βs = |S|+1
for each scan, weights according to
lemma 4, and solving using figure 5, not utilising sequential modification (i.e., the final line of the algorithm); we
refer to this as the convex variational (CV) algorithm
2) The method using the sequential modification of section
III-D, introducing a new scan at each step with γs = 0.55,
for past scans setting γs = 0, and with βs = κf,s +κs,1,s ,
again solving using figure 5; we refer to this as the CVsequential (CVS) algorithm

The proposed method seeks to estimate marginal distributions
of target states. This provides a complete summary of the information available when considering each target separately, and
is useful, for example, when seeking to provide a confidence
region for the target location, or when deciding whether it
is necessary to execute sensor actions which will provide
clarifying information. The experiments in the following
section demonstrate that the proposed methods provide a
scalable approach for solving problems of this type, addressing
limitations experienced using existing methods.
In some tracking problems, a particular type of uncertainty
arises, in which multiple modes appear in the joint distribution
which essentially correspond to exchanges of target identity; the
coalescence problem in JPDA is a well-known example of this
(e.g., [60]). In such instances, multiple modes will appear in the
estimates of the marginal distributions; this is by design, and
is a correct summary of the uncertainty which exists. In these
cases, extracting point estimates from marginal distributions
is not straight-forward, but can be performed using methods
such as the variational minimum mean optimal subpattern A. Illustrative example
assignment (VMMOSPA) estimator [61]. Alternatively, if a
The result in figure 6 illustrates the behaviour on a simplified
point estimate is all that is required, MAP-based methods can
version
of the problem, with three sensors, two targets, and
be used. Likewise, MF or hybrid MF-BP methods (e.g., [20],
a
low
false
alarm rate (10−6 ). The top row (a)-(e) shows the
applied to tracking in [38], [39]) tend to provide good estimates
of a particular mode (and hence point estimates), at the expense results for different algorithms utilising the first two sensors,
of not characterising the full multi-modal uncertainty which where the first sensor initialises the distribution for each
target (drawing particles along each bearing line), and the
exists (see [39]).
Past experiments (e.g., [32]) have shown that the accuracy second permits triangulation. The sensor locations are shown
of the beliefs provided by BP is highest when SNR is low, as triangles, while true target locations are shows as crosses.
e.g., high false alarm rate and/or low probability of detection. Measurements are illustrated as dotted grey lines. The marginal
Conversely, accuracy is lowest in high SNR conditions, e.g., low distributions of the two targets, as estimated by the various
false alarm rate, high probability of detection. As demonstrated algorithms studied, are shown in the background image.
• Due to the low false alarm rate, JPDA-BP (shown in (a))
in the next section, this can now be mitigated through the
essentially provides a MAP association. The marginal
use of FFE. This behaviour is complementary to traditional
distribution estimates indicate high confidence for each
solution techniques, which perform very well in high SNR
target, in an incorrect (ghost) location, assigning near-zero
conditions (when ambiguity is the least, permitting tractable,
probability density to the true target location.
exact solution) but fail in low SNR conditions where many
• JPDA-FBP(0.55) (shown in (b)) utilises the fractional
targets are interdependent.
BP method, based on figure 4 (for which convergence
IV. E XPERIMENTS
is proven as theorem 1), with γ = 0.55. The figure
The proposed method is demonstrated through a simulation
demonstrates correct characterisation of the uncertainty
which seeks to estimate the marginal distribution of several
in the problem, with each marginal distribution estimate
targets using bearings only measurements. The region of interest
showing significant probability density in both the true
is the square [−100, 100]2 ⊂ R2 . Tracks are initialised using a
target location, and the ghost location.
single accurate bearing measurement from one sensor, corrupted
• Since tracks are initialised using sensor 1, and updated
by Gaussian noise with 0.1◦ standard deviation (e.g., as may
using sensor 2, the two sensor problem is a single scan
be provided if there is an accurate, low false alarm rate
association problem, and JPDA-PBP and MSBP (shown in
sensor providing bearing measurements from a single location);
(c) and (d)) is identical to JPDA-BP, similarly indicating
particle filter representations of each track are initialised by
high confidence for each target in a ghost location, and
randomly sampling from the posterior calculated by combining
assigning near-zero probability density to the true target
these measurements with a uniform prior on the region of
location.
interest. The proposed algorithm is then utilised to refine the
• Similarly, CV (shown in (e)) is identical to JPDA-FBP,
sensor positions. Target-originated bearing measurements are
and again correctly characterises the uncertainty in the
corrupted with 1◦ standard deviation Gaussian noise. False
problem.

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

(Two sensors)

(a) JPDA-BP

(c) JPDA-PBP

(e) CV

(d) MSBP

100

100

100

100

100

0

0

0

0

0

1
−100
−100

2
0

100

0

0

100

3

2
0

−100
−100

100

−100
−100

−100
−100
100

0

(h) JPDA-PBP
3

2
0

100

2
0

100

−100
−100

100

0

−100
−100

100

(j) CV
3

0
2

1

2
0

−100
−100

100

3

0
1

2

1

(i) MSBP
100

0
1

1

2

100

0
1

−100
−100

−100
−100

1

2

100

3

100

1

(g) JPDA-FBP(0.55)

(f ) JPDA-BP
(Three sensors)

(b) JPDA-FBP(0.55)

13

0

100

1
−100
−100

2
0

100

Fig. 6. Example problem involving two targets and two or three sensors. Targets are marked as ‘+’, and sensors as ‘4’, and bearing measurements are
illustrated as dotted grey lines. Tracks are initialised with measurements from sensor 1, and updated with two measurements from sensor 2 (top row, (a)-(e)).
Bottom row (f)-(j), incorporates an additional scan from sensor 3, in which a single measurement was received. Background image shows marginal distribution
estimates for the two tracks overlaid.

The bottom row (f)-(j) shows the results for the same algorithms introducing a third sensor, which receives a measurement
on one of the two targets.
•

•

•

•

•

Utilising the measurement from the third sensor, JPDA-BP
(shown in (f)) correctly localises one of the two targets,
but the second remains invalid, indicating high confidence
in a ghost location, and assigning near-zero probability
to the true target location.
JPDA-FBP is shown in (g) to correctly localise one of
the two targets, but a bimodal distribution remains for the
second. This is unnecessary: localisation of the first target
effectively confirms that the upper measurement from
sensor 2 belongs to that target, which in turn confirms that
the lower measurement belongs to the other target. Thus
the distribution exhibits unnecessarily high uncertainty as
a consequence of not enforcing past association feasibility
constraints.
Because data from sensor 2 is not used when interpreting
data from sensor 3, JPDA-PBP (shown in (h)), still
indicates high confidence in a ghost location for both
targets, and near-zero probability in the true location.
By simultaneously optimising over multiple scans of data,
MSBP (shown in (i)) is able to recover from the incorrect
solution in (d) and arrive at the correct solution.
Likewise, by retaining association feasibility constraints
from previous scans, CV (shown in (j)) is able to utilise
the confirmation of the location of one target to resolve the
bimodal uncertainty in the other target, correctly localising
both targets.

Of the methods shown, JPDA-BP, JPDA-PBP and MSBP
exhibit false confidence in ghost locations in figures (a), (c),
(d), (f) and (h), and JPDA-FBP fails to localise the distribution
to the extent possible in (g). CV is the only method shown
that is able to correctly characterise the uncertainty present in
both cases.

B. Quantitative analysis
As illustrated through the example in figure 6, the goal in
this work is to produce a faithful estimate of the marginal
probability distributions. This is quite different to problems in
which the aim is to produce a point estimate of target location,
for which the standard performance measure is mean square
error (MSE). In 200 Monte Carlo trials of the scenario in figure
6(a)-(e), JPDA-BP, JPDA-PBP and MSBP resolve uncertainty
to a single mode for each target, which is correct 76% of the
time, and incorrect (as in figure 6(a), (c) and (d)) 24% of the
time. This incorrect resolution of uncertainty could lead to dire
outcomes if the decision is made to take an action based on the
incorrect characterisation of uncertainty (which indicates high
confidence in a single, incorrect mode, as in figure 6(a), (c) and
(d)), rather than wait until further information is obtained (as
the characterisation in figure 6(b) and (e) would direct). Since
MSE is a measure only of the proximity of the single point
estimate to the true location, it does not capture the correctness
of the uncertainty characterised in the marginal probability
distribution, and it is not an adequate measure in this class of
problem.
Instead, we utilise two performance measures, which directly
measure occurrences of the undesirable outcomes in figure
6(a), (c), (d), (f), (g) and (h), and which are known to behave
consistently in the presence of multi-modal uncertainty. The
performance measures are the entropy of the beliefs produced
by each method, and the high probability density (HPD) value
in which the true location lies. The HPD value is defined as
the total probability under a distribution that is more likely
than a given point; e.g., given a point x∗ (the true location of
the target) and a distribution p, the HPD value is
Z
∗
HPD(p, x ) =
p(x)dx.
(87)
x|p(x)≥p(x∗ )

Three examples of this are illustrated in figure 7; in each case,
the true target location is marked as x∗ , and the HPD value is

14

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

HPD(p, x∗ ) = 0.21

HPD(p, x∗ ) = 0.79

Results for twelve sensor problem

HPD(p, x∗ ) = 0.98

1

C:8.26%
B:10.6%
A:12.4%

CDF

0.8
0.6
0.4
0.2

x

∗

x

∗

x

0

∗

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

HPD

Fig. 7. Example of HPD value for different true target locations (x∗ ). HPD
is the area of probability that is more likely than a given point, i.e., the area
of the shaded region in each diagram.

1
0.8

CDF

the area of the shaded region. Thus if HPD(p, x∗ ) ≈ 1 then x∗
0.6
is in the distant tails of p(x) and is assigned very low likelihood,
0.4
as in figure 6(a), (c), and (d) while if HPD(p, x∗ ) = 0 then x∗
0.2
falls on the most likely value of p(x) (i.e., it is a MAP estimate).
0
Under mild conditions, it can be shown that if x∗ ∼ p(x)
−1
−0.5
0
0.5
1
1.5
2
(i.e., if p(x) correctly characterises the uncertainty in x∗ ) then
Entropy
∗
HPD(p, x ) ∼ U{0, 1} [62, section 9.7.2]. If the distribution
JPDA-BP
JPDA-PBP
CV
CVS(0.55)
MSBP
of HPD values is concentrated at the lower end, then the beliefs
generated are conservative, i.e., they overestimate uncertainty Fig. 8. CDF of HPD value of true location and entropy for beliefs of each
in such a way that the true value rarely falls in the tails. If the target over 200 Monte Carlo trials. Points A, B and C mark the percent of
distribution of HPD values is concentrated at the higher end, cases in which the true target location is less likely than 99% of the marginal
estimates produced by MSBP, JPDA-BP and JPDA-PBP, i.e., it is
then the beliefs are non-conservative, i.e., they underestimate distribution
in the tails in a similar manner to figure 6(a), (c), (d), (f) and (h).
uncertainty, and the true value is often falling in the tails.
The entropy of the distribution characterises its uncertainty;
for example, the entropy of a multivariate Gaussian distribution of the MSBP beliefs are significantly smaller than the other
with covariance P is 0.5 log |2πeP|. Entropy is often preferred methods, the point labelled as ‘A’ in the top figure reveals
over variance for multi-modal distributions as it is not affected that the true location is less likely than 99% of the belief for
by the distance between well-spaced modes (whereas the 12.4% of targets (treating each target in each Monte Carlo
distance between the modes will dominate the variance). The simulation as a sample). This indicates that if the MSBP belief
HPD value is not sufficient, e.g., since an estimator based is used to construct a 99% confidence region for the location
purely on the prior distribution should produce a uniform HPD of a particular target, the target does not lie within the region
distribution, but this would have a much higher entropy than a 12.4% of the time. An instance of this is illustrated in figure
solution which utilises all available measurement data. Likewise 6(d), where the beliefs effectively rule out the location of the
entropy is not sufficient, since one could devise a method true target. These results may be useful in applications where
of approximating beliefs which reports an arbitrarily small smaller entropy is desirable, but when consistency and accuracy
uncertainty; this would report large HPD values. There is not of the beliefs is essential, it is unacceptable.
JPDA-BP and JPDA-PBP also produce non-conservative
a single measure which adequately characterises performance
in this class of problem. This pair of values is necessary results, since the CDF of the HPD value consistently lies
to capture the undesirability of the behaviour in figure 6(a), below the x = y line (i.e., the CDF of a uniform distribution).
(c) and (d) (assigning near-zero likelihood to the true target The points labelled as ‘B’ and ‘C’ reveal that the true location
location, and producing a HPD(p, x∗ ) ≈ 1), as well as the is less likely than 99% of the belief for 10.6% of targets for
undesirability of the behaviour in figure 6(g) (not resolving JPDA-BP, and 8.26% of targets for JPDA-PBP. Again, this
uncertainty when adequate information exists to do so, and indicates that if the respective beliefs are used to construct
99% confidence regions for a particular target, the target does
thus increasing entropy).
For the quantitative experiment, twelve sensors are spaced not lie within the region for 10.6% or 8.26% of the time, as
equally around a circle with radius 100 units. False alarms illustrated in figures 6(a) and (c).
The convex variational (CV) method is shown to produce
follow a Poisson distribution with one per scan on average, and
targets are detected with probability 0.9. The number of targets conservative beliefs, since the CDF of the HPD value consisis Poisson distributed with expected value of twelve (but in tently lies above the x = y line. The cost of this conservatism
each simulation the true number is known by the estimator; the is beliefs with higher entropy; in many applications this may
method can be extended to accommodate an unknown number be preferable. The CVS method with γs = 0.55 is also shown
to significantly reduce instances where the target is in a very
of targets using [10]).
The results in figures 8–9 show the cumulative distribution low likelihood area of the belief. By tuning the value of γs a
function (CDF) of the HPD value (top) and entropy (bottom) trade-off between conservatism and entropy can be obtained.
for the various methods. Figure 8 shows that MSBP produces
A large family of heuristic methods can be developed by
significantly non-conservative results. Although the entropies employing the algorithm in figure 5 with weights that do not

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

Results for heuristic algorithms

proportion of the time. The result is a scalable, reliable
algorithm for estimating marginal probability distributions using
multiple scans.

1
0.8

CDF

15

0.6

ACKNOWLEDGEMENTS

0.4

The authors would like to thank the anonymous reviewers
for suggestions that helped to clarify many points.

0.2
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

HPD
1

CDF

0.8
0.6
0.4
0.2
0
−1

0

−0.5

0.5

1

1.5

2

Entropy
CV

CVH

CVS(0.55)

CVSH(0.55)

CVSH1(0.55)

Fig. 9. CDF of HPD value of true location and entropy for beliefs of each
target over 200 Monte Carlo trials.

ensure that each component (62), (63) is convex. We consider
an instance of this, which sets κf,x = κf,s = 1, κs,1,x =
−1,
P and κs,1,s = κs,2,x = κs,2,s = 0. As long as κf,x +
s∈S [κs,1,x + κs,2,x ] = −|S| + 1 and κf,s + κs,1,s + κs,2,s =
1 ∀ s ∈ S, the original objective remains unchanged, and if the
algorithm converges, the result is optimal (assuming convexity).
Experimentally, convergence appears to be both reliable and
rapid, though not guaranteed. In the rare case that convergence
is not obtained, weights may be reverted (immediately or via
a homotopy) to the form in lemma 4, for which convergence
is guaranteed, but slower in practice. With γs = 1 ∀ s ∈ S,
this can be seen to be
Pequivalent to multiple scan BP (which
is non-convex since s∈S γs > 1).
Figure 9 shows the results of the heuristic approach. The
algorithms with guaranteed convergence are marked as CV and
CVS(0.55); the heuristic equivalents are CVH and CVSH(0.55)
respectively. CVSH1(0.55) uses a single backward-forward
sweep after introducing each new sensor. The slight difference
between the method with guaranteed performance and the
heuristic method is caused by the different values of βs used
(since we must ensure that β̃s > 0.5; in each case we select
βs to set β̃s = 1). The results demonstrate that very similar
performance can be obtained with a single sweep.
V. C ONCLUSION
This paper has shown how the BP data association method
of [32] can be extended to multiple scans in a manner which
preserves convexity, using convex optimisation alongside a
convergent, BP-like method for optimising the FFE. In doing
so, it was demonstrated that the conservative beliefs can be
obtained, whereas the estimates provided by MSBP and JPDABP are significantly non-conservative, and can provide beliefs
which effectively rule out the true target location a significant

R EFERENCES
[1] S. S. Blackman and R. Popoli, Design and Analysis of Modern Tracking
Systems. Norwood, MA: Artech House, 1999.
[2] D. B. Reid, “An algorithm for tracking multiple targets,” IEEE Trans.
Autom. Control, vol. AC-24, no. 6, pp. 843–854, December 1979.
[3] T. Kurien, “Issues in the design of practical multitarget tracking
algorithms,” in Multitarget-Multisensor Tracking: Advanced Applications,
Y. Bar-Shalom, Ed. Norwood, MA: Artech-House, 1990, pp. 43–83.
[4] S. S. Blackman, “Multiple hypothesis tracking for multiple target
tracking,” IEEE Aerospace and Electronic Systems Magazine, vol. 19,
pp. 5–18, Jan. 2004.
[5] T. Fortmann, Y. Bar-Shalom, and M. Scheffe, “Sonar tracking of multiple
targets using joint probabilistic data association,” IEEE J. Ocean. Eng.,
vol. 8, no. 3, pp. 173–184, Jul 1983.
[6] L. Y. Pao, “Multisensor multitarget mixture reduction algorithms for
tracking,” Journal of Guidance, Control, and Dynamics, vol. 17, no. 6,
pp. 1205–1211, 1994.
[7] J. Vermaak, S. Maskell, and M. Briers, “A unifying framework for multitarget tracking and existence,” in Proc. 8th International Conference on
Information Fusion, July 2005.
[8] P. Horridge and S. Maskell, “Real-time tracking of hundreds of targets
with efficient exact JPDAF implementation,” in Proc. 9th International
Conference on Information Fusion, July 2006.
[9] D. Musicki and R. Evans, “Multiscan multitarget tracking in clutter
with integrated track splitting filter,” IEEE Trans. Aerosp. Electron. Syst.,
vol. 45, no. 4, pp. 1432–1447, October 2009.
[10] J. L. Williams, “Marginal multi-Bernoulli filters: RFS derivation of MHT,
JIPDA and association-based MeMBer,” IEEE Trans. Aerosp. Electron.
Syst., vol. 51, no. 3, July 2015.
[11] J. Roecker, “Multiple scan joint probabilistic data association,” IEEE
Trans. Aerosp. Electron. Syst., vol. 31, no. 3, pp. 1204–1210, Jul. 1995.
[12] K. Pattipati, R. Popp, and T. Kirubarajan, “Survey of assignment
techniques for multitarget tracking,” in Multitarget-Multisensor Tracking:
Applications and Advances, Y. Bar-Shalom and W. D. Blair, Eds.
Norwood, MA: Artech-House, 2000, vol. 3, ch. 2, pp. 77–159.
[13] A. B. Poore and S. Gadaleta, “Some assignment problems arising from
multiple target tracking,” Mathematical and Computer Modelling, vol. 43,
no. 9–10, pp. 1074–1091, 2006.
[14] M. J. Wainwright and M. I. Jordan, “Graphical models, exponential
families, and variational inference,” Foundations and Trends in Machine
Learning, vol. 1, no. 1–2, pp. 1–305, 2008.
[15] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles
and Techniques. Cambridge, MA, USA: MIT Press, 2009.
[16] J. Pearl, Probabilistic Reasoning in Intelligent Systems. San Francisco,
CA: Morgan Kaufmann, 1988.
[17] J. S. Yedidia, W. T. Freeman, and Y. Weiss, “Understanding belief
propagation and its generalizations,” Exploring artificial intelligence in
the new millennium, pp. 239–269, 2003.
[18] T. S. Jaakkola, “Tutorial on variational approximation methods,” in
Advanced mean field methods: theory and practice. MIT Press, 2000,
pp. 139–160.
[19] G. E. Kirkelund, C. N. Manchón, L. P. Christensen, E. Riegler, and B. H.
Fleury, “Variational message-passing for joint channel estimation and
decoding in MIMO-OFDM,” in Proc. 2010 Global Telecommunications
Conference, 2010.
[20] E. Riegler, G. E. Kirkelund, C. N. Manchón, M.-A. Badiu, and B. H.
Fleury, “Merging belief propagation and the mean field approximation:
A free energy approach,” IEEE Trans. Inf. Theory, vol. 59, no. 1, pp.
588–602, 2013.
[21] M. Wainwright, T. Jaakkola, and A. Willsky, “A new class of upper
bounds on the log partition function,” IEEE Trans. Inf. Theory, vol. 51,
no. 7, pp. 2313–2335, July 2005.
[22] T. Hazan and A. Shashua, “Norm-product belief propagation: Primal-dual
message-passing for approximate inference,” IEEE Trans. Inf. Theory,
vol. 56, no. 12, pp. 6294–6316, December 2010.

16

[23] R. McEliece, D. MacKay, and J.-F. Cheng, “Turbo decoding as an
instance of Pearl’s “belief propagation” algorithm,” IEEE J. Sel. Areas
Commun., vol. 16, no. 2, pp. 140–152, Feb 1998.
[24] L. Chen, M. J. Wainwright, M. Çetin, and A. S. Willsky, “Multitargetmultisensor data association using the tree-reweighted max-product
algorithm,” in Proc SPIE Signal Processing, Sensor Fusion, and Target
Recognition, vol. 5096, August 2003, pp. 127–138.
[25] L. Chen, M. Çetin, and A. S. Willsky, “Distributed data association
for multi-target tracking in sensor networks,” in Proc. 8th International
Conference on Information Fusion, July 2005.
[26] L. Chen, M. J. Wainwright, M. Çetin, and A. S. Willsky, “Data association
based on optimization in graphical models with application to sensor
networks,” Mathematical and Computer Modelling, vol. 43, no. 9–10,
pp. 1114–1135, 2006.
[27] A. Gning and L. Mihaylova, “Dynamic clustering and belief propagation
for distributed inference in random sensor networks with deficient links,”
in Proc. 12th International Conference on Information Fusion, July 2009,
pp. 656–663.
[28] M. Chertkov, L. Kroc, and M. Vergassola, “Belief propagation and beyond
for particle tracking,” arXiv, e-print arXiv:0806.1199v1, June 2008.
[29] B. Huang and T. Jebara, “Approximating the permanent with belief
propagation,” arXiv, e-print arXiv:0908.1769v1, August 2009.
[30] J. L. Williams and R. A. Lau, “Data association by loopy belief
propagation,” in Proc. 13th International Conference on Information
Fusion, Edinburgh, UK, July 2010.
[31] M. Chertkov, L. Kroc, F. Krzakala, M. Vergassola, and L. Zdeborov,
“Inference in particle tracking experiments by passing messages between
images,” Proceedings of the National Academy of Sciences, vol. 107,
no. 17, pp. 7663–7668, 2010.
[32] J. L. Williams and R. A. Lau, “Approximate evaluation of marginal
association probabilities with belief propagation,” IEEE Trans. Aerosp.
Electron. Syst., vol. 50, no. 4, October 2014.
[33] P. O. Vontobel, “The Bethe permanent of a non-negative matrix,” in Proc.
48th Allerton Conference on Communication, Control, and Computing,
Urbana-Champaign, IL, September/October 2010, pp. 341–346.
[34] J. L. Williams and R. A. Lau, “Convergence of loopy belief propagation
for data association,” in Proc. 6th International Conference on Intelligent Sensors, Sensor Networks and Information Processing, Brisbane,
Australia, December 2010, pp. 175–180.
[35] P. Vontobel, “The Bethe permanent of a nonnegative matrix,” IEEE Trans.
Inf. Theory, vol. 59, no. 3, pp. 1866–1901, 2013.
[36] K. P. Murphy, Y. Weiss, and M. I. Jordan, “Loopy belief propagation for
approximate inference: An empirical study,” in Proc. 15th Conference
on Uncertainty in Artificial Intelligence, 1999, pp. 467–476.
[37] R. A. Lau and J. L. Williams, “Tracking a coordinated group using
expectation maximisation,” in Proc. 8th International Conference on
Intelligent Sensors, Sensor Networks and Information Processing, Melbourne, Australia, April 2013.
[38] R. D. Turner, S. Bottone, and B. Avasarala, “A complete variational
tracker,” in Advances in Neural Information Processing Systems 27, 2014,
pp. 496–504.
[39] R. A. Lau and J. L. Williams, “A structured mean field approach for
existence-based multiple target tracking,” in Proc. 19th International
Conference on Information Fusion, July 2016.
[40] A. S. Rahmathullah, R. Selvan, and L. Svensson, “A batch algorithm for
estimating trajectories of point targets using expectation maximization,”
IEEE Trans Signal Process, vol. 64, no. 18, pp. 4792–4804, Sept 2016.
[41] H. Lan, Q. Pan, F. Yang, S. Sun, and L. Li, “Variational Bayesian
approach for joint multitarget tracking of multiple detection systems,” in
Proc 19th International Conference on Information Fusion, July 2016,
pp. 1260–1267.
[42] F. Meyer, P. Braca, P. Willett, and F. Hlawatsch, “Scalable multitarget
tracking using multiple sensors: A belief propagation approach,” in Proc.
18th International Conference on Information Fusion, July 2015.
[43] ——, “A scalable algorithm for tracking an unknown number of targets
using multiple sensors,” IEEE Trans Signal Process, vol. 65, no. 13, pp.
3478–3493, July 2017.
[44] A. Frank, P. Smyth, and A. Ihler, “Beyond MAP estimation with the
track-oriented multiple hypothesis tracker,” IEEE Trans. Signal Process.,
vol. 62, no. 9, pp. 2413–2423, May 2014.
[45] J. L. Williams, “Interior point solution of fractional Bethe permanent,”
in Proc. IEEE Workshop on Statistical Signal Processing, Gold Coast,
Australia, July 2014.
[46] M. Chertkov and A. B. Yedidia, “Approximating the permanent with
fractional belief propagation,” Journal of Machine Learning Research,
vol. 14, pp. 2029–2066, 2013.

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

[47] S. L. Lauritzen, Graphical Models. Oxford, UK: Clarendon Press, 1996,
vol. 17.
[48] R. E. Kalman, “A new approach to linear filtering and prediction
problems,” Transactions of the ASME Journal of Basic Engineering,
vol. 82, no. Series D, pp. 35–45, 1960.
[49] L. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition,” Proceedings of the IEEE, vol. 77, no. 2,
pp. 257–286, Feb 1989.
[50] A. Viterbi, “Error bounds for convolutional codes and an asymptotically
optimum decoding algorithm,” IEEE Trans. Inf. Theory, vol. 13, no. 2,
pp. 260–269, April 1967.
[51] T. M. Cover and J. A. Thomas, Elements of Information Theory. New
York, NY: John Wiley and Sons, 1991.
[52] A. Globerson and T. Jaakkola, “Approximate inference using conditional entropy decompositions,” Journal of Machine Learning Research:
Workshop and Conference Proceedings, vol. 2, pp. 131–138, 2007.
[53] H. A. Bethe, “Statistical theory of superlattices,” Proceedings of the
Royal Society of London, Series A—Mathematical and Physical Sciences,
vol. 150, no. 871, pp. 552–575, 1935.
[54] J. L. Williams and R. A. Lau, “Multiple scan data association by convex
variational inference (extended version),” arXiv, e-print arXiv:1607.07942,
January 2018.
[55] D. Musicki and R. J. Evans, “Joint integrated probabilistic data association: JIPDA,” IEEE Trans. Aerosp. Electron. Syst., vol. 40, no. 3, pp.
1093–1099, July 2004.
[56] R. T. Rockafellar, Convex Analysis. Princeton, NJ: Princeton University
Press, 1970.
[57] P. Tseng, “Dual coordinate ascent methods for non-strictly convex
minimization,” Mathematical Programming, vol. 59, no. 1–3, pp. 231–
247, 1993.
[58] A. L. Yuille, “CCCP algorithms to minimize the Bethe and Kikuchi
free energies: Convergent alternatives to belief propagation,” Neural
Computation, vol. 14, no. 7, pp. 1691–1722, June 2002.
[59] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
Numerical Recipes: The Art of Scientific Computing, 3rd ed. New York,
NY, USA: Cambridge University Press, 2007.
[60] H. A. Blom and E. A. Bloem, “Probabilistic data association avoiding
track coalescence,” IEEE Trans. Autom. Control, vol. 45, no. 2, pp.
247–259, February 2000.
[61] J. L. Williams, “An efficient, variational approximation of the best fitting
multi-Bernoulli filter,” IEEE Trans. Signal Process., vol. 63, no. 1, pp.
258–273, January 2015.
[62] S. Davey, N. Gordon, I. Holland, M. Rutten, and J. Williams, Bayesian
Methods in the Search for MH370, ser. SpringerBriefs in Electrical and
Computer Engineering. Singapore: Springer, 2016.
[63] S. Julier and J. Uhlmann, “Unscented filtering and nonlinear estimation,”
Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, Mar. 2004.
[64] N. Gordon, D. J. Salmond, and A. Smith, “Novel approach to non-linear
and non-Gaussian Bayesian state estimation,” IEE Proceedings F: Radar
and Signal Processing, vol. 140, pp. 107–113, 1993.
[65] D. P. Bertsekas, Nonlinear Programming, 2nd ed. Belmont, MA: Athena
Scientific, 1999.
[66] E. Kohlberg and J. W. Pratt, “The contraction mapping approach to
the Perron-Frobenius theory: Why Hilbert’s metric?” Mathematics of
Operations Research, vol. 7, no. 2, pp. 198–210, 1982.

A PPENDIX A
A SSOCIATION HISTORY MODEL
In section II-B, the problem of interest is formulated to
incorporate both continuous states xi and discrete association
hypothesis variables ais . Alternatively, we may formulate the
problem by defining association history hypotheses aiS =
(ai1 , . . . , aiS ), which detail which measurement corresponds to
the target in each scan. The role of the variational algorithm is to
determine the marginal association distribution pi (aiS ) for each
target. Calculation of the kinematic distribution conditioned on
i
an association hypothesis, pi,aS (xi ), can utilise well-studied
methods such as the Kalman filter (KF) [48], extended Kalman
filter (EKF), unscented Kalman filter (UKF) [63], or the particle
filter (PF) [64].

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

17

replaced by S 0 = {1, . . . , S 0 }, where S 0 = S + 1. If the new
scan represents a new time step, each hypothesis-conditioned
i
PDF pi,aS (xi ) first undergoes prediction according to standard
KF/EKF/UKF/PF expressions. A new single-target hypothesis
aiS 0 = (aiS , aiS 0 ) is generated for each combination of an old
single-target hypothesis aiS , and choice of event in the new
Definition 2. A single target association history hypothesis
scan, aiS 0 , where aiS 0 = 0 denotes a missed detection, and
(or single target hypothesis for short) is a sequence of timeaiS 0 = j ∈ {1, . . . , mS 0 } indicates that target i corresponded
stamped measurements that are hypothesised to correspond to
to measurement j. The parameters for the hypothesis aiS 0 =
the same target.
(aiS , 0) can be calculated using the expression:
Z
A global hypothesis for the scans in set S may be represented
i
i
i
wi,aS 0 = wi,aS [1 − PSd0 (xi )]pi,aS (xi )dxi ,
(94)
as aS = (a1S , . . . , anS ). Each single target hypothesis is
i,aiS
equipped with a hypothesis weight w
(utilised in the
i
i
pi,aS 0 (xi ) ∝ [1 − PSd0 (xi )]pi,aS (xi ).
(95)
calculation of the probability of the global hypotheses), and the
i
i
target state probability density function (PDF) conditioned on The hypothesis aS 0 = (aS , j), which updates the old singlei
the hypothesis pi,aS (xi ). Under this model, prediction steps target hypothesis aiS with measurement z jS 0 , is calculated using
may be easily introduced to incorporate a stochastic state the expressions:
model.
i R
i
i
wi,aS pS 0 (z jS 0 |xi )PSd0 (xi )pi,aS (xi )dxi
We denote by AiS the set of feasible single-target hypotheses
wi,aS 0 =
,
j
λfa
for target i ∈ {1, . . . , n} in the scans in S. The set of all
S 0 (z S 0 )
(96)
feasible global hypotheses (i.e., those in which no two targets
i
i
j
i,a
i
i
d
i
i,a
i
utilise the same measurement) can be written as:
p S 0 (x ) ∝ pS 0 (z S 0 |x )PS 0 (x )p S (x ).
(97)

The extension of these steps to accommodate an unknown,
AS = (a1S , . . . , anS ) aiS = (ai1 , . . . , aiS ) ∈ AiS ,
time-varying number of targets can be found in [10].

The association history model can be written in a graphical
ais 6= ajs ∀ s, i, j s.t. i 6= j, ais 6= 0 . (88) model form as:
Definition 1. A global association history hypothesis (or
global hypothesis for short) is a hypothesis for the origin of
every measurement received so far, i.e., for each measurement
it specifies from which target it originated, or if it was a false
alarm.

The joint distribution of states and hypotheses conditioned
p(X, aS , bS |ZS ) ∝



on measurements may be written as:
ms
n 

Y
Y
( n
)(
) Y
i
i
i
i,j i j 
ψsi (aiS , ais )
Y
Y
ψ
(x
,
a
)ψ
(a
)
ψ
(a
,
b
)
,
i
i
S
S
s
s
s

p(X, aS , bS |ZS ) ∝
wi,aS pi,aS (xi )
ψs (as , bs ) . i=1 
j=1
s∈S
i=1

s∈S

(89)
Marginalising the kinematic states, the probability of a global
hypothesis aS = (a1S , . . . , anS ) ∈ AS (and the corresponding
bS ) can be written in the form:
( n
)(
)
Y
Y
i,aiS
p(aS , bS |ZS ) ∝
w
ψs (as , bs ) .
(90)
i=1

s∈S

The joint PDF of all targets can be represented through a
total probability expansion over all global hypotheses:
p(X) =

X

p(aS )

aS ∈AS

n
Y

i

pi,aS (xi ),

(91)

i=1

where, for notational simplicity, we drop the explicit conditioning on ZS from p(X) and p(aS ). It is of interest to obtain
the marginal distributions of global hypothesis probabilities:
X
pi (aiS ) =
p(ãS ).
(92)
i
i
ãS =(ã1S ,...,ãn
S )∈AS |ãS =aS

From these marginal association distributions, we can find the
marginal state PDF of each target:
X
i
pi (xi ) =
pi (aiS )pi,aS (xi ).
(93)
aiS ∈AiS

We now describe the updates which occur when a new scan
of measurements is received, i.e., when S = {1, . . . , S} is

(98)
where ψ i (xi , aS ) =
ψsi (aiS , ais ) ensures that
(
1,
i
i
i
ψs (aS , as ) =
0,

i,aiS

i,aiS

p
(xi ), ψ i (aS ) = w
, and
i
aS and ais are in agreement:
aiS = (ãi1 , . . . , ãiS ), ãis = ais
otherwise.

(99)

This graph is illustrated in figure 10. Marginalising the
kinematic states xi (which can be done simply since they
are leaves), we arrive at the representation
p(aS , bS |ZS ) ∝



ms
n 

Y
Y
Y
ψsi (aiS , ais )
ψ i (aS )
ψsi,j (ais , bjs ) . (100)


i=1

s∈S

j=1

The derivation of section III may be applied directly to the
model in (100), substituting aiS in place of xi . The kinematic
distribution can then be recovered as
i

q(xi , aiS ) = q(aiS )ψ i (xi , aS ) = q(aiS )pi,aS (xi ).

(101)

A PPENDIX B
D ERIVATION OF B ETHE FREE ENERGY FORM
A. Single scan
In this section, we present two formulations for the single
scan problem, (102)-(105) and (117)-(120), and show they

18

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

function for the partial minimisation can be written as:
minimise
FB ([qsi (ais )], [qsi,j (ais , bjs )], [qsj (bjs )])
i,j
j
i
qs (as ,bs )


ms X
ms
n X
n
X
X
i 
+
λi,j
qsi,j (ais , bjs ) − qsi (ais )
s (as )
i=1 j=1 ais =0

+

bjs =0



ms X
n X
n
X

j 
λi,j
s (bs )

i=1 j=1 bjs =0

Fig. 10. Graphical model of association history formulation of multiple scan
problem.

are equivalent to the formulation of (27)-(30) after partial
minimisation. The formulation of (27)-(30) is similar to the
formulation in [35]; the difference is our formulation includes
the belief that target i is not detected, qsi,0 , and the belief
that measurement zsj is not used by any target, qs0,j , whereas
the formulation in [35] is for the matrix permanent problem,
which excludes qsi,0 and qs0,j (and is thus constrained such that
n = ms ).
The Bethe variational problem in section II-C involving
random variables ais and bjs can be solved by minimising:

i=1

−

qsi,j (ais , bjs ) = 1,

qsi,j (ais , bjs ) ≥ 0,

bjs =0

i
i,j j
where λi,j
s (as ) and λs (bs ) are dual variables. Solving the
dual function yields the solution:

qsi,j (ais , bjs ) =

1
ci,j

i
i,j j
ψsi,j (ais , bjs ) exp{−λi,j
s (as ) − λs (bs )},

(106)
where ci,j is the normalisation constant.
Using the marginalisation constraints (104), the pairwise
joint (106) and the definition of ψsi,j (ais , bjs ) in (26) (which
ensures that qsi,j (ais , bjs ) = 0 if ais = j, bjs 6= i or bjs = i,
ais 6= j), we find that qsi,j = qsi,j (ais = j, bjs = i), which is
related to the dual variables by
qsi,j =

FB ([qsi (ais )], [qsi,j (ais , bjs )], [qsj (bjs )]) =
ms
n
X
X

−
H(ais ) + E[log ψsi (ais )] −
H(bjs )
ms
n X
X


ais =0

qsi,j (ais , bjs ) − qsj (bjs )

ais =0

ms X
n
X

subject to



ms
X

1
ci,j

i
i,j j
exp{−λi,j
s (as = j)} exp{−λs (bs = i)}. (107)

Secondly, for i0 6= i and j 0 6= j, we find that qsi (ais = j 0 ) = qsi,j
0
and qsj (bjs = i0 ) = qsi ,j are related to the dual variables by

j=1

−I(ais ; bjs ) + E[log ψsi,j (ais , bjs )] , (102)

0

qsi,j =

i=1 j=1

1
ci,j

i
0
exp{−λi,j
s (as = j )}

n
X

0

j
0
exp{−λi,j
s (bs = i )},

i0 =0
0
i 6=i

(108)

subject to the constraints:
0

qsi,j (ais , bjs ) ≥ 0, qsi (ais ) ≥ 0, qsj (bjs ) ≥ 0,
(103)
m
n
s
X
X
qsi,j (ais , bjs ) = qsi (ais ),
qsi,j (ais , bjs ) = qsj (bjs ),
bjs =0
n
X

ais =0 bjs =0

1
j
0
exp{−λi,j
s (bs = i )}
ci,j

qsi,j (ais , bjs ) = 1,

ms
X
ais =0

qsi (ais ) = 1,

n
X

qsj (bjs ) = 1,

j =0
j 0 6=j

(109)
As we have already seen, qsi,j = qsi,j (ais = j, bjs = i). In the
remaining case, ais = j 0 6= j, bjs = i0 6= i, we again exploit the
structure of ψsi,j (ais , bjs ) to obtain

bjs =0

(105)
where ψsi (ais ) is defined by (32), I(ais ; bjs ) is defined in (6),
and ψsi,j (ais , bjs ) is in (26). Note that there is some redundancy
in these constraints, which is retained to reinforce that the
constraints in (103) and (105) are retained when the marginal
constraints (104) are relaxed in the next step.
Let the marginals qsi (ais ) and qsj (bjs ) be fixed and feasible.
Because the marginals are fixed, the Bethe variational problem
(102)-(105) is convex with respect to qsi,j (ais , bjs ). In addition,
since the marginals are feasible, then qsi (ais = j) = qsj (bjs =
i) , qsi,j . Relaxing the marginal constraints (104), the dual

i
0
exp{−λi,j
s (as = j )}.

0

ais =0

(104)
ms
X

qsi ,j =

ms
X

1

i
0
i,j j
0
exp{−λi,j
s (as = j )} exp{−λs (bs = i )}.
(110)
Substituting (107)–(110) into the pairwise normalisation constraint (105) yields:

qsi,j (ais , bjs ) =

qsi,j +

ci,j

n
1 X

ci,j

j
0
exp{−λi,j
s (bs = i )}

i0 =0
0
i 6=i

×

ms
X
j 0 =0
j 0 6=j

i
0
exp{−λi,j
s (as = j )} = 1. (111)

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

Subsequently, we find for i0 6= i, j 0 6= j,

rangement of the Bethe free energy (33) is:
0

qsi,j (ais

=j

0

, bjs

0

=i)=

19

0

qsi,j qsi ,j
1 − qsi,j

.

(112)

FB ([qsi (ais )], [qsi (xi , ais )], [qsi,j ]) =
n
X

−
H(ais ) + H(xi |ais ) + E[log ψ i (xi )ψsi (xi , ais )]
i=1

Substituting the marginals and the pairwise joint into the
entropies H(ais ), H(bjs ) and H(ais , bjs ) yields:

+

ms
X

qs0,j log qs0,j −

j=1

−H(ais ) =

ms
X

qsi,j log qsi,j ,

(113)

qsi,j log qsi,j ,

(114)

j=0

−H(bjs ) =

n
X

−H(ais , bjs ) = qsi,j log qsi,j +

ms
n X
X
qsi,j qsi ,j
i0 =0 j 0 =0
i0 6=i j 0 6=j

1−

0

qsi,j

0

log

0

qsi,j qsi ,j
1−

qsi,j

= qsi,j log qsi,j − (1 − qsi,j ) log (1 − qsi,j )
ms
n
X
X
0
0
0
0
+
qsi,j log qsi,j +
qsi ,j log qsi ,j (115)
j 0 =0
j 0 6=j

i0 =0
i0 6=i

qsi,j

Let
and
be fixed and feasible. Minimising the
Bethe free energy (121) with respect to qsi (xi , ais ) subject to
the constraints (34)-(38) yields the solution (using (13)):
q i (ai )ψ i (xi )ψ i (xi , ai )
qsi (xi , ais ) = Ps s i i 0 is i 0 si .
xi 0 ψ (x )ψs (x , as )

Substituting qsi (xi , ais ) into the Bethe free energy (121) results
in the equivalent
Bethe variational problem (27)-(30) where
P
wsi,j = xi ψ i (xi )ψsi (xi , j), and qsi (ais = j) = qsi,j .
B. Multiple scans
The Bethe variational problem in section II-B, which is
represented by figure 3(c), can be solved by minimising:

i=1

I(ais ; bjs ) = −qsi,j log qsi,j − (1 − qsi,j ) log (1 − qsi,j ). (116)
−
Substituting (113)-(116) into the single scan formulation (102)(105), we arrive at the equivalent Bethe variational problem
(27)-(30) where wsi,j = ψsi (ais = j).
The Bethe variational problem in section II-C involving
random variables xi , ais and bjs (illustrated in figure 3(a)) can
be solved by minimising:
FB ([q i (xi )], [qsi (xi , ais )], [qsi (ais )], [qsi,j (ais , bjs )], [qsj (bjs )]) =
ms
n
n
X
X
X

−
H(xi ) + E[log ψ i (xi )] −
H(ais ) −
H(bjs )
i=1

−

n
X

j=1


−I(xi ; ais ) + E[log ψsi (xi , ais )]

i=1

−

ms
n X
X


−I(ais ; bjs ) + E[log ψsi,j (ais , bjs )] , (117)

i=1 j=1

subject to the constraints (103)-(105) and

ms
X

qsi (xi , ais ) ≥ 0, qsi (xi ) ≥ 0,
(118)
X
qsi (xi , ais ) = qsi (xi ),
qsi (xi , ais ) = qsi (ais ), (119)

ais =0

xi

X

ms
X

xi ais =0

qsi (xi , ais )

= 1,

(122)

FB ([q i (xi )], [qsi (xi , ais )], [qsi (ais )], [qsi,j (ais , bjs )], [qsj (bjs )]) =
n
X

−
H(xi ) + E[log ψ i (xi )]

so that the mutual information (6) is:

i=1

(1 − qsi,j ) log(1 − qsi,j ). (121)

i=1 j=1

qsi (ais )

i=0
0

ms
n X
X

n
XX

H(ais ) −

s∈S i=1

ms
XX

H(bjs )

s∈S j=1

n
XX

−
−I(xi ; ais ) + E[log ψsi (xi , ais )]
s∈S i=1
ms
n X
XX

−
−I(ais ; bjs ) + E[log ψsi,j (ais , bjs )] , (123)
s∈S i=1 j=1

subject to the constraints (103)-(105) and (118)-(120) for s ∈
S. Partial minimisation over qsi,j (ais , bjs ) and rearrangement
produces the equivalent variational problem (51)-(57).
A PPENDIX C
P ROOF OF ALGORITHMS FOR MINIMISING PDCA
SUB - PROBLEMS
In this section, we prove lemmas 5 and 6, i.e., we derive
algorithms for minimising the blocks utilised in the PDCA
algorithm. Before we begin, we prove the preliminary result in
lemma 3, which shows that the block hs,1 is convex (convexity
of f and hs,2 is straight-forward).
Proof of lemma 3: If κs,1,x ≥ 0, then convexity with
respect to q i (xi ) is immediate. Otherwise, −κs,1,x > 0; let
κ̃ = κs,1,s + κs,1,x ≥ 0, and rewrite the first line of (62) as:
n
X

−
−κs,1,x H(ais |xi ) + κ̃H(xi |ais ) + κ̃H(ais ) .
i=1

X

qsi (xi )

= 1.

(120)

xi

Partial minimisation over the pairwise joint qsi,j (ais , bjs )
arrives at the Bethe variational problem (33)-(38). A rear-

The first two terms are convex by definition of conditional
entropy, as is the second line in (62), so we focus on the
remainder of the expression:
ms
n
n X
X
X
− κ̃
H(ais ) − γs
(1 − qsi,j ) log(1 − qsi,j ), (124)
i=1

i=1 j=1

20

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

Pms i,j
qs log qsi,j . Recognising (124) as:
where −H(ais ) = j=0


ms
ms
n X
n X
X
X
γs 
qsi,j log qsi,j −
(1 − qsi,j ) log(1 − qsi,j )
i=1 j=0

+(κ̃−γs )

while qsi (xi , ais ) can be recovered via (70). Dividing by (κf,s +
κs,1,s ), we obtain (74). Finally, since
∇qi (xi ) f = κf,x log q i (xi ) + κf,x − log ψ i (xi ),

ms
n X
X

qsi,j log qsi,j +γs

i=1 j=0

n
X
(1−qsi,0 ) log(1−qsi,0 ),
i=1

we find that the update in (48) reduces to
λis,1,x (xi ) = −µi (xi ) − κf,x log q i (xi ) − κf,x + log ψ i (xi )

(125)
we obtain the desired result; the first line is convex by
theorem
20 inP
[35], which shows that the function S(ξ) =
P
ξ
−
j
j ξj log
j (1 − ξj ) log(1 − ξj ) is convex on the domain
P
ξj ≥ 0, j ξj = 1; the second line is convex by convexity of
x log x.
Proof of lemma 5: Collecting terms, the objective to be
minimised is:
µ
Fs,1
([q i (xi )], [qsi (xi , ais )])
n
X

=−
(κf,x + κs,1,x )H(xi ) + E[φi (xi )]

= φi (xi ) − κf,x log q i (xi ) − κf,x ,

−

−

µ
Fs,2
([q i (xi )], [qτi (xi , aiτ )])
n
X

=−
(κf,x + κs,2,x )H(xi ) + E[φi (xi )]
i=1
n
XX

−
κf,τ H(xi , aiτ ) + E[φiτ (xi , aiτ )]

n
X

κf,τ H(xi , aiτ ) + E[φiτ (xi , aiτ )]

τ ∈S i=1

−

τ ∈S\{s} i=1
n
X
i=1

+ βs

qs0,j log qs0,j

ms
n X
X
− γs
(1 − qsi,j ) log(1 − qsi,j ).

j=1

i=1 j=1

(126)
For τ ∈ S\{s}, hs,1 is constant with respect to qτi (xi , aiτ ), so
λis,1,τ (xi , aiτ ) = 0. If we define
φ̃is (xi , ais ) = φi (xi ) + φis (xi , ais ),
ms
X
q i (xi ) =
qsi (xi , ais ),

(127)

(κf,s + κs,1,s )[H(xi |ais ) + H(ais )]

i=1

−

n
X

E[φ̃is (xi , ais )] + βs

i=1

ms
X

ms
n X
X
(1 − qsi,j ) log(1 − qsi,j ). (129)

n
X


κ̃H(xi ) + E[φi (xi )]

i=1
n
X

−
(κf,s + κs,2,s )H(ais |xi ) + E[φis (xi , ais )]
i=1

X

n
X


κf,τ H(xi , aiτ ) + E[φiτ (xi , aiτ )] , (134)

τ ∈S\{s} i=1

where κ̃ = κf,x + κf,s + κs,2,x + κs,2,s . Using lemma 2
to perform a partial minimisation of (134) with respect to
qτi (xi , aiτ ), holding q i (xi ) fixed, we find the result in (80) and
the remaining problem:

qs0,j log qs0,j

j=1

− γs

µ
Fs,2
([q i (xi )], [qτi (xi , ais )]) = −

−

then the terms in (126) that depend on qsi (xi , ais ) can be written
as
−

For τ ∈ S\{s}, (54) is not enforced, and hs,2 is constant
with respect to qτi (xi , aiτ ), so λis,2,τ (xi , aiτ ) = 0. Since the
constraint (54) is enforced for time s, (133) can be written
equivalently as

(128)

ais =0

n
X

n
X

κs,2,s H(xi , ais ) . (133)
i=1


(κf,s + κs,1,s )H(xi , ais ) + E[φis (xi , ais )]

ms
X

(132)

which is the result in (75). Following identical steps for
qsi (xi , ais ) gives the result in (76).
Proof of lemma 6: Collecting terms, the objective to be
minimised is:

i=1

X

(131)

i=1 j=0

c
µ
Fs,2
([q i (xi )]) =

n n
io
h
X
−
κ̃H(xi ) + E φi (xi ) + φ̃is (xi ) .
i=1

(135)

i=1 j=1

Using lemma 2 to minimise with respect to qsi (xi , ais ) while
holding qsi (ais ) fixed, we find that the optimisation becomes
"
n n
o
X
(κf,s + κs,1,s ) −
H(ais ) + E[φ̃is (ais )]
i=1

+ β̃s

ms
X
j=1

qs0,j

log qs0,j

#
ms
n X
X
i,j
i,j
− γ̃s
(1 − qs ) log(1 − qs ) ,
i=1 j=1

(130)

Using lemma 1, we obtain the result in (79).
Following similar steps to (132) gives the updates for
λis,2,x (xi ) and λis,2,s (xi , ais ).
A PPENDIX D
P ROOF OF CONVERGENCE OF SINGLE SCAN ITERATION
In this section, we prove convergence of an iterative
algorithm for solving the single scan block, illustrated in figure

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

Lemma 7. The solution of (136)-(139) lies in the relative
interior, i.e., 0 < qi,j < 1 ∀ (i, j) s.t. wi,j > 0.

2:
minimise

n X
m
X

qi,j log

i=1 j=1
n
X

+α

qi,0 log

i=1

qi,j
wi,j

Proof. Rewrite the objective in (136) in the form:

m
X
q0,j
qi,0
+β
q0,j log
wi,0
w
0,j
j=1

n X
m
X
−γ
(1 − qi,j ) log(1 − qi,j )

(1 − γ)

n X
m
X

subject to

j=0
n
X

qi,j log

i=1 j=1

+ (α − γ)

(136)

i=1 j=1
m
X

21

n
X

qi,0 log

i=1

"

qi,j = 1 ∀ i ∈ {1, . . . , n}

(137)

qi,j = 1 ∀ j ∈ {1, . . . , m}

(138)

qi,j
wi,j
m
X
qi,0
q0,j
+β
q0,j log
wi,0
w0,j
j=1

#
n m
qi,j X X
+γ
−
qi,j log
(1 − qi,j ) log(1 − qi,j ) .
wi,j i=1 j=1
i=0 j=1
n X
m
X

(140)

i=0

0 ≤ qi,j ≤ 1,

(139)

Consider two feasible points q 0 and q 1 , where q 0 is on the
boundary and q 1 is in the relative interior (such a point exists
by assumption 1). Let q λ = λq 1 + (1 − λ)q 0 , and denote the
objective evaluated at q λ by

where γ ∈ [0, α] ∩ [0, β] ∩ [0, 1), α ∈ (0.5, ∞), β ∈ (0.5, ∞).
In our analysis, we permit values wi,j = 0, maintaining a finite
f (λ) = g(λ) + h(λ),
objective by fixing the corresponding qi,j = 0, and defining
qi,j /wi,j , 1; since these take on fixed values, we do not where g(λ) is the first two lines of (140) evaluated at q λ ,
consider them to be optimisation variables. While we state and h(λ) is the final line. Lemma 3 shows that h(λ) is
the algorithm more generally, we prove convergence for three convex, therefore its gradient is monotonically non-decreasing.
cases:
Consequently it must be the case that:
1) n = m and wi,0 = w0,j = 0 (i.e., no missed
lim h0 (λ) = c < ∞.
(141)
detection/false alarm events)
λ↓0
2) wi,0 > 0 ∀ i, w0,j = 0 ∀ j, α = 1 (i.e., missed detections
The derivative of g(λ) is given by:
but no false alarms)
"
#
3) wi,0 > 0 ∀ i, w0,j > 0 ∀ j, α = 1 (i.e., missed detections
n X
m
1
0
X
λq
+
(1
−
λ)q
i,j
i,j
0
1
0
and false alarms)
g (λ) = (1−γ)
(qi,j −qi,j ) log
+1
wi,j
In case 1 above, α and β have no effect since qi,0 = 0 and
i=1 j=1
"
#
n
1
0
q0,j = 0. Similarly, in case 2, β has no effect since q0,j = 0.
X
λqi,0
+ (1 − λ)qi,0
1
0
(qi,0 − qi,0 ) log
+ (α − γ)
+1
Assumption 1 ensures that the problem has a relative interior
wi,0
i=1
(again, we exclude the qi,j variables for which wi,j = 0, since
"
#
m
1
0
X
they are fixed to zero).
λq0,j
+ (1 − λ)q0,j
1
0
+β
(q0,j − q0,j ) log
+1 .
w0,j
Assumption 1. There exists a feasible point in the relative
j=1
interior, i.e., there exists qi,j satisfying the constraints (137)(139) such that 0 < qi,j < 1 ∀ (i, j) s.t. wi,j > 0.
Assumption 2. The graph is connected, i.e., we can travel from
any left-hand side vertex ai , i ∈ {1, . . . , n} to any right-hand
side vertex bj , j ∈ {1, . . . , m} by following a path consisting
of edges (i0 , j 0 ) with wi0 ,j 0 > 0.
Assumption 1 can easily be shown to be satisfied if wi,0 >
0 ∀ i and w0,j > 0 ∀ j (i.e., missed detection and false alarm
likelihoods are non-zero). In problems without false alarms or
missed detections, the condition excludes infeasible problems
(e.g., where two measurements can only be associated with
a single target), and problems with trivial components (e.g.,
where a measurement can only be associated with one target,
so that measurement and target can be removed and the smaller
problem solved via optimisation). Assumption 2 ensures that
the problem is connected; this property is utilised in the proof
of case 1. Any problem in which the graph is not connected
can be solved more efficiently by solving each connected
component separately.

Since q 0 is on the boundary and q 1 is not, we must have:
lim g 0 (λ) = −∞.
λ↓0

(142)

By (141) and (142), we thus have that f 0 (λ) < 0 ∀ λ ∈
(0, ) for some  > 0. Thus the optimum cannot lie on the
boundary.
Lemma 8. The Karush-Kuhn-Tucker (KKT) optimality conditions [65] for the problem in (136) are:
qi,j
log
+ γ log(1 − qi,j ) + 1 + γ − λi − µj = 0 ∀ i, j > 0,
wi,j
(143)
qi,0
α log
+ α − λi = 0 ∀ i > 0,
(144)
wi,0
q0,j
β log
+ β − µj = 0 ∀ j > 0,
(145)
w0,j
as well as the primal feasibility conditions (137)-(139). The
conditions are necessary and sufficient for optimality. In case
1 (where wi,0 = 0 ∀ i and w0,j = 0 ∀ j) the solution is unique

22

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

up to a constant c being added to λi ∀ i and subtracted from
µj ∀ j. In other cases, the solution is unique.
Proof. One complication is that the objective is not convex
on Rn+1×m+1 but rather only on the subspace in which
either (137) or (138) is satisfied. We show that the regular
KKT conditions are still necessary and sufficient in this case.
Relaxing the non-negativity condition,9 the problem can be
expressed as:

The optimisation methodology we adopt is to define an
iterative method and prove that it converges to a point that
satisfies the KKT conditions, motivated by analysis of the BP
iteration in [32], [35]. Defining λ̄i = λi − α, µ̄i = µi − β and
κ = −1 − γ + α + β, the KKT conditions in (143)-(145) can
be rewritten as:

qi,0

minimise f (q)
subject to A1 q = b1 ,

q0,j =

A2 q = b2 .

The KKT conditions for this problem are:
∇f (q) − AT1 λ − AT2 µ = 0,
A1 q = b1 ,

A2 q = b2 .

(146)
(147)

Given a solution q 0 that satisfies A1 q 0 = b1 , we can
express any feasible q as q 0 + P(q − q 0 ) where P =
I − AT1 (A1 AT1 )−1 A1 is the matrix that projects onto the
null-space of A1 . Thus we can equivalently solve

A2 q = b2 ,

(148)
(149)

where, after taking the gradient of f in (148), we substitute
q 0 + P(q − q 0 ) = q since the point must satisfy the constraints
(149). The projection of the gradient is:

−1
m
X
wi,j exp{µ̄j + κ} 
,
exp λ̄i = wi,0 exp{( α1 − 1)λ̄i } +
(1 − qi,j )γ
j=1


"

(155)
where the RHS values of qi,j , λ̄i and µ̄j refer to the previous
iterates. The updates in (154) and (155) are applied alternately.
After each update, the values of qi,j are recalculated using
(151)-(153).
The iteration may be written equivalently in terms of the
parameterisation xi,j and yi,j , where

Similarly, given a point satisfying the KKT conditions for the
original problem, we can find a corresponding point satisfying
the modified KKT conditions (148)-(149) by inverting (150).
Thus points satisfying the KKT conditions for the original
problem (146)-(147) and the modified problem (148)-(149) are
in direct correspondence.
The expressions in (143)-(145) are found by forming the
Lagrangian and taking gradients. Uniqueness of the solution
comes from strict convexity of f . The freedom to choose
a constant offset is the result of linear dependence of the
constraints
in case 1 (since each set of constraints implies that
P
q
=
n).
i,j i,j

exp{µ̄j + κ}
,
(1 − qi,j )γ
= exp{( α1 − 1)λ̄i }, x0,j = exp{µ̄j },
exp{λ̄i + κ}
yi,j =
,
(1 − qi,j )γ
= exp{λ̄i }, y0,j = exp{( β1 − 1)µ̄j }.
xi,j =

xi,0

(150)

#−1

(1 − qi,j )γ

i=1

∗

λ̃ = λ∗ + (A1 AT1 )−1 A1 ∇f (q).

n
X
wi,j exp{λ̄i + κ}

exp µ̄j = w0,j exp{( β1 − 1)µ̄j } +

Thus a point (q , λ , µ ) satisfying the KKT conditions for
the modified problem (148)-(149) corresponds to a point
∗
(q ∗ , λ̃ , µ∗ ) in the KKT conditions for the original problem
(146)-(147), where
∗

(153)

While these expressions do not permit us to immediately solve
for qi,j , they permit application of an iterative method, in which
we repeatedly calculate new LHS values of qi,j by updating
either λ̄i via the equation:

P∇f (q) = ∇f (q) − AT1 (A1 AT1 )−1 A1 ∇f (q).
∗

∀ j > 0.

A2 q = b2 .

P∇f (q) − AT1 λ − AT2 µ = 0,

∗

(152)

or µ̄j via the equation:

Since the argument of f lies in the feasible subspace for the
first constraint, this problem is convex, and under Assumption 1
the Slater condition [65] is satisfied, so the KKT conditions are
necessary and sufficient. The KKT conditions for this modified
problem are:

A1 q = b1 ,

w0,j exp{ β1 µ̄j }

(151)

(154)

minimise f (q 0 + P(q − q 0 ))
subject to A1 q = b1 ,

wi,j exp{λ̄i + µ̄j + κ}
∀ i, j > 0,
(1 − qi,j )γ
= wi,0 exp{ α1 λ̄i } ∀ i > 0,

qi,j =

yi,0

(156)

(157)

Algebraic manipulation yields equivalent iterations in terms of
xi,j and yi,j as:
!−(1−γ)
(k+1)
xi,j

= ri,j (y

(k)

(k)
w0,j y0,j

),

+

X

(k)
wi0 ,j yi0 ,j

i0

−γ


(k)

× w0,j y0,j +

X

(k)

wi0 ,j yi0 ,j 

× eκ , (158)

i0 6=i
(k+1)

xi,0

(k)

1

= ri,0 (y (k) ) , (yi,0 ) α −1 ,

(159)
!−1

(k+1)

x0,j

= r0,j (y (k) ) ,

(k)

w0,j y0,j +

X

(k)

wi0 ,j yi0 ,j

,

i0
9 Alternatively,

define f (q) = ∞ for points violating the constraint.

(160)

,

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

23

thus

Fig. 11. Structure of iterative solution, alternating between half-iterations
x(k+1) = r(y (k) ) and y (k+1) = s(x(k+1) ).

and

i0 6=i

(k+1)

(k+1)

= si,j (x(k+1) ) , wi,0 xi,0

+

X

(k+1) 

Raising this to the power γ and multiplying by exp{(1−γ)µ̄∗j +
κ}, we obtain (158). Similar steps show (161) and (162), while
(159) and (163) are immediate.

wi,j 0 xi,j 0

j0

−γ


(k+1)

× wi,0 xi,0

+

X

(k+1) 

wi,j 0 xi,j 0

× eκ ,

j 0 6=j

(161)
−1


(k+1)

yi,0

(168)

−(1−γ)


yi,j

"
#−1
exp{µ̄∗j }
wi,j exp{λ̄∗i + κ}
∗
(167)
= exp{−µ̄j } −
∗
∗ )γ
1 − qi,j
(1 − qi,j

−1
n
∗
X wi0 ,j exp{λ̄ 0 + κ} 

i
= w0,j exp{( β1 − 1)µ̄∗j } +

∗ )γ
(1
−
q
0
i ,j
i0 =1

(k+1)

= si,0 (x(k+1) ) , wi,0 xi,0

+

X

(k+1) 

wi,j 0 xi,j 0

Our goal in what follows is to prove that the composite
operator r(s(·)) is a contraction, as defined below. We utilise
the same distance metric as in [32]:

,

d(x, x̃) = max log
i,j

j0

(162)
(k+1)

y0,j

(k+1)

= s0,j (x(k+1) ) , (x0,j

1

) β −1 .

(163)

0
The shorthand
i0Prepresents the sum over the set i ∈
{1, . . . , n}, while
the same summation,
i0 6=i represents P
excluding the i-th element. Similarly,
P j 0 represents the sum
over the set j 0 ∈ {1, . . . , m}, while j 0 6=j represents the same
summation, excluding the j-th element. The structure of this
iterative method is illustrated in figure 11. Note that if γ = 1,
this reduces to the BP iteration of [32].
At this point, we have stated but not derived the iteration
(158)–(163). The validity of the expressions is established by
proving that the solution of the KKT conditions is a fixed point
of the iteration (in lemma 9), and then showing that repeated
application of the expressions yields a contraction, which is
guaranteed to converge to the unique fixed point.

P

∗
Lemma 9. Let (qi,j
, λ̄∗i , µ̄∗j ) be the solution of the KKT
∗
conditions in lemma 8, and let x∗i,j and yi,j
be the values
∗
∗
∗
calculated from (qi,j , λ̄i , µ̄j ) using (156)–(157). Then x∗i,j and
∗
yi,j
are a fixed point of r(·) and s(·) in (158)-(163).
Pn
Proof. Feasibility implies that i=0 qi,j = 1 ∀ j. Therefore:

w0,j exp{ β1 µ̄∗j } +

n
X
wi,j exp{λ̄∗i + µ̄∗j + κ}

(1 −

i=1

∗ )γ
qi,j

= 1,

(164)

where we define 0/0 = 1.
Definition 3. An operation g(x) is a contraction with respect
to d(·, ·) if there exists α ∈ [0, 1) such that for all x, x̃
d[g(x), g(x̃)] ≤ αd(x, x̃).

(169)

If the expression is satisfied for α = 1, then g(x) is a nonexpansion.
Lemma 10. Let g(·) be the operator taking the weighted
combination with non-negative weights wi,j ≥ 0:
X
gk,l (x) =
wi,j,k,l xi,j .
i,j

g(·) is non-expansive with respect to d.
Proof. Let L = exp{d(x, x̃)} < ∞ (otherwise there is nothing
to prove), so that L1 xi,j ≤ x̃i,j ≤ Lxi,j . Then
X
X
gk,l (x̃) =
wi,j,k,l x̃i,j ≤ L
wi,j,k,l xi,j = Lgk,l (x),
i,j

i,j

(170)
and
gk,l (x̃) =

X
i,j

wi,j,k,l x̃i,j ≥

1
L

X

wi,j,k,l xi,j =

1
L gk,l (x).

i,j

(171)

or
"
exp{µ̄∗j }

xi,j
,
x̃i,j

=

w0,j exp{( β1

−

1)µ̄∗j }

+

n
X
wi,j exp{λ̄∗ + κ}

#−1

i

∗ )γ
(1 − qi,j
(165)
Equating terms, this proves (160), and shows that the first
factor in (158) is exp{(1 − γ)µ̄∗j }. To prove (158), note that
i=1

1
∗ =
1 − qi,j
1−

1
∗
wi,j exp{λ̄∗
i +µ̄j +κ}
∗ )γ
(1−qi,j

,

(166)

Lemma 11. Let f (·) be formed from two operators g(·) and
h(·) as
fi,j (x) = gi,j (x)ρg hi,j (x)ρh .
Suppose that g(·) and h(·) are contractions or non-expansions
with coefficients αg and αh . If αf = αg |ρg | + αh |ρh | < 1,
then f (·) is a contraction with respect to d. If αf = 1 then
f (·) is a non-expansion.

24

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

(k−1)

∗
maxi,j [yi,j /yi,j
] = L for some L with 1 < L < ∞ (if
L = 1 then convergence has occurred, and L = ∞ will only
d[f (x), f (x̃)]
∗
occur if yi,j
= 0, which contradicts lemma 7). The result is
fi,j (x)
similar
to
that
obtained by changing the distance to Hilbert’s
= max log
(172)
i,j
fi,j (x̃)
projective metric (e.g., [66]). We emphasise that this rescaling
does not need to be performed in the online calculation; rather
gi,j (x)ρg hi,j (x)ρh
(173) we are exploiting an equivalence to aid in proving convergence.
= max log
i,j
gi,j (x̃)ρg hi,j (x̃)ρh
Lemmas 13 and 14 establish an induction which shows
gi,j (x)
hi,j (x)
≤ |ρg | max log
+ |ρh | max log
(174) that after n steps, we are guaranteed to have reduced
i,j
i,j
gi,j (x̃)
hi,j (x̃)
(k+n)
∗
maxi,j [yi,j /yi,j
]. The induction commences with a single
≤ (|ρg |αg + |ρh |αh )d(x, x̃).
(175)
(k−1)
∗
edge with yi,j /yi,j
= 1, setting T (k−1) = {(i, j)}, and
v (k−1) = 1. As the induction proceeds, the set T (k) (or,
alternately, S (k) ) represents the edges for which improvement
The following results immediately from lemmas 10 and 11.
in the bound L is guaranteed, and v (k) < L (or, alternately,
Corollary 1. The operators r(·) and s(·) defined in (158)- u(k) ) represents the amount of improvement that is guaranteed.
The induction proceeds by alternately visiting the left-hand
(162) are non-expansions.
vertices and right-hand vertices (e.g., in figure 2), at each stage
Lemma 12. If d(x, x̃) ≤ log L̄ < ∞, then the operator s(·) adding to the set S (k) edges (i, j) for which w > 0, and an
i,j
is a contraction in cases 2 and 3 with a contraction factor edge that is incident on the vertex j is in T (k) (or, alternatively,
dependent on L̄ and wi,j .
adding to T (k) edges that could be traversed by starting from
a vertex i represented by an edge in S (k) ).
Proof. Lemma 2 in [32] shows that the update:
w
13. At iteration (k − 1), suppose that 1 ≤
P i,j
(176) Lemma
yi,j =
(k−1)
(k−1)
∗
∗
1 + j 0 6=j wi,j 0 xi,j 0
yi,j /yi,j
≤ L < ∞ ∀ i, j, and yi,j /yi,j
≤ v (k−1) <

Proof.

(k)

is a contraction. The proof of lemma 2 in [32] may be trivially
modified to show that the updates:
c
P 1
,
(177)
c2 + j 0 6=j wi,j 0 xi,j 0
c
P 1
(178)
c2 + j 0 wi,j 0 xi,j 0

L ∀ (i, j) ∈ T (k−1) . Then 1/L ≤ xi,j /x∗i,j ≤ 1 ∀ (i, j), and
(k)
1/u(k) ≤ xi,j /x∗i,j ≤ 1 ∀ (i, j) ∈ S (k) , where

are also contractions for any c1 > 0, c2 > 0. In cases 2 and
3, α = 1, so xi,0 = 1. Therefore, these results combined with
lemma 11, show that (161) and (162) are contractions. Since
| β1 − 1| < 1, (163) is also a contraction.

and

This is adequate to prove convergence in cases 2 and 3: s(·)
is a contraction, and r(·) is a non-expansion, so the composite
operator is a contraction. Combined with lemma 8, this shows
that the iteration converges to the unique solution of the KKT
conditions. The fact that the contraction factor in lemma 12
depends on an upper bound on the distance L̄ is not of concern;
since the combined operation is a contraction, the contraction
factor for the upper bound L̄ that we begin with will apply
throughout.
The final step is to prove convergence in case 1. For this,
we prove that n successive iterations of applying operators r(·)
and s(·) collectively form a contraction. The proof is based on
[35], but adapts it to address γ < 1, and to admit cases where
some edges have wi,j = 0.
As discussed in lemma 8, the solution in case 1 is not
changed by adding any constant c to λ̄i ∀ i and subtracting
it from µ̄j ∀ j; this is clear from (151), and was termed
message gauge invariance in [35, remark 30]. Incorporating
any such constant simply offsets all future iterations by the
value, having no impact on the qi,j iterates produced. Thus,
for the purpose of proving convergence, when analysing r(·),
(k−1)
∗
∗
we scale yi,j
such that mini,j [yi,j /yi,j
] = 1, and we denote

S (k) =

(i, j) ∈ {1, . . . , n}2 ∃i0 s.t. (i0 , j) ∈ T (k−1) , wi,j > 0 ,
(179)
(k−1) (k−1)

u(k) = max[θj
j

v

(k−1)

+(1−θj

)L]1−γ Lγ < L, (180)

where
(k−1)

θj

X

=

∗
wi,j yi,j

.X

∗
wi,j yi,j
.

(181)

i

i|(i,j)∈T (k−1)

Proof. Consider the sum in the first factor in (158) (remembering that w0,j = 0):
X
X
(k)
(k−1)
∗
σj∗ =
wi,j yi,j
, σj =
wi,j yi,j ,
(182)
i

i

(k)

so that σj /σj∗ ≥ 1, and
P
P
(k)
∗
∗
v (k−1) i|(i,j)∈T (k−1) wi,j yi,j
+ L i|(i,j)∈T
σj
/ (k−1) wi,j yi,j
P
≤
∗
σj∗
i wi,j yi,j
(k−1) (k−1)

= θj

v

(k−1)

+ (1 − θj

)L.

While a similar analysis could be applied to the second factor
of the expression for r(y) (as in [35]), to prove convergence
in the case with γ < 1, it is adequate to simply bound it by:
P
(k−1)
0
i0 6=i wi ,j yi0 ,j
1≤ P
≤ L.
(183)
∗
0
i0 6=i wi ,j yi0 ,j
Substituting these bounds into (158) gives the desired result.

WILLIAMS AND LAU: MULTIPLE SCAN DATA ASSOCIATION BY CONVEX VARIATIONAL INFERENCE

(k)

Lemma 14. At iteration k, suppose that 0 < L1 ≤ xi,j /x∗i,j ≤
(k)
1 ∀ i, j, and 1/L < 1/u(k) ≤ xi,j /x∗i,j ∀ (i, j) ∈ S (k) . Then
(k)
(k)
∗
∗
1 ≤ yi,j /yi,j
≤ L ∀ (i, j) and 1 ≤ yi,j /yi,j
≤ v (k) ∀ (i, j) ∈
(k)
T , where

T (k) = (i, j) ∈ {1, . . . , n}2 ∃j 0 s.t. (i, j 0 ) ∈ S (k) , wi,j > 0
(184)
and
"
#1−γ
(k)
(k)
ωi
(1 − ωi )
1
1
1
= min (k) +
> ,
(185)
γ
(k)
i
L
L
L
v
u
where
(k)

ωi

X

=

wi,j x∗i,j

.X

wi,j x∗i,j .

(186)

j

j|(i,j)∈S (k)

Proof. Taking either the left or right derivative of α(L) in
(189):
L
1
∂v(L) · log
v(L) − L log v(L)
.
(190)
∂α(L) =
log2 L
Thus it suffices to show that (omitting the iteration index
superscript from v)

j

so that

(k)
τi /τi∗
1

P

ũj

∂+ v(L) ≥

(k+l−1) (k−l−1)

(L) = θj

v

(k+l−1)

(L) + (1 − θj
(k+l)

(k+l)

ũ(k+l) (L) = max ũj
j

≤ 1, and
j|(i,j)∈S (k)

)L.

By the second result in lemma 16, ũj
will satisfy the
property (191) for each j. The pointwise maximum in (180):

j

P
∗
wi,j x∗i,j + L1 j|(i,j)∈S
/ (k) wi,j yi,j
u(k)
P
≥
∗
j wi,j yi,j
1
(k) 1
(k)
+ (1 − ωi ) .
= ωi
L
u(k−1)
The second factor in (161) can be bounded by the expression
P
(k)
0
1
j 0 6=j wi,j xi,j 0
≤P
≤ 1.
(188)
∗
0
L
j 0 6=j wi,j xi,j 0
(k)
τi
τi∗

v(L) log v(L)
,
L log L

v(L) log v(L)
.
L log L
(191)
We prove by induction, showing that the property in (191) is
maintained through the recursion in (180) and (185). The base
case is established by noting that if v (k−1) = 1, then (191)
holds. Now suppose that the property (191) is held for some
iteration (k + l − 1) with upper bound v (k+l−1) (L). Then let
∂− v(L) ≥

(k+l)

Proof. Following similar steps to the proof of lemma 13, we
define:
X
X
(k)
(k)
τi∗ =
wi,j x∗i,j , τi =
wi,j xi,j ,
(187)

25

(L)

will introduce a finite number of points where ũ(k+l) (L) is
continuous but the derivative is discontinuous. However, the
one-sided derivative at any of these points will satisfy (191)
since each component in the pointwise maximum satisfied it.
Finally, the result of (180) is:
u(k+l) (L) = [ũ(k+l) (L)]1−γ Lγ .

The first result in lemma 16 shows that this will satisfy (191).
Now we need to prove that the other half-iteration (185)
maintains the property (191). First, note that if L̄ = 1/L and
Substituting these bounds into (161) gives the desired result. ū(k+l) = 1/u(k+l) then the final result in lemma 16 shows
that if u(k+l) (L) satisfies (191), then:
We employ lemmas 13 and 14 by setting v (k−1) = 1 (scaling

∗
yi,j
accordingly), and T (k−1) to contain the edge(s) with
(k−1)
∗
yi,j /yi,j
= 1. Iteratively applying the lemmas for n steps,
(k+n)
∗
we find that 1 ≤ yi,j /yi,j
≤ v (k+n) < L ∀ (i, j) since
(k+n)

we will have T
containing all edges (since the graph is
connected). To prove linear convergence, we first need to show
that the distance is reduced by at least a constant, or that
α(L) ,

log v (k+n) (L)
< 1.
log L

ū(k+l) (L̄) log ū(k+l) (L̄)
,
L̄ log L̄
ū(k+l) (L̄) log ū(k+l) (L̄)
∂+ ū(k+l) (L̄) ≥
.
L̄ log L̄

∂− ū(k+l) (L̄) ≥

(192)

Subsequently, if:
(k+l)

ṽi

(k+l) (k−l)

(L̄) = ωi

ū

(k+l)

(L̄) + (1 − ωi

)L̄,
(k+l)

(189)

then the second result in lemma 16 establishes that ṽi
will
satisfy the property (192) for each i. As with the pointwise
maximum in the the previous case, the pointwise minimum
(185):
(k+l)
ṽ (k+l) (L̄) = min ṽi
(L̄)

where v (k+n) (L) depends on L through the recursion in (180)
and (185). The inequality in (189) can be established simply
i
by commencing from v (k−1) = 1, and observing that if
(k+l−1)
(k+l)
(k+l)
v
< L then u
< L and v
< L for any l > 0. will introduce a finite number of points where ṽ (k+l) (L̄) is
This is not adequate to prove convergence; we further need continuous but the derivative is discontinuous. However, the
to show that the constant α(L) is non-decreasing in L. This one-sided derivative at any of these points will satisfy (192).
ensures that the initial contraction rate (for the first n iterations) The first result in lemma 16 shows that the composition:
applies, at least, in all subsequent n-step iteration blocks.
v̄ (k+l) (L̄) = [ṽ (k+l) (L̄)]1−γ L̄γ
Lemma 15. α(L) is continuous and non-decreasing in L, i.e.,
its left and right derivatives everywhere satisfy
will also satisfy (192). Finally, the result of (185) is:
∂− α(L) ≥ 0,

∂+ α(L) ≥ 0.

v (k+l) (L) = 1/v̄ (k+l) (1/L)

26

ACCEPTED FOR PUBLICATION, IEEE TRANSACTIONS ON SIGNAL PROCESSING

which will satisfy (192) by the final result in lemma 16.
du(x)
dx

u(x) log u(x)
.
x log x

Lemma 16. Suppose that
≥
Then if v(x)
is given by any of the following:
1) v(x) = u(x)1−γ xγ ,
2) v(x) = θu(x) + (1 − θ)x
v(x) log v(x)
then dv(x)
. Finally, if y = 1/x and
dx ≥
x log x
v(y) = 1/u(x) = 1/u(1/y),
then

dv(y)
dy

≥

where
φis (j) = ∆γs [1 + log(1 − qsi,j )] − ∆βs [1 + log qs0,j ]. (194)
Consider the KKT conditions relating to [qsi,j ], since all
modifications relate to these variables. The conditions in the
original problem are:
νsi,j + ρis + σsj + γs log(1 − qsi,j ) + γs = 0,
νs0,j + σsj + βs log qs0,j + βs = 0,

v(y) log v(y)
.
y log y

(195)
(196)

where νsi,j is the dual variable for the
in (55), ρis is
Pconstraint
ms i,j
the dual variable for the constraint j=0 qs = 1, and σsj is
the dual variable for (57). For the modified problem, the same
two KKT conditions are:

Proof. For the first case:

dv(x)
du(x)
= (1 − γ)
u(x)−γ xγ + γu(x)1−γ xγ−1
dx
dx
u(x) log u(x)
≥ (1 − γ)
u(x)−γ xγ + γu(x)1−γ xγ−1
ν̄si,j + ρ̄is + σ̄sj + γ̄s log(1 − q̄si,j ) + γ¯s − φis (j) = 0, (197)
x log x
ν̄s0,j + σ̄sj + β̄s log q̄s0,j + β̄s = 0.
(198)
u(x)1−γ xγ [(1 − γ) log u(x) + γ log x]
=
x log x
Substituting in (194) and expanding γ̄s and β̄s , we find:
v(x) log v(x)
.
=
ν̄si,j + ρ̄is + σ̄sj + [γs + ∆γs ] log(1 − q̄si,j ) + γs + ∆γs
x log x
− ∆γs [1 + log(1 − qsi,j )] + ∆βs [1 + log qs0,j ] = 0, (199)
For the second case:
dv(x)
du(x)
Subsequently, by setting q̄si,j = qsi,j , q̄s0,j = qs0,j , ν̄si,j = νsi,j ,
=θ
+ (1 − θ)
dx
dx
ρ̄is = ρis and
u(x) log u(x)
+ (1 − θ)
≥θ
σ̄sj = σsj − ∆βs [1 + log qs0,j ],
(200)
x log x
θu(x) log u(x) + (1 − θ)x log x
we find a primal-dual solution (with identical primal values
=
x log x
[qsi,j ]) that satisfies the KKT conditions for the modified
v(x) log v(x)
problem, providing a certificate of optimality.
,
≥
x log x
where the final inequality is the result of convexity of x log x.
For the final result, let f (x) = 1/u(x) and g(y) = 1/y and
apply the chain rule:
dv(y)
= f 0 (g(y)) × g 0 (y)
dy
u0 (1/y)
1
=−
×− 2
2
u(1/y)
y
1
u(1/y) log u(1/y)
≥
×
u(1/y)2 y 2
(1/y) log(1/y)
[1/u(1/y)] log[1/u(1/y)]
v(y) log v(y)
=
=
.
y log y
y log y

A PPENDIX E
P ROOF OF SEQUENTIAL MODIFICATION
This section proves theorem 3, i.e., that the solution of the
problem in (58) is the same as the solution of the modified
problem of the same form, changing γs to γ̄s = γs + ∆γs , βs
to β̄s = βs + ∆βs , and ψsi (xi , ais ) as described in (86).
Proof of theorem 3: Let FBγ,β be the original problem
(in (58)) and F̄Bγ̄,β̄ be the modified problem. Note that the
modifying term in (86) depends only on ais , so we can
equivalently implement the modification by retaining the
unmodified ψsi (xi , ais ) and incorporating an additive term
− E[φis (ais )] = −

ms
X
j=0

qsi,j φis (j),

(193)

Jason L. Williams (S’01–M’07–SM’16) received
degrees of BE(Electronics)/BInfTech from Queensland University of Technology in 1999, MSEE from
the United States Air Force Institute of Technology
in 2003, and PhD from Massachusetts Institute of
Technology in 2007. He worked for several years as
an engineering officer in the Royal Australian Air
Force, before joining Australia’s Defence Science
and Technology Group in 2007. He is also an
adjunct associate professor at Queensland University
of Technology. His research interests include target
tracking, sensor resource management, Markov random fields and convex
optimisation.

Roslyn A. Lau (S’14) received the degrees of BE
(Computer Systems)/BMa&CS (Statistics) in 2005,
and MS (Signal Processing) in 2009, all from the
University of Adelaide, Adelaide, Australia. She is
currently a PhD candidate at the Australian National
University. She is also a scientist at the Defence
Science and Technology Group, Australia. Her research interests include target tracking, probabilistic
graphical models, and variational inference.

