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.
Similar content being viewed by others
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) [35, 38, 39], have almost exclusively focused on the squared error loss function of the form
with \(g:{\mathbb {R}}^D \rightarrow {\mathbb {R}}^N\) a forward map, \(t\in {\mathbb {R}}^N\) the data, \(\varGamma \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 [10, 11, 19] 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
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,9,10, 17], 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 Sect. 2. The homotopy approach to Bayesian inference is summarised in Sect. 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 Sect. 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 Sect. 7 followed by a summary statement in Sect. 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 Sect. 7.
The posterior distribution for class \(C_1\) is defined as a logistic sigmoid
acting on a linear combination of x-dependent and vector-valued features \(\phi _x \in {\mathbb {R}}^D\) so that
The model has D adjustable parameters \(\theta \in {\mathbb {R}}^{D}\). The probability of the complimentary 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
where \(y_n := y_{x_n}(\theta )\) with
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
and, using
obtain the gradient of the cross-entropy function
We rewrite this gradient more compactly as as
using
and
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
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 \(\varPsi _{\mathrm{data}}(\theta )\) by \(\theta _{\mathrm{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
with
that is, the posterior parameter distribution \(\pi _{\mathrm{post}}\) is given by
Here the prior mean \(m_{\mathrm{prior}} \in {\mathbb {R}}^D\) and covariance matrix \(\varSigma _{\mathrm{prior}} \in {\mathbb {R}}^{D\times D}\) need to be chosen appropriately.
It is known that \(\varPsi _{\mathrm{post}}\) is a strongly convex function since \(D^2_\theta \varPsi _{\mathrm{data}}(\theta ) \ge 0\) and, hence, the MAP estimator, that is,
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 _{\mathrm{MAP}}\) and covariance matrix
where
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 _{\mathrm{post}}\), we focus in this paper on ensemble transform methods based on the homotopy approach [6, 26]
with \(\tau \in [0,1]\), that is, \(\pi _0 = \pi _{\mathrm{prior}}\) and \(\pi _{\mathrm{post}} = \pi _1\). The PDFs \(\pi _\tau \) satisfy the evolution equation
where \(\pi _\tau [f]\) denotes the expectation value of a function \(f:{\mathbb {R}}^D \rightarrow {\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 \rightarrow {\mathbb {R}}\), which are defined by the elliptic partial differential equation (PDE)
where \(\varSigma _\tau \) denotes the covariance matrix of \(\pi _\tau \), that is,
Hence the evolution equations are given by
such that
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
\(\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 \(\varPsi _{\mathrm{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 [15, 17] 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
will lead to an associated evolution equation
such that (11) holds. The specific choice \(F_\tau = \varSigma _\tau \nabla V_\tau \) in (9) is motivated by the concept of affine-invariance [10, 11], which we explain next.
Let
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 \(\varPsi _{\mathrm{data}}\) and \(\varPsi _{\mathrm{prior}}\) and the corresponding problem for the transformed variable \({\overline{\theta }}\) by
Hence the PDFs \({\overline{\pi }}_\tau \) for \({\overline{\theta }}\) satisfy
Definition 1
A homotopy formulation (12) is called affine-invariant if the associated vector fields \(F_\tau \) and \({\overline{F}}_\tau \) satisfy
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
leads to an affine-invariant homotopy formulation and the transformed vector field is given by
with the potential \({\overline{V}}_\tau \) satisfying the elliptic PDE
Proof
Since the covariance matrices satisfy
the PDE (9) transforms into the PDE (16) under the affine transformation (13). Hence we have
and obtain
which verifies (14). \(\square \)
We note that the affine invariance of (10) is maintained under the forward Euler time discretisation
with step-size \(\varDelta \tau >0\) and \(\tau _k = k\varDelta \tau \), \(k = 0,\ldots ,K-1\) such that \(K\varDelta \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 \(\varSigma _\tau \nabla _\theta U_\tau (\theta )\) with the potential \(U_\tau \) chosen such that
A possible choice is
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 ) = \text{ const. }\) and (17) follows trivially. Upon formally replacing the drift term \(\varSigma _\tau \nabla _\theta U_\tau (\theta )\) by its EKS/ALDI stochastic mean-field representation, (10) turns into the affine-invariant stochastic evolution equation
where \(W_\tau \) denotes standard D-dimensional Brownian motion. The associated marginal PDF \(\pi _\tau \) still satisfies (8) provided \(\pi _0 = \pi _{\mathrm{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 \(\varDelta \) denote a randomly selected subset of \(\{1,\ldots ,N\}\) without replacement of cardinality \(N' \le N\). Given such a mini-batch \(\varDelta \), the full gradient (5) is replaced by its approximation
More abstractly, data sub-sampling in terms of mini-batches gives rise to random negative log-likelihood functions \(\varPsi _{\mathrm{data}}^\gamma \), where \(\gamma \) characterises the chosen mini-batch \(\varDelta ^\gamma \) and is the realisation of an appropriate random variable \(\varGamma \) such that
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 \(\varPsi _{\mathrm{data}}^\gamma \) is now being used in (9) instead of \(\varPsi _{\mathrm{data}}\), then we denote the resulting potential by \(V_\tau ^\gamma \), which is also unbiased, that is,
due to (19). Let us now time-step the implied mean-field ODEs (10) by the forward Euler method, that is,
where the \(\gamma _k\)’s are independent realisations of the random variable \(\varGamma \). A formal application of the modified equation analysis of [20] leads to the stochastic differential equation (SDE)
where \(W_\tau \) denotes standard D-dimensional Brownian motion and
Here \(\text{ cov}_\varGamma [f^\gamma ]\) denotes the covariance matrix of a vector-valued random function \(f^\gamma \in {\mathbb {R}}^D\) with respect to the random variable \(\varGamma \). 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 \(\varDelta \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 Sect. 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 _{\mathrm{prior}}\). The three methods considered here differ in the evolution equations for the ensemble \(\{\theta _\tau ^{(i)}\}\). We introduce the empirical mean
and the empirical covariance matrix
for \(\tau \ge 0\). We also introduce the matrix \(\varTheta _\tau \in {\mathbb {R}}^{D \times M}\) of ensemble deviations
which leads to the compact representation
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,
More generally,
for an observable \(f(\theta )\). We describe an affine-invariant extension of the FPF next. The key building stone 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
with covariance matrix
and (9) can be rewritten as
Following the FPF methodology, as carefully described in [34], we regularise the elliptic PDE problem (9) through the fixed-point equation
with regularisation parameter \(\epsilon >0\) and semigroup \(P_\epsilon = e^{\epsilon \varDelta _{\pi _\tau }}\). One then replaces \(V_\tau \) by \({{\tilde{V}}}_\tau \) in the evolution Eq. (10).
We note that \(P_\epsilon \) is the semi-group associated with the mean-field SDE
where \(W_\epsilon \), \(\epsilon \ge 0\), is standard Brownian motion in \({\mathbb {R}}^D\). The SDE (23) 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:
where \(n_\epsilon (\theta ) := \int k_\epsilon (\theta ,\theta ')\pi _\tau (\theta '){\mathrm{d}}\theta '\) is the normalisation constant,
and
is the Gaussian kernel associated with the pure diffusion in (23).
The next step is to replace \(\pi _\tau \) by the ensemble-induced empirical measure \({{\hat{\pi }}}_\tau \), which results in the empirical approximation
where \(n_\epsilon ^{(M)}(\theta ) := \sum _i k_\epsilon ^{(M)} (\theta ,\theta ^{(i)}_\tau )\) is the normalisation constant, and
Here \(g^{(M)}_\epsilon \) denotes the Gaussian kernel (24) with the covariance matrix \(\varSigma _\tau \) replaced by its empirical estimate \({{\hat{\varSigma }}}_\tau \). One finally introduces the \(M\times M\) Markov matrix \({\mathrm{T}}\) with entries
and the finite-dimensional fixed-point equation
with
Given the solution vector \(\tilde{\mathrm{V}} \in {\mathbb {R}}^M\) with entries \(\tilde{\mathrm{V}}^{(j)}\), \(j=1,\ldots ,M\), the approximation to the drift term in the ODE (10) is now provided as follows [34]:
Lemma 2
The affine-invariant FPF drift term for the ensemble members \(\theta _\tau ^{(i)}\), \(i=1,\ldots ,M\), is given by
with coefficients
Proof
The calculations leading to (26) are identical to those for the standard FPF [34] taking note that the covariance matrix \({\hat{\varSigma }}_\tau \) cancels out. \(\square \)
The resulting interacting particle system can be propagated in time using the forward Euler method with step size \(\varDelta \tau = 1/K\), that is,
\(k=0,\ldots ,K-1\), with the right-hand side defined by (26). Also note that \(M\rightarrow \infty \) leads formally to the mean-field forward Euler method (20).
Remark 4
We note that the Markov matrix \({\mathrm{T}}\), defined by (25), is non-reversible contrary to the underlying stochastic process (23). 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\rightarrow \infty \) and \(\epsilon \rightarrow 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 \(\varSigma _\tau \) of \(\pi _\tau \) from Eq. (8). Furthermore, if one assumes that \(\pi _\tau \) is well approximated by a Gaussian PDF with mean \(m_\tau \) and covariance matrix \(\varSigma _\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 \(\varSigma _\tau \), then their evolution equations are given by
and
respectively.
Proof
The evolution equations are derived using the following well-known identity for Gaussian PDFs \(\pi _\tau \) (see, for example, [23]):
as well as
for any scalar-valued smooth function \(f:{\mathbb {R}}^D \rightarrow {\mathbb {R}}\). \(\square \)
In order to derive corresponding equations for the ensemble \(\{\theta _\tau ^{(i)}\}\), we introduce
as well as
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
We finally obtain the evolution equations
for the ensemble mean and
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
and
one formally obtains the approximation
to the second-order Eqs. (30)–(31).
Lemma 4
The EnKBF equations (32) are of affine-invariant gradient structure
with potential
Proof
The lemma can be verified by direct calculation and taking note of
in particular. \(\square \)
The Eq. (33) also possesses an affine-invariant gradient-flow structure in the space of probability measures in the mean-field limit \(M\rightarrow \infty \). See [9, 10, 28] for details.
Remark 5
The authors of [17] proposed a different modification of the EnKF for classification problems. In our notation, it corresponds to
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 \(\varPsi _{\mathrm{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 \(\varTheta _\tau \) satisfy the evolution Eq. (31). The same propagation of the first and second-order moments is achieved by the stochastic equations
in the mean-field limit \(M\rightarrow \infty \), where \(W_\tau \) denotes standard N-dimensional Brownian motion. This suggests to replace the finite ensemble size formulation (32) by the interacting system of SDEs
Here \(W_\tau ^{(i)}\), \(i=1,\ldots ,M\), denote independent N-dimensional Brownian motions.
Again, the forward Euler method (27) can be used to solve (32) or (33) 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{\varSigma }}}_\tau \) for \(M\le 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 \le D\), then the empirical covariance matrix \({{\hat{\varSigma }}}_\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{\varTheta }}}_\tau \). The empirical covariance matrix is now defined by
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 (34) underestimates the cross-correlations. An unbiased estimator requires an additional rescaling of the off-diagonal entries by \(1/(1-\mu )\).
Remark 7
There are alternative methods for breaking the sub-space property of ensemble transform methods in case of \(M\le 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 (27). 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 (27) for the EnKBF formulation from Sect. 4.3. Starting from (32) we employ
with step size \(\varDelta \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 (30)–(31).
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 \rightarrow \infty \). Let us therefore turn the transport-based differential equation (10) into the McKean–Vlasov SDE
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 \varPsi _{\mathrm{prior}}\) of the prior in the drift term while (36) contains
instead, and (ii) in a similar vain, the potential \(\varPsi _{\mathrm{data}}(\theta )\) is replaced in (36) 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 (36) and (ii) that the SDE (36) 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 \(\varPsi _{\mathrm{data}}\) is a potential of the form
then the marginal PDFs \(\pi _\tau \) of (36) remain Gaussian for all \(\tau >0\), and the SDE (36) has a stationary measure which is given by \(\pi _\infty (\theta ) \propto \exp (-\varPsi _{\mathrm{post}}(\theta ))\).
Proof
Under the stated assumption, we have
as well as
Equating both expressions leads us to
The evolution equations arising from (36) for the mean \(m_\tau \) and the covariance matrix \(\varSigma _\tau \) are therefore given by
and
respectively. The associated equilibrium solutions are
and
as desired. We also note that the linearity of the gradient (37) implies that the PDFs \(\pi _\tau \) remain Gaussian. Finally, (37) is indeed the gradient of the potential
which concludes the proof. \(\square \)
Remark 8
One can replace the McKean–Vlasov SDE (36) by the following formulation
where the potential \(U_\tau \) now satisfies the elliptic PDE
instead of (9). Both formulations are equivalent for Gaussian PDFs \(\pi _\tau \). While (36) is closer to the homotopy approach of Sect. 4, formulation (38) has a closer resemblance to overdamped Langevin dynamics.
5.1 Asymptotic Behaviour
The case of a general negative log-likelihood function \(\varPsi _{\mathrm{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,
then following the arguments from Sect. 4.2, the first two moments satisfy
and
The equation for the covariance matrix is stable if \({{\tilde{\pi }}}_\tau [D_\theta ^2 \varPsi _{\mathrm{post}} ]\) is positive-definite for all \(\tau \ge 0\). This is, for example, the case for logistic regression.
We now discuss the limiting \(\tau \rightarrow \infty \) behavior of (36) in some more detail. Its equilibrium distributions are characterised by the following lemma.
Lemma 6
Stationary measures \(\pi _\infty \) of the mean-field equation (36) satisfy the PDE
or, in terms of expectation values with respect to test functions \(\phi \),
Proof
The lemma follows from the Fokker–Planck equation associated to (36) and the PDE (9). \(\square \)
Again making the Gaussian approximation (39), formally setting \(\tau = \infty \) and recalling Remark 8, we first obtain from (42) the identity
Then,
for \(\phi (\theta ) = \theta -m_\infty \), and
for \(\phi (\theta ) = (\theta -m_\infty )(\theta -m_\infty )^{\mathrm{T}}\), which provide a self-consistent system of equations for the mean \(m_\infty \) and the covariance matrix \(\varSigma _\infty \) of the Gaussian approximation (39). Here we have again used (28) and (29), respectively.
Note that (43) and (44) are compatible with (40) and (41), respectively. However, (43) and (44) 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 (36) or (38), 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 (36) utilising methods from sequential data assimilation, which are well suited to deal with the \(-\varSigma _\tau \nabla _\theta V_\tau \) term in the SDE (36). In other words, we numerically solve the SDE (36) with a step size \(\varDelta \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
and (ii) the SDE
over the same time-interval \(\varDelta \tau \). More precisely, given \(\theta _{\tau _k} \sim \pi _{\tau _k}\) we first find
followed by solving (45) with initial \({\tilde{\theta }}_{\tau _{k+1}}\) over the time-interval \(\varDelta \tau \) in order to obtain \(\theta _{\tau _{k+1}} \sim \pi _{\tau _{k+1}}\) for \(k=0,\ldots ,K-1\). Note that (45) can be solved robustly by the tamed forward Euler-Maruyama method
with \(\varXi _{k+1}\) a D-dimensional standard Gaussian random variable. As for ALDI, the \(\varSigma _\tau ^{1/2}\) multiplying the Brownian motion in (45) 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 (46) in terms of ensembles \(\{\theta _{\tau _k}^{(i)}\}\) rely on ensemble transform filters of the form [29]
for \(j =1,\ldots ,M\) with appropriate coefficients \(s_{ij}\), which satisfy
We collect the coefficients \(s_{ij}\) in the \(M\times M\) matrix S.
Lemma 7
An ensemble transform filter is affine-invariant (compare Sect. 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
at \(\tau _k\), then
and we conclude that
at \(\tau _{k+1}\) for all \(i=1,\ldots ,M\). \(\square \)
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
where \(w_{\tau _k} \in {\mathbb {R}}^M \) is the vector of normalised importance weights with entries
\(\mathrm{1}_M \in {\mathbb {R}}^M\) is the vector of ones, and \(\text{ 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].
Remark 9
The ensemble transform particle filter (ETPF) of [27, 29] can also be made affine-invariant by using the cost function
subject to the constraints \(s_{ij}\ge 0\), (48), and
in the optimal transport definition of the coefficients \(s_{ij}\) in (47). 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.
Example 1
We follow the example from [5] Section 4.2.1. Specifically, consider the Gaussian likelihood functions
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^{\mathrm{T}},1)^{\mathrm{T}}\), that is \(D = J+1\), and the true parameter value is
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
and
respectively. See Section 4.2.2 in [5] for details. In this example, we will however directly estimate the true parameter value (51) using Bayesian inference.
We consider the special case \(J = 2\) and use \(B = \mathrm{I}\), \(\mathrm{I}\) the identity matrix, \(\mu _1 = (-1,-1)^{\mathrm{T}}\), \(\mu _2 = (2,2)^{\mathrm{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
We now test the proposed algorithms using (i) an informative Gaussian prior with mean (52) and covariance matrix \(\varSigma _{\mathrm{prior}} = \mathrm{I}\), and (ii) a less informative Gaussian prior with mean \(m_{\mathrm{prior}} = 0\) and covariance matrix \(\varSigma _{\mathrm{prior}} = 4\,\mathrm{I}\). We time step the affine-invariant FPF, the second-order filter (30)–(31), and the EnKBF (32) with the forward Euler method using a step size \(\varDelta \tau = 10^{-3}\). A further reduction of the step size did not change the results in a statistically significant manner. The McKean–Vlasov SDE (36) 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 (38) were run up to time \(\tau = 10\) with step-size \(\varDelta \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 (36), explicitly involve the evaluation of the negative log-likelihood function (3). We found that we needed to replace (4) by
in order to avoid numerical instabilities due to exceedingly small values of \(\ln y_n\). The band-width, \(\epsilon \), in (24) 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 performance 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 (36) 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 (52) 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\).
We now consider a high-dimensional extension of the previous example, for which we further investigate the performance of the EnKBF (33) in terms of the tamed time-stepping method (35) for ensemble sizes \(M\le D\), localisation via dropouts, and data sub-sampling.
Example 2
We consider the simple feature map
that is \(J=D\) such that
for \(D=50\). The data is generated by first drawing a \(\theta = \theta _{\mathrm{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 _{\mathrm{ref}})\), and \(t_n = 0\) otherwise. We use the EnKBF tamed time-stepping method (35) with step size \(\varDelta \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_{\mathrm{prior}}=0\) and the covariance matrix \(\varSigma _{\mathrm{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{\varSigma }}}_{\tau _n}\) following (34). 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 _{\mathrm{ref}}\) and the ensemble mean \({{\hat{m}}}_\tau \) at final time \(\tau = 1\). We also state the spectral norm of \({{\hat{\varSigma }}}_\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' = 100\). We also compare the results to those from the standard EnKBF formulation (32). 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 \(\varDelta \tau \) in (35). However changing the threshold value \(\mu = 0.5\) in the dropout localisation (34) 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.
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 Sect. 2 to the more general logistic regression problem with \(y_x(\theta ) := \sigma (f_x(\theta ))\) for given functions \(f_x:{\mathbb {R}}^D \rightarrow {\mathbb {R}}\). Hence the gradient \(\nabla _\theta y_x(\theta )\) becomes
where the previously considered feature maps \(\phi _x\) now also depend on the parameters \(\theta \), that is,
The gradient of \(\varPsi _{\mathrm{data}}\) is thus explicitly given by
and the Hessian matrix of second-order derivatives becomes
We now discuss a derivative-free implementation of the EnKBF (32) for nonlinear logistic regression. The key difference to linear logistic regression is that the matrix \(\varPhi \in {\mathbb {R}}^{D\times N}\) now depends on \(\theta \) and is given by
In order to avoid the computation of the gradients of \(f_x\), one replaces \(\varSigma _\tau \varPhi (\theta _\tau )\) in (32) by the empirical correlation matrix
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 [13, 17] 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].
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
where the activation is given by
The adjustable parameters are now \(\theta = (\theta _1^{\mathrm{T}},\ldots ,\theta _L^{\mathrm{T}})^{\mathrm{T}} \in {\mathbb {R}}^{L D}\). The negative log-likelihood function becomes
with
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
\(l=1,\ldots ,L\). The Hessian matrix of second-order derivatives got \(D\times D\) block entries
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.
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}}\rightarrow {\mathbb {R}}^+\) which we represent using a sigmoidal random feature model, that is,
where \(\lambda ^*>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
Following the notation of Sect. 3.2 and using I uniformly distributed samples \({{\hat{x}}}_i\), \(i=1,\ldots ,I\), on \({\mathcal {X}}\), we obtain the random unbiased estimator
for this intractable likelihood \(\varPsi _{\mathrm{data}}\). Here we also used the previously introduced notation (4), that is, \(\lambda _x(\theta ) = \lambda ^*y_x(\theta )\).
We now apply the homotopy approach (10) to the associated Bayesian inference problem. Following the discussion in Sect. 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 (53) instead of \(\varPsi _{\mathrm{data}}\) in (9) and takes the limit \(\varDelta \tau \rightarrow 0\).
Example 3
Let us briefly discuss a related problem in a strictly Gaussian setting. We start from the intractable negative log-likelihood function
Its gradient is provided by
and the associated mean-field EnKBF becomes
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 \mathrm{U}[0,1]\), \(i=1,\ldots ,I\), and introduce a time step \(\varDelta \tau \) leading to
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 \(\varDelta \tau \rightarrow 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 \(\varDelta \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 (36) and (38), 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.
References
Adams, R., Murray, I., MacKay, D.: Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. In: Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, pp. 9–16. Association for Computing Machinery, New York, NY, USA (2009). https://doi.org/10.1145/1553374.1553376
Amezcua, J., Kalnay, E., Ide, K., Reich, S.: Ensemble transform Kalman–Bucy filters. Q.J.R. Meteor. Soc. 140, 995–1004 (2014). https://doi.org/10.1002/qj.2186
Bergemann, K., Reich, S.: A localization technique for ensemble Kalman filters. Quarterly Journal of the Royal Meteorological Society 136, 701–707 (2010). https://doi.org/10.1002/qj.591
Beskos, A., Crisan, D., Jasra, A.: On the stability of sequential Monte Carlo methods in high dimensions. The Annals of Applied Probability 24, 1396–1445 (2014). https://doi.org/10.1214/13-AAP951
Bishop, C.M.: Pattern Recognition and Machine Learning. Springer, New York (2006)
Daum, F., Huang, J., Noushin, A.: Exact particle flow for nonlinear filters. In: I. Kadar (ed.) Signal Processing, Sensor Fusion, and Target Recognition XIX, vol. 7697, pp. 92–110. International Society for Optics and Photonics, SPIE (2010). https://doi.org/10.1117/12.839590
Donner, C., Opper, M.: 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
Evensen, G.: Data Assimilation: The Ensemble Kalman Filter, second edn. Springer, Berlin (2009). https://doi.org/10.1007/978-3-642-03711-5
Garbuno-Inigo, A., Hoffmann, F., Li, W., Stuart, A.: Interacting Langevin diffusions: Gradient structure and ensemble Kalman sampler. SIAM J. Appl. Dyn. Syst. 19(1), 412–441 (2020). https://doi.org/10.1137/19M1251655
Garbuno-Inigo, A., Nüsken, N., Reich, S.: Affine invariant interacting Langevin dynamics for Bayesian inference. SIAM J. Appl. Dyn. Syst. 19, 1633–1658 (2020). https://doi.org/10.1137/19M1304891
Goodman, J., Weare, J.: Ensemble samplers with affine invariance. Comm. Appl. Math. and Comput. Science 5, 65–80 (2010). https://doi.org/10.2140/camcos.2010.5.65
Gottwald, G., Reich, S.: Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation. Physica D 423, 132911 (2021). https://doi.org/10.1016/j.physd.2021.132911
Haber, E., Lucka, F., Ruthotto, L.: Never look back—A modified EnKF method and its application to the training of neutral networks without back propagation (2018). arXiv:1805.08034
Heber, F., Ž. Trst’anová, Leimkuhler, B.: 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
Huang, D., Schneider, T., Stuart, A.: Uncented Kalman inversion (2021). arXiv:2102.01580
Kass, R., Raftery, A.: Bayes factors. Journal of the American Statistical Association 90, 773–795 (1995). https://doi.org/10.1080/01621459.1995.10476572
Kovachki, N., Stuart, A.: Ensemble Kalman inversion: A derivative-free technique for machine learning tasks. Inverse Probl. 35(9), 095005 (2019). https://doi.org/10.1088/1361-6420/ab1c3a
Lee, J., Xiao, L., Schoenholz, S., Bahri, Y., Novak, R., Sohl-Dickstein, J., Pennington, J.: Wide neural networks of any depth evolve as linear models under gradient descent. Journal of Statistical Mechanics: Theory and Experiment 2020(12), 124002 (2020). https://doi.org/10.1088/1742-5468/abc62b
Leimkuhler, B., Matthews, C., Weare, J.: Ensemble preconditioning for Markov chain Monte Carlo simulation. Statistics and Computing 28, 277–290 (2018). https://doi.org/10.1007/s11222-017-9730-1
Li, Q., Tai, C., E, W.: 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
Liu, F., Huang, X., Chen, Y., Suykens, A.: Random features for kernel approximation: A survey on algorithms, theory and beyond (2020). arXiv:2004.11154
Nadeem, A., Potthast, R.: Transformed and generalized localization for ensemble methods in data assimilation. Mathematical Methods in the Applied Sciences 39, 619–634 (2016). https://doi.org/10.1002/mma.3496
Opper, M.: A Bayesian approach to online learning. In: D. Saad (ed.) Online Learning in Neural Networks, pp. 363–378. Cambridge University Press (1999). https://doi.org/10.1017/CBO9780511569920.017
Pavliotis, G., Stuart, A., Vaes, U.: Derivative-free Bayesian inversion using multiscale dynamics (2021). arXiv:2102.00540
Rabovsky, M., McClelland, J.: 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). https://doi.org/10.1098/rstb.2019.0313
Reich, S.: A dynamical systems framework for intermittent data assimilation. BIT Numer. Math. 51, 235–249 (2011). https://doi.org/10.1007/s10543-010-0302-4
Reich, S.: A nonparametric ensemble transform method for Bayesian inference. SIAM J. Sci. Comput. 35, A2013–A2024 (2013). https://doi.org/10.1137/130907367
Reich, S., Cotter, C.: Ensemble filter techniques for intermittent data assimilation. In: M. Cullen, F.M. A., S. Kindermann, R. Scheichl (eds.) Large Scale Inverse Problems. Computational Methods and Applications in the Earth Sciences, Radon Ser. Comput. Appl. Math., vol. 13, pp. 91–134. Walter de Gruyter, Berlin (2013). https://doi.org/10.1515/9783110282269
Reich, S., Cotter, C.: Probabilistic forecasting and Bayesian data assimilation. Cambridge University Press, New York (2015). https://doi.org/10.1017/CBO9781107706804
Reich, S., Weissmann, S.: Fokker–Planck particle systems for Bayesian inference: Computational approaches. SIAM/ASA J. Uncertainty Quantification 9, 446–482 (2021). https://doi.org/10.1137/19M1303162
Sakov, P., Oliver, D., Bertino, L.: An iterative EnKF for strongly nonlinear systems. Monthly Weather Review 140, 1988–2004 (2012). https://doi.org/10.1175/MWR-D-11-00176.1
Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: 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
Sullivan, T.: Introduction to Uncertainty Quantification, Texts in Applied Mathematics, vol. 63. Springer (2015). https://doi.org/10.1007/978-3-319-23395-6
Taghvaei, A., Mehta, P., Meyn, S.: Diffusion map-based algorithm for gain function approximation in the feedback particle filter. SIAM/ASA J. Uncertain. Quantif. 8(3), 1090–1117 (2020). https://doi.org/10.1137/19M124513X
Taghvaei, A., de Wiljes, J., Mehta, P., Reich, S.: Kalman filter and its modern extensions for the continuous-time nonlinear filtering problem. ASME. J. Dyn. Sys., Meas., Control. 140, 030904–030904–11 (2017). https://doi.org/10.1115/1.4037780
Tödter, J., Kirchgessner, P., Nerger, L., Ahrens, B.: Assessment of a nonlinear ensemble transform filter for high-dimensional data assimilation. Monthly Weather Review 144, 409–427 (2016). https://doi.org/10.1175/MWR-D-15-0073.1
Wormell, C., Reich, S.: Spectral convergence of diffusion maps: Improved error bounds and an alternative normalisation. SIAM J. Numer. Anal. 59, 1687–1734 (2021). https://doi.org/10.1137/20M1344093
Yang, T., Blom, H.A.P., Mehta, P.G.: The Continuous-Discrete Time Feedback Particle Filter. In: 2014 American Control Conference, pp. 648–653. American Automatic Control Council (2014). https://doi.org/10.1109/ACC.2014.6859259
Yang, T., Mehta, P.G., Meyn, S.P.: Feedback particle filter. IEEE Trans. Automat. Control 58, 2465–2480 (2013). https://doi.org/10.1109/TAC.2013.2258825
Open Access
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Communicated by Teresa Krick and Hans Munthe-Kaas.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Invited paper based on the FoCM 2021 Online Seminar lecture Statistical inverse problems and affine-invariant gradient flow structures in the space of probability measures presented by Sebastian Reich in June 2021.
SR has been partially funded by Deutsche Forschungsgemeinschaft (DFG) - Project-ID 318763901 - SFB1294..
Appendix
Appendix
We provide an analysis of the SDE
which is obtained from the modified SDE (21) by setting \(\varOmega _\tau (\theta ) = \varOmega _\tau \), that is, by ignoring the state dependence of the diffusion term, and using the EnKBF drift term from linear regression [3], that is,
Conditioned on \(W_\tau \), the evolution equation for the (conditional) mean becomes
and that for the (conditional) covariance matrix
We can further decompose (55) into an equation for the mean, which we denote by \(\mu _\tau \), and the covariance \(P_\tau \). We obtain the two evolution equations
and
respectively. The initial conditions are \(\mu _0 = m_0\) and \(P_0 = 0\).
We note that (57) and (56) are equivalent to the equations arising from the standard Kalman–Bucy filter with solutions
and
respectively, where \(K_\tau \) denotes the Kalman gain matrix
If we now consider the fully observed case, that is \(G = I\), under small measurement error covariance \(\varGamma = \epsilon I\), \(\epsilon \ll 1\), then we find that
If we define
for simplicity, where \(\pi _\tau \) is the exact filtering distribution for the Kalman–Bucy filter, then \(\varOmega _0 = {\mathcal {O}}(\epsilon ^{-2})\) while \(\varOmega _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, \(\varSigma _\tau \), at \(\tau = 1\) provided the step-size \(\varDelta \tau \) in (54) is chosen small enough.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Pidstrigach, J., Reich, S. Affine-Invariant Ensemble Transform Methods for Logistic Regression. Found Comput Math 23, 675–708 (2023). https://doi.org/10.1007/s10208-022-09550-2
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10208-022-09550-2
Keywords
- Logistic regression
- Bayesian inference
- Interacting particle systems
- Affine invariance
- Ensemble Kalman filter
- Langevin dynamics