Markov Chain Monte Carlo Data Association for
General Multiple Target Tracking Problems
Songhwai Oh, Stuart Russell, Shankar Sastry
Abstract—In this paper, we consider the general multiple
target tracking problem in which an unknown number of
targets appears and disappears at random times and the goal
is to find the tracks of targets from noisy observations. We
propose an efficient real-time algorithm that solves the data
association problem and is capable of initiating and terminat-
ing a varying number of tracks. We take the data-oriented,
combinatorial optimization approach to the data association
problem but avoid the enumeration of tracks by applying a
sampling method called Markov chain Monte Carlo (MCMC).
The MCMC data association algorithm can be considered as a
deferred logic since its decision about forming a track is based
on the current and past observations. But, at the same time, it
can be considered as an approximation to the optimal Bayesian
filter. The algorithm shows remarkable performance compared
to the greedy algorithm and the multiple hypothesis tracker
(MHT) under the extreme conditions, such as a large number
of targets in a dense environment, low detection probabilities,
and a large number of false alarms.
I. INTRODUCTION
The multiple target tracking plays an important role in
many areas of engineering such as surveillance, computer
vision, and signal processing [1], [4]. Under the most
general setup, a varying number of indistinguishable targets
is moving around in a region with continuous motions and
the positions of moving targets are sampled at random
intervals. The measurements about the positions are noisy,
with detection probability less than one, and there is a
noise background of spurious position reports, i.e., false
alarms. Targets arise at random in space and time. Each
target persists independently for a random length of time
and ceases to exist. A track of a target is defined as a path
in space-time traveled by the target. The essence of the
multiple target tracking problem is to find tracks from the
noisy observations and it requires solutions to both data
association and state estimation problems [16].
The data association problem in multiple target tracking
is described as a problem of finding a partition of observa-
tions such that each element of a partition is a collection
of observations generated by a single target or clutter
[16]. However, due to the noises in state transitions and
observations, we cannot expect to find the exact solution.
This data-oriented view of data association has been applied
and extended by many authors [9], [17], [15], [8], [5],
[14]. The most successful multiple target tracking algorithm
This work was supported by DARPA F30602-00-2-0538 and ARO
MURI DAAD 19-OZ-1-0383.
The authors are with the Department of Electrical Engineering and
Computer Sciences, University of California, Berkeley, CA 94720,
{sho,russell,sastry}@eecs.berkeley.edu.
based on this view is the multiple hypothesis tracker (MHT)
[15]. In MHT, each hypothesis associates past observations
with a target and, as a new set of observations arrives, a new
set of hypotheses is formed from the previous hypotheses.
Each hypothesis is scored by its posterior and the algorithm
returns a hypothesis with the highest score as a solution.
MHT is categorized as a deferred logic [14] in which the
decision about forming a new track or removing an existing
track is delayed until enough observations are collected.
Hence, MHT is capable of initiating and terminating a
varying number of tracks and suitable for surveillance
applications in which an autonomous tracker is required.
However, the construction of new hypotheses requires an
enumeration of all possibilities and the size of hypotheses
grows exponentially. The initial implementation and later
extensions proposed several heuristics, such as pruning,
gating, clustering and N-scan-back logic, to reduce the
complexity of the problem [15], [8]. However, the heuristics
are used at the expense of the optimality and the algorithm
can still suffer in a dense environment. Furthermore, the
running time at each step of the algorithm cannot be
bounded easily, making it difficult to be deployed in a
real-time surveillance system. As a method of pruning, an
efficient method of finding k-best hypotheses based on the
algorithm by Murty [10] is developed in [5].
A different approach to the data association problem is
the joint probabilistic data association filter (JPDAF) [1].
JPDAF is a suboptimal single-stage approximation to the
optimal Bayesian filter. Given a fixed number of targets,
JPDAF enumerates all possible associations between the
latest set of observations and the known tracks and clutter
and computes each association weight. For each association,
the conditional expectation of the state of a target is
estimated by a filtering algorithm. Then, the state of a target
is estimated by summing over the conditional expectations
weighted by association weights. JPDAF is a sequential
tracker in which the associations between the known targets
and the latest observations are made sequentially and the
associations made in the past are not reversible [14]. Since
only the current set of observations is considered, JPDAF
cannot initiate or terminate tracks. Also JPDAF assumes
a fixed number of targets and requires a good initial state
for each target. There are restricted extensions to JPDAF to
allow the formation of a new track (see [4] and references
within). Other multiple target tracking algorithms include
the multisensor multitarget mixture reduction (MTMR) [12]
and the probabilistic multi-hypothesis tracker (PMHT) [18].
But they also assume a fixed number of targets and cannot
initiate or terminate tracks. Recently, a Bayesian model-
based approach for tracking an unknown number of targets
which can initiate and terminate tracks is presented in [11].
The sequential trackers are more efficient than deferred
logic trackers such as MHT but they are prone to make
erroneous associations [14]. In addition, it is inferred that
the exact calculation at each stage is NP-hard [3] since the
related problem of finding the permanent of a 0-1 matrix
is #P-complete [19]. In [6], a single-stage data association
problem is considered and a leave-one-out heuristic is
developed to avoid the enumeration of all possible asso-
ciations. Later, the approach is extended to a multi-stage
data association problem using Markov chain Monte Carlo
(MCMC) [13].
The data association problem of multiple target tracking
formulated under the data-oriented view is also known
to be NP-hard [14]. Hence, we cannot expect to find an
optimal solution in polynomial time unless P = N P .
An optimization approach to data association has been
applied as a 0-1 integer programming problem [9] and as
a multidimensional assignment problem [14]. In both cases
one needs to find a feasible set of tracks from all possible
tracks to prevent the exponential explosion and compute
the cost of each feasible track, such as the negative log
likelihood. Then the optimization routine finds a subset
from the feasible tracks such that the combined costs are
minimized while satisfying the constraints, i.e., each track
has at most one observation at each time and no two tracks
share the same observation. The gating method similar to
the ones described in [17], [15] is used to find a feasible
set of tracks. However, in a dense environment, the size of
the feasible tracks can be very large and the complexity
of the optimization routine increases dramatically, since the
number of parameters in the optimization routine depends
on the number of feasible tracks.
The main contribution of this paper is the development
of an efficient real-time algorithm that solves the data
association problem and is capable of initiating and ter-
minating a varying number of tracks. We take the data-
oriented, combinatorial optimization approach to the data
association problem but avoid the enumeration of tracks by
applying a sampling method called Markov chain Monte
Carlo (MCMC). The immediate benefit of using MCMC is
the low memory requirement. The MCMC data association
(MCMCDA) algorithm can be considered as a deferred
logic since its decision about forming a track is based on
the current and past observations. But, at the same time,
it can be considered as an approximation to the optimal
Bayesian filter if it is used to approximate the association
probabilities or expectations such as the average link travel
time as done in [13]. So, from the Bayesian point of view,
the algorithm can be considered as a generalization of
[13] to handle an unknown number of objects, missing
observations and false alarms. MCMCDA shows remarkable
performance compared to the greedy algorithm and MHT
under the extreme conditions such as a large number of
targets in a dense environment, low detection probabilities,
and a large number of false alarms.
The remainder of this paper is structured as follows. We
formally state the (discrete-time) general multiple target
tracking problem in Section II. In Section III, we present
a general purpose MCMCDA algorithm for multiple target
tracking. The algorithm is applied in simulation to extreme
situations and its performance is compared with the greedy
algorithm and MHT in Section IV.
II. GENERAL MULTIPLE TARGET TRACKING
A. Problem
Let T ∈ Z
+
be the duration of surveillance. Let K
be the unknown number of objects moving around the
surveillance region R for some duration [t
k
i
, t
k
f
] ⊂ [1, T ]
for k = 1, . . . , K. Let V be the volume of R. Each object
arises at a random position in R at t
k
i
, moves independently
around R until t
k
f
and disappears. At each time, an existing
target persists with probability (1 − p
z
) and disppears with
probability p
z
. The number of objects arising at each time
over R has a Poisson distribution with a parameter (λ
b
V )
where λ
b
is the birth rate of new objects per unit time,
per unit volume. The initial position of a new object is
uniformly distributed over R.
Let F
k
: R
d
→ R
d
be the discrete-time dynamics of the
object k, where d is the dimension of the state variable,
and let x
k
t
∈ R
d
be the state of the object k at time t for
k = 1, 2, . . . , K. The object k moves according to
x
k
t+1
= F
k
(x
k
t
) + w
k
t
for t = t
k
i
, . . . , t
k
f
− 1,
where w
k
t
∈ R
d
are white noise processes. The noisy
observation about the state of the object is measured with
the detection probability p
d
which is less than unity. There
are also false alarms and the number of false alarms has a
Poisson distribution with a parameter (λ
f
V ) where λ
f
is the
false alarm rate per unit time, per unit volume. Let n
t
be
the number of observations at time t which includes both
noisy observations and false alarms. Let y
j
t
∈ R
m
be the
j-th observation at time t for j = 1, . . . , n
t
, where m is
the dimensionality of each observation vector. Each object
generates a unique observation at each sampling time if it
is detected. Let H
j
: R
d
→ R
m
be the observation model.
Then the observations are generated as follows:
y
j
t
=
H
j
(x
k
t
) + v
j
t
if j-th observation is from x
k
t
u
t
otherwise,
where v
j
t
∈ R
m
are white noise processes and u
t
∼
Unif(R) is a random process for false alarms. Notice that,
with probability 1 − p
d
, the object is not detected and we
call this a missing observation. We assume that targets are
indistinguishable in this paper. But, if observations include
target type or attribute information, the state variable can
be extended to include target type information.
Under the data-oriented approach, the multiple target
tracking problem is to partition the observations such that
the posterior is maximized, i.e., the maximum a posteriori
(MAP) estimate. Under the Bayesian approach, if we are
given a function defined on Ω, the collection of all partitions
of observations (see below for its definition), we seek the
expected value of the function given the observations. The
MAP estimate found under the data-oriented approach may
not be robust in the Bayesian sense. But it is sometimes
more convenient when estimating parameters whose dimen-
sion is dependent on the number of tracks, such as the states
of targets. Since the size of x
t
= (x
1
t
T
, . . . , x
K
t
T
)
T
depends
on the number of tracks K, the estimation of x
t
without
fixing the number of tracks is not meaningful. Hence, under
the Bayesian approach, if a single set of state estimation is
required, we might first estimate the most likely number of
targets and then estimate the expected values of states given
the estimated number of targets.
B. Probabilistic Model
Let us first specify the dynamic and measurement models.
Here we use the usual linear system model but the method
can be easily extended to non-linear models coupled with a
non-linear regression algorithm. If an object is observed k
times at t
1
, t
2
, . . . , t
k
, its dynamic and measurement models
can be expressed as:
x
t
i+1
= A(t
i+1
− t
i
)x
t
i
+ G(t
i+1
− t
i
)w
t
i
y
t
i
= Cx
t
i
+ v
t
i
for i = 1, . . . , k,
(1)
where w
t
i
and v
t
i
are white Gaussian noises with zero mean
and covariance Q and R, respectively. A(·), G(·), and C are
matrices with appropriate sizes. The entries of the matrix
A(t
i+1
−t
i
) and G(t
i+1
−t
i
) are determined by the sampling
interval t
i+1
− t
i
for each i. For clarity, the subsequence
notation for the time index is suppressed for now. Let ¯x
t
be the expected value of x
t
given y
1
, . . . , y
t−1
;
¯
P
t
be the
covariance of x
t
given y
1
, . . . , y
t−1
; ˆx
t
be the expected
value of x
t
given y
1
, . . . , y
t
; and
ˆ
P
t
be the covariance of
x
t
given y
1
, . . . , y
t
.
Let y
t
= {y
j
t
: j = 1, . . . , n
t
} and Y =
S
t∈{1,,...,T }
y
t
.
Let Ω be a collection of partitions of Y such that, for ω ∈ Ω,
1) ω = {τ
0
, τ
1
, . . . , τ
K
};
2)
S
K
k=0
τ
k
= Y and τ
i
∩ τ
j
= ∅ for i 6= j;
3) τ
0
is a set of false alarms;
4) |τ
k
∩ y
t
| ≤ 1 for k = 1, . . . , K and t = 1, . . . , T ; and
5) |τ
k
| > 1 for k = 1, . . . , K.
Here, K is the number of tracks for the given partition
ω ∈ Ω. We call τ
k
a track when there is no confusion
although the actual track is the set of estimated states
from the observations τ
k
. However, we assume there is a
deterministic function that returns a set of estimated states
given a set of observations, so no distinction is required. We
denote by τ
k
(t) the observation at time t that is assigned to
the track τ
k
. Notice that τ
k
(t) can be empty if it is a missing
observation. The fourth requirement says that a track can
have at most one observation at each time, but, in the case
of multiple sensors, we can easily relax this requirement to
allow multiple observations per track. A track is assumed to
contain at least two observations since we cannot distinguish
a track with a single observation from a false alarm. Once
a partition ω ∈ Ω is chosen, the tracks τ
1
, . . . , τ
K
∈ ω and
a set of false alarms τ
0
∈ ω are completely determined.
Hence, for each track, we can estimate the states of an
object independently since each object moves independently
from the other objects. For each track τ ∈ ω, we apply the
Kalman filter to estimate the states ¯x
t
(τ) and covariances
B
t
(τ), where B
t
(τ) = C
¯
P
t
(τ)C
T
+ R is the conditional
observation covariance at time t for the track τ .
Let e
t
be the number of targets from time t −1 and a
t
be
the number of new targets at time t. Let z
t
be the number of
targets terminated at time t and c
t
= e
t
− z
t
. Let d
t
be the
number of detections at time t and u
t
= e
t
−z
t
+a
t
−d
t
be
the number of undetected targets. Finally, let f
t
= n
t
− d
t
be the number of false alarms. It can be shown that the
posterior of ω is:
P (ω|Y ) =
1
Z
Q
T
t=1
p
z
t
z
(1 − p
z
)
c
t
p
d
t
d
(1 − p
d
)
u
t
λ
a
t
b
λ
f
t
f
×
Q
τ ∈ω\{τ
0
}
Q
|τ |−1
i=1
N (τ(t
i+1
)|¯x
t
i+1
(τ), B
t
i+1
(τ)),
(2)
where Z is a normalizing constant and N (·|µ, Σ) is the
Gaussian density function with mean µ and covariance
matrix Σ. Now under the data-oriented, combinatorial op-
timization approach, our goal is to find a partition of
observations such that P (ω|Y ) is maximized.
III. MCMC DATA ASSOCIATION ALGORITHM
In this section, we develop an MCMC sampler to solve
the multiple target tracking problem. Solving complex prob-
lems by sampling methods such as Markov chain Monte
Carlo (MCMC) has become more tractable, due to the
increased computational power. MCMC-based algorithms
play a significant role in many fields such as physics,
statistics, economics, and engineering [2]. In some cases,
MCMC is the only known general algorithm which finds
a good approximate solution to a complex problem in
polynomial time [7]. MCMC techniques have been applied
to the complex probability distribution integration problems,
counting problems such as #P-complete problems, and
combinatorial optimization problems [7], [2]. The MCMC
approach applied to the combinatorial optimization prob-
lems is generally known as simulated annealing.
The set Ω becomes a state space of the MCMC sampler
and we sample from Ω such that its stationary distribution
is P (ω|Y ). If we are at state ω ∈ Ω, we propose ω
0
∈ Ω
following the proposal distribution q(ω, ω
0
). The move is
accepted with an acceptance probability A(ω, ω
0
) where
A(ω, ω
0
) = min
1,
P (ω
0
|Y )q(ω
0
, ω)
P (ω|Y )q(ω, ω
0
)
, (3)
otherwise the sampler stays at ω, so that the detailed balance
is satisfied. If we make sure that the chain is irreducible
and aperiodic, then the chain converges to its stationary
distribution. The sampler consists of five types of moves.
They are
1) birth/death move pair;
2) split/merge move pair;
3) extension/reduction move pair;
4) track update move; and
5) track switch move.
We index each move by an integer such that m = 1 for a
birth move, m = 2 for a death move and so on. The move
m is chosen randomly from the distribution ξ
K
(m) where
K is the number of tracks of the current partition ω. When
there is no track, we can only propose a birth move, so we
set ξ
0
(m = 1) = 1 and 0 for all other moves. When there
is only a single target, we cannot propose a merge or track
switch move, so ξ
1
(m = 4) = ξ
1
(m = 8) = 0. For other
values of K and m, we assume ξ
K
(m) > 0. The MCMC
data association (MCMCDA) algorithm is described in
Algorithm 1. The inputs are the set of all observations
Y , the number of samples n
mc
, and the initial state ω
init
.
At each step of the algorithm, ω is the current state of
the Markov chain. The acceptance probability A(ω, ω
0
) is
defined in (3) where the posterior (2) is used.
Algorithm 1 (MCMC Data Association (MCMCDA)):
Input: Y, n
mc
, ω
init
Output: ˆω = arg max p(ω (n)|Y )
ω ← ω
init
for n = 1 to n
mc
sample m from ξ
K
(·)
propose ω
0
based on m and ω (described below)
sample U from Unif[0, 1]
if U < A(ω, ω
0
)
ω ← ω
0
end
ω(n) ← ω
end
In Algorithm 1, we use MCMC to find a solution to a
combinatorial optimization problem. So it can be considered
as simulated annealing at a constant temperature. No burn-in
samples are used since we are simply looking for a partition
which maximizes the posterior. In addition, the memory
requirement of the algorithm is at its bare minimum. Instead
of keeping all {ω(n)}
n
mc
n=1
, we can simply keep the partition
with the maximum posterior. If the algorithm is used
to estimate E
P (ω|Y )
f(ω) for some bounded function f,
we will need burn-in samples and need to maintain the
sufficient statistics for the desired expectation.
In order to make the algorithm more efficient, we make
two additional assumptions: (1) the maximal directional
speed of any target in R is less than ¯v; and (2) the number of
consecutive missing observations of any track is less than
¯
d.
The first assumption is reasonable in a surveillance scenario
since, in many cases, the maximal speed of a vehicle is
generally known based on its type and terrain conditions.
The second assumption is a user defined parameter and it
can be used as one of the criteria to distinguish an event of
a new object’s appearance from an event of a continuation
of an existing object. We will now assume that these two
new conditions are added to the definition of Ω so each
element ω ∈ Ω satisfies these two additional assumptions.
We now introduce a data structure which is used to
propose a new partition ω
0
in Algorithm 1. We define a
neighborhood tree of observations as
L
d
(y
j
t
) = {y
k
t+d
∈ y
t+d
: ky
j
t
− y
k
t+d
k ≤ d · ¯v}
for d = 1, . . . ,
¯
d, j = 1, . . . , n
t
and t = 1, . . . , T − 1. Here
k · k is the usual Euclidean distance. This neighborhood
tree groups temporally separated observations based on their
distances. The parameter d allows missing observations.
The use of this neighborhood tree makes the algorithm
more scalable since distant observations will be considered
separately and makes the computations of the proposal dis-
tribution easier. It is similar to the clustering technique used
in MHT but L
d
(·) is fixed for a given set of observations
Y .
We now describe each move of the sampler in detail.
First, let ζ(d) be a distribution of a random variable d taking
values from {1, 2, . . . ,
¯
d}. We assume the current state of
the chain is ω = ω
0
∪ ω
1
∈ Ω, where ω
0
= {τ
0
} and
ω
1
= {τ
1
, . . . , τ
K
}. The proposed partition is denoted by
ω
0
= ω
00
∪ ω
01
∈ Ω. Note the abuse of notation below with
indexing of time, i.e., when we say τ(t
i
), t
i
means the time
at which a target corresponding to the track τ is observed
i times.
A. Birth and Death Moves
For a birth move, we increase the number of tracks from
K to K
0
= K +1 and select t
1
uniformly at random (u.a.r.)
from {1, . . . , T − 1} as an appearance time of a new track.
Let τ
K
0
be the track of this new object. Then we choose
d
1
from the distribution ζ. Let L
1
d
1
= {y
j
t
1
: L
d
1
(y
j
t
1
) 6=
∅, y
j
t
1
6∈ τ
k
(t
1
), j = 1, . . . , n
t
1
, k = 1, . . . , K}. L
1
d
1
is
a set of observations at t
1
such that, for any y ∈ L
1
d
1
,
y does not belong to other tracks and y has at least one
descendant in L
d
1
(y). We choose τ
K
0
(t
1
) u.a.r. from L
1
d
1
.
If L
1
d
1
is empty, the move is rejected since the move is
not reversible. Once the initial observation is chosen, we
then choose the subsequent observations for the track τ
K
0
.
For i = 2, 3, . . ., we choose d
i
from ζ and choose τ
K
0
(t
i
)
u.a.r. from L
d
i
(τ
K
0
(t
i−1
))\ {τ
k
(t
i−1
+d
i
) : k = 1, . . . , K}
unless this set is empty. But, for i = 3, 4, . . ., the processe
of adding observations to τ
K
0
terminates with probability
γ, where 0 < γ < 1. If |τ
K
0
| ≤ 1, the move is rejected. We
then propose this modified partition where ω
01
= ω
1
∪{τ
K
0
}
and ω
00
= {τ
0
\τ
K
0
}. For a death move, we simply choose k
u.a.r. from {1, . . . , K} and delete the k-th track and propose
a new partition where ω
01
= ω
1
\{τ
k
} and ω
00
= {τ
0
∪τ
k
}.
B. Split and Merge Moves
For a split move, we select τ
s
(t
r
) u.a.r. from {τ
k
(t
i
) :
|τ
k
| ≥ 4, i = 2, . . . , |τ
k
| − 2, k = 1, . . . , K}. Then we split
the track τ
s
into τ
s
1
and τ
s
2
such that τ
s
1
= {τ
s
(t
i
) :
i = 1, . . . , r} and τ
s
2
= {τ
s
(t
i
) : i = r + 1, . . . , |τ
s
|}.
The modified track partition becomes ω
01
= (ω
1
\ {τ
s
}) ∪
{τ
s
1
} ∪ {τ
s
2
} and the false alarm partition ω
00
is updated
accordingly. For a merge move, we consider the set
M = {(τ
k
1
(t
f
), τ
k
2
(t
1
)) : τ
k
2
(t
1
) ∈ L
t
1
−t
f
(τ
k
1
(t
f
)),
f = |τ
k
1
| for k
1
6= k
2
, 1 ≤ k
1
, k
2
≤ K}.
We select a pair (τ
s
1
(t
f
), τ
s
2
(t
1
)) u.a.r. from M . The tracks
are combined into a single track τ
s
= τ
s
1
∪ τ
s
2
. Then we
propose a new partition where ω
01
= (ω
1
\({τ
s
1
}∪{τ
s
2
}))∪
{τ
s
} and ω
00
with appropriate rearrangements.
C. Extension and Reduction Moves
In a track extension move, we select a track τ u.a.r. from
K available tracks in ω. We reassign observations for τ after
the disappearance time t
|τ |
as done in the track birth move.
For a track reduction move, we select a track τ u.a.r. from
K available tracks in ω and r u.a.r. from {2, . . . , |τ| − 1}.
We shorten the track τ to {τ (t
1
), . . . , τ (t
r
)} by removing
the observations assigned to τ after the time t
r+1
.
D. Track Update Move
In a track update move, we select a track τ u.a.r. from
K available tracks in ω. Then we pick r u.a.r. from
{1, 2, . . . , |τ|} and reassign observations for τ after the time
t
r
as done in the track birth move.
E. Track Switch Move
For a track switch move, we select a pair of observa-
tions (τ
k
1
(t
p
), τ
k
2
(t
q
)) from two different tracks such that,
τ
k
1
(t
p+1
) ∈ L
d
(τ
k
2
(t
q
)) and τ
k
2
(t
q+1
) ∈ L
d
0
(τ
k
1
(t
p
)),
where d = t
p+1
− t
q
, d
0
= t
q+1
− t
p
and 0 < d, d
0
≤
¯
d.
Then we let
τ
k
1
= {τ
k
1
(t
1
), . . . , τ
k
1
(t
p
), τ
k
2
(t
q+1
), . . . , τ
k
2
(t
|τ
k
2
|
)}
τ
k
2
= {τ
k
2
(t
1
), . . . , τ
k
2
(t
q
), τ
k
1
(t
p+1
), . . . , τ
k
1
(t
|τ
k
1
|
)}.
Theorem 1: Suppose that 0 < p
z
, p
d
< 1 and λ
b
, λ
f
> 0.
If ζ(d) > 0, for all d ∈ {1, . . . ,
¯
d}, then the Markov chain
designed by Algorithm 1 is irreducible.
Proof: The birth and death moves are sufficient to
illustrate the irreducibility of the chain. Since 0 < p
z
, p
d
< 1
and λ
b
, λ
f
> 0, P (ω|Y ) > 0 for all ω ∈ Ω. Take an arbitrary
partition ω ∈ Ω, say ω = {τ
0
, τ
1
, . . . , τ
K
}. Now consider
the partition ω
0
∈ Ω, such that ω
0
= {τ
0
0
}, i.e., ω
0
assigns all
observations as false alarms. Since ω is arbitrary, the chain
is irreducible if the chain can move from ω
0
to ω and from ω
to ω
0
. For the move from ω
0
to ω, consider K consecutive
birth moves: ω
0
= ω
0
, ω
1
= {{τ
0
0
\ τ
1
}, τ
1
}, . . . , ω
K
=
{{τ
0
0
\ {∪
K
i=1
τ
i
}}, τ
1
, . . . , τ
K
} = ω. Since ω ∈ Ω, all tracks
τ
k
are legal, i.e., τ
k
satisfies the constraints described in Sec-
tion II-B and, for i = 1, . . . , |τ
k
|−1, τ
k
(t
i+1
) ∈ L
d
(τ
k
(t
i
))
for 1 ≤ d = t
i+1
− t
i
≤
¯
d. Thus, ω
k
∈ Ω for all k. Because
ζ(d) > 0 and all tracks τ
k
are legal, the probability of
proposing τ
k
at ω
k−1
by the birth move is positive and
q(ω
k
, ω
k+1
) > 0. For the move from ω to ω
0
, consider
K consecutive death moves: ω
K
= ω, ω
K−1
, . . . , ω
0
= ω
0
.
The probability of removing the track τ
k
at ω
k
by the death
move is positive and q(ω
k+1
, ω
k
) > 0. Since P (ω
k
|Y ) > 0
for all k, the chain can move from ω
0
to ω and from ω to
ω
0
. Hence, the chain is irreducible.
The Markov chain designed by Algorithm 1 is irreducible
(Theorem 1) and aperiodic since there is always a positive
probability of staying at the current state in the track update
move. In addition, the transitions described in Algorithm 1
satisfy the detailed balance. Hence, by the ergodic theorem,
the chain converges to its stationary distribution. Notice that
the other moves are designed to improve the performance
of the algorithm.
Algorithm 2 (Greedy Multiple Target Tracking):
Input: Y, σ(threshold function)
Output: ω = ω
0
∪ ω
1
ω
1
← ∅
for t = 1 to T − 1
repeat
ω
t
← ∅
G ← {(y
i
t
, y
j
s
) : 1 ≤ s − t ≤
¯
d and y
i
t
, y
j
s
6∈ τ ∈ ω
1
}
foreach (y
i
t
, y
j
s
) in G
τ
t
← ∅
estimate an initial state from (y
i
t
, y
j
s
)
τ
t
(t
1
) ← y
i
t
and τ
t
(t
2
) ← y
j
s
r ← s
while r < T
for d = 1 to
¯
d
B ← {y ∈ L
d
(τ
t
(r)) : y 6∈ τ ∈ ω
1
}
if B 6= ∅
τ
t
(r + d) ← arg min
y∈B
ky − ¯x
r+d
(τ
t
)k
break
end
end
r ← r + d
end
if p(τ
t
|Y ) ≥ σ(|τ
t
|)
ω
t
← ω
t
∪ {τ
t
}
end
end
ω
1
← ω
1
∪ {arg max
τ ∈ω
t
p(τ |Y )}
until ω
t
6= ∅
end
ω
0
← {Y \ (∪
τ ∈ω
1
τ )}
IV. SIMULATION RESULTS
For the simulations we consider the surveillance over a
rectangular region on a plane, R = [0, L]×[0, L] ⊂ R
2
. The
state vector is x = [x, y, ˙x, ˙y]
T
where (x, y) is a position
on R along the usual x and y axes and ( ˙x, ˙y) is a velocity
vector. The linear system model (1) is used where δ is an
interval between observations and
A(δ) =
1 0 δ 0
0 1 0 δ
0 0 1 0
0 0 0 1
G(δ) =
δ
2
2
0
0
δ
2
2
δ 0
0 δ
C =
1 0
0 1
0 0
0 0
T
The covariance matrices are Q = diag(100, 100) and R =
diag(25, 25).
The complexity of the multiple target tracking problems
can be measured by several metrics: (1) the intensity of
the false alarm rate λ
f
; (2) the detection probability p
d
;
(3) the number of tracks K; and (4) the density of tracks.
