1
PiVoT: A Variational Solution for Real-time Large-scale Multi-object Detection and Tracking under Heavy Clutter
arXiv:2607.13891v1 [cs.LG] 15 Jul 2026
Runze Gan, Member, IEEE, Qing Li, Member, IEEE, Simon J. Godsill, Fellow, IEEE, Mike E. Davies, Fellow, IEEE, and James R. Hopgood, Senior Member, IEEE
Abstract—Multi-object detection and tracking from noisy point clouds remain challenging in many data-scarce radar applications. Current Bayesian trackers based on Poisson measurement models offer a training-free solution but struggle to achieve accuracy and efficiency under severe clutter, large object populations, and full-resolution Doppler point clouds. We address this with PiVoT, a fast, clutter-resilient multi-object tracker for both positional and Doppler measurements. PiVoT performs end-toend detection and tracking of a large and time-varying number of objects without external clustering or detectors, through joint inference of object states, shapes, existence probabilities, data association, and measurement rates. Its efficiency is driven by several variational inference innovations, such as theoretically justified birth pruning, quadratic-to-linear complexity reductions for exact updates, and a computationally efficient Doppler Poisson model. Experiments show that PiVoT substantially outperforms existing Bayesian trackers in challenging scenes, while also demonstrating exceptional scalability to a thousand objects, robustness to clutter visually inseparable from objects, and realtime operation on full-scale modern automotive radar datasets, where it attains performance comparable to a deep-learning detection benchmark as a training-free joint detector and tracker. Index Terms—Multi-object tracking, radar point clouds, variational inference, object detection, Doppler measurements, Poisson measurement model, data association, Bayesian inference.
I. I NTRODUCTION
D
ETECTING and tracking multiple objects directly from noisy point-cloud-like measurements [1]–[4] is increasingly relevant for modern sensing systems. A prominent example is high-resolution radar that provides many spatial returns with Doppler velocity over long detection ranges (see Fig. 1 for an illustration), even under poor weather conditions [5]. While recent advances in point-cloud deep learning have enabled effective detection and subsequent tracking in datarich automotive settings [1], [3], labelled data scarcity remains This work was supported through the SIGNeTS project (W911NF-20-20225) and the Engineering and Physical Sciences Research Council (EPSRC) under grant EP/X025365/1, in collaboration with Leonardo, as part of the ‘Smart Products Made Smarter’ project. R. Gan is with the Institute for Imaging, Data and Communications (IDCOM), University of Edinburgh, Edinburgh EH9 3FG, U.K., and also with the Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, U.K. (e-mail: [email protected]; [email protected]). Q. Li is with the School of Mathematics, University of Edinburgh, Edinburgh EH9 3FD, U.K., and also with the Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, U.K. (e-mail: [email protected]; [email protected]). S. J. Godsill is with the Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, U.K. (e-mail: [email protected]). M. E. Davies, and J. R. Hopgood are with IDCOM, University of Edinburgh, Edinburgh EH9 3FG, U.K. (e-mail: {mike.davies, james.hopgood}@ed.ac.uk).
a fundamental challenge in many radar-dependent applications [6], [7], and is especially restrictive in surveillance radar, where controlled data collection and annotation require costly and time-consuming field trials. In such a data-constrained setting, multi-object tracking has historically been formulated within a model-based Bayesian framework [9], [10], where object kinematics and measurement generation are modelled probabilistically, and detection and tracking are performed as approximate Bayesian inference without the need for training data. To account for the point-cloud-like measurements from modern high-resolution sensors, the non-homogeneous Poisson process (NHPP) measurement model was introduced in [11], characterising a generation process where each object produces a number of measurements drawn from a Poisson distribution. Many Bayesian multi-object trackers have been actively developed based on this NHPP measurement model, commonly referred to as extended-object trackers [2]. Examples include random-finite-set-based [2], [12], [13] and vectorbased methods [14]–[16], employing inference techniques such as multiple-hypothesis tracking [12], belief propagation [13], [14], and Monte Carlo methods [15], [16]. While effective in ideal scenarios, these trackers often struggle with accuracy and efficiency in challenging real-world settings with heavy clutter and a large number of objects. Such clutter may arise from noisy backgrounds, as commonly observed in automotive radar data (see Fig. 1(b)), and from sea or rain clutter [17]. Large-scale tracking further involves many objects and massive data, for example when tracking drone swarms or groups of insects or animals. A further limitation of existing Bayesian extended-object trackers is that they often neglect or inefficiently exploit valuable Doppler point cloud information, which is routinely provided by many automotive and surveillance radars [5], [8]. To this end, we introduce PiVoT (Poisson Measurementsbased Variational Multi-object Detection and Tracking), an efficient variational inference framework for detecting and tracking a large, time-varying number of objects directly from point-cloud measurements, robust to dense clutter and supporting optional incorporation of Doppler information. PiVoT does not require training data and jointly estimates object kinematic states, shapes, Poisson rates (the mean measurement counts), existence and detectability probabilities within a unified probabilistic model. Its performance on challenging real automotive
2
Fig. 1: PiVoT joint detection and tracking on real radar data, with associated workflow. (a) Camera image from the RadarScenes dataset [8] (not used by PiVoT). (b) Corresponding radar point cloud with Doppler shown by arrows, where coloured points and arrows mark ground-truth moving-object returns. PiVoT estimates are overlaid as ellipses (shapes) and curves (tracks), with coloured estimates for moving objects and light-grey estimates for stationary objects or clutter, reliably recovering the true movers under heavy clutter. (c) Satellite view of the maritime surveillance area. (d) Maritime radar measurements (blue point clusters), overlaid by PiVoT estimates (red ellipses: shapes; curves: tracks). (e) Workflow of PiVoT at each time step. and maritime radar data§ is illustrated in Fig. 1. A. Related work The efficiency of PiVoT arises from a novel variational inference routine tailored to the NHPP measurement model with uncertain object existence, admitting closed form, parallelisable coordinate ascent variational inference (CAVI) updates [19], [20] and a linear-time global optimum, see Fig. 1(e). Although CAVI has been applied to NHPP-based tracking in [21], [22], these methods assume a fixed number of objects, and extending them to a time-varying cardinality introduces substantial technical difficulties (see Section III-A). Previous works [23]–[25] have also applied CAVI and expectationmaximisation to joint detection and tracking; however, they are under the point object model rather than the NHPP, and [23], [24] further require iterative belief propagation at each CAVI iteration. Extending these methods to NHPP is also problematic, [24] loses closed-form updates, while [25] retains the original issue of exponential scaling with the number of objects. In contrast to [23]–[25], PiVoT preserves closedform structure and scalability via a two-stage variational inference that targets distinct marginal posteriors to exploit the model’s tractable structure. It also addresses the computational bottlenecks in existing NHPP-based trackers [12]–[16]. To handle a varying number of objects, PiVoT uses existence indicators and birth models, as in [12]–[15], [23]–[25]. A key innovation to improve efficiency is adaptive birth removal (Fig. § Video demonstrations of PiVoT on real automotive and maritime radar data, together with simulated heavy-clutter and large-scale tracking scenes, are provided in [18].
1(e), Stage 1), which progressively removes ineffective births during CAVI and is formally supported by Theorem IX.1. PiVoT extends naturally to Doppler point clouds through our Doppler-augmented NHPP model, which, unlike existing nonlinear Doppler models [26]–[28], retains a linear Gaussian per-object likelihood and thus preserves the efficient closedform update structure of the positional NHPP model [11]. PiVoT directly processes all point cloud measurements without an external detector or clustering. Instead, its end to end joint inference inherently provides clutter robust, parallelisable clustering to detect an unknown number of objects while incorporating previous tracking results. This distinguishes PiVoT from existing NHPP-based trackers, e.g. [12]–[14] that benefit from external clustering for data association, track initiation, or measurement censoring, and from automotive radar trackers, e.g. [28], [29] that require clustering or pretrained networks as detectors. For Doppler point clouds, PiVoT does not prune near-zero Doppler measurements, as they may arise from tangential motion. Unlike conventional clustering methods [30], [31] that group measurements by similar Doppler values, PiVoT accounts for geometry and kinematics, enabling detection of objects that may generate measurements with quite different Doppler values. See Section VI-C for details. B. Contributions Our key contribution is to present PiVoT, a fast, trainingfree, clutter resilient method for large-scale multi-object detection and tracking, applicable to both positional and Doppler point cloud measurements. Compared to existing NHPP-based trackers, PiVoT delivers substantial performance gains in challenging settings, as shown in Section VIII. It also demonstrates
3
the following capabilities (illustrated in supplementary videos [18]) beyond existing Bayesian trackers: 1) robustness in dense clutter, even when objects are difficult to distinguish by human inspection; 2) detecting and tracking a thousand objects in under a second on a standard laptop without gating; 3) to our knowledge, being the first fully training-free§ detector and tracker to run efficiently on full-scale modern automotive radar datasets, with performance comparable to a deep-learning detection benchmark as shown in Section VIII. At the same time, PiVoT additionally provides estimates of tracks, shapes, data association, detectability and existence probabilities. The central methodological contribution of PiVoT is a novel two-stage variational inference framework (Section III-B) that targets distinct marginal posteriors to exploit the NHPP model’s tractable structure. This design addresses several inherent challenges of applying standard variational inference directly to the full posterior, as detailed in Section III-A. Within this two-stage framework (see Fig. 1(e)), three further methodological innovations enable PiVoT’s efficiency and strong empirical performance: 1) a procedure for early removal of ineffective birth objects that would otherwise converge to having no associated measurements, with theoretical guarantees developed in Theorem IX.1 in Section IX, thereby greatly accelerating Stage 1 inference; 2) an exact global optimiser for object existence and detectability probabilities in Stage 2, which reduces the standard quadratic complexity to linear time, as detailed in Section V; and 3) a Doppler-augmented NHPP model (Section VI-A) that, to our knowledge, is the first NHPP model to incorporate Doppler information while preserving a linear Gaussian per-object likelihood, providing a convenient building block for efficient handling of Doppler point clouds in both PiVoT and future trackers. A preliminary version of PiVoT was presented in [32]. The current paper substantially extends it with broader experiments and technical developments, including deeper CAVI theory and proofs, a numerically stable and efficient existence evaluation, a principled treatment of mean-field limitations, and a Doppler-augmented NHPP model (Doppler PiVoT) validated on automotive radar point clouds from RadarScenes [8].
C. Paper layout The rest of the paper is organised as follows. Section II formulates the problem as an inference task under a general NHPP model. Section III discusses the main inference challenges and introduces the proposed two-stage variational inference framework, with detailed Stages 1 and 2 inference steps presented in Sections IV and V. Section VI presents the Doppler-augmented NHPP model and Doppler PiVoT. Section VII describes the implementation details. Section VIII presents the experimental results, while the theoretical analysis is deferred to Section IX. Section X concludes the paper. Tables IV–VI in Appendix A summarise the main notation used throughout the paper. All appendices cited in the main paper are provided in the supplementary material. § No supervised learning is used in either detection or tracking.
II. P ROBLEM FORMULATION We assume that existing objects may cease to exist (die) and new objects may emerge (birth) at each time step, as in [12]–[15], [23]. At time step n, consider that there are Kn = Kn−1 + Knb potentially existing objects. These Kn objects include Kn−1 legacy objects from the previous time step n − 1, and Knb newly born objects at the current time step n. We assign each of the Kn objects an index k ∈ {1, ..., Kn }. Legacy objects from time n − 1 are indexed as k = 1, ..., Kn−1 , while newly born objects at time n are indexed as k = Kn−1 + 1, ..., Kn . For each object k = 1, ..., Kn , we define the following variables. The kinematic state Xn,k is a vector containing position, velocity, and any other motion parameters. The scalar Poisson rate Λn,k > 0 equals the expected number of measurements generated by object k at time n. The precision matrix Pn,k has an inverse that describes the elliptical shape of object k. The binary indicators Dn,k ∈ {0, 1} and En,k ∈ {0, 1} represent detectability and existence respectively. More details about these features will be introduced in Section II-A. The clutter Poisson rate (the expected number of clutter measurements) is denoted by Λn,0 . Finally, we denote Xn , Λn , Pn , Dn and En as the respective collections of these variables for all objects k = 1, ..., Kn , for example, Λn = [Λn,1 , ..., Λn,Kn ]. The measurements at time step n are denoted by Yn = [Yn,1 , Yn,2 , ..., Yn,Mn ], where Mn is the total number of measurements. We define the measurement-oriented association θn = [θn,1 , θn,2 , ..., θn,Mn ], where each θn,j ∈ {0, 1, ..., Kn } for j = 1, 2, ..., Mn indicates the origin of the measurement Yn,j . Specifically, θn,j = k indicates that Yn,j was generated by object k if k = 1, ..., Kn , or by clutter if k = 0. A. Model for measurements and association We assume that existing objects produce measurements following an NHPP measurement model [11], [22] with a detection probability pdn,k ∈ (0, 1] [12]. Specifically, an object k = 1, . . . , Kn at time n is detectable only if it exists, and given existence En,k = 1 it is detectable, Dn,k = 1, with probability pdn,k . An object can generate measurements only when Dn,k = 1, so the detectability indicator is useful to model occasional occlusion. Mathematically, YKn p(Dn |En ) = p(Dn,k |En,k ), (1) k=1
where p(Dn,k |En,k = 1) is pdn,k if Dn,k = 1 and 1 − pdn,k otherwise, while p(Dn,k |En,k = 0) = δ[Dn,k = 0], where δ[·] is Kronecker delta function. Finally, we assume pdn,k = 1 for all newly born objects k = Kn−1 +1, ..., Kn at current time n. Conditional on Dn,k = 1, object k generates measurements according to an NHPP. Specifically, the number of measurements is Poisson distributed with rate Λn,k , and each measurement is independent and identically distributed with density ℓ(· | Xn,k , Pn,k ). When Dn,k = 0, object k produces no measurements. Consequently, conditional on Dn , the full set of measurements Yn forms an NHPP obtained by superposing conditionally independent NHPPs from clutter and the Kn objects with rates Λn,0 and Dn,k Λn,k respectively. The likelihood is given by [11]:
4
The joint prior for θn , Mn in (4) can be verified to be [22] XKn Dn,k Λn,k ) p(θn , Mn |Λn , Dn ) = M1n ! exp(−Λn,0 − k=1
× h(Λn , Dn , θn ) PMn
h(Λn , Dn , θn ) = Λn,0j=1 =
Mn Y j=1
Fig. 2: Simulated NHPP measurements from three objects over three time steps. Black dots denote clutter, and coloured dots are object-generated measurements with Poisson rate Λn,k and shape Pn,k (shown as ellipses), generated only when Dn,k = 1 (which implies En,k = 1). In particular, object k = 3 no longer exists from time n onwards, while object k = 2 still exists at time n but is not detectable (e.g. due to occlusion).
δ[θn,j =0]
Λn,0 δ[θn,j = 0] +
Kn X
Kn Y
(6) P Mn
(Λn,k Dn,k )
j=1 δ[θn,j =k]
(7)
k=1
Λn,k Dn,k δ[θn,j = k] ,
(8)
k=1
where (7) and (8) provide alternative but equivalent expressions, each useful in different contexts encountered later. B. Model for object dynamics and other features
1) Model for legacy objects: At each time n, we assume that transitions are independent across legacy objects k = 1, 2, ..., Kn−1 and also independent among each object’s PKn features Xn,k , Λn,k , Pn,k , En,k . Specifically, each legacy exp(−Λn,0 − k=1 Dn,k Λn,k ) object survives from time n−1 to n with a survival probability p(Yn , Mn |Dn , Xn , Λn , Pn ) = Mn ! psn,k ∈ (0, 1], i.e. p(En,k |En−1,k = 1) = psn,k if En,k = 1 M K n n Y X s Λn,k Dn,k ℓ(Yn,j |Xn,k , Pn,k ) , (2) and 1 − pn,k if En,k = 0. A legacy object that does not × Λn,0 ℓ0 (Yn,j )+ exist (En−1,k = 0) remains nonexistent in future time steps j=1 k=1 where ℓ0 (Yn,j ) and ℓ(Yn,j |Xn,k , Pn,k ) denote the clutter and (En,k = 0). We assume linear Gaussian state transition with transition matrix Fn and covariance Qn for each legacy object: single-object measurement likelihoods, respectively. Depending on the sensor, each measurement Yn,j may conp(Xn,k |Xn−1,k ) = N (Xn,k ; Fn Xn−1,k , Qn ). (9) tain only position or position together with Doppler velocity. We first consider the classical positional-only model, where Moreover, we adopt the commonly assumed heuristic transiYn,j reduces to a positional measurement yn,j ∈ RdY with tion (e.g. in [2], [12], [13]) for rate Λn,k and shape Pn,k , which dY = 2 for planar and dY = 3 for spatial tracking. For this preserves the Gamma and Wishart distribution forms of Λn,k positional-only model, i.e. Yn,j = yn,j , the clutter likelihood and Pn,k in the prediction step while slightly increasing the −1 ℓ0 in (2) may be uniform over the surveillance area or uncertainty of Λn,k and Pn,k , as detailed in Section VII. 2) Model for newly born objects: At each time n, we follow a predefined spatial map, while the object measurement assume a birth model with a maximum of Knb new births. likelihood ℓ is assumed to follow a linear Gaussian form: −1 ℓ(Yn,j |Xn,k , Pn,k ) = N (yn,j ; HXn,k , Pn,k ) (3) For each newly born object k = Kn−1 + 1, ..., Kn , its birth prior is assumed to be independent across features Xn,k , Λn,k , where H is the observation matrix mapping the object state P , E , as well as independent of all other newly born n,k n,k −1 to the positional measurement space, and Pn,k is the co- and legacy objects’ features. Specifically, the birth priors for variance matrix (with Pn,k being its precision) that captures existence E , state X , rate Λ , and shape (precision n,k n,k n,k the object’s elliptical shape. Example measurements from this matrix) P n,k follow Bernoulli, Gaussian, Gamma, Wishart positional-only NHPP model, together with the corresponding distributions, respectively: Λn,k , Pn,k , Dn,k and En,k are illustrated in Fig. 2. The Doppler NHPP model will be presented in Section VI. p(En,k ) = Ber(pbn,k ), p(Xn,k ) = N (µbn,k , Σbn,k ), (10) b The PiVoT inference framework in Sections III-V applies p(Λn,k ) = G(ηn,k , ρbn,k ), p(Pn,k ) = W(Φbn,k , ϕbn,k ), to both settings, whereas the final update expressions in Section IV are specific to the positional model introduced here; where the Bernoulli parameter pbn,k ∈ (0, 1) denotes the the corresponding Doppler case is given in Section VI. probability of birth, i.e. En,k = 1. µbn,k , Σbn,k denote the mean b By introducing the data association θn , the joint likelihood and covariance, ηn,k , ρbn,k are the shape and scale parameters, of Yn , Mn , θn , Dn can be expressed as follows, which can Φb is the scale matrix, and ϕb is the degrees of freedom. n,k n,k be verified using (2) following a similar procedure in [22]: For reliable detection of all potentially born objects in a cluttered environment, the maximum birth count Knb should p(Yn , Mn , θn , Dn |Xn , Λn , Pn , En ) (4) be sufficiently large, often far exceeding the actual number of = p(Yn |θn , Pn , Xn )p(θn , Mn |Λn , Dn )p(Dn |En ), births. However, PiVoT can adaptively reduce Knb , removing The association-conditioned likelihood in (4) is given by ineffective birth objects during inference to improve efficiency YMn p(Yn |θn , Pn , Xn ) = p(Yn,j |θn,j , Pn , Xn ) (5) and refine the birth model, see Section IV-C and IV-D. For j=1 efficient detection/clustering over large areas, PiVoT typically ( ℓ(Yn,j |Xn,θn,j , Pn,θn,j ), θn,j ̸= 0 uses identical non-informative birth priors p(Xn,k ) (see Secp(Yn,j |θn,j , Pn , Xn ) = tion IV-D), unlike the usual multi-Bernoulli birth setting. ℓ0 (Yn,j ), θn,j = 0
5
C. Inference goal: filtering posterior and predictive prior The goal of inference at time step n is to obtain an accurate and efficient approximation of the following filtering posterior: p̂n (Xn , Λn , Pn , θn , En , Dn |Yn ) ∝ p̂n (Xn , Λn , Pn , En ) × p(Yn , Mn , θn , Dn |Xn , Λn , Pn , En ), (11) where the likelihood term (second line) follows the exact NHPP joint likelihood given in (4)-(6) and (1). The prior term p̂n (Xn , Λn , Pn , En ) consists of the exact birth prior (Section II-B2) and the predictive prior for legacy objects, constructed using the transition model (Section II-B1) and the approximated posterior from the previous time step. In PiVoT, this predictive prior retains the same independent factorisation and distribution forms as the birth prior, leading to: Kn Y p̂n (Xn , Λn , Pn , En ) = p̂n (Xn,k )p̂n (Λn,k )p̂n (Pn,k)p̂n (En,k ) k=1
p̂n (Xn,k ) = N (µ′n,k , Σ′n,k ),
′ p̂n (Λn,k ) = G(ηn,k , ρ′n,k ),
p̂n (Pn,k ) = W(Φ′n,k , ϕ′n,k ),
p̂n (En,k ) = Ber(pe′ n,k ), (12) For newly born objects k = Kn−1 +1, ..., Kn , the prior parameters in (12) are exactly identical to those in the birth prior (10), b b b b i.e. [µ′n,k , Σ′n,k , ..., ϕ′n,k , pe′ n,k ] = [µn,k , Σn,k , ..., ϕn,k , pn,k ]. For legacy objects k = 1, ..., Kn−1 , the prior is the predictive prior, with parameters computed from the prediction step, as detailed in (46) in Section VII. In the rest of this paper, we omit the subscript n in the notation for clarity, unless necessary to avoid confusion. III. P I VOT: I NFERENCE CHALLENGES AND SOLUTION This section first outlines the challenges of approximating the filtering posterior in (11). We then present our solution, a two-stage variational inference framework that serves as the backbone of PiVoT and summarises the method in a nutshell. A. Inference challenges and motivation The key inference routine of PiVoT is CAVI, which demonstrated superior efficiency and reliability for tracking a fixed number of objects under the NHPP model [21], [22]. However, extending it to our models in Section II for a varying number of objects is challenging. A standard mean-field assumption: q(θ, X, Λ, P, E, D) = q(θ)q(X)q(Λ)q(P )q(E, D), (13) fails when approximating p̂ in (11) by minimising the Kullback-Leibler divergence (KLD) KL(q∥p̂). The reason is as follows. Since p̂ is zero whenever a measurement is associated to a non-existent object, q cannot place probability mass there. This forces variational inference (VI) into one of two failures: 1) forcing q(Ek = 0) = 0, meaning that object k must exist with no uncertainty, or 2) preventing object k from associating with measurements, making updates impossible. Either outcome renders the naive mean-field assumption in (13) fundamentally flawed. Detailed derivations are given in Appendix B. Possible remedies include modifying the model to allow non-existent objects to generate measurements [23], but inference performance is highly sensitive to the extra Poisson rates introduced for non-existent objects. Another approach replaces q(θ) in (13) with q(θ|D), which preserves detectability dependence but introduces computationally prohibitive evaluations
of q(θ|D) over 2Kn configurations as in [25]. While general inference methods like Monte Carlo and gradient-based VI could be used, they compromise efficiency. Instead, PiVoT introduces alternative objectives and suitable approximations to retain efficient closed-form CAVI updates, as shown below. B. PiVoT: Two-stage variational inference framework PiVoT follows a mean-field approximation similar to (13), but employs two separate inference stages with distinct, well motivated objectives to address inference challenges. Specifically, it approximates p̂ in (11) as Kn Y q(X, Λ, P, θ, E, D) = q1 (X)q1 (Λ)q1 (P )q1 (θ) q2 (Ek , Dk ), (14) k=1 where q1 and q2 are evaluated in the first and second stages, respectively, as introduced below and summarised in Fig. 1(e). 1) Detection and tracking stage: Stage 1 aims to use VI to approximate the marginal posterior p̂(X, Λ, P, θ|Y ) from (11) with q1 (X)q1 (Λ)q1 (P )q1 (θ):X q1 (X)q1 (Λ)q1 (P )q1 (θ) ≈ p̂(X, Λ, P, θ, E, D|Y ) (15) E,D X ∝ p̂(X, Λ, P )p(Y |θ, P, X) p(θ, M |Λ, D)p̂(D), (16) D
Here, (16) follows from (4) and (11). The factors p̂(X, Λ, P ), p(Y |θ, P, X) and p(θ, M |Λ, D) are specified in (12), (5) and (6), respectively. The shorthand p̂(D) is defined as X YKn p̂(D) := p(D|E)p̂(E) = p̂(Dk ), k=1 (17) XE d p̂(Dk ) = p(Dk |Ek )p̂(Ek ) = Ber(pe′ k pk ), Ek
where the simplified factorised form follows from (1) and (12). Since E, D are marginalised out, the resulting posterior is not necessarily zero for all possible associations θ, avoiding the mean-field issue discussed in Section III-A. Parallelisable CAVI then iteratively updates each variational distribution q1 . This stage forms the core of PiVoT, jointly detecting new objects, tracking legacy objects, and estimating rates and shapes. Notably, the design in Section IV enables automatic clutter-robust clustering for birth detections while improving efficiency through early identification and removal of ineffective births, as supported by Theorem IX.1 in Section IV-C. 2) Existence evaluation stage: In Stage 2, the main objective is to independently evaluate q2 (Ek , Dk ) in (14) for each object k = 1, 2, ..., Kn , by targeting the following object-wise marginal approximation with VI: q2 (Ek ,Dk )q2 (θ|Ek , Dk )q1 (X)q1 (Λ)q1 (P ) X ≈ p̂(X, Λ, P, θ, E, D|Y ) (18) Ek− ,Dk−
∝ p̂(X, Λ, P )p̂(Ek )p(Dk |Ek )p(Y |θ, P, X) X × p(θ, M |Λ, D)p̂(Dk− ), (19) Dk−
where q1 (X)q1 (Λ)q1 (P ) is obtained from the first stage and remains fixed in the second stage. The subscript k− denotes all indexed variables except the k-th one, e.g. Ek− includes all components of E except Ek . (19) is derived similarly to (16), with p̂(Dk− ) obtained directly from (17). The rationale behind this stage is that q1 (X)q1 (Λ)q1 (P ) from the first stage should already approximate p̂(X, Λ, P |Y )
6
well, so keeping it fixed allows efficient evaluation of q2 (Ek , Dk )q2 (θ|Ek , Dk ) to approximate p̂(Ek , Dk , θ|Y ), though only q2 (Ek , Dk ) is required in (14). Here q2 (θ|Ek , Dk ) depends on Ek , Dk to avoid the meanfield issue, and to avoid the 2Kn configurations when conditioning on all E, both discussed in Section III-A. Notably, Section V presents a method that directly computes the globally optimal q2 (Ek , Dk ) for all k = 1, ..., Kn without evaluating q2 (θ|Ek , Dk ). This makes the second stage surprisingly fast, with the total cost of evaluating all q2 (Ek , Dk ) comparable to just two iterations of q1 (θ) updates in the first stage; thus Stage 2 takes much less time than Stage 1. For each object k, Stage 2 can be extended with an optional refinement step to improve the accuracy of inferred object features, as detailed in Appendix D-2. C. Approximation for efficient two-stage inference The two-stage VI requires the unnormalised marginal posteriors in (16) and (19), wherein the marginal associations in each are respectively given as follows using (6), (7): Kn PM X e−Λ0 δ[θj =0] Y Λ0 j=1 Ξk (Λk , θ), p(θ, M |Λ, D)p̂(D) = M! D k=1 PM X e−Λ0 −Dk Λk δ[θj =0] (20) p(θ, M |Λ, D)p̂(Dk−) = Λ0 j=1 M! Dk− PM Y × (Λk Dk ) j=1 δ[θj =k] Ξk (Λk , θ), k−
h i P Mn Ξk (Λk , θ) := Ep̂(Dk ) e−Dk Λk (Λk Dk ) j=1 δ[θj =k] .
(21)
Although Ξk (Λk , θ) can be expressed analytically, their product complicates efficient optimisation in both stages. To enable efficient two-stage inference in PiVoT, we introduce one approximation: e′ d
PMn
d j=1 δ[θj =k] . Ξk (Λk , θ) ≈ e−pk pk Λk (pe′ (22) k pk Λk ) This approximation replaces Dk in Ξk with its mean d Ep̂(Dk ) [Dk ] = pe′ k pk using (17), corresponding to a first-order Taylor approximation of the expectation Ξk [33]. An alternative and perhaps more principled interpretation of (22) is that it replaces a marginalised point process within the superposed measurement process by an NHPP with intensity d pe′ k pk Λk ℓ(Yj |Xk , Pk ). This approximation is detailed and justified in Appendix C, where it is shown to be the unique NHPP minimising the KLD from the original marginalised process. Furthermore, Appendix D analyses the impact of the approximation in (22) on inference accuracy. A key takeaway is that, for an object sufficiently far from others, the approximation has minimal impact on accuracy, and residual error can be removed through the additional Stage 2 refinement steps.
IV. S TAGE 1: D ETECTION AND TRACKING This section first presents the Stage 1 evidence lower bound (ELBO) [19], then the corresponding CAVI updates and an early identification and removal procedure for ineffective birth objects, both under the Stage 1 block in Fig. 1(e). This removal step, for which we develop a formal guarantee, enables efficient discovery of birth objects over large surveillance areas.
A. Stage 1 ELBO Recall from Section III-B1 that Stage 1 aims to approximate the marginal posterior in (15) with variational inference. This amounts to optimising q1 to maximise the following ELBO: F1 = Eq1 (X)q1 (Λ)q1 (P )q1 (θ) log
p̄(X, Λ, P, θ, Y ) , q1 (X)q1 (Λ)q1 (P )q1 (θ)
where the joint density p̄(X, Λ, P, θ, Y ) is defined as p̄(X, Λ, P, θ, Y ) := p̂(X, Λ, P )p(Y |θ, P, X) PK
× M1 ! e−Λ0 −
e′ d k=1 pk pk Λk
d e′ d h(Λ, D = [pe′ 1 p1 , ..., pK pK ], θ), (23)
where h(·) isP defined in (7). The last line in (23) approximates the marginal D p(θ, M |Λ, D)p̂(D) in (20) using (22), so that the resulting CAVI updates admit closed forms. This ELBO F1 enables efficient CAVI to achieve the Stage 1’s inference goal in (15). Specifically, CAVI minimises KL(q1 (X)q1 (Λ)q1 (P )q1 (θ)∥p̄(X, Λ, P, θ|Y )) by iteratively updating each of q1 (X), q1 (Λ), q1 (P ), q1 (θ) to its optimal value while keeping others fixed. Each update ensures a non-increasing KLD, and subsequently, convergence to a stationary point is guaranteed. B. Coordinate ascent updates This section presents the updates for each q1 (X), q1 (Λ), q1 (P ), q1 (θ), using standard CAVI update formula [19], [20]: q1 (zi ) ∝ exp Eq1 (zi− ) log p̄(X, Λ, P, θ, Y ) , (24) where zi represents each variable in {X, Λ, P, θ}, zi− denotes all remaining variables, and p̄ is given in (23). The final update expressions given below are for the positional-only likelihood in (3); the corresponding Doppler-augmented updates are presented in Section VI. The specific structure of p̄ in (23) implies that the updated q1 obtained via the CAVI formula in (24) factorises as Kn M Y Y q1 (X)q1 (Λ)q1 (P )q1 (θ) = q1 (Xk )q1 (Λk )q1 (Pk ) q1 (θj ), k=1
j=1
q1 (Xk ) = N (µk , Σk), q1 (Λk ) = G(ηk , ρk), q1 (Pk ) = W(Φk , ϕk) (25) where (µk , Σk ), (ηk , ρk ), (Φk , ϕk ) are the corresponding variational parameter pairs for the Gaussian, Gamma, and Wishart distributions. Each parameter pair is updated in turn while the others are held fixed and used in the required expectations. The following updates will apply to all legacy and newlyborn objects, differing only in the respective model parameters. State update: Setting zi = X in (24) and substituting (25), (23), (12), (5), the CAVI update for q1 (X) reduces to updating each q1 (Xk ) for k = 1, ..., Kn independently according to YM q1 (Xk ) ∝ p̂(Xk ) exp(q1 (θj = k)Eq1 (Pk ) log ℓ(Yj |Xk , Pk )) j=1 (26) Using (25) and the positional-only likelihood ℓ in (3), each µk , Σk in (25) can be obtained independently from a Kalman filter with prior p̂(Xk ) in (12) and a linear Gaussian likelihood:
7
q1 (Xk ) ∝ p̂(Xk )N y k ; HXk , Rk , (27) PM −1 Φk j=1 yj q1 (θj = k) . , y k = PM Rk = PM ϕk j=1 q1 (θj = k) j=1 q1 (θj = k) Rate update: Similarly, setting zi = Λ in (24) and substituting (25), (23), (12), the update of q1 (Λ) reduces to independently evaluating each q1 (Λk ) for k = 1, ..., Kn as
their current variational distribution (i.e. without restarting the optimisation, as illustrated in Fig. 1(e)). This procedure is formally justified by Theorem IX.1 in Section IX, informally stated below. In short, any removed birth is provably redundant at CAVI convergence. Theorem IX.1 (informal). Under mild modelling conditions, one can compute a sufficient threshold L such that, if ςk := P M PM j=1 q1 (θj = k) < L at any CAVI iteration for any newly born j=1 q1 (θj =k) e′ d q1 (Λk ) = G(ηk , ρk) ∝ p̂(Λk ) exp(−pk pk Λk )Λk , object k, then subsequent CAVI updates drive ςk monotonically M to 0. Equivalently, at convergence, that birth object assigns X d ′ ηk = ηk′ + q1 (θj = k), ρk = ρ′k /(pe′ p ρ + 1). (28) zero probability of associating with any measurement. k k k j=1 The formal theorem and full analysis are deferred to SecShape update: For the CAVI update of q1 (P ), using (24) with tion IX. To our knowledge, existing CAVI analyses do not (25), (23), (12), (5), it becomes evaluating each q1 (Pk ) via provide a computable sufficient condition certifying that a YM component will be driven to a zero-association fixed point. We q1 (Pk ) ∝ p̂(Pk ) exp(q1 (θj = k)Eq1 (Xk ) log ℓ(Yj |Xk , Pk )). derive such a condition by bounding how the total association j=1 mass propagates through subsequent coordinate updates, and With the prior p̂(Pk ) = W(Φ′k , ϕ′k ) in (12) and the positionalby using the Lambert-W characterisation in Corollary H.2.1 only likelihood ℓ in (3), this gives q1 (Pk ) = W(Φk , ϕk ) with in Appendix H to obtain an explicit sufficient threshold. parameters updated independently for each k = 1, ..., Kn as: As Theorem IX.1 gives only a sufficient condition, in M h i X practice we set L heuristically, for example L = 0.5 for speed −1 ′ −1 ⊤ ⊤ Φk = Φ k + q1 (θj = k) (yj −Hµk )(yj −Hµk ) +HΣk H or L = 0.05 for more conservative pruning, noting that an j=1 M overly large L may discard valid births. X q1 (θj = k). (29) ϕk = ϕ′k + The modelling conditions required by Theorem IX.1 are j=1 given in Assumption 1 in Section IX. Notably, they are satisAssociation update: Finally, for the update of q1 (θ), using fied by the efficient Model Scenario in the following Section (24) with (25), (23), (12), (5), the association for each data Yj IV-D, together with an extra uniform clutter assumption. (i.e. q1 (θj )) is independently updated for j = 1, ..., M as Kn D. Efficient clutter-robust birth detection: parameter suggesX d k pe′ q1 (θj ) ∝ Λ0 ℓ0 (Yj )δ[θj = 0] + tions and demo k pk Sj δ[θj = k], (30) k=1 Denote by Xkp and Xkv the positional and remaining compok Sj := exp(Eq1 (Λk ) log Λk +Eq1 (Xk )q1 (Pk ) log ℓ(Yj |Xk , Pk )). nents of the kinematics state X . We recommend the following k model scenario for efficient PiVoT implementation: Substituting the positional-only likelihood ℓ in (3), we have Model Scenario: For each newly-born object k, assume the Sjk = exp(Tjk ), (31) birth prior in (10) is factorisable as p(Xk ) = p(X p )p(X v ), k k h i p Tjk := − 12 (yj −Hµk )⊤ϕk Φk (yj −Hµk )+Tr(ϕk Φk HΣk H ⊤) and that p(Xk ) is a flat Gaussian approximating a uniform positional birth prior. + ψ(ηk )+log ρk + 12 ψdY ( ϕ2k )+log |Φk |−dY log(π) (32) In this model scenario, the uninformative positional birth prior greatly simplifies birth state update: (27) reduces to where dY is the dimension of yj , | · | denotes the determinant, q1 (Xk ) = N (Xkp ; y k , Rk )p(Xkv ); see Remark 1 under Asand Tr(·) denotes the trace. ψ(·) and ψdY (·) are the digamma sumption 1 in Section IX. Hence, the updated birth state is and multivariate digamma functions, respectively. obtained directly without using a Kalman filter in any CAVI iteration, which reduces the per birth computational cost. C. Early identification of ineffective births For efficient clutter-robust birth detection, we suggest: Exhaustive birth detection under heavy clutter requires (i) The above Model Scenario holds, i.e. a factorised prior initially setting a birth count Knb far above the true num- with a flat Gaussian p(Xkp ) is used for all newly-born objects. ber, increasing feature updates and association dimensionality. (ii) Constrain each newly born object’s initial positional However, many ineffective births that ultimately fail to asso- variational distribution q1 (Xkp ) to a distinct small region ciate with measurements upon CAVI convergence are naturally (crucial for updated q1 (Xk ) to localise the true object within irrelevant to new object detection. Identifying and removing it under dense clutter [22]). This can be achieved by setting a them early, by eliminating them from the prior and reducing small covariance Σpk with a distinct mean µpk , where µpk and Σpk the association dimension, greatly improves CAVI efficiency. denote the corresponding positional subvector and submatrix PiVoT PMimplements this by monitoring the total association of the variational parameters µk and Σk in (25). mass, j=1 q1 (θj = k), in each CAVI iteration. If this count (iii) Tile the entire surveillance area with these small regions falls below a predefined threshold L for newly born objects for exhaustive detection, which implicitly determines an often k = Kn−1 + 1, ..., Kn , those births are deemed ineffective very large Knb . However, ineffective birth removal in PiVoT and removed from the model, and subsequent CAVI iterations Stage 1 inference (Section IV-C) then adaptively reduces Knb , update only the remaining objects’ X, Λ, P , continuing from improving efficiency and refining the model.
8
Fig. 3: Clutter-robust detection using PiVoT Stage 1. The data is from Dataset 3 (first time step) in Section VIII-A, showing a zoomed-in region of the surveillance area. Blue and red dots represent clutter and true object measurements. Black ellipses denote the positions in q1 (Xk ) and shapes q1 (Pk ) in that iteration. (a) Initial CAVI iteration with 1312 births across the full surveillance area. (b) 3rd iteration with 392 remaining births after ineffective birth removal (Section IV-C). (c) 24th iteration (converged) with 47 remaining births. Red-highlighted ellipses mark objects with existence probability above 0.1, computed in Stage 2 (Section V) and displayed next to them. Notably, these red ellipses precisely correspond to true objects. An example of such an initial q1 (Xk ) is shown in Fig. 3(a), where each initial q1 (Xk ) covers a distinct small region, and collectively spans the full surveillance area (extending beyond the figure). This results in an initial Knb of 1312. Fig. 3(b) and 3(c) further illustrate the reliable clustering performance of PiVoT Stage 1 under dense clutter, using the same parameter settings as in Section VIII-A. The proposed early ineffective birth removal progressively reduces Knb to 952, 631, 392 in the first three iterations, ultimately converging at 47 by iteration 24, greatly improving efficiency. Each iteration refines cluster positions while gradually eliminating those containing only sparse measurements. Notably, the relatively high existence probabilities (computed in PiVoT Stage 2) accurately distinguish true births from clutter-induced clusters in Fig. 3(c). V. S TAGE 2: E XISTENCE EVALUATION This section formulates the Stage 2 ELBO and its optimal q2 (Ek , Dk ) over existence and detectability (see Fig. 1(e), Stage 2). While direct evaluation scales as O(Kn2 M ), we derive an exact closed-form expression that reduces the cost to O(Kn M ). We also provide a numerically stable implementation that preserves the O(Kn M ) complexity, making Stage 2 comparable in cost to a few updates of q1 (θ) in (30). A. Stage 2 ELBO and the optimal variational solution 1) Stage 2 ELBO: Recall from Section III-B2 that Stage 2 inference targets (18). Specifically, holding Stage 1 solutions q1 (X)q1 (Λ)q1 (P ) fixed, Stage 2 seeks q2 (θ, Ek , Dk ) by maximising, for each object k = 1, ..., Kn , the ELBO
p̄(θ, M |Λ, Dk ) := M1 ! e−Λ0 −Dk Λk −
PKn
e′ d s=1,s̸=k ps ps Λs
(34)
d e′ d e′ d e′ d ×h(Λ, D = [pe′ 1 p1 , ..., pk−1 pk−1 , Dk , pk+1 pk+1 , ..., pK pK ], θ),
where h is given in (8). The only approximation in F2k , relative to the ELBO under the exact Stage 2 unnormalised posterior in (19), isP the use of p̄(θ, M |Λ, Dk ) in (34) to approximate the marginal Dk− p(θ, M | Λ, D)p̂(Dk− ) in (19) and (20). This approximation follows directly from (22), in the same manner as the Stage 1 ELBO formulation in Section IV-A. 2) Optimal variational distribution: Maximising the Stage 2 ELBO F2k with Stage 1 solution q1 (X)q1 (Λ)q1 (P ) fixed admits a closed-form global optimum q2∗ (Ek , Dk )q2∗ (θ|Ek , Dk ) for each object k = 1, 2, ..., Kn . Using the optimal update formula for conditional variational distributions [20], [34], this optimum is given by log q2∗ (θ|Ek , Dk ) = Eq1 (X,Λ,P ) log p̄(X, Λ, P, θ, Ek , Dk , Y ) + c1 ∗ log q2 (Ek , Dk ) = c2 − Eq2∗ (θ|Ek ,Dk ) log q2∗ (θ|Ek , Dk )
(35) (36)
+ Eq1 (X,Λ,P )q2∗ (θ|Ek ,Dk ) log p̄(X, Λ, P, θ, Ek , Dk , Y ), where c1 , c2 are constants. Recall that only q2∗ (Ek , Dk ) is needed as the filtering posterior in (14). However, from (36), direct computation requires first evaluating q2∗ (θ|Ek , Dk ), which alone entails M evaluations of the categorical distribution q2∗ (θj |Ek , Dk ) of size Kn + 1, leading to a total of O(Kn2 M ) exponential/logarithm operations for all k = 1, ..., Kn . This makes direct computation infeasible for large Kn and M . Next, we introduce a surprisingly efficient method to compute all q2∗ (Ek , Dk ) directly.
B. Efficient evaluation of the optimal q2∗ (Ek , Dk ) p̄(X, Λ, P, θ, Ek , Dk , Y ) q2 (θ, Ek , Dk)q1(X)q1(Λ)q1(P ) Here we derive an exact expression for q2∗ (Ek , Dk ) with where the structured variational distribution q2 (θ, Ek , Dk ) = overall complexity O(Kn M ) for all k = 1, . . . , Kn . The on θ by substituting q2 (Ek , Dk )q2 (θ|Ek , Dk ) is used to fully preserve conditional key idea is to collapse the dependence ∗ the normalised form of log q (θ|E , D ) into (36), and then k k 2 dependence, and the joint density p̄ is defined as rearranging the resulting expression until the summation over p̄(X, Λ, P, θ, Ek , Dk ,Y ) = p̂(X, Λ, P )p̂(Ek )p(Dk |Ek ) θ becomes tractable. The final form reveals that the expensive sum ∗ × p(Y |θ, P, X)p̄(θ, M |Λ, Dk ), (33) part (Sj in (39)) for computing q2 (Ek , Dk ) is shared across
F2k = Eq2 (θ,Ek ,Dk)q1(X)q1(Λ)q1(P) log
9
different k and therefore needs to be computed only once. This shared structure is not evident from the direct update in (36), yet is crucial for efficient evaluation. We defer the detailed derivation to Appendix E and present the final form below: q2∗ (Ek , Dk ) ∝ p̂(Ek )p(Dk |Ek ) exp(g(Dk )), (37) M X d k g(Dk ) = −ηk ρk Dk + log Sjsum + (Dk − pe′ (38) k pk )Sj j=1
Sjsum = Λ0 ℓ0 (Yj ) +
Kn X
d k pe′ k pk Sj ,
(39)
k=1 where Sjk is defined in (30), and Sjsum is computed once for all k. If q1 (θ) is the last updated distribution in Stage 1, all Sjk have already been computed in (30) along with Sjsum as the normalisation constant. The other quantities p̂(Ek ) and pe′ k are defined in (12), p(Dk |Ek ) and pdk in (1), and ηk , ρk in (25). The globally optimal q2∗ (Ek , Dk ) in (37) and (38) is now
in a highly simplified form, requiring only O(Kn M ) exponential/logarithm operations (the main computational cost) for updating all Kn objects, compared to at least O(Kn2 M ) via direct computation in (36). This is much more efficient and matches the complexity of updating q1 (θ) in (30). Finally, the object k’s existence probability is X q2∗ (Ek ) = Ber(pek ) ∝ p̂(Ek ) p(Dk |Ek ) exp(g(Dk )) (40) Dk
pek = le /(le + ln ), ln = (1 − pe′ k ) exp(g(Dk = 0)), e′ d le = pk pk exp(g(Dk = 1)) + (1 − pdk ) exp(g(Dk = 0)) , C. Numerically stable computation of g(Dk ) The numerically stable evaluation of g(Dk ) using (38) requires care, especially for g(Dk = 0), which involves leaveone-out log-sums. A naive stable evaluation would recompute these sums using log-sum-exp for each k, reintroducing quadratic complexity. We avoid this by reusing stable global sums where possible, with a separate treatment only when the removed term is the largest one in the sum. Appendix F gives the detailed numerically stable computation of g(Dk ), specified in (65)–(68), which preserves the O(Kn M ) complexity. VI. P I VOT WITH D OPPLER V ELOCITY M EASUREMENTS PiVoT can be readily extended to incorporate Doppler velocity measurements, which supply crucial radial motion information commonly exploited in radar systems to greatly improve tracking accuracy. This section first introduces a new Doppler-incorporated Poisson measurement model that retains the efficient linear Gaussian single-object likelihood. We then derive the corresponding PiVoT updates, where most of the existing updates remain unchanged. A. Doppler-augmented Poisson measurement model Doppler velocity measures the relative radial velocity of a detection point, namely the component of its actual velocity along the line of sight to the sensor. Fig. 4 illustrates this geometric effect. To make the model concrete, we consider a 2D Cartesian sensing setting. Each object state Xk contains planar velocity
components vkx , vky for the x and y axis, and we let G extract them so that GXk = [vkx , vky ]⊤ . For any positional measurement yj ∈ R2 , let αj be its bearing angle relative to the sensor and define the line-of-sight unit vector ej := [cos αj , sin αj ]. The true radial velocity of object k at this bearing is therefore vkj = ej GXk . A visual illustration of vkx , vky , vkj , αj for a detection point yj is provided in Fig. 4(b). Since Doppler is naturally attached to each positional detection, we model each measurement Yj = (yj , vj ) in a marked NHPP, where each positional measurement yj is equipped with a mark vj ∈ R. This retains the NHPP likelihood in (2), with ℓ0 and ℓ defined below. Assuming the observed Doppler velocity vj is a Gaussian measurement of the true radial velocity vkj , the single-object likelihood still admits a convenient linear Gaussian form: ℓ(Yj | Xk , Pk ) = N (vj ; ej GXk , σv2 ) N (yj ; HXk , Pk−1 ), ej := [cos αj , sin αj ],
(41)
where σv2 is the Doppler noise variance. As defined above, αj is the bearing angle of yj (see Fig. 4(b)), and G and H are the velocity and positional observation matrices, respectively. For instance, when Xk = [pxk , vkx , pyk , vky ]⊤ represents the position and velocity in the x and y axes, we have 0 1 0 0 1 0 0 0 G= , H= . (42) 0 0 0 1 0 0 1 0 The bearing angle αj in (41) can be computed from the known sensor position that receives yj . For example, if the positional measurement yj = [yjx , yjy ]⊤ and the corresponding sensor is located at [sxj , syj ]⊤ , then αj = atan2(yjy − syj , yjx − sxj ). The clutter likelihood ℓ0 (Yj ) in (2) may be defined using a prelearned spatial–Doppler map or any suitable distribution over part or all of Yj = (yj , vj ), analogous to Section II-A. For the Doppler component vj , a common choice is a mixture of a uniform distribution over a prescribed velocity range and a normal distribution centred at zero, placing additional mass on stationary clutter. This completes the formulation of the Doppler-augmented Poisson measurement model. All joint priors and likelihood terms in (4)–(8) remain valid, with the single-object likelihood ℓ replaced by (41). Finally, we note that the convenient linear Gaussian singleobject likelihood (41) in the Doppler-augmented Poisson measurement model relies on the assumption that all detection points from the same object share the same true velocity (see Fig. 4(b)). This is valid for a rigid object moving approximately in a straight line (including accelerated motion), but not for turning motion that induces velocity variation across the rigid body, or for nonrigid objects exhibiting micro Doppler effects. In such cases, the Doppler noise in (41) is expected to absorb part of the resulting model mismatch. More accurate Doppler models that incorporate yaw rate as part of the object state exist [26], [27] and can be extended to the NHPP setting, but these lose the convenient linear Gaussian form of (41). B. Doppler-augmented PiVoT updates Since the two-stage inference framework and the corresponding updates in Sections III to V are formulated for the
10
Fig. 4: (a) Automotive radar data from [8], with PiVoT estimates of moving objects overlaid (ellipses: shape, curves: tracks). Coloured points and arrows denote positional measurements and Doppler velocities from ground truth moving objects, whereas black points and arrows denote static background and clutter. (b) Simulated measurements under the Doppler NHPP model. Ellipses show object shapes Pk . Crosses denote positional measurements yj from clutter (black), stationary objects (yellow) and moving objects (other colours). Measurements from the same object share the same true velocity, shown by blue arrows, while their radial velocities relative to the sensor (yellow star) are indicated by red arrows. One example detection yj is annotated with its bearing angle αj , true velocity components vkx , vky and radial velocity vkj (k being the true generating object’s index). general NHPP likelihood in (2), which includes the Doppleraugmented model in Section VI-A, the derived update rules remain valid with only minor changes to the final expressions. 1) Stage 1 updates: Stage 1 performs CAVI to iteratively update q1 (X), q1 (Λ), q1 (P ), and q1 (θ), which retain the same factorised structure and distributional forms as in (25). The state update still follows the rule in (26), but with single-object likelihood ℓ now given by (41). This yields the following independent state updates for each object k
2) Stage 2 updates: Recall that Stage 2 evaluates the existence and detectability distribution q2 (Ek , Dk ) for each object k, using the q1 (X, Λ, P, θ) obtained in Stage 1. All Stage 2 results in Sections V-B, V-C depend on the singleobject likelihood ℓ only through Sjk . Therefore, the optimal q2∗ (Ek , Dk ), q2∗ (Ek ) in (37), (40), the expression for g(Dk ) in (38), (39), and its numerically stable and efficient evaluations in Appendix F all remain unchanged. We simply note that Sjk is now evaluated using its Doppler-augmented form in (45).
q1 (Xk ) ∝ p̂(Xk ) N y k ; HXk , Rk N v k ; GXk , U k , (43)
C. Properties of Doppler-augmented PiVoT The Doppler-augmented PiVoT enables efficient detection and tracking using all received Doppler measurements through principled model-based inference. We highlight two advantages over common Doppler pruning and clustering heuristics. First, PiVoT associates Doppler measurements using object shape, kinematics, and sensor geometry jointly through inference, rather than grouping similar Doppler values as in distance-based clustering methods [30], [31]. In practice, detections from the same object can exhibit markedly different Doppler velocities due to their relative geometry to the sensor, as seen from the arrows on the blue points in Fig. 4(a) and the green crosses in Fig. 4(b). Distance-based clustering may therefore incorrectly separate them, whereas PiVoT accounts for the geometric model and maintains a consistent association. Second, for detection of moving objects, measurements with Doppler velocity close to zero are still useful and are fully exploited in PiVoT. As illustrated in Fig. 4(b), the purple points have almost zero radial velocity but arise from a moving object, and these measurements provide evidence for the existence of objects moving approximately perpendicular to the sensor’s line of sight. In contrast, a simple preprocessing step in radar pipelines is to prune such near-zero Doppler returns by treating them as stationary clutter.
where p̂(Xk ) is given in (12), and y k and Rk are as defined in (27). The matrices H and G are exemplified in (42). The Doppler related terms U k and v k take the form XM −1 U k = σv2 q1 (θj = k) e⊤ e , (44) j j j=1 XM −1 XM vk = q1 (θj = k) e⊤ q1 (θj = k) e⊤ j ej j vj , j=1
j=1
where σv2 , ej are given in (41). This update in (43) can be implemented by applying two standard Kalman filter updates in sequence, one to incorporate y k as in the positional update in (27), and one to incorporate v k for the Doppler measurements. The rate and shape updates remain unchanged: each q1 (Λk ) still follows (28) and each q1 (Pk ) follows (29), as these do not involve the additional Doppler term N (vj ; ej GXk , σv2 ) that augments the single-object likelihood ℓ in (41) relative to (3). The data association still follows the independent updates in (30) for each q1 (θj ). The clutter likelihood ℓ0 (Yj ) in (30) now incorporates the Doppler measurement. With the Doppleraugmented likelihood ℓ in (41), the Sjk in (30) becomes 2 Sjk = exp(Tjk − 2σ1 2 ej GΣk G⊤ e⊤ j )N (vj ; ej Gµk , σv ), (45) v
where Tjk is defined in (32); vj , σv2 , ej , G are as in the Doppler likelihood (41), and µk , Σk are variational parameters in (25).
VII. A LGORITHMIC R EALISATION OF P I VOT Algorithm 1 summarises the generic PiVoT routine for both positional-only and Doppler-augmented models. This section
11
Algorithm 1: Generic PiVoT at time step n Prediction: Evaluate the predictive prior p̂n in (12) using the prediction step (46), with birth prior (10) set as in points (i) and (iii) of Section IV-D. Stage 1 Initialisation: Initialise q1 according to Appendix G-1. Stage 1 CAVI updates (Section IV-B): while not converged do For k = 1, 2, ..., Kn , update kinematic states q1 (Xn,k ) using (26). For k = 1, 2, ..., Kn , update measurement rates q1 (Λn,k ) using (28). For k = 1, 2, ..., Kn , update object shapes q1 (Pn,k ) using (29). For j = 1, 2, ..., Mn , update data association q1 (θn,j ) using (30). Ineffective PMn birth removal (Section IV-C): if j=1 q1 (θn,j = k) < L for any k = Kn−1 + 1, ..., Kn then Remove components {X, Λ, P, E, D}n,k for all such k from p̂n and q1 , restrict each θn,j to the remaining indices, then update Kn and relabel indices accordingly such that k ∈ {1, ..., Kn }. Stage 2 (Section V): For k = 1, 2, . . . , Kn , compute the globally optimal object existence and detectability distributions q2∗(En,k , Dn,k ), q2∗(En,k ) using (37), (40), with g(Dk ) evaluated from (65)–(68) in Appendix F.
Remark: This generic algorithm covers both positional-only PiVoT and Doppler-augmented PiVoT. In the q1 (Xn,k ) update, (26) reduces to (27) and (43), respectively. Similarly, the Sjk terms (used in (30) and (65)–(66)) reduce to (31) and (45), respectively.
gives implementation details for the prediction step in Algorithm 1; Stage 1 initialisation is detailed in Appendix G. We then describe the extraction of state estimates and a Doppler specific post-processing step for reporting confidently moving objects in automotive radar data. A. Prediction After obtaining q1 (Xn−1 ), q1 (Λn−1 ), q1 (Pn−1 ) in (25) from Stage 1, and q2∗ (En−1,k ) = Ber(pen−1,k ) in (40) from Stage 2 at time step n − 1, the prediction step for parameters in (12) at time step n for legacy object k = 1, ..., Kn−1 is: µ′n,k = Fn µn−1,k , Σ′n,k = Fn Σn−1,k Fn⊤ + Qn , ′ ηn,k = ηn−1,k γΛ,n ,
ρ′n,k = ρn−1,k /γΛ,n ,
ϕ′n,k = max{ϕn−1,k γP,n , dY − 1}, Φ′n,k = Φn−1,k ϕn−1,k /ϕ′n,k ,
(46)
e s pe′ n,k = pn−1,k pn,k ,
where Fn , Qn are from the state transition in (9), γΛ,n , γP,n ≤ 1 are forgetting factors for Λ and P , and psn,k ∈ (0, 1] is the survival probability defined in Section II-B1. The state and existence predictions follow from the linear Gaussian transition and the survival process defined in Section II-B1. The rate and shape predictions preserve the mean while inflating uncertainty via the forgetting factors. Alternative prediction models for Λ and P are also available, e.g. [12], [35]. For a moving platform, such as automotive radar on an ego vehicle, measurements and object states may need to be transformed between time steps to account for ego motion.
Fig. 5: PiVoT tracking results on the RadarScenes [8] dataset (a) before and (b) after moving object post-processing. The ego vehicle is located at (0, 0) and moves to the right. Black and coloured dots/arrows denote measurements from ground truth stationary background and moving objects, respectively. Red trajectories and ellipses in (a) show all estimated tracks and shapes, including stationary and moving objects. All objects shown in (b) are regarded as moving. Among those, orange, cyan, and dashed ellipses indicate tracks rejected due to Doppler uninformativeness, ghost suppression, and the concrete object shape constraint, respectively. previously reported objects will not be reported in the current time step if q2∗ (En,k ) < Pstp . It is important to note that the thresholds Prep , Pstp affect only the final selection of estimates from the set of all processed objects. PiVoT still performs variational inference over the full set of potential objects with the same computational efficiency, regardless of the choice of reporting thresholds. C. Post-processing for moving object selection
B. Estimation extraction To prevent unbounded growth of Kn , a pruning step is applied at each time step, removing objects with existence probability q2∗ (En,k ) < Ppru from the posterior. This is necessary because, in principle, Kn can grow to infinity over time, as every potential object born at previous time steps, including those with zero existence posterior, remains in the system and continues to be processed. For reporting a reliable estimation output at the current time step, PiVoT introduces another two thresholds: when reporting an unreported object for the first time, we require q2∗ (En,k ) > Prep to select only convincing estimates from heavy clutter;
PiVoT processes the full measurement stream to jointly detect and track all persistent structures, including genuinely moving objects, stationary infrastructure, and clutter induced artefacts. For the automotive radar experiments, we therefore apply a lightweight deterministic post-processing stage to report a reliable subset of moving objects. This stage operates only on inferred track statistics and does not modify the variational inference. In brief, tracks are retained only when they show 1) consistent kinematic motion over a short temporal window, 2) informative Doppler support along the estimated direction of motion when required to disambiguate motion, 3) plausible road user shape and sufficient dominance of their
12
Fig. 6: True tracks (black), shapes, one-run PiVoT estimates (red dashed tracks and ellipses), and final-time-step measurements for trajectory/dataset pairs (a) T1/DS1, (b) T2/DS3, (c) T3/DS5, and (d) T4/DS6. True shapes (grey-to-black for earlier-to-later time steps) and shape estimates are shown at multiple steps in (a)–(b), and only at the final step in (c)–(d). TABLE I: Performance comparison on DS1–DS6. Time and GOSPA are averaged over all time steps and 50 Monte Carlo runs† . Each dataset header gives the trajectory set, detection probability pdk , and clutter rate Λ0 . Object numbers, tracks, shapes, and measurement rates are defined by T1–T4, see Section VIII-A1. DS2 additionally sets the T1 measurement rate to Λk = 5. Method
Time (s) GOSPA Location Missed False Time (s) GOSPA Location Missed DS1: T1, pdk = 1, Λ0 = 47.25
PiVoT PMBMF PMBMC SPA3000 SPA15000
PiVoT PMBMC SPA100
0.0077 82.8745 79.6145 1.3000 1.9600 0.4747 108.0524 89.4524 17.200 1.4000 0.1837 109.6324 90.4324 18.400 0.8000 1.6411 133.4078 118.2878 8.1800 6.9400 7.3063 93.0783 91.5983 0.8800 0.6000
False
DS2: T1 with Λk = 5, pdk = 0.9, Λ0 = 47.25 0.0106 0.6243 0.2927 1.3807 6.0615
126.4433 105.1433 8.4400 141.2471 114.7471 12.500 141.5576 114.7576 13.000 170.7621 140.6221 16.260 134.4843 118.8843 8.6000
12.860 14.000 13.800 13.880 7.0000
Time (s) GOSPA Location Missed False DS3: T2, pdk = 0.9, Λ0 = 47.25 0.0091 1.6159 0.6221 5.8976 –
191.1176 159.5576 16.860 14.700 253.0662 190.9662 17.900 44.200 261.1853 193.6853 19.900 47.600 463.1671 252.0071 160.14 51.020 – – – –
DS4: T2, pdk = 1, Λ0 = 1575
DS5: T3, pdk = 1, Λ0 = 472.5
DS6: T4, pdk = 1, Λ0 = 1575
0.2013 223.3302 168.1702 26.540 28.620 9.3477 355.4586 208.2586 120.90 26.300 143.44 869.7601 325.1601 445.70 98.900
0.0550 395.5997 361.5197 24.460 9.6200 4.7251 621.1152 458.9152 115.80 46.400 16.638 1930.081 774.9614 782.28 372.84
0.3008 1000.147 817.8874 160.08 22.180 39.983 1931.305 883.2056 925.00 123.10 216.10 3537.380 1396.480 1655.7 485.20
† The only exceptions are that SPA 100 uses 10 runs on DS4 and DS6, and SPA15000 is omitted on DS3 due to excessive runtime.
associated detections within the estimated extent, and 4) no evidence of being a multipath ghost induced by stationary structures. Fig. 5 illustrates how these checks refine the full PiVoT output into the final reported moving object set. This automotive radar specific stage is used only to report a reliable and interpretable subset of moving road users from the richer set of PiVoT tracks, without changing inference; further procedural details are therefore given in Appendix G-2. VIII. E XPERIMENTS This section first evaluates PiVoT against existing modelbased trackers on various simulated tracking scenarios in Section VIII-A, followed by a large-scale simulation demo tracking 1000 objects in Section VIII-B. Finally, Section VIII-C compares PiVoT with a deep learning detection benchmark on a real-world automotive radar dataset. A. Evaluation on simulated scenes This section compares PiVoT’s detection and tracking performance with existing Bayesian NHPP-based trackers, including the Poisson multi-Bernoulli mixture (PMBM) filter [12] and the sum-product algorithm (SPA) [14]. Since these methods do not yet accommodate Doppler measurements, we restrict this comparison to positional NHPP model in (3). The evaluation uses 6 simulated datasets of increasing complexity.
1) Experimental setting: The six datasets are built from four trajectory sets, T1–T4, shown in Fig. 6. We use two −1 object types: small objects with (Λk , Pn,k ) = (5, 100I2 ) and −1 large objects with (Λk , Pn,k ) = (8, 800I2 ). T1 contains 4– 10 randomly appearing and disappearing large objects. T2 contains 11–25 objects, with 60% small and 40% large. T3 and T4 contain 21–50 and 42–99 objects, respectively, both with 30% small and 70% large. Each trajectory set has 50 time steps with sampling interval τ = 1 s. Objects follow a 2D constant-velocity model with positional NHPP measurements. In (3) and (9), we set Fn = diag(Fn1 , Fn2 ), Qn = diag(Q1n , Q2n ), and Hn = diag(Hn1 , Hn2 ), where for each d = 1, 2, 3 1 τ τ /3 τ 2 /2 Fnd = , Qdn = 25 2 , Hd = 1 0 . 0 1 τ /2 τ (47) For each dataset DS1–DS6, the true tracks are fixed from one of T1–T4, while 50 independent measurement sets are generated from the model, each over 50 time steps. The detailed configurations are given in Table I. DS1–DS2 contain only large objects, and DS3–DS6 contain both object types. The detailed parameter settings for all compared methods are provided in Appendix J. In summary, PiVoT uses uninformative parameterisations across DS1–DS6 by design, to assess robustness across diverse scenes. PMBM and SPA are
13
configured more informatively to improve performance; for example, in DS1–DS2, their shape priors concentrate more on the single object type, while PiVoT uses uninformative priors spanning both types. SPA also uses the ground truth measurement rate, since it is not estimated in [14]. In Table I, PMBMC and PMBMF use 16 and 100 DBSCAN distance values, respectively, evenly spaced over [1, 100], a range found to give strong performance. The subscript in SPA denotes the number of particles, e.g., SPA3000 uses 3000 particles. All methods are implemented in MATLAB, with runtime measured as the average execution time, in seconds, per time step (System: Apple M1, 16GB RAM). Tracking performance is evaluated using the GOSPA metric [36], a standard metric in multi-object tracking literature, with cut-off distance 100, penalty parameter 1, and α = 2. The total GOSPA error decomposes into localisation error for matched objects, missedobject, and false-object errors, as reported in Table I. 2) Results: The results in Table I establish PiVoT as a promising advancement in multi-object tracking, in both accuracy and efficiency. Across all datasets, PiVoT consistently achieves the lowest GOSPA scores and the fastest runtime. Its advantages are especially clear in DS4–DS6, which involve heavy clutter and/or large object numbers. In these challenging cases, PiVoT achieves substantially higher tracking accuracy while being around 50–100 times faster than the second-best method. In the simpler DS1–DS2 scenarios, SPA achieves lower missed object and false object components, but remains second to PiVoT in overall GOSPA and requires more than 500 times the computation time. As the tracking challenge intensifies in DS4–DS6, SPA deteriorates substantially. PMBM then becomes the second-best method, although its accuracy also degrades in highly cluttered scenes with many objects. In contrast, PiVoT maintains robust tracking performance, with GOSPA values increasing roughly proportionally to the object number. SPA’s deterioration in DS4–DS6 may reflect its sensitivity to uncertainty in the shape prior as two object types are present: it is observed that SPA improves clearly when the shape prior is sufficiently concentrated around the true shape. PiVoT’s tracking and shape estimates are visualised in Fig. 6 for representative datasets. The estimated tracks and shapes closely follow the ground truth, with few missed and false detections, even under the heavy clutter and frequent coalescence in scenes (c) and (d). The initially imperfect shape estimates in Fig. 6(a) and (b) are due to the uninformative prior and limited early observations, but are soon corrected to reasonable accuracy. In Fig. 6(d), and more generally in DS6, severe clutter makes true objects with sparse measurements difficult to distinguish by human inspection. A video demonstrating PiVoT’s accurate detection and tracking in this challenging setting is provided in [18].
Fig. 7: True and estimated object number.
Fig. 8: Tracking results for up to 1036 objects. Black curves are true tracks and red dotted curves are PiVoT estimates. The right panel enlarges the region marked by the green rectangle. steps. Objects follow a constant-velocity motion model with parameters in (47) and the true tracks are shown in Fig. 8. −1 Object elliptical shape covariance is Pn,k = 50 I, and the measurement rates for all objects are 5. Detection probability is pdn,k = 1. The clutter density is 10−6 per unit area, resulting in a clutter rate of 15.75. The object number varies as in Fig. 7. The challenge of this case is efficiently handling a large volume of data while managing many birth and death events. Further, the close proximity of objects leads to frequent coalescence, making data association particularly difficult. PiVoT uses the same parameter setting as in Section VIII-A, except that the ground-truth object rates and shapes are assumed known. Fig. 7 and Fig. 8 present the estimated trajectories and estimated object counts over time alongside the ground truth, demonstrating PiVoT’s ability to track a thousand objects simultaneously, with timely detection of births and deaths and no excessive track loss. Furthermore, PiVoT achieves this with remarkable computational efficiency, requiring only 0.67 s per time step on average to track over 1000 objects in MATLAB 2024b (System: Apple M4 Pro, 24GB RAM). This implementation uses only vectorised operations, without explicit parallelisation or any gating technique. A video demonstration of this large-scale scene is provided in [18]. Given its performance, further speed improvements could be expected through dedicated CPU/GPU parallelisation and gating over a large surveillance area, making PiVoT a highly promising solution for real-time multi-object tracking in large-scale, cluttered environments.
B. Demo of large-scale multi-object tracking in clutter This section demonstrates PiVoT’s capability to detect and track 1000 objects in clutter. Since running PMBM and SPA in this large-scale setting would be computationally impractical, we report only PiVoT results as a scalability demonstration. We consider a densely populated tracking environment with up to 1036 objects appearing and disappearing over 50 time
C. Evaluation on RadarScenes dataset To demonstrate PiVoT’s capability to operate on real-world automotive radar data with Doppler information, we evaluate it on RadarScenes [8], a dataset with typical radar clutter and diverse road users. RadarScenes is widely used for assessing machine learning methods for automotive radar-based multi-
14
TABLE II: RadarScenes results on the validation and test splits. Time is measured per 1 s of radar data. All object matching and GOSPA metrics are macro-averaged over frames, except APglobal , which pools object estimates over all frames. Split
Method
Time (s)
Object matching, IoU 0.3 F1
Object matching, IoU 0.5
Precision Recall APmacro APglobal
F1
GOSPA error
Precision Recall APmacro APglobal GOSPA Location Missed False
RadarGNNO 4.895 0.6943 RadarGNNB 4.899 0.7407 Val. RadarGNNP 4.896 0.7886 PiVoT 0.830 0.8016
0.6383 0.7085 0.8291 0.8547
0.8259 0.8279 0.7897 0.7993
0.7729 0.7797 0.7585 0.7825
0.7482 0.7415 0.6916 0.7106
0.6017 0.6447 0.6931 0.7829
0.5551 0.6190 0.7309 0.8313
0.7096 0.7144 0.6892 0.7815
0.6458 0.6550 0.6456 0.7597
0.5901 0.5862 0.5535 0.6739
4.9084 4.2197 3.4500 2.8114
1.5588 1.5153 1.3384 0.8704
0.7112 2.6384 0.8098 1.8946 1.2850 0.8266 1.2400 0.7009
RadarGNNO 5.135 0.6643 RadarGNNB 5.139 0.7109 Test RadarGNNP 5.133 0.7538 PiVoT 0.861 0.7723
0.6157 0.6884 0.8103 0.8372
0.7826 0.7826 0.7396 0.7573
0.7211 0.7302 0.7072 0.7357
0.6921 0.6827 0.6263 0.6748
0.5688 0.6117 0.6571 0.7554
0.5282 0.5933 0.7050 0.8186
0.6669 0.6706 0.6438 0.7407
0.5964 0.6074 0.5992 0.7131
0.5316 0.5255 0.4907 0.6469
5.6744 4.8692 4.0263 3.6252
1.7170 1.6411 1.3731 0.9633
0.9971 2.9603 1.1585 2.0697 1.8314 0.8218 1.7185 0.9434
TABLE III: Point-wise moving foreground segmentation results on RadarScenes. Macro metrics are averaged over frames, while global metrics are computed by pooling all points. Split
Method
F1
Macro Precision Recall
F1
Global Precision Recall
Val.
RadarGNN 0.8258 PiVoT 0.8551
0.7710 0.8824
0.9108 0.8887 0.8541 0.8785
0.8250 0.8880
0.9631 0.8692
Test
RadarGNN 0.7861 PiVoT 0.8119
0.7311 0.8500
0.8801 0.8607 0.8076 0.8410
0.7980 0.8516
0.9342 0.8307
object detection. We therefore compare with RadarGNN [37], a strong deep learning baseline on this dataset with a publicly available implementation and released trained checkpoints. We note that the two methods differ in their design goals, temporal aggregation windows, and input features. PiVoT is evaluated in its full tracking configuration with moving object selection (Section VII-C), where each frame accumulates radar data over 0.25 s. This shorter frame length reduces latency and enables more timely responses in real-time tracking. By contrast, RadarGNN uses frames accumulated over 0.5 s and performs frame-wise detection without temporal state propagation, consistent with its original design and training in [37]. The input information also differs: PiVoT currently uses only positional and Doppler measurements, whereas RadarGNN additionally uses radar cross-section information. RadarGNN’s object classification capability is not considered here, since this experiment focuses on the object detection and tracking performance across all moving objects. We emphasise that this comparison is not intended to claim state-of-the-art performance on RadarScenes, but to demonstrate that a fully training-free, model-based joint detector and tracker can achieve efficient and competitive performance on a full-scale modern automotive radar benchmark. To our knowledge, PiVoT is the first such method to report quantitative results comparable to a deep learning detection benchmark. 1) Experimental setting: We evaluate all methods on the RadarScenes validation and test splits, following [37]. For RadarGNN, we use the best-performing translation-invariant model reported in [37] with the released checkpoints trained on the RadarScenes training split. Besides the original postprocessing configuration, denoted by RadarGNNO , we consider two more configurations. Since RadarGNNO tends to produce many false alarms, we tune the confidence thresholds used to remove background points and objects to obtain a balanced precision-recall configuration, denoted by
RadarGNNB , and a precision-oriented configuration, denoted by RadarGNNP . These configurations are selected based on detection metrics evaluated on the validation split. For PiVoT, the parameterisation is selected only through visual inspection of tracking videos from a few sequences in the validation split. All runtime evaluations are conducted on a standard laptop with a 2.4 GHz 8-core Intel Core i9 CPU and 32GB RAM. PiVoT is implemented in MATLAB, while RadarGNN is run in Python through Docker with access to 7 of the 8 CPU cores. Object-level detection and tracking results are reported in Table II. We use object-matching metrics, including precision, recall, F1 score, and average precision (AP), under point-wise IoU thresholds of 0.3 and 0.5 as in [3], together with GOSPA metrics [36] using cut-off distance 3 m, penalty parameter 1, and α = 2. For GOSPA, unavailable ground-truth object positions are computed by averaging the corresponding groundtruth points. Table III further reports point-wise movingforeground segmentation performance, where all RadarGNN configurations share the same point-wise results. 2) Computational time: Table II reports the average processing time per 1 s radar data interval. For RadarGNN, the full inference chain takes 4.897 s on the Val. split and 5.136 s on the Test split on average. In contrast, PiVoT takes 0.830 s and 0.861 s, respectively, both below the 1 s data interval, satisfying the real-time processing requirement on average and running much faster than RadarGNN under the evaluated CPU setting. We note that RadarGNN is expected to run much faster on a GPU, while PiVoT may also be further accelerated. 3) Detection and Tracking Results: From Table II, PiVoT achieves the best frame-wise macro-averaged overall performance across all RadarGNN configurations, including F1 score, APmacro , and GOSPA. This improvement is mainly due to its higher precision and, similarly, fewer false alarms, which reflects PiVoT’s ability to recover true moving objects from dense clutter. This precision advantage remains even compared with RadarGNNP , whose thresholds are deliberately tuned towards higher precision at the cost of recall, with recall already lower than that of PiVoT. By contrast, RadarGNN often achieves higher APglobal and recall, and similarly fewer missed detections in Table II. This is mainly because PiVoT uses a conservative moving-object selection rule, which requires moving evidence over two to three frames before declaring a moving object. As shown in columns 1, 3, and 4 of Fig. 9, this can cause PiVoT to miss objects that have only recently appeared in the field of view.
15
Fig. 9: Visualisation on representative RadarScenes frames. The results are matched to the camera images. Coloured points denote ground-truth moving-object measurements, with colours indicating object identities. RadarGNNB contains more radar points from longer accumulation. PiVoT shows estimated shapes and tracks using ellipses and curves, while RadarGNNB shows detections using bounding boxes. Dashed cyan ellipses and curves in PiVoT’s first column indicate a ghost target. The same trend is also reflected in Table III, which reports point-wise moving-foreground segmentation between movingobject and background measurements. PiVoT is stronger in precision and macro-averaged F1, whereas RadarGNN has an advantage in recall and global F1. An exception appears in Table II under the stricter IoU threshold of 0.5, where RadarGNN performs worse than PiVoT across all object-level metrics. This is mainly because RadarGNN bounding boxes can miss many ground-truth points and sometimes split the same object into multiple boxes. Both effects are visible in Fig. 9. The third column particularly illustrates poor box alignment when the object is diagonally oriented with respect to the coordinate axes, possibly due to limited representation of such orientations in the training data. Fig. 9 further supports the preceding analysis through representative visual comparisons between RadarGNNB and PiVoT. In particular, PiVoT produces fewer false alarms and better shape estimates in terms of measurement inclusion. Its missed detections mainly occur for newly appearing objects or objects with abnormal measurements. For example, in column 4, the missed object has Doppler arrows inconsistent with its motion, possibly due to upstream radar signal-processing artefacts. In column 1, PiVoT also correctly suppresses a ghost target that RadarGNNB falsely detects as a moving object, highlighting the benefit of retaining stationary structures as evidence for ghost-target rejection (see also Fig. 5). We note that PiVoT also provides multi-object tracks directly, as shown in Fig. 9, along with estimates such as object kinematics and existence probabilities. In contrast, RadarGNN outputs require additional post-processing to form tracks. Finally, [18] provides frame-by-frame comparison videos between PiVoT and RadarGNNB on three complete radar sequences, together with additional PiVoT demonstrations on varied RadarScenes clips.
IX. T HEORETICAL ANALYSIS This section presents Theorem IX.1, which justifies PiVoT’s early identification of ineffective births in Section IV-C. Further details and the proof are provided in Appendices H and I. We first state the modelling assumption and introduce the admissible bound used in the theorem. Assumption 1. The single-object measurement likelihood ℓ in (2) follows (3). The clutter Poisson intensity λc = Λ0 ℓ0 (Yj ) is uniform over the surveillance area. The birth prior p(Xn,k ) in (10) is uninformative in the observation space, i.e. HΣbk H ⊤ = cIdY with c → ∞, where IdY is the identity matrix, and H is in (3). This induces a uniform prior over RdY for HXk . Remark 1. A convenient consequence of this assumption is that the updated state q1 (Xk ) in (25), (27) satisfies HΣk H ⊤ = Rk (see Lemma H.3 in Appendix H-2) and Hµk = y k with Rk , y k given in (27). This greatly simplifies the state update for newly-born objects, as used in Section IV-C’s Model Scenario. Definition 1 (Admissible bound). For any s > 0 and any newly born object k, define the birth-common constants c(s) := exp(0.5[ψdY ((ϕ′k + s)/2) + log |Φ′k | − dY log π]) d ′ 1 × pe′ k pk ρk exp(ψ(ηk + s)) λc ,
t(s) :=c(s) exp(− d2sY ),
(48) (49)
′ where ϕ′k , Φ′k , pe′ k , ηk are the (birth) prior parameters in (12); dY , ψdY , ψ are as in (32); ρk is given in (28) and fixed across CAVI iterations; and λc is the clutter intensity in Assumption 1. For a newly born object k, a scalar Bk (s) is called an admissibleP bound if, for any current association distribution M satisfying j=1 q1 (θj = k) < s, one round of CAVI updates (29), (28), and (27) produces parameters µk , ϕk , Φk such that
∆k (s) ≤ Bk (s)
(50)
16
where ∆k (s) is defined as ∆k (s) :=
XM
X. C ONCLUSION
1 + t(s)
j=1 dk (yj ) + t(s)
,
(51)
dk (yj ) := exp(0.5(yj − Hµk )⊤ϕk Φk (yj − Hµk )),
(52)
Remark 2. The universal choice Bk (s) = M is always admissible, since dk (yj ) ≥ 1 for all j. Tighter admissible bounds are possible; Appendix I-2 gives one such Bk (s) for the known-shape case. We now present the main theorem. Theorem IX.1. Let Assumption 1 hold. For any s > 0 and any newly born object k = Kn−1 +1, . . . , Kn , let Bk (s) be an admissible bound in Definition 1, with c(s) defined in (48). Denote by W−1 : [−1/e, 0) → (−∞, −1] the −1 branch of the Lambert W function, and define ( /2 k (s) exp(− BdYk (s) )<e s, if c(s)B dY /2 (53) Lk (s) := min{s, Vk (s)}, otherwise Vk (s) :=
1 1 1 Bk (s) − dY /2 W−1
dY /2 − c(s)B exp k (s)
dY /2 Bk (s)
(54)
Suppose CAVI iteratively PMapplies the updates (30)–(27). If, at any CAVI iteration, j=1 q1 (θj = k) < Lk (s), then PM j=1 q1 (θj = k) decreases monotonically to 0 in subsequent iterations, with consecutive updates satisfying PM new ∆k (s) j=1 q1 (θj = k) ≤ 1, (55) < PM old Bk (s) j=1 q1 (θj = k) where ∆k (s) is evaluated from (51) after the corresponding CAVI round starting from q1old (θ). Remark 1. A larger Lk (s) enables earlier identification of guaranteed ineffective births, improving inference efficiency. Appendix I-1 shows that a smaller admissible bound Bk (s) only increases the threshold Lk (s), or leaves it unchanged. Hence a tighter Bk (s) is preferred when available; see Remark 2 under Definition 1 for admissible choices of Bk (s). Remark 2. Lk (s) and Vk (s) remain fixed across CAVI iterations for fixed s and a chosen admissible bound Bk (s). If Bk (s) is also independent of k, the same threshold Lk (s) can be used for all birth components sharing the same birth prior and for all iterations. More refined control is possible by recomputing Lk (s) with different choices of s or Bk (s). Remark 3. The threshold Lk (s) in Theorem IX.1 provides only a sufficient condition for monotonic decrease. In practice, we use a simple heuristic threshold L, such as 0.5 for higher inference efficiency or 0.05 for more reliable identification of ineffective births. This heuristic is convenient, although an overly large L may risk discarding valid births. Further discussion on how the threshold Lk (s) and related functions depend on Bk (s) and s, together with a tighter construction of Bk (s), is provided in Appendix I. The proof of Theorem IX.1, supporting lemmas, and detailed CAVI update analysis are provided in Appendix H.
We present PiVoT, a fast, training-free framework for largescale multi-object detection and tracking in heavy clutter from point clouds. Experiments show that PiVoT significantly outperforms existing Bayesian trackers in cluttered or largescale scenes and scales to a thousand objects without gating. PiVoT runs efficiently on a full-scale automotive radar dataset, with performance comparable to a deep-learning detection benchmark, while requiring no training and additionally providing tracking, shape, and existence probability estimation. For automotive radar applications, future work includes a formal Bayesian treatment of moving-object indication, 4D Doppler radar modelling, the incorporation of radar crosssection information, and object classification. More broadly, PiVoT could be extended towards decentralised sensor fusion [38] and nonlinear Gaussian measurement models. R EFERENCES [1] H. Wang and Y. Tian, “Sequential point clouds: A survey,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 46, no. 8, pp. 5504–5523, 2024. [2] K. Granström and M. Baum, “A tutorial on multiple extended object tracking,” Authorea Preprints, 2022. [3] N. Scheiner, F. Kraus, N. Appenrodt, J. Dickmann, and B. Sick, “Object detection for automotive radar point clouds–a comparison,” AI Perspectives, vol. 3, no. 1, p. 6, 2021. [4] A. Cavagna, S. Melillo, L. Parisi, and F. Ricci-Tersenghi, “SpaRTA tracking across occlusions via partitioning of 3D clouds of points,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 43, no. 4, pp. 1394–1403, 2019. [5] S. S. Kari, A. A. B. Raj, and B. K, “Evolutionary developments of today’s remote sensing radar technology—right from the telemobiloscope: A review,” IEEE Geoscience and Remote Sensing Magazine, vol. 12, no. 1, pp. 67–107, 2023. [6] B. I. Ahmad, C. Rogers, S. Harman, H. Dale, M. Jahangir, M. Antoniou, C. Baker, M. Newman, and F. Fioranelli, “A review of automatic classification of drones using radar: key considerations, performance evaluation, and prospects,” IEEE Aerospace and Electronic Systems Magazine, vol. 39, no. 2, pp. 18–33, 2023. [7] N. Scheiner, F. Weishaupt, J. F. Tilly, and J. Dickmann, “New challenges for deep neural networks in automotive radar perception: An overview of current research trends,” Automatisiertes Fahren 2020: Von der Fahrerassistenz zum autonomen Fahren 6. Internationale ATZFachtagung, pp. 165–182, 2021. [8] O. Schumann, M. Hahn, N. Scheiner, F. Weishaupt, J. F. Tilly, J. Dickmann, and C. Wöhler, “RadarScenes: A real-world radar point cloud data set for automotive applications,” in 2021 IEEE 24th International Conference on Information Fusion (FUSION). IEEE, 2021, pp. 1–8. [9] Y. Bar-Shalom and X.-R. Li, Multitarget-multisensor tracking: principles and techniques. YBs Storrs, CT, 1995, vol. 19. [10] I. J. Cox and S. L. Hingorani, “An efficient implementation of Reid’s multiple hypothesis tracking algorithm and its evaluation for the purpose of visual tracking,” IEEE Transactions on pattern analysis and machine intelligence, vol. 18, no. 2, pp. 138–150, 1996. [11] K. Gilholm, S. Godsill, S. Maskell, and D. Salmond, “Poisson models for extended target and group tracking,” in Signal and Data Processing of Small Targets 2005, vol. 5913. SPIE, 2005, pp. 230–241. [12] K. Granström, M. Fatemi, and L. Svensson, “Poisson multi-Bernoulli mixture conjugate prior for multiple extended target filtering,” IEEE Transactions on Aerospace and Electronic Systems, vol. 56, no. 1, pp. 208–225, 2019. [13] Y. Xia, Á. F. Garcı́a-Fernández, F. Meyer, J. L. Williams, K. Granström, and L. Svensson, “Trajectory PMB filters for extended object tracking using belief propagation,” IEEE Transactions on Aerospace and Electronic Systems, vol. 59, no. 6, pp. 9312–9331, 2023. [14] F. Meyer and J. L. Williams, “Scalable detection and tracking of geometric extended objects,” IEEE Transactions on Signal Processing, vol. 69, pp. 6283–6298, 2021.
17
[15] Q. Li, R. Gan, and S. Godsill, “A scalable Rao-Blackwellised sequential MCMC sampler for joint detection and tracking in clutter,” in 2023 26th International Conference on Information Fusion (FUSION). IEEE. [16] Q. Li, R. Gan, J. Liang, and S. J. Godsill, “An adaptive and scalable multi-object tracker based on the non-homogeneous Poisson process,” IEEE Transactions on Signal Processing, vol. 71, pp. 105–120, 2023. [17] A. Bole, A. Wall, and A. Norris, “Chapter 3 - target detection,” in Radar and ARPA Manual, 3rd ed. Butterworth-Heinemann, 2014, pp. 139–213. [18] PiVoT project page and demos: https://runzegan.github.io/projects/pivot/. [19] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” Journal of the American statistical Association, vol. 112, no. 518, pp. 859–877, 2017. [20] C. Bishop, Pattern Recognition and Machine Learning. Springer, 2006. [21] R. Gan, Q. Li, and S. Godsill, “A variational Bayes association-based multi-object tracker under the non-homogeneous Poisson measurement process,” in 2022 25th International Conference on Information Fusion (FUSION). IEEE, 2022, pp. 1–8. [22] R. Gan, Q. Li, and S. J. Godsill, “Variational tracking and redetection for closely-spaced objects in heavy clutter,” IEEE Transactions on Aerospace and Electronic Systems, vol. 60, no. 4, pp. 5286–5311, 2024. [23] R. D. Turner, S. Bottone, and B. Avasarala, “A complete variational tracker,” Advances in Neural Information Processing Systems, 2014. [24] R. A. Lau and J. L. Williams, “A structured mean field approach for existence-based multiple target tracking,” in 2016 19th International Conference on Information Fusion (FUSION). IEEE, 2016. [25] S. J. Davey and D. A. Gray, “Integrated track maintenance for the PMHT via the hysteresis model,” IEEE transactions on Aerospace and Electronic Systems, vol. 43, no. 1, pp. 93–111, 2007. [26] C. Knill, A. Scheel, and K. Dietmayer, “A direct scattering model for tracking vehicles with high-resolution radars,” in 2016 IEEE Intelligent Vehicles Symposium (IV). IEEE, 2016, pp. 298–303. [27] K. Thormann and M. Baum, “Incorporating range rate measurements in ekf-based elliptical extended object tracking,” in 2021 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI). IEEE, 2021, pp. 1–6. [28] A. Scheel and K. Dietmayer, “Tracking multiple vehicles using a variational radar model,” IEEE Transactions on Intelligent Transportation Systems, vol. 20, no. 10, pp. 3721–3736, 2018. [29] J. F. Tilly, S. Haag, O. Schumann, F. Weishaupt, B. Duraisamy, J. Dickmann, and M. Fritzsche, “Detection and tracking on automotive radar data with deep learning,” in 2020 IEEE 23rd International Conference on Information Fusion (FUSION). IEEE, 2020, pp. 1–7. [30] N. Scheiner, N. Appenrodt, J. Dickmann, and B. Sick, “A multistage clustering framework for automotive radar data,” in 2019 IEEE Intelligent Transportation Systems Conference (ITSC). IEEE, 2019, pp. 2060–2067. [31] C. Malzer and M. Baum, “Constraint-based hierarchical cluster selection in automotive radar data,” Sensors, vol. 21, no. 10, p. 3410, 2021. [32] R. Gan, Q. Li, J. R. Hopgood, M. E. Davies, and S. Godsill, “PiVoT: Poisson measurements-based variational multi-object detection and tracking,” in 2025 28th International Conference on Information Fusion (FUSION). IEEE, 2025, pp. 1–8. [33] A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes, 4th ed. McGraw-Hill, 2002, Section 5-4. [34] R. Gan and S. J. Godsill, “Conditionally factorized variational Bayes with importance sampling,” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022. [35] J. W. Koch, “Bayesian approach to extended object and cluster tracking using random matrices,” IEEE Transactions on Aerospace and Electronic Systems, vol. 44, no. 3, pp. 1042–1059, 2008. [36] A. S. Rahmathullah, Á. F. Garcı́a-Fernández, and L. Svensson, “Generalized optimal sub-pattern assignment metric,” in 2017 20th International Conference on Information Fusion (FUSION). IEEE, 2017, pp. 1–8. [37] F. Fent, P. Bauerschmidt, and M. Lienkamp, “RadarGNN: Transformation invariant graph neural network for radar-based perception,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2023, pp. 182–191. [38] Q. Li, R. Gan, and S. J. Godsill, “Decentralised variational inference frameworks for multi-object tracking on sensor networks,” IEEE Transactions on Signal Processing, vol. 73, pp. 2753–2767, 2025.
1
Supplementary Materials for PiVoT: A Variational Solution for Real-time Large-scale Multi-object Detection and Tracking under Heavy Clutter A PPENDIX A N OTATION Tables IV–VI summarise the notation for variables, parameters, probability laws, and functions used throughout the paper. TABLE IV: Main variable and parameter notation. Note that the time index subscript n is omitted from Section III onwards to simplify notation, unless needed to avoid ambiguity. Notation n k, j Kn b Kn Mn Yn,j , Yn = [Yn,1 , . . . , Yn,Mn ] yn,j dY RdY IdY Xn,k , Xn = [Xn,1 , . . . , Xn,Kn ] Λn,k , Λn = [Λn,1 , . . . , Λn,Kn ] Pn,k , Pn = [Pn,1 , . . . , Pn,Kn ] En,k , En = [En,1 , . . . , En,Kn ] Dn,k , Dn = [Dn,1 , . . . , Dn,Kn ] Λn,0 θn,j , θn = [θn,1 , . . . , θn,Mn ] Xk− , Λk− , Pk− Ek− , Dk− Fn , Qn pdn,k psn,k , pbn,k µbn,k , Σbn,k b , ρb ηn,k n,k Φbn,k , ϕbn,k pe′ n,k µ′n,k , Σ′n,k ′ , ρ′n,k ηn,k Φ′n,k , ϕ′n,k γΛ,n , γP,n Xkp , Xkv L H G vn,j σv2 αn,j , en,j j vn,k x , vy vn,k n,k ȳk R̄k ūk Ūk Ppru , Prep , Pstp e
Description Discrete time index. Object index and measurement index, respectively. Number of potentially existing objects at time n. Number of birth objects at time n. Number of received measurements at time n. The j-th measurement and the full measurement set at time n. Positional component of Yn,j in the positional-only and Doppler-augmented models. Dimension of the positional measurement space. dY -dimensional real vector space Identity matrix of dimension dY . Kinematic state of object k and the collection of all object states. Poisson rate of object k and the collection of all object rates. Shape precision matrix of object k and the collection of all object shape precision matrices. Existence indicator of object k and the collection of all existence indicators. Detectability indicator of object k and the collection of all detectability indicators. Clutter Poisson rate. Association variable for measurement Yn,j and the collection of all measurement oriented associations. Collections of all object states, rates, and shape precision matrices except those of object k. Collections of all existence and detectability indicators except those of object k. State transition matrix and covariance for legacy objects at time n; see (9). Detection probabilities of object k at time n; see (1). Survival and birth probabilities of object k at time n; see Section II-B1 and (10). Mean and covariance of the Gaussian birth prior for Xn,k ; see (10). Shape and scale parameters of the Gamma birth prior for Λn,k ; see (10). Scale matrix and degrees of freedom of the Wishart birth prior for Pn,k ; see (10). Predictive prior existence probability of object k; see (12). Mean and covariance of the Gaussian predictive prior for Xn,k ; see (12). Shape and scale parameters of the Gamma predictive prior for Λn,k ; see (12). Scale matrix and degrees of freedom of the Wishart predictive prior for Pn,k ; see (12). Forgetting factors for the rate and shape prediction step; see (46). Positional and remaining components of the kinematic state Xk , respectively; used in Section IV-D. Practical threshold for early removal of ineffective births. Observation matrix mapping the object state to the positional measurement space, as exemplified in (42). Matrix that extracts the velocity component from the object state, as exemplified in (42). Scalar Doppler velocity component of Yn,j in the Doppler-augmented model, where Yn,j = (yn,j , vn,j ). Variance of the scalar Doppler velocity measurement noise in the Doppler-augmented model; see (41). Bearing angle of measurement Yn,j and the corresponding line of sight unit vector; see (41) and Fig. 4(b). True radial velocity of object k at the bearing of measurement j; see Fig. 4(b). Cartesian velocity components of object k along the x and y axes in the Doppler-augmented model; see Fig. 4(b). Pseudo measurement for object k constructed from all positional measurements for state update; see (27). Covariance matrix associated with the positional pseudo measurement ȳk ; see (27). Pseudo measurement for object k constructed from all Doppler measurements for state update; see (44). Covariance matrix associated with the Doppler pseudo measurement ūk ; see (44). Existence probability thresholds for pruning, reporting, and stopping reporting an object, respectively; see Section VII-B. Euler’s number, the base of the natural exponential, used in Theorem IX.1.
2
TABLE V: Probability law and variational notation. The time index subscript n is omitted from Section III onwards for clarity. Notation p p̂n p̄ q1 q2 q2∗
Description Exact probability law specified by the generative model; see Sec. II. Filtering probability law at time n, formed by combining the predictive prior constructed from the filtered result at time n − 1 with the exact current measurement likelihood; see (11). Approximate filtering probability law obtained after applying the approximation in (22); defined separately for Stage 1 and Stage 2 in (23) and (33). Stage 1 variational law used for detection, tracking, rate estimation, shape estimation, and association. Stage 2 object-wise variational law used for existence and detectability evaluation. Optimal Stage 2 variational law under the corresponding Stage 2 objective.
TABLE VI: Main named functions. Notation N (·) G(·) W(·) Ber(·) δ[·] ψ(·) ψdY (·) exp(·) |·| Tr(·) W−1 (·) ℓ0 (Yj ) ℓ(Yj |Xk , Pk ) h(Λn , Dn , θn ) Ξk (Λk , θ) Sjk g(Dk ) c(s), t(s) dk (yj ) ∆k (s) Bk (s) Lk (s) Vk (s)
Description Gaussian distribution, parameterised by a mean and covariance. Gamma distribution, parameterised by shape and scale parameters. Wishart distribution, parameterised by a scale matrix and degrees of freedom. Bernoulli distribution, parameterised by its success probability. Kronecker delta indicator, equal to 1 when the condition is true and 0 otherwise. Digamma function. Multivariate digamma function of dimension dY . Exponential function. Determinant of a square matrix. Trace of a square matrix. The −1 branch of the Lambert W function. Clutter measurement likelihood. Single-object measurement likelihood, defined separately for the positional-only and Doppler-augmented models. Product term in the association prior. Object specific term obtained after marginalising Dk in the association prior. Expected contribution of object k to measurement Yj in the association update. Detectability dependent term used in the Stage 2 update of q2 (Ek , Dk ). Birth common constants used to construct the admissible threshold. Exponential distance score between measurement yj and object k. Aggregate score controlling contraction of the association mass of object k. Admissible upper bound on ∆k (s) for birth object k. Sufficient threshold of birth object k for association mass to decay to zero. Auxiliary threshold function used in the construction of Lk (s).
Defined in Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Standard notation Sections II-A and VI-A (3) and (41), respectively. (7), (8) (21) (30) (38) (48) and (49) (52) (51) (50) (53) (54)
3
A PPENDIX B L IMITATIONS OF THE NAIVE M EAN -F IELD A SSUMPTION This appendix examines the inevitable failure of using the naive mean-field structure in (13) for the variational distribution q(θ, X, Λ, P, E, D) to approximate the true posterior p̂(θ, X, Λ, P, E, D|Y ) in (11) by minimising KL(q∥p̂). Applying the standard CAVI update formula to this variational inference problem, we can verify that the mean-field structure q in (13) needs to satisfy the following finer mean-field structure to achieve a lower KL(q∥p̂): q(θ, X, Λ, P, E, D) =
M Y
q(θj )
j=1
K Y
q(Xk )q(Λk )q(Pk )q(Ek , Dk)
k=1
(56) A fundamental issue arises because the true posterior p̂(θ, X, Λ, P, E, D|Y ) in (11) involves p(θ, M |Λ, D) in (6) and p(Dk |Ek ) in (1) as multiplication factors. Specifically, for any k = 1, ..., K and j = 1, ..., M , the joint factor p(θj = k, θj− , M |Λ, D)p(Dk |Ek = 0) is always zero. Consequently, p̂(θj = k, Ek = 0, θj− , Ek− , X, Λ, P, D) is always zero, i.e. p̂ is zero whenever a measurement Yj is associated with a non-existent object k. Since a variational distribution cannot assign probability mass where p̂ = 0 (as otherwise KL(q∥p̂) would be infinite), we must have: q(θj = k, Ek = 0, θj− , Ek− , X, Λ, P, D) = 0, for all j, k, and for all values of θj− , Ek− , X, Λ, P, D. Since q follows the structure in (56), this constraint implies
1) Point process interpretation of the approximation: Recall from Section II-A that the overall measurement process is constructed by superposing a clutter process with intensity Λ0 ℓ0 (Yj ) and K independent NHPPs, where the k-th NHPP has intensity Dk Λk ℓ(Yj |Xk , Pk ), dependent on the detectability variable Dk . From (4)-(7), the joint association and likelihood density for the overall measurement process is QK p(Y, M, θ|D, X, Λ, P ) = M1 ! p(Y |θ, P, X) Ω0 s=1 Ωs (58) PM
Ω0 = e−Λ0 Λ0 j=1
δ[θj =0]
PM
, Ωs = e−Ds Λs (Ds Λs )
where each Ωs is contributed solely by the s-th NHPP and is independent of the other superposed components. Now consider marginalising Dk under prior p̂(Dk ), similarly to (16) and (19). The global joint density P p̂(D k )p(Y, M, θ|D, X, Λ, P ) retains the same form as Dk (58), except that Ωk is replaced by its expectation under p̂(Dk ), i.e. Ep̂(Dk ) Ωk = Ξk (Λk , θ), which is exactly the quantity defined in (21). From a modelling perspective, this corresponds to replacing the k-th component of the superposed measurement process, which was originally an NHPP conditioned on Dk , with a marginalised point process. The other components remain unchanged due to independence. To examine this Dk -marginalised point process more closely, let Z = [Z1 , . . . , Zm ] denote the measurement set generated by object k alone, with cardinality m. Its likelihood density pk (Z, m|Xk , Λk , Pk ) is obtained by integrating out Dk from the likelihood of the original k-th NHPP, i.e. Ym pk (Z, m|Xk , Λk , Pk ) = pk (m|Λk ) ℓ(Zj |Xk , Pk ), j=1
pk (m|Λk ) := Ep̂(Dk ) Pois(m; Dk Λk ). q(θj = k)q(Ek = 0, Dk )q(θj− , Ek− , X, Λ, P, Dk− ) = 0. (57) To allow an object k even a small probability of associating with at least one measurement, i.e. q(θj = k) ̸= 0, the only way to ensure the constraint (57) for all possible values of θj− , Ek− , X, Λ, P, Dk− and Dk is to enforce q(Ek = 0, Dk ) = 0 for all values of Dk , which in turn implies q(Ek = 0) = 0. This leads to an inevitable failure: 1) Forcing q(Ek = 0) = 0 means the object k must always exist, eliminating any meaningful uncertainty quantification of existence. 2) Otherwise, the object k cannot associate with any measurement, preventing it from ever being updated. Either outcome is problematic, making the naive mean-field assumption in (13) fundamentally flawed. A PPENDIX C O PTIMAL NHPP A PPROXIMATION V IEW OF (22) This appendix provides a detailed interpretation of the approximation in (22) from the perspective of point process approximation. We first show that approximating Ξk (Λk , θ) in (21) by (22) is equivalent to replacing a marginalised point process with an NHPP. We then justify this replacement by showing that the approximating NHPP minimises the KLD from the original marginalised process.
j=1 δ[θj =s]
(59)
where Pois(m; Dk Λk ) denotes the Poisson distribution with d rate Dk Λk , and p̂k (Dk ) = Ber(pe′ k pk ) as defined in (17). This marginalised process is generally not an NHPP. Rather, it is a mixture between an empty process and an NHPP, with the original spatial density ℓ(·|Xk , Pk ) retained. The approximation replaces Ωk in (58) by e′ d
P Mn
d e−pk pk Λk (pe′ k pk Λk )
j=1 δ[θj =k]
.
Comparing this expression with the original form of Ωk in (58), we see that this corresponds to replacing the k-th NHPP d with rate Dk Λk by a new NHPP with rate pe′ k pk Λk . Since this substitution does not alter the spatial distribution encoded in the likelihood term p(Y |θ, P, X) in (58), this approximating NHPP inherits the same spatial distribution ℓ(·|Xk , Pk ). Hence, the intensity function of the approximating NHPP is: d pe′ k pk Λk ℓ(·|Xk , Pk ). Consequently, replacing Ξk (Λk , θ) in (21) with the approximation in (22) is equivalent to replacing the Dk -marginalised point process (i.e. the k-th component superposed into the overall measurement process) with an NHPP of intensity d pe′ k pk Λk ℓ(·|Xk , Pk ). 2) Justification as KLD minimisation: We justify the employed NHPP approximation by showing that it arises naturally as the solution to a KLD minimisation problem. Specifically, we consider the following result:
,
4
Lemma C.1. Let q(Z, m) denote the likelihood density of an NHPP. Among all NHPP likelihoods, the q(Z, m) that minimises the KLD from the Dk -marginalised point process, i.e. KL(pk (Z, m|Xk , Λk , Pk )||q(Z, m)),
(60)
where pk (Z, m|Xk , Λk , Pk ) is defined in (59), is given by q(Z, m) = Pois(m; Λk Ep̂(Dk) Dk )
Ym j=1
ℓ(Zj |Xk , Pk ). (61)
Proof. As a general form, the likelihood of an NHPP can be expressed as q(Z, m) = Pois(m; λ)
Ym j=1
qℓ (Zj ),
(62)
where λ > 0 is the rate, and qℓ (Zj ) is the spatial distribution. Subsequently, the KLD KL(pk ||q) in (60) can be expressed as KL(pk ||q) = KL(pk (m|Λk )||Pois(m; λ)) (63) Xm + Epk (m|Λk ) KL(ℓ(Zj |Xk , Pk )||qℓ (Zj )). j=1
The second term is minimised to zero by setting qℓ (·) = ℓ(·|Xk , Pk ), which does not affect the first term. Therefore, the optimal spatial distribution is qℓ (Zj ) = ℓ(Zj |Xk , Pk ). The remaining KLD minimisation task is to find the λ that minimises KL(pk (m|Λk )||Pois(m; Λ)). Since the Poisson distribution belongs to the exponential family, its KLD is minimised by matching the expectation of its sufficient statistic, which is m. That is, the optimal Λ satisfies Epk (m|Λk ) T (m) = EPois(m;Λ) T (m) = λ. Using the definition of pk (m|Λk ) in (59), we compute X λ= Ep̂(Dk ) Pois(m; Dk Λk )m = Λk Ep̂(Dk) Dk . m
Substituting this λ and qℓ (Zj ) = ℓ(Zj |Xk , Pk ) into (62) yields the (61) stated in the lemma, completing the proof. Using Lemma C.1 and the definition of p̂(Dk ) in (17), we find that the NHPP that minimises the KLD in (60) from the Dk -marginalised point process has intensity d pe′ k pk Λk ℓ(·|Xk , Pk ), which exactly matches the employed approximate NHPP. The resulting approximation in (22) for Ξk (Λk , θ) is hence justified. A PPENDIX D A PPROXIMATION ACCURACY A NALYSIS AND R EFINEMENT FOR I SOLATED O BJECTS This appendix provides justification for why the approximation introduced in (22) has minimal impact on the inference accuracy of an object k that is well-separated from others, and also explains how to remove this approximation error through an additional refinement step in Stage 2.
1) Justification for well-separated object k: Specifically, we show that the approximation error does not affect the Stage 2 inference for object k where q1 (X, Λ, P ) is fixed. The argument is to show that, under the same gating reduction that is safe for a sufficiently isolated object, the Stage 2 update without the approximation in (22) and the employed Stage 2 update with the approximation in (22) yield the same expression for g(Dk ). It is well known that for an object k sufficiently isolated from others, gating techniques can be employed to restrict consideration to only a subset of measurements that lie within a sufficiently large gate around this object. Denote the set of labels for these measurements as Gk . In this scenario, associations can be safely assumed to involve either object k or clutter, as the remaining objects are too distant to contribute meaningfully. The Stage 2 updates derived in (37) can then be applied for this single-object case, and the optimal log q2∗ (Ek , Dk ) is given by (37), where g(Dk ) takes the form: X log Dk Sjk +Λ0 ℓ0 (Yj ) . (64) g(Dk ) = −ηk ρk Dk + j∈Gk
Note that this expression involves no approximation via Ξ Qk (Λk , θ) from (22), as only object k is considered and the k− Ξk (Λk , θ) in (20), to which the approximation applies, is absent. For comparison, consider the derived expression for log q2∗ (Ek , Dk ) in the multi-object case with Kn objects, where the approximation in (22) is applied. The corresponding g(Dk ) becomes: g(Dk ) = −ηk ρk Dk +
M X j=1
Kn X k d s log Dk Sj +Λ0 ℓ0 (Yj )+ pe′ s ps Sj . s=1 s̸=k
For j ∈ / Gk , the log-summands can be ignored as constants independent of Dk since Sjk ≈ 0. This is due to the very large value of (Yj −Hµk )⊤ϕk Φk (Yj −Hµk ) in the definition of Sjk in (30), as measurements outside the gate are far from object k. Similarly, for j ∈ Gk , Sjs ≈ 0 for all s ̸= k, since those objects s P are far from the measurements within the gate. e′ d s Consequently, s̸=k ps ps Sj can be neglected. Neglecting these terms, the resulting expression for g(Dk ) under the approximation coincides with the exact form in (64), showing that the approximation in (22) does not impact the Stage 2 inference for object k. Hence, the remaining influence of the approximation comes from the q1 (X, Λ, P ) computed from the Stage 1 inference. 2) Approximation mitigation and Stage 2 refinement: The approximation error from (22) introduced in Stage 1 for q1 (X, Λ, P ) can be mitigated or removed by extending the Stage 2 objective to include a refinement step. Specifically, we may incorporate Xk , Λk , Pk in q2 and exclude them from the fixed q1 . In doing so, we continue to approximate the same target distribution p̂ as in (18), but with a refined and richer variational form: q2 (Ek , Dk )q2 (θ|Ek , Dk )q2 (Xk )q2 (Λk )q2 (Pk ) × q1 (Xk− )q1 (Λk− )q1 (Pk− ). In this refined Stage 2 objective, the detectability variable Dk of the object being evaluated is retained explicitly. Hence,
5
unlike in Stage 1 where the approximation in (22) is applied to every object, this approximation is not applied to the objectspecific contribution of k. The approximation can only enter through the marginalised components of the other objects. For a sufficiently well-separated object, these other-object contributions are negligible within the local gate of object k, following the same reasoning as above. Thus, the refinement can reduce, and in the isolated gated case remove, the residual impact of the Stage 1 approximation on the inferred object features. Importantly, this refinement procedure is flexible. One may choose to refine only a subset of variables or a subset of objects, rather than updating all of Xk , Λk , and Pk for every k. Furthermore, structured dependencies can be incorporated, such as q2 (Xk |Ek ), where q2 (Xk |Ek = 1) can approximate a more accurate conditional posterior p̂(Xk | Ek = 1, Y ), as advocated in works such as [14], [24]. Nevertheless, all such refinement steps require additional computation and iterative updates, and are therefore omitted in this paper.
A PPENDIX E D ERIVATION OF EFFICIENT EVALUATION OF q2∗ (Ek , Dk ) Here we derive the optimal q2∗ (Ek , Dk ) in (37) and (38), which enables an efficient Stage 2 implementation. First, observing from (35), the normalised q2∗ (θ|Ek , Dk ) is
PKn
k=0 f (bk )δ[θn,j = k] for f equal to exp or log. Therefore, we arrive at c
log q2∗ (Ek , Dk ) = log p̂(Ek )p(Dk |Ek ) − Dk Eq1 (Λ) Λk YM XKn h + log exp(Eq1 (X,P ) log p(Yj |θj , P, X)) j=1 θj =0 XKn i e′ d × Λ0 δ[θj = 0]+ Λ̃k Dk δ[θj = k]+ p p Λ̃ δ[θ = s] s j s=1 s s s̸=k
where Λ̃k := exp(Eq1 (Λk ) log Λk ) for all k = 1, ..., Kn . Finally, performing the summation over θj , substituting p(Yj |θj , P, X) in (5) and the expectation Eq1 (Λ) Λk using (25), we obtain the final expression for the optimal q2∗ (Ek , Dk ): c
log q2∗ (Ek , Dk ) = log p̂(Ek )p(Dk |Ek ) + g(Dk ), g(Dk ) = −ηk ρk Dk +
M X
Kn X d s log Dk Sjk +Λ0 ℓ0 (Yj )+ pe′ s ps Sj
s=1 j=1 s̸=k k where Sj is defined in (30). This yields the log q2∗ (Ek , Dk )
expression in (37). Since each logarithm differs by only one summand for different k, we simplify g(Dk ) using a shared total sum: M X d k log Sjsum + (Dk − pe′ g(Dk ) = −ηk ρk Dk + k pk )Sj j=1
Sjsum = Λ0 ℓ0 (Yj ) +
Kn X
d k pe′ k pk Sj .
k=1
exp(Eq1 (X,Λ,P ) log p̄(X, Λ, P, θ, Ek , Dk , Y )) P . θ exp(Eq1 (X,Λ,P ) log p̄(X, Λ, P, θ, Ek , Dk , Y )) Substituting this into the log q2∗ (θ|Ek , Dk ) term in (36), the last line of (36) cancels, yielding log q2∗ (Ek , Dk ) up to an additive constant as
This matches (38) and completes the derivation. A PPENDIX F N UMERICALLY STABLE COMPUTATION OF g(Dk )
While (38) and (39) offer an efficient way to compute g(Dk ), they may encounter numerical issues in practice. This P appendix presents a numerically stable yet efficient evaluation log θ exp(Eq1 (X,Λ,P ) log p̄(X, Λ, P, θ, Ek , Dk , Y )). of g(Dk = 1) and g(Dk = 0) that mitigates overflow/underflow, and quantifies its computational cost. For all j = 1, ..., M Next, substituting (33), (5), and discarding terms constant and k = 1, ..., Kn , the following quantities are assumed to be ∗ with respect to Ek , Dk , the optimal log q2 (Ek , Dk ) is precomputed and stored with sufficient precision. These values are also required in the log-sum-exp computation of Sjsum in c log q2∗ (Ek , Dk ) = log p̂(Ek )p(Dk |Ek ) (39), and thus should already be available: P + log θ exp(Eq1 (X,Λ,P ) log p(Y |θ, P, X)p̄(θ, M |Λ, Dk )) sum k • Log-domain terms log Sj , along with log Sj , Xh = log p̂(Ek )p(Dk |Ek )+log exp(Eq1 (Λ) log p̄(θ, M |Λ, Dk )) d k hkj = log(pe′ k pk Sj ) θ YM i × exp Eq1 (X,P ) log p(Yj |θj , P, X) . for k = 1, ..., Kn , and h0j = log(Λ0 ℓ0 (Yj )). j=1 Kn h X YM • The maximum value mj = maxk=0 hkj and its index Ij ∈ c = log p̂(Ek )p(Dk |Ek )−Dk Eq1 (Λ) Λk +log Ij {0, 1, ..., Kn } such that mj = hj . θ j=1 XKn k sum • The stabilised terms exp(hj −mj ), and their sum Uj = e′ d exp Eq1 (Λ) log PKn s=1 ps ps Λs δ[θj = s] +Λk Dk δ[θj = k] k s̸=k exp(h − m ). j j k=0 i 1) g(Dk = 1): For Dk = 1, the two summands inside each + Λ0 δ[θj = 0] exp Eq1 (X,P ) log p(Yj |θj , P, X) , logarithm in (38) are positive, so we can apply the log-sumexp trick to obtain: c where = denotes equality up to additive constant terms with M respect to Ek and Dk . The X Q last lineQfollows from (34), (8) g(D = 1) = −η ρ + max{aj , bkj }+log 1+exp(−|aj −bkj |) k k k and the identity exp E log( j aj ) = j exp E log aj . Furtherj=1 more, since the delta indicators selectPexactly one intensity Kn sum d k aj = log Sj , bkj = log(1−pe′ (65) term for each fixed θn,j , we have f ( k=0 bk δ[θn,j = k]) = k pk )+log Sj .
6
The quantities aj and bkj are either precomputed or can be efficiently evaluated from stored quantities. The main computational cost lies in evaluating log(1 + exp(−|aj − bkj |)) for all j and k = 1, ..., Kn , which requires Kn M exponential and logarithmic operations in total to compute all g(Dk = 1). 2) g(Dk = 0): For Dk = 0 the computation is slightly more complex. From (38) we have g(Dk = 0) =
XM j=1
gjk ,
d k gjk = log(Sjsum − pe′ k pk Sj ). (66)
To evaluate gjk numerically stably for all k = 1, . . . , Kn and j = 1, . . . , M , thereby ensuring sufficient precision for all g(Dk = 0), we consider two scenarios: a) For k ̸= Ij : We use gjk = mj + log(Ujsum − exp(hkj − mj )),
(67)
where mj , Ujsum , exp(hkj −mj ) are all precomputed and stored. Since k ̸= Ij , the largest stabilised term remains in the sum after subtraction, so Ujsum − exp(hkj − mj ) ∈ [1, Kn ]. This keeps the logarithm numerically stable. b) For k = Ij : The subtraction in the expression above may be close to zero, risking underflow before the logarithm is taken. In this case, we apply the log-sum-exp trick using the second-highest value cj from {h0j , h1j , ..., hK j }: XKn s gjk =cj + log exp(h0j − cj ) + s=1 exp(hj − cj ) . (68) s̸=k
Since cj is the largest value after excluding hkj , the largest exponentiated term in the sum is 1, while all remaining terms lie in (0, 1]. This gives a stable log-sum-exp evaluation of gjk . 3) Cost analysis: We now assess the cost in terms of exponential and logarithmic operations. Computing all g(Dk = 0) for k = 1, ..., Kn requires evaluating all gjk via either (67) or (68). Specifically: 1) Each use of (67) involves one exponential and one logarithm, and is applied up to Kn M times. 2) Each use of (68) involves Kn exponentials and one logarithm, and is applied at most M times. Hence, the total number of exponentials and logarithms required is bounded by 2Kn M . Recall that computing all g(Dk = 1) for k = 1, ..., Kn via (65) mainly requires Kn M exponential and logarithmic operations in total. The numerically stable evaluation of all g(Dk ) therefore retains O(Kn M ) complexity. In practice, the total number of exponential/logarithmic evaluations is bounded by 3Kn M , making it comparable to performing two to three updates of q1 (θ) in (30). This remains substantially more efficient than the direct computation via (36), which requires at least O(Kn2 M ) operations. A PPENDIX G A DDITIONAL IMPLEMENTATION DETAILS This appendix supplements Section VII with additional implementation details on variational initialisation and movingobject selection.
1) Variational initialisation and closely spaced birth pruning: Points (i)–(iii) of Section IV-D suggest using a flat factorised birth prior, initialising newly born components over distinct small positional regions, and tiling these regions over the surveillance area. Building on this birth initialisation, we recommend first running two to three CAVI iterations for existing objects k = 1, . . . , Kn−1 with births disabled, so that existing tracks are prioritised in associating with measurements. The resulting q1 (Xn,k ) for existing objects is then combined with the birth initialisation above to start the subsequent CAVI inference with new births included. This reduces the risk that newborn objects explain measurements belonging to existing tracks. In addition, when updating q1 (θn ) in the first CAVI iteration via (30), we recommend using a clutter likelihood ℓ0 (Yn,j ) lower than the expected true level. This biases the early association update towards explaining measurements by objects, which helps q1 (Xn,k ) localise within its prescribed region. After CAVI converges, multiple variational components may represent a single physical object. This is most likely under non elliptical measurement patterns or gradual target entry, where an existing component explains the already observed portion while a birth component captures the newly observed part. Such duplication inflates the estimated object count. A simple remedy is to prune newly born components whose converged positional means lie very close to other objects, and then rerun CAVI to restore convergence. 2) Post-processing details for Doppler PiVoT moving object selection: This appendix details the post-processing stage for moving-object selection used in the Doppler PiVoT automotive radar experiments. The procedure is applied after inference, operates only on inferred track statistics, and does not modify the variational updates. Its purpose is to extract a reliable subset of moving road users from the broader set of scene elements tracked by PiVoT. Tracks are first nominated as moving candidates using minimum speed gates and velocity Mahalanobis distance over short time windows, together with basic directional and trajectory smoothness checks. This stage removes stationary tracks and artefacts induced by flickering clutter and ego motion. An example of surviving candidates are shown in Fig. 5(b). Second, candidates are examined for Doppler informativeness. Tracks whose associated Doppler velocities are largely perpendicular to the estimated motion direction, or have weak projected components along that direction, are treated as Doppler uninformative and flagged as unreliable. This suppresses objects for which stationary and moving hypotheses are unidentifiable under the adopted model. For example, the orange ellipses in Fig. 5(b) correspond to stationary road barrier segments for which the radar line of sight is nearly orthogonal to the ego motion. These segments consistently appear at the side of the ego vehicle with near-zero Doppler and can therefore be misinterpreted as moving consistently with the ego vehicle, rather than being stationary. Such tracks are rejected at this stage unless supported by smooth trajectories, kinematic consistency, and evidence of independent motion. Third, a concrete object shape constraint is enforced. Objects are rejected if (i) estimated shapes are implausible for
7
road participants (e.g. are excessively large or have unrealistic aspect ratios), or (ii) their associated detections fail to dominate the measurements within the estimated 95 percent shape ellipse (e.g. excessive stationary measurements inside). Tracks removed at this stage are shown as dashed ellipses in Fig. 5(b). Finally, a ghost suppression step removes objects arising from multipath reflections, leveraging PiVoT’s ability to track stationary structures. Potential ghost and true target pairs are first identified based on symmetric motion directions about the perpendicular bisector of their positions, with ghosts characterised by longer sensor range, shorter track duration, and fewer associated measurements. These candidates are then validated by checking whether the inferred reflection points lie within stationary objects. Confirmed ghost tracks are removed. Tracks rejected by this step are shown as cyan ellipses in Fig. 5(b), whose reflection points fall within the elongated stationary ellipses corresponding to road barriers in Fig. 5(a). Only tracks that satisfy all consistency checks are reported as moving objects (the remaining red ellipses in Fig. 5 (b)). All criteria and window lengths are fixed across experiments. This post-processing stage is computationally negligible compared to inference and serves solely to improve the reliability and interpretability of the reported outputs.
T HEORETICAL FOUNDATIONS AND PROOF OF T HEOREM IX.1 This appendix presents the proof of Theorem IX.1, supported by auxiliary lemmas analysing the CAVI updates. 1) Preliminary algebraic lemmas: We begin by presenting algebraic lemmas that will be used in the subsequent analysis.
1 This inequality always holds if − pq < − 1e , establishing 1 condition 1 of Lemma H.2. Otherwise, when − pq ≥ − 1e , 1 ) or the inequality holds if and only if t < W−1 (− pq 1 1 1 t > W0 (− pq ). Substituting back t = p − x , this becomes 1 1 1 1 1 1 x > p − W−1 (− pq ) or x < p − W0 (− pq ). Since both 1 1 W−1 (− pq ) and W0 (− pq ) are negative, rearranging yields condition 2 of Lemma H.2, completing the proof.
Corollary H.2.1. Let p, q, r, x > 0. Then, 1 x>p 1− 1 + q exp(− xr )
(72)
if and only if one of the following conditions holds: 1 1) pq pq r r r exp(− p ) < e, or 2) r exp(− p ) ≥ e and either x < p − −1 −1 r r r 1 or x > p1 − 1r W0 (− pq exp( pr )) . r W−1 (− pq exp( p )) Proof. Let x′ = xr , p′ = pr , q ′ = q exp(− pr ). Then, inequality 1 (72) is equivalent to x′ > p′ 1 − 1+q′ exp( . Applying 1 1 − ) p′ x′ Lemma H.2 completes the proof.
Lemma H.3. Suppose the prior covariance Σ′k in (12), used in the state update (27), satisfies HΣ′k H ⊤ = cIdY , where c > 0 is a scalar and IdY is the identity matrix. Then the updated 1 −1 covariance Σk in (27) satisfies HΣk H ⊤ = (R−1 . k + c IdY ) Proof. Using the Gaussian update formula and Woodbury identity, the updated covariance Σk in (27) is
(69)
Proof. Since a + 1 ≤ a + 1/x, adding α − a ≥ 0 to both the a+1 numerator and denominator of a+1/x does not decrease the a+1 α+1 ratio; hence, a+1/x ≤ α+1/x . Therefore, we have 1 a+1 1 1 α+1 ≥ 1− = 1− , 1− 1 + a α + 1/x 1 + a a + 1/x 1 + ax hence completing the proof. Lemma H.2. Let p, q, x > 0, and let W−1 : [− 1e , 0) → (−∞, −1], W0 : [− 1e , ∞) → [−1, ∞) denote the −1 and 0 branches of the Lambert W function, respectively. Then, 1 (70) x>p 1− 1 + q exp( p1 − x1 ) if and only if one of the following conditions holds: 1) pq < e, −1 1 or 2) pq ≥ e and either x < p1 − W−1 (− pq or x > 1 −1 1 . p − W0 (− pq p Proof. Let t = p1 − x1 , such that x = 1−pt . Then, (70) becomes
p 1 pqet > p(1 − ) = 1 − pt 1 + qet 1 + qet
1 + qet > qet − pqtet ⇐⇒ tet > −(pq)−1 .
2) Analysis of CAVI update: Next, we present two lemmas and a proposition concerning the employed CAVI updates and their bounding behaviour.
A PPENDIX H
Lemma H.1. Let 0 < a ≤ α, and 0 < x ≤ 1. Then, 1 α+1 1 ≤ 1− . 1− 1 + ax 1 + a α + 1/x
Since 1 − pt = 1 − (1 − p/x) > 0, multiplying both sides of (71) by (1 − pt)(1 + qet )/p yields
(71)
Σk = (Σ′k
−1
−1 + H ⊤ R−1 k H)
= Σ′k − Σ′k H ⊤ (Rk + HΣ′k H ⊤ )−1 HΣ′k . Then, applying H and H ⊤ to both sides, we obtain HΣk H ⊤ = HΣ′k H ⊤ − HΣ′k H ⊤ (Rk + HΣ′k H ⊤ )−1 HΣ′k H ⊤ = cIdY − c2 (Rk + cIdY )−1 = c(IdY − ( 1c Rk + IdY )−1 ) −1 1 −1 = c[IdY − (IdY − (cR−1 )] = (R−1 , k + IdY ) k + c IdY )
where the second last equality follows from the Woodbury identity. This completes the proof. PM Lemma H.4. Let j=1 q1 (θj = k) < s. Then the ηk , ϕk , Φk obtained from updates (28), (29) using such a q1 (θ) satisfy: ψ(ηk ) < ψ(ηk′ + s),
ψdY (ϕk /2) < ψdY ((ϕ′k + s)/2), (73)
log |Φk | ≤ log |Φ′k |,
(74)
where dY is the dimension of yj . Proof. The inequalities in (73) follow from the strict monotonicity of the digamma function ψ(η) for η > 0 and the multivariate digamma function ψdY (ϕ/2) for ϕ > dY − 1. The inequalities in (74) follow from the fact that Φ′k − Φk is a positive semi definite matrix, as implied by (29).
8
Proposition H.5. Define c(s) as in (48). Let Assumption 1 hold, and suppose CAVI applies one round of updates sequentially via (29)-(27) with q1 (θ) = q1− (θ), followed by updating q1 (θ) to q1+ (θ) using (31)-(30). Then for any newlyborn k = Kn−1 + 1, ..., Kn and any s > 0 such that PM object − + j=1 q1 (θj = k) < s, q1 (θj = k) satisfies q1+ (θj = k) < 1 −
1 dY /2 P 1+ dkc(s) (yj ) exp(− M q − (θ =k) ) j=1
,
(75)
j
1
where dk (yj ) is evaluated from (52) using the parameters ϕk , Φk , µk obtained from that round of updates (29) and (27). Proof. From (30) and Assumption 1, we have q1+ (θj = k) =
d k d k pe′ pe′ k pk Sj k pk Sj PKn e′ d l ≤ e′ λc + pk pdk Sjk λc + l=1 pl pl Sj 1 =1− (76) 1 e′ d k 1 + λc pk pk Sj
where Sjk > 0 is defined in (30). Using (31), (48), and Lemma d k H.4, the term λ1c pe′ k pk Sj can be upper bounded as 1 e′ d k c(s) exp(−0.5 Tr(ϕk Φk HΣk H ⊤)) pk pk Sj < λc exp(0.5(yj −Hµk )⊤ϕk Φk (yj −Hµk )) c(s) = exp(−0.5 Tr(ϕk Φk HΣk H ⊤)), (77) dk (yj ) where dk (yj ) is defined in (52), and µk , Σk are obtained from state update (27). From Lemma H.3 and Assumption 1, ⊤ HΣk H ⊤ equals PM the− Rk in (27), and thus Tr(ϕk Φk HΣk H ) equals dY / j=1 q1 (θj = k). Subsequently, (77) becomes 1 e′ d k c(s) dY /2 . pk pk Sj < exp − PM − λc dk (yj ) j=1 q1 (θj = k) Substituting this inequality into (76) for a further upper bound yields (75), completing the proof. 3) Proof of Theorem IX.1: With Proposition H.5, Corollary H.2.1, and Lemma H.1, we are now ready to prove Theorem IX.1. PM Proof. Suppose j=1 q1− (θj = k) < s. Proposition H.5 first gives the one-step upper bound (75) on q1+ (θj = k). We further 1 , bound its right-hand side using Lemma H.1 with x = dk (y j) a = c(s) exp(− PM dqY−/2 ), and α = c(s) exp(− d2sY ) = (θ =k) j=1
1
j
t(s), where t(s) is as defined in (49). This gives 1 q1+ (θj = k) < 1 − 1 + c(s) exp(− PM dqY−/2 ) (θ =k)
where Bk (s) is an admissible bound from Definition 1. We now seek a sufficient condition under which X M 1 Bk (s) 1− < q1− (θj = k), 1 + c(s) exp(− PM dqY−/2 ) j=1 (θ =k) j=1
1
j
so thatP applying (79) yields the desired bound (55) in Theorem M ∆k (s) PM − + IX.1: q (θ = k) < j j=1 q1 (θj = k). This j=1 1 Bk (s) sufficient condition can be established via Corollary PM H.2.1, by setting p = Bk (s), q = c(s), r = dY /2, x = j=1 q1− (θj = k). Corollary H.2.1 then implies that the inequality holds if either of the following two conditions PMis satisfied: dY /2 − k (s) exp(− ) ≥ e and 1) c(s)B j=1 q1 (θj = k) < Vk (s), dY /2 Bk (s) where Vk (s) is defined in (54), or /2 k (s) exp(− BdYk (s) ) < e. 2) c(s)B dY /2 PM Recall that we have also assumed j=1 q1− (θj = k) < s. Therefore, combining the above gives the threshold PM conditions − Lk (s) in (53). Hence, if q (θ = k) < Lk (s), then j j=1 1 PM + ∆k (s) PM − j=1 q1 (θj = k) < Bk (s) j=1 q1 (θj = k). Since ∆k (s) ≤ B (s) by the admissibility in (50), PMk + PM defined − it follows that q (θ = k) < q (θ = k) < j j j=1 1 j=1 1 Lk (s), so the sufficient condition remains satisfied at the next iteration. By induction, the same bound continues to hold, and PM q (θ j=1 1 j = k) decreases monotonically across iterations, k (s) with the ratio of successive values bounded by ∆ Bk (s) ≤ 1. PM It remains to show that j=1 q1 (θj = k) converges to zero. PM (r) Let βr := j=1 q1 (θj = k) denote the association mass after the r-th subsequent update from any iteration satisfying βr < Lk (s), and define 1 . f (x) := Bk (s) 1 − 1 + c(s) exp(−dY /(2x)) The preceding analysis shows that βr is monotonically decreasing and bounded below by zero, and therefore βr → β∞ ≥ 0. Suppose β∞ > 0. Then since β∞ < Lk (s), the construction of Lk (s) gives f (β∞ ) < β∞ . However, (79) and ∆k (s) ≤ Bk (s) give βr+1 < f (βr ). Taking limits in this inequality, using βr+1 → β∞ and continuity of f , yields β∞ ≤ f (β∞ ), contradicting f (β∞ ) < β∞ . Hence, β∞ = 0, completing the proof.
A PPENDIX I A DDITIONAL D ETAILS ON THE T HRESHOLD Lk (s) IN T HEOREM IX.1
This appendix provides additional details on the threshold in Theorem IX.1. We first analyse how the threshold Lk (s) 1 + t(s) depends on the admissible bound Bk (s) and the parameter s, × . (78) showing that a tighter Bk (s) gives a no smaller threshold for dk (yj ) + t(s) fixed s. We then give a practical construction of such a bound Summing over j and using ∆k (s) in (51) yields for the known-shape case. M X 1) Threshold Lk (s) behaviour with respect to s and Bk (s): 1 q1+ (θj = k) < 1 − ∆k (s) Here we examine how the threshold Lk (s) in Theorem IX.1 dY /2 1 + c(s) exp(− PM q− (θ =k) ) j=1 j changes with s and the admissible bound Bk (s). The main j=1 1 message is twofold. First, for fixed s, a tighter/lower ad1 ∆k (s) = Bk (s) 1− , (79) missible bound B (s) gives a no smaller threshold L (s). dY /2 k k B (s) k 1 + c(s) exp(− PM q− (θ =k) ) j j=1 1 Second, the dependence on s is more subtle: decreasing s j=1
1
j
9
improves the internal bound terms, in particular increasing Vk (s) and making the first-case condition in (53) easier to satisfy. However, Lk (s) is also capped by s, so Lk (s) is not monotone in s in general. We first consider the dependence on Bk (s) for fixed /2 s. Since Bk (s) exp(− BdYk (s) ) decreases with smaller Bk (s), decreasing Bk (s) makes the first case in (53), namely c(s)Bk (s) /2 exp(− BdYk (s) ) < e, easier to satisfy. Moreover, in dY /2 the second case of (53), it will be shown below that Vk (s) increases as Bk (s) decreases. This gives the claimed monotonicity of Lk (s) with respect to Bk (s). Indeed, fix s and compare two admissible bounds B1 < B2 . If B2 satisfies the first case in (53), then so does B1 , and both thresholds equal s. If B2 is in the second case but B1 moves to the first case, then the new threshold is s, which is no smaller than min{s, Vk (s)}. Finally, if both remain in the second case, then Vk (s) is larger for B1 , and hence min{s, Vk (s)} is no smaller. Therefore, for fixed s, tightening/lowering Bk (s) can only increase, or leave unchanged, the sufficient threshold Lk (s). We next consider the role of s. Decreasing s > 0 decreases c(s) in (48) and hence also decreases t(s) in (49). Since dk (yj ) ≥ 1, each term in (51) is nondecreasing in t(s), so ∆k (s) also decreases. Therefore a smaller admissible bound Bk (s) may be obtained. In addition, the condition term /2 c(s)Bk (s) exp(− BdYk (s) ) in (53) decreases with a smaller s, so dY /2 the first case is easier to satisfy. The analysis below also shows that Vk (s) increases as s decreases. However, since Lk (s) is either s or min{s, Vk (s)}, these favourable changes do not imply monotonicity of Lk (s) with respect to s. Finally, we show that Vk (s) increases when either s or Bk (s) decreases. Recall from Section H-3 that Vk (s) is constructed using Corollary H.2.1. Specifically, define the function −1 1 1 r r fl (p, q, r) = − W−1 − exp . (80) p r pq p Then V (s) = fl (p, q, r) with p = Bk (s), q = c(s), r = dY /2. Smaller Bk (s) implies smaller p, and smaller s implies smaller c(s) and hence smaller q. Since fl (p, q, r) increases as p or q decreases, as formalised in Lemma I.1 presented below, this establishes the desired monotonicity: Vk (s) increases as either s or Bk (s) decreases. Lemma I.1. Let p, q, r > 0. The function fl (p, q, r) defined r in (80), when well-defined (i.e., pq r exp(− p ) ≥ e), increases as either p or q decreases, or as r increases. 1 Proof. Define gx (p, q, r) = p 1 − 1+q exp(− . When r x) fl (p, q, r) is well defined, Corollary H.2.1 gives that the solution set of the inequality x > gx (p, q, r) is −1 S r x ∈ (−∞, fl (p, q, r)) ( p1 − 1r W0 (− pq exp( pr )) , ∞). Observe that gx (p, q, r) decreases when either p or q decreases, or when r increases. Hence, any x satisfying x > gx (p, q, r) for a given (p, q, r) still satisfies the inequality after any such parameter change. Therefore, the original solution set is contained in the updated solution set. Since the two intervals above are separated by a nonempty gap except at the boundary case, the lower-branch interval (−∞, fl (p, q, r)) cannot shrink
under such a parameter change. Consequently, fl (p, q, r) increases as either p or q decreases, or as r increases. 2) Construction of Bk (s) for the known-shape case: We consider the case where the object shape is known or the shape prior is sufficiently strong that the shape update is negligible. In this case, we assume that the expected precision of q1 (Pk ) remains fixed at its prior value, namely ϕk Φk = ϕ′k Φ′k . We now construct a tighter admissible bound Bk (s) in (50) for this setting. From (51), we have ∆k (s) =
M X
1 + t(s) , ⊤ ′ ′ exp(0.5(y −Hµ ) j k ϕk Φk (yj −Hµk )) + t(s) j=1
where t(s) is defined in (49). By Remark 1 under Assumption 1, Hµk = y k , where y k in (27) is a weighted average of the measurements. Hence Hµk lies in the surveillance region S containing all measurements, given that S is convex, e.g. rectangular. Therefore, ∆k (s) ≤ max y∈S
M X
1 + t(s) . ⊤ϕ′ Φ′ (y −y)) + t(s) exp(0.5(y −y) j k k j j=1
This maximum can be bounded by a grid based construction. Let G1 , . . . , GNG be a partition of S. Then a practical and relatively tight admissible bound Bk (s) can be obtained as ∆k (s) ≤ max
max
i=1,...,NG y∈Gi
M X
1 + t(s) ⊤ϕ′ Φ′ (y −y)) + t(s) exp(0.5(y −y) j k k j j=1
M X
1 + t(s) i=1,...,NG miny∈Gi exp(0.5(yj −y)⊤ϕ′k Φ′k (yj −y)) + t(s) j=1
≤ max
:= Bk (s), where the last inequality becomes exact in the limiting case where the grid is arbitrarily fine, i.e. each Gi consists of individual points in S. The inner minimisation miny∈Gi exp(0.5(yj − y)⊤ ϕ′k Φ′k (yj − y)) for each j and Gi can be efficiently computed for diagonal Φ′k and rectangular grid cells. For non diagonal Φ′k , one can instead build the grid in transformed coordinates where the quadratic form becomes diagonal. This gives the same type of computation and remains cheap for dY = 2 or 3. A PPENDIX J PARAMETERISATION OF COMPARED METHODS IN S ECTION VIII-A Here we list the parameterisations of PMBM, SPA, and PiVoT used in the simulation experiments in Section VIII-A. For all methods, the kinematic transition model uses the ground truth parameterisation in (47), and the positional birth prior covers the entire surveillance region. For PMBM, the parameterisation follows [12]. PMBM reports objects with existence probability above 0.5 under the highest weight global hypothesis. The gate probability is set to 0.999, the maximum number of global hypotheses is set to 100, and Murty’s algorithm is capped at 20 assignments. The Poisson birth rate is set to 0.5, and survival probability is set to 0.99. The clutter rate and detection probability are given
10
as the ground truth for all DS1–DS6. The rate prior follows a Gamma distribution with shape parameter 2 and mean set to 8 for DS3–DS6, giving a prior that accommodates both measurement rates well. For DS1–DS2, its mean is instead set to the ground truth. The shape prior follows an inverse Wishart distribution with degrees of freedom 10. Its mean is set to the ground truth 800I2 for DS1–DS2, and to 400I2 for DS3– DS6 to better accommodate both object types. The forgetting factors for shape and rate prediction are set to τ = 10 and η = 1.2, respectively, see Table III of [12]. These factors are defined differently from those in PiVoT, but the chosen values inflate the uncertainty to a comparable level. Finally, PMBMC and PMBMF use 16 and 100 DBSCAN distance values, respectively, evenly spaced over [1, 100], a range found to give strong overall performance. For SPA, the object declaration threshold is set to 0.5, the pruning threshold to 10−3 , and the number of iterations to 3. The mean number of births is set to 0.5, the survival probability to 0.99, and the clutter rate to the ground truth. Since SPA does not handle rate estimation or detection probability, we first take the measurement rate as the ground truth rate for DS1–DS2 and the larger rate 8 for DS3–DS6, then multiply it by the ground truth detection probability (not equal to 1 only for DS2 and DS3) to account for missed detections. The shape prior is set as an inverse Wishart distribution with degrees of freedom 10. Its mean is set to the ground truth 800I2 for DS1– DS2, and to 400I2 for DS3–DS6 to better accommodate both object types. The degrees of freedom for shape prediction is set to 100. SPA is implemented with different particle sizes, indicated by the subscripts in Table I. For PiVoT, identical uninformative birth priors are used as suggested in Section IV-D, with birth probability pbn,k = 0.2 and ineffective birth removal threshold L = 0.5. The thresholds Ppru , Prep , and Pstp in Section VII-B are set to 0.03, 0.8, and 0.3, respectively. To assess PiVoT’s robustness across different scenes, we assume an unknown clutter rate and use the same uninformative shape and birth priors for all DS1– DS6. The shape (precision matrix) prior is set as a Wishart distribution with mean I2 /400 and degree of freedom 3, and the rate prior is set as a Gamma distribution with shape parameter 2 and mean 8. The unknown clutter rate is set to the current number of measurements at each time step, which is a crude upper bound but works well under the current parameterisation. The transition parameters psn,k , γΛ,n , and γP,n in (46) are all set to 0.9. The initial q1 is shown in Fig. 3(a) in Section IV-D.