Variational Gaussian filtering
via Wasserstein gradient flows

Adrien Corenflos is funded by the Academy of Finland project 321891 (ADAFUME). Hany Abdulsamad is funded by the Finnish Center for Artificial Intelligence (FCAI).

Adrien Corenflos1, Hany Abdulsamad2

Department of Electrical Engineering and Automation

Aalto University, Finland

Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
We present a novel approach to approximate Gaussian and mixture-of-Gaussians filtering. Our method relies on a variational approximation via a gradient-flow representation. The gradient flow is derived from a Kullback–Leibler discrepancy minimization on the space of probability distributions equipped with the Wasserstein metric. We outline the general method and show its competitiveness in posterior representation and parameter estimation on two state-space models for which Gaussian approximations typically fail: systems with multiplicative noise and multi-modal state distributions.

Kalman filtering, variational inference, state-space models, Wasserstein gradient flow.

§ INTRODUCTION

State-space models (or hidden Markov models) are a class of models widely used to represent latent dynamics that are partially or indirectly observed. They typically arise in ecological, economical, and tracking applications <cit.>. Formally, state-space models are given by a set of dynamics and noisy observations, often depending on a set of parameters $\theta$

\begin{equation}\label{eq:ssm}
\begin{aligned}
X_0 & \sim p_0(\cdot \given \theta), \\
X_{k+1} & \sim p_k(\cdot \given X_k, \theta),\\
Y_k & \sim h_k(\cdot \given X_k, \theta).
\end{aligned}
\end{equation}

While the problem of inference in such models is generally intractable, computing the filtering distribution $p(x_k \given y_{0:k})$ can typically be done exactly if the state-space is finite ($x_k$ can only take a finite number of values) or when all the (conditional) densities in (<ref>) are Gaussian using the celebrated Kalman filter [2]. When this is not the case, relying on approximations becomes necessary. Two important types of approximations are approximate Gaussian filters <cit.>, and Monte Carlo filters <cit.>.

Although the standard filtering problem is important, one may also be interested in system identification, which, in the parametric context, refers to learning $\theta$ from a sequence of observations $\{y_{0}, \ldots, y_{K}\}$.

In this article, we pay particular attention to two classes of models, typically neglected in approximate Gaussian filtering. The first class is that of models with multiplicative noise, for which stochastic volatility models are an illustrative example, often used in economics to model financial returns [5]. These are usually given as an auto-regressive latent state $x_k$, and observations $y_k$ following

\begin{equation}\label{eq:stochvol}
\begin{aligned}
Y_k & = \exp(X_k/2) \, \eta_k, \\
X_{k+1} & = \mu + \alpha(X_{k} - \mu) + \sigma \, \epsilon_{k},
\end{aligned}
\end{equation}

where the noise processes are correlated

\begin{equation}\label{eq:corr-noise}
\begin{pmatrix}
\epsilon_k \\ \eta_k
\end{pmatrix}
\sim \N\left(
\begin{pmatrix}
0 \\ 0
\end{pmatrix},
\begin{pmatrix}
1 & \rho \\ \rho & 1
\end{pmatrix}
\right).
\end{equation}

The second class we consider are systems for which the state (filtering) posterior is a multi-modal distribution. A simple example of this form can be given by constructing a latent state $x_k$ with random walk dynamics while the observations $y_k$ are a function of  the absolute value of $x_k$:

\begin{equation}\label{eq:multimodal}
\begin{aligned}
X_{k+1} & = X_{k} + \epsilon_k, && \, \epsilon_k \sim \N(0, 1), \\
Y_k & = |X_k| + \eta_k, && \, \eta_k \sim \N(0, 1).
\end{aligned}
\end{equation}

If, for example, $X_0 \sim \N(0, \delta^2)$, it is straightforward to see that the state filtering distribution will be bi-modal and fully symmetric with respect to the x-axis.

§.§ Contributions

Existing approximate Gaussian filtering methods suffer from several drawbacks. The linearization methods of [6, 7] for example require computing conditional expectations $m_k(x) = \mathbb{E}\left[Y_k \given X_k=x\right]$. In the case of the stochastic volatility model (<ref>), this quantity will unequivocally be null (at least for $\rho = 0$, see Section <ref> for details). Consequently, applying these methods to (<ref>) will lead to a filtering (or a smoothing) solution independent of the gathered observations, which is problematic.

Additionally, classical linearization methods do not extend directly to multi-modal distributions, and they do not, to the best of our knowledge, enable handling mixtures of Gaussians. In view of this, our contributions are the following:

* We rephrase the filtering problem as an iterative distribution fitting problem.

* We apply the variational inference method from [8] to propagate Gaussian approximations between time steps and formulate it as a fixed point iteration for efficient gradient calculation.

* We use our method for parameter estimation in stochastic volatility models and filtering of multi-model target distributions.

§ METHODOLOGY

§.§ Variational Inference via Wasserstein Gradient-flows

Let $\pi(x) \propto \exp(-V(x))$ be an arbitrary target distribution known up to a normalizing constant. Given a variational family of distributions, $q_{\phi}(x)$, $\phi \in \Phi$, and a measure of discrepancy $L(\phi) = D(\pi, q_{\phi})$ between $\pi$ and $q_{\phi}$, it is natural to try and approximate $\pi$ with $q_{\phi}$ by minimizing $L$.

A typical discrepancy is given by the Kullback–Leibler (KL) divergence  [9]

\begin{equation}
L(\phi) = \kl(q_{\phi} \given \pi) \coloneqq \int q_{\phi}(x) \log \frac{q_{\phi}(x)}{\pi(x)} \, \d x.
\end{equation}

This divergence presents a number of attractive properties for statistical inference: (i) it is positive, (ii) it only requires evaluating $V(x)$ and does not necessitate knowing the normalizing constant of $\pi$, and (iii), it is exact in the sense that $\kl(q_{\phi} \given \pi) = 0$ if and only if $q_{\phi} = \pi$. For more details, we refer the reader to [10].

For example, if $q_{\phi}(x) = \N(x \given \mu, \Sigma)$ is in the family of well-defined Gaussians, parameterized by their mean and covariance, then the inference problem can be cast as a minimization of $L(\cdot)$ with respect to $(\mu, \Sigma)$. Furthermore, when ignoring the positivity constraints on $\Sigma$, it is possible to define a gradient flow on $\phi = (\mu, \Sigma)$, akin to a gradient descent on $L(\cdot)$ in continuous time

\begin{equation}\label{eq:euclidean-flow}
\frac{\d \phi_t}{\d t} = -\nabla_{\!\phi} L(\phi_t),
\end{equation}

which, under an assumption of convexity, will converge to the minimizer of the objective $L(\cdot)$ [11].

While this procedure is correct in essence, it targets the problem indirectly by first assuming an arbitrary parameterization of the model which may or may not respect convexity. A more direct approach is to fit $q_{\phi}$ to $\pi$ in terms of a minimization problem over the space of probability distributions, where we want to minimize $L(q) = D(\pi, q)$. In this case, it is possible to define an analog to the gradient flow (<ref>) by equipping the space of probability distributions with the Wasserstein distance <cit.>. Under this metric, we can define a trajectory of probability distributions $q_t(x) \in \mathcal{P}(\mathbb{R}^d)$ via a partial differential equation

\begin{equation}\label{eq:wasserstein-flow}
\frac{\partial q_t(x)}{\partial t} = \nabla \cdot \left[ q_t(x) \nabla \log \frac{q_t(x)}{\pi(x)} \right],
\end{equation}

where $\nabla \cdot$ is the divergence operator, expressed in Euclidean coordinates.

Interestingly, by restricting $q_t$ to represent a Gaussian distribution, it was shown in [8], following [13], that (<ref>) can be reformulated into coupled ordinary differential equations (ODEs) on the mean $\mu_t$ and covariance $\Sigma_t$ of $q_t$

\begin{equation}\label{eq:mean-cov-flow}
\begin{aligned}
\frac{\d \mu_{t}}{\d t} & = -\mathbb{E} \left[\nabla V(Z_t)\right] \\
\frac{\d \Sigma_t}{\d t} & =
\begin{aligned}[t]
& 2 \, I - \mathbb{E} \left[\nabla V(Z_t) \otimes (Z_t - \mu_t) \right] \\
& - \mathbb{E} \left[(Z_t - \mu_t) \otimes \nabla V (Z_t)\right],
\end{aligned}
\end{aligned}
\end{equation}

where $I$ is the identity matrix with dimensions $d \times d$ and $Z_{t} \sim \N(\mu_t, \Sigma_t)$ is Gaussian. Now, provided that we can compute or approximate the expectations arising in (<ref>), we can find a minimizer $q^{*}(x) = \N(x \given m, P)$ of $L(\cdot)$ by integrating these coupled ODEs until convergence.

§.§ Filtering as Variational Inference

The filtering problem involves inferring the posterior distribution $p(x_k \given y_{0:k})$ for each time step $k$. To do so, it is often possible to rely on the familiar innovation-prediction decomposition <cit.>

\begin{equation}\label{eq:bayes-filtering}
\begin{aligned}
p(x_k \given y_{0:k}) & \propto p(y_k \given x_k) p(x_k \given y_{0:k-1}), \\
p(x_k \given y_{0:k-1}) & = \! \int \! p(x_k \given x_{k-1}) p(x_{k-1} \given y_{0:k-1}) \, \d x_{k-1}.
\end{aligned}
\end{equation}

For simplicity, we will focus on state-space models with affine Gaussian transition models

\begin{equation}
p(x_k \given x_{k-1}) = \N(x_k \given A_{k-1} \, x_{k-1} + b_{k-1}, Q_{k-1}),
\end{equation}

covering those previously presented in (<ref>) and (<ref>).

In that case, assuming a Gaussian approximation is available for the previous time step $k-1$,

\begin{equation}
p(x_{k-1} \given y_{0:k-1}) \approx \N(x_{k-1} \given m_{k-1}, P_{k-1}),
\end{equation}

then the prediction step of the filter leads to another Gaussian distribution $p(x_k \given y_{0:k-1}) \approx \N(x_k \given \bm_{k}, \bP_{k})$ with mean and covariance

\begin{equation}
\begin{aligned}
\bm_k & = A_{k-1} \, m_{k-1} + b_{k-1}, \\
\bP_k & = A_{k-1} \, P_{k-1} \, A_{k-1}^{\top} + Q_{k-1}.
\end{aligned}
\end{equation}

Given this predictive distribution, the innovation step computes the approximate filtering distribution

\begin{equation}
p(x_k \given y_{0:k}) \propto \exp(-V(x_k)),
\end{equation}

where the potential function $V(x_k)$ is given as

\begin{equation}\label{eq:potential}
V(x_k) = -\log p(y_k \given x_k) - \log \N(x_k \given \bm_k, \bP_k).
\end{equation}

Consequently, in order to find the parameters of a variational Gaussian $\N(x_k \given m_k, P_k)$ that approximates the posterior $p(x_k \given y_{0:k})$, it suffices to follow the recipe from Section <ref> and integrate (<ref>) up to stationarity, starting from the predictive parameters $(\bm_k, \bP_k)$, or any other approximation of $(m_k, P_k)$.

When the transition $p(x_k \given x_{k-1})$ is not Gaussian, it is possible to use a similar variational approach to propagate the Gaussian approximation $p(x_{k-1} \given y_{0:k-1}) \approx \N(x_{k-1} \given m_{k-1}, P_{k-1})$ through the non-Gaussian dynamics and approximate the marginal  $p(x_{k} \given y_{0:k-1})$ with another Gaussian $\N(x_k \given \bm_k, \bP_k)$. This method was used in [13] to propagate the Gaussian approximation through dynamics defined by a stochastic differential equation. Combining this with our approach is, therefore, de facto possible. However, due to the limited scope, we leave this aspect for future work and only consider Gaussian dynamics in the remainder of this article.

Finally, because the likelihood of the observations is given by $p(y_{0:k}) = p(y_{0:k-1})\int p(y_k \given x_k) p(x_k \given y_{0:k-1}) \, \d x_k$, it is easy to derive an approximation of the marginal log-likelihood of the model by recursion. That is because the quantity $\ell_{k} = \int p(y_k \given x_k) p(x_k \given y_{0:k-1}) \, \d x_k$ is evaluated as part of (<ref>), and can therefore be used to provide an online estimate of the log-likelihood increments, which, in turn, can be used, for example, in a system identification scenario. We return to this point in Sections <ref> and <ref>.

§.§ The Multi-modal Filtering Case

To generalize the variational technique from the previous section, let us suppose that the filtering distribution at time $k-1$ is instead given by a mixture density

\begin{equation}
p(x_{k-1} \given y_{0:k-1}) = \sum_{i=1}^{N_i} w^{(i)} \, \N(x_{k-1} \given m_{k-1}^{(i)}, P_{k-1}^{(i)}),
\end{equation}

with $w^{(i)} = 1/N_i, \forall i \in [1, N_i]$. In this case, when the transition dynamics are Gaussian, it is easy to show that

\begin{equation}\label{eq:pred-mixture}
p(x_{k} \given y_{0:k-1}) = \sum_{i=1}^{N_i} w^{(i)} \, \N(x_k \given \bm_{k}^{(i)}, \bP_{k}^{(i)}),
\end{equation}

where $\forall i \in [1, N_i]$

\begin{equation}
\begin{aligned}
\bm_k^{(i)} & = A_{k-1} \, m_{k-1}^{(i)} + b_{k-1}, \\
\bP_k^{(i)} & = A_{k-1} \, P_{k-1}^{(i)} \, A_{k-1}^{\top} + Q_{k-1}.
\end{aligned}
\end{equation}

These updates correspond to a tractable prediction step. As a consequence, we only need to understand how to perform the $\emph{innovation}$ step to arrive at the posterior $p(x_{k} \given y_{0:k})$.

Conveniently, [8] also shows that the duality between the gradient flow (<ref>) and the coupled ODEs (<ref>) extends to the case of the finite variational mixture of Gaussians

\begin{equation}
q_t(x) = \sum_{i=1}^{N_i} \omega^{(i)} q_t^{(i)}(x)= \sum_{i=1}^{N_i} \omega^{(i)} \, \N(x \given \mu_t^{(i)}, \Sigma_t^{(i)}),
\end{equation}

with $\omega^{(i)} = 1/N_i, \forall i \in [1, N_i]$. In this case, rather than a pair of ODEs, we obtain a system of such ODEs

\begin{equation} \label{eq:mixture-mean-cov-flow}
\begin{aligned}
\frac{\d \mu_t^{(i)}}{\d t} & = - \mathbb{E} \left[\nabla \log \frac{q_t}{\pi}(Z_t^{(i)}) \right] \\
\frac{\d \Sigma_t^{(i)}}{\d t} & =
\begin{aligned}[t]
& - \mathbb{E} \left[\nabla^2 \log \frac{q_t}{\pi}(Z_t^{(i)})\right] \Sigma_t^{(i)} \\
& - \Sigma_t^{(i)} \mathbb{E} \left[\nabla^2 \log \frac{q_t}{\pi}(Z_t^{(i)}) \right],
\end{aligned}
\end{aligned}
\end{equation}

where, for all $i$, $Z_t^{(i)} \sim \N(\mu_t^{(i)}, \Sigma_t^{(i)})$ is Gaussian, and $\nabla^2$ denotes the Hessian operator.

This result means that, given a Gaussian mixture approximation of the predictive $p(x_k \given y_{1:k-1})$, we can obtain a Gaussian mixture approximation of $p(x_k \given y_{0:k})$ by integrating (<ref>) and following a similar approach to Section <ref>.

§.§ Numerical Considerations and Implementation

In practice, the integrals arising in (<ref>) and (<ref>) are not available in closed form. Therefore, we need to resort to numerical integration. This can be done by using any form of deterministic or stochastic Gaussian quadrature, for example, Monte Carlo <cit.> or sigma-points <cit.> methods. The two approaches have their pros and cons: Monte Carlo will give the correct solution on average, while deterministic quadrature is bound to be biased. Nonetheless, deterministic rules are sometimes more practical when no stochasticity should be allowed in the system. Whichever is chosen, we will obtain an approximation

\begin{equation}\label{eq:mean-cov-flow-approx}
\frac{\d \mu_t}{\d t} \approx F_m(\mu_t, \Sigma_t), \quad \frac{\d \Sigma_t}{\d t}\approx F_P(\mu_t, \Sigma_t),
\end{equation}

of (<ref>) for the choice of $V(x_k)$ given by (<ref>). Moreover, the same approximation scheme can be used to compute log-likelihood increments $\log \mathbb{E} \left[p(y_k \given x_k)\right]$ under the predictive $\N(x_k \given \bm_{k}, \bP_{k})$. In Section <ref>, we use Gauss–Hermite <cit.> quadrature integration rules with order five.

Uni-modal Wasserstein Gradient-flow Filter

[1]

* input: Measurements $y_{0:K}$ and prior $(m_0, P_0)$

* output: Filtering posterior distributions $(m_{1:K}, P_{1:K})$ and marginal log-likelihood $\ell = \log p(y_{0:K})$

Set $\ell \gets 0$
$k \gets 1$ to $K$
Set $\ell_k \gets \log \mathbb{E}_{x \sim \N(\bm_k, \bP_k)}p(y_k \given x_{k})$
Set $\ell \leftarrow \ell + \ell_k$ Log-likelihood
$(m_k, P_k) \gets (\bm_k, \bP_k)$
Not converged Innovation
$(m_k, P_k) \gets I(m_k, P_k)$
Set $\bm_{k+1} \gets A_{k} \, m_k + b_k$ Prediction
Set $\bP_{k+1} \gets A_{k} \, P_k \, A_{k}^{\top} + Q_{k}$

$(m_{1:K}, P_{1:K})$, $\ell$

We now assume that we have chosen an integration method $I$ for which the stationary solution of (<ref>) is a fixed point: $(m_k, P_k) = I(m_k, P_k)$. This can, for example, be an Euler integration step with a small step size. This fixed-point perspective allows us to leverage the implicit function theorem and bypass the loop when computing gradients for system identification. For the sake of brevity, we omit the details here and refer to [15] for the method and to our code[Implementation of the fixed-point iteration is available under <https://github.com/hanyas/wasserstein-flow-filter/blob/master/wasserstein_filter/utils.py>] for a Python implementation.

Therefore, the final algorithm for variational uni-modal filtering from Section <ref> is given by Algorithm <ref>. Due to space constraints, we do not reproduce the algorithm for the multi-modal case from Section <ref>: the procedure follows the same steps, albeit for a larger system of ODEs.

§ EMPIRICAL EVALUATION

As discussed in Section <ref>, introducing our method is motivated by the problems posed by multiplicative noise and multi-modality in approximate Gaussian filtering. Consequently, we aim to demonstrate its effectiveness on these two problems. We compare the variational Wasserstein filter (VWF) to a bootstrap particle filter (PF, [16]) using $500$ particles and the continuous resampling scheme from [17]. This resampling allows using the particle filter in a standard parameter estimation scenario. The code to reproduce these experiments can be found at <https://github.com/hanyas/wasserstein-flow-filter>.

§.§ Stochastic Volatility with Leverage

First, we consider the stochastic volatility model as given by (<ref>) and (<ref>).
Because the noises $\epsilon_k$ and $\eta_k$ are correlated, we need to construct the joint distribution over $X_k$ and its generating noise $\epsilon_k$ and form the augmented two-dimensional state
$\zeta_k =
\begin{bmatrix}
X_k & \epsilon_{k}
\end{bmatrix}^\top$.
The dynamics of $\zeta_k$ are then given by

\begin{equation}\label{eq:stacked-vol}
\begin{pmatrix}
X_{k+1} \\ \epsilon_{k+1}
\end{pmatrix} =
\begin{pmatrix}
\alpha & \sigma \\
0 & 0
\end{pmatrix}
\begin{pmatrix}
X_{k} \\ \epsilon_{k}
\end{pmatrix}
\begin{pmatrix}
\mu(1 - \alpha)\\
\end{pmatrix}
+ q_{k} \begin{pmatrix}
\end{pmatrix},
\end{equation}

where $q_{k}$ is a standard Gaussian random variable. With these dynamics, we have the following observation model

\begin{equation}
Y_k = \exp(X_k / 2) (\rho \epsilon_k + \sqrt{1 - \rho^2} \, r_k),
\end{equation}

with $r_k$ being a standard Gaussian random variable. Contrary to the original form of the model, the noise processes are now de-correlated so that we can apply our method to the 2-dimensional system defined by $\zeta_k$.

To perform an empirical comparison, we simulate the model and record three trajectories with $K=1\,000$, $K=1\,500$, and $K=2\,000$ observations. We use parameters $\alpha_*=0.975$, $\mu_* = 0.5$, $\sigma_*^2 = 0.02$, and $\rho_*=-0.8$. These values correspond to those in [17], typical of a standard stock market.

To illustrate the problem of using methods based on linearizing the conditional observation mean $\mathbb{E}\left[Y_k \given X_k, \epsilon_k\right]$, we evaluate the marginal log-likelihood of the data approximated by an extended Kalman filter (EKF) targeting (<ref>) with parameters $(\alpha_*, \mu_*, \sigma_*)$ while varying levels of correlation $\rho$ and number of observations $K$. The resulting curves for the particle filter of [17], our method, and the extended Kalman filter are shown in Figure <ref>. As the (model) correlation $\rho$ between the two noise-generating processes decreases, the predictive value of the observations becomes negligible. As a consequence, the marginal likelihood, as approximated by the EKF, will struggle to capture the right level of correlation. This means that calibrating the model using an extended (or similar) Kalman filter will result in an inconsistent estimate, at least for the parameter $\rho$.

To further confirm this hypothesis, we follow [17] and perform joint maximum likelihood parameter estimation given $K=1\,000$ observations. We compare our method to theirs and to an extended Kalman filter. The likelihood is optimized in all cases using the gradient of the log-likelihood, delivered by automatic reverse differentiation. In the particular case of the variational Wasserstein filter, differentiation is made efficient by our fixed point formulation of the variational calibration.

For a statistical comparison, we repeat this experiment over ten trials associated with different random seeds and, therefore, ten distinct sets of observations. We report the mean plus or minus one standard deviation of the parameter estimates for each algorithm. However, because the particle filter is inherently a stochastic method, even for a fixed set of observations, we perform $25$ sub-trials per set of observations and use the median as a parameter estimate per trial. The results are given in Table <ref> and confirm our intuition. The extended Kalman filter provides inconsistent solutions, while our method delivers estimates comparable to the estimates of the particle filter.

Comparing the (normalized) marginal log-likelihood as a function of the leverage parameter $\rho$ in a stochastic volatility model. We plot the estimates as returned by an extended Kalman filter (red), a bootstrap particle filter (blue), and a variational Wasserstein filter (green). The dotted, dashed, and solid lines correspond to different numbers of observations $K=1\,000$, $K=1\,500$, and $K=2\,000$. The true value is $\rho=-0.8$ (vertical line). The variational Wasserstein and particle filters deliver consistent approximations independent of the data size, while the extended Kalman filter does not.

Statistics for parameter estimation of the stochastic volatility model. Results are averaged over ten trials.

max width=

$\mu$     $\alpha$     $\sigma$     $\rho$

True     0.50     0.975     0.14     -0.80

PF [17] 	    $0.56 \, (\pm 0.08)$     $0.972 \, (\pm 0.007)$     $0.15 \, (\pm 0.02)$     $-0.85 \, (\pm 0.06)$

VWF 							    $0.56 \, (\pm 0.07)$     $0.972 \, (\pm 0.009)$     $0.15 \, (\pm 0.02)$     $-0.80 \, (\pm 0.04)$

EKF 							    $0.69 \, (\pm 0.33)$     $0.780 \, (\pm 0.590)$     $0.21 \, (\pm 0.25)$     $-0.58 \, (\pm 0.54)$

§.§ Multi-modal Example

We now turn to the multi-modal motivating example (<ref>). We simulate $K=500$ observations from the model, where we have taken the variance at origin to be $\delta^2 = 1$. We then perform variational filtering using the mixture version of Algorithm <ref> with $N_i=2$ mixture components.

Our goal is a correct representation of the filtering posterior. Thus we will not assess the performance in terms of root mean square error, which is inappropriate for such problems. Instead, we qualitatively report the resulting filtering distributions in Figure <ref>. There is visually hardly any difference between the particle filtering estimate and our method.

One caveat to this positive result is that the bi-modal variational Wasserstein filter tends to collapse when the two modes are too close to each other. It is still unclear whether this is a feature of the general method or the ODE solver, and further investigations are warranted.

Performing filtering on a multi-modal dynamical system. We compare the filtering result of a particle filter (top) with that of a variational Wasserstein filter equipped with a mixture-of-Gaussians posterior representation (bottom). Both filters capture the bi-modal posterior distribution (red and blue) inferred from the observations (green) induced by the true states (black).

§ CONCLUSION

In this article, we have presented a novel approach for approximate Gaussian filtering, which avoids the reliance on enabling assumptions [6], which are usually not amenable to multiplicative noise and not extendable to multi-modal distributions. Several questions remain open:

Assumptions We have assumed, for simplicity, that the transition dynamics are governed by affine Gaussian densities. Overcoming this limitation is an interesting research question. In fact, in [13], the coupled ODEs (<ref>) were originally introduced to propagate the Gaussian through a nonlinear stochastic differential equation.

Numerics The numerical scheme chosen here is, primarily for ease of exposition, different from that of [8], which uses an iterative method called JKO [11] after [18], instead of integrating the differential equation (<ref>) directly. It is unclear which approach is the best fit for filtering applications.

Mixture weights The weights of the mixture in the ODE (<ref>) are not allowed to vary. This modeling restriction needs to be relaxed. However, that may lead to identifiability issues (for example, a single Gaussian can be represented with a mixture of two Gaussians with different weights in infinite ways). Furthermore, introducing time-varying weights should be done carefully.

§ INDIVIDUAL CONTRIBUTIONS

The original idea and redaction of the article are due to Adrien Corenflos. Both authors contributed to the design of the methodology. The implementation and experiments are due to Hany Abulsamad.

[1]
A. H. Jazwinski, Stochastic Processes and Filtering Theory.1em
plus 0.5em minus 0.4emCourier Corporation, 2007.

[2]
R. E. Kalman, “A new approach to linear filtering and prediction problems,”
Journal of Basic Engineering, 1960.

[3]
S. Särkkä, Bayesian Filtering and Smoothing.1em plus
0.5em minus 0.4emCambridge University Press, 2013.

[4]
N. Chopin and O. Papaspiliopoulos, An Introduction to Sequential Monte
Carlo.1em plus 0.5em minus 0.4emSpringer, 2020.

[5]
S. L. Heston, “A closed-form solution for options with stochastic volatility
with applications to bond and currency options,” The Review of
Financial Studies, 1993.

[6]
A. F. Garcìa-Fernàndez, L. Svensson, M. R. Morelande, and S. Särkkä,
“Posterior linearization filter: Principles and implementation using sigma
points,” IEEE Transactions on Signal Processing, 2015.

[7]
F. Tronarp, Á. F. García-Fernández, and S. Särkkä, “Iterative
filtering and smoothing in nonlinear and non-Gaussian systems using
conditional moments,” IEEE Signal Processing Letters, 2018.

[8]
M. Lambert, S. Chewi, F. Bach, S. Bonnabel, and P. Rigollet, “Variational
inference via Wasserstein gradient flows,” arXiv preprint
arXiv:2205.15902, 2022.

[9]
J. M. Joyce, “Kullback–Leibler divergence,” in International
Encyclopedia of Statistical Science.1em plus 0.5em minus 0.4emSpringer, 2011.

[10]
K. P. Murphy, Probabilistic Machine Learning: An Introduction.1em plus 0.5em minus 0.4emMIT Press, 2022.

[11]
F. Santambrogio, “Euclidean, metric, and Wasserstein gradient flows: An
overview,” Bulletin of Mathematical Sciences, 2017.

[12]
C. Villani, Optimal Transport: Old and New.1em plus 0.5em minus
0.4emSpringer Berlin Heidelberg, 2009.

[13]
S. Särkkä, “On unscented Kalman filtering for state estimation of
continuous-time nonlinear systems,” IEEE Transactions on Automatic
Control, 2007.

[14]
G. Pagès, Numerical Probability.1em plus 0.5em minus
0.4emSpringer, 2018.

[15]
B. Christianson, “Reverse accumulation and attractive fixed points,”
Optimization Methods & Software, 1994.

[16]
N. J. Gordon, D. J. Salmond, and A. F. Smith, “Novel approach to
nonlinear/non-Gaussian Bayesian state estimation,” in IEEE
Proceedings of Radar and Signal Processing, 1993.

[17]
S. Malik and M. K. Pitt, “Particle filters for continuous likelihood
evaluation and maximisation,” Journal of Econometrics, 2011.

[18]
R. Jordan, D. Kinderlehrer, and F. Otto, “The variational formulation of the
Fokker–Planck equation,” SIAM Journal on Mathematical Analysis,