# Affine-invariant ensemble transform methods for logistic regression

Jakiw Pidstrigach and Sebastian Reich  Institut für Mathematik, Universität
Potsdam, Karl-Liebknecht-Str. 24/25, 14476 Potsdam

###### Abstract.

We investigate the application of ensemble transform approaches to Bayesian
inference of logistic regression problems. Our approach relies on appropriate
extensions of the popular ensemble Kalman filter and the feedback particle
filter to the cross entropy loss function and is based on a well-established
homotopy approach to Bayesian inference. The arising finite particle evolution
equations as well as their mean-field limits are affine-invariant.
Furthermore, the proposed methods can be implemented in a gradient-free manner
in case of nonlinear logistic regression and the data can be randomly
subsampled similar to mini-batching of stochastic gradient descent. We also
propose a closely related SDE-based sampling method which again is affine-
invariant and can easily be made gradient-free. Numerical examples demonstrate
the appropriateness of the proposed methodologies.

##  1\. Introduction

Statistical inference for logistic regression and classification problems is a
well studied problem from both parameter optimisation and Bayesian inference
perspectives [5]. While many classical optimisation and Markov chain Monte
Carlo (MCMC) methods are applicable, an efficient inference in the context of
semi-parametric models still poses computational challenges. With this paper,
we address this challenge by investigating the extension of coupling-of-
measures and ensemble transform methodologies [29] from the squared error loss
function to the cross entropy loss function typically used for logistic
regression and classification.

More specifically, previous work on ensemble methods for Bayesian inference,
such as the ensemble Kalman filter (EnKF) [8] and the feedback particle filter
(FPF) [39, 38, 34], have almost exclusively focused on the squared error loss
function of the form

(1) $\Psi_{\rm data}(\theta)=\frac{1}{2}(g(\theta)-t)^{\rm
T}\Gamma^{-1}(g(\theta)-t)\,,$

with $g:{\mathbb{R}}^{D}\to{\mathbb{R}}^{N}$ a forward map,
$t\in{\mathbb{R}}^{N}$ the data, $\Gamma\in{\mathbb{R}}^{N\times N}$ the
measurement error covariance matrix, and $\theta\in{\mathbb{R}}^{D}$ the
parameters to be estimated. Notable exceptions include the application of
ensemble Kalman inversion (EKI) [17] and the modified EnKF formulation of [13]
to the training of neural networks with a cross entropy loss function. While
these methods seek to minimise a regularised cross entropy loss function, the
more recent work [14] has also investigated ensemble-based sampling methods
for Bayesian inference in the context of logistic regression and
classification. This work relies on the time-stepping of appropriate
stochastic differential equations (SDEs) and is in line with several other
ensemble-based sampling methods such as the ensemble preconditioned MCMC
methods of [19], the ensemble Kalman sampler (EKS) [9], and affine invariant
Langevin dynamics (ALDI) [10]. Contrary to such invariance-of-measures based
SDE methodologies, the approach taken in this paper is instead founded on the
homotopy approach to Bayesian inference, as first formulated in a time-
continuous framework in [6, 26], and which is close in spirit to the iterative
application of the EnKF [31] and parameter estimation methods based on
sequential Monte Carlo (SMC) methods [4].

In addition to expanding the homotopy-based approach to logistic regression,
we address two further important concepts, namely affine-invariance of the
proposed methods [11, 19, 10] and sub-sampling of data points, as widely used
in stochastic gradient descent [20]. Both concepts are used to improve the
computational efficiency of optimisation and sampling methods. Take, for
example, a two dimensional Gaussian random variable with mean zero and
covariance matrix

$\Sigma=\left(\begin{array}[]{cc}\epsilon^{2}&0\\\ 0&1\end{array}\right).$

A random walk Metropolis–Hastings algorithm will sample inefficiently whenever
$\epsilon^{2}$ is vastly different from one. An affine-invariant modification,
on the other hand, will sample this problem as efficiently as if $\epsilon$
were set to $\epsilon=1$ [11]. Sub-sampling replaces the exact gradient of a
cost functional by a cheaper to evaluate stochastic approximation, which
agrees with the exact gradient in expectation.

We also discuss the possibility for derivative-free implementations [8, 17, 9,
10], localisation [8, 29] via dropouts [32], and efficient linearly implicit
time-stepping methods [2].

Our main contributions with regard to Bayesian homotopy methods for logistic
regression are:

  * •

an extension of affine-invariant ensemble transform approaches, such as the
EnKF, to logistic regression,

  * •

an affine-invariant generalisation of the FPF and its application to logistic
regression,

  * •

an extension of data sub-sampling (mini-batches) to Bayesian homotopy methods.

In a further step, we combine these homotopy approaches with SDE-based
sampling methods in order to derive

  * •

a derivative-free and affine-invariant SDE-based sampling methods for logistic
regression.

We demonstrate the appropriateness of the proposed methods by means of a set
of numerical experiments. Extensions to nonlinear and multi-class logistic
regression [5] are also discussed. We also briefly discuss the application to
sigmoidal Cox processes [1].

The layout of this paper is as follows. The required mathematical background
material on both logistic regression is collected in Section 2. The homotopy
approach to Bayesian inference is summarised in Section 3. There we also
present an affine-invariant formulation of the homotopy approach and discuss
and analyse data sub-sampling in the spirit of [20]. Section 4 develops three
different algorithmic approaches for the implementation of homotopy-based
Bayesian inference for logistic regression. More specifically, we propose an
affine-invariant modification of the FPF and two extensions of the EnKF to
logistic regression. We also discuss robust numerical implementations
combining dropouts [32] with localisation [8, 29] and linearly implicit time-
stepping methods [2]. Section 5 combines SDE-based sampling methods with a
homotopy-based drift term in order to derive a gradient-free and affine-
invariant algorithm for Bayesian inference. While the proposed method is only
exact for Gaussian measures, it can be used for faster equilibration and
approximate inference provided the posterior distribution is close to
Gaussian. Numerical results are presented and evaluated in Section 6 and
provide a proof-of-concept while a more detailed exploration is left to
follow-up work. Several possible extensions of the proposed methods, namely
multi-class classification, nonlinear logistic regression, and sigmoidal Cox
processes, are discussed in Section 7 followed by a summary statement in
Section 8.

##  2\. Mathematical problem formulation

This paper is motivated by the classical logistic regression problem arising
from classification into $L>1$ classes, denoted by $C_{l}$, over an input
space $x\in{\mathbb{R}}^{J}$ [5]. In order to simplify the exposition, we
start with the binary classification problem, that is, $L=2$ and focus on the
linear case. The extension to nonlinear and multi-class logistic regression is
discussed in Section 7.

The posterior probability for class $C_{1}$ is defined as a logistic sigmoid

$\sigma(a)=\frac{1}{1+\exp(-a)}$

acting on a linear combination of $x$-dependent and vector-valued features
$\phi_{x}\in{\mathbb{R}}^{D}$ so that

(2) $\pi(C_{1}|\phi_{x})=\sigma\left(\theta^{\rm T}\phi_{x}\right).$

The model has $D$ adjustable parameters $\theta\in{\mathbb{R}}^{D}$. The
probability of the complementary class $C_{2}$ is given by
$\pi(C_{2}|\phi_{x})=1-\pi(C_{1}|\phi_{x})$. Given a data set
$\\{(x_{n},t_{n})\\}_{n=1}^{N}$ with labels $t_{n}\in\\{0,1\\}$, the negative
log-likelihood function is given by the cross-entropy function

(3) $\Psi_{\rm data}(\theta)=-\sum_{n=1}^{N}\left\\{t_{n}\ln
y_{n}+(1-t_{n})\ln(1-y_{n})\right\\},$

where $y_{n}:=y_{x_{n}}(\theta)$ with

(4) $y_{x}(\theta):=\sigma\left(\theta^{\rm T}\phi_{x}\right).$

Here $t_{n}=1$ for samples $x_{n}$ which are assigned to class $C_{1}$ and
$t_{n}=0$ for samples from the complementary class $C_{2}$.

We introduce the further shorthand

$\phi_{n}=\phi_{x_{n}}$

and, using

$\frac{\rm d}{{\rm d}a}\sigma(a)=\sigma(a)(1-\sigma(a)),$

obtain the gradient of the cross-entropy function

(5) $\nabla_{\theta}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}(y_{n}-t_{n})\phi_{n}\,.$

We rewrite this gradient more compactly as as

$\nabla_{\theta}\Psi_{\rm data}(\theta)=\Phi\,(y(\theta)-t)$

using

$y(\theta)=(y_{1},\ldots,y_{n})^{\rm T}\in{\mathbb{R}}^{N},\quad
t=(t_{1},\ldots,t_{N})^{\rm
T}\in{\mathbb{R}}^{N},\quad\Phi=\left(\phi_{1},\ldots,\phi_{N}\right)\in{\mathbb{R}}^{D\times
N}.$

The term $y(\theta)-t$ is often called the innovation in the Kalman filter
literature. The Hessian matrix of second-order derivative is given by

$D^{2}_{\theta}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}\phi_{n}y_{n}(1-y_{n})\phi_{n}^{\rm T}=\Phi
R(\theta)\Phi^{\rm T}$

where $R(\theta)\in{\mathbb{R}}^{N\times N}$ denotes the diagonal matrix with
entries $r_{nn}=y_{n}(1-y_{n})$.

Let us denote a minimiser of $\Psi_{\rm data}(\theta)$ by $\theta_{\rm MLE}$.
We note that the maximum likelihood estimators (MLEs) are not uniquely
determined in the over-parametrised setting, that is, when $D\gg N$. For that
reason and in order to avoid overfitting when $D\approx N$, we introduce a
Gaussian prior probability density function (PDF) over the parameters via

$\pi_{\rm prior}(\theta)\propto\exp(-\Psi_{\rm prior}(\theta)),\qquad\Psi_{\rm
prior}(\theta):=\frac{1}{2}(\theta-m_{\rm prior})^{\rm T}\Sigma_{\rm
prior}^{-1}(\theta-m_{\rm prior}),$

that is, the posterior parameter distribution $\pi_{\rm post}$ is given by

(6) $\pi_{\rm post}(\theta)\propto\exp(-\Psi_{\rm
post}(\theta)),\qquad\Psi_{\rm post}(\theta):=\Psi_{\rm
data}(\theta)+\Psi_{\rm prior}(\theta)\,.$

Here the prior mean $m_{\rm prior}\in{\mathbb{R}}^{D}$ and covariance matrix
$\Sigma_{\rm prior}\in{\mathbb{R}}^{D\times D}$ need to be chosen
appropriately.

It is known that $\Psi_{\rm post}$ is a strongly convex function since
$D^{2}_{\theta}\Psi_{\rm data}(\theta)\geq 0$ and, hence, the MAP estimator,
that is,

$\theta_{\rm MAP}=\arg\min_{\rm\theta\in{\mathbb{R}}^{D}}\Psi_{\rm
post}(\theta),$

is uniquely defined. Furthermore, assuming the availability of sufficiently
many data points, that is $N\gg 1$, the posterior is well approximated by a
Gaussian with mean $\theta_{\rm MAP}$ and covariance matrix

$\Sigma_{\rm MAP}=\left(D^{2}_{\theta}\Psi_{\rm post}(\theta_{\rm
MAP})\right)^{-1},$

where

$D^{2}_{\theta}\Psi_{\rm post}(\theta)=\Phi R(\theta)\Phi^{\rm T}+\Sigma_{\rm
prior}^{-1}\in{\mathbb{R}}^{D\times D}$

denotes the Hessian matrix of second-order derivatives. This is an implication
of the Bernstein–von Mises theorem [33].

##  3\. The homotopy approach to Bayesian inference

While there are many Monte Carlo methods available for sampling from the
posterior $\pi_{\rm post}$, we focus in this paper on ensemble transform
methods based on the homotopy approach [6, 26]

(7) $\pi_{\tau}(\theta)\propto\exp(-\tau\Psi_{\rm data}(\theta))\,\pi_{\rm
prior}(\theta)$

with $\tau\in[0,1]$, that is, $\pi_{0}=\pi_{\rm prior}$ and $\pi_{\rm
post}=\pi_{1}$. The PDFs $\pi_{\tau}$ satisfy the evolution equation

(8) $\partial_{\tau}\pi_{\tau}=-(\Psi_{\rm data}-\pi_{\tau}\left[\Psi_{\rm
data}\right])\,\pi_{\tau},$

where $\pi_{\tau}[f]$ denotes the expectation value of a function
$f:{\mathbb{R}}^{D}\to{\mathbb{R}}$ under the PDF $\pi_{\tau}$.

In order to derive the associated evolution equations for random variables
$\theta_{\tau}\sim\pi_{\tau}$, we introduce potentials
$V_{\tau}:{\mathbb{R}}^{D}\to{\mathbb{R}}$, which are defined by the elliptic
partial differential equation (PDE)

(9)
$\nabla_{\theta}\cdot(\pi_{\tau}\Sigma_{\tau}\nabla_{\theta}V_{\tau})=-(\Psi_{\rm
data}-\pi_{\tau}\left[\Psi_{\rm data}\right])\,\pi_{\tau},$

where $\Sigma_{\tau}$ denotes the covariance matrix of $\pi_{\tau}$, that is,

$\Sigma_{\tau}=\pi_{\tau}\left[(\theta-\pi_{\tau}[\theta])(\theta-\pi_{\tau}[\theta])^{\rm
T}\right].$

Hence the evolution equations are given by

(10) $\frac{\rm d}{{\rm
d}\tau}\theta_{\tau}=-\Sigma_{\tau}\nabla_{\tau}V_{\tau}(\theta_{\tau}),\qquad\theta_{0}\sim\pi_{\rm
prior},$

such that

(11) $\theta_{1}\sim\pi_{\rm post}.$

See [26] for further details.

######  Remark 1.

If different models $\mathcal{M}_{i}$, $i=1,\ldots,I_{M}$, need to be compared
via their Bayes factors $\mathbb{P}[\mathcal{M}_{i}]$ [16], then the homotopy
approach (8) can also be used and gives rise to the evolution equation

$\frac{\rm d}{{\rm d}\tau}\mathbb{P}_{\tau}=-\pi_{\tau}\left[\Psi_{\rm
data}\right],$

$\tau\in[0,1]$, which is to be solved for each model $\mathcal{M}_{i}$ with
$\mathbb{P}_{0}$ its prior probability, $\pi_{\tau}$ its transformed parameter
PDFs, as defined by (8), and $\Psi_{\rm data}$ the model’s negative log-
likelihood function. It then holds that
$\mathbb{P}[\mathcal{M}_{i}]=\mathbb{P}_{1}$.

######  Remark 2.

It is, of course, feasible to use a stopping time different from $\tau=1$ in
the homotopy approach (10). This has already been hinted at in in the context
of optimisation and Tikhonov regularisations in [26]; while subsequently being
fully explored by the EKI methodology. See [17, 15] for a detailed description
in the context of machine learning as well as for an extensive literature
survey.

###  3.1. Affine-invariance

Note that any vector field $F_{\tau}(\theta)$, which satisfies

$\nabla_{\theta}\cdot(\pi_{\tau}F_{\tau})=-(\Psi_{\rm
data}-\pi_{\tau}\left[\Psi_{\rm data}\right])\,\pi_{\tau},$

will lead to an associated evolution equation

(12) $\frac{\rm d}{{\rm
d}\tau}\theta_{\tau}=-F_{\tau}(\theta_{\tau}),\qquad\theta_{0}\sim\pi_{\rm
prior}$

such that (11) holds. The specific choice $F_{\tau}=\Sigma_{\tau}\nabla
V_{\tau}$ in (9) is motivated by the concept of affine-invariance [11, 10],
which we explain next.

Let

(13) $\theta=A\overline{\theta}+b$

denote an invertible affine transformation between two variables
$\theta\in{\mathbb{R}}^{D}$ and $\overline{\theta}\in{\mathbb{R}}^{D}$.
Furthermore, the Bayesian inference problem for the variable $\theta$ is
defined by $\Psi_{\rm data}$ and $\Psi_{\rm prior}$ and the corresponding
problem for the transformed variable $\overline{\theta}$ by

$\overline{\Psi}_{\rm data}(\overline{\theta}):=\Psi_{\rm
data}(A\overline{\theta}+b),\qquad\overline{\Psi}_{\rm
prior}(\overline{\theta}):=\Psi_{\rm prior}(A\overline{\theta}+b).$

Hence the PDFs $\overline{\pi}_{\tau}$ for $\overline{\theta}$ satisfy

$\overline{\pi}_{\tau}(\overline{\theta})=\pi_{\tau}(A\overline{\theta}+b)\,|A|,\qquad\tau\in[0,1]\,.$

######  Definition 1.

A homotopy formulation (12) is called affine-invariant if the associated
vector fields $F_{\tau}$ and $\overline{F}_{\tau}$ satisfy

(14) $F_{\tau}(A\overline{\theta}+b)=A\overline{F}_{\tau}(\overline{\theta})$

for any suitable choice of $A\in{\mathbb{R}}^{D\times D}$ and
$b\in{\mathbb{R}}^{D}$, which implies
$\theta_{\tau}=A\overline{\theta}_{\tau}+b$ for all $\tau>0$ if it holds at
$\tau=0$.

######  Lemma 1.

The choice

(15) $F_{\tau}(\theta)=\Sigma_{\tau}\nabla_{\theta}V_{\tau}(\theta)$

leads to an affine-invariant homotopy formulation and the transformed vector
field is given by

$\overline{F}_{\tau}(\overline{\theta})=\overline{\Sigma}_{\tau}\nabla_{\overline{\theta}}\overline{V}_{\tau}(\overline{\theta})$

with the potential $\overline{V}_{\tau}$ satisfying the elliptic PDE

(16)
$\nabla_{\overline{\theta}}\cdot(\overline{\pi}_{\tau}\overline{\Sigma}_{\tau}\nabla_{\overline{\theta}}\overline{V}_{\tau})=-(\overline{\Psi}_{\rm
data}-\overline{\pi}_{\tau}\left[\overline{\Psi}_{\rm
data}\right])\,\overline{\pi}_{\tau}.$

###### Proof.

Since the covariance matrices satisfy

$\Sigma_{\tau}=A\overline{\Sigma}_{\tau}A^{\rm T},$

the PDE (9) transforms into the PDE (16) under the affine transformation (13).
Hence we have

$\overline{V}_{\tau}(\overline{\theta})=V_{\tau}(A\overline{\theta}+b)$

and obtain

$\displaystyle\Sigma_{\tau}\nabla_{\theta}V_{\tau}(\theta)$
$\displaystyle=A\overline{\Sigma}_{\tau}A^{\rm
T}\nabla_{\theta}V_{\tau}(A\overline{\theta}+b)$
$\displaystyle=A\overline{\Sigma}_{\tau}\nabla_{\overline{\theta}}V_{\tau}(A\overline{\theta}+b)$
$\displaystyle=A\overline{\Sigma}_{\tau}\nabla_{\overline{\theta}}\overline{V}_{\tau}(\overline{\theta})\,,$

which verifies (14). ∎

We note that the affine invariance of (10) is maintained under the forward
Euler time discretisation

$\theta_{\tau_{k+1}}=\theta_{\tau_{k}}-\Delta\tau\Sigma_{\tau_{k}}\nabla_{\theta}V_{\tau_{k}}(\theta_{\tau_{k}})$

with step-size $\Delta\tau>0$ and $\tau_{k}=k\Delta\tau$, $k=0,\ldots,K-1$
such that $K\Delta\tau=1$. This immediately follows from property (14).

######  Remark 3.

The choice of the vector field $F_{\tau}(\theta)$ in (15) is not unique. It
is, for example, possible to add any drift term
$\Sigma_{\tau}\nabla_{\theta}U_{\tau}(\theta)$ with the potential $U_{\tau}$
chosen such that

(17) $\nabla_{\theta}\cdot(\pi_{\tau}\Sigma_{\tau}\nabla_{\theta}U_{\tau})=0.$

A possible choice is

$U_{\tau}(\theta)=\log\pi_{\tau}+\tau\Psi_{\rm data}(\theta)-\log\pi_{\rm
prior}(\theta)$

and the associated operator on the left hand side of (17) corresponds to the
nonlinear Fokker–Planck equation arising from the EKS/ALDI mean-field
evolution equations [9, 10] for sampling from the PDF $\pi_{\tau}$. Because of
(7), $U_{\tau}(\theta)=\mbox{const.}$ and (17) follows trivially. Upon
formally replacing the drift term
$\Sigma_{\tau}\nabla_{\theta}U_{\tau}(\theta)$ by its EKS/ALDI stochastic
mean-field representation, (10) turns into the affine-invariant stochastic
evolution equation

(18) ${\rm
d}\theta_{\tau}=-\Sigma_{\tau}\nabla_{\theta}\left\\{V_{\tau}(\theta)+\tau\nabla_{\theta}\Psi_{\rm
data}(\theta)-\log\pi_{\rm prior}(\theta)\right\\}{\rm
d}t+\Sigma^{1/2}_{\tau}{\rm d}W_{\tau},$

where $W_{\tau}$ denotes standard $D$-dimensional Brownian motion. The
associated marginal PDF $\pi_{\tau}$ still satisfies (8) provided
$\pi_{0}=\pi_{\rm prior}$ at initial time $\tau=0$.

###  3.2. Data sub-sampling

If the data sets are large, then it makes sense to randomly sub-sample the
data and to only use these subsets in the computation of the gradient (5).
More precisely, let $\Delta$ denote a randomly selected subset of
$\\{1,\ldots,N\\}$ without replacement of cardinality $N^{\prime}\leq N$.
Given such a mini-batch $\Delta$, the full gradient (5) is replaced by its
approximation

$\nabla_{\theta}\Psi_{\rm
data}(\theta)\approx\frac{N}{N^{\prime}}\sum_{n\in\Delta}(y_{n}-t_{n})\phi_{n}\,.$

More abstractly, data sub-sampling in terms of mini-batches gives rise to
random negative log-likelihood functions $\Psi_{\rm data}^{\gamma}$, where
$\gamma$ characterises the chosen mini-batch $\Delta^{\gamma}$ and is the
realisation of an appropriate random variable $\Gamma$ such that

(19) $\mathbb{E}_{\Gamma}[\Psi_{\rm data}^{\gamma}(\theta)]=\Psi_{\rm
data}(\theta)$

for all parameter values $\theta\in\mathbb{R}^{D}$. See [20] for a detailed
investigation of data sub-sampling in the context of stochastic gradient
descent and related methods. If $\Psi_{\rm data}^{\gamma}$ is now being used
in (9) instead of $\Psi_{\rm data}$, then we denote the resulting potential by
$V_{\tau}^{\gamma}$, which is also unbiased, that is,

$\mathbb{E}_{\Gamma}[\nabla_{\theta}V_{\tau}^{\gamma}(\theta)]=\nabla_{\theta}V_{\tau}(\theta)$

due to (19). Let us now time-step the implied mean-field ODEs (10) by the
forward Euler method, that is,

(20)
$\theta_{\tau_{k+1}}=\theta_{\tau_{k}}-\Delta\tau\Sigma_{\tau_{k}}\nabla_{\theta}V_{\tau_{k}}^{\gamma_{k}}(\theta_{\tau_{k}})\,,$

where the $\gamma_{k}$’s are independent realisations of the random variable
$\Gamma$. A formal application of the modified equation analysis of [20] leads
to the stochastic differential equation (SDE)

(21) ${\rm
d}\theta_{\tau}=-\Sigma_{\tau}\left\\{\nabla_{\theta}V_{\tau}(\theta_{\tau})\,{\rm
d}\tau-(\Delta\tau\Omega_{\tau}(\theta_{\tau}))^{1/2}{\rm
d}W_{\tau}\right\\}\,,$

where $W_{\tau}$ denotes standard $D$-dimensional Brownian motion and

$\Omega_{\tau}(\theta):=\mbox{cov}_{\Gamma}\left[\nabla_{\theta}V_{\tau}^{\gamma}(\theta)\right]\,.$

Here $\mbox{cov}_{\Gamma}[f^{\gamma}]$ denotes the covariance matrix of a
vector-valued random function $f^{\gamma}\in{\mathbb{R}}^{D}$ with respect to
the random variable $\Gamma$. We conclude from (21) that the approximation
error inflicted on the posterior distribution at $\tau=1$ can be made
arbitrarily small by reducing the time-step $\Delta\tau$ of the forward Euler
method (20). While a rigorous mathematical analysis is left to future work, we
discuss an application to Cox point processes in Section 7.3 below and analyse
a Gaussian approximation in the Appendix in more detail.

##  4\. Ensemble transform algorithms

We now discuss several Monte Carlo implementations of the proposed affine-
invariant homotopy approach to Bayesian inference for logistic regression. All
methods start with a set of $M$ samples $\theta_{0}^{(i)}$, $i=1,\ldots,M$,
from the prior distribution $\pi_{0}=\pi_{\rm prior}$. The three methods
considered here differ in the evolution equations for the ensemble
$\\{\theta_{\tau}^{(i)}\\}$. We introduce the empirical mean

$\hat{m}_{\tau}=\frac{1}{M}\sum_{i=1}^{M}\theta_{\tau}^{(i)}$

and the empirical covariance matrix

$\hat{\Sigma}_{\tau}=\frac{1}{M-1}\sum_{i=1}^{M}(\theta_{\tau}^{(i)}-\hat{m}_{\tau})(\theta_{\tau}^{(i)}-\hat{m}_{\tau})^{\rm
T}$

for $\tau\geq 0$. We also introduce the matrix
$\Theta_{\tau}\in{\mathbb{R}}^{D\times M}$ of ensemble deviations

(22)
$\Theta_{\tau}=\left(\theta_{\tau}^{(1)}-\hat{m}_{\tau},\theta_{\tau}^{(2)}-\hat{m}_{\tau},\ldots,\theta_{\tau}^{(M)}-\hat{m}_{\tau}\right),$

which leads to the compact representation

(23) $\hat{\Sigma}_{\tau}=\frac{1}{M-1}\Theta_{\tau}\Theta_{\tau}^{\rm T}$

of the empirical covariance matrix. The empirical measure associated with the
ensemble $\\{\theta_{\tau}^{(i)}\\}$ at time $\tau$ is denoted by
$\hat{\pi}_{\tau}$ implying, for example,

$\hat{m}_{\tau}=\hat{\pi}_{\tau}[\theta].$

More generally,

$\hat{\pi}_{\tau}[f]=\frac{1}{M}\sum_{i=1}^{M}f(\theta_{\tau}^{(i)})$

for an observable $f(\theta)$. We describe an affine-invariant extension of
the FPF next. The key building block is a diffusion map approximation to the
PDE (9).

###  4.1. Diffusion maps and the FPF

The right-hand side of the elliptic PDE (9) leads to the elliptic operator

$\Delta_{\pi}\phi:=\frac{1}{\pi}\nabla_{\theta}\cdot(\pi\Sigma\nabla_{\theta}\phi)$

with covariance matrix

$\Sigma=\pi\left[(\theta-\pi[\theta])(\theta-\pi[\theta])^{\rm T}\right]$

and (9) can be rewritten as

$-\Delta_{\pi_{\tau}}V_{\tau}=\Psi_{\rm data}-\pi_{\tau}\left[\Psi_{\rm
data}\right].$

Following the FPF methodology, as carefully described in [35], we regularise
the elliptic PDE problem (9) through the fixed-point equation

$\tilde{V}_{\tau}=P_{\epsilon}\tilde{V}_{\tau}+\epsilon(\Psi_{\rm
data}-\pi_{\tau}\left[\Psi_{\rm data}\right])$

with regularisation parameter $\epsilon>0$ and semigroup
$P_{\epsilon}=e^{\epsilon\Delta_{\pi_{\tau}}}$. One then replaces $V_{\tau}$
by $\tilde{V}_{\tau}$ in the evolution equation (10).

We note that $P_{\epsilon}$ is the semi-group associated with the SDE

(24) ${\rm
d}S_{\epsilon}=\Sigma_{\tau}\nabla_{\theta}\log\pi_{\tau}(S_{\epsilon})\,{\rm
d}\epsilon+\sqrt{2}\Sigma_{\tau}^{1/2}{\rm d}W_{\epsilon},$

where $W_{\epsilon}$, $\epsilon\geq 0$, is standard Brownian motion in
${\mathbb{R}}^{D}$. The SDE (24) is affine invariant and is closely related to
the EKS/ALDI mean-field equations [9, 10].

The diffusion map approximation $T_{\epsilon}$ of $P_{\epsilon}$ is defined as
follows:

$T_{\epsilon}f(\theta)=\frac{1}{n_{\epsilon}(\theta)}\int_{R^{D}}k_{\epsilon}(\theta,\theta^{\prime})f(\theta^{\prime})\pi_{\tau}(\theta^{\prime}){\rm
d}\theta^{\prime}\,,$

where $n_{\epsilon}(\theta):=\int
k_{\epsilon}(\theta,\theta^{\prime})\pi_{\tau}(\theta^{\prime}){\rm
d}\theta^{\prime}$ is the normalisation constant,

$k_{\epsilon}(\theta,\theta^{\prime}):=\frac{g_{\epsilon}(\theta,\theta^{\prime})}{\sqrt{\int
g_{\epsilon}(\theta,\theta^{\prime\prime})\pi_{\tau}(\theta^{\prime\prime}){\rm
d}\theta^{\prime\prime}}\sqrt{\int
g_{\epsilon}(\theta^{\prime},\theta^{\prime\prime})\pi_{\tau}(\theta^{\prime\prime}){\rm
d}\theta^{\prime\prime}}}$

and

(25)
$g_{\epsilon}(\theta,\theta^{\prime}):=\exp\left(-\frac{1}{4\epsilon}(\theta-\theta^{\prime})^{\rm
T}\Sigma_{\tau}^{-1}(\theta-\theta^{\prime})\right)$

is the Gaussian kernel associated with the pure diffusion in (24).

The next step is to replace $\pi_{\tau}$ by the ensemble-induced empirical
measure $\hat{\pi}_{\tau}$, which results in the empirical approximation

$T_{\epsilon}^{(M)}f(\theta):=\frac{1}{n_{\epsilon}^{(M)}(\theta)}\sum_{i=1}^{M}k_{\epsilon}^{(M)}(\theta,\theta^{(i)}_{\tau})f(\theta^{(i)}_{\tau})\,,$

where
$n_{\epsilon}^{(M)}(\theta):=\sum_{i}k_{\epsilon}^{(M)}(\theta,\theta^{(i)}_{\tau})$
is the normalisation constant, and

$k_{\epsilon}^{(M)}(\theta,\theta^{\prime}):=\frac{g^{(M)}_{\epsilon}(\theta,\theta^{\prime})}{\sqrt{\sum_{i}g^{(M)}_{\epsilon}(\theta,\theta^{(i)}_{\tau})}\sqrt{\sum_{i}g^{(M)}_{\epsilon}(\theta^{\prime},\theta^{(i)}_{\tau})}}\,.$

Here $g^{(M)}_{\epsilon}$ denotes the Gaussian kernel (25) with the covariance
matrix $\Sigma_{\tau}$ replaced by its empirical estimate
$\hat{\Sigma}_{\tau}$. One finally introduces the $M\times M$ Markov matrix
${\rm T}$ with entries

(26) ${\rm
T}_{ij}=\frac{1}{n_{\epsilon}^{(M)}(\theta_{\tau}^{(i)})}k_{\epsilon}^{(M)}(\theta_{\tau}^{(i)},\theta_{\tau}^{(j)})$

and the finite-dimensional fixed-point equation

$\tilde{\rm V}={\rm T}\tilde{\rm V}+\epsilon\,\Delta{\rm\Psi}_{\rm data}$

with

$\Delta{\rm\Psi}_{\rm data}:=\left(\Psi_{\rm
data}(\theta_{\tau}^{(1)})-\hat{\pi}_{\tau}[\Psi_{\rm data}],\ldots,\Psi_{\rm
data}(\theta_{\tau}^{(M)})-\hat{\pi}_{\tau}[\Psi_{\rm data}]\right)^{\rm T}.$

Given the solution vector $\tilde{\rm V}\in{\mathbb{R}}^{M}$ with entries
$\tilde{\rm V}^{(j)}$, $j=1,\ldots,M$, the approximation to the drift term in
the ODE (10) is now provided as follows [35]:

$\hat{\Sigma}_{\tau}\nabla_{\theta}\tilde{V}_{\tau}(\theta):=\hat{\Sigma}_{\tau}\nabla_{\theta}\left\\{\frac{1}{n_{\epsilon}^{(M)}(\theta)}\sum_{j=1}^{M}k_{\epsilon}^{(M)}(\theta,\theta^{(j)}_{\tau})\left(\tilde{\rm
V}^{(j)}+\epsilon\,\Delta\Psi_{\rm data}^{(j)}\right)\right\\}.$

######  Lemma 2.

The affine-invariant FPF drift term for the ensemble members
$\theta_{\tau}^{(i)}$, $i=1,\ldots,M$, is given by

(27)
$\hat{\Sigma}_{\tau}\nabla_{\theta}\tilde{V}_{\tau}(\theta_{\tau}^{(i)})=\sum_{j=1}^{M}s_{ij}\theta_{\tau}^{(j)}$

with coefficients

$s_{ij}=\frac{1}{2\epsilon}{\rm T}_{ij}\left(r_{j}-\sum_{k=1}^{M}{\rm
T}_{ik}r_{k}\right),\qquad r_{j}:=\tilde{\rm
V}^{(j)}+\epsilon\,\Delta\Psi_{\rm data}^{(j)}\,.$

###### Proof.

The calculations leading to (27) are identical to those for the standard FPF
[35] taking note that the covariance matrix $\hat{\Sigma}_{\tau}$ cancels out.
∎

The resulting interacting particle system can be propagated in time using the
forward Euler method with step size $\Delta\tau=1/K$, that is,

(28)
$\theta_{\tau_{k+1}}^{(i)}=\theta_{\tau_{k}}^{(i)}-\Delta\tau\hat{\Sigma}_{\tau_{k}}\nabla_{\theta}\tilde{V}_{\tau_{k}}(\theta_{\tau_{k}}^{(i)}),\qquad
i=1,\ldots,M,$

$k=0,\ldots,K-1$, with the right-hand side defined by (27). Also note that
$M\to\infty$ leads formally to the mean-field forward Euler method (20).

######  Remark 4.

We note that the Markov matrix ${\rm T}$, defined by (26), is non-reversible
contrary to the underlying stochastic process (24). This problem can be
addressed by a bi-stochastic approximation, which also leads to improved
convergence rates. See [37] for more details.

While the FPF leads to a consistent approximation in the limit $M\to\infty$
and $\epsilon\to 0$, its application is limited to low-dimensional problems.
We next describe two approximations that extend to high-dimensional inference
problems.

###  4.2. Second-order methods

One can derive evolution equations for the mean $m_{\tau}$ and the covariance
matrix $\Sigma_{\tau}$ of $\pi_{\tau}$ from equation (8). Furthermore, if one
assumes that $\pi_{\tau}$ is well approximated by a Gaussian PDF with mean
$m_{\tau}$ and covariance matrix $\Sigma_{\tau}$, the following evolution
equations for the mean and the covariance matrix arise.

######  Lemma 3.

If the temporal PDF $\pi_{t}$ in (9) is Gaussian with mean $m_{\tau}$ and
covariance matrix $\Sigma_{\tau}$, then their evolution equations are given by

$\frac{\rm d}{{\rm
d}\tau}m_{\tau}=-\Sigma_{\tau}\,\pi_{\tau}[\nabla_{\theta}\Psi_{\rm
data}]=-\Sigma_{\tau}\Phi(\pi_{\tau}[y]-t),$

and

$\frac{\rm d}{{\rm
d}\tau}\Sigma_{\tau}=-\Sigma_{\tau}\,\pi_{\tau}[D^{2}_{\theta}\Psi_{\rm
data}]\Sigma_{\tau}=-\Sigma_{\tau}\Phi\,\pi_{\tau}[R]\Phi^{\rm
T}\,\Sigma_{\tau},$

respectively.

###### Proof.

The evolution equations are derived using the following well-known identity
for Gaussian PDFs $\pi_{\tau}$ (see, for example, [23]):

(29)
$\pi_{\tau}[\theta\,(f(\theta)-\pi_{\tau}[f])]=\Sigma_{\tau}\pi_{\tau}[\nabla_{\theta}f]$

as well as

(30) $\pi_{\tau}[(\theta-\pi_{\tau}[\theta])(\theta-\pi_{\tau}[\theta])^{\rm
T}(f(\theta)-\pi_{\tau}[f])]=\Sigma_{\tau}\pi_{\tau}[D^{2}_{\theta}f]\Sigma_{\tau}$

for any scalar-valued smooth function $f:{\mathbb{R}}^{D}\to{\mathbb{R}}$. ∎

In order to derive corresponding equations for the ensemble
$\\{\theta_{\tau}^{(i)}\\}$, we introduce

$\hat{y}_{\tau}:=\hat{\pi}_{\tau}[y]=\frac{1}{M}\sum_{j=1}^{M}y_{\tau}^{(j)},\qquad
y_{\tau}^{(j)}:=y\left(\theta_{\tau}^{(j)}\right)\in{\mathbb{R}}^{N}\,,$

as well as

$\hat{R}_{\tau}:=\hat{\pi}_{\tau}[R]=\frac{1}{M}\sum_{j=1}^{M}R_{\tau}^{(j)},$

where the $R_{\tau}^{(j)}\in{\mathbb{R}}^{N\times N}$ are diagonal matrices
with entries $r_{nn}^{(j)}=y_{\tau,n}^{(j)}\left(1-y_{\tau,n}^{(j)}\right)$
for $n=1,\ldots,N$ with

$y_{\tau,n}^{(j)}=\sigma\left((\theta_{\tau}^{(j)})^{\rm
T}\phi_{x_{n}}\right)\in{\mathbb{R}}\,.$

We finally obtain the evolution equations

(31) $\frac{\rm d}{{\rm
d}\tau}\hat{m}_{\tau}=-\hat{\Sigma}_{\tau}\Phi(\hat{y}_{\tau}-t)$

for the ensemble mean and

(32) $\frac{\rm d}{{\rm
d}\tau}\Theta_{\tau}=-\frac{1}{2}\hat{\Sigma}_{\tau}\Phi\hat{R}_{\tau}\Phi^{\rm
T}\Theta_{\tau}$

for the ensemble deviations (22). Note that these equations can be used to
propagate the ensemble $\\{\theta_{\tau}^{(i)}\\}$ regardless of the Gaussian
assumption made for their derivation. In the following section, we will
introduce further approximations which allow us to introduce an affine-
invariant gradient flow structure.

###  4.3. Ensemble Kalman–Bucy filter

We now turn to a formulation with gradient flow structure in the spirit of the
ensemble Kalman–Bucy filter (EnKBF) for quadratic loss functions (1) [3, 26,
28]. Provided that

$\hat{y}_{\tau}=\hat{\pi}_{\tau}[y]\approx
y(\hat{m}_{\tau})=y(\hat{\pi}_{\tau}[\theta])$

and

$y_{\tau}^{(i)}=y\left(\theta_{\tau}^{(i)}\right)\approx\hat{y}_{\tau}+\hat{R}_{\tau}\Phi^{\rm
T}(\theta_{\tau}^{(i)}-\hat{m}_{\tau})\,,$

one formally obtains the approximation

(33) $\frac{\rm d}{{\rm
d}t}\theta_{\tau}^{(i)}=-\frac{1}{2}\hat{\Sigma}_{\tau}\Phi\,(y_{\tau}^{(i)}+y(\hat{m}_{\tau})-2t)\,.$

to the second-order equations (31)–(32).

######  Lemma 4.

The EnKBF equations (33) are of affine-invariant gradient structure

(34) $\frac{\rm d}{{\rm
d}t}\theta_{\tau}^{(i)}=-\hat{\Sigma}_{\tau}\nabla_{\theta^{(i)}}\Psi_{\rm
KBF}(\\{\theta_{\tau}^{(j)}\\})\,,$

with potential

(35) $\Psi_{\rm KBF}(\\{\theta^{(j)}\\})=\frac{1}{2}\sum_{j=1}^{M}\Psi_{\rm
data}(\theta^{(j)})+\frac{M}{2}\Psi_{\rm data}(\hat{m}_{\tau})\,.$

###### Proof.

The lemma can be verified by direct calculation and taking note of

$\nabla_{\theta^{(i)}}\Psi_{\rm
data}(\hat{m}_{\tau})=\frac{1}{M}\Phi\,(y(\hat{m}_{\tau})-t)$

in particular. ∎

The equation (34) also possesses an affine-invariant gradient-flow structure
in the space of probability measures in the mean-field limit $M\to\infty$. See
[28, 10, 9] for details.

######  Remark 5.

The authors of [17] proposed a different modification of the EnKF for
classification problems. In our notation, it corresponds to

$\frac{\rm d}{{\rm
d}t}\theta_{\tau}^{(i)}=-\sum_{n=1}^{N}\left\\{\frac{1}{M-1}\sum_{j=1}^{M}\left(y_{\tau,n}^{(j)}-\frac{1}{M}\sum_{l=1}^{M}y_{\tau,n}^{(l)}\right)\left(\frac{t_{n}}{y_{\tau,n}^{(i)}}-\frac{1-t_{n}}{1-y_{\tau,n}^{(i)}}\right)\theta_{\tau}^{(j)}\right\\}.$

The modified EnKF formulation of [13] uses an implementation that is closer to
ours but does not actually propagate ensembles. Also recall that both methods
are to be used for reducing the loss function $\Psi_{\rm data}(\theta)$ via
minimisation instead of sampling from the posterior PDF (6), which is the
subject of this paper.

######  Remark 6.

In addition to the EnKBF formulation considered so far, there exists another
variant which is based on stochastic perturbations in the innovation [8, 26].
We briefly explain how to extend this alternative formulation to logistic
regression. We recall that the ensemble deviations $\Theta_{\tau}$ satisfy the
evolution equation (32). The same propagation of the first and second-order
moments is achieved by the stochastic equations

${\rm
d}(\theta_{\tau}-m_{\tau})=-\Sigma_{\tau}\Psi\left\\{\pi_{\tau}[R]\Psi^{\rm
T}(\theta_{\tau}-m_{\tau})\,{\rm d}\tau+\pi_{\tau}[R]^{1/2}{\rm
d}W_{\tau}\right\\}$

in the mean-field limit $M\to\infty$, where $W_{\tau}$ denotes standard
$N$-dimensional Brownian motion. This suggests to replace the finite ensemble
size formulation (33) by the interacting system of SDEs

${\rm d}\theta_{\tau}^{(i)}=-\hat{\Sigma}_{\tau}\Psi\left(y_{\tau}^{(i)}\,{\rm
d}\tau+\hat{R}_{\tau}\,{\rm d}W_{\tau}^{(i)}-t\,{\rm d}\tau\right).$

Here $W_{\tau}^{(i)}$, $i=1,\ldots,M$, denote independent $N$-dimensional
Brownian motions.

######  Remark 7.

One can replace the potential (35) by

$\Psi_{\rm
EnKBF}(\\{\theta^{(j)}\\})=\frac{2-\alpha}{2}\sum_{j=1}^{M}\Psi_{\rm
data}(\theta^{(j)})+\frac{\alpha\,M}{2}\Psi_{\rm data}(\hat{m}_{\tau})\,.$

with parameter $\alpha\in[1,2)$. This modification entails a form of variance
inflation [8, 29] for $\alpha>1$ since the implied time evolution of the
ensemble deviations (22) becomes less contractive.

Again, the forward Euler method (28) can be used to solve (33) or (34) in
time. However, the inherent stiffness of the equations of motion can lead to
very small time-steps [2]. This issue and the rank deficiency of
$\hat{\Sigma}_{\tau}$ for $M\leq D$ are addressed in the following subsection.

###  4.4. Dropout and time-stepping

We now discuss two numerical issue that are common to all ensemble-base
methods considered in this section.

First, we observe that whenever $M\leq D$, then the empirical covariance
matrix $\hat{\Sigma}_{\tau}$ is rank-deficient and the ensemble
$\\{\theta_{\tau}^{(i)}\\}$ is propagated in the subspace defined by the
ensemble at initial time. The concept of covariance localisation has been
pioneered in the geoscience community in order to lift this restriction [8,
29]. In the context of logistic regression we suggest to instead utilise the
idea of dropout training as pioneered in the machine learning community [32].
We employ the concept of dropout as follows: In each time step, a randomly
chosen number of entries in $\theta_{\tau}^{(i)}$, $i=1,\ldots,M$, is set to
zero. The resulting modified matrix of ensemble deviations (22) is denoted by
$\tilde{\Theta}_{\tau}$. The empirical covariance matrix is now defined by

(36)
$\hat{\Sigma}_{\tau}=\frac{1}{(1-\mu)(M-1)}\tilde{\Theta}_{\tau}\tilde{\Theta}_{\tau}^{\rm
T}\,,$

where $\mu\in(0,1)$ denotes the fraction of entries randomly set to zero. All
other aspects of the previously considered algorithms remain unaltered. The
estimator (36) underestimates the cross-correlations. An unbiased estimator
requires an additional rescaling of the off-diagonal entries by $1/(1-\mu)$.

######  Remark 8.

The random dropout approximation (36) does not provide an unbiased
approximation to the empirical covariance matrix (23). An unbiased
approximation would require that the off-diagonal entries of
$\tilde{\Theta}_{\tau}\tilde{\Theta}_{\tau}^{\rm T}$ get scaled by
$1/((1-\mu)^{2}(M-1))$ in (36). While such an unbiased approximation fits into
the modified equation analysis framework of Section 3.2, we found that (36)
works well for practical applications as its puts a stronger emphasis on the
diagonal entries of the covariance matrix, which is in line with standard
covariance localisation techniques [8, 29].

######  Remark 9.

There are alternative methods for breaking the sub-space property of ensemble
transform methods in case of $M\leq D$. For example, the authors of [17]
suggest to randomly perturb the ensemble members after each time step of EKI.
In our context, this approach can be viewed as combining the EnKBF with
diffusion at low temperature. However, it is not clear that diffusion is
sufficient to effectively explore the full parameter space in the context of
the homotopy approach considered in this paper.

Second, we return to the issue of efficient time-stepping of the interacting
particle systems. We had previously considered the forward Euler method (28).
However, the method can encounter severe step-size restrictions due to its
restricted domain of stability. Following [2], we therefore consider the
following tamed version of the forward Euler discretisation (28) for the EnKBF
formulation from Section 4.3. Starting from (33) we employ

(37)
$\theta_{\tau_{k+1}}^{(i)}=\theta_{\tau_{k}}^{(i)}-\frac{\Delta\tau}{2}\hat{\Sigma}_{\tau_{k}}\Phi\left(I+\Delta\tau\hat{R}_{\tau_{k}}\Phi^{\rm
T}\hat{\Sigma}_{\tau_{k}}\Phi\right)^{-1}(y_{\tau_{k}}^{(i)}+y(\hat{m}_{\tau_{k}})-2t)$

with step size $\Delta\tau>0$ for $k=0,\ldots,K-1$. If the necessary matrix
inversions become too computationally expensive, one can diagonalise the
inverse as proposed in [2]. The same modification can be applied to the
second-order formulation (31)-(32).

##  5\. Ensemble transform Langevin dynamics

While we have focused on the homotopy approach in the previous section, we now
combine the homotopy approach with overdamped Langevin dynamics in order to
sample from the desired target distribution (6) as $\tau\to\infty$. Let us
therefore turn the transport-based differential equation (10) into the
McKean–Vlasov SDE

(38) ${\rm
d}\theta_{\tau}=-\Sigma_{\tau}\left\\{\nabla_{\theta}V_{\tau}(\theta_{\tau})+\frac{1}{2}\Sigma_{\rm
prior}^{-1}(\theta_{\tau}+m_{\tau}-2m_{\rm prior})\right\\}{\rm
d}\tau+\Sigma_{\tau}^{1/2}{\rm d}W_{\tau}$

by adding diffusion and a prior-related drift term to (10). This SDE is
motivated by the affine-invariant EKS [9] and ALDI [10] variants of standard
Langevin dynamics. Note that the crucial differences are that (i) EKS/ALDI
uses the gradient $\nabla_{\theta}\Psi_{\rm prior}$ of the prior in the drift
term while (38) contains

$\frac{1}{2}\Sigma_{\rm prior}^{-1}(\theta_{\tau}+m_{\tau}-2m_{\rm prior})$

instead, and (ii), in a similar vein, the potential $\Psi_{\rm data}(\theta)$
is replaced in (38) by a time-dependent potential $V_{\tau}$ which is defined
as the solution of the elliptic PDE (9). We will find that these modifications
are motivated by the two facts that (i) one can find efficient gradient-free
time-stepping methods for the SDE (38) and (ii) that the SDE (38) samples from
the posterior PDF in the case of Gaussian distributions as we demonstrate
next.

######  Lemma 5.

If the initial $\pi_{0}$ is Gaussian and $\Psi_{\rm data}$ is a potential of
the form

$\Psi_{\rm data}(\theta)=\frac{1}{2}(G\theta-t)^{\rm
T}\Gamma^{-1}(G\theta-t),$

then the marginal PDFs $\pi_{\tau}$ of (38) remain Gaussian for all $\tau>0$,
and the SDE (38) has a stationary measure which is given by
$\pi_{\infty}(\theta)\propto\exp(-\Psi_{\rm post}(\theta))$.

###### Proof.

Under the stated assumption, we have

$-\frac{1}{\pi_{\tau}(\theta)}\nabla_{\theta}\cdot(\pi_{\tau}(\theta)\Sigma_{\tau}\nabla_{\theta}V_{\tau}(\theta))=(\theta-
m_{\tau})^{\rm
T}\nabla_{\theta}V_{\tau}(\theta)-\nabla_{\theta}\cdot\Sigma_{\tau}\nabla_{\theta}V_{\tau}(\theta)$

as well as

$\displaystyle\Psi_{\rm data}(\theta)-\pi_{\tau}[\Psi_{\rm data}]$
$\displaystyle=\frac{1}{2}(G\theta-t)^{\rm T}\Gamma^{-1}(G\theta-t)\,\,-$
$\displaystyle\qquad\qquad\frac{1}{2}\pi_{\tau}[\theta^{\rm T}G^{\rm
T}\Gamma^{-1}G\theta]+m_{\tau}^{\rm T}G^{\rm T}\Gamma^{-1}t-\frac{1}{2}t^{\rm
T}\Gamma^{-1}t$ $\displaystyle=\frac{1}{2}(\theta-m_{\tau})^{\rm T}G^{\rm
T}\Gamma^{-1}(G\theta+Gm_{\tau}-2t)\,\,-$
$\displaystyle\qquad\qquad\frac{1}{2}\pi_{\tau}[\theta^{\rm T}G^{\rm
T}\Gamma^{-1}G\theta]+\frac{1}{2}m_{\tau}^{\rm T}G^{\rm
T}\Gamma^{-1}Gm_{\tau}.$

Equating both expressions leads us to

(39) $\nabla_{\theta}V_{\tau}(\theta)=\frac{1}{2}G^{\rm
T}\Gamma^{-1}(G\theta+Gm_{\tau}-2t).$

The evolution equations arising from (38) for the mean $m_{\tau}$ and the
covariance matrix $\Sigma_{\tau}$ are therefore given by

$\displaystyle\frac{\rm d}{{\rm d}\tau}m_{\tau}$
$\displaystyle=-\Sigma_{\tau}\left(G^{\rm
T}\Gamma^{-1}(Gm_{\tau}-t)+\Sigma_{\rm prior}^{-1}(m_{\tau}-m_{\rm
prior})\right)$ $\displaystyle=-\Sigma_{\tau}\Sigma_{\rm
post}^{-1}(m_{\tau}-m_{\rm post})$

and

$\displaystyle\frac{\rm d}{{\rm d}\tau}\Sigma_{\tau}$
$\displaystyle=-\Sigma_{\tau}\left(G^{\rm T}\Gamma^{-1}G+\Sigma_{\rm
prior}^{-1}\right)\Sigma_{\tau}+\Sigma_{\tau}$
$\displaystyle=-\Sigma_{\tau}\Sigma_{\rm
post}^{-1}\Sigma_{\tau}+\Sigma_{\tau},$

respectively. The associated equilibrium solutions are

$m_{\infty}=m_{\rm post}:=m_{\rm prior}-\Sigma_{\rm post}G^{\rm
T}\Gamma^{-1}(Gm_{\rm prior}-t)$

and

$\Sigma_{\infty}=\Sigma_{\rm post}:=\left(G^{\rm T}\Gamma^{-1}G+\Sigma_{\rm
prior}^{-1}\right)^{-1},$

as desired. We also note that the linearity of the gradient (39) implies that
the PDFs $\pi_{\tau}$ remain Gaussian. Finally, (39) is indeed the gradient of
the potential

$V_{\tau}(\theta)=\frac{1}{4}(G\theta-t)^{\rm
T}\Gamma^{-1}(G\theta-t)+\frac{1}{2}\theta^{\rm T}G^{\rm
T}\Gamma^{-1}(Gm_{\tau}-t)\,,$

which concludes the proof. ∎

######  Remark 10.

One can replace the McKean–Vlasov SDE (38) by the following formulation

(40) ${\rm
d}\theta_{\tau}=-\Sigma_{\tau}\nabla_{\theta}U_{\tau}(\theta_{\tau})\,{\rm
d}\tau+\Sigma_{\tau}^{1/2}{\rm d}W_{\tau},$

where the potential $U_{\tau}$ now satisfies the elliptic PDE

$\nabla_{\theta}\cdot(\pi_{\tau}\Sigma_{\tau}\nabla_{\theta}U_{\tau})=-(\Psi_{\rm
post}-\pi_{\tau}\left[\Psi_{\rm post}\right])\,\pi_{\tau}$

instead of (9). Both formulations are equivalent for Gaussian PDFs
$\pi_{\tau}$. While (38) is closer to the homotopy approach of Section 4,
formulation (40) has a closer resemblance to overdamped Langevin dynamics.

###  5.1. Asymptotic behaviour

The case of a general negative log-likelihood function $\Psi_{\rm data}$ is
more delicate. However, if the associated posterior (6) is non-Gaussian but is
well approximated by a Gaussian, and this also applies to all intermediate
PDFs $\pi_{\tau}$, that is,

(41)
$\pi_{\tau}(\theta)\approx\tilde{\pi}_{\tau}(\theta)\propto\exp\left(-\frac{1}{2}(\theta-
m_{\tau})^{\rm T}\Sigma_{\tau}^{-1}(\theta-m_{\tau})\right),$

then following the arguments from Section 4.2, the first two moments satisfy

(42) $\frac{\rm d}{{\rm
d}\tau}m_{\tau}=-\Sigma_{t}\tilde{\pi}_{\tau}[\nabla_{\theta}\Psi_{\rm post}]$

and

(43) $\frac{\rm d}{{\rm
d}t}\Sigma_{\tau}=-\Sigma_{t}\tilde{\pi}_{\tau}[D_{\theta}^{2}\Psi_{\rm
post}]\Sigma_{t}+\Sigma_{t}\,.$

The equation for the covariance matrix is stable if
$\tilde{\pi}_{\tau}[D_{\theta}^{2}\Psi_{\rm post}]$ is positive-definite for
all $\tau\geq 0$. This is, for example, the case for logistic regression.

We now discuss the limiting $\tau\to\infty$ behavior of (38) in some more
detail. Its equilibrium distributions are characterised by the following
lemma.

######  Lemma 6.

Stationary measures $\pi_{\infty}$ of the mean-field equation (38) satisfy the
PDE

$\frac{1}{2}\nabla_{\theta}\cdot\left(\pi_{\infty}\Sigma_{\infty}\left\\{\nabla_{\theta}\log\pi_{\infty}+\Sigma_{\rm
prior}^{-1}(\theta+m_{\infty}-2m_{\rm
prior})\right\\}\right)=\pi_{\infty}(\Psi_{\rm data}-\pi_{\infty}[\Psi_{\rm
data}])$

or, in terms of expectation values with respect to test functions $\phi$,

(44a)
$\displaystyle\frac{1}{2}\pi_{\infty}\left[\nabla_{\theta}\phi\cdot\Sigma_{\infty}\left\\{\nabla_{\theta}\log\pi_{\infty}+\Sigma_{\rm
prior}^{-1}(\theta+m_{\infty}-2m_{\rm prior})\right\\}\right]$ (44b)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad=\pi_{\infty}\left[\phi\left(\Psi_{\rm
data}-\pi_{\infty}\left[\Psi_{\rm data}\right]\right)\right]\,.$

###### Proof.

The lemma follows from the Fokker–Planck equation associated to (38) and the
PDE (9). ∎

Again making the Gaussian approximation (41), formally setting $\tau=\infty$
and recalling Remark 10, we first obtain from (44) the identity

$\frac{1}{2}\tilde{\pi}_{\infty}\left[\nabla_{\theta}\phi\cdot\Sigma_{\infty}\nabla_{\theta}\log\tilde{\pi}_{\infty}\right]=\tilde{\pi}_{\infty}\left[\phi\left(\Psi_{\rm
post}-\tilde{\pi}_{\infty}\left[\Psi_{\rm post}\right]\right)\right]\,.$

Then,

(45)
$\frac{1}{2}\tilde{\pi}_{\infty}[\nabla_{\theta}\log\tilde{\pi}_{\infty}]=0=\tilde{\pi}_{\infty}[\nabla_{\theta}\Psi_{\rm
post}]$

for $\phi(\theta)=\theta-m_{\infty}$, and

(46)
$\tilde{\pi}_{\infty}[D^{2}_{\theta}\log\tilde{\pi}_{\infty}]=\Sigma_{\infty}^{-1}=\tilde{\pi}_{\infty}[D^{2}_{\theta}\Psi_{\rm
post}]$

for $\phi(\theta)=(\theta-m_{\infty})(\theta-m_{\infty})^{\rm T}$, which
provide a self-consistent system of equations for the mean $m_{\infty}$ and
the covariance matrix $\Sigma_{\infty}$ of the Gaussian approximation (41).
Here we have again used (29) and (30), respectively.

Note that (45) and (46) are compatible with (42) and (43), respectively.
However, (45) and (46) only require Gaussianity (or a near-Gaussianity) for
the equilibrium distribution $\pi_{\infty}$.

###  5.2. Numerical implementation

The most appealing aspect of our modified overdamped Langevin dynamics is the
possibility of implementing the SDEs (38) or (40), respectively, in a
derivative-free manner, that is, without explicit knowledge of the potential
$V_{\tau}$ and its gradient.

We describe a suitable time-stepping method in detail for the formulation (38)
utilising methods from sequential data assimilation, which are well suited to
deal with the $-\Sigma_{\tau}\nabla_{\theta}V_{\tau}$ term in the SDE (38). In
other words, we numerically solve the SDE (38) with a step size $\Delta\tau$
by alternating in each time step between (i) a data assimilation step (such as
SMC or an ensemble transform filter [29]) with negative log-likelihood

$l(\theta)=\Delta\tau\Psi_{\rm data}(\theta)$

and (ii) the SDE

(47) ${\rm d}\theta_{t}=-\frac{1}{2}\Sigma_{\tau}\Sigma_{\rm
prior}^{-1}(\theta_{\tau}+m_{\tau}-2m_{\rm prior})\,{\rm
d}\tau+\Sigma_{t}^{1/2}{\rm d}W_{t}$

over the same time-interval $\Delta\tau$. More precisely, given
$\theta_{\tau_{k}}\sim\pi_{\tau_{k}}$ we first find

(48) $\tilde{\theta}_{\tau_{k+1}}\sim\exp(-\Delta\tau\Psi_{\rm
data})\,\pi_{\tau_{k}}$

followed by solving (47) with initial $\tilde{\theta}_{\tau_{k+1}}$ over the
time-interval $\Delta\tau$ in order to obtain
$\theta_{\tau_{k+1}}\sim\pi_{\tau_{k+1}}$ for $k=0,\ldots,K-1$. Note that (47)
can be solved robustly by the tamed forward Euler-Maruyama method

$\theta_{\tau_{k+1}}=\tilde{\theta}_{\tau_{k+1}}-\frac{\Delta\tau}{2}\tilde{\Sigma}_{\tau_{k+1}}\left(\Sigma_{\rm
prior}+\Delta\tau\tilde{\Sigma}_{\tau_{k+1}}\right)^{-1}(\tilde{\theta}_{\tau_{k+1}}+\tilde{m}_{\tau_{k+1}}-2m_{\rm
prior})+\sqrt{\Delta\tau}\tilde{\Sigma}_{\tau_{k+1}}^{1/2}\Xi_{k+1}$

with $\Xi_{k+1}$ a $D$-dimensional standard Gaussian random variable. As for
ALDI, the $\Sigma_{\tau}^{1/2}$ multiplying the Brownian motion in (47)
requires a correction term in order to sample from the correct distribution
for finite ensemble sizes $M$. See [10] for details.

Our numerical implementations of (48) in terms of ensembles
$\\{\theta_{\tau_{k}}^{(i)}\\}$ rely on ensemble transform filters of the form
[29]

(49)
$\tilde{\theta}_{\tau_{k+1}}^{(j)}=\sum_{i=1}^{M}\theta_{\tau_{k}}^{(i)}s_{ij}$

for $j=1,\ldots,M$ with appropriate coefficients $s_{ij}$, which satisfy

(50) $\sum_{i=1}^{M}s_{ij}=1.$

We collect the coefficients $s_{ij}$ in the $M\times M$ matrix $S$.

######  Lemma 7.

An ensemble transform filter is affine-invariant (compare Section 3.1) if the
associated coefficients $\overline{s}_{ij}$ for the filter in the transformed
variable $\overline{\theta}$ satisfy $s_{ij}=\overline{s}_{ij}$ for all
$i,j=1,\ldots,M$.

###### Proof.

If

$\theta_{\tau_{k}}^{(i)}=A\overline{\theta}_{\tau_{k}}^{(i)}+b$

at $\tau_{k}$, then

$\theta_{\tau_{k+1}}^{(j)}=\sum_{i=1}^{M}\theta_{\tau_{k}}^{(i)}s_{ij}=\sum_{i=1}^{M}(A\overline{\theta}_{\tau_{k}}^{(i)}+b)s_{ij}=A\sum_{i=1}^{M}\overline{\theta}_{\tau_{k}}^{(i)}\overline{s}_{ij}+b=A\overline{\theta}_{\tau_{k+1}}^{(j)}+b\,,$

and we conclude that

$\theta_{\tau_{k+1}}^{(i)}=A\overline{\theta}_{\tau_{k+1}}^{(i)}+b$

at $\tau_{k+1}$ for all $i=1,\ldots,M$. ∎

ensemble size/method | $M=50$ | $M=100$ | $M=200$ | $M=400$
---|---|---|---|---
McKean–Vlasov SDE | $\left(\begin{array}[]{c}-3.34\\\ -3.39\\\ 3.22\end{array}\right)$ | $\left(\begin{array}[]{c}-3.35\\\ -3.39\\\ 3.21\end{array}\right)$ | $\left(\begin{array}[]{c}-3.35\\\ -3.39\\\ 3.21\end{array}\right)$ | $\left(\begin{array}[]{c}-3.35\\\ -3.39\\\ 3.22\end{array}\right)$
FPF | $\left(\begin{array}[]{c}-3.31\\\ -3.35\\\ 3.20\end{array}\right)$ | $\left(\begin{array}[]{c}-3.33\\\ -3.37\\\ 3.21\end{array}\right)$ | $\left(\begin{array}[]{c}-3.33\\\ -3.37\\\ 3.21\end{array}\right)$ | $\left(\begin{array}[]{c}-3.33\\\ -3.37\\\ 3.21\end{array}\right)$
second-order | $\left(\begin{array}[]{c}-3.36\\\ -3.40\\\ 3.22\end{array}\right)$ | $\left(\begin{array}[]{c}-3.36\\\ -3.41\\\ 3.22\end{array}\right)$ | $\left(\begin{array}[]{c}-3.36\\\ -3.41\\\ 3.22\end{array}\right)$ | $\left(\begin{array}[]{c}-3.36\\\ -3.41\\\ 3.22\end{array}\right)$
EnKBF | $\left(\begin{array}[]{c}-3.27\\\ -3.31\\\ 3.20\end{array}\right)$ | $\left(\begin{array}[]{c}-3.27\\\ -3.32\\\ 3.20\end{array}\right)$ | $\left(\begin{array}[]{c}-3.27\\\ -2.31\\\ 3.19\end{array}\right)$ | $\left(\begin{array}[]{c}-3.27\\\ -3.31\\\ 3.19\end{array}\right)$
Table 1. Informative prior: We display the ensemble mean averaged over
$L=1000$ repeated experiments for the four proposed methods as a function of
the ensemble size $M$. The reference value from an ALDI simulation is
$(-3.32,-3.36,3.20)^{\rm T}$. We find that all methods reproduce the reference
solution provided by ALDI within small variations.

The EnKF and the nonlinear ensemble transform filter (NETF) of [36], for
example, are affine-invariant. More specifically, the NETF leads to the
transformation matrix

$S=w_{\tau_{k}}{\rm 1}^{\rm T}_{M}+\sqrt{M}\left({\rm
diag}\,(w_{\tau_{k}})-w_{\tau_{k}}w_{\tau_{k}}^{\rm T}\right)^{1/2}$

where $w_{\tau_{k}}\in{\mathbb{R}}^{M}$ is the vector of normalised importance
weights with entries

$w_{\tau_{k}}^{(i)}=\frac{\exp(-\Delta\tau\Psi_{\rm
data}(\theta_{\tau_{k}}^{(i)}))}{\sum_{j=1}^{M}\exp(-\Delta\tau\Psi_{\rm
data}(\theta_{\tau_{k}}^{(j)}))},$

${\rm 1}_{M}\in{\mathbb{R}}^{M}$ is the vector of ones, and
$\mbox{diag}\,(w_{\tau_{k}})\in{\mathbb{R}}^{M\times M}$ is the diagonal
matrix with diagonal entries $w_{\tau_{k}}$. The affine invariance follows
from the invariance of the importance weights $w_{\tau_{k}}$. The affine
invariance of the EnKF has been previously demonstrated, for example, in [22].

ensemble size/method | $M=50$ | $M=100$ | $M=200$ | $M=400$
---|---|---|---|---
McKean–Vlasov SDE | $0.89$ | $0.83$ | $0.80$ | $0.78$
FPF | $0.98$ | $0.90$ | $0.86$ | $0.83$
second-order | $0.75$ | $0.75$ | $0.74$ | $0.74$
EnKBF | $0.80$ | $0.79$ | $0.78$ | $0.78$
Table 2. Informative prior: We display the spectral norm of the final ensemble
covariance matrices averaged over $L=1000$ repeated experiments for the four
proposed methods as a function of the ensemble size $M$. The ALDI reference
value is $0.82$. We find that the FPF slightly overestimates the variance for
smaller ensemble sizes, while both the EnKBF and the second-order formulation
slightly underestimate it across all ensemble sizes.

######  Remark 11.

The ensemble transform particle filter (ETPF) of [27, 29] can also be made
affine-invariant by using the cost function

(51)
$V(S)=\frac{1}{2}\sum_{i,j=1}^{M}s_{ij}\,(\theta_{\tau_{k}}^{(i)}-\theta_{\tau_{k}}^{(j)})^{\rm
T}\hat{\Sigma}_{\tau_{k}}^{-1}(\theta_{\tau_{k}}^{(i)}-\theta_{\tau_{k}}^{(j)})$

subject to the constraints $s_{ij}\geq 0$, (50), and

(52) $\frac{1}{M}\sum_{j=1}^{M}s_{ij}=w_{\tau_{k}}^{(i)}$

in the optimal transport definition of the coefficients $s_{ij}$ in (49).
Again, the affine invariance follows from the affine invariance of both the
cost function $V(S)$ and the importance weights $w_{\tau_{k}}$.

##  6\. Numerical example

In this section, we provide a couple of relatively simple numerical
illustrations for the methods proposed in this paper. We start with a low-
dimensional toy problem for which all proposed methods can easily be
implemented and tested.

ensemble size/method | $M=50$ | $M=100$ | $M=200$ | $M=400$
---|---|---|---|---
McKean–Vlasov SDE | $\left(\begin{array}[]{c}-2.61\\\ -2.65\\\ 2.15\end{array}\right)$ | $\left(\begin{array}[]{c}-2.60\\\ -2.64\\\ 2.15\end{array}\right)$ | $\left(\begin{array}[]{c}-2.59\\\ -2.62\\\ 2.15\end{array}\right)$ | $\left(\begin{array}[]{c}-2.59\\\ -2.62\\\ 2.15\end{array}\right)$
FPF | $\left(\begin{array}[]{c}-2.16\\\ -2.15\\\ 1.64\end{array}\right)$ | $\left(\begin{array}[]{c}-2.28\\\ -2.30\\\ 1.83\end{array}\right)$ | $\left(\begin{array}[]{c}-2.39\\\ -2.41\\\ 1.95\end{array}\right)$ | $\left(\begin{array}[]{c}-2.44\\\ -2.47\\\ 2.00\end{array}\right)$
second-order | $\left(\begin{array}[]{c}-2.29\\\ -2.32\\\ 1.82\end{array}\right)$ | $\left(\begin{array}[]{c}-2.30\\\ -2.32\\\ 1.83\end{array}\right)$ | $\left(\begin{array}[]{c}-2.31\\\ -2.33\\\ 1.83\end{array}\right)$ | $\left(\begin{array}[]{c}-2.31\\\ -3.34\\\ 1.84\end{array}\right)$
EnKBF | $\left(\begin{array}[]{c}-2.14\\\ -2.16\\\ 1.73\end{array}\right)$ | $\left(\begin{array}[]{c}-2.16\\\ -2.18\\\ 1.74\end{array}\right)$ | $\left(\begin{array}[]{c}-2.17\\\ -2.19\\\ 1.75\end{array}\right)$ | $\left(\begin{array}[]{c}-2.17\\\ -2.19\\\ 1.76\end{array}\right)$
Table 3. Less informative prior: We display the ensemble mean averaged over
$L=1000$ repeated experiments for the four proposed methods as a function of
the ensemble size $M$. We find that all methods qualitatively reproduce the
reference value $(-2.56,-2.59,2.15)^{\rm T}$ provided by ALDI. While the
second-order method and the EnKBF suffer from a systematic bias, the FPF
approaches the reference solution as the ensemble size increases. The
McKean–Vlasov SDE formulation yields good approximations for all ensemble
sizes.

######  Example 1.

We follow the example from [5] Section 4.2.1. Specifically, consider the
Gaussian likelihood functions

$\pi(x|C_{i})\propto\exp\left(-\frac{1}{2}(x-\mu_{i})^{\rm
T}B^{-1}(x-\mu_{i})\right)$

for inputs $x$ to belong to classes $C_{i}$, $i=1,2$. The implied posterior
PDF $\pi(C_{1}|\phi_{x})$ is given by (2) with $\phi_{x}=(x^{\rm T},1)^{\rm
T}$, that is $D=J+1$, and the true parameter value is

(53) $\theta=\left(\begin{array}[]{c}B^{-1}(\mu_{1}-\mu_{2})\\\
-\frac{1}{2}\mu_{1}^{\rm T}B^{-1}\mu_{1}+\frac{1}{2}\mu_{2}^{\rm
T}B^{-1}\mu_{2}+\ln\frac{\pi(C_{1})}{\pi(C_{2})}\end{array}\right)\,.$

Furthermore, given data points $(x_{n},t_{n})$, $n=1,\ldots,N$, the maximum
likelihood estimators for $\pi(C_{1})$, $\mu_{i}$, and $B$ are given by

$\hat{\pi}(C_{1})=\frac{N_{1}}{N},\qquad N_{1}=\sum_{n=1}^{N}t_{n},$
$\hat{\mu}_{1}=\frac{1}{N_{1}}\sum_{n=1}^{N}t_{n}x_{n},\qquad\hat{\mu}_{2}=\frac{1}{N-N_{1}}\sum_{n=1}^{N}(1-t_{n})x_{n}\,,$

and

$\hat{B}=\frac{1}{N}\left\\{\sum_{n=1}^{N}t_{n}(x_{n}-\hat{\mu}_{1})(x_{n}-\hat{\mu}_{1})^{\rm
T}+\sum_{n=1}^{N}(1-t_{n})(x_{n}-\hat{\mu}_{2})(x_{n}-\hat{\mu}_{2})^{\rm
T}\right\\}\,,$

respectively. See Section 4.2.2 in [5] for details. In this example, we will
however directly estimate the true parameter value (53) using Bayesian
inference.

We consider the special case $J=2$ and use $B={\rm I}$, ${\rm I}$ the identity
matrix, $\mu_{1}=(-1,-1)^{\rm T}$, $\mu_{2}=(2,2)^{\rm T}$, and
$\pi(C_{1})=\pi(C_{2})=1/2$ to generate $N=100$ data points. We note that the
chosen values for $\mu_{1}$, $\mu_{2}$, and $B$ lead to

(54) $\theta=\left(\begin{array}[]{c}-3\\\ -3\\\ 3\end{array}\right).$

We now test the proposed algorithms using (i) an informative Gaussian prior
with mean (54) and covariance matrix $\Sigma_{\rm prior}={\rm I}$, and (ii) a
less informative Gaussian prior with mean $m_{\rm prior}=0$ and covariance
matrix $\Sigma_{\rm prior}=4\,{\rm I}$. We time step the affine-invariant FPF,
the second-order filter (31)-(32), and the EnKBF (33) with the forward Euler
method using a step size $\Delta\tau=10^{-3}$. A further reduction of the step
size did not change the results in a statistically significant manner. The
McKean–Vlasov SDE (38) is also implemented with the NETF as the inner data
assimilation step. Furthermore, we implemented the ALDI method [10] for
comparison since it is known to exactly sample from the posterior PDF (6).
Both ALDI and the McKean–Vlasov SDE (40) were run up to time $\tau=10$ with
step-size $\Delta\tau=10^{-2}$ at which point the interacting particle systems
were considered to be in equilibrium.

The implementation of the FPF, as well as the SDE (38), explicitly involve the
evaluation of the negative log-likelihood function (3). We found that we
needed to replace (4) by

$y_{n}=0.99\,\sigma(\theta^{\rm T}\phi_{x_{n}})+0.005$

in order to avoid numerical instabilities due to exceedingly small values of
$\ln y_{n}$. The band-width, $\epsilon$, in (25) has been set to
$\epsilon=0.1$.

All experiments were repeated $L=1000$ times in order to average out random
sampling effects. The numerically computed posterior means and spectral norm
of the posterior covariance matrices can be found in Tables 1 and 2,
respectively, for the informative prior, and in Tables 3 and 4, respectively,
for the less informative prior. Reference values from an ALDI simulation are
provided in the captions. While all methods reproduce those reference values
very well in case of the informative prior, this picture changes for the less
informative prior, where one expects a stronger nonlinear and non-Gaussian
behavior of the associated Bayesian inference problem. One finds that the
accuracy of the FPF improves as the ensemble size increases while the second-
order method and the EnKBF suffer from a systematic bias. We also find that
the McKean–Vlasov SDE formulation (38) behaves rather well over the full range
of ensemble sizes. We conclude that all tested methods are able to
qualitatively reproduce the exact parameter value (54) even in the case of the
less informative prior. The results also indicate that the posterior
distribution is close to Gaussian, while the intermittent PDFs $\pi_{\tau}$
must deviate significantly from a Gaussian distribution.

In terms of computational cost, the EnKBF and the second-order method perform
comparable and are the least expensive by far for all ensemble sizes. While
the McKean–Vlasov SDE is twice as expensive as the EnKBF at $M=50$, this ratio
goes up to a factor of five at $M=400$. The FPF, on the other hand, is ten
times as expensive as the EnKFB at $M=50$ which increases steeply to a factor
of one hundred at $M=400$.

ensemble size/method | $M=50$ | $M=100$ | $M=200$ | $M=400$
---|---|---|---|---
McKean–Vlasov SDE | $1.36$ | $1.14$ | $1.03$ | $0.97$
FPF | $3.61$ | $2.58$ | $2.01$ | $1.66$
second-order | $0.46$ | $0.45$ | $0.45$ | $0.45$
EnKBF | $0.60$ | $0.59$ | $0.59$ | $0.59$
Table 4. Less informative prior: We display the spectral norm of the final
ensemble covariance matrices averaged over $L=1000$ repeated experiments for
the four proposed methods as a function of the ensemble size $M$. The ALDI
reference value is $1.18$.We find that the FPF overestimates the variance,
while both the EnKBF and the second-order formulation underestimate it. While
the FPF appears to approach the correct reference value for increased ensemble
sizes, both the EnKBF and the second-order filter appear to suffer from a
systematic bias. The McKean–Vlasov SDE approach leads to a slightly degraded
performance in terms of ensemble variations as the ensemble size increases.

We now consider a high-dimensional extension of the previous example, for
which we further investigate the performance of the EnKBF (34) in terms of the
tamed time-stepping method (37) for ensemble sizes $M\leq D$, localisation via
dropouts, and data sub-sampling.

######  Example 2.

We consider the simple feature map

$\phi_{x}=x,$

that is $J=D$ such that

$y_{x}(\theta)=\sigma(\theta^{\rm T}x)$

for $D=50$. The data is generated by first drawing a $\theta=\theta_{\rm ref}$
from the standard normal distribution in ${\mathbb{R}}^{D}$. Next, $N=1000$
data points are generated by randomly drawing $x_{n}$, $n=1,\ldots,N$, again
from a standard normal distribution in ${\mathbb{R}}^{D}$. Those points are
assigned the label $t_{n}=1$ with probability $y_{x_{n}}(\theta_{\rm ref})$,
and $t_{n}=0$ otherwise. We use the EnKBF tamed time-stepping method (37) with
step size $\Delta\tau=1/200$ for parameter inference. The initial ensemble
$\theta_{0}^{(i)}$, $i=1,\ldots,M$ is drawn from the prior Gaussian
distribution with mean $m_{\rm prior}=0$ and the covariance matrix
$\Sigma_{\rm prior}$ equal to the identity matrix. We vary the ensemble size
$M$ between $M=20$ and $M=100$ and use dropout localisation for the
computation of the empirical covariance matrices $\hat{\Sigma}_{\tau_{n}}$
following (36). More specifically, an entry of $\theta_{\tau_{n}}^{(i)}$ is
set to zero if $\eta<0.5$ where the $\eta$’s are i.i.d. uniform random
variables from the interval $[0,1]$. We report in Table 5 the
$l_{2}$-difference between $\theta_{\rm ref}$ and the ensemble mean
$\hat{m}_{\tau}$ at final time $\tau=1$. We also state the spectral norm of
$\hat{\Sigma}_{\tau}$ at the final time in Table 6. In a second set of
experiments, we apply mini-batching to the localised EnKBF with random batches
of size $N^{\prime}=100$. We also compare the results to those from the
standard EnKBF formulation (33). Furthermore, ALDI is run for a sample size of
$M=100$ in order to provide the numerical benchmark values listed in the
captions of Tables 5 and 6. All results have been averaged over $L=1000$
independent simulations. We also report the standard deviations in addition to
the averaged values. Our numerical results indicate that localisation by
dropout is effective for small ensemble sizes, and that mini-batching does not
degrade performance while providing significant computational savings. As for
the previous example, the EnKBF is able to capture the posterior mean quite
well while the posterior variances are systematically underestimated. Ensemble
inflation, as widely used for the EnKF [8, 29], could help to further improve
the variance estimates. We found that our results are insensitive to the
choice of the step size $\Delta\tau$ in (37). However changing the threshold
value $\mu=0.5$ in the dropout localisation (36) affects the results
significantly. Smaller values of $\mu$ lead in particular to reduced
$l_{2}$-differences for larger ensemble sizes, $M$, while being less effective
for smaller $M$. See Table 7.

ensemble size/method | $M=20$ | $M=40$ | $M=60$ | $M=80$ | $M=100$
---|---|---|---|---|---
EnKBF | 6.26$\pm 0.75$ | 4.55$\pm$0.83 | 2.67$\pm$0.68 | 1.99$\pm$0.56 | 1.69$\pm$0.48
EnKBF | 1.29$\pm$0.29 | 1.19$\pm$0.23 | 1.28$\pm$0.30 | 1.35$\pm$0.34 | 1.39$\pm$0.46
(with dropout) |  |  |  |  |
EnKBF | 2.14$\pm$0.46 | 1.56$\pm$0.30 | 1.42$\pm$0.25 | 1.37$\pm$0.22 | 1.35$\pm 0.21$
(with dropout & mini-batch) |  |  |  |  |
Table 5. We display the $l_{2}$-difference between the posterior ensemble
means and the true parameter value averaged over $L=1000$ repeated experiments
and its standard deviation for three different implementations of the EnKBF as
a function of the ensemble size $M$. The reference value provided by ALDI is
$1.10$. We find that localisation by dropout dramatically improves the
behavior of the EnKBF especially for small ensemble sizes. Using mini-batches
reduces the approximation quality slightly while providing significant
computational savings.

##  7\. Generalisations and extensions

In this section, we discuss possible extensions of the proposed methods to
nonlinear logistic regression including derivative-free implementations,
multi-class logistic regression, and sigmoidal Cox processes.

###  7.1. Nonlinear logistic regression

We generalise linear logistic regression from Section 2 to the more general
logistic regression problem with $y_{x}(\theta):=\sigma(f_{x}(\theta))$ for
given functions $f_{x}:{\mathbb{R}}^{D}\to{\mathbb{R}}$. Hence the gradient
$\nabla_{\theta}y_{x}(\theta)$ becomes

$\nabla_{\theta}y_{x}(\theta)=\phi_{x}(\theta)\,(1-y_{x}(\theta))\,y_{x}(\theta),$

where the previously considered feature maps $\phi_{x}$ now also depend on the
parameters $\theta$, that is,

$\phi_{x}(\theta):=\nabla_{\theta}f_{x}(\theta).$

The gradient of $\Psi_{\rm data}$ is thus explicitly given by

$\nabla_{\theta}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}\phi_{x_{n}}(\theta)\,(y_{x_{n}}(\theta)-t_{n})$

and the Hessian matrix of second-order derivatives becomes

$D^{2}_{\theta}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}\left\\{\phi_{x_{n}}\,(\theta)\,y_{x_{n}}(\theta)\,(1-y_{x_{n}}(\theta))\,\phi_{x_{n}}(\theta)^{\rm
T}+D_{\theta}\phi_{x_{n}}(\theta)\,(y_{x_{n}}(\theta)-t_{n})\right\\}.$

We now discuss a derivative-free implementation of the EnKBF (33) for
nonlinear logistic regression. The key difference to linear logistic
regression is that the matrix $\Phi\in{\mathbb{R}}^{D\times N}$ now depends on
$\theta$ and is given by

$\Phi(\theta)=(\nabla_{\theta}f_{x_{1}}(\theta),\ldots,\nabla_{\theta}f_{x_{N}}(\theta)).$

In order to avoid the computation of the gradients of $f_{x}$, one replaces
$\Sigma_{\tau}\Phi(\theta_{\tau})$ in (33) by the empirical correlation matrix

$C_{\tau}:=\frac{1}{M-1}\sum_{i=1}^{M}(\theta_{\tau}^{(i)}-\hat{m}_{\tau})\left(f_{x_{1}}(\theta_{\tau}^{(i)}),\ldots,f_{x_{N}}(\theta_{\tau}^{(i)})\right)\in{\mathbb{R}}^{D\times
N}$

between $\theta_{\tau}$ and $f_{x_{n}}(\theta_{\tau})$ for $n=1,\ldots,N$. See
[8] for the basic idea in the context of the EnKF and also [17, 13] for the
specific application of this idea in the context of machine learning. See also
the recent work [24, 30] on tuneable derivative-free approximations. We also
note that this substitution becomes exact in the case of linear functions
$f_{x}(\theta)$. Such a linear approximation is, for example, justified for
functions $f_{x}(\theta)$ defined by sufficiently wide neural networks [18].

ensemble size/method | $M=20$ | $M=40$ | $M=60$ | $M=80$ | $M=100$
---|---|---|---|---|---
EnKBF | 0.014$\pm$0.004 | 0.034$\pm$0.009 | 0.058$\pm$0.012 | 0.079$\pm$0.014 | 0.097$\pm$0.015
EnKBF | 0.043$\pm$0.012 | 0.071$\pm$0.019 | 0.088$\pm$0.022 | 0.100$\pm$0.024 | 0.109$\pm$0.025
(with dropout) |  |  |  |  |
EnKBF | 0.041$\pm$0.012 | 0.066$\pm$0.018 | 0.084$\pm$0.021 | 0.096$\pm$0.024 | 0.105$\pm$0.025
(with dropout and mini-batch) |  |  |  |  |
Table 6. We provide the spectral norm of the computed posterior covariance
matrices in the setting described in Table 5. The reference value provided by
ALDI is $0.22$. We find that both all three EnKBF implementations
underestimate the variance. However, dropout localisation increases the
ensembles spread which in turn improves the estimated ensemble means as
displayed in Table 5.

###  7.2. Multi-class logistic regression

The extension from two-class to multi-call classification is straightforward,
see [5]. The posterior probabilities for $L$ classes $C_{l}$, $l=1,\ldots,L$,
are given by

$\pi(C_{l}|\phi_{x})=\frac{\exp\left(a_{x,l}\right)}{\sum_{j=1}^{L}\exp\left(a_{x,j}\right)},$

where the activation is given by

$a_{x,l}:=\theta_{l}^{\rm T}\phi_{x}.$

The adjustable parameters are now $\theta=(\theta_{1}^{\rm
T},\ldots,\theta_{L}^{\rm T})^{\rm T}\in{\mathbb{R}}^{LD}$. The negative log-
likelihood function becomes

$\Psi_{\rm data}(\theta)=-\sum_{n=1}^{N}\sum_{l=1}^{L}t_{nl}\ln
y_{nl}(\theta)$

with

$y_{nl}(\theta):=\frac{\exp\left(a_{x_{n},l}\right)}{\sum_{j=1}^{L}\exp\left(a_{x_{n},j}\right)}$

and data $(x_{n},\\{t_{nl}\\}_{l=1}^{L})$, where $t_{nl}\in\\{0,1\\}$ subject
to $\sum_{l=1}^{L}t_{nl}=1$ for all $n=1,\ldots,N$. The gradients satisfy

$\nabla_{\theta_{l}}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}(y_{nl}-t_{nl})\phi_{n},$

$l=1,\ldots,L$. The Hessian matrix of second-order derivatives got $D\times D$
block entries

$D_{\theta_{l}}D_{\theta_{j}}\Psi_{\rm
data}(\theta)=\sum_{n=1}^{N}y_{nl}(\delta_{lj}-y_{nj})\phi_{n}\phi_{n}^{\rm
T}$

with $\delta_{lj}=1$ if $l=j$ and $\delta_{lj}=0$ otherwise. We conclude that
these formulas closely resemble the corresponding expressions for binary
classification and, hence, all previously considered algorithms easily
generalise to multi-class regression.

ensemble size/method | $M=20$ | $M=40$ | $M=60$ | $M=80$ | $M=100$
---|---|---|---|---|---
EnKBF | 3.30$\pm$1.01 | 1.73$\pm$0.60 | 1.26$\pm$0.23 | 1.14$\pm$0.31 | 1.12$\pm$0.36
(with dropout) |  |  |  |  |
EnKBF | 3.39$\pm$0.86 | 1.86$\pm$0.61 | 1.38$\pm$0.25 | 1.24$\pm$0.23 | 1.19$\pm$0.27
(with dropout and mini-batch) |  |  |  |  |
Table 7. Same results as displayed in Table 5 except that the dropout rate has
been reduced from $\mu=0.5$ to $\mu=0.2$. While the estimation errors increase
for the smallest ensemble sizes the reduced rate improves the performance for
larger ensemble sizes.

###  7.3. Sigmoidal Cox processes

We consider Cox processes (CPs) on a domain
$\mathcal{X}\subset{\mathbb{R}}^{J}$. Following [1], the CP is parametrised by
an intensity function $\lambda_{x}:\mathcal{X}\to{\mathbb{R}}^{+}$ which we
represent using a sigmoidal random feature model, that is,

$\lambda_{x}(\theta)=\lambda^{\ast}\sigma(\theta^{\rm T}\phi_{x}),$

where $\lambda^{\ast}>0$ provides the upper threshold value for the intensity
function. Note that a Gaussian process is used by [1] while we rely here on
the closely related random feature maps $\phi_{x}$ [21], which easier
integrate into the homotopy approach and the EnKF in particular [12].

See [7] for an efficient Bayesian inference for CPs using variable
augmentation. Here we suggest an alternative approach using an associated
unbiased estimator of the negative log-likelihood in the forward Euler time-
stepping (20).

More specifically, given $N$ events $x_{n}$, $n=1,\ldots,N$, the negative log-
likelihood function is provided by

$\Psi_{\rm data}(\theta)=\int_{\mathcal{X}}\lambda_{x}(\theta){\rm
d}x-\sum_{n=1}^{N}\ln\lambda_{x_{n}}(\theta)\,.$

Following the notation of Section 3.2 and using $I$ uniformly distributed
samples $\hat{x}_{i}$, $i=1,\ldots,I$, on $\mathcal{X}$, we obtain the random
unbiased estimator

(55) $\Psi_{\rm
data}^{\gamma}(\theta)=\frac{\lambda^{\ast}}{I}\sum_{i=1}^{I}y_{\hat{x}_{i}}(\theta)-\sum_{n=1}^{N}\ln
y_{x_{n}}(\theta)-N\ln\lambda^{\ast}$

for this intractable likelihood $\Psi_{\rm data}$. Here we also used the
previously introduced notation (4), that is,
$\lambda_{x}(\theta)=\lambda^{\ast}y_{x}(\theta)$.

We now apply the homotopy approach (10) to the associated Bayesian inference
problem. Following the discussion in Section 3.2, the evolution converges to
the correct posterior if in each time step (20) one uses a different
realisation $\gamma_{\tau_{k}}$ of the unbiased estimator (55) instead of
$\Psi_{\rm data}$ in (9) and takes the limit $\Delta\tau\to 0$.

######  Example 3.

Let us briefly discuss a related problem in a strictly Gaussian setting. We
start from the intractable negative log-likelihood function

$\Psi_{\rm data}(\theta):=\frac{1}{2}\int_{0}^{1}\left(\theta^{\rm
T}\phi_{x}-t\right)^{2}{\rm d}x.$

Its gradient is provided by

$\nabla_{\theta}\Psi_{\rm data}(\theta)=\int_{0}^{1}\phi_{x}\left(\theta^{\rm
T}\phi_{x}-t\right){\rm d}x$

and the associated mean-field EnKBF becomes

$\frac{\rm d}{{\rm
d}\tau}\theta_{t}=-\frac{1}{2}\Sigma_{\tau}\int_{0}^{1}\phi_{x}\left(\left(\theta_{\tau}+m_{\tau}\right)^{\rm
T}\phi_{x}-2t\right){\rm d}x\,.$

So far, everything is exact and no approximation has been made. We now replace
the integral by its Monte Carlo approximation using $I$ randomly sampled
$\hat{x}_{i,k}\sim{\rm U}[0,1]$, $i=1,\ldots,I$, and introduce a time step
$\Delta\tau$ leading to

$\theta_{\tau_{k+1}}=\theta_{\tau_{k}}-\frac{\Delta\tau}{2I}\Sigma_{\tau_{k}}\sum_{i=1}^{I}\phi_{\hat{x}_{i,k}}\left(\left(\theta_{\tau_{k}}+m_{\tau_{k}}\right)^{\rm
T}\phi_{\hat{x}_{i,k}}-2t\right),$

where we re-sample the $\hat{x}_{i,k}$’s in each time-step $\tau_{k}$,
$k=0,\ldots,K-1$. A theoretical investigation now boils down to studying the
limit $\Delta\tau\to 0$, which can be performed along the modified equation
analysis of [20]. We note that the approximation error can be reduced by
either increasing $I$ or decreasing the step size $\Delta\tau$.

##  8\. Conclusions

We have presented affine-invariant extensions of the popular EnKF and related
methods such as the FPF to logistic regression. In addition, we have proposed
a novel SDE-based affine-invariant and derivative-free sampling method which
utilises the robust ensemble transform approach for Bayesian inference as an
inner time-stepping method. We have also introduced a number of algorithmic
choices such as localisation via dropout, linearly implicit time-stepping
methods, and mini-batching the data, which help to improve the efficiency of
the proposed methods. While we have only considered two simple numerical
examples, future work will consider more challenging applications such as the
training of the Sentence Gestalt model of natural language comprehension [25].
A numerical exploration of the affine-invariant stochastic homotopy
formulation (18) would also be of practical and theoretical interest.

A number of theoretical questions also remain to be answered such as the
convergence of the ensemble methods to their mean-field limits, the impact of
the various approximations on the computed posterior representations, and the
ergodicity and convergence to equilibrium properties of (38) and (40),
respectively. Finally, an embedding of the ensemble transform Langevin
dynamics into a Markov chain Monte Carlo framework would allow for an exact
sampling from the posterior.

### Acknowledgements

SR is supported by Deutsche Forschungsgemeinschaft (DFG) - Project-ID
318763901 - SFB1294.

## References

  * [1] R.P. Adams, I. Murray, and D.J.C. MacKay.  Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities.  In Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, pages 9–16, New York, NY, USA, 2009. Association for Computing Machinery.  doi:10.1145/1553374.1553376.
  * [2] J. Amezcua, E. Kalnay, K. Ide, and S. Reich.  Ensemble transform Kalman–Bucy filters.  Q.J.R. Meteor. Soc., 140:995–1004, 2014.  doi:10.1002/qj.2186.
  * [3] K. Bergemann and S. Reich.  A localization technique for ensemble Kalman filters.  Quarterly Journal of the Royal Meteorological Society, 136:701–707, 2010.  doi:10.1002/qj.591.
  * [4] A. Beskos, D. Crisan, and A. Jasra.  On the stability of sequential Monte Carlo methods in high dimensions.  The Annals of Applied Probability, 24:1396–1445, 2014.  doi:10.1214/13-AAP951.
  * [5] Ch. M. Bishop.  Pattern Recognition and Machine Learning.  Springer-Verlag, New York, 2006.
  * [6] F. Daum, J. Huang, and A. Noushin.  Exact particle flow for nonlinear filters.  In Ivan Kadar, editor, Signal Processing, Sensor Fusion, and Target Recognition XIX, volume 7697, pages 92–110. International Society for Optics and Photonics, SPIE, 2010.  doi:10.1117/12.839590.
  * [7] Ch. Donner and M. Opper.  Efficient Bayesian inference of sigmoidal Gaussian Cox processes.  J. Mach. Learn. Research, 19(67):1–34, 2018.  URL: http://jmlr.org/papers/v19/17-759.html.
  * [8] G. Evensen.  Data Assimilation: The Ensemble Kalman Filter.  Springer-Verlag, Berlin, second edition, 2009.  doi:10.1007/978-3-642-03711-5.
  * [9] A. Garbuno-Inigo, F. Hoffmann, W. Li, and A.M. Stuart.  Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler.  SIAM J. Appl. Dyn. Syst., 19(1):412–441, 2020.  doi:10.1137/19M1251655.
  * [10] A. Garbuno-Inigo, N. Nüsken, and S. Reich.  Affine invariant interacting Langevin dynamics for Bayesian inference.  SIAM J. Appl. Dyn. Syst., 19:1633–1658, 2020.  doi:10.1137/19M1304891.
  * [11] J. Goodman and J. Weare.  Ensemble samplers with affine invariance.  Comm. Appl. Math. and Comput. Science, 5:65–80, 2010.  doi:10.2140/camcos.2010.5.65.
  * [12] G. Gottwald and S. Reich.  Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation.  Physica D, 423:132911, 2021.  doi:10.1016/j.physd.2021.132911.
  * [13] E. Haber, F. Lucka, and L. Ruthotto.  Never look back - A modified EnKF method and its application to the training of neutral networks without back propagation, 2018.  arXiv:1805.08034.
  * [14] F. Heber, Ž. Trst’anová, and B. Leimkuhler.  Posterior sampling strategies based on discretized stochastic differential equations for machine learning applications.  Journal of Machine Learning Research, 21(228):1–33, 2020.  URL: http://jmlr.org/papers/v21/19-339.html.
  * [15] D.Z. Huang, T. Schneider, and A.M. Stuart.  Uncented Kalman inversion, 2021.  arXiv:2102.01580.
  * [16] R.E. Kass and A.E. Raftery.  Bayes factors.  Journal of the American Statistical Association, 90:773–795, 1995\.  doi:10.1080/01621459.1995.10476572.
  * [17] N.B. Kovachki and A.M. Stuart.  Ensemble Kalman inversion: A derivative-free technique for machine learning tasks.  Inverse Probl., 35(9):095005, aug 2019.  doi:10.1088/1361-6420/ab1c3a.
  * [18] J. Lee, L. Xiao, S.S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  Journal of Statistical Mechanics: Theory and Experiment, 2020(12):124002, 2020.  doi:10.1088/1742-5468/abc62b.
  * [19] B. Leimkuhler, C. Matthews, and J. Weare.  Ensemble preconditioning for Markov chain Monte Carlo simulation.  Statistics and Computing, 28:277–290, 2018.  doi:10.1007/s11222-017-9730-1.
  * [20] Q. Li, C. Tai, and W. E.  Stochastic modified equations and dynamics of stochastic gradient algorithms I: Mathematical foundations.  Journal of Machine Learning Research, 20:1–47, 2019.  URL: http://jmlr.org/papers/v20/17-526.html.
  * [21] F. Liu, X. Huang, Y. Chen, and A.K. Suykens.  Random features for kernel approximation: A survey on algorithms, theory and beyond, 2020.  arXiv:2004.11154.
  * [22] A. Nadeem and R. Potthast.  Transformed and generalized localization for ensemble methods in data assimilation.  Mathematical Methods in the Applied Sciences, 39:619–634, 2016\.  doi:10.1002/mma.3496.
  * [23] M. Opper.  A Bayesian approach to online learning.  In DavidEditor Saad, editor, Online Learning in Neural Networks, pages 363–378. Cambridge University Press, 1999.  doi:10.1017/CBO9780511569920.017.
  * [24] G.A. Pavliotis, A.M. Stuart, and U. Vaes.  Derivative-free Bayesian inversion using multiscale dynamics, 2021.  arXiv:2102.00540.
  * [25] M. Rabovsky and J.L. McClelland.  Quasi-compositional mapping from form to meaning: A neural network based approach to capturing neural responses during human language comprehension.  Phil. Trans. Royal Soc. B: Biological Sciences, 375:20190313, 2019\.  doi:10.1098/rstb.2019.0313.
  * [26] S. Reich.  A dynamical systems framework for intermittent data assimilation.  BIT Numer. Math., 51:235–249, 2011.  doi:10.1007/s10543-010-0302-4.
  * [27] S. Reich.  A nonparametric ensemble transform method for Bayesian inference.  SIAM J. Sci. Comput., 35:A2013–A2024, 2013.  doi:10.1137/130907367.
  * [28] S. Reich and C.J. Cotter.  Ensemble filter techniques for intermittent data assimilation.  In M. Cullen, Freitag M. A., S. Kindermann, and R. Scheichl, editors, Large Scale Inverse Problems. Computational Methods and Applications in the Earth Sciences, volume 13 of Radon Ser. Comput. Appl. Math., pages 91–134. Walter de Gruyter, Berlin, 2013.  doi:10.1515/9783110282269.
  * [29] S. Reich and C.J. Cotter.  Probabilistic forecasting and Bayesian data assimilation.  Cambridge University Press, New York, 2015.  doi:10.1017/CBO9781107706804.
  * [30] S. Reich and S. Weissmann.  Fokker–Planck particle systems for Bayesian inference: Computational approaches.  SIAM/ASA J. Uncertainty Quantification, 9:446–482, 2021.  doi:10.1137/19M1303162.
  * [31] P. Sakov, D.S. Oliver, and L. Bertino.  An iterative EnKF for strongly nonlinear systems.  Monthly Weather Review, 140:1988–2004, 2012.  doi:10.1175/MWR-D-11-00176.1.
  * [32] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(56):1929–1958, 2014.  URL: http://jmlr.org/papers/v15/srivastava14a.html.
  * [33] T.J. Sullivan.  Introduction to Uncertainty Quantification, volume 63 of Texts in Applied Mathematics.  Springer, 2015.  doi:10.1007/978-3-319-23395-6.
  * [34] A. Taghvaei, J. de Wiljes, P.G. Mehta, and S. Reich.  Kalman filter and its modern extensions for the continuous-time nonlinear filtering problem.  ASME. J. Dyn. Sys., Meas., Control., 140:030904–030904–11, 2017\.  doi:10.1115/1.4037780.
  * [35] A. Taghvaei, P.G. Mehta, and S.P. Meyn.  Diffusion map-based algorithm for gain function approximation in the feedback particle filter.  SIAM/ASA J. Uncertain. Quantif., 8(3):1090–1117, 2020.  doi:10.1137/19M124513X.
  * [36] J. Tödter, P. Kirchgessner, L. Nerger, and B. Ahrens.  Assessment of a nonlinear ensemble transform filter for high-dimensional data assimilation.  Monthly Weather Review, 144:409–427, 2016.  doi:10.1175/MWR-D-15-0073.1.
  * [37] C.L. Wormell and S. Reich.  Spectral convergence of diffusion maps: Improved error bounds and an alternative normalisation.  SIAM J. Numer. Anal., 59:1687–1734, 2021.  doi:10.1137/20M1344093.
  * [38] Tao Yang, Henk A P Blom, and Prashant G Mehta.  The Continuous-Discrete Time Feedback Particle Filter.  In 2014 American Control Conference, pages 648–653. American Automatic Control Council, 2014.  doi:10.1109/ACC.2014.6859259.
  * [39] Tao Yang, Prashant G. Mehta, and Sean P. Meyn.  Feedback particle filter.  IEEE Trans. Automat. Control, 58:2465–2480, 2013.  doi:10.1109/TAC.2013.2258825.

## Appendix

We provide an analysis of the SDE

(56) ${\rm
d}\theta_{\tau}=-\Sigma_{\tau}\left\\{\nabla_{\theta}V_{\tau}(\theta_{\tau})\,{\rm
d}\tau-(\Delta\tau\Omega_{\tau})^{1/2}{\rm d}W_{\tau}\right\\},$

which is obtained from the modified SDE (21) by setting
$\Omega_{\tau}(\theta)=\Omega_{\tau}$, that is, by ignoring the state
dependence of the diffusion term, and using the EnKBF drift term from linear
regression [3], that is,

$\nabla_{\theta}V_{\tau}(\theta)=\frac{1}{2}G^{\rm
T}\Gamma^{-1}(G\theta+Gm_{\tau}-2t).$

Conditioned on $W_{\tau}$, the evolution equation for the (conditional) mean
becomes

(57) ${\rm d}m_{\tau}=-\Sigma_{\tau}\left\\{G^{\rm
T}\Gamma^{-1}(Gm_{\tau}-t){\rm d}\tau-(\Delta\tau\Omega_{\tau})^{1/2}{\rm
d}W_{\tau}\right\\}$

and that for the (conditional) covariance matrix

(58) $\frac{{\rm d}}{{\rm d}\tau}\Sigma_{\tau}=-\Sigma_{\tau}G^{\rm
T}\Gamma^{-1}G\Sigma_{\tau}\,.$

We can further decompose (57) into an equation for the mean, which we denote
by $\mu_{\tau}$, and the covariance $P_{\tau}$. We obtain the two evolution
equations

(59) $\frac{\rm d}{{\rm d}\tau}\mu_{\tau}=-\Sigma_{\tau}G^{\rm
T}\Gamma^{-1}(G\mu_{\tau}-t)$

and

(60) $\frac{{\rm d}}{{\rm d}\tau}P_{\tau}=-\Sigma_{\tau}G^{\rm
T}\Gamma^{-1}GP_{\tau}-P_{\tau}G^{\rm
T}\Gamma^{-1}G\Sigma_{\tau}+\Delta\tau\Sigma_{\tau}\Omega_{\tau}\Sigma_{\tau},$

respectively. The initial conditions are $\mu_{0}=m_{0}$ and $P_{0}=0$.

We note that (59) and (58) are equivalent to the equations arising from the
standard Kalman–Bucy filter with solutions

$\mu_{\tau}=\mu_{0}-K_{\tau}(G\mu_{0}-t)$

and

$\Sigma_{\tau}=\Sigma_{0}-K_{\tau}G\Sigma_{0},$

respectively, where $K_{\tau}$ denotes the Kalman gain matrix

$K_{\tau}=\tau\Sigma_{0}G^{\rm T}(\Gamma+\tau G\Sigma_{0}G^{\rm T})^{-1}.$

If we now consider the fully observed case, that is $G=I$, under small
measurement error covariance $\Gamma=\epsilon I$, $\epsilon\ll 1$, then we
find that

$\mu_{1}\approx\theta_{\rm MLE}=t,\qquad\Sigma_{1}\approx\Gamma=\epsilon
I,\qquad
P_{1}<\frac{\Delta\tau}{2}\Sigma_{1}\Omega_{1}\Sigma_{1}\approx\epsilon^{2}\frac{\Delta\tau}{2}\Omega_{1}.$

If we define

$\Omega_{\tau}=\pi_{\tau}[\Omega_{\tau}(\theta)]$

for simplicity, where $\pi_{\tau}$ is the exact filtering distribution for the
Kalman–Bucy filter, then $\Omega_{0}=\mathcal{O}(\epsilon^{-2})$ while
$\Omega_{1}=\mathcal{O}(\epsilon^{-1})$. Hence, the variance, $P_{\tau}$, of
the posterior conditional mean, $m_{\tau}$, is much smaller than the posterior
(as well as the frequentist) variance, $\Sigma_{\tau}$, at $\tau=1$ provided
the step-size $\Delta\tau$ in (56) is chosen small enough.