Discrete gradients for computational Bayesian inference

TL;DR

This paper explores the use of discrete gradient methods to efficiently solve stiff differential equations arising in continuous-time Bayesian inference, specifically in ensemble Kalman-Bucy filters and Fokker-Planck discretizations.

Contribution

It introduces discrete gradient techniques for Bayesian inference, providing a comparative analysis with existing numerical methods for stiff differential equations.

Findings

Discrete gradient methods improve numerical stability.

They outperform semi-implicit methods in certain scenarios.

Enhanced efficiency in Bayesian filtering applications.

Abstract

In this paper, we exploit the gradient flow structure of continuous-time formulations of Bayesian inference in terms of their numerical time-stepping. We focus on two particular examples, namely, the continuous-time ensemble Kalman-Bucy filter and a particle discretisation of the Fokker-Planck equation associated to Brownian dynamics. Both formulations can lead to stiff differential equations which require special numerical methods for their efficient numerical implementation. We compare discrete gradient methods to alternative semi-implicit and other iterative implementations of the underlying Bayesian inference problems.

Tables

Table 1. Table 1. RMSE for sequential data assimilation applied to the Lorenz-63 model. The parameter α = 1 𝛼 1 \alpha=1 corresponds to a standard square root ensemble Kalman filter, while α < 1 𝛼 1 \alpha<1 results in a Gaussian mixture approximation to the prior and posterior distributions. Small improvements over the ensemble Kaman filter can be found for ensemble sizes M ≥ 35 𝑀 35 M\geq 35 .

$\alpha \backslash M$	$15$	$20$	$25$	$30$	$35$	$40$	$45$	$50$
$0.8$	$2.5145$	$2.4726$	$2.4471$	$2.4406$	$2.4221$	$2.4173$	$2.4162$	$2.4112$
$0.85$	$2.5048$	$2.4585$	$2.4357$	$2.4330$	$2.4186$	$2.4161$	$2.4020$	$2.4026$
$0.9$	$2.4989$	$2.4493$	$2.4396$	$2.4259$	$2.4177$	$2.4103$	$2.4106$	$2.4075$
$0.95$	$2.4998$	$2.4605$	$2.4399$	$2.4327$	$2.4242$	$2.4132$	$2.4162$	$2.4131$
$1.0$	$2.4804$	$2.4597$	$2.4518$	$2.4354$	$2.4334$	$2.4247$	$2.4244$	$2.4203$

Equations

$$\mathrm{d}{z}_{\tau }=-A(z_{\tau }){\nabla }_{z}V(z_{\tau })\mathrm{d}\tau ,$$

$$\mathrm{d}V(z_{\tau })=-{\nabla }_{z}V(z_{\tau })\cdot A(z_{\tau }){\nabla }_{z}V(z_{\tau })\mathrm{d}\tau \le 0$$

$${\pi }_{\ast }(x)\propto \pi (y|x){\pi }_0(x),$$

$$\mathbb{E}[g(X_{\ast })]\approx \frac{1}{M}\sum _{i=1}^{M}g(x_{\tau }^i).$$

$$\mathrm{d}{x}_{\tau }^i=-P_{\tau }^{xx}{\nabla }_{x^i}V(\{x_{\tau }^j\})\mathrm{d}\tau ,i=1,\dots ,M,$$

$$P_{\tau }^{xx}=\frac{1}{M-1}\sum _{j=1}^M(x_{\tau }^j-{\overline{x}}_{\tau })(x_{\tau }^j-{\overline{x}}_{\tau }{)}^{\mathrm{T}},{\overline{x}}_{\tau }=\frac{1}{M}\sum _{j=1}^M{x}_{\tau }^j$$

$$V(\{x_{\tau }^j\})=\frac{M}{2}\left\{S({\overline{x}}_{\tau })+\frac{1}{M}\sum _{i=1}^{M}S(x_{\tau }^i)\right\}.$$

$$S(x)=-\log \pi (y|x).$$

$$y=h(x)+R^{1/2}\Xi ,\Xi \sim \mathrm{N}(0,I),$$

$$S(x)=\frac{1}{2}(h(x)-y{)}^{\mathrm{T}}{R}^{-1}(h(x)-y).$$

$$z_{\tau }=\left(\begin{array}{c}{x}_{\tau }^1\\ x_{\tau }^2\\ ⋮\\ x_{\tau }^M\end{array}\right)$$

$$\mathrm{d}{x}_{\tau }^i=-P_{\tau }^{xh}{R}^{-1}\left(\frac{1}{2}(h(x_{\tau }^i)+{\overline{h}}_{\tau })-y\right)$$

$$P_{\tau }^{xh}=\frac{1}{M-1}\sum _{j=1}^M(x_{\tau }^j-{\overline{x}}_{\tau })(h(x_{\tau }^j)-{\overline{h}}_{\tau }{)}^{\mathrm{T}},{\overline{h}}_{\tau }=\frac{1}{M}\sum _{j=1}^{M}h(x_{\tau }^j).$$

$${\tilde{\pi }}_{\tau }(x)=\frac{1}{M}\sum _{j=1}^M\psi (x-x_{\tau }^j)=\int \psi (x-x^{\prime }){\pi }_{\tau }(x^{\prime })\mathrm{d}{x}^{\prime },$$

$${\pi }_{\tau }(x)=\frac{1}{M}\sum _{i=1}^M\delta (x-x_{\tau }^i)$$

$$\int \psi (x)\mathrm{d}x=1.$$

$$V(\{x_{\tau }^j\})$$

$$=\frac{1}{M}\sum _{j=1}^M\left\{\log {\tilde{\pi }}_{\tau }(x_{\tau }^j)-\log {\pi }_{\ast }(x_{\tau }^j)\right\}$$

$$=\frac{1}{M}\sum _{j=1}^M\left\{\log \left(\frac{1}{M}\sum _{l=1}^M\psi (x_{\tau }^j-x_{\tau }^l)\right)-\log {\pi }_{\ast }(x_{\tau }^j)\right\}$$

$$\mathrm{d}{x}_{\tau }^i=-M{\nabla }_{x^i}V(\{x_{\tau }^j\})\mathrm{d}\tau ,$$

$${\nabla }_{x_{\tau }^i}V(\{x_{\tau }^j\})=\frac{1}{M}\left\{{\nabla }_{x_{\tau }^i}\log {\tilde{\pi }}_{\tau }(x_{\tau }^i)+\frac{1}{M}\sum _{j\ne i}\frac{1}{{\tilde{\pi }}_{\tau }(x_{\tau }^j)}{\nabla }_{x_{\tau }^i}\psi (x_{\tau }^i-x_{\tau }^j)-{\nabla }_{x_{\tau }^i}\log {\pi }_{\ast }(x_{\tau }^i)\right\}.$$

$${\rm KL}_{\mathcal{H}}(\pi_{\tau}|\pi^{\ast})=\Bigl{\langle}\pi_{\tau},\log\left(\frac{\widetilde{\pi}_{\tau}}{\pi_{\ast}}\right)\Bigr{\rangle},$$

$$\frac{\delta {\mathrm{KL}}_{\mathcal{H}}({\pi }_{\tau }|{\pi }^{\ast })}{\delta {\pi }_{\tau }}(x)=\log {\tilde{\pi }}_{\tau }(x)-\log {\pi }_{\ast }(x)+\int \psi (x-x^{\prime })\frac{{\pi }_{\tau }(x^{\prime })}{{\tilde{\pi }}_{\tau }(x^{\prime })}\mathrm{d}{x}^{\prime }$$

$$M{\nabla }_{x_{\tau }^i}V(\{x_{\tau }^j\})={\nabla }_x\left(\frac{\delta {\mathrm{KL}}_{\mathcal{H}}({\pi }_{\tau }|{\pi }^{\ast })}{\delta {\pi }_{\tau }}\right)(x_{\tau }^i).$$

$${\nabla }_x\left(\frac{\delta {\mathrm{KL}}_{\mathcal{H}}({\pi }_{\infty }|{\pi }^{\ast })}{\delta {\pi }_{\infty }}\right)(x_{\infty }^i)=0$$

$${\pi }^{\ast }(x)\propto {\tilde{\pi }}_{\infty }(x)e^{l(x)},l(x)=\int \psi (x-x^{\prime })\frac{{\pi }_{\infty }(x^{\prime })}{{\tilde{\pi }}_{\infty }(x^{\prime })}\mathrm{d}{x}^{\prime },$$

$$\frac{\delta {\mathrm{KL}}_{\mathcal{H}}({\pi }_{\infty }|{\pi }^{\ast })}{\delta {\pi }_{\tau }}(x)=\text{mbox}const.$$

$$\mathrm{d}{X}_{\tau }={\nabla }_x\log {\pi }_{\ast }(X_{\tau })\mathrm{d}\tau +\sqrt{2}\mathrm{d}{W}_{\tau }$$

$$\underset{\tau \to \infty }{\lim }{p}_{\tau }={\pi }_{\ast }$$

$$\psi (x-x^{\prime }):=\mathrm{n}(x;x^{\prime },B)\propto \tilde{\psi }(x-x^{\prime }):=\exp \left(-\frac{1}{2}(x-x^{\prime }{)}^{\mathrm{T}}{B}^{-1}(x-x^{\prime })\right)$$

$${\pi }_{\mathrm{em}}(x)=\frac{1}{M}\sum _{i=1}^M\delta (x-{\hat{x}}_0^i)$$

Full text

1. Introduction

A number of algorithmic approaches to Bayesian inference (Robert, 2001; Tarantola, 2005; Kaipio & Somersalo, 2005) can be rephrased in terms of evolution equations of the form

[TABLE]

$\tau\geq 0$, with $A(z)$ being symmetric positive semi-definite and $V:\mathbb{R}^{N_{z}}\to\mathbb{R}$ being an appropriate potential. It holds that

[TABLE]

along solutions $z_{\tau}\in\mathbb{R}^{N_{z}}$ of (1). ODEs of the form (1) are referred to as gradient dynamics or of gradient flow structure (Stuart & Hymphries, 1996). Generally speaking, numerical approximations of (1) will not preserve the decay property (2) unless the time-step $\Delta\tau>0$ is chosen sufficiently small. Such time-step restrictions can be problematic if the ODE system (1) is stiff (Ascher, 2008). Discrete gradient methods (McLachlan et al., 1999) are attractive in this context, as they have been devised such that (2) holds regardless of the chosen step-size $\Delta\tau$. However, these methods are implicit, which implies an increased computational cost per time-step. In this paper, we explore the pros and cons of discrete gradient methods in the context of computational Bayesian inference. More specifically, we consider two particular examples of such gradient dynamics arising from particle methods for sequential data assimilation (Law et al., 2015; Reich & Cotter, 2015) and from a transformation approach to posterior sampling (Crisan & Xiong, 2010; Daum & Huang, 2011; Reich, 2011). In both cases, the decay property (2) plays a crucial role in ensuring the statistical consistency of the associated inference method when the state variable $z_{\tau}$ is treated as a random variable with the initial $z_{0}$ following a given prior distribution $\pi_{0}$.

The paper is structured as follows. Section 2 provides the necessary background for the two considered continuous-time gradient flow systems, namely, the ensemble Kalman–Bucy filter (EnKBF) and a particle discretisation of the Fokker–Planck equations associated to Brownian dynamics. Section 3 summarises the discrete gradient method and its numerical implementation. We also discuss the semi-implicit Euler method as an alternative time-stepping method. Numerical results and a comparison between different methods is provided in Section 4. The paper closes with a summary in Section 5.

2. Examples of gradient systems for Bayesian inference

The following two examples of gradient flow systems are both concerned with dynamically transforming a set of Monte Carlo samples $\{x^{i}_{0}\}$, $i=1,\ldots,M$, from a prior distribution $\pi_{0}(x)$ into samples $\{x^{i}_{\tau}\}$ from a posterior distribution

[TABLE]

which can be used to approximate expectation values of a function $g$ of a random variable $X_{\ast}$ with distribution given by (3), that is:

[TABLE]

The first example, namely, the ensemble Kalman–Bucy filter (Bergemann & Reich, 2010a, b, 2012), achieves this for $\tau=1$, while the second example, based on a particle discretisation of the Fokker–Planck equation (Russo, 1990; Degond & Mustieles, 1990; Reich, 2019), requires $\tau\to\infty$. We now describe both approaches in more detail.

2.1. Ensemble Kalman–Bucy filter

We have $M$ particles $x_{\tau}^{i}\in\mathbb{R}^{N_{x}}$ that evolve in time according to the gradient dynamics (Bergemann & Reich, 2010a, b)

[TABLE]

with the symmetric positive semi-definite (covariance) matrix $P_{\tau}^{xx}$ given by

[TABLE]

and the potential $V:\mathbb{R}^{N_{x}\times M}\to\mathbb{R}$ by111Different variants of $V$ can be used when the EnKBF is used for minimisation only. See, for example, Kovachki & Stuart (2018); Reich (2019).

[TABLE]

Here $S(x)$ denotes the negative log-likelihood function, that is,

[TABLE]

We are primarily interested in data models of the form

[TABLE]

and, hence,

[TABLE]

The ODE system (5) is typically solved over a unit time interval, that is, $\tau\in[0,1]$. The initial particles $\{x_{0}^{i}\}$ represents the prior distribution while the particles $\{x_{\tau}^{i}\}$ at final time approximate the posterior distribution $\pi_{\ast}$ and can be used in (4) with $\tau=1$.

The ODE system (5) is of the form (1) with

[TABLE]

and $A(z_{\tau})$ being a block diagonal matrix with repeated entries $P_{\tau}^{xx}$. Numerical problems arising from the stiffness of the gradient dynamics have been discussed by Amezcua et al. (2014).

We mention that there are derivative-free extensions of the EnKBF to non-linear forward maps $h$ (Bergemann & Reich, 2012; Reich & Cotter, 2015). Consider, for example,

[TABLE]

with

[TABLE]

This formulation is no longer of gradient flow structure unless $h$ is linear. There are also stochastic formulations of the EnKBF (Bergemann & Reich, 2012; Law et al., 2015; Reich & Cotter, 2015).

We finally mention that EnKBF is closely related to iterative implementations of the ensemble Kalman filter as considered, for example, by Emerik & Reynolds (2012); Sakov et al. (2012); Chen & Oliver (2013). Those implementations can also be thought of particular time-stepping methods for the continuous-time EnKBF. We will return to this aspect in Section 4.1. Iterative implementations of the Kalman–Bucy filter are also related to the natural gradient approach from statistical learning (Ollivier, 2018).

2.2. Particle flow Fokker–Planck dynamics

Given $M$ particles $x_{\tau}^{i}$ and a kernel function $\psi(x)\geq 0$ of a Reproducing Kernel Hilbert Space (RKHS) $\mathcal{H}$, we approximate their density in $\mathcal{H}$ by

[TABLE]

where $\pi_{\tau}$ stands for the empirical measure

[TABLE]

and we have assumed that $\psi(x)=\psi(-x)$ as well as

[TABLE]

Let us denote the standard $L_{2}$-inner product of two functions $f$ and $g$ by $\langle f,g\rangle$ and the associated inner product in $\mathcal{H}$ by $\langle f,g\rangle_{\mathcal{H}}$. The Kullback–Leibler divergence between $\widetilde{\pi}_{\tau}(x)$ and the posterior target PDF $\pi_{\ast}(x):=\pi(x|y)$ is then approximated in $\mathcal{H}$ by

[TABLE]

and the particle flow Fokker–Planck dynamics is given by

[TABLE]

for $i=1,\ldots,M$. See Reich (2019) for more details. The gradient of $V$ with respect to $x_{\tau}^{i}$ is given by

[TABLE]

We note that one can replace the normalised kernel $\psi$ by an unnormalised version $\widetilde{\psi}$ without changing the gradient dynamics. Furthermore, the ODE system (20) is of the form (1) with $z_{\tau}$ defined by (11) and $A(z_{\tau})=MI$ where $I$ is the identity matrix.

Let us look at the evolution equations (20) and their geometric properties in some more detail. The variational derivative of the Kullback–Leibler divergence in $\mathcal{H}$,

[TABLE]

with respect to $\pi_{\tau}$ is given by

[TABLE]

and, hence, one finds that

[TABLE]

Furthermore, as $\tau\to\infty$, the particle flow Fokker–Planck dynamics (20) leads to

[TABLE]

and the asymptotic particle positions $\{x_{\infty}^{i}\}$ can be used to approximate expectation values with respect to $X_{\ast}\sim\pi_{\ast}$ based on

[TABLE]

which follows from (23) for an equilibrium density $\widetilde{\pi}_{\infty}$, that is,

[TABLE]

We use the simpler, less accurate in the limit $M\to\infty$ empirical approximation (4) for $\tau>0$ sufficiently large in the numerical experiments of Section 4.

For $M\to\infty$ and the kernel function $\psi(x)$ approaching a Dirac delta function, the evolution equations (20) provide a mesh-free particle approximation to the Fokker–Planck equation for the marginal PDFs $p_{t}$ of Brownian dynamics

[TABLE]

(Reich & Cotter, 2015). It holds under suitable conditions on $\pi_{\ast}$ that the marginal distributions $p_{\tau}$ of the solutions $X_{\tau}$ to the SDE (28) satisfy

[TABLE]

in a weak sense (Pavliotis, 2014).

As a specific example, let us consider the RKHS $\mathcal{H}$ with Gaussian kernel

[TABLE]

for given covariance matrix $B$. Here ${\rm n}(x,\overline{x},B)$ denotes the probability density function of a Gaussian random variable with mean $\overline{x}$ and covariance matrix $B$. We now provide a construction for the required covariance matrix $B$ and the initial particle positions $x_{0}^{i}$ that yields a kernel approximation to the given prior distribution $\pi_{0}$. Given a set of $M$ samples $\widehat{x}^{i}_{0}$ from the prior distribution $\pi_{0}$, the associated empirical measure is given by

[TABLE]

with empirical mean $\overline{x}_{0}$ and empirical covariance matrix $P^{xx}_{0}$. Setting

[TABLE]

for $\alpha\in(0,1]$ in (30), the associated measure (14) in $\mathcal{H}$ is defined as follows:

[TABLE]

with particle positions

[TABLE]

In other words, the covariance matrix $B$ and the particle positions $x^{i}_{0}$ are chosen such that the mean and covariance matrix under $\widetilde{\pi}_{0}$ are identical to $\overline{x}_{0}$ and $P_{0}^{xx}$, respectively, that is,

[TABLE]

and

[TABLE]

The initial particle positions $\{x_{0}^{i}\}$ are now evolved under the ODEs (20). Note that $\alpha=1$ leads to $x^{i}_{0}=\overline{x}_{0}$ and points of $\nabla_{x^{i}}V=0$, that is, stationary points of (20) agree with critical points of the posterior distribution $\pi_{\ast}$. It is possible to extend the dynamics in the particles $\{x_{\tau}^{i}\}$ by an evolution equation in the kernel matrix $B$ (Yserentant, 1997).

Furthermore, one can put (20) into the more general form of (1), that is,

[TABLE]

for $i=1,\ldots,M$. A particular choice is given by $\mathcal{A}(\{x_{\tau}^{j}\})=M\,P_{t}^{xx}$. A related preconditioned Brownian dynamics formulation (28) has been proposed recently by Garbuno-Inigo et al. (2019).

We finally mention that the Stein variational approach (Liu & Wang, 2016) has recently become popular as an alternative dynamical approach for transforming a set of Monte Carlo samples from a prior distribution into samples from a posterior distribution. However, it is currently unknown whether the resulting interacting particle systems possess a gradient flow structure. See Detommaso et al. (2018) for a preconditioned Stein variational formulation in the spirit of (37).

3. Discrete gradient time-stepping

Discrete gradient methods are numerical time-stepping methods for ODEs of type (1) that allow one to maintain the decay property (2) under arbitrary choices of the step-size $\Delta\tau>0$. See McLachlan et al. (1999) for an overview of such methods. Here we focus on the following particular form of discrete gradient dynamics:

[TABLE]

$\theta\in[0,1]$, with the discrete gradient

[TABLE]

and symmetric positive definite matrix $A(z)$. It follows that

[TABLE]

as desired. Note that (40) is a rather natural discrete representation of the continuous (2).

There exists several variants of the method, for example, the time-symmetric formulation of Gonzalez (1996) (that is, with $\theta=1/2$) and the local update version of Reich (1996).

We rewrite (38) in the form

[TABLE]

with the scalar factor $\gamma_{n}$ defined by

[TABLE]

Hence one can think of the discrete gradient method (38) as a $\theta$-method (Stuart & Hymphries, 1996) with an implicitly determined effective step-size

[TABLE]

Formulation (41) with $\theta>0$ leads naturally to the following iterative solution procedure at each time-step $\tau_{n}$:

• •

Set $l=0$, $z_{n+\theta}^{l}=z_{n}$, and $\gamma_{n}^{l}=1$.

• •

For $l\geq 0$ solve till convergence the minimisation problem

[TABLE]

and set

[TABLE]

The minimisation problem (44) can be solved by the Gauss–Newton method (Nocedal & Wright, 2006) when applied to the EnKBF with data model (9). Hence only first order derivatives of the nonlinear forward operator $h$ are required. Alternatively, one can use a quasi-Newton method (Nocedal & Wright, 2006) directly on the nonlinear equations arising from the chosen time-step method. In case $A(z)$ is only positive semi-definite, the inverse in (44) should be replaced by the pseudo-inverse and the update is performed in the subspace spanned by the image of $A(z_{n+\theta}^{l})$.

In the subsequent numerical examples, we will compare the discrete gradient formulation (38) to the standard semi-implicit Euler method

[TABLE]

which can also be reformulated as a minimisation problem, that is,

[TABLE]

We note that (48) is structurally similar to the minimisation problem (44). However, the semi-implicit Euler method (47) requires the solution of such a minimisation problem only once per time-step while the discrete gradient method leads to a more elaborate implementation. Furthermore, iteration (47) does not obey a discrete decay property of the form (40).

While the discrete decay property (40) is sufficient to reach a critical point of (20) at $\tau\to\infty$, a sufficiently accurate approximation of the ensemble Kalman–Bucy filter equations is required in order to reproduce the posterior distribution. This statement is confirmed by the numerical experiments in the following section. We note in this context that the discrete gradient dynamics (38) is first order for all $\theta\in(0,1]$ except for $\theta=1/2$ when the method becomes second order.

4. Numerical results

We study the behaviour of the proposed gradient flow systems and their discretisation for a linear and a nonlinear one-dimensional inference problem, as well as for a partially observed nonlinear three-dimensional system. We start with the ensemble Kalman–Bucy filter formulation.

4.1. Ensemble Kalman–Bucy filter

4.1.1. Linear problem

We consider the scalar inference problem with Gaussian prior $X\sim{\rm N}(m_{0},\sigma_{0})$ and forward model

[TABLE]

The posterior distribution is also Gaussian with mean

[TABLE]

Kalman gain matrix $K=\sigma_{0}/(\sigma_{0}+r)$, and variance

[TABLE]

The experimental parameters are set to $m_{0}=1/2$, $\sigma_{0}=1$, $r=0.02$, $y=0.1$.

We implement the EnKBF formulation with $M=2$ ensemble members. Their initial positions $x_{0}^{i}$ are uniquely determined by the requirement that

[TABLE]

Furthermore, the analytic solutions $x_{\tau}^{i}$, $\tau\in(0,1]$, can be recovered from the associated solutions for the mean and variance, that is,

[TABLE]

with time-dependent Kalman gain

[TABLE]

The potential (7) reduces to

[TABLE]

with gradients

[TABLE]

The matrix $A(x^{1},x^{2})$ is given by

[TABLE]

We implement the discrete gradient method (38) with $\theta=1$ and the semi-implicit Euler (47) method for $\Delta\tau\in\{0.1,0.2,0.5,1\}$. The relative small value of $r$ renders (5) stiff and the explicit Euler method is unstable for any of the four chosen step-sizes $\Delta\tau$. The numerical results are compared to the analytic solutions as given by (53) and can be found in Figure 1. We find that the discrete gradient method systematically overestimates the variance while the semi-implicit Euler method has the opposite effect on the variance. The semi-implicit Euler method also performs remarkably well with respect to the mean while the discrete gradient method leads to relatively large errors in the final mean for step-sizes $\Delta\tau=0.5$ and $\Delta\tau=1.0$.

4.1.2. Nonlinear problem

We consider the following nonlinear forward map

[TABLE]

and a measurement error $\Xi\sim{\rm N}(0,r)$ with $r=1$. The observed value is $y=2$ and the prior distribution is Gaussian with $X\sim{\rm N}(-2,1/2)$. The posterior target distribution is given by

[TABLE]

The posterior mean and variance are given by $\overline{x}^{\ast}\approx 0.2095$ and $\sigma^{\ast}\approx 0.0211$, respectively. See Example 5.9 in Reich & Cotter (2015) for more details.

Given $M$ particles $x_{n}^{i}$ and a state vector $z_{n}=(x_{n}^{1},\ldots,x_{n}^{M})^{\rm T}$ the required solution $z_{n+\theta}^{l+1}\in\mathbb{R}^{M}$ to the minimisation problem (44) can be recast in the following form:

[TABLE]

where $f:\mathbb{R}^{M}\to\mathbb{R}^{2M+1}$ is defined by

[TABLE]

and the norm $\|\cdot\|_{C^{l}}$ by

[TABLE]

The Gauss–Newton method replaces the nonlinear minimisation problem (60) by the linearised problem

[TABLE]

Here the Jacobian matrix $Df(z)\in\mathbb{R}^{2M+1\times M}$ is provided by

[TABLE]

where $h^{\prime}(x)=7x^{2}/4-7x+8$,

[TABLE]

and $D(c)\in\mathbb{R}^{M\times M}$ is a diagonal matrix with the entries of $c$ on its diagonal.

The semi-implicit Euler method (47) leads to a rather similar Gauss–Newton iteration

[TABLE]

with $C^{l}$ in (60) being replaced by

[TABLE]

and $\theta=1$.

We compare the results from the semi-implicit Euler and discrete gradient method with $\theta=1$ to the gradient-free, explicit time-stepping method (IEnKF)

[TABLE]

which has been discussed, for example, by Amezcua et al. (2014); Blömker et al. (2018); de Wiljes et al. (2018). The discretisation (68) is also related to the DEnKF formulation of the ensemble Kalman as proposed by Sakov & Oke (2008).

We used $M=100$ particles in our numerical experiments and step-sizes of $\Delta\tau\in\{0.01,0.1,0.2,0.5\}$. An explicit Euler discretisation is stable only for the smallest of these step-sizes and is used to provide the reference solution. See Figure 2 for details.

We find that both the semi-implicit Euler and the discrete gradient method approximate the exact time-evolution of the gradient flow system, as represented by the explicit Euler method, rather well for step-sizes $\Delta\tau\in\{0.01,0.1\}$. The IEnKF method leads to systematically different results in the variance which, however, correspond rather well to the exact posterior value for this particular problem. We also note that the IEnKF method also performs well for the larger step-size of $\Delta\tau=0.2$. Overall the IEnKF method (68) emerges as the computationally cheapest and most accurate method among all the four tested time-stepping procedures.

4.2. Particle flow Fokker–Planck dynamics

We repeat the experiments from the previous subsection using particle flow Fokker–Planck dynamics. While its implementation is more involved, particle flow Fokker–Planck dynamics have the advantage of being more generally applicable than the Kalman-based formulations. We also implement the particle flow Fokker–Planck dynamics for a sequential data assimilation procedure for the Lorenz-63 model (Lorenz, 1963).

4.2.1. Linear problem

We use exactly the same experimental set up as already used for the EnKBF, that is, the measurement error is Gaussian with variance $r=0.02$ the the prior is Gaussian with mean $m_{0}=1/2$ and variance $\sigma_{0}=1$. We now use $M=10$ particles to approximate the (Gaussian) posterior distribution with a Gaussian mixture (33). The parameter $\alpha$ in (32) is set to $\alpha=0.005$. The particle system reaches fairly quickly a stationary configuration which approximates the first two moments of the posterior distribution rather well. See Figure 3. We also find that the discrete gradient method with $\theta=1$ relaxes to the stationary configuration within a single time-step even for large step-sizes $\Delta\tau$. The semi-implicit Euler method also performs rather well.

4.2.2. Nonlinear problem

We now come to a more challenging example with a non-Gaussian posterior distribution and consider the same nonlinear forward map and experimental setup as described in Section 4.1.2. Again, we use $M=100$ particles and consider the discrete gradient method with $\theta=1$ and the semi-implicit Euler method. The reference solution is also provided by the explicit Euler with small step size $\Delta\tau=2.5\times 10^{-6}$. Here we consider step-sizes of $\Delta\tau\in\{0.002,0.005,0.02,0.05\}$, see Figure 4 for results. The minimisation problem (44) is solved using the fminunc Matlab routine, which is based on a quasi-Newton method. The parameter $\alpha$ in (32) is set to $\alpha=0.01$.

As for the linear example, we find that the discrete gradient method with $\theta=1$ converges rapidly to the stationary configuration with both the mean and the variance of the posterior distribution well approximated. The semi-implicit method also behaves well with a only slightly slower convergence rate at large step-sizes.

Remark 4.1**.**

One could also implement the preconditioned formulation (37) with $\mathcal{A}(\{x_{n}^{j}\})$ being the empirical covariance matrix of the particle system times the ensemble size. This can imply significant computational savings for higher-dimensional problems if the kernel matrix $B$ in (30) is also provided by an empirical covariance matrix. Furthermore, the explicit computation of gradients can be avoided similar to what is done in the IEnKF implementation (68). More precisely, consider a negative log-likelihood function of the form (10), then

[TABLE]

See also Garbuno-Inigo et al. (2019) for related ideas in the context of Brownian dynamics. However, the gradient flow structure is lost under this approximation.

4.2.3. Lorenz-63

We use the experimental setting of Chustagulprom et al. (2016). Specifically, the Lorenz-63 differential equations (Lorenz, 1963) with the standard parameters $\sigma=10$, $\rho=28$, and $\beta=8/3$ are solved numerically with the implicit midpoint rule and step-size $\Delta t=0.01$. We only observe the first component of the state vector in observation intervals of $\Delta t_{\rm obs}=0.12$ with observation error variance $R=8$. A total of $K=200,000$ assimilation cycles are performed. The ensemble size varies between $M=15$ and $M=50$ and particle rejuvenation is applied with $\beta=0.2$ (Chustagulprom et al., 2016). The quality of the data assimilation method is assessed using the root mean squared error (RMSE)

[TABLE]

with $N_{x}=3$ and $t_{k}=0.12k$. Here $x_{\rm ref}(t_{k})\in\mathbb{R}^{3}$ denotes the reference solution and $\overline{x}_{\rm a}(t_{k})$ the analysis mean at observation time $t_{k}$.

At any given time in the sequential data assimilation procedure, a forecast ensemble $\{x_{\rm f}^{i}\}$ is generated by propagating the posterior sample (also referred to as analysis ensemble) $\{x_{\rm a}^{i}\}$ from the previous time step through the discretised Lorenz-63 equations. The transformation of this forecast ensemble to the posterior for the next time step via Gaussian mixture approximations and the particle flow Fokker–Planck dynamics is conducted as follows. One first computes the empirical mean $\overline{x}_{\rm f}$ and the empirical covariance matrix $P_{\rm f}$ from the forecast ensemble. Given a parameter value $\alpha\in(0,1]$, a Gaussian mixture approximation to the forecast distribution is provided by

[TABLE]

with covariance matrix $B_{\rm f}=(2\alpha-\alpha^{2})P_{\rm f}$ and

[TABLE]

as already discussed in Section 2.2. Given an observation model of the form (9) with linear forward map $h(x)=Hx$, the analysis (or posterior) distribution is also a Gaussian mixture

[TABLE]

with normalised weights

[TABLE]

new centers

[TABLE]

and new covariance matrix

[TABLE]

We next use the particle flow Fokker–Planck dynamics (20) with $\pi_{\ast}=\pi_{\rm a}$, kernel functions $\psi(x)={\rm n}(x;0,B_{\rm a})$ and initial particle positions $x_{0}^{i}=\overline{x}_{\rm a}^{i}$ in order to find an equally weighted Gaussian mixture approximation to (73). Let us denote the resulting ensemble at a suitably chosen final time $\tau_{\rm end}$ by $\{x_{\ast}^{i}\}$. Then the analysis ensemble is finally given by

[TABLE]

Note that $\alpha=1$ leads back to the standard ensemble Kalman filter (Evensen, 2006; Reich & Cotter, 2015), while one formally approaches a standard particle filter for $\alpha\to 0$. However, it should be noted that the particle flow Fokker–Planck dynamics is not well-defined for $\alpha\to 0$ and fixed ensemble size $M$.

The RMSEs for ensemble sizes $M\in\{15,\ldots,50\}$ and parameters $\alpha\in\{0.8,\ldots,1\}$ in (72) can be found in Table 1. One finds that $\alpha<1$ leads to some improvements especially for larger ensemble sizes. Overall these improvements are, however, not as impressive as those obtained by the hybrid methods considered by Chustagulprom et al. (2016).

5. Summary

We have discussed two instances of gradient dynamics arising from Bayesian inference. The gradient flow structure lends itself naturally to discrete gradient time-stepping method. However, such methods are implicit and potentially costly to implement. We found that discrete gradient dynamics does not improve the behaviour of iterated ensemble Kalman filters such as the derivative-free formulation (68). The situation is different for the particle flow Fokker–Planck dynamics, where the discrete gradient method was shown to converge rapidly. This difference in behaviour can be explained by the fact that ensemble Kalman filter requires an accurate representation of its dynamics while the Fokker–Planck dynamics only requires a rapid transition to an equilibrium solution. Applications of particle flow Fokker–Planck dynamics to a Lorenz-63 sequential data assimilation problem showed small improvements over a standard ensemble Kalman filter. Further applications of these methods may arise in the context of machine learning. See, for example, Ollivier (2018); Kovachki & Stuart (2018).

Acknowledgement. This research has been partially funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1294 ‘Data Assimilation’  (projects A02 and B04).