Poisson Multi-Bernoulli Approximations for Multiple Extended Object Filtering

1 Poisson Multi-Bernoulli Approximations for Multiple Extended Object Filtering Y uxuan Xia, Karl Granstr ¨ om, Lennart Svensson, Maryam F atemi, ´ Angel F . Garc ´ ıa-Fern ´ andez, Jason L. W illiams Abstract —The Poisson multi-Bernoulli mixtur e (PMBM) is a multi-object conjugate prior for the closed-f orm Bayes random finite sets filter . The extended object PMBM filter provides a closed-form solution for multiple extended object filtering with standard models. This paper considers computationally lighter alternativ es to the extended object PMBM filter by propagating a Poisson multi-Bernoulli (PMB) density thr ough the filtering recursion. A new local hypothesis repr esentation is presented where each measurement creates a new Bernoulli component. This facilitates the developments of methods f or efficiently approximating the PMBM posterior density after the update step as a PMB. Based on the new hypothesis repr esentation, two approximation methods are pr esented: one is based on the track- oriented multi-Bernoulli (MB) approximation, and the other is based on the variational MB approximation via K ullback–Leibler diver gence minimisation. The perf ormance of the pr oposed PMB filters with gamma Gaussian inv erse-Wishart implementations are evaluated in a simulation study . Index T erms —Multi-object filtering, extended object, random finite sets, multi-object conjugate prior , multi-Bernoulli, Kull- back–Leibler divergence, Gaussian in verse-Wishart I . I N T RO D U C T I O N Multi-object tracking (MO T) is the processing of sequences of noisy measurements obtained from multiple sources to esti- mate the sequences of object states [1]–[3]. Con ventional MO T algorithms are usually tailored to the point object assumption: each object is modelled as a point without spatial extent and giv es rise to at most one measurement at each time step. Howe ver , the point object assumption becomes unrealistic if multiple resolution cells of the sensor are occupied by a single object such that each object may giv e rise to multiple measurements at each time step. T ypical examples of such scenarios include, for example, vehicle and pedestrian tracking using automotiv e radar and laser range sensors. This paper considers the multiple extended object filtering problem where the objectiv e is to estimate both the kinematic and extent states of the current set of objects. An overvie w of the extended object tracking literature can be found in [4]. The detections of an extended object are commonly modelled by an inhomogeneous Poisson point process (PPP), proposed in [5]. In the PPP measurement model, the number of detections is Y . Xia, K. Granstr ¨ om and L. Svensson are with the Department of Electrical Engineering, Chalmers University of T echnology , G ¨ oteborg, Sweden. E- mail: firstname.lastname@chalmers.se. M. Fatemi is with Zenseact, G ¨ oteborg, Sweden. E-mail: maryam.fatemi@zenseact.com. A. F . Garc ´ ıa-Fern ´ andez is with the Department of Electrical Engineering and Electronics, Univ ersity of Liverpool, Li verpool, United Kingdom and also with the ARIES research centre, Universidad Antonio de Nebrija, Madrid, Spain. Email: angel.garcia- fernandez@liv erpool.ac.uk. J. L. Williams is with the Commonwealth Sci- entific and Industrial Research Organization, Brisbane, Australia. Email: jason.williams@data61.csiro.au. Poisson distributed, and the detections are spatially distrib uted around the object. The PPP measurement model has been integrated into a number of multiple e xtended object filters including, for example, the extended object joint probabilistic data association (JPD A) filter [6]–[10], the extended object multiple hypothesis tracker [11], [12], filters based on belief propagation [13], [14] and sev eral random finite sets (RFSs) based extended object filters [15]–[20]. It is also possible to model the object detections using a cluster process [21], and it has been integrated into the probability hypothesis density filter [22]–[28]. RFSs facilitate an elegant Bayesian formulation of the MO T problem, and a significant trend in RFSs-based MO T is the dev elopment of multi-object conjugate priors. W ith PPP birth model, the Poisson multi-Bernoulli mixture (PMBM) is a multi-object conjugate prior for general measurement models [29], including, for example, point objects [30] and extended objects with PPP measurement model [20], in the sense that the PMBM form of the multi-object density is preserved through the filtering recursion. The simulation studies in [20], [31] hav e shown that the extended object PMBM filter , in general, performs well in comparison to other RFSs-based extended object filters. The PMBM recursion can also handle a multi-Bernoulli (MB) birth model by setting the PPP intensity to zero and adding new Bernoulli components in the prediction step, resulting in the MB mixture (MBM) filter [32], [33]. The MBM filter can be further extended to consider MBs with deterministic object existence, referred to as MBM 01 , at the expanse of an exponential growth in the number of MBs. The labelled MBM 01 filtering recursion is analogous to the generalised labelled multi-Bernoulli (GLMB) filter [19], [34]. The reasoning ov er the correspondence of measurements and objects, also known as data association, is an inherent challenge in MOT . The unknown data associations lead to an intractably large number of terms in the PMBM density . T o achieve computational tractability of the extended object PMBM filter , it is necessary to reduce the number of PMBM parameters, see [20, Sec. V -B] for different reduction methods. The computational complexity of the PMBM filter, howe ver , still scales with the number of MBs. It is therefore of interest to consider a special variant of the extended object PMBM filter that only propagates a Poisson MB (PMB) density through the filtering recursion. This paper focuses on developing ef ficient PMB approxima- tions of the PMBM multi-object posterior density for multiple extended object filtering. W e first identify the major challenges in extended object PMB approximation: 1) determining the 2 number Bernoulli components in the approximate PMB; 2) de- termining which local hypothesis densities should be merged; 3) determining how to merge the selected local hypothesis densities into one. W e then proceed to present a more efficient local hypothesis representation of the PMBM conjugate prior where each measurement creates a new Bernoulli component, whereas for the hypothesis structure in [20], [29] each non- empty subset of measurements creates a ne w Bernoulli compo- nent. This facilitates the de velopment of PMB approximation methods tailored to extended objects. Importantly , the new hypothesis representation results in an equiv alent PMBM representation of the multi-object density but with far fe wer Bernoulli components. W e present two extended object PMB filters that approx- imate the extended object PMBM filter using two different PMB approximation methods. The first PMB filter is based on the track-oriented MB (TO-MB) approximation, which was first presented in [30] for point objects and shown to be similar to the JPD A approximation. The idea of using TO- MB approximation for extended object filtering has recently been extended in [29] to handle scenarios with coexisting point and extended objects, b ut with a less efficient hypothesis representation where each non-empty subset of measurements creates a new Bernoulli component. The second PMB filter is based on variational approximation that aims to minimise the Kullback–Leibler di ver gence (KLD) between the PMBM and the approximate PMB. T wo different formulations of the minimisation objective are studied: one is based on finding the most likely permutation of local hypotheses for each global hypothesis, and the other follo ws the linear programming formulation in [35]. Implementations of the proposed extended object PMB filters, based on gamma Gaussian in verse-W ishart (GGIW) single object models [36], are also presented. Simu- lation results show that the extended object PMB filters offer an appealing trade-off between computational complexity and estimation performance. The remainder of the paper is organised as follows. In Section II, we overvie w the e xtended object PMBM conjugate prior with standard models and present the problem formula- tion. The extended object PMB filtering recursion is presented in Section III. The new local hypothesis representation of the PMBM conjugate prior and the TO-MB approximation are presented in Section IV . The v ariational MB approximation and its two different formulations are presented in Section V . The GGIW implementations of the PMB filters are giv en in Section VI. The simulation results are shown in Section VII and the conclusions are drawn in Section VIII. I I . B A C K G RO U N D A N D P RO B L E M F O R M U L A T I O N In this section, we first briefly introduce Bayesian MO T and the standard models for multiple extended object tracking. W e then o vervie w the e xtended object PMBM conjugate prior b uilt upon these models. Last, the problem formulation is presented. A. Bayesian multi-object filtering The set of objects at time step k is denoted X k . By mod- elling X k as a RFS, both object set cardinality | X k | and single object state x k ∈ X k are random variables. In e xtended object tracking, the object state models both the kinematic properties, and the extent, of the object. The set of measurements at time step k is denoted Z k =  z 1 k , . . . , z m k k  , including clutter and object measurements. Further , Z k = ( Z 1 , . . . , Z k ) denotes the sequence of measurements from time step 1 to k . The multi-object filtering density is denoted f k | k  X k | Z k  . The multi-object Bayes filter propagates in time the multi- object density f k − 1 | k − 1  X k − 1 | Z k − 1  using the Chapman- K olmogorov prediction f k | k − 1  X k | Z k − 1  = Z f k,k − 1 ( X k | X k − 1 ) f k − 1 | k − 1  X k − 1 | Z k − 1  δ X k − 1 , (1) and then updates the density using the Bayes’ update f k | k ( X k | Z k ) = f k ( Z k | X k ) f k | k − 1 ( X k | Z k − 1 ) R f k ( Z k | X k ) f k | k − 1 ( X k | Z k − 1 ) δ X k , (2) where f k,k − 1 ( X k | Z k − 1 ) is the multi-object transition density , f k ( Z k | X k ) is the multi-object measurement set density , and R f ( X ) δ X is set inte gral, defined in [37, Sec. 11.3.3]. B. Multi-object dynamic and measur ement models 1) Multi-object dynamic model: Each object surviv es from time step k − 1 to time step k with a probability of surviv al p S ( x k − 1 ) . Each surviv ed object ev olves independently with a Markovian transition density π k,k − 1 ( x k | x k − 1 ) . New objects appear independently of the objects that already exist. The object birth is assumed to be a PPP with intensity D b k ( x ) . 2) Extended object measurement model: The set of mea- surements Z k is a union of a set of clutter measurements and a set of object measurements. The clutter is modelled as a PPP with intensity κ ( z ) = λc ( z ) , where λ is the Poisson clutter rate and c ( z ) specifies the spatial density of clutter measurements. An extended object with state x k is detected with detection probability p D ( x k ) , and if detected, the object measurements are modelled as a PPP with Poisson rate γ ( x k ) and spatial distribution φ ( ·| x k ) , independent of any other objects and their corresponding generated measurements. The conditional extended object measurement set likelihood for a non-empty set of measurements is the product of the probability of detection and the PPP density , ` Z ( x k ) = p D ( x k ) e − γ ( x k ) Y z ∈ Z γ ( x k ) φ ( z | x k ) . (3) For an e xtended object state x k , the effecti ve probability of detection is the product of the probability of detection and the probability that object generates at least one measurement 1 − e − γ ( x k ) . Therefore, the effecti ve probability that the object x k is not detected is q D ( x k ) = 1 − p D ( x k ) + p D ( x k ) e − γ ( x k ) . (4) Note that the conditional likelihood for an empty measurement set, ` ∅ ( x k ) , is also described by (4). 3 C. PMBM conjugate prior For the abov e multi-object dynamic and measurement mod- els, the multi-object density f k | k 0 ( · ) of X k giv en the sequence of measurements Z k 0 , where k 0 ∈ { k − 1 , k } , is a PMBM [20]. The PMBM density is of the form [30], [38] f pmbm k | k 0 ( X k ) = X X u k U X d k = X k f ppp k | k 0 ( X u k ) f mbm k | k 0  X d k  , (5a) f ppp k | k 0 ( X u k ) = e − R D u k | k 0 ( x ) d x Y x ∈ X u k D u k | k 0 ( x ) , (5b) f mbm k | k 0  X d k  = X a ∈A k | k 0 w a k | k 0 X U n k | k 0 l =1 X l = X d k n k | k 0 Y i =1 f i,a i k | k 0  X i  , (5c) where U denotes the disjoint union, D u k | k 0 ( · ) is the intensity of the PPP f ppp k | k 0 ( · ) represent objects that ha ve ne ver been detected (referred to as undetected objects), and f mbm k | k 0 ( · ) is an MBM representing potential objects that ha ve been detected at some point up to time step k 0 . The set of global hypotheses is de- noted A k | k 0 , and each global hypothesis a =  a 1 , . . . , a n k | k 0  is giv en by selecting a local hypothesis a i ∈ n 1 , . . . , h i k | k 0 o for each Bernoulli component i ∈ { 1 , . . . , n k | k 0 } subject to certain association constraints [20], [38]. The global hypothe- sis weight w a k | k 0 is w a k | k 0 ∝ n k | k 0 Y i =1 w i,a i k | k 0 , (6) where w i,a i k | k 0 is the weight of local hypothesis a i , and the proportionality denotes that normalisation is required to en- sure that P a ∈A k | k 0 w a k | k 0 = 1 . In the MBM (5c), there are n k | k 0 Bernoulli components, and for each Bernoulli there are h i k | k 0 possible local hypotheses. The i -th Bernoulli with local hypothesis a i has density f i,a i k | k 0 ( X ) =        1 − r i,a i k | k 0 X = ∅ , r i,a i k | k 0 p i,a i k | k 0 ( x ) X = { x } , 0 otherwise , (7) where r i,a i k | k 0 is the probability of existence and p i,a i k | k 0 ( · ) is the existence-conditioned single object state density . D. Pr oblem F ormulation W e aim to de velop an extended object filter that propagates a PMB density over time, i.e., a special case of the PMBM density (5) that has an MBM with a single MB. W e denote the PMB predicted/filtering density at time step k , with k 0 ∈ { k − 1 , k } , as f pmb k | k 0 ( X k ) = X U n k | k 0 l =1 X l U Y = X k f ppp k | k 0 ( Y ) n k | k 0 Y i =1 f i k | k 0  X i  , (8) where f ppp k | k 0 ( · ) is of the form (5b) with intensity D u k | k 0 ( · ) , n k | k 0 is the number of Bernoulli components, and the i -th Bernoulli component has probability of e xistence r i k | k 0 and e xistence- conditioned single object state density p i k | k 0 ( · ) . Giv en the PMB filtering density f pmb k − 1 | k − 1 ( X k − 1 ) , the pre- dicted density f pmb k | k − 1 ( X k ) retains the PMB form [30]. The Bayes’ update of f pmb k | k − 1 ( X k ) , howe ver , results in a PMBM f pmbm k | k ( X k ) due to the unknown data associations. T o obtain a PMB filtering recursion, the true filtering density f pmbm k | k ( X k ) needs to be approximated by a PMB f pmb k | k ( X k ) , which should retain as much information from f pmbm k | k ( X k ) as possible. A natural choice for solving this problem is to minimise the KLD between the PMBM and the approximate PMB, min f pmb D KL  f pmbm k | k ( X k )    f pmb k | k ( X k )  , min f pmb Z f pmbm k | k ( X k ) log f pmbm k | k ( X k ) f pmb k | k ( X k ) ! δ X k . (9) Analytically minimising (9) is intractable, and thus we hav e to use approximations to obtain an efficient algorithm. The above problem can be simplified by noting that the PPP representing undetected objects is not hypothesis-dependent, and thus does not require approximation [30]. It is therefore sufficient to consider approximating the MBM f mbm k | k  X d k  representing detected objects in f pmbm k | k ( X k ) as an MB: arg min f mb D KL  f mbm k | k  X d k     f mb k | k  X d k   . (10) Lastly , following the recycling method in [39], we note that Bernoulli components with a small probability of existence are well approximated as a PPP . Practically , this allows these Bernoulli components to be modelled through the PPP , reduc- ing the complexity of subsequent data association problems. I I I . E X T E N D E D O B J E C T P M B F I L T E R In this section, we present the track-oriented e xtended object PMB filtering recursion, which includes prediction, update, MB approximation and rec ycling. Compared to the hypothesis- oriented filtering recursion [20], the track-oriented filtering recursion [38] has more efficient hypothesis management and facilitates the de velopment of MB approximation methods. A. Pr ediction W e denote the inner product of a ( x ) and b ( x ) as h a ; b i , R a ( x ) b ( x ) dx for notational con venience. Lemma 1 (PMB prediction) . Given the PMB filtering density at time step k − 1 of the form (8) and the multi-object dynamic model in Section II-B1, the pr edicted density at time step k is a PMB of the form (8), with n k | k − 1 = n k − 1 | k − 1 and D u k | k − 1 ( x ) = D b k ( x ) + D D u k − 1 | k − 1 ; π k,k − 1 ( x |· ) p S ( · ) E , r i k | k − 1 = r i k − 1 | k − 1 D p i k − 1 | k − 1 ; p S E , p i k | k − 1 ( x ) = D p i k − 1 | k − 1 ; π k,k − 1 ( x |· ) p S ( · ) E D p i k − 1 | k − 1 ; p S E . 4 Lemma 1 is a special case of the PMBM prediction [20], [30] as a PMB is a PMBM with only one mixture component. B. Update W e refer to measurement z j k using the index pair ( k , j ) , and the set of all such index pairs at time step k is denoted M k . W e also let M i,a i k ⊆ M k be the set of index pairs associated to local hypothesis a i of the i -th Bernoulli component. In the measurement update step of an extended object PMB filter , ev ery measurement should be assigned to one and only one local hypothesis in a global hypothesis, and more than one measurement can be assigned to the same local hypothesis at the same time step. The set of global hypotheses is then A k | k = (  a 1 , . . . , a n k | k  : a i ∈ n 1 , . . . , h i k | k o ∀ i, n k | k ] i =1 M i,a i k = M k ) . (12) This means that each global hypothesis, represented as an MB, corresponds to a unique association of measurements Z k to n k | k − 1 Bernoulli components for previously detected objects and n k | k − n k | k − 1 new Bernoulli components for newly detected objects. Lemma 2 (PMB update) . Given the PMB pr edicted density at time step k of the form (8), the multi-object measure- ment model in Section II-B2, and a measur ement set Z k =  z 1 k , . . . , z m k k  , the updated distribution is a PMBM of the form (5), with n k | k = n k | k − 1 + 2 m k − 1 and D u k | k ( x ) = q D ( x ) D u k | k − 1 ( x ) . (13) F or each pre-e xisting Bernoulli component f i k | k − 1 ( · ) , i ∈ { 1 , . . . , n k | k − 1 } , there are h i k | k = 2 m k local hypotheses, corr esponding to a misdetection and an update using a non- empty subset of Z k . The misdetection hypothesis for Bernoulli component i ∈ { 1 , . . . , n k | k − 1 } is M i, 1 k = ∅ , (14a) w i, 1 k | k = 1 − r i k | k − 1 + r i k | k − 1 D p i k | k − 1 ; q D E , (14b) r i, 1 k | k = r i k | k − 1 D p i k | k − 1 ; q D E 1 − r i k | k − 1 + r i k | k − 1 D p i k | k − 1 ; q D E , (14c) p i, 1 k | k ( x ) = q D ( x ) p i k | k − 1 ( x ) D p i k | k − 1 ; q D E . (14d) Let Z 1 k , . . . , Z 2 m k − 1 k be the non-empty subsets of Z k . The local hypothesis for Bernoulli component i ∈ { 1 , . . . , n k | k − 1 } and measur ement set Z j k is M i,j k = n ( k , l ) : z l k ∈ Z j k o , (15a) w i,j k | k = r i k | k − 1 D p i k | k − 1 ; ` Z j k E , (15b) r i,j k | k = 1 , (15c) p i,j k | k ( x ) = ` Z j k ( x ) p i k | k − 1 ( x ) D p i k | k − 1 ; ` Z j k E . (15d) F or the new Bernoulli component ( i = n k | k − 1 + j ) initiated by measur ement set Z j k , ther e are h i k | k = 2 local hypotheses M i, 1 k = ∅ , w i, 1 k | k = 1 , r i, 1 k | k = 0 , (16a) M i, 2 k = n ( k , l ) : z l k ∈ Z j k o , (16b) w i, 2 k | k =    κ ( z ) + D D u k | k − 1 ; ` Z j k E | Z j k | = 1 , D D u k | k − 1 ; ` Z j k E | Z j k | > 1 , (16c) r i, 2 k | k =         D u k | k − 1 ; ` Z j k  κ ( z )+  D u k | k − 1 ; ` Z j k  | Z j k | = 1 , 1 | Z j k | > 1 , (16d) p i, 2 k | k ( x ) = ` Z j k ( x ) D u k | k − 1 ( x ) D D u k | k − 1 ; ` Z j k E . (16e) Local hypothesis (16a) corresponds to the case that a non- empty subset of Z j k is associated with another Bernoulli, hence this local hypothesis has probability of e xistence equal to zero and non-valid single object density . Local hypothesis (16b) corresponds to the case that the i -th Bernoulli component is created by the set Z j k of measurements. Lemma 2 is a special case of the PMBM update [20] by considering a PMB prior . C. MB appr oximation and its challenges T o complete the PMB filtering recursion, we would like to approximate the MBM (5c) with k 0 = k in the PMBM filtering density as an MB of the form f mb k | k ( X k ) = X U e n k | k l =1 X l = X k e n k | k Y i =1 e f i k | k ( X i ) . (17) There are mainly three challenges in this MB approximation problem. The first challenge is to determine the number of Bernoulli components e n k | k in the approximate MB f mb k | k ( X k ) . Each mixture component in (5c) corresponds to a unique global hypothesis with n k | k local hypotheses, but some of these local hypotheses may correspond to non-existent objects. Considering that 2 m k − 1 ne w Bernoullis are created at each time step, it may be inefficient to set e n k | k = n k | k − 1 + 2 m k − 1 . Since the Bernoulli components in (5c) are order-in variant, the second challenge is to determine which local hypothesis densities f i,a i k | k ( · ) ( a ∈ A k | k , i ∈ { 1 , . . . , n k | k } ) should be merged to form e f l k | k ( · ) for each l ∈  1 , . . . , e n k | k  . The third challenge is about how to merge the selected local hypothesis densities into one, such that each Bernoulli component in (5c) only corresponds to a single local hy- pothesis. The single Bernoulli approximation of a Bernoulli mixture can be obtained by taking the weighted sum of different parameters. The approximate Bernoulli component, without further approximation, then contains a single object density collecting mixture components that can be arbitrarily dissimilar . It is usually more desirable to represent the single 5 1 2 Predicted Bernoullis New Bernoullis 3 4 5 6 7 (a) Hypothesis structure (b) Global hypothesis a 1 (c) Global hypothesis a 2 (d) Global hypothesis a 3 Fig. 1. Illustrativ e example with three global hypotheses a 1 , a 2 and a 3 where a i denotes the index of the local hypothesis (index ed from left to right) of the i -th Bernoulli component. The local hypothesis densities are Gaussian and described by ellipses. Local hypotheses marked in red and blue represent the updates of predicted Bernoullis, whereas local hypotheses marked in green and magenta represent Bernoullis for the newly detected objects. object density using a unimodal distribution, which simplifies computation and permits closed-form ev aluation of objecti ves, such as the KLD. In extended object filtering, the merging of different single object densities is typically achieved by KLD minimisation, see, e.g., [40], [41]. W e illustrate these challenges with the following example. Example 1. Consider the scenario shown in Fig . 1, with an MBM r epr esenting thr ee global hypotheses. F or simplicity , we assume that ther e is no misdetection, and that each Bernoulli r epresenting an object has pr obability of existence equal to one. The two objects that are detected for the fir st time are closely-spaced, so it is difficult to tell whether ther e ar e two closely-spaced new objects or one larger new object. The first challenge is determining how to choose the number of Bernoullis in the appr oximate MB. In this example, the MBM contains seven Bernoulli components, but each global hypothesis contains local hypotheses corresponding to non- existent objects. Thus, it is non-trivial to select the number of Bernoulli components in the appr oximate MB. The second c hallenge is determining how to select the local hypothesis densities to be merged. F or example, when creating appr oximate Bernoulli components for pr eviously detected objects, we may mer ge f 1 , 1 ( · ) and f 2 , 2 ( · ) into one, but we may also mer ge f 1 , 1 ( · ) and f 1 , 2 ( · ) into one. The pr oblem becomes more complicated when we consider how to form appr oximate Bernoulli components corresponding to the newly detected objects. The thir d c hallenge is determining how to mer ge the selected local hypothesis densities. F or example , if we choose to mer ge f 1 , 1 ( · ) and f 2 , 2 ( · ) , we would like the information loss due to mer ging to be minimised. D. Recycling For the approximate MB (17), we can apply recycling [39] to Bernoulli components with low probability of existence. The recycled Bernoulli components are first approximated as a PPP in the sense of KLD minimisation, and then the approximate PPP is added to the PPP of the PMB filter . After recycling, we set the probability of existence of the recycled Bernoulli components in (17) to zero, and the PPP intensity of undetected objects is e D u k | k ( x ) = D u k | k ( x ) + X i ∈ { 1 ,..., e n k | k } : e r i k | k <τ r e r i k | k e p i k | k ( x ) , where τ r is the recycling threshold, and e r i and e p i ( · ) are the probability of e xistence and the single object density of Bernoulli e f i k | k ( · ) , respectiv ely . By having fe wer Bernoulli components in the MB for detected objects, the comple xity of the data association problem in next time step can be reduced. I V . T R A C K - O R I E N T E D M U LT I - B E R N O U LL I A P P RO X I M A T I O N In this section, we first present a class of possible local hypothesis representations in the PMBM conjugate prior for extended objects. This allows us to find a compact local hy- pothesis structure for Bernoulli components. Based on an im- prov ed local hypothesis representation for the ne wly detected objects of (16), we then present the TO-MB approximation of the MBM (5c) in the PMBM filtering density . W e refer to the resulting filter as track-oriented PMB (T O-PMB). The relations of TO-PMB to the labelled MB (LMB) and JPD A for extended objects are also discussed. 6 A. Local hypothesis structur e for newly detected objects Let F ( M k ) denote the set of all the subsets of M k . The following theorem describes a family of hypothesis represen- tations of the updated PMBM posterior that only differ in the representations of the new Bernoulli components; its proof is provided in Appendix A. Theorem 1. The updated PMBM density (5), whose set of global hypotheses meets (12), can be equivalently r epr esented by the same local hypotheses for pre viously detected objects (14), (15), and n k ≥ m k new Bernoulli components with local hypotheses whose associations satisfy [ i,j n M i,j k o = F ( M k ) , (18a) ∀ i ∃ j : M i,j k = ∅ , (18b) ∀ i, j, j 0 : j 6 = j 0 → M i,j k 6 = M i,j 0 k , (18c) ∀ i, i 0 , j, j 0 : i 6 = i 0 , j 6 = j 0 , M i,j k 6 = ∅ → M i,j k 6 = M i 0 ,j 0 k , (18d) ∀ i, j, j 0 : j 6 = j 0 , M i,j k 6 = ∅ , M i,j 0 k 6 = ∅ → M i,j k ∩ M i,j 0 k 6 = ∅ , (18e) wher e the range of variable i is { n k | k − 1 + 1 , . . . , n k | k } with n k | k = n k | k − 1 + n k . The density of a new Bernoulli component associated with non-empty M i,j k ( j ∈ n 1 , . . . , h i k | k o ) wher e M i,j k = n ( k , l ) : z l k ∈ M i,j k o is obtained as in (16c)-(16e) by letting Z j k = M i,j k . Compared to the representations of the new Bernoulli com- ponents (16) in Lemma 2 where each non-empty subset of measurements creates a new Bernoulli component, here we consider a more general case where several non-empty subsets may be assigned to the same new Bernoulli component. Con- straint (18a) is necessary since any non-empty subset of M k can correspond to the first detection of an undetected object. Constraint (18b) means that each new Bernoulli component may represent a non-existent object. Constraints (18c) and (18d) ensure that there do not exist two local hypotheses with the same non-empty subset of M k . At last, constraint (18e) ensures that two disjoint non-empty subsets of M k cannot be assigned to the same Bernoulli component. This is intuiti ve as local hypotheses created by measurement clusters with non- empty intersections can never co-exist in a global hypothesis. In other words, for every Bernoulli component i , where the range of v ariable i is { n k | k − 1 + 1 , . . . , n k | k } , there should be at least one measurement which is shared by all the local hypotheses assigned to the component. The constraints (18) can be interpreted as assigning 2 m k − 1 local hypotheses with non-empty subsets of M k to n k ≥ m k new Bernoulli components such that local hypotheses of the same Bernoulli component satisfy (18e). B. An alternative PMB update In the T O-MB approximation, all the Bernoulli components are forced to be independent, and local hypothesis densities of the same Bernoulli component are mer ged [32]. According to Lemma 2, for a predicted PMB with n k | k − 1 Bernoulli components and a set Z k of m k measurements, the TO-MB approximation of the updated distribution contains e n k | k = n k | k − 1 + 2 m k − 1 Bernoulli components. Howe ver , creating as man y new Bernoulli components as there are non-empty measurement subsets may yield intractably many Bernoulli components with low probability of existence. Because m k measurements can correspond to at most m k newly detected objects, the most efficient class of local hy- pothesis representations has n k | k = n k | k − 1 + m k . That is, each measurement creates a new Bernoulli component, as in the case of point object filtering. T o construct such a local hypothesis representation that satisfies (18), one strategy is to choose M i, 1 k = ∅ ∀ i > n k | k − 1 and ( k , i ) ∈ M i,a i k ∀ i > n k | k − 1 , a i ≥ 2 . Let M i,j − 1 k = n ( k , l ) : z l k ∈ M i,j − 1 k o ∀ j ≥ 2 where M i,j k is the j -th subset of measurements of the i - th Bernoulli component 1 . The set S i = n M i, 1 k , . . . , M i,i k o of subsets of measurements associated to local hypotheses of the i -th Bernoulli component can be built recursively as S i =  z i k  ∪    [ M ∈∪ i − 1 j =1 S j  z i k  ∪ M     with S 1 =  z 1 k  . Note that S i is a set of sets and  z i k  represents a set of single-measurement set  z i k  [42, Appendix D]. This results in an alternativ e PMB update, which gi ves the same PMBM posterior (5) as giv en by the PMB update in Lemma 2. Lemma 3 (Alternative PMB update) . Given the PMB pr e- dicted density at time step k of the form (8), the multi-object measur ement model in Section II-B2, and a measur ement set Z k =  z 1 k , . . . , z m k k  , the updated distribution can be written as a PMBM of the form (5), with n k | k = n k | k − 1 + m k and updated PPP intensity D u k | k ( x ) (13). F or each pr e-existing Bernoulli component f i k | k − 1 ( · ) , i ∈ { 1 , . . . , n k | k − 1 } , ther e ar e h i k | k = 2 m k local hypotheses, corr esponding to a misdetection and an update using a non-empty subset of Z k , see (14), (15). F or the new Bernoulli component ( i = n k | k − 1 + j ), there ar e h i k | k = 2 j − 1 + 1 local hypotheses ( ι ∈  1 , . . . , 2 j − 1  ) M i, 1 k = ∅ , w i, 1 k | k = 1 , r i, 1 k | k = 0 , (19a) M i,ι +1 k = n ( k , l ) : z l k ∈ M i,ι k o , (19b) w i,ι +1 k | k =    κ ( z ) + D D u k | k − 1 ; ` M i,ι k E | M i,ι k | = 1 , D D u k | k − 1 ; ` M i,ι k E | M i,ι k | > 1 , (19c) r i,ι +1 k | k =         D u k | k − 1 ; ` M i,ι k  κ ( z )+  D u k | k − 1 ; ` M i,ι k  | M i,ι k | = 1 , 1 | M i,ι k | > 1 , (19d) 1 Note that a dif ferent and varying number of local hypotheses is assigned to each Bernoulli component. This is necessary to correctly represent each association event. For example, if each Bernoulli component was assigned ev ery hypothesis containing the corresponding measurement, many local hypotheses w ould be duplicated over multiple Bernoulli components, leading to double-counting of association events. 7 1 2 3 4 5 6 7 (a) Hypothesis structure (Lemma 2) 1 2 3 (b) Hypothesis structure (Lemma 3) Fig. 2. An illustration of the two hypothesis structures in Example 2, given by Lemma 2 with 7 Bernoulli components (top) and Lemma 3 with 3 Bernoulli components (bottom), respectively , where the time index k is omitted from the local hypothesis representation. Global hypotheses a 1 and b 1 correspond to data partition {{ 1 } , { 2 } , { 3 }} ; global hypotheses a 2 and b 2 correspond to data partition {{ 1 } , { 2 , 3 }} ; global hypotheses a 3 and b 3 correspond to data partition {{ 1 , 2 } , { 3 }} ; global hypotheses a 4 and b 4 correspond to data partition {{ 1 , 3 } , { 2 }} ; and global hypotheses a 5 and b 5 correspond to data partition {{ 1 , 2 , 3 }} . p i,ι +1 k | k ( x ) = ` M i,ι k ( x ) D u k | k − 1 ( x ) D D u k | k − 1 ; ` M i,ι k E . (19e) For a better understanding of the new hypothesis structure, let us consider the following simple example. Example 2. Consider the case wher e n k | k − 1 = 0 and m k = 3 . Ther e ar e 5 differ ent ways to partition three measur ements, so the number of global hypotheses is 5. According to Lemma 2, 7 Bernoulli components ar e created, and the local hypotheses under these components are: 1) ∅ and { ( k , 1) } , 2) ∅ and { ( k , 2) } , 3) ∅ and { ( k , 3) } , 4) ∅ and { ( k , 1) , ( k , 2) } , 5) ∅ and { ( k , 1) , ( k , 3) } , 6) ∅ and { ( k , 2) , ( k , 3) } , 7) ∅ and { ( k , 1) , ( k , 2) , ( k , 3) } . The 5 global hypotheses r epr esented using these local hypotheses ar e: a 1 = (2 , 2 , 2 , 1 , 1 , 1 , 1) , a 2 = (2 , 1 , 1 , 1 , 1 , 2 , 1) , a 3 = (1 , 1 , 2 , 2 , 1 , 1 , 1) , a 4 = (1 , 2 , 1 , 1 , 2 , 1 , 1) and a 5 = (1 , 1 , 1 , 1 , 1 , 1 , 2) . According to Lemma 3, 3 Bernoulli components are created, and the local hypotheses under these components ar e: 1) ∅ and { ( k , 1) } , 2) ∅ , { ( k , 2) } and { ( k , 1) , ( k , 2) } , 3) ∅ , { ( k , 3) } , { ( k, 1) , ( k , 3) } , { ( k , 2) , ( k , 3) } and { ( k , 1) , ( k , 2) , ( k , 3) } . The 5 global hy- potheses repr esented using these local hypotheses corr espond- ing to a 1 , a 2 , a 3 , a 4 and a 5 ar e: b 1 = (2 , 2 , 2) , b 2 = (2 , 1 , 4) , b 3 = (1 , 3 , 2) , b 4 = (1 , 2 , 3) and b 5 = (1 , 1 , 5) , r espectively . Note that, compar ed to the original hypothesis r epresentation, the number of Bernoulli components r educes fr om 7 to 3, see also F ig. 2. C. Multi-Bernoulli appr oximation Giv en the PMBM filtering density resulting from Lemma 3, the TO-MB approximation is [30] f mb k | k ( X k ) = X U e n k | k l =1 X l = X k e n k | k Y i =1 e f i k | k ( X i ) , (20a) e n k | k = n k | k − 1 + m k , (20b) w a k | k ∝ e n k | k Y i =1 w i,a i k | k , (20c) e f i k | k ( X ) = h i X a i =1 X b ∈A k | k : b i = a i w b k | k f i,a i k | k ( X ) . (20d) The TO-MB approximation (20) can be seen to preserve the probability hypothesis density of PMBM filtering density [32]. It can be also shown that (20) minimises the KLD on a single object space augmented with an auxiliary v ari- able 2 , which indicates if the object remains undetected or corresponds to the i -th Bernoulli component [43]. Finally , it should be noted that the approximate MB filtering density (20) depends on the indices of the individual measurements z j k , with j ∈ { 1 , . . . , m k } . A different measurement indexing may result in a different approximate MB, but local hypothesis densities created by disjoint measurement sets will never be merged, regardless of the measurement index. As will be discussed in Section V -B, the problem of finding the local hypothesis representation, that leads to the most accurate TO- MB approximation in the sense of KLD minimisation, can be solved using variational approximation. Example 3. Consider the same scenario and the measur ement indexing shown in F ig. 1 and suppose that global hypotheses a 1 , a 2 and a 3 have weights w a 1 , w a 2 and w a 3 , r espectively . The PMBM filtering density resulting fr om Lemma 3 has two new Bernoulli components: the first one contains local hy- potheses f 4 , 2 ( · ) and f 6 , 2 ( · ) , whereas the second one contains local hypotheses f 3 , 2 ( · ) , f 5 , 2 ( · ) and f 7 , 2 ( · ) . The appr oximate MB obtained using (20) contains four Bernoulli components: e f 1 ( X ) = ( w a 1 + w a 2 ) f 1 , 1 ( X ) + w a 3 f 1 , 2 ( X ) , e f 2 ( X ) = ( w a 1 + w a 2 ) f 2 , 1 ( X ) + w a 3 f 2 , 2 ( X ) , e f 3 ( X ) = w a 1 f 4 , 2 ( X ) + w a 2 f 6 , 2 ( X ) , e f 4 ( X ) = w a 1 f 3 , 2 ( X ) + w a 2 f 5 , 2 ( X ) + w a 3 f 7 , 2 ( X ) . D. Relations to JPD A and LMB The TO-MB approximation was shown to be exactly the approximation made in the joint integrated probabilistic data association filter [44], a variant of the JPD A filter by incorpo- rating a probability of existence for each object. That is, the joint association distribution is approximated by the product of its marginals. The main difference is related to the modelling of undetected objects and objects detected for the first time; see [30] for more discussions. The PMBM filtering recursion can also handle an MB birth model by setting the PPP intensity for undetected objects to zero and adding new Bernoulli components in the prediction step, resulting in the MBM filter [32], [33]. Performing the TO-MB approximation after the update step of the MBM filter will then give a TO-MB filter . If the Bernoulli component is uniquely labelled in the MB birth, then the TO-MB filter can be considered as a fast implementation of the LMB filter without MBM 01 expansion [45], [46]. Compared to the MB birth model, the Poisson birth model is measurement-driven. The extended object LMB filter can also use an adaptive object birth model that is measurement-driven and allows the 2 The KLD on the space of sets of objects with auxiliary variables is an upper bound on the KLD for sets of objects without auxiliary variables [43]. 8 newborn objects to appear anywhere in the state space [19]. Howe ver , such adaptiv e birth model requires careful choice of a few design parameters which are application dependent. For the LMB filter, object trajectories can be formed by linking object state estimates with the same label. For the PMBM and the TO-PMB filters, object trajectories may be obtained based on filter meta-data (or auxiliary variables [43]). For filters based on PMBM conjugate prior, the optimal track building approach is by computing multi-object densities on sets of trajectories [47], which has led to, e.g., the trajectory PMBM filter [38], [48] and the trajectory PMB filter [43]. V . V A R I A T I O N A L M U LT I - B E R N O U L L I A P P RO X I M A T I O N The idea of using variational inference for obtaining an MB approximation that minimises the KLD from the true posterior was first proposed in [35] for point object filtering. This section presents the variational MB approximation of the MBM in the extended object PMBM posterior density resulted from the PMB update in Lemma 3 where each measurement creates a new Bernoulli component, i.e., n k | k = n k | k − 1 + m k . T wo dif ferent formulations of the minimisation objecti ve in variational approximations are studied: one is based on finding the most likely permutation of local hypotheses in each global hypothesis, and the other follo ws the linear programming formulation proposed in [35]. A. Appr oximate solution of KLD minimisation The objective is to find an MB of the form (17) with n k | k = n k | k − 1 + m k Bernoulli components that minimises the KLD argmin [ e f i ] Z f mbm k | k ( X k ) log f mbm k | k ( X k ) f mb k | k ( X k ) ! δ X k = argmin [ e f i ] − Z f mbm k | k ( X k ) log f mb k | k ( X k ) δ X k , (21) where f mbm k | k ( · ) is an MBM of the form (5c), f mb k | k ( · ) is an MB of the form (17) consisting of Bernoulli components e f i k | k ( · ) , and h e f i i denotes a collection of Bernoulli densities e f i k | k ( · ) with i ∈ { 1 , . . . , n k | k } . W e further let Π N denote the set of permutation functions on I N = { 1 , ..., N } : Π N , { π : I N → I N | i 6 = j ⇒ π ( i ) 6 = π ( j ) } . After simplifying the MB set integral in (21) into a series of Bernoulli set integrals, the solution of (21) is the same as the solution of the optimisation [35, Theorem 5] argmin [ e f i ] − X a ∈A k | k w a k | k Z · · · Z n k | k Y i =1 f i,a i k | k  X i  × log X π ∈ Π n k | k n k | k Y i =1 e f π ( i ) k | k  X π ( i )  δ X 1 · · · δ X n k | k , (22) which is suitable for applying expectation-maximisation (EM). Due to the set representation of the multi-object density , we should describe the correspondence between the underlying local hypothesis density f i,a i k | k ( · ) in global hypothesis a and the approximate Bernoulli component e f i k | k ( · ) using a missing data distribution q a ( π ) , which probabilistically determines which local hypothesis densities in f mbm k | k ( · ) should be merged to obtain each approximate Bernoulli component in f mb k | k ( · ) . Im- portantly , this relaxes the independence assumption of TO-MB approximation where only local hypothesis densities under the same Bernoulli component can be merged. By incorporating the missing data distribution q a ( π ) for each global hypothesis a and using the log-sum inequality , an upper bound of the objectiv e function in (22) is giv en by [35, Sec. III-A] J h e f i i , [ q a ( π )]  = − X a ∈A k | k ,π ∈ Π n k | k w a k | k q a ( π ) × − log q a ( π ) + n k | k X i =1 Z f i,a i k | k  X i  log e f π ( i ) k | k  X i  δ X i ! , (23) where [ q a ( π )] denotes a collection of missing data distrib ution q a ( · ) with a ∈ A k | k , q a ( π ) ≥ 0 ∀ a, π and P π ∈ Π n k | k q a ( π ) = 1 ∀ a . The optimisation problem (22) can then be approxi- mately solved by minimising the upper bound (23): min [ e f i ] , [ q a ( π )] J h e f i i , [ q a ( π )]  . (24) An approximate solution of (24) can be obtained using EM by coordinate descent 3 , which proceeds by alternating between estimating the missing data distribution q a ( π ) with a ∈ A k | k (E-step), and optimising e f i k | k ( · ) with i ∈ { 1 , . . . , n k | k } that minimises (23) (M-step). These two steps can be solved as q a ( π ) ∝ n k | k Y i =1 exp  Z f i,a i k | k ( X ) log e f π ( i ) k | k ( X ) δ X  , (25) e f i k | k ( X ) = X a ∈A k | k ,π ∈ Π n k | k w a q a ( π ) f π − 1 ( i ) ,a π − 1 ( i ) k | k ( X ) , (26) where π − 1 ( · ) is the in verse of the permutation function π ( · ) . B. F inding the most likely assignment Exact solution of the missing data distrib ution q a ( π ) in (25) may be computed by enumerating all the possible permutations of local hypothesis densities contained in global hypothesis a in a brute force fashion. Approximate solutions of (25) may also be obtained by only considering permutations with high weights using e.g., Murty’ s algorithm [49] or Gibbs sampling [50]. Empirically , we hav e found that it is more ef ficient from a computational perspectiv e to only consider the most likely permutation of local hypothesis densities for each global hypothesis. In this case, the missing data distrib ution under 3 In most applications of EM, the missing data distribution is estimated for each finite number of training samples. In turn, this guarantees that the log- sum inequality is tight at the optimum, and hence that the procedure will con ver ge to a local minimum of (24) [35]. 9 each global hypothesis becomes a point mass, and the E-step (25) and the M-step (26) reduce to argmin e π a − n k | k X i =1 Z f i,a i k | k ( X ) log e f e π a ( i ) k | k ( X ) δ X , (27) e f i k | k ( X ) = X a ∈A k | k w a f e π − 1 a ( i ) ,a e π − 1 a ( i ) k | k ( X ) , (28) where the most likely assignment e π a for each global hy- pothesis a ∈ A k | k can be obtained by solving a 2D assign- ment problem with computational complexity O  n 3 k | k  [51]. Hence, the total computational complexity of solving (27) for all the global hypotheses is O  |A k | k | n 3 k | k  . In practical implementations, dif ferent approximation methods can be used to keep the number of global hypotheses |A k | k | at a tractable lev el; this is discussed in Section VI-C. Because the MB distribution is in variant to the indices of the Bernoulli components it contains, the selection of the per- mutation function e π a for each global hypothesis a ∈ A k | k will not change the MBM f mbm k | k ( · ) , but it determines which local hypothesis densities should be merged. In TO-MB approx- imation, local hypothesis densities f i,a i k | k ( · ) corresponding to the same Bernoulli component, i.e., with the same superscript i , are merged. The simplified E-step (27) can be understood as constructing an alternative local hypothesis representation where local hypothesis densities f e π − 1 a ( i ) ,a e π − 1 a ( i ) k | k ( · ) with a ∈ A k | k correspond to the i -th Bernoulli component, whereas the simplified M-step (28) is equiv alent to applying TO-MB approximation to f mbm k | k ( · ) using the alternati ve local hypothesis representation determined by (27). This closely connects to the problem of searching for the local hypothesis representation for ne wly detected objects that yields the most accurate TO- MB approximation in the sense of KLD minimisation; the difference here is that the optimisation problem is now jointly solved for all the detected objects. It was also shown in [35] that the missing data distribution in the above formulation acts in an equi valent manner to the selection of an ordered distribution in the same unordered family that can be approximated by the desired distribution family in the Kullback–Leibler JPDA filter [52]. Howe ver , note that the related optimisation problem in [52] is solved using importance sampling, which is more computationally demanding than solving 2D assignment problems. C. Efficient appr oximation of the feasible set It was observed in [35] that the EM algorithm with the objectiv e function e J h e f i i , [ q a ( π )]  = − X a ∈A k | k ,π ∈ Π n k | k w a k | k q a ( π ) × n k | k X i =1 Z f i,a i k | k  X i  log e f π ( i ) k | k  X i  δ X i (29) tends to yield solution with lower KLD. For notational con- venience, let the local hypotheses f h k | k ( · ) for all the Bernoulli components be index ed through the set H k | k with |H k | k | = P n k | k i =1 h i k | k . The minimisation of (29) can be solved equiv a- lently as [35, Theorem 6] min q ( h,j ) ∈P − n k | k X j =1 Z   X h ∈H k | k q ( h, j ) f h k | k ( X )   log e f j k | k ( X ) δ X , (30) where the polytope P is P =      q ( h, j ) = n k | k X i =1 X a ∈A k | k | a i = h w a X π ∈ Π n k | k | π ( i )= j q a ( π )       q a ( π ) ≥ 0 , X π ∈ Π n k | k q a ( π ) = 1    , (31) a =  a 1 , . . . , a n k | k  and q ( h, j ) specifies the weight of local hypothesis f h k | k ( · ) in Bernoulli component e f j k | k ( · ) . An approximate solution of (30) can be obtained by con- sidering a relaxation of the polytope P M =    q ( h, j ) ≥ 0       X h ∈H k | k q ( h, j ) = 1 , n k | k X j =1 q ( h, j ) = n k | k X i =1 X a ∈A k | k | a i = h w a    . (32) Compared to (31), the missing data distribution in (32) is not constrained to depend on q a ( π ) . In this case, the E-step (25) and the M-step (26) reduce to argmin e q ( h,j ) ∈M − X h ∈H k | k n k | k X j =1 e q ( h, j ) Z f h k | k ( X ) log e f j k | k ( X ) δ X , (33) e f j k | k ( X ) = X h ∈H k | k e q ( h, j ) f h k | k ( X ) , (34) where the optimisation problem in volv ed in (33) is equiv alent to the minimum-cost flow problem that can be solved by linear programming in polynomial time [53, Chap. 7]. The computational complexity of solving (33) depends on the number of Bernoulli components n k | k and the number of local hypotheses H k | k , b ut not the number of global hypotheses |A k | k | , as opposed to the optimisation problem in (27). D. An illustrative example The v ariational MB approximation in volves an iterative update of the missing data distrib ution and the approximate MB. When initialised using the mar ginal association proba- bilities, the variational MB approximation may be considered an impro vement on the TO-MB approximation in the sense that an upper bound of the KLD (21) between the MBM and the approximate MB is minimised. W e illustrate how the missing data distributions change from one iteration of the variational MB approximation to the ne xt for the two different optimisation methods, respectiv ely , in the following example. 10 Example 4. Consider the same scenario shown in F ig. 1. Suppose that we have used the TO-MB appr oximation to obtain an appr oximate MB with four Bernoulli components, as described in Example 3. F or the variational MB approximation using 2D assign- ment, we can find the most likely permutation of Bernoulli components for each MB accor ding to (27). In this case, the sum of the cr oss entr opy between f 2 , 2 ( · ) and e f 1 ( · ) and the cr oss entropy between f 1 , 2 ( · ) and e f 2 ( · ) (cf. (27)) is smaller than the sum of the cross entropy between f 1 , 2 ( · ) and e f 1 ( · ) and the cr oss entr opy between f 2 , 2 ( · ) and e f 2 ( · ) . In order to minimise the objective (24), the or der of the first two Bernoulli components in the thir d MB should be swapped. After the r eordering , the first two approximate Bernoulli components is e f 1 ( X ) = ( w a 1 + w a 2 ) f 1 , 1 ( X ) + w a 3 f 2 , 2 ( X ) , e f 2 ( X ) = ( w a 1 + w a 2 ) f 2 , 1 ( X ) + w a 3 f 1 , 2 ( X ) . Under the same assumption, for the variational MB appr ox- imation using linear pr ogramming, in or der to minimise the objective (29), a pr oportion of the weight of assigning f 1 , 2 ( X ) to e f 1 ( X ) should be shifted to e f 2 ( X ) , and accor dingly a pr oportion of the weight of assigning f 2 , 2 ( X ) to e f 2 ( X ) should be shifted to e f 1 ( X ) . F or example, the first two appr oximate Bernoulli components after weight shifting may become e f 1 ( X ) = ( w a 1 + w a 2 ) f 1 , 1 ( X ) + ( w a 3 − ∆ w ) f 2 , 2 ( X ) + ∆ w f 1 , 2 ( X ) , e f 2 ( X ) = ( w a 1 + w a 2 ) f 2 , 1 ( X ) + ( w a 3 − ∆ w ) f 1 , 2 ( X ) + ∆ w f 2 , 2 ( X ) , wher e both the weight of f 1 , 2 ( X ) in e f 1 ( X ) and the weight of f 2 , 2 ( X ) in e f 2 ( X ) are 0 ≤ ∆ w ≤ w a 3 , given by the optimal ˜ q ( · , · ) via solving (33). V I . G G I W I M P L E M E N TA T I O N Solving the multiple extended object filtering problem re- quires not only an MO T filtering recursion, but also a single extended object model. There are sev eral single extended object models av ailable in the literature, including, for ex- ample, the random matrix model [54]–[56], the multiplicative error model [57], the random hypersurface model [58] and its variant Gaussian process model [59]. In this paper, we hav e chosen the random matrix model in [55], in which the object shape is approximated as an ellipse. The random matrix model is simple to use, yet flexible enough to be applicable to many real scenarios, e.g., pedestrian tracking using laser scanner [17], [60] and ship tracking using marine radar [36], [61]. The random matrix model has also been used in sev eral multiple extended object filters, see, e.g., [17]–[20]. A. Single object model In the random matrix model, the e xtended object state x k is the combination of a vector ξ k describing the object kinematic state, a matrix χ k describing the object extent and a scalar γ k being the measurement model Poisson rate. The measurement likelihood for a single measurement z is φ ( z k | x k ) = N ( z k ; H ξ k , z χ k + R k ) , (35) where H selects the object position from ξ k , z is a scaling factor , and R k is the covariance of the Gaussian measurement noise. The single object density for the PPP model (5b) with single measurement likelihood (35) is a product of gamma, Gaussian and in verse-W ishart distributions [36], [55] p ( x ) = G AM ( γ ; a, b ) N ( ξ ; m, P ) I W d ( χ ; v , V ) , (36a) , G G I W ( x ; ζ ) , (36b) where G AM ( γ ; a, b ) is the gamma distribution defined on γ > 0 with shape a and rate b , I W d ( χ ; v , V ) is the in verse- W ishart distrib ution defined on positiv e definite d × d matrix with degrees of freedom v and d × d scale matrix V , and ζ = ( a, b, m, P , v , V ) is the GGIW density parameters. The motion models are giv en by ξ k = g ( ξ k − 1 ) + w k − 1 , (37a) χ k = M ( ξ k − 1 ) χ k − 1 M ( ξ k − 1 ) T , (37b) γ k = γ k − 1 , (37c) where g ( · ) is a kinematic motion model, w k − 1 is Gaussian process noise with zero mean and covariance Q , and M ( · ) is a transformation matrix. B. Bernoulli-GGIW distribution T o solve the optimisation problems (27) and (33), the cross entropy between two Bernoulli-GGIW distributions needs to be computed. Because the gamma distribution, the Gaussian distribution and the in verse-W ishart distribution in (36a) are mutually independent, analytical expression of such cross entropy exists; see Appendix D for details. T o approximate (20d), (28) and (34) as a single Bernoulli component, we also need expressions for merging a Bernoulli- GGIW mixture distribution. Suppose that we hav e a mixture of Bernoulli components index ed by h ∈ H and that the h -th Bernoulli component has weight w h , probability of existence r h and GGIW density G G I W ( x h ; ζ h ) . According to the result in Appendix B, the best-fitting Bernoulli-GGIW in the sense of KLD minimisation has probability of existence e r = X h ∈H w h r h , and GGIW parameters argmin e ζ D KL P h ∈H w h r h G G I W ( x h ; ζ h ) P h ∈H w h r h      G G I W ( x ; e ζ ) ! . The merging of multi variate Gaussian mixture that minimises the KLD can be achie ved using moment matching. The numer - ical methods for mer ging gamma mixture and inv erse-Wishart mixture via analytical KLD minimisation are presented in [40] and [41], respectiv ely . 11 C. Appr oximations for computational tractability The pseudo code for the prediction and the update of PPP- GGIW and Bernoulli-GGIW can be found in [20], so we focus on approximation methods for keeping the computational complexity of the GGIW -PMB filter at a tractable level. First, gating is performed to remov e unlikely measurement- to-object associations. In extended object filtering, the gates need to take into account both the position and extent of the object, as well as the state uncertainties. Second, clustering algorithms with different hyperparameter settings are used to obtain sev eral different partitions of the measurements. Last, for each measurement partition we only consider cluster-to- object assignments with relativ ely high likelihoods, which can be obtained using, e.g., ranked assignment [49] or Gibbs sampling [50]. The abov e clustering and assignment procedure usually results in a small subset of the cluster-to-object assign- ments with relati vely high likelihoods, and therefore greatly reduces the number of possible data associations in the update step. Alternativ e approaches to handling the data association problem in extended object filtering include the random sam- pling methods [60] and the sum-product algorithms [13], [14]. When performing MB approximation, it is computationally lighter to only consider MBs with non-negligible normalised weight. Thus, we prune MBs whose updated weight fall below a threshold before MBM merging. After MB approximation, Bernoulli components with probability of existence belo w a threshold are recycled and added to the updated PPP intensity representing undetected objects, see Section III-D. The PPP intensity is typically represented by a mixture. T o reduce the computational complexity of the measurement update of PPP , mixture components with weight below a threshold are pruned from the PPP intensity , and similar mixture components, in the sense of small KLD, can be merged. W e emphasise that the approximations are applied to the elements a of the outer sum in (5c) the elements of D u ( · ) in (5b). The summations mapping elements of X k on the LHS of (5a) to X u k and X d k , and mapping the elements of X d k on the LHS to X 1 k , . . . , X n k ensure symmetry of the respecti ve distributions regardless of the pruning and recycling approximations that we apply . These symmetrising summations do not need to be e xplicitly ev aluated in prediction, update or estimation; they remain implicit throughout and are not subject to approximation. V I I . S I M U L A T I O N S A N D R E S U LT S This section presents the results from a Monte Carlo simu- lation with 100 runs where the performance of eleven dif ferent extended object filtering implementations 4 are ev aluated: 1) PMBM filter with track-oriented implementation [20]. 2) MBM filter with track-oriented implementation. 3) TO-PMB filter with PMB update in Lemma 2 where each measurement cluster creates a ne w Bernoulli com- ponent, referred to as TO-PMB-C. 4) TO-PMB filter with PMB update in Lemma 3 where each measur ement creates a ne w Bernoulli component, referred to as TO-PMB-M. 4 MA TLAB code is available at https://github.com/yuhsuansia. 5) PMB filter with v ariational MB approximation based on 2D assignment (27), referred to as V -PMB-A. 6) PMB filter with v ariational MB approximation based on linear pr ogramming (33), referred to as V -PMB-LP . 7) TO-MB filter , see Section IV -D. 8) MB filter with variational MB approximation based on 2D assignment , referred to as V -MB-A. 9) MB filter with variational MB approximation based on linear pr ogramming , referred to as V -MB-LP . 10) LMB filter [19] with joint prediction and update [62]. 11) LMB filter with adaptive birth model [19] as well as joint prediction and update, referred to as LMB-AB. W e ev aluate these implementations in three different simulated scenarios, with true object trajectories and extents illustrated in Fig. 3. For all the scenarios, the object’ s surviv al probability is p S = 0 . 99 , the object’ s Poisson measurement rate is γ = 10 and its initial extent is randomly generated from an in verse- W ishart distribution with mean χ = 4 I 2 where I 2 is an identity matrix with size 2. In the first scenario, 96 randomly generated objects are born around four locations, and they appear in and disappear from the surveillance area [ − 200 , 200] × [ − 200 , 200] at different time steps. The probability of detection is p D = 0 . 9 , and the clutter is uniformly distrib uted in the surveillance area with Poisson rate λ = 60 . In the second scenario, 26 objects first get close to each other and then separate. The probability of detection is p D = 0 . 3 and the Poisson clutter rate is λ = 100 . In the third scenario, 8 objects are all born around the same location at time step 1 . The probability of detection is p D = 0 . 7 and the Poisson clutter rate is λ = 30 . The ellipsoidal gate size in probability is 0 . 999 . For solving the data association problem, we first apply DBSCAN [63] with different hyperparameters to obtain a set of measurement partitions. DBSCAN is also used to identify the clusters of new Bernoulli components in PMBM and TO-PMB. W e then use Murty’ s algorithm 5 to find the d 20 · w a e best cluster- to-object assignments for each measurement partition. The updated MBs are pruned to only contain high-weight MBs that correspond to 99% of the likelihood. For filters without MB approximation, the number of updated MBs is also capped at 100 . For filters with MB birth model, Bernoulli components with probability of existence smaller than 0 . 001 are pruned. For filters with Poisson birth model, Bernoulli components with probability of existence smaller than 0 . 1 are recycled, and PPP mixture components with weight smaller than 0 . 001 are pruned. For filters with v ariational MB approximation, the con vergence threshold for the EM algorithm is 0 . 1 . For V - PMB-A and V -MB-A, the 2D assignment problem is solved using a modified Jonker-V olgenant algorithm [51]. For V - PMB-LP and V -MB-LP , the linear programming is addressed using Gurobi optimiser [65]. The kinematic state is ξ k = [ p k , v k ] T , and it describes the object’ s position p k ∈ R 2 and velocity v k ∈ R 2 . The sym- metric positi ve definite random matrix χ k is tw o-dimensional. 5 The C++ implementation in the T racker Component Library [64] is used. 12 (a) Scenario with 96 objects (b) Scenario with 26 objects (c) Scenario with 8 objects Fig. 3. True object trajectories for the simulated scenarios. The object extent is shown for every 10 time steps. In the first scenario (100 time steps), objects are born around four possible locations. In the second scenario (101 time steps), at most one object is born per time step, and objects move in proximity around the origin. In the third scenario (50 time steps), objects are all born around the same location at time step 1. Objects mo ve with a nearly constant velocity model (37) with g ( ξ k ) =  I 2 T s I 2 0 2 I 2  ξ k , Q = q 2 " T 3 s 3 I 2 T 2 s 2 I 2 T 2 s 2 I 2 T s I 2 # where T s = 1 is the time interval and q = 0 . 3 is the standard acceleration de viation. The transformation matrix is M ( ξ k ) = I 2 . For the measurement model (35), the scaling factor is z = 0 . 25 and the measurement noise cov ariance is R k = 0 . 5 I 2 . All the filters extract state estimates by first determining the maximum a posteriori estimate of the cardinality . For filters with MB approximation, mean values of the object states are reported from the corresponding number of Bernoulli compo- nents with the highest probability of existence. F or PMBM and MBM, mean values of the object states are reported from the global hypothesis with corresponing cardinality with highest weight in a similar manner , see [32, Sec. VI.B] for details. For performance ev aluation of extended object estimates with ellipsoidal extent, a comparison study in [66] has shown that a reasonable choice is the Gaussian W asserstein distance (GWD) metric, which is defined as [67] d GW ( x , ˆ x ) =    H  ξ − ˆ ξ     2 + T r  χ + ˆ χ − 2  χ 1 2 ˆ χχ 1 2  1 2  , where Tr ( X ) is the trace of matrix X . T o ev aluate the multi- object filtering performance, the GWD metric is integrated into the generalised optimal sub-pattern assignment (GOSP A) metric [68] with parameters α = 2 , c = 60 and p = 1 , which allows for the decomposition of the estimation error into localisation errors and costs for missed and false objects. For the first scenario, the Poisson birth intensity is a GGIW mixture with four components, and each of them has weight 0 . 25 and Gaussian density with mean [ ± 50 , ± 50 , 0 , 0] T and cov ariance diag (100 , 100 , 4 , 4) cov ering a potential birth area, whereas the MB birth density has four Bernoulli component, and each of them has probability of existence 0 . 25 and single object GGIW density covering a potential birth area. For the second scenario, the Poisson birth intensity in this case is a single GGIW component with weight 0 . 25 and Gaussian density with mean [ − 100 , − 100 , 0 , 0] T and cov ari- ance diag (10000 , 10000 , 16 , 16) cov ering a large surveillance area, whereas the MB birth density only contains a single Fig. 4. Percentage of time steps for which there exists at least one measurement cluster that simultaneously falls within the ellipsoidal gates of at least n Bernoulli or PPP components for PMBM, V -PMB-LP and LMB in three different scenarios. Bernoulli with probability of existence 0 . 25 and the same GGIW component as it is in the Poisson birth intensity . For birth models used in the third scenario, the only dif ference compared to the second scenario is that the Gaussian density now has mean [0 , 0 , 0 , 0] T and cov ariance diag (16 , 16 , 4 , 4) . For a quantitativ e analysis of the data association uncer- tainty in the three different scenarios, we also computed the percentage of time steps for which there exists at least one measurement cluster that simultaneously falls within the ellipsoidal gates of at least n Bernoulli or PPP components, see Fig. 4 for the results of PMBM, V -PMB-LP and LMB. As can be seen, V -PMB-LP has the lo west data association uncertainty as it maintains fewer Bernoulli components at each time step. For the three dif ferent scenarios, the data association problem is most challenging in the scenario with 26 objects where many objects move in proximity . The mean GOSP A error and its decomposition, along with the mean cycle time per time step 6 , for the ele ven different filters in two different scenarios are presented in T able I. Also, the mean GOSP A error and its decomposition for the PMBM, MBM, T O-PMB-M, V -PMB-LP and LMB filters in three 6 MA TLAB implementation on 3 GHz 6-Core Intel Core i5. 13 T ABLE I M E AN G O S P A E R RO R A N D C Y CL E T I M E ( S E CO N D S ) P E R T I M E S T EP , A V E R AG ED OV E R 1 0 0 M O N TE C A R LO R UN S . L E GE N D : T O L .– G O S P A E R RO R ; L O C . – L O CA L I S A T I ON E R RO R ; M I S .– M I S SE D O B J E CT E R RO R ; F A L . – FA LS E O B JE C T E R RO R ; C Y C .– C Y C LE T I M E Scenario with 96 objects Scenario with 26 objects Scenario with 8 objects Filter T O L . L OC . M I S . F A L . C Y C . T O L . L O C . M I S . F AL . C YC . T OL . L O C . M I S . F A L . C Y C . TO - P MB - C 75.38 35.68 19.79 19.91 1.04 156.99 48.12 85.98 22.88 0.69 21.62 5.76 9.77 6.08 0.18 TO - P MB - M 75.37 35.69 19.73 19.95 0.73 156.71 48.25 85.86 22.60 0.57 21.36 5.80 9.57 5.99 0.14 V - P M B- A 75.30 35.68 19.65 19.98 0.82 152.58 47.63 81.06 23.90 0.64 21.09 5.55 9.53 6.02 0.15 V - P M B- L P 75.30 35.67 19.65 19.98 0.81 152.16 47.30 80.68 24.17 0.62 21.11 5.53 9.54 6.04 0.15 P MB M 74.89 35.61 19.31 19.97 2.74 149.02 57.62 60.00 31.40 1.81 20.42 5.29 9.49 5.65 0.40 TO - M B 82.83 41.48 21.11 20.23 1.17 188.53 61.07 100.88 26.58 0.88 54.35 10.85 39.77 3.72 0.14 V - M B -A 76.89 36.86 19.03 21.00 1.47 168.06 56.61 80.91 30.54 1.00 51.22 10.49 36.72 4.01 0.14 V - M B -L P 76.71 36.72 18.98 21.01 1.35 167.42 56.01 81.21 30.20 0.93 51.19 10.53 36.64 4.03 0.14 M BM 75.56 35.98 18.67 20.91 4.64 169.86 72.27 54.54 43.04 3.20 50.17 10.30 35.58 4.29 0.37 L MB 83.03 40.42 19.19 23.42 2.07 249.15 64.81 98.37 85.97 0.79 57.48 11.13 42.73 3.61 0.23 L MB - A B 165.39 56.19 77.03 32.17 1.55 265.74 65.70 145.23 54.81 0.22 50.23 15.03 29.83 5.38 0.13 different scenarios are shown in Fig. 5, 6 and 7, respectiv ely . It can be seen that the PMBM filter has the best estimation performance for all the scenarios. LMB-AB shows the worst performance in the first scenario as it does not know the region where new objects will appear . LMB-AB outperforms LMB in the third scenario where there is a birth model mismatch. All filters with Poisson birth model outperform their coun- terparts with MB birth model in terms of both estimation error and computational time. This is mainly because the PMBM conjugate prior has a more efficient hypothesis structure than the MBM conjugate prior [32]. TO-PMB-C has slightly better estimation performance than TO-PMB-M, and it is also faster . The noticable difference in mean cycle time is because in T O- PMB-C many more new Bernoulli components are created in the update step, which increases the complexity of the PMB filtering recursion in subsequent time steps. All filters with v ariational MB approximation outperform their counterparts with TO-MB approximation and the differ - ence is most noticable in the second scenario with coalescence. The multi-object density becomes strongly multimodal when objects are closely-spaced, and in this case it is less accu- rate to consider the TO-MB approximation. The PMB filters with v ariational MB approximation not only ha ve very close estimation performance to the PMBM filter in all the three scenarios, b ut they are also significantly faster . For the two variants of the variational MB approximation, V -PMB-LP has similar performance compared to V -PMB-A, whereas V -MB- LP has slightly better performance than V -MB-A. V I I I . C O N C L U S I O N S W e have proposed two different PMB approximations of the multi-object posterior density for e xtended object filtering. They are both based on a new and more efficient local hy- pothesis representation for newly detected objects. T wo PMB filters and their GGIW implementations are presented: the TO- PMB filter is based on the TO-MB approximation and the V -PMB filter is based on the variational MB approximation. The PMB filters offer a good comprise between computational complexity and accuracy . Specifically , they are generally faster than the PMBM filter with slightly w orse performance. The two V -PMB filters hav e better performance than the T O-PMB filter , especially in scenarios with a high risk of coalescence. Fig. 5. GOSP A error and its decomposition versus time for PMBM, TO- PMB-M, V -PMB-LP , MBM and LMB in scenario with 96 objects. Fig. 6. GOSP A error and its decomposition versus time for PMBM, TO- PMB-M, V -PMB-LP , MBM and LMB in scenario with 26 objects. R E F E R E N C E S [1] S. Blackman and R. Popoli, Design and analysis of modern trac king systems . Norwood, MA: Artech House, 1999. 14 Fig. 7. GOSP A error and its decomposition versus time for PMBM, TO- PMB-M, V -PMB-LP , MBM and LMB in scenario with 8 objects. [2] S. Challa, M. R. Morelande, D. Mu ˇ sicki, and R. J. Evans, Fundamentals of object tracking . Cambridge University Press, 2011. [3] F . Meyer , T . Kropfreiter, J. L. W illiams, R. Lau, F . Hlawatsch, P . Braca, and M. Z. Win, “Message passing algorithms for scalable multitarget tracking, ” Pr oceedings of the IEEE , vol. 106, no. 2, pp. 221–259, 2018. [4] K. Granstr ¨ om, M. Baum, and S. Reuter , “Extended object tracking: Intro- duction, overvie w and applications, ” Journal of Advances in Information Fusion , vol. 12, no. 2, 2017. [5] K. Gilholm and D. Salmond, “Spatial distribution model for tracking extended objects, ” IEE Pr oceedings-Radar , Sonar and Navigation , vol. 152, no. 5, pp. 364–371, 2005. [6] L. Hammarstrand, L. Svensson, F . Sandblom, and J. Sorstedt, “Extended object tracking using a radar resolution model, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 48, no. 3, pp. 2371–2386, 2012. [7] B. Habtemariam, R. Tharmarasa, T . Thayaparan, M. Mallick, and T . Kirubarajan, “ A multiple-detection joint probabilistic data association filter , ” IEEE Journal of Selected T opics in Signal Processing , vol. 7, no. 3, pp. 461–471, 2013. [8] R. Streit, “JPD A intensity filter for tracking multiple extended objects in clutter , ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2016, pp. 1477–1484. [9] G. V iv one and P . Braca, “Joint probabilistic data association tracker for extended target tracking applied to X-band marine radar data, ” IEEE Journal of Oceanic Engineering , v ol. 41, no. 4, pp. 1007–1019, 2016. [10] S. Y ang, K. Thormann, and M. Baum, “Linear-time joint probabilistic data association for multiple extended object tracking, ” in IEEE 10th Sensor Array and Multichannel Signal Pr ocessing W orkshop (SAM) . IEEE, 2018, pp. 6–10. [11] M. Wieneke and W . K och, “ A PMHT approach for extended objects and object groups, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 48, no. 3, pp. 2349–2370, 2012. [12] S. P . Coraluppi and C. A. Carthel, “Multiple-hypothesis tracking for targets producing multiple measurements, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 54, no. 3, pp. 1485–1498, 2018. [13] F . Meyer and M. Z. W in, “Scalable data association for extended object tracking, ” IEEE T ransactions on Signal and Information Pr ocessing over Networks , vol. 6, pp. 491–507, 2020. [14] F . Meyer and J. L. Williams, “Scalable detection and tracking of extended objects, ” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) . IEEE, 2020, pp. 8916–8920. [15] R. Mahler , “PHD filters for nonstandard targets, I: Extended targets, ” in 12th International Confer ence on Information Fusion . IEEE, 2009, pp. 915–921. [16] K. Granstr ¨ om, C. Lundquist, and O. Orguner , “Extended target tracking using a Gaussian-mixture PHD filter, ” IEEE T ransactions on Aerospace and Electr onic Systems , vol. 48, no. 4, pp. 3268–3286, 2012. [17] K. Granstr ¨ om and U. Orguner , “ A PHD filter for tracking multiple extended targets using random matrices, ” IEEE T ransactions on Signal Pr ocessing , v ol. 60, no. 11, pp. 5657–5671, 2012. [18] C. Lundquist, K. Granstr ¨ om, and U. Or guner , “ An extended target CPHD filter and a gamma Gaussian in verse Wishart implementation, ” IEEE Journal of Selected T opics in Signal Pr ocessing , vol. 7, no. 3, pp. 472– 483, 2013. [19] M. Beard, S. Reuter , K. Granstr ¨ om, B.-T . V o, B.-N. V o, and A. Scheel, “Multiple extended target tracking with labeled random finite sets, ” IEEE T ransactions on Signal Processing , vol. 64, no. 7, pp. 1638–1653, 2016. [20] K. Granstr ¨ om, M. Fatemi, and L. Svensson, “Poisson multi-Bernoulli mixture conjugate prior for multiple extended target filtering, ” IEEE T ransactions on Aerospace and Electr onic Systems , v ol. 56, no. 1, pp. 208–225, 2020. [21] A. Swain and D. Clark, “Extended object filtering using spatial indepen- dent cluster processes, ” in 13th International Conference on Information Fusion . IEEE, 2010, pp. 1–8. [22] R. P . Mahler , “Detecting, tracking, and classifying group targets: a unified approach, ” in Signal Pr ocessing, Sensor Fusion, and T arget Recognition X , vol. 4380. International Society for Optics and Photonics, 2001, pp. 217–228. [23] ——, “Extended first-order Bayes filter for force aggregation, ” in Signal and Data Pr ocessing of Small T ar gets 2002 , v ol. 4728. International Society for Optics and Photonics, 2002, pp. 196–207. [24] A. Swain and D. E. Clark, “First-moment filters for spatial independent cluster processes, ” in Signal Pr ocessing, Sensor Fusion, and T ar get Recognition XIX . International Society for Optics and Photonics, 2010, p. 76970I. [25] A. Swain and D. Clark, “Bayesian estimation of the intensity for independent cluster point processes: An analytic solution, ” Pr ocedia En vir onmental Sciences , vol. 7, pp. 56–61, 2011. [26] ——, “The PHD filter for extended target tracking with estimable extent shape parameters of varying size, ” in 15th International Conference on Information Fusion . IEEE, 2012, pp. 1111–1118. [27] D. Clark and R. Mahler, “Generalized PHD filters via a general chain rule, ” in 15th International Conference on Information Fusion . IEEE, 2012, pp. 157–164. [28] A. J. Swain, “Group and extended target tracking with the probability hypothesis density filter , ” Ph.D. dissertation, Heriot-W att University , 2013. [29] ´ A. F . Garc ´ ıa-Fern ´ andez, J. L. W illiams, L. Svensson, and Y . Xia, “ A Poisson multi-Bernoulli mixture filter for coexisting point and extended targets, ” IEEE T ransactions on Signal Processing , 2021. [30] J. L. W illiams, “Marginal multi-Bernoulli filters: RFS derivation of MHT, JIPDA, and association-based member , ” IEEE Tr ansactions on Aer ospace and Electr onic Systems , vol. 51, no. 3, pp. 1664–1687, 2015. [31] K. Granstr ¨ om, M. Fatemi, and L. Svensson, “Gamma Gaussian inv erse- Wishart Poisson multi-Bernoulli filter for extended target tracking, ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2016, pp. 893–900. [32] ´ A. F . Garc ´ ıa-Fern ´ andez, J. L. Williams, K. Granstr ¨ om, and L. Svensson, “Poisson multi-Bernoulli mixture filter: direct derivation and imple- mentation, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 54, no. 4, pp. 1883–1901, 2018. [33] ´ A. F . Garc ´ ıa-Fern ´ andez, Y . Xia, K. Granstr ¨ om, L. Sv ensson, and J. L. W illiams, “Gaussian implementation of the multi-Bernoulli mixture fil- ter , ” in 22th International Conference on Information Fusion (FUSION) . IEEE, 2019, pp. 1–8. [34] B.-T . V o and B.-N. V o, “Labeled random finite sets and multi-object conjugate priors, ” IEEE T ransactions on Signal Processing , vol. 61, no. 13, pp. 3460–3475, 2013. [35] J. L. Williams, “ An efficient, variational approximation of the best fitting multi-Bernoulli filter, ” IEEE Tr ansactions on Signal Pr ocessing , vol. 63, no. 1, pp. 258–273, 2015. [36] K. Granstr ¨ om, A. Natale, P . Braca, G. Ludeno, and F . Serafino, “Gamma Gaussian inverse Wishart probability hypothesis density for extended target tracking using X-band marine radar data, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 53, no. 12, pp. 6617–6631, 2015. [37] R. P . Mahler, Statistical Multisour ce-Multitar get Information Fusion . Artech House Norwood, MA, 2007. [38] Y . Xia, K. Granstr ¨ om, L. Sv ensson, ´ A. F . Garc ´ ıa-Fern ´ andez, and J. L. W illiams, “Extended target Poisson multi-Bernoulli mixture trackers based on sets of trajectories, ” in 22th International Confer ence on Information Fusion (FUSION) . IEEE, 2019, pp. 1–8. [39] J. L. Williams, “Hybrid Poisson and multi-Bernoulli filters, ” in Pr oceed- ings of International Conference on Information Fusion . IEEE, 2012, pp. 1103–1110. [40] K. Granstr ¨ om and U. Orguner , “Estimation and maintenance of mea- surement rates for multiple extended target tracking, ” in Proceedings 15 of International Confer ence on Information Fusion . IEEE, 2012, pp. 2170–2176. [41] ——, “On the reduction of Gaussian in verse Wishart mixtures, ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2012, pp. 2162–2169. [42] R. P . Mahler , Advances in Statistical Multisource-Multitar get Informa- tion Fusion . Artech House Norwood, MA, 2014. [43] ´ A. F . Garc ´ ıa-Fern ´ andez, L. Svensson, J. L. W illiams, Y . Xia, and K. Granstr ¨ om, “Trajectory Poisson multi-Bernoulli filters, ” IEEE T rans- actions on Signal Processing , vol. 68, pp. 4933–4945, 2020. [44] D. Musicki and R. Evans, “Joint integrated probabilistic data associa- tion: JIPD A, ” IEEE transactions on Aerospace and Electronic Systems , vol. 40, no. 3, pp. 1093–1099, 2004. [45] J. Olofsson, C. V eib ¨ ack, and G. Hendeby , “Sea ice tracking with a spatially indexed labeled multi-Bernoulli filter , ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2017, pp. 1–8. [46] T . Kropfreiter , F . Meyer , and F . Hlawatsch, “ A fast labeled multi- Bernoulli filter using belief propagation, ” IEEE Tr ansactions on Aer ospace and Electr onic Systems , vol. 56, no. 3, pp. 2478–2488, 2020. [47] ´ A. F . Garc ´ ıa-Fern ´ andez, L. Svensson, and M. R. Morelande, “Multiple target tracking based on sets of trajectories, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 56, no. 3, pp. 1685–1707, 2020. [48] K. Granstr ¨ om, L. Svensson, Y . Xia, J. Williams, and ´ A. F . Garc ´ ıa-Fem ´ andez, “Poisson multi-Bernoulli mixture trackers: continuity through random finite sets of trajectories, ” in Pr oceedings of Interna- tional Conference on Information Fusion . IEEE, 2018, pp. 1–5. [49] M. L. Miller , H. S. Stone, and I. J. Cox, “Optimizing Murty’ s ranked assignment method, ” IEEE T ransactions on Aer ospace and Electronic Systems , vol. 33, no. 3, pp. 851–862, 1997. [50] M. Fatemi, K. Granstr ¨ om, L. Svensson, F . J. Ruiz, and L. Hammarstrand, “Poisson multi-Bernoulli mapping using Gibbs sampling, ” IEEE T rans- actions on Signal Processing , vol. 65, no. 11, pp. 2814–2827, 2017. [51] D. F . Crouse, “On implementing 2D rectangular assignment algorithms, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 52, no. 4, pp. 1679–1696, 2016. [52] L. Svensson, D. Svensson, M. Guerriero, and P . Willett, “Set JPD A filter for multitarget tracking, ” IEEE T ransactions on Signal Pr ocessing , vol. 59, no. 10, pp. 4677–4691, 2011. [53] C. H. Papadimitriou and K. Steiglitz, Combinatorial optimization: algorithms and complexity . Courier Corporation, 1998. [54] J. W . Koch, “Bayesian approach to extended object and cluster tracking using random matrices, ” IEEE T ransactions on Aer ospace and Electr onic Systems , vol. 44, no. 3, 2008. [55] M. Feldmann, D. Franken, and W . K och, “Tracking of extended objects and group targets using random matrices, ” IEEE T ransactions on Signal Pr ocessing , v ol. 59, no. 4, pp. 1409–1420, 2011. [56] B. T uncer and E. Ozkan, “Random matrix based extended target tracking with orientation: A new model and inference, ” IEEE T ransactions on Signal Pr ocessing , 2021. [57] S. Y ang and M. Baum, “Tracking the orientation and axes lengths of an elliptical extended object, ” IEEE T ransactions on Signal Pr ocessing , vol. 67, no. 18, pp. 4720–4729, 2019. [58] M. Baum and U. D. Hanebeck, “Extended object tracking with random hypersurface models, ” IEEE T ransactions on Aerospace and Electronic systems , vol. 50, no. 1, pp. 149–159, 2014. [59] N. W ahlstr ¨ om and E. ¨ Ozkan, “Extended target tracking using Gaussian processes, ” IEEE T ransactions on Signal Pr ocessing , vol. 63, no. 16, pp. 4165–4178, 2015. [60] K. Granstr ¨ om, L. Svensson, S. Reuter , Y . Xia, and M. Fatemi, “Likelihood-based data association for extended object tracking using sampling methods, ” IEEE Tr ansactions on Intelligent V ehicles , vol. 3, no. 1, pp. 30–45, 2018. [61] M. Schuster , J. Reuter, and G. W anielik, “Probabilistic data association for tracking extended targets under clutter using random matrices, ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2015, pp. 961–968. [62] B.-N. V o, B.-T . V o, and H. G. Hoang, “ An efficient implementation of the generalized labeled multi-Bernoulli filter , ” IEEE T ransactions on Signal Pr ocessing , vol. 65, no. 8, pp. 1975–1987, 2016. [63] M. Ester , H.-P . Kriegel, J. Sander , X. Xu et al. , “ A density-based algorithm for discovering clusters in large spatial databases with noise. ” in Kdd , vol. 96, no. 34, 1996, pp. 226–231. [64] D. F . Crouse, “The tracker component library: free routines for rapid prototyping, ” IEEE Aer ospace and Electr onic Systems Magazine , v ol. 32, no. 5, pp. 18–27, 2017. [65] L. Gurobi Optimization, “Gurobi optimizer reference manual, ” 2020. [Online]. A vailable: http://www .gurobi.com [66] S. Y ang, M. Baum, and K. Granstr ¨ om, “Metrics for performance ev aluation of elliptic extended object tracking methods, ” in Proceedings of International Confer ence on Multisensor Fusion and Integr ation for Intelligent Systems (MFI) . IEEE, 2016, pp. 523–528. [67] C. R. Gi vens, R. M. Shortt et al. , “ A class of Wasserstein metrics for probability distributions. ” The Michigan Mathematical Journal , vol. 31, no. 2, pp. 231–240, 1984. [68] A. S. Rahmathullah, ´ A. F . Garc ´ ıa-Fern ´ andez, and L. Svensson, “Gen- eralized optimal sub-pattern assignment metric, ” in Pr oceedings of International Conference on Information Fusion . IEEE, 2017, pp. 1–8. 16 Supplementary Materials A P P E N D I X A A. Pr oof of Theorem 1 Because the dif ferent local hypothesis representations only differ in how the ne wly detected objects are represented, for simplicity we only demonstrate the case with only newly detected objects, i.e., n k | k − 1 = 0 . Using the local hypothesis representation in Lemma 2, the update of a PPP prior density with a set Z k =  z 1 k , . . . , z m k k  of measurements is a PMBM where the MBM density describ- ing the set X k of newly detected objects is f k | k ( X k ) = X a ∈A k | k w a k | k X U 2 m k − 1 l =1 X l = X k 2 m k − 1 Y i =1 f i,a i k | k  X i  (38) where the set A k | k = n a 1 , . . . , a 2 m k − 1  : a i ∈ { 1 , 2 } ∀ i, 2 m k − 1 ] i =1 M i,a i k = M k ) (39) of global hypotheses is made up of local hypotheses (16). By construction, there is a one-to-one correspondence between the set of partitions of M k and the set A k | k . Now consider a local hypothesis representation with n k ≥ m k Bernoulli components, and the MBM density describing the set X k of newly detected objects is instead expressed as f 0 k | k ( X k ) = X b ∈B k | k w b k | k X U n k l =1 X l = X k n k Y i =1 f i,b i k | k  X i  (40) where the set B k | k =  b 1 , . . . , b n k  : b i ∈ n 1 , . . . , h i k | k o ∀ i, n k ] i =1 M i,b i k = M k ) (41) of global hypotheses is made up of local hypotheses whose associations M i,j k satisfy (18). Note that for a local hypothesis with M i,j k its weight w i,j k | k and density f i,j k | k ( · ) are computed in the same manner (cf. (16)) in both representations (38) and (40). The objecti ve is to prove that (38) and (40) represent the same MBM density . W e use the follo wing two lemmas, which are proved later in Appendix A-B and A-C, to proceed with the proof. Lemma 4. If global hypotheses a and b repr esent the same partition of M k , the two corr esponding weights and MB densities given by a ∈ A k | k in (38) and b ∈ B k | k in (40), r espectively , ar e identical. That is, it holds that w a k | k = w b k | k , (42a) X U 2 m k − 1 l =1 X l = X k 2 m k − 1 Y i =1 f i,a i k | k  X i  = X U n k l =1 X l = X k n k Y ι =1 f ι,b ι k | k ( X ι ) (42b) wher e the global hypothesis weights w a k | k and w b k | k can be computed according to (6) using the local hypothesis weights n w i,a i k | k o and n w ι,b ι k | k o , r espectively . Lemma 5. F or a local hypothesis repr esentation with n k ≥ m k Bernoulli components satisfying (18), ther e is a one-to-one corr espondence between its corr esponding set B k | k of global hypotheses, given by (41), and the set of partitions of M k . Since there is also a one-to-one correspondence between the set of partitions of M k and the set A k | k , for ev ery a ∈ A k | k Lemma 5 ensures that there is precisely one b ∈ B k | k that corresponds to the same partition of M k . According to Lemma 4, we can conclude that (38) and (40) represent the same MBM density . This completes the proof. B. Pr oof of Lemma 4 The constraints (39) and (41), on the global hypotheses, ensure that each global hypothesis is a collection of local hypotheses, one for each Bernoulli component, and that the set of local hypotheses with non-empty M i,j k in a global hypothesis corresponds to a partition of M k . Therefore, if global hypotheses a and b represent the same partition of M k , their corresponding sets of local hypotheses with non-empty M i,j k are the same, i.e., n M i,a i k , w i,a i k | k , f i,a i k | k ( · ) o M i,a i k 6 = ∅ = n M ι,b ι k , w ι,b ι k | k , f ι,b ι k | k ( · ) o M ι,b ι k 6 = ∅ . Local hypotheses associated with empty M i,j k hav e local hypothesis weight w i,j k | k = 1 and Bernoulli density f i,j k | k ( ∅ ) = 1 (cf. 16a). Therefore, adding these local hypothese to a global hypothesis does not change its weight and density . This means that, for global hypotheses a ∈ A k | k and b ∈ B k | k , as long as their corresponding sets of local hypotheses with non-empty M i,j k represent the same partition of M k , (42) holds, i.e., they hav e the same weight (42a) and MB density (42b). This completes the proof. C. Pr oof of Lemma 5 Let us first establish that there is at least one b ∈ B k | k for ev ery partition of M k . T o this end, given an arbitrary partition P of M k , we will now show that there is a corresponding b ∈ B k | k . Constraint (18a) means that the local hypothesis representation contains all subsets in P . This ensures that for each element of P we can find an identical set M i,j k . Constraint (18e) means that two different Bernoulli compo- nents cannot ha ve local hypotheses whose M i,j k hav e non- empty intersection. This ensures that two elements in P can nev er be associated to local hypotheses of the same Bernoulli component. Finally , (18b) ensures that for ev ery Bernoulli component i , there is a local hypothesis with empty M i,j k . T o construct a global hypothesis b ∈ B k | k that corresponds to P , we can first go through the elements of P , and, for each element, we find an identical set M i,j k and set b i = j . 17 T ABLE II G G IW P R E DI C T I ON Input: GGIW parameter ζ , motion model g ( · ) , process noise cov ariance Q , transformation matrix M ( · ) , time interval T s , maneuvering correlation constant τ , measurement rate parameter η , dimension d of the extent state Output: Predicted GGIW parameter ζ + . ζ + =                  a + = a η , b + = b η , m + = g ( m ) , P + = GP G T + Q, v + = 2 d + 2 + e − T s /τ ( v − 2 d − 2) , V + = e − T s /τ M ( m ) V M ( m ) T , where G = ∇ ξ g ( ξ ) | ξ = m . At last, for each of the remaining Bernoulli components, we find its local hypothesis with empty M i,j k and set b i = j . T o summarise, constraints (18a), (18b) and (18e) together ensure that there is at least one b ∈ B k | k for ev ery partition of M k . T o complete the proof, it remains to establish that no par- tition P is represented by two different b ∈ B k | k . Constraints (18c), (18d) imply that each Bernoulli component has only one local hypothesis with empty M i,j k and that there does not exist two local hypotheses with the same non-empty subset of M k . Therefore, the uniqueness of each global hypothesis b ∈ B k | k is ensured. This completes the proof. A P P E N D I X B In this appendix, we present the mer ging of a mixture of Bernoulli densities, in the sense of KLD minimisation. Let f ( X ) = P h ∈H w h f h ( X ) be a mixture of Bernoulli components, where the single object density of Bernoulli component f h ( X ) is from a family of distributions F , i.e., p h ( x ) ∈ F ∀ h ∈ H . The objecti ve is then to approximate the Bernoulli mixture f ( X ) by a single Bernoulli distribution e f ( X ) , whose single object density is from the the same family of distributions, i.e., e p ( x ) ∈ F . The approximate Bernoulli density e f ( X ) that minimises the KLD D KL  f ( X ) k e f ( X )  has probability of existence and single object density [30]: e r = X h ∈H w h r h , (43a) e p ( x ) = argmin p ( x ) ∈F D KL P h ∈H w h r h p h ( x ) P h ∈H w h r h      p ( x ) ! . (43b) For distributions from the exponential family , the KLD min- imisation problem (43b) can be analytically solved by match- ing the expected sufficient statistics. A P P E N D I X C In this appendix, we present the pseudo code of factorised GGIW prediction and update with non-negligible measurement noise co variance [20], [36], see T able II and T able III, respec- tiv ely . T ABLE III G G IW U P DATE Input: GGIW parameter ζ + , set W of detections, measurement model H , measurement noise cov ariance R , scaling factor z , dimension d of the extent state Output: Updated GGIW parameter ζ and predicted likelihood ` . ζ =                  a = a + + | W | , b = b + + 1 , m = m + + K , P = P + − K H P + , v = v + + | W | , V = V + + N + ˆ Z , ` = π − d | W | 2 | W | − d 2 | V + | v + − d − 1 2 Γ d  v − d − 1 2     ˆ X    | W | 2 Γ( a )( b + ) a + | V | v − d − 1 2 Γ d  v + − d − 1 2     ˆ R    | W |− 1 2 | S | 1 2 Γ( a + )( b ) a where ¯ z = 1 | W | X z i ∈ W z i , Z = X z i ∈ W  z i − ¯ z   z i − ¯ z  T , ˆ X = V + ( v + − 2 d − 2) − 1 ,  = ¯ z − H m + , ˆ R = z ˆ X + R, S = H P + H T + ˆ R | W | , K = P + H T S − 1 , N = ˆ X 1 / 2 ˆ S − 1 / 2  T ˆ S − 1 / 2 ˆ X T / 2 , ˆ Z = ˆ X 1 / 2 ˆ R − 1 / 2 Z ˆ R − 1 / 2 ˆ X T / 2 . A P P E N D I X D In this appendix, we present the cross entropy between two Bernoulli-GGIWs. Suppose f h ( · ) and e f i ( · ) are two Bernoulli densities with the following form f h ( X ) = ( 1 − r h X = ∅ , r h G G I W  x h ; ζ h  X = { x } , (44a) e f i ( X ) = ( 1 − r i X = ∅ , r i G G I W  x i ; ζ i  X = { x } , (44b) where ζ = ( a, b, m, P , v , V ) is the GGIW density parameters. The negati ve cross entropy between f h ( · ) and e f i ( · ) is Z f h ( X ) log e f i ( X ) δ X =  1 − r h  log  1 − r i  + r h log r i + r h  Z N  ξ ; m h , P h  log N  ξ ; m i , ˆ P i  dξ + Z G AM  γ ; a h , b h  log G AM  γ ; a i , b i  dγ + Z I W  χ ; v h , V h  log I W  χ ; v i , V i  dχ  , (45a) Z N  ξ ; m h , P h  log N  ξ ; m i , P i  dξ = 18 − d ξ 2 log(2 π ) − 1 2 log   P i   − 1 2 T r  P h +  m h − m i   m h − m i  T   P i  − 1  , (45b) Z G AM  γ ; a h , b h  log G AM  γ ; a i , b i  dγ = a i log b i − log Γ  a i  +  a i − 1   ψ 0  a h  − log b h  − b i a h b h , (45c) Z I W  χ ; v h , V h  log I W  χ ; v i , V i  dχ = − d χ + 1 2 log     V h 2     + v i 2 d χ X j =1 ψ 0  v h − d χ − j 2  − log Γ d χ  v i − d χ − 1 2  − v h − d χ − 1 2 T r   V h  − 1 V i  + v i − d χ − 1 2 log     V h  − 1 V i    , (45d) where d ξ is the dimension of ξ , d χ is the dimension of χ , Γ( · ) is the gamma function, Γ d χ ( · ) is the multiv ariate gamma function, and ϕ 0 ( · ) is the digamma function.