A Multi-Scan Labeled Random Finite Set Model for Multi-object State Estimation

State space models in which the system state is a finite set--called the multi-object state--have generated considerable interest in recent years. Smoothing for state space models provides better estimation performance than filtering by using the full posterior rather than the filtering density. In multi-object state estimation, the Bayes multi-object filtering recursion admits an analytic solution known as the Generalized Labeled Multi-Bernoulli (GLMB) filter. In this work, we extend the analytic GLMB recursion to propagate the multi-object posterior. We also propose an implementation of this so-called multi-scan GLMB posterior recursion using a similar approach to the GLMB filter implementation.


I. INTRODUCTION
In Bayesian estimation for state-space models, smoothing yields significantly better estimates than filtering by using the history of the states rather than the most recent state [18], [6], [10]. Conditional on the observation history, filtering only considers the current state via the filtering density, whereas smoothing considers the sequence of states up to the current time via the posterior density. Numerical methods for computing the filtering and posterior densities have a long history and is still an active area of research, see for example [1], [25], [9], [10]. Recursive computation of the posterior density is also known as smoothing-while-filtering [6].
A generalisation of state-space models that has attracted substantial interest in recent years is Mahler's Finite Set Statistics (FISST) framework for multi-object system [14], [15], [16], [17]. Instead of a vector, the state of a multi-object system at each time, called the multi-object state, is a finite set of vectors. Since its inception, a host of algorithms have been developed for multi-object state estimation [16], [17]. By incorporating labels (or identities), multi-object state estimation provides a state-space formulation of the multi-object tracking problem where the aim is to estimate the number of objects and their trajectories [5], 1 [2], [16]. Numerically, this problem is far more complex than standard state estimation due to additional challenges such as false measurements, misdetection and data association uncertainty.
In multi-object state estimation, the labeled multi-object filtering recursion admits an analytic solution known as the Generalized Labeled Multi-Bernoulli (GLMB) filter [28], [30]. Moreover, this recursion can be implemented with linear complexity in the number of measurements and quadratic in the number of hypothesized objects [31]. Since the filtering density only considers information on the current multi-object state, earlier estimates cannot be updated with current data. Consequently, apart from poorer performance compared to smoothing, an important drawback in a multi-object context is track fragmentation, where terminated trajectories are picked up again as new evidence from the data emerges.
In this paper, we extend the GLMB filtering recursion to a (labeled) multi-object posterior recursion.
Such posterior captures all information on the set of underlying trajectories and eliminates track fragmentation as well as improving general tracking performance. Specifically, by introducing the multiscan GLMB model, an analytic multi-object posterior recursion is derived. Interestingly, the multi-scan GLMB recursion takes on an even simpler and more intuitive form than the GLMB recursion. In implementation, however, the multi-scan GLMB recursion is far more challenging. Like the GLMB filter, the multi-scan GLMB filter needs to be truncated, and as shown in this article, truncation by retaining components with highest weights minimizes the L 1 truncation error. Unlike the GLMB filter, finding the significant components of a multi-scan GLMB filter is an NP-hard multi-dimensional assignment problem.
To solve this problem, we propose an extension of the Gibbs sampler for the 2-D assignment problem in [31] to higher dimensions. The resulting technique can be applied to compute the GLMB posterior off-line in one batch, or recursively as new observations arrive, thereby performing smoothing-while-filtering.
The remainder of this article is divided into 5 Sections. Section II summarizes relevant concepts in Bayesian multi-object state estimation and the GLMB filter. Section III introduces the multi-scan GLMB model and the multi-scan GLMB posterior recursion. Section IV presents an implementation of the multi-scan GLMB recursion using Gibbs sampling. Numerical studies are presented in Section V and conclusions are given in Section VI.

II. BACKGROUND
Following the convention in [28], the list of variables X m , X m+1 , ..., X n is abbreviated as X m:n , and the inner product f (x)g(x)dx is denoted by f, g . For a given set S, 1 S (·) denotes the indicator function of S, and F(S) denotes the class of finite subsets of S. For a finite set X, its cardinality (or number of elements) is denoted by |X|, and the product x∈X f (x), for some function f , is denoted by the multi-object exponential f X , with f ∅ = 1. In addition we use for a generalization of the Kroneker delta that takes arbitrary arguments.

A. Trajectories and Multi-object States
This subsection summarizes the representation of trajectories via labeled multi-object states.
At time k, an existing object is described by a vector x ∈ X and a unique label ℓ = (s, α), where s is the time of birth, and α is a unique index to distinguish objects born at the same time (see Fig. 1 in [30]). Let B s denote the label space for objects born at time s, then the label space for all objects up to time k (including those born prior to k) is given by the disjoint union A trajectory is a sequence of labeled states with a common label, at consecutive times [28], i.e. a trajectory with label ℓ = (s, α) and kinematic states x s , x s+1 , ..., x t ∈ X, is the sequence A labeled multi-object state at time i is a finite subset X of X×L i with distinct labels. More concisely, let L : X×L k → L k be the projection defined by L((x, ℓ)) = ℓ, then X has distinct labels if and only if the distinct label indicator ∆(X) δ |X| [|L(X)|] equals one. The labeled states, at time i, of a set S of trajectories (with distinct labels) is the labeled multi-object state denotes the labeled state of trajectory τ at time i.
Consider a sequence X j:k of labeled multi-object states in the interval {j : k}. Let x i , ℓ) denote the element of X i with label ℓ ∈ L(X i ). Then the trajectory in X j:k with label ℓ ∈ ∪ k i=j L(X i ) is the sequence of states with label ℓ: where is the start time of label ℓ in the interval {j : k}, and is the latest time in {s(ℓ) : k} such that label ℓ still exists.
The multi-object state sequence X j:k can thus be equivalently represented by the set of all such trajectories, i.e.
The left and right hand sides of (5) are simply different groupings of the labeled states on the interval {j : k}. The multi-object state sequence groups the labeled states according to time while the set of trajectories groups according to labels (see also figure 1 of [30]).
For the rest of the article, single-object states are represented by lowercase letters (e.g. x, x), while multi-object states are represented by uppercase letters (e.g. X, X), symbols for labeled states and their distributions are bolded to distinguish them from unlabeled ones (e.g. x, X, π, etc).
Given the observation history Z 1:k = (Z 1 , ..., Z k ), all information on the set of objects (and their trajectories) is captured in the multi-object posterior density, π 0:k (X 0:k ) π 0:k (X 0:k |Z 1:k ). Note that the dependence on Z 1:k is omitted for notational compactness. Similar to standard Bayesian state estimation [6], [10], the (multi-object) posterior density can be propagated forward recursively by where g k (·|·) is the multi-object likelihood function at time k, f k|k−1 (·|·) is the multi-object transition density to time k, and h k (Z k |Z 1:k−1 ) is the normalising constant, also known as the predictive likelihood.
A valid f k|k−1 (·|·) ensures each surviving object keeps the same label and dead labels never reappear [28], so that the multi-object history X 0:k represents a set of trajectories.
Markov Chain Monte Carlo (MCMC) approximations of the posterior have been proposed in [29] and [8] for detection and image measurements respectively. Combining MCMC with the generic multi-object particle filter [27] has also been suggested in [13].
A cheaper alternative is the multi-object filtering density, π k (X k ) π 0:k (X 0:k )δX 0:k−1 , which can be propagated by the multi-object Bayes filter [14], [16] Under the standard multi-object system model, the filtering recursion (7) admits an analytic solution known as the Generalized labeled Multi-Bernoulli (GLMB) filter [28], [30]. For a general system model, the generic multi-object particle filter can be applied, see for example [19].

C. Multi-object System model
Given a multi-object state X k−1 (at time k − 1), each state x k−1 = (x k−1 , ℓ k−1 ) ∈ X k−1 either survives with probability P S,k−1 (x k−1 ) and evolves to a new state (x k , ℓ k ) with probability density . Further, for each ℓ k in a (finite) birth label space B k at time k, either a new object with state (x k , ℓ k ) is born with probability P B,k (ℓ k ) and density f B,k (x k , ℓ k ), or unborn with probability Q B,k (ℓ k ) = 1 − P B,k (ℓ k ). The multi-object state X k (at time k) is the superposition of surviving states and new born states, and the multi-object transition density f k|k−1 (X k |X k−1 ) is given by equation (6) in [30]. An alternative form (using multi-scan exponential notation introduced in the next section) is given in subsection III-A.
Given a multi-object state X k , each x k ∈ X k is either detected with probability P D,k (x k ) and generates a detection z with likelihood g D,k (z|x k ) or missed with probability Q D,k (x k ) = 1−P D,k (x k ). The multiobject observation Z k is the superposition of the observations from detected objects and Poisson clutter with intensity κ k . Assuming that, conditional on X k , detections are independent of each other and clutter, the multi-object likelihood function is given by [28], [30] , Θ k denotes the set of positive 1-1 maps (i.e. those that never assign distinct arguments to the same positive value) from L k to {0:|Z k |}, and Θ k (I) denotes the subset of Θ k with domain I. The map θ k assigns a detected label ℓ to measurement z θk(ℓ) ∈ Z k , while for an undetected label θ k (ℓ) = 0.
The GLMB family is closed under the Bayes multi-object filtering recursion (7) and an explicit expression relating the filtering density at time k to that at time k − 1 is given by (14) of [31]. This recursion can be expressed in the following form 1 , which complies with, and facilitates the generalisation to the posterior recursion: Given the GLMB filtering density at time k − 1, GLMB filtering density at time k is given by The number of components of the GLMB filtering density grows super-exponentially in time.
Truncation by discarding components with small weights minimizes the L 1 approximation error in the multi-object density [30]. This can be achieved by solving the ranked assignment problems using Murty's algorithm or Gibbs sampling [31].

III. GLMB POSTERIOR RECURSION
In this section we extend the GLMB model to the multi-scan case, and subsequently derive an analytic recursion for the multi-scan GLMB posterior.

A. Multi-scan GLMB
Recall the equivalence between the multi-object state sequence X j:k and the set (5). For any function h taking the trajectories to the non-negative reals we introduce the following so-called multi-scan exponential notation: Note from (2) that a trajectory x (ℓ) s(ℓ):t(ℓ) is completely characterised by ℓ and the kinematic states x The multi-scan exponential notation is quite suggestive since properties). It also provides an intuitive expression for the multi-object transition density in [28].
Proposition 1: For the multi-object dynamic model described in subsection II-C, the multi-object transition density is given by For completeness the proof is given in Appendix VII-C.
It is clear that the multi-scan GLMB (density) reduces to a GLMB (density) when j = k.
Similar to the GLMB, the multi-scan GLMB (22) can be expressed in the so-called δ-form: where ξ ∈ Ξ, I j:k ∈ F(L j )×...×F (L k ). Each term or component of a multi-scan GLMB consists of a weight w (ξ) (I j:k ) and a multi-scan exponential p (ξ) Xj:k with label history that matches I j:k . The weight w (ξ) (I j:k ) can be interpreted as the probability of hypothesis (ξ, I j:k ), and for each ℓ ∈ ∪ k i=j L(I i ), p (ξ) (x s(ℓ):t(ℓ) ; ℓ) is the joint probability density of its kinematic states, given hypothesis (ξ, I j:k ).
By setting f to 1 in the above proposition, the multi-scan GLMB integrates to 1, and hence, is a FISST density. Some useful statistics from the multi-scan GLMB follows from the above proposition for suitably defined functions of the labels.

Corollary 4:
The cardinality distribution, i.e. distribution of the number of trajectories is given by Corollary 5: The joint probability of existence of a set of trajectories with labels L is given by As a special case, the probability of existence of trajectory with label ℓ is Corollary 6: The distribution of trajectory lengths is given by and the distribution of the length of trajectory with label ℓ is Similar to its single-scan counterpart, a number of estimators can be constructed for a multi-scan GLMB. The simplest would be to find the multi-scan GLMB component with the highest weight w (ξ) (I j:k ) and compute the most probable or expected trajectory estimate from p (ξ) (·; ℓ) for each ℓ ∈ ∪ k i=j L(I i ). Alternatively, instead of the most significant, we can use the most significant amongst components with the most probable cardinality n * (determined by maximizing the cardinality distribution (27)).
Another class of estimators, based on existence probabilities, can be constructed as follows. Find the set of labels L * with highest joint existence probability by maximizing (28). Then for each ℓ ∈ L * determine the most probable length m * by maximizing (31) and compute the trajectory density from which the most probable or expected trajectory estimate can be determined. Alternatively, instead of the label set with highest joint existence probability, we can use the label set of cardinality n * with highest joint existence probability. Another option is to find the n * labels with highest individual existence probabilities and use the same strategy for computing the trajectory estimates.

B. Multi-scan GLMB Posterior Recursion
Just as the GLMB is closed under the filtering recursion (7), the multi-scan GLMB is closed under the posterior recursion (6). Moreover, the multi-scan GLMB posterior recursion is, in essence, the GLMB filtering recursion without the marginalization of past labels and kinematic states. This is stated more concisely in the following Proposition (see Appendix VII-E proof).

Proposition 7:
Under the standard multi-object system model, if the multi-object posterior at time k −1 is a multi-scan GLMB of the form where ξ ∈ Ξ, then the multi-object posterior at time k is the multi-scan GLMB: Note that p (ξ,θk) The multi-scan GLMB posterior recursion (33)-(34) bears remarkable resemblance to the GLMB filtering recursion (9)- (10). Indeed, the weight increments for multi-scan GLMB and GLMB components are identical. Arguably, the multi-scan GLMB recursion is more intuitive because it does not involve marginalization over previous label sets nor past states of the trajectories.
The multi-scan GLMB recursion initiates trajectories for new labels, update trajectories for surviving labels, terminates trajectories for disappearing labels, and stores trajectories that disappeared earlier.
, the update of trajectories for surviving labels is the same as the GLMB filter, but without marginalization of past kinematic states. On the other hand, termination/storing of trajectories for disappearing/disappeared labels are not needed in the GLMB filter.

C. Cannonical Multi-scan GLMB Posterior
Without summing over the labels nor integrating the probability densities of the trajectories, the canonical expression for the multi-scan GLMB posterior takes on a rather compact form. To accomplish this, we represent each θ k ∈ Θ k by an extended association map γ k : Let Γ k denote the set of positive 1-1 maps from L k to {−1:|Z k |}, and (with a slight abuse of notation) denote the live labels of γ k , i.e. the domain D(θ k ), by Then for any γ k ∈ Γ k , we can recover θ k ∈ Θ k by θ k (ℓ) = γ k (ℓ) for each ℓ ∈ L(γ k ). It is clear that there is a bijection between Θ k and Γ k , and hence θ 1:k can be completely represented by γ 1:k .
Starting with an empty initial posterior π 0 (X 0 ) = δ 0 [L(X 0 )], by iteratively applying Proposition 7, the posterior at time k is given by where L(γ 0 ) ∅, Note that the multi-scan GLMB (40) is completely parameterized by the components (w L p -norm (of the difference) or information theoretic divergences such as Kullback-Leibler [20], Renyi, Cauchy-Schwarz [12], [3] and so on, can be extended to the multi-scan case by replacing the set integral with multiple set integrals.
Computing the multi-scan GLMB posterior first and then approximating it (according to a certain criterion) is not tractable. The main challenge is: given a prescibed number of terms, what is the best multi-scan GLMB approximation of the posterior, without having to compute all of its terms.
In this work we consider the L 1 -norm criterion, for which an optimal multi-scan GLMB approximation with a prescribed number of terms can be determined (see next section). Using the same lines of arguments as Proposition 5 of [30], the L 1 -error between a multi-scan GLMB and its truncation is given by the following result.
Hence, given a multi-scan GLMB, the minimum L 1 -norm approximation for a prescribed number of terms can be obtained by keeping only those with highest weights. Furthermore, this can be accomplished, without exhaustively computing all terms of the posterior, by solving the multi-dimensional or multi-frame assignment problem [21]. However, this problem is NP-hard for more than two scans. The most popular approximate solutions are based on Lagrangian relaxation [22], [23], [24], which is still computationally intensive. In this work we propose a more efficient algorithm by extending the Gibbs sampler of [31] to the multi-dimensional case.

IV. IMPLEMENTATION
In this section we consider problem of finding significant components (characterized by their association history) of the multi-scan GLMB (39) by sampling from some discrete probability distribution π. To ensure that mostly high-weight components are selected, π should be constructed so that association histories with high weights are more likely to be chosen than those with low weights. A natural choice to set with L(γ j ) = ∅, so that where w (γ0:k) 0:k is given by (40).
We present two techniques for sampling from (45). The first is based on sampling from the factors (44), i.e., γ j ∼ π (j) (·|γ 0:j−1 ), for j = 1 : k, using the Gibbs sampler of [31]. The alternative is a full Gibbs sampler with (45) as the stationary distribution.
end This approach ensures that the sample γ 1:k is a valid association history. However, we have to run each Gibbs sampler sufficiently long to ensure that sample γ j is close enough to being a sample from π (j) (·|γ 0:j−1 ) so that γ 1:k is close enough to being a sample from (45). Nonetheless, sampling from the factors can be use to generate a good starting point for the full Gibbs sampler.
Observe from (44) and (45) that for a valid γ 1:k , i.e. π(γ 1:k ) > 0, it is necessary that 1 Γi (γ i ) = 1 (i.e. γ i is positive 1-1), and 1 F (Bi⊎L(γi−1)) (L(γ i )) = 1 (i.e. dead labels at i − 1 cannot be live at i, or equivalently, a live label at i cannot be dead at i − 1) for i = 1 : k. Thus, in addition to being positive 1-1, consecutive elements of a valid γ 1:k must be such that dead labels remains dead at the next time.
Closed form expressions for the conditionals are given in the following Proposition (see Appendix VII-F for proof).
Proposition 10: Starting from any valid initial state, the Gibbs sampler defined by the family of conditionals (46) converges to the target distribution (45) at an exponential rate. More concisely, let π j denote the jth power of the transition kernel, then where, h = k + 1, β min γ1:k,γ ′ 1:k ∈Γk π h (γ ′ 1:k |γ 1:k ) > 0 is the least likely h-step transition probability.
The proof follows along the same line as Proposition 4 of [31], with the 2-step transition probability replaced by the (k + 1)-step transition probability. Instead of going from one arbitrary state of the chain to another via the all-zeros state in 2 steps as in [31], in this case we go to the all-negative state (consists of all -1) in k steps or less, and from this state to the other state in one additional step.
Similar to Gibbs sampling for the 2D assignment problem, the proposed Gibbs sampler has a relatively short burn-in period. For the purpose of approximating the multi-scan GLMB posterior density, it is not necessary to wait for samples from the stationary distribution since each distinct sample constitutes one term in the approximant, and reduces the L 1 approximation error by an amount proportional to its weight.
Recall that one of the proposed estimator is based on the most significant γ 1:k . The full Gibbs sampler above can be used in the simulated annealing setting to find the best γ 1:k more efficiently.

C. Computing Multi-Scan GLMB Components
This subsection details the computations the multi-scan GLMB (39) parameters for linear Gaussian multi-object models, i.e.
where the single-object state x is a d-dimensional vector, N (·; m, P ) denotes a Gaussian density with mean m and covariance P .
It follows from (42), (43) that the single object densities p Otherwise ℓ ∈ L(γ k ), and we have the following cases.
Case: s(ℓ) < k, then p  Note that the corresponding prediction density to time k is given by the marginal

V. NUMERICAL EXPERIMENTS
This section presents a preliminary numerical study demonstrating the improvements in tracking performance of the proposed multi-scan GLMB tracker over the single-scan counterpart [31]. For this purpose we use the linear Gaussian scenario shown in Section IV(A) of [31]. This scenario involves an unknown and time varying number of targets (up to 10 in total) over 100 time steps with births, deaths and crossings. Individual object kinematics are described by a 4D state vector of position and velocity that follows a constant velocity model with sampling period of 1s, and process noise standard deviation σ ν = 5m/s 2 . The survival probability P S = 0.99, and the birth model is an LMB with parameters  The standard GLMB tracker [31] is run with a maximum 10000 components. The proposed multi-scan GLMB tracker is run for 100 iterations of the multi-scan Gibbs sampler. The multi-target filtering errors are evaluated using the OSPA metric with parameters c = 100m and p = 1 as shown in Figure 1.
In addition the multi-target tracking errors are evaluated using the OSPA (2) metric [4] with the same parameters and a window length of 10 time steps. Observe from Figure 1 that there is a significant

VI. CONCLUSIONS
By introducing the multi-scan GLMB model, we extended the GLMB filtering recursion to propagate the labeled multi-object posterior density, i.e. multi-object smoothing-while-filtering. Conceptually the multi-scan GLMB recursion is more intuitive, but numerically it is far more challenging than the (singlescan) GLMB recursion. We showed that computing the multi-scan GLMB posterior with minimal L 1 -error (from its exact value) requires solving a multi-dimensional assignment problem with very high dimensions.
Further, we developed an efficient and highly parallelizable algorithm for solving such multi-dimensional assignment problems using Gibbs sampling, and subsequently a novel multi-object smoothing-whilefiltering algorithm. Numerical multi-object tracking examples demonstrated that the proposed algorithm significantly improves tracking performance as well as eliminating track fragmentation, a problem often found in multi-object filters. This is possible because the labeled multi-object posterior contains all information on the entire history of the underlying trajectories.

A. Properties of Multi-scan Exponentials
To present relevant properties of multi-scan exponentials, we introduce some useful partitionings for the labels of the multi-object state sequence X j:k . Given a time i in {j : k}, a label ℓ ∈ ∪ k r=j L(X r ) is alive at i iff ℓ ∈ L(X i ), terminates at t(ℓ) < i (before i) iff ℓ ∈ L(X t ) − L(X t+1 ), and born at time s(ℓ) > i (after i) iff ℓ ∈ L(X s ) ∩ B s . The set of labels in X j:k can be partitioned into labels terminated before i, live labels at i, and labels born after i, i.e. k r=j L(X r ) = When i = k, − −−− → L(X k ) = ∅ and the set of labels in X j:k can be partitioned into labels terminated before k and live labels at k, i.e. (49) becomes In addition, if j = k − 1, then ← −−− − L(X k ) = L(X k−1 ) − L(X k ), and using the decomposition L(X k ) = (L(X k−1 ) ∩ L(X k ))⊎(B k ∩ L(X k )), the set of labels in X k−1:k can be partitioned into labels terminated at k − 1, labels survived to k, and labels born at k, i.e. (51) becomes When i = j, ← −−− − L(X j ) = ∅ and the set of labels in X j:k can be partitioned into live labels at j and labels born after j, i.e. (50) becomes The following Lemma summarizes some useful properties of multi-scan exponentials.
Lemma A.1: Let X j:k be a sequence of multi-object states (generated by a set of trajectories) and g, h be two functions taking trajectories to the reals. Then:  (iv) For any i in {j : k} Proof: (i) and (ii) follows straight from the definition of multi-scan exponential.
To prove (iii) noting from (48) that the set of labels ∪ k r=j L(X r ) can be partitioned into those terminated before i, live at i, and born after i, we partition the set X j:k of trajectories accordingly, i.e.
Hence, using (ii) gives (iv) Using the corollaries to part (iii), we partition [g] Xj:i and [h] Xi:k , and then combine them as follows

C. Proof of Proposition 1 (Multi-object Transition)
Using (52) to partition L(X k−1 ) ∪ L(X k ) into disappearing labels at time k − 1, surviving labels at time k, and new born labels at time k, and noting that we have k , ℓ) to denote the element of the multi-object state X k at time k, with label ℓ ∈ L(X k ), then the multi-object transition density given in [28], [30] can be rewritten as where the last step follows from Lemma A.1 part (ii).

F. Proof of Proposition 9
We note the following conditions.