Newton-Flow Particle Filters based on Generalized Cram\'er Distance

TL;DR

This paper introduces a novel recursive particle filter that uses a generalized Cramér distance and Newton flow to efficiently update particles in high-dimensional spaces without degeneracy, improving stability and simplicity.

Contribution

The paper presents a new particle filtering method based on a differentiable generalized Cramér distance and Newton flow, avoiding density estimation and ensuring non-degeneracy in high dimensions.

Findings

The filter is simple to implement and efficient.

It never degenerates in high-dimensional problems.

It can replace classic particle filters as a plugin.

Abstract

We propose a recursive particle filter for high-dimensional problems that inherently never degenerates. The state estimate is represented by deterministic low-discrepancy particle sets. We focus on the measurement update step, where a likelihood function is used for representing the measurement and its uncertainty. This likelihood is progressively introduced into the filtering procedure by homotopy continuation over an artificial time. A generalized Cram\'er distance between particle sets is derived in closed form that is differentiable and invariant to particle order. A Newton flow then continually minimizes this distance over artificial time and thus smoothly moves particles from prior to posterior density. The new filter is surprisingly simple to implement and very efficient. It just requires a prior particle set and a likelihood function, never estimates densities from samples, and can be used as a plugin replacement for classic approaches.

Equations

$${\underline{\mathbf{x}}}_{k+1}={\underline{a}}_k({\underline{\mathbf{x}}}_k,{\underline{u}}_k,{\underline{\mathbf{w}}}_k)$$

$${\underline{\mathbf{y}}}_k={\underline{h}}_k({\underline{\mathbf{x}}}_k)+{\underline{\mathbf{v}}}_k$$

$$f_k^L(\underline{x})=f_k({\hat{\underline{y}}}_k|{\underline{x}}_k).$$

$$f_e(\underline{x})=\frac{1}{c_e}\cdot f_p(\underline{x})\cdot f_L(\underline{x})$$

$$c_e={\int }_{{\mathrm{I}\mathrm{R}}^N}{f}_p(\underline{x})\cdot f_L(\underline{x})\mathrm{d}\underline{x}$$

$$f_p(\underline{x})=\sum _{i=1}^L{w}_{p,i}\cdot \delta (\underline{x}-{\underline{x}}_{p,i}),$$

$${\tilde{f}}_e(\underline{x})\propto f_L(\underline{x})\cdot f_p(\underline{x}),$$

$${\tilde{f}}_e(\underline{x})$$

$$=\sum _{i=1}^L\underset{{\bar{w}}_{e,i}}{\underbrace{w_{p,i}\cdot f_L({\underline{x}}_{p,i})}}\cdot \delta (\underline{x}-{\underline{x}}_{p,i}).$$

$${\tilde{w}}_{e,i}=\frac{{\bar{w}}_{e,i}}{{\sum }_{i=1}^L{\bar{w}}_{e,i}}.$$

$${\tilde{f}}_e(\underline{x})=\sum _{i=1}^L{\tilde{w}}_{e,i}\cdot \delta (\underline{x}-{\tilde{\underline{x}}}_{e,i}).$$

$$f_L(\underline{x},\gamma )=\{\begin{array}{ll}1& \gamma =0\\ \text{convenient interpolation}& \gamma \in (0,1)\\ f_L(\underline{x})& \gamma =1\end{array}.$$

$${\tilde{f}}_e(\underline{x},\gamma )\propto f_L(\underline{x},\gamma )f_p(\underline{x}).$$

$$f_L(\underline{x},\gamma )={\left[f_L(\underline{x})\right]}^{\gamma }.$$

$$f_L(\underline{x},\gamma )={\left[f_L(\underline{x})\right]}^{g(\gamma )}.$$

$${\dot{f}}_L(\underline{x},\gamma )$$

$$=\frac{\partial }{\partial \gamma }\left\{\exp (g(\gamma )\log (f_L(\underline{x})\left)\right)\right\}$$

$$=\log (f_L(\underline{x}))\dot{g}(\gamma )f_L(\underline{x},\gamma ).$$

$$f_L(\underline{x},\gamma )=\exp \left(-\gamma \frac{1}{2}(\underline{y}-\underline{h}(\underline{x}){)}^{\top }{\mathbf{C}}_v^{-1}(\underline{y}-\underline{h}(\underline{x}))\right).$$

$$D(\underline{\eta }(\gamma ),\gamma )=D\left({\tilde{f}}_e(\underline{x},\gamma ),f_e\left(\underline{x},\underline{\eta }(\gamma )\right)\right).$$

$$\underline{G}(\underline{\eta }(\gamma ),\gamma )=\underset{D_{\underline{\eta }}(\underline{\eta }(\gamma ),\gamma )}{\underbrace{\frac{\partial D(\underline{\eta }(\gamma ),\gamma )}{\partial \underline{\eta }(\gamma )}}}=\underline{0}.$$

$$\underset{{\underline{G}}_{\underline{\eta }}(\underline{\eta }(\gamma ),\gamma )}{\underbrace{\frac{\partial \underline{G}(\underline{\eta }(\gamma ),\gamma )}{\partial {\underline{\eta }}^T(\gamma )}}}\cdot \frac{\partial \underline{\eta }(\gamma )}{\partial \gamma }+\underset{{\underline{G}}_{\gamma }(\underline{\eta }(\gamma ),\gamma )}{\underbrace{\frac{\partial \underline{G}(\underline{\eta }(\gamma ),\gamma )}{\partial \gamma }}}=\underline{0}$$

$$\underset{D_{\underline{\eta }\underline{\eta }}(\underline{\eta }(\gamma ),\gamma )}{\underbrace{\frac{{\partial }^{2}D(\underline{\eta }(\gamma ),\gamma )}{\partial \underline{\eta }(\gamma )\partial {\underline{\eta }}^T(\gamma )}}}\cdot \underset{\dot{\underline{\eta }}(\gamma )}{\underbrace{\frac{\partial \underline{\eta }(\gamma )}{\partial \gamma }}}+\underset{D_{\underline{\eta }\gamma }(\underline{\eta }(\gamma ),\gamma )}{\underbrace{\frac{{\partial }^{2}D(\underline{\eta }(\gamma ),\gamma )}{\partial \underline{\eta }(\gamma )\partial \gamma }}}=\underline{0}.$$

$$\underline{J}(\underline{\eta }(\gamma ),\gamma )=D_{\underline{\eta }\gamma }(\underline{\eta }(\gamma ),\gamma )={\underline{G}}_{\gamma }(\underline{\eta }(\gamma ),\gamma )$$

$$\mathbf{H}(\underline{\eta }(\gamma ),\gamma )=D_{\underline{\eta }\underline{\eta }}(\underline{\eta }(\gamma ),\gamma )={\underline{G}}_{\underline{\eta }}(\underline{\eta }(\gamma ),\gamma ).$$

$$\mathbf{H}(\underline{\eta }(\gamma ),\gamma )\cdot \dot{\underline{\eta }}(\gamma )+\underline{J}(\underline{\eta }(\gamma ),\gamma )=\underline{0}.$$

$$\dot{\underline{\eta }}(\gamma )=-{\mathbf{H}}^{-1}(\underline{\eta }(\gamma ),\gamma )\cdot \underline{J}(\underline{\eta }(\gamma ),\gamma )\stackrel{\mathrm{def}}{=\text{joinrel}=}\underline{R}(\underline{\eta }(\gamma ),\gamma ).$$

$$\dot{\underline{\eta }}(\gamma )=\underline{R}(\underline{\eta }(\gamma ),{\tilde{f}}_e(\underline{x},\gamma ),{\dot{\tilde{f}}}_e(\underline{x},\gamma )).$$

$$f_L^{\text{eff}}(\underline{x},\gamma )=\frac{f_L(\underline{x},\gamma )}{f_L(\underline{x},{\gamma }_{\ast })}.$$

$$\dot{\underline{\eta }}(\gamma )=\underline{R}(\underline{\eta }(\gamma ),f_L^{\text{eff}}(\underline{x},\gamma ),{\dot{f}}_L^{\text{eff}}(\underline{x},\gamma )).$$

$${\tilde{f}}_e(\underline{x},\tilde{\underline{w}}(\gamma ),\tilde{\underline{X}})=\sum _{i=1}^L{\tilde{w}}_i(\gamma )\cdot \delta (\underline{x}-{\tilde{\underline{x}}}_i)$$

$$\tilde{\underline{X}}={\left[\begin{array}{c}{\tilde{\underline{x}}}_1^T,{\tilde{\underline{x}}}_2^T,\dots ,{\tilde{\underline{x}}}_L^T\end{array}\right]}^T$$

Full text

Note

The first version of this manuscript was prepared as a basis for discussions during my plenary talk at the 2025 IEEE International Conference on Multisensor Fusion and Integration (MFI 2025) at Texas A&M University, College Station, Texas on September 3, 2025.

This slightly modified version was prepared for my plenary talk at the 17th Symposium Sensor Data Fusion in Bonn, Germany, on November 25, 2025. It includes changes based on comments by Ryne Beeson, Ondřej Straka, Daniel Frisch, and Yiqing Wang. In addition, the distance measure was modified to naturally avoid singularities.

The next version will include sections on implementation details, state-of-the-art, and an evaluation of the proposed filtering method.

Acronyms

DMD Dirac mixture density GMD Gaussian mixture density NFPF Newton-flow particle filter ODE ordinary differential equation pdf probability density function SDE stochastic differential equation

Notation

Vectors will be indicated by small underlined letters, $\displaystyle\underline{0}$ and $\displaystyle\underline{1}$ are a vector of zeros and a vector of ones, respectively. Random values will be denoted by small boldface letters, so random vectors are small boldface underlined letters. Matrices are denoted by capital boldface letters. The identity matrix of dimension $\displaystyle N$ is $\displaystyle{\mathbf{I}}_{N}$.

I Introduction

I-A Estimation Problem to be Solved

We consider recursive estimation of the state of a nonlinear dynamic system based on measurements sequentially arriving at discrete time steps. The evolution of the system state is described by a process model either in continuous time by a stochastic differential equation (SDE) or in discrete time by a stochastic difference equation, both usually in state-space form. The relation of states to measurements is described by a stochastic measurement model. The goal is to continually provide an estimate of the true system state with each new measurement.

Due to the uncertainties in both process and measurement model, only an uncertain estimate of the state can be provided, which is described by a probability density function (pdf). This pdf is propagated through the process model (prediction step or time update) and updated with a new measurement (filter step or measurement update).

For performing time and measurement update, specific density representations have to be selected. As performing the prediction step for significant nonlinearities is demanding for continuous densities, particle representations are typically selected. These simplify the prediction step, but lead to difficulties in performing the filter step. The focus of this paper is a new method to execute the filter step.

I-B Challenges of Particle-based Estimation

The prediction step entails propagating the particles forward through the system model, which is relatively straightforward. The filtering step, however, although looking almost trivial at first glance, is far more complex to actually implement. It requires finding a particle set compatible with both the prediction and a given measurement. In a Bayesian setting, this results in a multiplication of the prior density and a likelihood function (resulting from measurement equation and actual measurement). Naïvely performing the Bayesian update step would result in a set of weighted particles for the posterior density. Especially for “narrow” likelihoods, many particles would receive zero weights, and not contribute to the density representation. This effect is called particle degeneration, and many methods for its avoidance have been published.

Many of the methods for restoring the health of a particle set propose resampling methods that rely on the subsequent prediction step. The bootstrap particle filter, for example, just replicates particles according to their relative weight in the filter step without changing their locations. Only the process noise of the following prediction step is then responsible for introducing particle variability. This leads to an undesired dependency of the filter step from the parameters of the process model, e.g., its noise covariance, and is also problematic for chaining several filter steps without having intermediate prediction steps.

One of our goals for the filter developed in this paper is the strict separation of the filter step from the prediction step: The filter step should be self-contained.

I-C New Filtering Method: Key Ideas

In this paper, we advocate a radically different method for particle-based filtering that inherently avoids sample degeneration issues.

It is based on the idea of progressive Bayesian estimation in [1, 2, 3], which employs continuous densities for representing the posterior. Progressive estimation for Dirac mixture densities (DMDs) was introduced in [4, 5], where direct discretization of the artificial time related to the homotopy was used, resulting in a discrete flow. In contrast, in this paper, ordinary differential equations (ODEs) for continuous particle flows will be derived. Compared to [6] the set distance measure from [7] is used.

The method consists of three main ingredients that will be explained in detail in the following: (i) progressively introducing the likelihood by homotopy continuation, (ii) deterministic sampling as density representation, and (iii) using a set distance measure for deriving a Newton particle flow from prior to posterior density.

Homotopy Continuation

Instead of performing the filter step at once with the given likelihood function, an artificial time $\displaystyle\gamma$ ranging between $\displaystyle 0$ and $\displaystyle 1$ is introduced, which is then used for parametrizing the likelihood. For $\displaystyle\gamma=0$, no update is performed. For $\displaystyle\gamma=1$, the original likelihood is recovered. The posterior density is also parametrized with $\displaystyle\gamma$ so that it starts as the prior density for $\displaystyle\gamma=0$ and moves towards the (unknown) posterior for $\displaystyle\gamma\rightarrow 1$. For a particle representation, this so-called homotopy continuation approach results in a particle flow over the artificial time $\displaystyle\gamma$.

Deterministic Sampling

For the density representation by particles, we use deterministic samples (also called low-discrepancy samples). Compared to random samples, (i) the density is covered more homogeneously, (ii) the number of particles required for coverage is reduced, (iii) statistical properties, e.g., moments, show a much better convergence, (iv) and reproducibility is ensured. In this paper, we use the sampling method from [8] also used in [9].

Newton Flow Induced by Set Distance Measure

The flow induced by the homotopy continuation very rarely produces exact results. In general, we have to be content with approximations. This means that we have to continually track the closest approximation to the true posterior over the artificial time. Closeness is measured by a set distance measure, the so called generalized Cramér-von Mises distance [7, 10] or Cramér distance for short, that allows the efficient comparison of DMDs and is differentiable w.r.t. to the particle locations and weights. The optimal particle approximation is tracked over artificial time $\displaystyle\gamma$ by a Newton flow induced by this generalized Cramér distance.

I-D New Filtering Method: Unique Features

Here, we will give a short overview of the advantages of the new filtering method and also of its differences w.r.t. the state-of-the-art.

Modest Assumptions

For performing the filter step, the new method only requires (i) the prior density in the form of a (possibly weighted) particle set and (ii) a likelihood function. Nothing else is required! In particular, the underlying continuous prior density is not required and no attempt towards its reconstruction is ever made. This is the natural setup encountered in recursive particle filtering, where the last estimate is a particle set used in the next filter step. Avoiding density reconstruction is especially important in higher-dimensional spaces as density estimation from samples is notoriously difficult here.

Mapping-free Flow

The new filtering method is different from other particle flow filters in that it does not calculate an explicit mapping for transporting particles. This can be interpreted as a direct online optimization of the distance measure versus finding a mapping policy first. Such a policy is often found by interpreting the homotopic version of the measurement update equation as a Fokker-Planck equation and using the associated SDE for moving the particles. This necessitates first finding the SDE and second numerically solving it for calculating the desired particle locations, which typically leads to numerical issues due to various reasons, e.g., stiffness. In contrast, the new method directly yields a closed-form ODE by calculating the gradient and Hessian of the generalized Cramér distance. This ODE can be solved with standard methods.

Generality

Methods relying on explicit mappings between the involved random variables are ultimately restricted to calculating posterior particle sets. The reason is that calculating the flow for Gaussian mixture densities (GMDs) or exponential densities based on the SDE for the posterior random variable (associated to the Fokker-Planck equation) would add another layer of complication. In fact, the SDE is a dynamic system describing the evolution of the posterior random variable and thus, not even weight changes of a particle set are possible (as the system mapping only changes sample locations, but not sample weights). In contrast, the new filtering method can be applied with different density representations. See [11] for an early approach involving GMDs.

II Problem Formulation

II-A Models

We consider a dynamic system with random state vector $\displaystyle{\bm{\underline{x}}}_{k}\in\mathds{R}^{N}$, defined at discrete time steps $\displaystyle k\in\mathds{N}\cup\{0\}$. The state evolves according to a nonlinear system equation

[TABLE]

with nonlinear (deterministic) function $\displaystyle\underline{a}_{k}(\cdot,\cdot,\cdot)$, known input $\displaystyle\underline{u}_{k}$, and system noise $\displaystyle{\bm{\underline{w}}}_{k}$.

The state is usually hidden, i.e., $\displaystyle{\bm{\underline{x}}}_{k}$ is a latent variable not directly available, and is indirectly observed via measurements related to the state via the measurement equation

[TABLE]

with nonlinear (deterministic) function $\displaystyle\underline{h}_{k}(\cdot)$ and measurement noise $\displaystyle\underline{v}_{k}$.

The models (1) and (2) relate random variables and their realizations, and we call them generative models. Corresponding probabilistic models can be derived that describe the system evolution and the measurement extraction via conditional densities $\displaystyle f_{k}(\underline{x}_{k+1}|\underline{x}_{k})$ and $\displaystyle f_{k}(\underline{y}_{k}|\underline{x}_{k})$. For a specific measurement $\displaystyle\hat{\underline{y}}_{k}$, the so-called likelihood function $\displaystyle f_{k}^{L}$ is obtained by

[TABLE]

II-B Representation of Uncertain Latent States

In the case of stochastic model and sensor uncertainties, the latent state is represented by a probability density function (pdf). Continuous pdfs such as Gaussian or Gaussian mixture densities can be used in the case of linear systems. For nonlinear systems, however, the required transformations are difficult to handle for continuous densities, unless some form of linearization is employed, which limits the achievable quality. Hence, often sample/particle representations are used. In this paper, we employ deterministic sample representations we call Dirac mixtures, where the weights and locations are obtained by optimization rather than by random sampling from a continuous density.

II-C Estimation and Challenges

We assume that the measurements $\displaystyle\hat{\underline{y}}_{1}$, $\displaystyle\hat{\underline{y}}_{2}$, $\displaystyle\cdots$ arrive sequentially and that we desire to provide a state estimate with each new measurement. This requires the alternating execution of a prediction and a filter step.

The prediction step, or time update, entails propagating the estimated state density $\displaystyle f_{k}^{e}(\underline{x}_{k})$ at time step $\displaystyle k$ to the next time step $\displaystyle k+1$ by using the generative description in equation (1) or the probabilistic description $\displaystyle f_{k}(\underline{x}_{k+1}|\underline{x}_{k})$. This results in the predicted state density $\displaystyle f_{k+1}^{p}(\underline{x}_{k+1})$.

The filter step, or measurement update, uses a measurement taken at time $\displaystyle k$ to correct the given prediction $\displaystyle f_{k}^{p}(\underline{x}_{k})$ either based on the generative measurement equation in (2) or the likelihood function (3).

Remark II.1**.**

Please note that in this paper, we focus on a single filter step, omit the time index $\displaystyle k$, and change $\displaystyle f_{k}^{p}(\cdot)$, $\displaystyle f_{k}^{e}(\cdot)$, $\displaystyle f_{k}^{L}(\cdot)$ to $\displaystyle f_{p}(\cdot)$, $\displaystyle f_{e}(\cdot)$, $\displaystyle f_{L}(\cdot)$.

II-D Filter Step for Continuous Densities

For a continuous prior density $\displaystyle f_{p}(\underline{x})$ and a continuous likelihood function $\displaystyle f_{L}(\underline{x})$, the exact Bayesian filter step for calculating the desired continuous posterior density $\displaystyle f_{e}(\underline{x})$ is given by

[TABLE]

with normalization constant

[TABLE]

so that $\displaystyle\textstyle\int_{{\rm I\hskip-1.35526ptR}^{N}}f_{e}(\underline{x})\,\mathrm{d}\underline{x}=1$ holds. Calculation of the normalization constant is in general a major computational burden.

II-E Challenges for General Density Representations

The filter step in (4) is exact and assumes that the posterior density can completely capture the shape of the true distribution. This is true, e.g., in linear systems. For nonlinear systems, it only works when the posterior $\displaystyle f_{e}(\underline{x})$ lives in a general infinite-dimensional distribution function space. For densities with limited representation capabilities, i.e., finite mixtures or finite particle sets, as required in practical implementations, the density type changes and/or the complexity goes up. For example, a GMD becomes a non-Gaussian mixture during the update or an equally weighted particle set becomes an unequally weighted particle set.

In recursive Bayesian estimation, we desire the densities describing the considered latent variables to be closed under the measurement update step, at least in an approximate sense. This means that for a prior density of a certain type, we want the posterior density to be of the same type111Prior and posterior are then called (approximately) conjugate distributions w.r.t. a likelihood.. In that case, we can use the resulting posterior as a prior for the next recursion step without changing its mechanism. This means that we have to find the best approximation of the resulting posterior lying within the space of permissible densities. This is equivalent to projecting the true posterior to the closest permissible approximation, which requires a distance measure.

The optimization problem of finding the best approximative posterior density among the permissible densities is typically a complicated numerical minimization task that involves local minima. It is exacerbated by the fact that the prior is not a good initial guess, but usually the only available one.

II-F Filter Step for DMDs

In this paper, we focus on the case that the prior density $\displaystyle f_{p}(\underline{x})$ is not a continuous density, but given as a set of $\displaystyle L$ weighted samples (or particles). It can be formally written as a DMD

[TABLE]

with weights $\displaystyle w_{p,i}>0$, $\displaystyle\textstyle\sum_{i=1}^{L}w_{p,i}=1$ and sample locations $\displaystyle\underline{x}_{p,i}$.

For prior DMDs, the posterior is also a DMD, and we can omit the normalization constant from the Bayesian filter step (4)

[TABLE]

as normalizing $\displaystyle\tilde{f}_{e}(\underline{x})$ becomes trivial. The tilde $\displaystyle\tilde{\phantom{\ast}}$ marks this as the true posterior that will be approximated later. Plugging in $\displaystyle f_{p}(\underline{x})$ from (6) gives

[TABLE]

Normalization to obtain the posterior weights is simply performed by

[TABLE]

The posterior DMD is now given as

[TABLE]

II-G Challenges for DMDs

The posterior $\displaystyle\tilde{f}_{e}(\underline{x})$ in (10) is derived from the straightforward application of the Bayesian filter step to Dirac mixtures. Please note that in this case, only the weights change. The posterior sample locations do not change w.r.t. the prior samples, i.e., $\displaystyle\tilde{\underline{x}}_{e,i}=\underline{x}_{p,i}$ for $\displaystyle i=1,\ldots,L$. This leads to a serious problem: The samples are not equally weighted anymore and do not equally contribute to the representation of the posterior. Often, some particle weights are (close to) zero, in fact dying out, leading to particle degeneracy. A typical scenario is large system noise, which spreads the particles during the prediction step combined with low measurement noise leading to “narrow” likelihoods.

Many solutions, some systematic, many of heuristic nature, have been proposed to solve the degeneracy problem, which is a fundamental and difficult problem.

Our goal is to derive a Bayesian filter that inherently avoids degeneracy without any heuristic approaches. It should be easy to understand, simple to implement, numerically stable, and robust.

III Homotopy Continuation

In order to circumvent the problems coming with the one-shot measurement update shown in the previous section, [1] proposed to gradually introduce the information from the measurement. According to [2], this is done by carrying out the measurement update progressively over an artificial time $\displaystyle\gamma$, where $\displaystyle\gamma$ ranges in the (arbitrarily selected) interval222Other intervals could be used when more convenient. This also includes $\displaystyle\gamma\in[0,\infty]$. $\displaystyle\gamma\in[0,1]$.

To achieve the desired effect, we parametrize the likelihood with the artificial time $\displaystyle\gamma$, resulting in $\displaystyle f_{L}(\underline{x},\gamma)$. For the initial value of $\displaystyle\gamma$, $\displaystyle\gamma=0$, an uninformative likelihood, e.g., $\displaystyle f_{L}(\underline{x},\gamma=0)=1$, is selected that simply returns the prior unchanged. For the final value of $\displaystyle\gamma$, $\displaystyle\gamma=1$, we want to recover the original likelihood function. In between the extreme values for $\displaystyle\gamma$, i.e., $\displaystyle\gamma\in(0,1)$, we desire an interpolation between the constant likelihood and the original likelihood that is convenient to use and leads to the desired estimation results.

In summary, see also [6, p. 8], the progressive likelihood function $\displaystyle f_{L}(\underline{x},\gamma)$ parametrized by the artificial time $\displaystyle\gamma$ is given by

[TABLE]

Plugging the parametrized likelihood into the Bayesian update equation (7) yields the true parametrized posterior

[TABLE]

There are infinitely many progression schedules and, strictly speaking, a good schedule depends on the specific filtering setup, i.e., the prior density and the original likelihood. An intuitive and generic progression schedule is attained by taking the power of the original likelihood as

[TABLE]

This obviously reduces to $\displaystyle f_{L}(\underline{x},\gamma=0)=1$ and $\displaystyle f_{L}(\underline{x},\gamma=1)=f_{L}(\underline{x})$. As typically most of the change in the posterior caused by changes in $\displaystyle\gamma$ is experienced for small $\displaystyle\gamma$, a generalization is to use nonlinear functions of $\displaystyle\gamma$ with $\displaystyle g(\gamma=0)=0$ and $\displaystyle g(\gamma=1)=1$ that start out with a smaller slope around zero such as $\displaystyle g(\gamma)=\gamma^{2}$. The exponentiation with such a function gives a parametrized likelihood of the form

[TABLE]

Its derivative is

[TABLE]

In many applications, we assume additive Gaussian measurement noise, which leads to likelihoods from the exponential family. For these likelihoods, we obtain simpler expressions due to the logarithm taken above. For additive zero-mean Gaussian noise with covariance matrix $\displaystyle{\mathbf{C}}_{v}$ and $\displaystyle g(\gamma)=\gamma$, the likelihood function corresponding to (2) is

[TABLE]

IV Newton Flow Induced by Distance Measure

In a progressive Bayesian particle flow, a given intermediate posterior density undergoes continual changes for increasing artificial time $\displaystyle\gamma$ under the influence of the homotopic changes of the likelihood function. Hence, we require a distance measure that matches a slightly (infinitesimally) changed posterior density with an approximation from the desired density class. For particle sets, infinitesimal weight changes have to be instantly compensated by suitable particle location changes as we want to keep all weights equal.

We start with the definition of an abstract distance measure. The ODE for describing the flow will be derived in two steps. In the first step, we will derive a so-called iterative flow where the given true posterior remains unchanged over the flow. Just the approximative posterior is changed by the flow. In the second step, a so-called recursive ODE is derived, where the reference density is replaced by the approximation at either discrete time steps or continually.

IV-A Abstract Distance Measure

We are given a true or reference posterior density $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ explicitly depending on the artificial time $\displaystyle\gamma$ and an approximation $\displaystyle f_{e}(\underline{x},\underline{\eta}(\gamma))$ of the reference density. The parameter vector $\displaystyle\underline{\eta}(\gamma)$ is adjusted over time $\displaystyle\gamma$ so that the approximate posterior follows the evolving reference density. As a result, $\displaystyle\underline{\eta}(\gamma)$ implicitly depends on $\displaystyle\gamma$.

To calculate the necessary changes in the parameter vector $\displaystyle\underline{\eta}(\gamma)$, an abstract distance measure for comparing $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ and $\displaystyle f_{e}(\underline{x},\underline{\eta}(\gamma))$ is used

[TABLE]

Changing $\displaystyle\gamma$ induces variations in the distance that have to be compensated by changes in $\displaystyle\underline{\eta}(\gamma)$ in order to stay in the minimum.

Finding the optimal approximate posterior density by minimizing the distance measure (17) is generally difficult due to local minima. Hence, we start at the known minimum for $\displaystyle\gamma=0$, i.e., $\displaystyle f_{e}(\underline{x},\gamma=0)=\tilde{f}_{e}(\underline{x},\gamma=0)=f_{p}(\underline{x})$, and track the optimal $\displaystyle\gamma$-dependent approximate posterior $\displaystyle f_{e}(\underline{x},\underline{\eta}(\gamma))$ as the minimum of (17) with infinitesimal changes of the artificial time $\displaystyle\gamma$. The desired final approximate posterior is obtained at $\displaystyle\gamma=1$.

IV-B Derivation of Iterative ODE

The necessary condition for a minimum is obtained by calculating the gradient vector $\displaystyle\underline{G}(\underline{\eta}(\gamma),\gamma)$, i.e., the partial derivative of $\displaystyle D(\underline{\eta}(\gamma),\gamma)$ w.r.t. parameter vector $\displaystyle\underline{\eta}(\gamma)$, and setting it to zero

[TABLE]

As we are interested in the change of the parameter vector $\displaystyle\underline{\eta}(\gamma)$ w.r.t. $\displaystyle\gamma$, we take the partial derivative of $\displaystyle\underline{G}(\underline{\eta}(\gamma),\gamma)$ w.r.t. $\displaystyle\gamma$

[TABLE]

which is equivalent to

[TABLE]

Using the abbreviations under the braces, we denote the derivative of the gradient vector w.r.t. $\displaystyle\gamma$ as

[TABLE]

and the Hesse matrix as

[TABLE]

We can now rewrite (20) and obtain the ODE

[TABLE]

It is an implicit ODE with mass matrix $\displaystyle{\mathbf{H}}(\cdot,\cdot)$, where we have to solve for $\displaystyle\dot{\underline{\eta}}(\gamma)$. It can be formally rewritten as an explicit ODE

[TABLE]

Please note that we desist from explicitly inverting the Hessian to convert the ODE into an explicit one, as numerical inversion has a high computational complexity and undesirable numerical properties. We rather solve for $\displaystyle\dot{\underline{\eta}}(\gamma)$ in the linear system (23) while exploiting the properties of the Hessian.

We call this a Newton flow for the parameter vector $\displaystyle\underline{\eta}(\gamma)$ as a continuous version of Newton’s method, in analogy to a gradient flow being a continuous version of a gradient-descent method.

We also call this an iterative flow, as the reference posterior $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ at time $\displaystyle\gamma$ is maintained as the product of the prior particle set times the progressive likelihood $\displaystyle f_{L}(\underline{x},\gamma)$. As the product of prior density times likelihood typically is not in the same density class as the prior, the reference posterior $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ strays into the territory of impermissible densities and, for significant changes in $\displaystyle\gamma$, creeps beyond the border of usefulness. To make the effect on the flow in (24) more visible, we note that changes in $\displaystyle\gamma$ enter the right-hand side $\displaystyle\underline{R}(\cdot,\cdot)$ via the reference posterior $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ and its derivative $\displaystyle\dot{\tilde{f}}_{e}(\underline{x},\gamma)$ w.r.t. $\displaystyle\gamma$. This means that we can write $\displaystyle\underline{R}(\cdot,\cdot)$ as an explicit function of $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ and $\displaystyle\dot{\tilde{f}}_{e}(\underline{x},\gamma)$ as

[TABLE]

In order to avoid the reference posterior becoming useless, we have to replace it with its approximation $\displaystyle f_{e}(\underline{x},\underline{\eta}(\gamma))$ either from time to time or continually. We will call the resulting ODE a recursive flow.

IV-C Derivation of Recursive ODE

For deriving the recursive flow, we start with finite integration intervals333This is different from directly taking discrete steps in $\displaystyle\gamma$ as in [12] as we integrate the ODE over the interval.. Imagine that we first integrate the ODE from $\displaystyle\gamma=0$ to some small $\displaystyle\gamma_{\ast}$ so that the posterior flows from the prior to the best permissible approximation $\displaystyle f(\underline{x},\underline{\eta}(\gamma_{\ast}))$ of the (impermissible) true posterior $\displaystyle\tilde{f}_{e}(\underline{x},\gamma_{\ast})\propto f_{p}(\underline{x})f_{L}(\underline{x},\gamma_{\ast})$. For DMD, the weights of the reference posterior have changed and are now not equal anymore. To keep it from degenerating further, we replace the original reference posterior $\displaystyle\tilde{f}_{e}(\gamma_{\ast})$ by its best approximation $\displaystyle f(\underline{x},\underline{\eta}(\gamma_{\ast}))$ at time $\displaystyle\gamma_{\ast}$. From this time on, we compare with this newly reset reference density, that we then simply call reference posterior for the next steps. Of course, we have to modify the likelihood function that we use from $\displaystyle\gamma_{\ast}$ on as the new reference density already includes the likelihood integrated up to that time point. This involves calculating an effective likelihood $\displaystyle f_{L}^{\text{eff}}(\underline{x},\gamma)$ by dividing the original likelihood $\displaystyle f_{L}(\underline{x},\gamma)$ by the likelihood $\displaystyle f_{L}(\underline{x},\gamma_{\ast})$ used to far

[TABLE]

We could now repeat this procedure for the next $\displaystyle\gamma$ interval (and in fact, the filtering method could be implemented this way). However, we want to avoid the reference density leaving the space of permissible densities at all. This is achieved by simply using infinitesimal integration intervals. This is equivalent to replacing the reference posterior by its best approximation for all $\displaystyle\gamma$ in the ODE. Now the ODE does not include the reference density at all anymore

[TABLE]

We call this a recursive flow.

V Generalized Cramér Distance

Instead of the abstract distance measure assumed in the previous section, we now use the specific distance measure in (48) from the Appendix that compares two DMDs that may have (i) different numbers of components $\displaystyle M$ and $\displaystyle L$, (ii) at mutually distinct locations, (iii) with non-equal weights. For the particle flow in this paper, we use a specific configuration of this distance measure that compares two particle sets having the same numbers of samples ($\displaystyle M=L$).

V-A Adaptation of General Distance Measure to Particle Flow

Let us define the given posterior density at artificial time $\displaystyle\gamma\in[0,1]$ as $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$. When $\displaystyle\tilde{f}_{e}(\underline{x},\gamma)$ represents a particle set, the dependency on $\displaystyle\gamma$ stems from its $\displaystyle\gamma$-dependent weights. This dependency is caused by the parametrized likelihood $\displaystyle f_{L}(\underline{x},\gamma)$ as defined in Sec. III used for weighting the samples in the Bayesian update step. From (8), by plugging in $\displaystyle f_{L}(\underline{x},\gamma)$, we obtain

[TABLE]

with

[TABLE]

and

[TABLE]

We have constant $\displaystyle\tilde{\underline{x}}_{i}=\underline{x}_{p,i}$ and normalized time-varying weights

[TABLE]

with

[TABLE]

It is obvious, that for $\displaystyle\gamma>0$, the weights $\displaystyle\tilde{w}_{i}(\gamma)$ become non-equal.

As we desire to maintain equal weights, we define an approximating posterior density

[TABLE]

with fixed and equal weights $\displaystyle w_{i}=1/L$ and $\displaystyle\gamma$-dependent locations $\displaystyle\underline{x}_{i}(\gamma)$. We collect the locations in a parameter vector $\displaystyle\underline{\eta}(\gamma)$ with

[TABLE]

and the $\displaystyle\gamma$-independent weights in a vector

[TABLE]

We define the distance $\displaystyle D$ between the true density $\displaystyle\tilde{f}_{e}$ and its approximation $\displaystyle f_{e}$ based on the distance measure (48) as

[TABLE]

The best approximation of the posterior that maintains equal weights is tracked as the minimum of $\displaystyle D$ over artificial time $\displaystyle\gamma$, starting with the posterior density for $\displaystyle\gamma=0$, i.e., $\displaystyle f_{e}\bigl(\underline{x},\underline{w},\underline{\eta}(0)\bigr)=f_{p}(\underline{x})$.

V-B ODE* for Iterative Particle Flow*

During the particle flow, the densities $\displaystyle\tilde{f}_{e}$ and $\displaystyle f_{e}$ have the same number of components, so we have $\displaystyle M=L$. In that case, the gradient in (54) from the Appendix simplifies to

[TABLE]

We note that $\displaystyle\tilde{w}_{i}$ is a function of $\displaystyle\gamma$, i.e., $\displaystyle\tilde{w}_{i}=\tilde{w}_{i}(\gamma)$, while $\displaystyle\tilde{\underline{x}}_{i}$ and $\displaystyle w_{i}$ are independent of $\displaystyle\gamma$, and $\displaystyle\underline{\eta}_{i}$ only implicitly depends on $\displaystyle\gamma$. The derivative of the gradient vector w.r.t. $\displaystyle\gamma$ is given by

[TABLE]

with elements

[TABLE]

The normalized weights $\displaystyle\tilde{w}_{i}(\gamma)$ induced by the parametrized likelihood $\displaystyle f_{L}(\underline{x},\gamma)$ are given in (31) and (32). The derivatives of the normalized weights w.r.t. artificial time $\displaystyle\gamma$ are

[TABLE]

for $\displaystyle i=1,2,\ldots,L$ with

[TABLE]

V-C ODE* for Recursive Particle Flow*

In analogy to the procedure in Subsec. IV-C, we now derive a recursive version of the particle flow. This leads to simplifications of the derivative of the gradient vector w.r.t. $\displaystyle\gamma$ from (38). Its elements are now given by

[TABLE]

The derivatives of the normalized weights w.r.t. artificial time $\displaystyle\gamma$ are simplified to

[TABLE]

for $\displaystyle i=1,2,\ldots,L$.

The Hesse matrix in (56) from the Appendix is now also simplified. The diagonal block matrices are given by

[TABLE]

for $\displaystyle k=1,\ldots,L$ and the off-diagonal block matrices are still given by

[TABLE]

for $\displaystyle k=1,\ldots,L$, $\displaystyle l=1,\ldots,L$, $\displaystyle l\neq k$.

VI Conclusions

The proposed Newton-flow particle filter (NFPF) is a systematic approach to recursively estimate the hidden state of a nonlinear dynamic system. No heuristics are involved. The NFPF is easy to understand and simple to use as almost no tuning factors are required. It is fast for several reasons: (i) The number of particles to be processed is much smaller than in a standard particle filter. (ii) The distance measure, its gradient, and the Hesse matrix can be calculated in closed form. (iii) No optimization is required, just the solution of an ODE over a finite time interval. The proposed estimator shows a high performance compared to state-of-the-art approaches, produces unbiased high-quality estimates, is robust, is degeneration-free by design, and works in high dimensions. The proposed plugin replacement for a particle filter step is fully differentiable. This allows for end-to-end learning of the measurement nonlinearities or the noise densities, which will be exploited in future work.

Appendix

For convenience, we summarize the set distance measure derived in [7, 10] and list all the required formulas in the slightly different notation used in this paper.

VI-A Distance Measure

We are given two DMDs that we want to compare. Usually, this comparison is used for the approximation of one DMD by the other. The reference DMD is called $\displaystyle\tilde{f}(\underline{x})$. It comprises $\displaystyle M$ components with weights $\displaystyle\tilde{w}_{1},\tilde{w}_{2},\ldots,\tilde{w}_{M}$ and locations

[TABLE]

for $\displaystyle j=1,\ldots,M$. The approximating DMD $\displaystyle f(\underline{x})$ has $\displaystyle L$ components with weights $\displaystyle w_{1},w_{2},\ldots,w_{L}$ and locations

[TABLE]

for $\displaystyle i=1,\ldots,L$.

The numbers of components $\displaystyle L$ and $\displaystyle M$ can be different. In reduction problems, the number of components of the approximating density is usually less than or equal to the number of components $\displaystyle M$ of the reference density. In some cases, we even have $\displaystyle L\ll M$. For the application of the distance in this paper, the numbers of components are equal, i.e., $\displaystyle M=L$.

The generalized Cramér-von Mises distance [7] is given by $\displaystyle\sqrt{D}$ with

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

$\displaystyle D_{E}$ can be viewed as a penalty term that ensures equal means and is given by the squared difference between the means of $\displaystyle\tilde{f}(\underline{x})$ and $\displaystyle f(\underline{x})$

[TABLE]

$\displaystyle K$ is a large constant for weighting the penalty term.

Please note that in contrast to [7], we added the small constant $\displaystyle d_{\text{min}}$ here to avoid singularities when two samples assume the same location.

VI-B Gradient

The gradient of $\displaystyle D$ in (48) w.r.t. $\displaystyle\underline{x}_{k}$ is given by444For minimization, we use the squared distance $\displaystyle D$ for simplicity.

[TABLE]

with

[TABLE]

VI-C Hesse Matrix

The Hesse matrix $\displaystyle{\mathbf{H}}$ of the distance measure $\displaystyle D$ is composed of blocks $\displaystyle{\mathbf{H}}_{kl}\in{\rm I\hskip-1.69998ptR}^{N\times N}$ for $\displaystyle k=1,\ldots,L$, $\displaystyle l=1,\ldots,L$ according to

[TABLE]

We have to distinguish between diagonal and off-diagonal block matrices. The diagonal block matrices are given by

[TABLE]

for $\displaystyle k=1,\ldots,L$ with

[TABLE]

where $\displaystyle{\mathbf{I}}_{N}$ is the $\displaystyle N$-dimensional identity matrix. The off-diagonal block matrices are given by

[TABLE]

for $\displaystyle k=1,\ldots,L$, $\displaystyle l=1,\ldots,L$, $\displaystyle l\neq k$. The Hesse matrix is symmetric and positive definite.