0.0.0.1 Keywords.

Ensemble Kalman inversion; second-order particle methods; inverse problems; derivative-free optimization; interacting particle systems.

0.0.0.2 MSC2020.

Primary 65J22; Secondary 65C35, 37N40, 34D20, 65K10.

1 Introduction↩︎

Inverse problems arise in many areas of applied mathematics, science and engineering, where one aims to recover unknown parameters or states from indirect and typically noisy observations. Classical approaches to inverse problems rely on regularization techniques, such as Tikhonov-type methods, iterative regularization, and variational formulations; see, for instance, [1]–[4]. In the Bayesian framework, the solution of the inverse problem is interpreted as a posterior probability distribution, combining the information contained in the data with prior knowledge on the unknown [3], [4].

Ensemble Kalman methods provide a derivative-free computational framework for data assimilation and inverse problems. Starting from the Ensemble Kalman Filter introduced in [5], ensemble-based Kalman methodologies have been widely used in geophysical applications, reservoir engineering and uncertainty quantification; see, among others, [6]–[8]. In the context of inverse problems, Ensemble Kalman Inversion (EKI) was introduced and analyzed as an iterative ensemble method in [9]. Its key feature is the use of empirical covariances computed from an ensemble of particles to define a Kalman-type update, thereby avoiding the explicit computation of derivatives of the forward map. This makes EKI particularly attractive for large-scale or black-box inverse problems.

The mathematical analysis of EKI has developed in several directions. Convergence properties and regularization effects have been studied in the linear and noisy settings in [10]–[12], while Tikhonov regularization within EKI has been investigated in [13]. Further variants include constrained formulations [14]–[16], regularized ensemble Kalman methods [17], hierarchical approaches [18], and applications to machine learning tasks [19], [20]. Mean-field, continuous-time and spectral formulations have also played an important role in clarifying the dynamical structure of the method [21]–[26].

From a dynamical systems viewpoint, continuous-time EKI can be interpreted as an interacting particle system driven by a covariance-preconditioned descent direction. In the linear inverse problem setting, this structure reveals both the strengths and the limitations of the method. On the one hand, the empirical covariance adapts the dynamics to the geometry of the ensemble. On the other hand, the covariance may collapse prematurely, so that the covariance matrix loses rank and the dynamics becomes confined to a possibly non-optimal subspace. This mechanism makes the method sensitive to the initial ensemble and may lead to stationary configurations which do not correspond to the desired solution. Stabilization mechanisms for continuous-time EKI have therefore been proposed using covariance inflation and relaxation [27].

Recent work has also explored acceleration mechanisms for EKI, inspired by classical inertial and accelerated optimization methods such as Polyak’s heavy-ball method [28] and Nesterov’s accelerated gradient scheme [29]. The connection between accelerated optimization methods and continuous-time second-order dynamics has been clarified, for instance, in [30]. In [31], Nesterov acceleration is incorporated into EKI and related Kalman inversion variants by means of a particle-level nudging step.

The particle interpretation of EKI connects naturally with interacting particle methods for optimization, such as Particle Swarm Optimization [32]–[34] and Consensus-Based Optimization [35], [36]. Related developments include stochastic and mean-field connections between PSO and CBO [37], memory and random selection mechanisms [38], polarized and mirror-type consensus dynamics [39], [40], and micro–macro decompositions of PSO [41].

In this paper we introduce an inertial Ensemble Kalman particle system for inverse and optimization problems. The method extends continuous-time EKI by assigning positions and velocities to the particles and by coupling the Kalman-type force with damping, attraction towards the ensemble mean, and short-range repulsion. Unlike Nesterov-type accelerations, which modify the discrete ensemble update, the proposed dynamics is a continuous-time second-order interacting particle system. The repulsive interaction is designed to destabilize fully collapsed configurations, while the attraction term prevents uncontrolled dispersion, thereby providing a mechanism to preserve ensemble diversity without abandoning the covariance-based structure of EKI. For linear inverse problems, we derive the mean and fluctuation dynamics, prove linear instability of fully collapsed configurations in a suitable parameter regime, characterize asymptotic equilibria through optimality on the subspace retained by the limiting covariance, and analyze the associated frozen-covariance mean dynamics. Numerical experiments on linear, nonconvex, and nonlinear inverse problems show improved robustness in adverse initialization regimes where first-order EKI-type methods may collapse prematurely.

2 From continuous-time EKI to an inertial particle system↩︎

2.1 Continuous-time EKI↩︎

We recall the continuous-time formulation of EKI and the main structural features that motivate the second-order dynamics introduced below. Let \[y = G(u^\dagger) + \eta\] be the observed data, where \(u^\dagger \in \mathbb{R}^d\) is the unknown parameter, \(G:\mathbb{R}^d \to \mathbb{R}^K\) is the continuous forward map, and \(\eta\) denotes observational noise. We assume that \(\eta\sim\mathcal{N}(0,\Gamma)\), where \(\Gamma \in \mathbb{R}^{K\times K}\) is a symmetric positive definite covariance matrix, and we consider the least-squares functional \[\Phi(u) = \frac{1}{2} \|y-G(u)\|_{\Gamma}^{2}. \label{eq:least-squares-functional}\tag{1}\] Here and in the following, \(\|z\|_\Gamma^2 := z^T\Gamma^{-1}z\).

Let \(U(t)=\{u_j(t)\}_{j=1}^J\) be an ensemble of \(J\) particles in \(\mathbb{R}^d\). We denote by \[\bar u(t) = \frac{1}{J}\sum_{j=1}^J u_j(t), \qquad \bar G(t) = \frac{1}{J}\sum_{j=1}^J G(u_j(t))\] the empirical means in parameter and observation space. The empirical cross-covariance between parameters and model outputs is defined by \[C^{uG}(U(t)) = \frac{1}{J}\sum_{j=1}^J (u_j(t)-\bar u(t))\otimes (G(u_j(t))-\bar G(t)) \in \mathbb{R}^{d\times K}. \label{eq:cross-covariance}\tag{2}\] The continuous-time EKI dynamics [10] is then given by \[\dot{u}_j(t) = C^{uG}(U(t))\Gamma^{-1}(y-G(u_j(t))), \qquad j=1,\dots,J. \label{eq:continuous-eki}\tag{3}\] This system can be interpreted as an interacting particle method in which the drift of each particle is determined by empirical covariance information extracted from the ensemble.

In the linear case \(G(u)=Gu\), with \(G\in\mathbb{R}^{K\times d}\), the least-squares functional satisfies \[\nabla\Phi(u)=G^T\Gamma^{-1}(Gu-y)=Au-b, \qquad A:=G^T\Gamma^{-1}G,\quad b:=G^T\Gamma^{-1}y . \label{eq:Phi-A-b-definitions}\tag{4}\] Moreover, \[C^{uG}(U(t))=C(U(t))G^T,\] where \[C(U(t)) = \frac{1}{J}\sum_{j=1}^J (u_j(t)-\bar u(t))\otimes(u_j(t)-\bar u(t))\] is the empirical covariance. Hence the continuous-time EKI dynamics reduces to \[\dot{u}_j(t) = -C(U(t))(Au_j(t)-b), \qquad j=1,\dots,J. \label{eq:linear-continuous-eki}\tag{5}\] Thus, in the linear setting, EKI is a covariance-preconditioned gradient flow. The dynamics is still nonlinear, since the preconditioner \(C(U(t))\) depends on the evolving ensemble.

In the linear case, averaging 5 gives the mean dynamics \[\dot{\bar u}(t) = -C(U(t))(A\bar u(t)-b). \label{eq:first-order-mean-dynamics}\tag{6}\] If \(e_j(t)=u_j(t)-\bar u(t)\), then the fluctuations satisfy \[\dot{e}_j(t) = -C(U(t))A e_j(t), \qquad j=1,\dots,J. \label{eq:first-order-fluctuation-dynamics}\tag{7}\] These identities show that the evolution of the mean is entirely mediated by the empirical covariance, while the fluctuations determine the directions retained by the ensemble.

A fundamental structural property of EKI is the subspace property: the dynamics remains confined to the affine space generated by the initial ensemble, \[u_j(t)\in \bar u(0)+\operatorname{span}\{e_k(0):k=1,\dots,J\}, \qquad j=1,\dots,J. \label{eq:subspace-property}\tag{8}\] Moreover, if \(C(U_*)=0\), then all particles coincide and the right-hand side of 5 vanishes. Thus every fully collapsed ensemble is a stationary configuration, independently of whether its common value minimizes \(\Phi\). More generally, if the covariance loses rank, the dynamics can only move along the directions retained by the ensemble covariance. At a stationary configuration \(U_*\), 6 gives \[C(U_*)(A\bar u_*-b)=0,\] so that the gradient of \(\Phi\) is required to vanish only on the range of the empirical covariance. This covariance-collapse mechanism motivates the inertial interacting particle system introduced below.

2.2 The inertial particle system↩︎

In the previous section the EKI particles were denoted by \(u_j(t)\), consistently with the inverse problem variable \(u\). In the second-order formulation it is convenient to use a mechanical notation. We therefore denote by \[X(t)=\{x_j(t)\}_{j=1}^J, \qquad V(t)=\{v_j(t)\}_{j=1}^J\] the particle positions and velocities, respectively. The position \(x_j(t)\in\mathbb{R}^d\) plays the role of the parameter particle \(u_j(t)\) in first-order EKI.

We introduce a second-order particle system of the form \[\begin{cases} \dot{x}_j(t)=v_j(t),\\[0.4em] \dot{v}_j(t) = a_j^{\rm EKI}(X(t),V(t)) + a_j^{\rm rep}(X(t)) + a_j^{\rm att}(X(t)), \end{cases} \qquad j=1,\dots,J. \label{eq:second-order-abstract}\tag{9}\] The three terms in the acceleration represent, respectively, a Kalman-type relaxation force, a short-range repulsive interaction, and an attraction force towards the ensemble mean.

As before, we define \[\bar x(t)=\frac{1}{J}\sum_{j=1}^J x_j(t), \qquad \bar G(t)=\frac{1}{J}\sum_{j=1}^J G(x_j(t)),\] and the empirical cross-covariance \[C^{xG}(X(t)) = \frac{1}{J}\sum_{j=1}^J (x_j(t)-\bar x(t))\otimes (G(x_j(t))-\bar G(t)). \label{eq:second-order-cross-covariance}\tag{10}\] The Kalman-type relaxation force is defined by \[a_j^{\rm EKI}(X(t),V(t)) = -\gamma v_j(t) + \beta C^{xG}(X(t))\Gamma^{-1}(y-G(x_j(t))), \label{eq:second-order-eki-force}\tag{11}\] where \(\gamma>0\) is a damping coefficient and \(\beta>0\) controls the strength of the EKI correction. The damping term dissipates kinetic energy, while the second term is the natural Kalman-type force inherited from continuous-time EKI.

The repulsive interaction is introduced to counteract premature collapse of the ensemble. We define \[a_j^{\rm rep}(X(t)) = k\sum_{i\neq j} f(\|x_i(t)-x_j(t)\|) (x_j(t)-x_i(t)), \label{eq:second-order-repulsion}\tag{12}\] where \(k\geq0\) is the repulsion strength and \(f:[0,\infty)\to(0,\infty)\) is a smooth non-increasing interaction kernel. In the numerical experiments and in the stability discussion below, we use the regularized inverse-power kernel \[f(r)=\frac{1}{(\varepsilon+r)^p}, \qquad \varepsilon>0,\quad p>1. \label{eq:repulsive-kernel}\tag{13}\] This choice produces a strong short-range repulsion while remaining bounded at \(r=0\). The parameter \(\varepsilon\) controls the regularization near the origin, whereas \(p\) controls the decay of the interaction at large distances.

The attraction force towards the ensemble mean is defined by \[a_j^{\rm att}(X(t)) = -\alpha(x_j(t)-\bar x(t)), \label{eq:second-order-attraction}\tag{14}\] with \(\alpha\geq0\). Its role is to prevent the ensemble from dispersing excessively. The proposed dynamics is therefore based on a competition between attraction and repulsion: the attraction term promotes collective coherence, while the repulsive term acts against complete collapse.

Combining 9 –14 , we obtain the second-order EKI system \[\begin{cases} \dot{x}_j(t)=v_j(t),\\[0.4em] \begin{align} \dot{v}_j(t) &= -\gamma v_j(t) + \beta C^{xG}(X(t))\Gamma^{-1}(y-G(x_j(t))) \\ &\quad + k\displaystyle\sum_{i\neq j} f(\|x_i(t)-x_j(t)\|) (x_j(t)-x_i(t)) - \alpha(x_j(t)-\bar x(t)), \end{align} \end{cases} \label{eq:second-order-eki}\tag{15}\] for \(j=1,\dots,J\).

The model can be interpreted as an inertial, heavy-ball-type reformulation of continuous-time EKI. In the overdamped regime \(\gamma\gg1\), with \(\beta/\gamma=\mathcal{O}(1)\), the velocity equation formally relaxes to a first-order Kalman-type drift. At the same time, the additional attraction–repulsion mechanism changes the internal geometry of the ensemble and is designed to reduce the tendency of first-order EKI to collapse prematurely.

The analysis below focuses on the linear case for which the Kalman force reduces to a covariance-preconditioned gradient term.

3 Analysis in the linear inverse problem setting↩︎

We analyze the inertial EKI system introduced in Section 2.2 in the linear inverse problem setting. Thus \(G(x)=Gx\), and \[\nabla\Phi(x)=Ax-b, \qquad A=G^T\Gamma^{-1}G,\quad b=G^T\Gamma^{-1}y.\] Since \(C^{xG}(X)=C(X)G^T\), the particle positions satisfy \[\begin{align} \ddot x_j(t) + \gamma \dot{x}_j(t) &= -\beta C(X(t))(Ax_j(t)-b) \\ &\quad + k\sum_{i\neq j} f(\|x_i(t)-x_j(t)\|) (x_j(t)-x_i(t)) \\ &\quad - \alpha(x_j(t)-\bar x(t)), \end{align} \label{eq:second-order-position-form}\tag{16}\] for \(j=1,\dots,J\). It is a covariance-preconditioned heavy-ball dynamics coupled with attraction–repulsion interactions within the ensemble. This is the form used throughout the analysis. The associated velocities are \(v_j=\dot{x}_j\).

We first record the equations satisfied by the ensemble mean and the fluctuations. Let \[\tilde{x}_j(t):=x_j(t)-\bar x(t), \qquad \frac{1}{J}\sum_{j=1}^J \tilde{x}_j(t)=0.\]

The mean dynamics 17 is an inertial EKI equation. The internal ensemble geometry is controlled by 18 , where Kalman contraction, mean attraction, and pairwise repulsion compete. This separation is one of the useful structural properties of the model.

For the regularized inverse-power kernel 13 the interaction term is nonsingular at particle collisions and satisfies the local Lipschitz assumption in Proposition 1. Thus local well-posedness holds without any additional structural assumption.

3.1 Instability of collapsed configurations and non-collapsed equilibria↩︎

We now study the local stability of fully collapsed configurations. These configurations are stationary states of the particle system. If all particles have the same position and zero velocity, the empirical covariance vanishes and both the attraction and repulsion terms are zero. The following result shows that, for a sufficiently strong repulsive interaction, such configurations are unstable with respect to zero-mean perturbations.

We relate the instability of collapsed configurations to the possible asymptotic states of the ensemble.

Recall that an ensemble configuration \(X\) is collapsed if \(C(X)=0\).

The ensemble spread is given by \[S(X(t)):=\frac{1}{J}\sum_{j=1}^J\|\tilde{x}_j(t)\|^2 =\operatorname{Tr}(C(X(t))),\] and \[\begin{align} \frac{\mathrm d^2}{\mathrm dt^2}S(X(t)) + \gamma\frac{\mathrm d}{\mathrm dt}S(X(t)) &= \frac{2}{J}\sum_{j=1}^J \|\dot{\tilde{x}}_j(t)\|^2 - 2\beta\operatorname{Tr}\!\left(C(X(t))AC(X(t))\right) \notag\\ &\quad - 2\alpha S(X(t)) + \frac{k}{J} \sum_{i\neq j} f(r_{ij}(t))r_{ij}(t)^2, \label{eq:spread-dynamics} \end{align}\tag{20}\] where \(r_{ij}(t)=\|\tilde{x}_i(t)-\tilde{x}_j(t)\|\). This identity shows the competition between the repulsive contribution, which increases the spread, and the attraction and Kalman terms, that tend to decrease it.

3.2 Limiting mean dynamics: optimality and frozen-covariance decay↩︎

We now discuss the behavior of the ensemble mean in the limiting covariance regime. First, we characterize the optimality condition satisfied by asymptotic equilibria. Then we consider the frozen-covariance dynamics, which provides a linear asymptotic model along the retained directions.

The limiting covariance may be nonzero but still rank-deficient. Therefore, the equilibrium condition for the ensemble mean does not necessarily imply full optimality in \(\mathbb{R}^d\). It only yields optimality along the directions retained by the limiting ensemble covariance.

We conclude the analysis with a decay estimate for the mean dynamics in the frozen-covariance regime.

Let \[\mathcal{V}_\infty:=\operatorname{Range}(C_\infty)\] be the subspace retained by the limiting covariance, and let \(x^\dagger\) be a minimizer of \(\Phi\) on the affine space associated with \(\mathcal{V}_\infty\). Equivalently, \(x^\dagger\) satisfies \[C_\infty(Ax^\dagger-b)=0. \label{eq:frozen-retained-minimizer}\tag{21}\] The frozen-covariance mean dynamics is \[\ddot{\bar x}(t) + \gamma \dot{\bar x}(t) = -\beta C_\infty(A\bar x(t)-b). \label{eq:frozen-mean-dynamics}\tag{22}\]

4 Numerical experiments↩︎

In this section we investigate the behavior of the proposed second-order EKI dynamics on a set of representative test problems. We illustrate the qualitative effect of inertia, damping, attraction and repulsion on the ensemble dynamics. In particular, we monitor the evolution of the data misfit, the error of the ensemble mean, and the ensemble spread, comparing the second-order dynamics with first-order EKI-type methods.

All experiments are performed using a time-discrete approximation of the particle system. We first describe the discretization used for the second-order model, emphasizing the implicit-explicit treatment of the damping term.

4.1 Time discretization↩︎

We discretize the second-order EKI system 15 by an implicit-explicit scheme [42]. The position equation is advanced explicitly, while the velocity equation is treated implicitly only in the linear damping term. Given \(X^n=\{x_j^n\}_{j=1}^J\) and \(V^n=\{v_j^n\}_{j=1}^J\), we set \[x_j^{n+1}=x_j^n+\Delta t\,v_j^n, \qquad j=1,\dots,J. \label{eq:numerics-position-update}\tag{25}\] Then the velocity is updated by \[\begin{align} v_j^{n+1} &= v_j^n + \Delta t \Big[ -\gamma v_j^{n+1} + \beta C^{xG}(X^{n+1})\Gamma^{-1} \bigl(y-G(x_j^{n+1})\bigr) \Big] \notag\\ &\quad + \Delta t \Bigg[ k\sum_{i\neq j} f(\|x_i^n-x_j^n\|) (x_j^n-x_i^n) - \alpha(x_j^n-\bar x^n) \Bigg]. \label{eq:numerics-velocity-imex-implicit} \end{align}\tag{26}\] Solving explicitly for \(v_j^{n+1}\), we obtain \[\begin{align} v_j^{n+1} &= \frac{1}{1+\gamma\Delta t}v_j^n + \frac{\beta\Delta t}{1+\gamma\Delta t} C^{xG}(X^{n+1})\Gamma^{-1} \bigl(y-G(x_j^{n+1})\bigr) \notag\\ &\quad + \frac{k\Delta t}{1+\gamma\Delta t} \sum_{i\neq j} f(\|x_i^n-x_j^n\|) (x_j^n-x_i^n) - \frac{\alpha\Delta t}{1+\gamma\Delta t} (x_j^n-\bar x^n). \label{eq:numerics-velocity-imex} \end{align}\tag{27}\] The scheme is explicit with respect to the particle interactions and does not require the solution of nonlinear systems. The only implicit contribution is the damping term, whose effect is accounted for through the scalar factor \((1+\gamma\Delta t)^{-1}\).

In the overdamped regime \(\gamma\gg1\), with \(\beta/\gamma=\mathcal{O}(1)\), the update 27 formally relaxes to a first-order Kalman-type correction evaluated at the updated positions. This provides a simple and robust discretization of the inertial dynamics while preserving the computational structure of ensemble-based methods.

In the following experiments, we use the implicit–explicit discretization 25 –27 for the second-order dynamics, whereas the first-order EKI-type dynamics are advanced by an explicit Euler scheme.

4.2 Diagnostics and stopping criteria↩︎

In all numerical experiments we monitor a common set of diagnostic quantities in order to compare the behavior of the different ensemble dynamics. Given an ensemble \[X^n=\{x_j^n\}_{j=1}^J, \qquad \bar x^n=\frac{1}{J}\sum_{j=1}^J x_j^n,\] we evaluate the weighted data misfit at the ensemble mean, \[\Phi(\bar x^n) = \frac{1}{2} \bigl(y-G(\bar x^n)\bigr)^T \Gamma^{-1} \bigl(y-G(\bar x^n)\bigr). \label{eq:numerical-misfit}\tag{28}\] When the true parameter \(x^\dagger\) is known, we also compute the relative error of the ensemble mean, \[e_x^n = \frac{\|\bar x^n-x^\dagger\|}{\|x^\dagger\|}. \label{eq:numerical-mean-error}\tag{29}\] To measure ensemble diversity, we use the normalized empirical spread \[S^n = \left( \frac{1}{Jd} \sum_{j=1}^J \|x_j^n-\bar x^n\|^2 \right)^{1/2}. \label{eq:numerical-spread}\tag{30}\]

For inverse problems with additive Gaussian observational noise \(\eta\sim\mathcal{N}(0,\Gamma)\), we use a discrepancy-type stopping criterion based on the weighted residual. More precisely, if \(K\) denotes the dimension of the data space, we stop the iteration when \[\Phi(\bar x^n) \leq \frac{\tau K}{2}, \label{eq:discrepancy-stopping}\tag{31}\] where \(\tau\geq1\) is a prescribed tolerance. This choice is motivated by the fact that, after whitening by \(\Gamma^{-1/2}\), the noise has approximately \(K\) independent unit-variance components, so that the expected squared weighted noise level is of order \(K\).

In practice, the discrepancy principle is combined with a maximum number of iterations, or equivalently a final integration time. In the experiments below we also impose a minimum number of iterations before activating 31 , in order to avoid premature stopping during the initial transient. If the discrepancy criterion is not met, the final output is taken at the final iteration. For multi-seed experiments, we report mean values and standard deviations of the relevant diagnostics.

Some test problems involve additional problem-specific quantities. In particular, for the Darcy flow inverse problem we also report the relative error on the reconstructed permeability field, \[e_a^n = \frac{\|a_{\rm rec}^n-a^\dagger\|}{\|a^\dagger\|}, \label{eq:numerical-field-error}\tag{32}\] where \(a^\dagger\) is the discretized true permeability field and \(a_{\rm rec}^n\) is the permeability field reconstructed from the ensemble mean at iteration \(n\).

4.3 A linear elliptic inverse problem↩︎

We first consider a linear inverse problem associated with a one-dimensional elliptic equation. This test is useful because it is consistent with the linear setting analyzed in Section 3, while still retaining the ill-posed character of inverse problems.

Let \(p\) solve \[-p''(s)+p(s)=u(s), \qquad s\in(0,\pi), \qquad p(0)=p(\pi)=0. \label{eq:elliptic-1d-model}\tag{33}\] The inverse problem consists in recovering the forcing term \(u\) from noisy observations of the state \(p\). After discretization on \(d\) interior grid points, the forward map is linear and can be written as \[p = G u,\] where \(G\in\mathbb{R}^{d\times d}\) is the inverse of the finite-difference discretization of the operator \(-\mathrm d^2/\mathrm ds^2+I\). Thus both the unknown forcing \(u\) and the observed state \(p\) are represented by vectors in \(\mathbb{R}^{d}\). Synthetic data are generated as \[y = G u^\dagger+\eta, \qquad \eta\sim\mathcal{N}(0,\sigma^2 I),\] with \(\sigma = 0.1\) and \[u^\dagger(s)=10\sin(8s).\] In the following numerical simulations, we consider \(d=K=63\).

We compare three ensemble methods: standard first-order EKI, Nesterov-accelerated EKI, and the proposed second-order EKI. The Nesterov-accelerated method is implemented according to the particle-level nudging strategy of [31]. The first-order method and the Nesterov-accelerated EKI are advanced explicitly, whereas the second-order dynamics is discretized by the implicit–explicit scheme described in Section 4.1.

The initial ensemble is chosen to be both concentrated and biased away from the truth. More precisely, \[x_j^0\sim \mathcal{N}(-4\cdot\mathbf{1},\sigma_0^2 I), \qquad \sigma_0=10^{-2}, \qquad j=1,\dots,J,\] and the initial velocities of the second-order method are set to zero. This initialization is deliberately unfavorable for standard EKI, since the initial ensemble has small spread and explores a region far from the true parameter. All ensemble methods are run with \(J=80\) particles.

The parameters used in the second-order method are \[\gamma=2,\qquad \beta=0.5,\qquad \alpha=0.4,\qquad k=0.1,\qquad \varepsilon=10^{-3},\qquad p=1.5.\] Final time is chosen as \(T=1000\) and the time step as \(\Delta t = 10^{-1}\). All methods are stopped according to the discrepancy principle 31 , with the additional safeguards described in Section 4.2. In particular, we take \(\tau = 1.05\) and the minimum number of iterations before stopping is \(20\).

Figure 1 reports the evolution of the weighted data misfit, the parameter relative error at the ensemble mean, and the ensemble spread. The Nesterov-accelerated method substantially improves the decay of the data misfit with respect to standard EKI. However, it also accelerates the collapse of the ensemble, as shown by the rapid decay of the spread. This behavior is consistent with the fact that Nesterov acceleration modifies the first-order update but does not introduce an explicit mechanism to preserve ensemble diversity.

By contrast, the second-order method reaches the discrepancy threshold much earlier but maintains a significantly larger spread throughout the computation. This allows to achieve the smallest parameter error in this test. This suggests that the improvement is not merely due to acceleration, but rather to the modified internal ensemble dynamics induced by the damping and attraction–repulsion terms. In particular, the repulsive interaction prevents the ensemble from collapsing too early, while the damping controls the inertial motion. In conclusion, acceleration alone does not prevent covariance collapse.

Figure 2 shows the reconstruction of the control \(u\) and of the observed state \(p=Gu\). Both Nesterov-accelerated EKI and the second-order method provide a more accurate reconstruction of the observed state and follows the data more closely than standard EKI. The recovered control exhibits more oscillations, reflecting the ill-posedness of the inverse map and the absence of an explicit regularization term. This behavior is expected: a better fit in the observation space does not necessarily imply a uniformly smoother reconstruction of the parameter.

Table 1 summarizes the final values of the main diagnostics. In addition to the final misfit, parameter relative error and spread, we report the time at which the discrepancy criterion is first reached and the time step used by each method.

To check that the behavior observed above is not an artifact of a single initial ensemble, we repeated the experiment for different values of the initial ensemble spread \(\sigma_0\). The results in Table 2 show that the second-order method remains competitive across the tested regimes and gives the smallest relative parameter error for the most concentrated and the most dispersed initial ensembles. The Nesterov-accelerated method reaches the discrepancy threshold faster, but this is not systematically reflected in a smaller parameter error, consistently with the ill-posedness of the inverse map.

We observed the same qualitative behavior for structured low-frequency initial ensembles, obtained by sampling random combinations of the first sine modes. This indicates that the improved performance of the second-order method is not merely an artifact of the full-dimensional Gaussian initialization. These additional tests are not reported here for brevity, but they are available in Appendix 6.

4.4 Nonconvex Ackley benchmark↩︎

We next consider the Ackley function as a nonconvex benchmark to illustrate the behavior of the different ensemble dynamics in the presence of multiple local minima. In order to keep the test within the inverse-problem formulation, we set \[G(x)=\operatorname{Ack}(x), \qquad x\in\mathbb{R}^d, \qquad y=0\in\mathbb{R},\] so that the least-squares functional becomes \[\Phi(x) = \frac{1}{2} \operatorname{Ack}(x)^2.\] The global minimizer is \(x^\dagger=0\in\mathbb{R}^d\), and no observational noise is added.

Since this benchmark is noise-free, no discrepancy principle is used. Instead, the iterations are stopped when either \[\Phi(\bar x)\leq 10^{-8} \qquad\text{or}\qquad \|\bar x\|\leq 10^{-3},\] or when the final integration time is reached.

We use the one-dimensional Ackley function to compare standard EKI, Nesterov-accelerated EKI, and the proposed second-order dynamics. We then use the two-dimensional Ackley function to isolate the effect of the repulsive interaction by comparing two second-order simulations with weak and strong repulsion.

4.4.1 One-dimensional test↩︎

We first consider the one-dimensional Ackley landscape. The initial ensemble is concentrated around a local basin far from the global minimizer, \[x_j^0 \sim \mathcal{N}(5,\,0.1^2), \qquad j=1,\dots,J,\] and the initial velocities of the second-order method are set to zero. All methods are run up to the final time \(T=20\) with time step \(\Delta t=10^{-2}\) and \(J=80\) particles. For the second-order dynamics we use \[\gamma=2,\qquad \beta=1,\qquad \alpha=0.4,\qquad k=0.05,\qquad \varepsilon=10^{-3},\qquad p=1.5.\]

Figure 3 reports the evolution of the misfit and of the ensemble spread. Standard EKI and Nesterov-accelerated EKI remain trapped near the initial basin and show almost no effective decrease of the misfit. Their ensemble spread also remains small, indicating that the ensemble does not explore enough of the nonconvex landscape. By contrast, the second-order method rapidly increases the ensemble spread and then drives the ensemble mean toward the global minimizer.

This behavior is confirmed in Figure 4, where the final particles are displayed on the Ackley landscape. The first-order methods remain concentrated near the initial local basin, whereas the second-order dynamics produces a wider ensemble distribution and reaches the global minimizer at the level of the ensemble mean. In this test, the second-order method stops at time \(t=16.86\), with \[\Phi(\bar x)=8.18\cdot 10^{-6}, \qquad |\bar x|=9.98\cdot 10^{-4},\] while the two first-order methods remain close to \(x\simeq 4.96\) until the final time. See Table 3.

4.4.2 Two-dimensional test↩︎

We now use the two-dimensional Ackley function to isolate the effect of the repulsive interaction in the second-order dynamics. The initial ensemble is concentrated near a local basin away from the global minimum: \[x_j^0 \sim \mathcal{N}\bigl((3.5,3.5)^T,\,0.15^2 I\bigr), \qquad j=1,\dots,J,\] with \(J=100\), and the initial velocities are set to zero. The final time is \(T=30\) and the time step is \(\Delta t=10^{-2}\).

We compare two second-order simulations with the same parameters \[\gamma=2,\qquad \beta=1,\qquad \alpha=0.5,\qquad \varepsilon=10^{-2},\qquad p=1.5,\] but different repulsion strengths. In the weakly repulsive regime we set \[k=10^{-6},\] whereas in the strongly repulsive regime we set \[k=0.05.\]

Figure 5 displays the final ensembles on the Ackley landscape, together with the trajectory of the ensemble mean. In the weakly repulsive case, the ensemble collapses near a non-optimal local basin and the mean remains far from the global minimizer. This test illustrates the role of the attraction–repulsion mechanism: attraction alone tends to concentrate the ensemble, while sufficiently strong repulsion prevents premature collapse and allows the inertial dynamics to explore the nonconvex landscape.

Figure 6 reports the corresponding time histories of the misfit and ensemble spread. In the weakly repulsive case, the spread rapidly decays to zero and the misfit remains large. In the strongly repulsive case, the spread is maintained at a positive level and the misfit decreases by several orders of magnitude. The final values for weak repulsion are \[\Phi(\bar x)=4.13\cdot 10^{1}, \qquad \|\bar x\|=4.27, \qquad S=1.08\cdot 10^{-3}.\] In contrast, for stronger repulsion the ensemble remains non-collapsed and the mean reaches the global minimizer with high accuracy: \[\Phi(\bar x)=3.00\cdot 10^{-4}, \qquad \|\bar x\|=8.06\cdot 10^{-3}, \qquad S=2.19.\] Table 3 summarizes the final diagnostics.

4.5 Darcy flow inverse problem↩︎

We finally consider a nonlinear inverse problem governed by a Darcy-type elliptic equation. This test is more challenging than the previous ones because the forward map is nonlinear and the unknown parameter controls the permeability field of the PDE.

Let \(D=(0,1)^2\). We consider \[-\nabla\cdot\bigl(a(s)\nabla u(s)\bigr)=f(s), \qquad s\in D, \label{eq:darcy-forward}\tag{34}\] with homogeneous Dirichlet boundary condition \[u(s)=0, \qquad s\in\partial D.\] Here \(u\) denotes the hydraulic head, \(a\) is the permeability field, and we take \[f(x)\equiv 1.\]

The permeability is parameterized on a \(3\times3\) partition of the domain. More precisely, \[a(s) = \sum_{\ell=1}^9 e^{\theta_\ell}\chi_{D_\ell}(s), \qquad \theta=(\theta_1,\dots,\theta_9)^T\in\mathbb{R}^9, \label{eq:darcy-parametrization}\tag{35}\] where \(\{D_\ell\}_{\ell=1}^9\) are the nine subdomains, ordered row-wise from top-left to bottom-right. The exponential parameterization guarantees positivity of the permeability.

The reference parameter used to generate synthetic data is \[\theta^\dagger = (-0.4,\;0.8,\;1.1,\;-1.0,\;0.2,\;0.7,\;-0.6,\;1.0,\;0.1)^T. \label{eq:darcy-theta-true}\tag{36}\]

We observe the solution at nine interior sensors located on the \(3\times3\) grid \[(0.25,0.25),\;(0.25,0.50),\;\dots,\;(0.75,0.75).\] Let \[\mathcal{G}(\theta)\in\mathbb{R}^K, \qquad K=9,\] denote the forward observation map obtained by solving 34 and extracting the solution values at the sensor locations. Synthetic data are generated as \[y=\mathcal{G}(\theta^\dagger)+\eta, \qquad \eta\sim\mathcal{N}(0,\Gamma), \qquad \Gamma=\sigma^2 I_K, \qquad \sigma=0.02. \label{eq:darcy-data}\tag{37}\] The inverse problem consists in recovering \(\theta\) from the noisy data \(y\).

The initial ensemble is chosen to be both concentrated and biased away from the truth: \[\theta_j^0 \sim \mathcal{N}(\mu_0,\sigma_0^2 I), \qquad \mu_0=-1.5\cdot\mathbf{1}, \qquad \sigma_0=0.1, \qquad j=1,\dots,J. \label{eq:darcy-initialization}\tag{38}\] This unfavorable initialization is designed to test whether the second-order dynamics can prevent premature ensemble collapse and recover useful directions even when the initial cloud is poorly informative.

All methods use the diagnostic quantities and stopping rules described in Section 4.2. In addition to the parameter error and the weighted misfit, we also report the relative permeability reconstruction error 32 .

4.5.1 Parameter sensitivity↩︎

Before comparing with standard EKI, we perform a short parameter sensitivity study for the second-order method in the initialization regime 38 . We vary the damping coefficient \(\gamma\), the strength of the data-driven Kalman force \(\beta\), and the repulsion strength \(k\), while keeping \[\alpha=0.2, \qquad \varepsilon=10^{-3}, \qquad p=1.5\] fixed. The tested grid is \[\gamma\in\{1,2\}, \qquad \beta\in\{0.5,1\}, \qquad k\in\{0.01,0.05,0.1\}. \label{eq:darcy-tuning-grid}\tag{39}\]

The purpose of this tuning step is not to identify a universally optimal parameter set, but to select a robust configuration for the difficult initialization regime considered here. The results indicate that a moderate Kalman forcing, stronger damping, and stronger repulsion are preferable in this setting. In particular, the choice \[\gamma=2, \qquad \beta=0.5, \qquad k=0.1 \label{eq:darcy-selected-parameters}\tag{40}\] provides the best compromise among data fit, parameter reconstruction, permeability reconstruction, and robustness across seeds. The complete tuning results are reported in Appendix 7.

4.5.2 Comparison with standard EKI↩︎

We now compare standard EKI with the second-order method using the parameter configuration 40 . Both methods are run with the same noisy observations, the same ensemble size, and the same multi-seed protocol.

Figure 7 reports the mean time evolution, over \(10\) seeds, of the weighted data misfit, the parameter error, and the ensemble spread. Standard EKI rapidly reduces the initial spread and then stagnates with a relatively large residual. By contrast, the second-order method maintains a substantially larger ensemble spread and achieves a much smaller data misfit. This behavior is consistent with the mechanism observed in the previous tests: the attraction–repulsion dynamics prevents premature collapse and allows the ensemble to keep exploring relevant directions.

To complement the mean trajectories, Figure 8 shows the distributions of the final parameter error, weighted misfit, and relative permeability error. The corresponding aggregate statistics are reported in Table 4.

The second-order method improves all the reported final metrics. It achieves a smaller parameter error, a significantly smaller weighted data misfit, and a smaller relative permeability reconstruction error. At the same time, it maintains a much larger ensemble spread. In this experiment, the larger spread is not a drawback: it is precisely the mechanism that prevents the ensemble from collapsing too early near the unfavorable initial configuration. The standard first-order dynamics, instead, quickly loses diversity and fails to recover from the poor initialization.

Overall, the Darcy experiment confirms the main numerical result. In a nonlinear inverse problem with a concentrated and biased initial ensemble, the second-order dynamics is more robust than standard EKI because it combines the Kalman-type descent mechanism with an internal attraction–repulsion dynamics that preserves useful ensemble diversity.

5 Conclusions↩︎

We introduced a second-order Ensemble Kalman particle system for inverse and optimization problems. The proposed dynamics augments the continuous-time EKI formulation with explicit velocities, damping, attraction towards the ensemble mean, and short-range repulsion between particles. The resulting model can be interpreted as a heavy-ball-type reformulation of EKI in which the internal geometry of the ensemble is directly controlled by competing attractive and repulsive mechanisms.

For linear inverse problems, we analyzed the mean and fluctuation dynamics of the second-order system. The analysis shows that fully collapsed configurations may be destabilized by a sufficiently strong repulsive interaction. We also characterized the optimality condition satisfied by asymptotic equilibria. Finally, we showed exponential decay along the retained directions under a natural positivity assumption.

The numerical experiments confirm the structural role of the second-order terms. The method is not uniformly faster, nor universally superior to first-order EKI-type methods. Its main advantage appears in adverse initialization regimes, where standard EKI and accelerated variants may collapse prematurely. In such cases, the attraction–repulsion mechanism maintains a nonzero ensemble spread and allows the particles to continue exploring informative directions. This behavior is visible in the linear elliptic inverse problem and in the Darcy flow inverse problem.

Acknowledgments↩︎

This work was carried out within the activities of the PRIN PNRR Project 2022 No. P2022JC95T, “Data-driven discovery and control of multiscale interacting artificial agent systems”, funded by MUR (Ministry of University and Research) and Next Generation EU – European Commission.

G.V. acknowledges the support of Sapienza University under Ateneo Project 2024 “Advanced Computational Methods for Real-World Applications: Data-Driven Models, Hyperbolic Equations and Optimal Control”.

G.V. is a member of the INdAM Research National Group of Scientific Computing (INdAM–GNCS) and of the SIMAI’s Activity Group on “Multiscale Modelling of Interacting Agents”.

6 Additional tests for the linear elliptic problem↩︎

We propose additional results on the linear elliptic problem discussed in Section 4.3. In particular, we consider a different initialization of the ensemble based on low-frequency sine modes. Let \[\mathcal{S}_\ell(s)=\sin(\ell s), \qquad \ell=1,\dots,L.\] We define the corresponding discrete basis vectors on the grid \(\{s_j\}_{j=1}^d\) and normalize these vectors in \(\mathbb{R}^d\), and use them to construct the initial particles.

More precisely, let \(\mathbf{S}\in\mathbb{R}^{d\times L}\) be the matrix whose columns are the normalized discrete sine modes, we define the initial ensemble by \[\label{eq:init95sine} x_j^0 = \sigma_0 \sum_{\ell=1}^L \xi_{\ell,j}\,\mathbf{S}_{\ell j}, \qquad j=1,\dots,J,\tag{41}\] where, \(\sigma_0>0\) controls the spread of the initialization and \(\xi_{\ell,j}\sim\mathcal{N}(0,1)\) are independent random coefficients. For the numerical experiments, we choose a moderate number of low-frequency modes, for example \(L=6\), and different values of \(\sigma_0\).

7 Additional results on the Darcy tuning study↩︎

For completeness, we report the parameter configurations and the aggregate tuning results used to select the second-order EKI parameters in Section 4.5.1. The configurations are obtained from the grid \[\gamma\in\{1,2\}, \qquad \beta\in\{0.5,1\}, \qquad k\in\{0.01,0.05,0.1\},\] with \(\alpha=0.2\), \(\varepsilon=10^{-3}\), and \(p=1.5\).

References↩︎