data assimilation, ensemble filtering, data-driven modeling, stochastic optimal control, adjoint matching

1 Introduction↩︎

Multiscale modeling in scientific computing concerns the relation between an underlying dynamical description and the quantities accessible to computation, experiment, or measurement [1], [2]. From a forward perspective, one derives, approximates, or learns reduced mesoscopic or macroscopic dynamics for relevant observables [3], [4]. A complementary inverse perspective seeks to recover the evolving underlying state from partial, noisy, or indirect observations rather than to construct a closed model for those observables. This is the setting of data assimilation (DA), which combines dynamical models with observational data to infer states, trajectories, and other latent quantities [5], [6]. Originating in numerical weather prediction [7], DA is now widely used in oceanography, geophysics, ecology, and other partially observed dynamical systems [8]–[11].

In this work, we focus on the sequential data-assimilation setting. Sequential DA tracks the evolving state by alternating between a forecast step, which propagates the current state uncertainty through the dynamical model, and an analysis step, which incorporates newly available observations [5], [6]. Among numerical realizations of this recursion, ensemble-based methods are the scalable approach most widely used in large-scale systems [12]. They represent the forecast uncertainty by a finite collection of sample states and update this ensemble whenever new data arrive. Once a forecast ensemble has been generated, the central numerical task considered here is the analysis transformation mapping from the forecast distribution \(\pi^{\mathrm{f}}\) to the Bayesian analysis distribution \(\pi^{\mathrm{a}}\propto\ell(\cdot;y)\pi^{\mathrm{f}}\) by taking into account the likelihood model \(\ell\) and the observation data \(y\).

Existing ensemble filters differ mainly in the computational interface through which they use the observation model. For additive-Gaussian observations, the EnKF constructs its analysis update from empirical first and second moments [12], [13]. The bootstrap particle filter evaluates the likelihood at forecast particles and then resamples according to the resulting importance weights [14], [15]. Among recent flow methods, EnSF evaluates the observation-likelihood score along a diffusion path, whereas EnFF uses likelihood-reweighted Monte-Carlo guidance or a localized likelihood-gradient approximation [16], [17]. These interfaces are effective when the observation model provides the required quantities, but many scientific sensing mechanisms do not. The map from state to observation can be non-smooth, many-to-one, implicit, or available only through a simulator. For example, a quantized sensor reports only the interval in which a value falls [18], an intensity measurement records magnitude but discards phase [19], and a radiative-transfer code may provide forward samples without a tractable likelihood [20]. These observations can still define informative analysis updates, but they may not supply a reliable additive residual, evaluable likelihood weights, or analytic likelihood guidance along a transport path. This motivates separating the distribution targeted by the analysis from the numerical mechanism used to realize the update.

We encode the observation information by a scalar energy \(J(x;y)\) and define the target analysis law by \(\pi^{\mathrm{a}}(x\mid y)\propto\pi^{\mathrm{f}}(x)\exp[-J(x;y)]\). When \(J(x;y)=-\log p(y\mid x)\), this is the standard Bayesian posterior. More generally, \(J\) may be a compatibility energy that is not a normalized observation density, or a surrogate learned from simulator samples. We call this formulation implicit data assimilation. Because the target is defined as a tilt of the forecast law rather than as a Gaussian moment update, it can represent multiple state configurations compatible with a many-to-one observation.

To approximate this target with a finite ensemble, the ensemble controlled-flow filter introduces an auxiliary transport time and steers a base diffusion whose terminal law is the forecast distribution. Under the memoryless construction and exact optimal control, the controlled terminal law equals the energy-based analysis law. In the practical method, the optimal control is approximated by adjoint matching [21]: terminal gradients \(\nabla J(Z_1;y)\) are propagated backward through the base dynamics to form regression targets for the control. EnCF therefore requires a differentiable analytic energy, but it does not prescribe an observation-likelihood gradient at every transport-time step. 1 summarizes this design. The observation reaches the analysis only through the energy interface, and one assimilation cycle realizes the tilt as a controlled transport of the forecast ensemble.

When the observation is represented by a differentiable analytic energy, EnCF uses that energy directly. When only simulator samples are available, when a Gaussian residual model gives the wrong conditional law, or when direct gradients are unavailable or uninformative, its likelihood-free variant EnCF-LF learns a smooth surrogate energy function and uses its gradient in the same adjoint-matching solver. We analyze EnCF as a numerical scheme for the Bayes filtering map. It is exact with an exact base flow and optimal control, its one-step error splits into separately refinable approximation terms, and these local errors do not accumulate when the exact filter is stable. The experiments indicate that the controlled-flow formulation is most useful in regimes with non-Gaussian, multimodal, many-to-one, or implicit observation models. The likelihood-free variant addresses the weaker-access case in which only forward simulator samples are available, by learning the energy used by the same controlled-flow solver.

Our contributions are summarized as follows.

1. Implicit data assimilation. We formulate the analysis update through a scalar observation energy \(J\), defining the target law as the exponential tilt \(\pi^{\mathrm{a}}\propto\pi^{\mathrm{f}}\mathrm{e}^{-J}\) of the forecast distribution. This formulation separates the specification of the observation model from the numerical analysis mechanism and accommodates explicit likelihoods, implicit constraints, joint models, and learned surrogate energies without requiring an explicit forward observation map.

2. Ensemble controlled-flow filter. We introduce the Ensemble Controlled-flow Filter, which realizes the energy-defined analysis as a stochastic optimal-control problem. The observation-dependent feedback control is approximated by adjoint matching from terminal energy gradients. Its likelihood-free variant, EnCF-LF, learns a conditional surrogate energy from simulator samples and applies the same controlled-flow solver.

3. Error analysis. We establish exactness in the ideal limit and derive a one-step total-variation error decomposition that separates the effects of the empirical base flow, control approximation, surrogate energy, and time discretization. Under exponential stability of the exact filtering recursion, we further show that uniformly bounded local errors yield a uniformly bounded multi-step filtering error.

4. Regime-based evaluation. We compare EnCF and EnCF-LF with classical and flow-based ensemble filters across additive-Gaussian, non-Gaussian, many-to-one, multimodal, and implicit observation models. The results empirically validate the intended regime of the method: tuned Kalman filters remain preferable for smooth additive-Gaussian observations, whereas controlled transport becomes advantageous when the analysis geometry is non-Gaussian, multimodal, or implicitly defined.

2 Background and related work↩︎

Consider the discrete-time state-space model \[\label{eq:ssm} x_k=\mathcal{M}_k(x_{k-1};\eta_k)\in\mathbb{R}^n,\qquad\eta_k\sim\mathcal{D}_k,\tag{1}\] where \(\mathcal{M}_k\) is a (possibly stochastic) prediction operator and \(\eta_k\) denotes the model noise distributed according to \(\mathcal{D}_k\). At each time step \(k\), we observe the system \(x_k\) through an observation \(y_k\in\mathbb{R}^m\) satisfying \(y_k\sim p(\cdot\mid x_k)\). We assume \(p(y_k\mid x_k,y_{<k})=p(y_k\mid x_k)\), i.e., the current observation is conditionally independent of previous observations given the current state. We denote by \[\pi^{\mathrm{f}}_k=\pi^{\mathrm{f}}_k(x_k)\mathrel{\vcenter{:}}= p(x_k\mid y_{<k})\quad\text{and}\quad\pi^{\mathrm{a}}_k=\pi^{\mathrm{a}}_k(x_k)\mathrel{\vcenter{:}}= p(x_k\mid y_{\le k})\] the forecast and analysis distributions of the system state \(x_k\), respectively, where \(y_{<k}\) and \(y_{\le k}\) collect the observations before and up to step \(k\).

Sequential data assimilation seeks to recursively track the analysis distribution \(\pi^{\mathrm{a}}_k\) by alternating between a prediction step and an analysis step. Specifically, for each step \(k\), the Chapman–Kolmogorov equation gives \[p(x_k\mid y_{<k})=\int q_k^\mathcal{M}(x_k\mid x_{k-1})p(x_{k-1}\mid y_{<k})\mathop{}\!\mathrm{d}x_{k-1},\] where \(q_k^\mathcal{M}\) denotes the transition kernel induced by \(\mathcal{M}_k\). Applying Bayes’ rule then yields \[p(x_k\mid y_{\le k})=\frac{p(y_k\mid x_k,y_{<k})p(x_k\mid y_{<k})}{p(y_k\mid y_{<k})}=\frac{p(y_k\mid x_k)p(x_k\mid y_{<k})}{p(y_k\mid y_{<k})}.\] The second equality follows from the conditional-independence assumption above. Consequently, sequential data assimilation can be summarized by the two-step recursion \[\pi^{\mathrm{f}}_k=(\mathcal{M}_k)_\sharp\pi^{\mathrm{a}}_{k-1}\quad(\text{prediction}),\qquad \pi^{\mathrm{a}}_k\propto \ell(\cdot;y_k)\pi^{\mathrm{f}}_k\quad(\text{analysis}), \label{eq:recursion}\tag{2}\] where \(\ell(x_k;y_k)\) denotes the likelihood density associated with \(p(y_k\mid x_k)\).

Classical ensemble-based filters realize this recursion by evolving an ensemble of \(N\) particles \(\{x_k^{(j)}\}_{j=1}^N\) over time. In the prediction step, each ensemble member is propagated independently according to \[x_k^{\mathrm{f},(j)}=\mathcal{M}_k(x_{k-1}^{\mathrm{a},(j)};\eta_k^{(j)}),\quad \eta_k^{(j)}\sim\mathcal{D}_k.\] If the previous analysis ensemble approximates \(\pi^{\mathrm{a}}_{k-1}\), independently propagating its members produces an ensemble approximation of \(\pi^{\mathrm{f}}_k\). Different ensemble filters mainly differ in how they transform the forecast ensemble members so as to approximate the analysis distribution \(\pi^{\mathrm{a}}_k\propto \ell(\cdot;y_k)\pi^{\mathrm{f}}_k\). For an additive-Gaussian likelihood \[\ell(x;y_k)=\mathcal{N}(y_k;\mathcal{H}(x),R_k),\] in which \(R_k\) is the observation-error covariance, the EnKF approximates the forecast distribution by a Gaussian whose moments match those of the forecast ensemble, and then shifts each member by the empirical Kalman gain, formulated as \[x_k^{\mathrm{a},(j)}=x_k^{\mathrm{f},(j)}+P_{xy}^\mathrm{f}\bigl(P_{yy}^\mathrm{f}+R_k\bigr)^{-1}\bigl[y_k^{(j)}-\mathcal{H}(x_k^{\mathrm{f},(j)})\bigr],\qquad y_k^{(j)}\sim\mathcal{N}(y_k,R_k). \label{eq:enkf}\tag{3}\] The matrices \(P^\mathrm{f}_{xy},P^\mathrm{f}_{yy}\) denote the cross-covariance and predicted observation covariance estimated from the ensembles, respectively. The EnKF recovers the exact analysis distribution as \(N\to\infty\) when the forecast is Gaussian and \(\mathcal{H}\) is linear; finite \(N\) introduces covariance-estimation and perturbed-observation sampling errors. Outside the linear–Gaussian setting, an update determined by empirical first and second moments need not reproduce a non-Gaussian or multimodal analysis distribution.

The bootstrap particle filter instead preserves the empirical forecast measure and reweights it by the likelihood as \[w_k^{(j)}=\frac{\ell(x_k^{\mathrm{f},(j)};y_k)}{\sum_{j'=1}^N \ell(x_k^{\mathrm{f},(j')};y_k)},\qquad \tilde{\pi}^{\mathrm{a}}_k\mathrel{\vcenter{:}}=\sum_{j=1}^N w_k^{(j)}\delta_{x_k^{\mathrm{f},(j)}}. \label{eq:pf}\tag{4}\] The analysis ensemble \(\{x_k^{\mathrm{a},(j)}\}_{j=1}^N\) is then obtained by resampling from \(\tilde{\pi}^{\mathrm{a}}_k\). For a fixed assimilation horizon, the bootstrap particle approximation converges to the filtering distribution as \(N\to\infty\) under the usual sequential-Monte-Carlo integrability assumptions, but it faces the curse of dimensionality and suffers from severe degeneracy in high-dimensional settings [15], [22].

Flow-based filters construct stochastic or deterministic transports whose terminal laws approximate the analysis distribution. EnSF constructs a score-based diffusion transport [16], whereas EnFF constructs a flow-matching transport [17]. Related work extends score- and flow-based assimilation to partially observed systems, learned latent dynamics, and stochastic-interpolant models [23]–[25]. EnSF evaluates \(\nabla_z\log p(y_k\mid z)\) along the diffusion path, whereas EnFF uses either likelihood-reweighted Monte-Carlo guidance or a localized likelihood-gradient approximation. Data-driven methods address unknown or unresolved dynamics by learning or replacing forecast-model components from training data [26]–[29]. Other methods perform assimilation in learned reduced-order spaces [30]–[33].

The comparison above exposes three limitations of existing analysis interfaces. The EnKF reduces the update to empirical first and second moments and can therefore miss non-Gaussian or multimodal analysis distributions. The bootstrap particle filter is asymptotically consistent, but its observation-likelihood weights can degenerate in high dimensions. EnSF and EnFF avoid explicit likelihood-weighted resampling, yet they require an evaluable likelihood or likelihood gradient to prescribe observation-dependent guidance along the transport path. The data-driven and reduced-order methods above primarily address unknown or expensive dynamics and representation cost; here we instead assume access to forecast samples and isolate the analysis update. We seek an update that does not impose a Gaussian approximation, does not weight and resample forecast particles by the observation likelihood, and separates the specification of the observation from the transport mechanism. The next section defines the target energy-based analysis density \(\pi^{\mathrm{a}}_k(x)\propto\pi^{\mathrm{f}}_k(x)\exp[-J_k(x;y_k)]\) and constructs a controlled-flow approximation that learns the control by adjoint matching from terminal energy gradients rather than prescribing a likelihood-based guidance rule at every transport-time step.

3 The ensemble controlled-flow filter↩︎

In this section, we formulate the analysis step as a controlled transport from the forecast law to an energy-tilted analysis law. We begin by introducing implicit data assimilation, in which the observation enters through a scalar energy rather than an additive-Gaussian residual or a prescribed analytic representation of the likelihood. We then cast the resulting analysis update as a stochastic optimal-control problem and use this formulation to place several existing flow-based filters within a common controlled-flow perspective [34]. Building on this viewpoint, we develop the ensemble controlled-flow filter, which approximates the optimal feedback control through adjoint matching. Finally, when an analytic observation energy is unavailable but forward simulation of the observation mechanism is possible, we introduce the likelihood-free variant EnCF-LF, which learns a surrogate energy from simulator samples and uses the same controlled-flow solver.

3.1 Implicit data assimilation via energy tilts↩︎

Recall that the classical analysis step in 2 reweights the forecast law by the likelihood \(\ell(\cdot;y_k)\). Particular numerical filters often require more restrictive structures or computational interfaces, such as an additive-Gaussian residual, a moment-based closure, or the likelihood score along a transport path. We instead use a scalar energy as the computational interface between the observation model and the analysis solver.

The energy \(J_k(x;y_k)\) measures the compatibility between a candidate state \(x\) and the observation \(y_k\), and is defined up to the addition of a term independent of \(x\). It accommodates several observation interfaces.

When an evaluable likelihood is available, taking \[J_k(x;y_k)=-\log \ell(x;y_k) \label{eq:likelihood-energy}\tag{7}\] recovers the standard Bayesian analysis. This includes non-Gaussian and many-to-one observation models and is not restricted to the additive-Gaussian setting.

The observation relation may instead be specified implicitly through \(\mathcal{R}_k(x,y_k)\approx0\), without requiring an explicit forward representation \(y_k=\mathcal{H}(x_k)+\varepsilon_k\). A soft enforcement of this relation is obtained from \[J_k(x;y_k) = \frac{\left\| \mathcal{R}_k(x,y_k) \right\|^2}{2\sigma_{\mathrm{tol}}^2}, \label{eq:constraint-energy}\tag{8}\] where \(\sigma_{\mathrm{tol}}>0\) determines the tolerance assigned to violations of the relation.

More generally, any strictly positive compatibility potential \(\Psi_k(x;y_k)\) can be represented by \[J_k(x;y_k)=-\log\Psi_k(x;y_k),\] provided that 5 holds. In this form, 6 defines a Gibbs-type measure with the forecast law as its reference measure, \(J_k\) as the energy, and \(\mathcal{Z}_k\) as the partition function [35]. When \(\Psi_k\) is the likelihood, this measure coincides with the Bayesian posterior. A general compatibility potential instead defines an energy-based analysis that need not correspond to a normalized observation law.

When the observation mechanism is accessible only through forward simulation, a differentiable surrogate energy \(E_\phi(x;y_k)\) may be learned from simulated state–observation pairs and substituted for \(J_k\). This construction is developed in 3.4.

In all these cases, the formulation reweights the full forecast law without selecting a particular inverse branch. It can therefore represent multiple well-separated state configurations that are compatible with the same observation.

3.2 The optimal controlled flow↩︎

To realize the energy-tilted analysis law in 6 , we introduce an auxiliary transport time \(\tau\in[0,1]\) within each assimilation cycle \(k\). Let \(\mathbb{P}\) denote the path law of the base diffusion \[\mathop{}\!\mathrm{d}Z_\tau = f(Z_\tau,\tau)\mathop{}\!\mathrm{d}\tau+\sigma(\tau)\mathop{}\!\mathrm{d}W_\tau, \qquad Z_0\sim\mu_0, \quad \tau\in[0,1], \label{eq:osde}\tag{9}\] whose terminal law is the forecast law, \(\mathrm{Law}_{\mathbb{P}}(Z_1)=\pi^{\mathrm{f}}_k\). The base drift generally depends on the forecast law at cycle \(k\), although this dependence is suppressed in the notation. We then steer the base flow toward the analysis law \(\pi^{\mathrm{a}}_k\) through an observation-dependent drift \(u\), \[\mathop{}\!\mathrm{d}Z_\tau = \left[ f(Z_\tau,\tau)+u(Z_\tau,\tau;y_k) \right]\mathop{}\!\mathrm{d}\tau +\sigma(\tau)\mathop{}\!\mathrm{d}W_\tau, \qquad Z_0\sim\mu_0, \quad \tau\in[0,1]. \label{eq:csde}\tag{10}\] We denote the path law of 10 by \(\mathbb{Q}^u\). This representation also includes deterministic flows when \(\sigma=0\). For the stochastic optimal-control construction below, we assume that \(\sigma(\tau)>0\).

The control is chosen by minimizing the free-energy objective \[\mathcal{J}_k(u) \mathrel{\vcenter{:}}= \mathbb{E}_{\mathbb{Q}^u} \left[ \frac{1}{2} \int_0^1 \left\| \sigma(\tau)^{-1}u(Z_\tau,\tau;y_k) \right\|^2 \mathop{}\!\mathrm{d}\tau + J_k(Z_1;y_k) \right]. \label{eq:soc}\tag{11}\] The running term penalizes deviation from the base path law, whereas the terminal term favors states compatible with the observation. Under the corresponding Girsanov change of measure, the running cost equals \(D_{\mathrm{KL}}\left(\mathbb{Q}^u\,\middle\|\,\mathbb{P}\right)\) for path laws sharing the initial distribution \(\mu_0\).

These are standard conditions for the Doob \(h\)-transform associated with a terminal energy [36], [37]. For the VP construction below, the Gaussian transition kernel provides interior regularization under suitable integrability conditions.

The memoryless condition is essential because otherwise the normalization of the terminal tilt generally depends on \(Z_0\), and the controlled terminal law need not coincide with the global analysis law in 6 .

Several existing flow-based ensemble filters admit the decomposition in 10 , including the score-based EnSF [16], the flow-matching EnFF [17], and the Schrödinger-bridge EnSBF [40]. They differ in their initial law, base dynamics, and observation-dependent drift. This common representation does not imply that all these methods solve the stochastic optimal-control problem 11 . Rather, it isolates how each method constructs its forecast transport and incorporates the observation. The resulting decompositions are summarized in 1.

3.3 The adjoint-matching control↩︎

We now construct EnCF by specifying a memoryless base flow and learning its observation-dependent control. Let \(\alpha_\tau\) and \(\beta_\tau\) be interpolation coefficients satisfying \((\alpha_0,\beta_0)=(0,1)\) and \((\alpha_1,\beta_1)=(1,0)\) so that the Gaussian kernel \(q_{\tau\mid0}(z\mid x)=\mathcal{N}(z;\alpha_\tau x,\beta_\tau^2 I_n)\) interpolates from the reference law \(\mu_0=\mathcal{N}(0,I_n)\) to a point mass at \(x\). The associated reverse variance-preserving (VP) diffusion [41], [42] has coefficients \[b(\tau) = \frac{\dot{\alpha}_\tau}{\alpha_\tau}, \qquad \sigma^2(\tau) = 2b(\tau)\beta_\tau^2 - \frac{\mathop{}\!\mathrm{d}\beta_\tau^2}{\mathop{}\!\mathrm{d}\tau} = \frac{ 2\beta_\tau \bigl( \dot{\alpha}_\tau\beta_\tau - \dot{\beta}_\tau\alpha_\tau \bigr) }{ \alpha_\tau }. \label{eq:memoryless}\tag{13}\] This choice makes the endpoints of the base process independent, \(Z_0\perp Z_1\), and hence satisfies the memoryless condition required in 1.

At assimilation cycle \(k\), let \(\{\hat{Z}_k^{(j)}\}_{j=1}^N\) denote the forecast ensemble approximating \(\pi^{\mathrm{f}}_k\). Then we may obtain an empirical base drift by Monte-Carlo approximation as \[\hat{f}_k(z,\tau) = b(\tau)z + \sigma^2(\tau) \sum_{j=1}^N w_j(z,\tau) \frac{ \alpha_\tau\hat{Z}_k^{(j)}-z }{ \beta_\tau^2 },\quad w_j(z,\tau)=\frac{q_{\tau\mid0}(z\mid\hat{Z}_k^{(j)})}{\sum_{j'=1}^Nq_{\tau\mid0}(z\mid\hat{Z}_k^{(j')})}. \label{eq:EnCF-f}\tag{14}\] The corresponding population drift transports \(\mu_0\) to \(\pi^{\mathrm{f}}_k\), while 14 replaces the forecast law by its ensemble approximation [43].

With the base flow fixed, it remains to approximate the optimal control in ?? . We parameterize the observation-dependent drift as \[u_\theta(z,\tau;y_k) = \sigma^2(\tau)a_\theta(z,\tau;y_k), \label{eq:control}\tag{15}\] where \(a_\theta\) approximates the logarithmic gradient in ?? . In the experiments, \(a_\theta\) is represented by a convolutional neural network, although the construction is independent of this choice.

Direct minimization of \(\mathcal{J}_k(u_\theta)\) in 11 requires differentiating through the sampled controlled trajectories. Adjoint matching [21], [44] instead constructs a least-squares target by solving a lean adjoint equation backward along each sampled path \(Z\sim\mathbb{Q}^{u_\theta}\), with the path treated as stop-gradient. Because the optimal terminal control in ?? is known explicitly at \(\tau=1\) as \(u^\star(z,1;y_k)=-\sigma^2(1)\nabla_z J_k(z;y_k)\), we start from this terminal information and define the lean adjoint by \[\widetilde{a}(1;Z) = \nabla_z J_k(Z_1;y_k), \qquad \frac{\mathop{}\!\mathrm{d}\widetilde{a}(\tau;Z)}{\mathop{}\!\mathrm{d}\tau} = - \nabla_z f(Z_\tau,\tau)^\mathsf{T} \widetilde{a}(\tau;Z). \label{eq:adjoint}\tag{16}\] The control is then trained using the adjoint-matching objective \[\mathcal{L}_{\mathrm{AM}}(\theta;Z) = \frac{1}{2} \int_0^1 \sigma^2(\tau) \left\| a_\theta(Z_\tau,\tau;y_k) + \widetilde{a}(\tau;Z) \right\|^2 \mathop{}\!\mathrm{d}\tau. \label{eq:am}\tag{17}\] Thus \(a_\theta\) is regressed onto \(-\widetilde{a}\), or equivalently, the controlled drift \(u_\theta\) is regressed onto \(-\sigma^2\widetilde{a}\). The observation enters the backward equation only through the terminal seed \(\nabla_zJ_k(Z_1;y_k)\). Moreover, the stop-gradient construction avoids retaining an automatic-differentiation graph through the complete controlled rollout.

3.4 The likelihood-free surrogate energy↩︎

When the observation law is unavailable but observation samples can be generated, we replace \(J_k\) by a surrogate energy \(E_\phi\). We use two parameterizations according to the available observation structure.

When a meaningful observation mapping exists, Gaussian-LF learns a pointwise conditional model with mean \(m_\phi(x)\) and scale \(s_\phi\). Its energy is \[E_\phi(x;y) = \frac{1}{2} \left\| \frac{m_\phi(x)-y}{s_\phi} \right\|^2,\] up to terms independent of \(x\). The parameters are fitted online by Gaussian maximum likelihood on simulator pairs \((x^{(j)},y^{(j)})\). This model is appropriate for smooth additive observations. For more general observations such as multimodal or implicit conditional laws, the energy is directly parameterized as a scalar network. The network is trained by noise-contrastive estimation, using matched simulator pairs as positives and mismatched observations together with perturbed states as negatives. At each assimilation cycle, \(E_\phi\) is first trained with the control fixed and is then held fixed during adjoint matching. In both variants, the controlled-flow solver uses the surrogate only through the terminal gradient \(\nabla_x E_\phi(Z_1;y_k)\).

4 Theory↩︎

Building on the controlled-flow formulation developed in 3, this section establishes structural a priori guarantees for the resulting filtering scheme along three complementary directions. We first study consistency in the ideal limit, showing that an exact base flow together with the optimal control recovers the exact filtering distribution. We then characterize local accuracy by decomposing the one-step filtering error into refinable contributions from the empirical base flow, control approximation, surrogate energy, and time discretization. Finally, using the perturbation framework for Markov filters [45], we establish stability of the approximate recursion by showing that, whenever the exact filter is stable, these locally injected errors remain uniformly controlled rather than accumulating over assimilation cycles.

4.1 Consistency in the ideal limit↩︎

We first consider the idealized setting in which the forecast propagation, the base flow, and the control approximation are exact. In this setting, the controlled-flow construction recovers the target analysis law at every assimilation cycle.

In practice EnCF departs from this ideal by four separately measurable amounts, which the next subsection isolates in 1 and bounds in 6, a base-drift Monte-Carlo error, the suboptimality of the learned control, the surrogate energy gap, and the time-discretization of the controlled flow.

4.2 Local accuracy and one-step error decomposition↩︎

We isolate four sources of one-step approximation error under two structural assumptions. 1 first decomposes the error into these contributions. We then bound the terms that admit analytic control and assemble the resulting a priori estimate in 6.

Formally, starting from the previous analysis \(\nu\mathrel{\vcenter{:}}=\pi^{\mathrm{a}}_{k-1}\), the exact one-step filtering operator \(\Phi_k\mathrel{\vcenter{:}}= \mathcal{B}_{k}\circ(\mathcal{M}_k)_\sharp\) consists of a predictor \((\mathcal{M}_k)_\sharp:\pi^{\mathrm{a}}_{k-1}\mapsto\pi^{\mathrm{f}}_k\) and an exact analysis operator \(\mathcal{B}_{k}(\nu)(\mathop{}\!\mathrm{d}x)\propto \mathrm{e}^{-J_k(x;y_k)}\nu(\mathop{}\!\mathrm{d}x)\). In this subsection, we study how the practical update operator \(\hat{\Phi}_k\) differs from the exact \(\Phi_k\).

Fix an assimilation cycle \(k\). Let \(\widehat{\mathbb{Q}}\) denote the path law of the practical continuous-time flow \[\mathrm dZ_\tau = \bigl[ \hat{f}(Z_\tau,\tau) + u_\theta(Z_\tau,\tau;y_k) \bigr]\mathrm d\tau + \sigma(\tau)\mathop{}\!\mathrm{d}W_\tau, \qquad Z_0\sim\mu_0,\] and let \(\mathbb{Q}^\star\) denote the path law of the ideal flow

\[\mathrm dZ_\tau = \bigl[ f(Z_\tau,\tau) + u^\star(Z_\tau,\tau;y_k) \bigr]\mathrm d\tau + \sigma(\tau)\mathop{}\!\mathrm{d}W_\tau, \qquad Z_0\sim\mu_0.\] We write \(\widehat{\mathbb{Q}}_1\) and \(\mathbb{Q}^\star_1\) for their terminal marginals, and \(\widehat{\mathbb{Q}}^{\mathrm{EM}}_1\) for the terminal law produced by the Euler–Maruyama discretization of the practical flow. Thus, the practical terminal law \(\hat{\mathbb{Q}}^{\mathrm{EM}}_1=\hat{\pi}^{\mathrm{a}}_k\mathrel{\vcenter{:}}=\hat{\Phi}_k(\pi^{\mathrm{a}}_{k-1})\) and the ideal terminal law \(\mathbb{Q}^\star_1=\pi^{\mathrm{a}}_k=\Phi_k(\pi^{\mathrm{a}}_{k-1})\).

We then discuss the four terms in turn. The decomposition identifies the numerical levers of the method, the ensemble representation and base drift for the base-drift term, the adjoint matching optimization error for the control term, the surrogate gradient mismatch for the energy term, and the number of integration steps \(N_\tau\) for the discretization term. These terms have different theoretical status. The base-flow and control errors are retained as approximation residuals, whereas the energy and discretization terms admit explicit analytic bounds under additional regularity assumptions.

For the base-drift term, 8 proves in the appendix the self-normalized importance-sampling rate \(\mathcal{O}(N^{-1/2})\), with a constant governed by the \(\chi^2\) dispersion of the base weights. The control term measures the optimization errors between the learned control and the nonparametric fixed point induced by the on-policy adjoint-matching regression as detailed in 2. We retain this term as a residual because a general convergence rate would additionally require approximation estimates for the control class and quantitative guarantees for the optimization procedure, which is beyond our scope in this paper.

Two consequences make this the pivot of the implicit-energy and likelihood-free design. First, although \(u_J^\star\) is nonlocal in \(J\), replacing \(J\) by a learned surrogate perturbs the filtering map only through \(\left\| \nabla E_\phi-\nabla J \right\|_\infty\) under the stated assumptions, so energy learning by adjoint matching is stable in the control metric. Besides, the control is invariant to shifts \(J\mapsto J+c\) in \(x\). EnCF-LF therefore need not estimate the tilt normalizer over the state, making the maximum-likelihood training reasonable.

Combining the estimate of each component in the decomposition 18 , we have the following aggregated result.

6 separates the local filtering error according to the four approximations made by the practical method. The first term depends on the accuracy with which the empirical base dynamics represent the forecast transport. The second reflects the expressiveness and optimization accuracy of the learned control. The third is governed by the \(C^1\)-accuracy of the surrogate energy, while the fourth decreases as the transport-time discretization is refined. The bound is therefore structural and a priori rather than an online error certificate, since the exact base drift \(f\) and optimal control \(u^\star\) are generally unavailable. The constant \(L_\mathrm{E}\) also indicates that surrogate-energy errors may be amplified when the tilted bridge measures become poorly conditioned, as reflected by large Poincaré or score-variance constants. This is consistent with the intuition that sharply concentrated observation models tend to be more challenging.

4.3 Multi-step stability and error propagation↩︎

The one-step estimate of 6 quantifies the local error introduced during a single assimilation cycle. To propagate this estimate over multiple cycles, we compare the practical recursion with the exact filtering recursion started from the same intermediate distributions. The key requirement is stability of the exact filtering operators with respect to perturbations of their input laws.

To study the accumulated filtering error, we write \(\Phi_{k\to l}\mathrel{\vcenter{:}}=\Phi_{l}\circ\cdots\circ\Phi_{k+1}\) for \(k\le l\) and set \(\Phi_{k\to k}\) as the identity. Then we define \[\label{eq:dk-ek-def} d_k\mathrel{\vcenter{:}}=\mathrm{TV}\left(\hat{\pi}^{\mathrm{a}}_k,\pi^{\mathrm{a}}_k\right),\quad e_k\mathrel{\vcenter{:}}=\mathrm{TV}\left(\hat{\pi}^{\mathrm{a}}_k,\Phi_k(\hat{\pi}^{\mathrm{a}}_{k-1})\right)=\mathrm{TV}\left(\hat{\Phi}_k(\hat{\pi}^{\mathrm{a}}_{k-1}),\Phi_k(\hat{\pi}^{\mathrm{a}}_{k-1})\right)\tag{24}\] as the accumulated filtering error and the one-step filtering error, respectively. Recall that we have actually provided an upper bound of \(e_k\) in 6.

Estimates of the form 25 are classical under suitable uniform mixing conditions on the forecast kernels and upper and lower bounds on the observation potentials [45], [48]. Note that 4 is a structural property of the exact filtering recursion independent of the choice of numerical implementation.

7 does not assert that the practical filtering map is itself contractive. Rather, it shows that local perturbations of a stable exact filtering recursion are geometrically attenuated as they propagate forward in time. Thus, uniformly bounded local errors produce a uniformly bounded global error, with amplification factor \(C_{\mathrm{stab}}/(1-\varrho)\).

The one-step estimate 6 remains an a priori structural bound because its base-flow and control terms involve the unavailable exact drift and optimal control. 7 complements this local estimate by describing how such errors propagate across assimilation cycles. When \(\varrho\) is close to one, the exact filter forgets perturbations slowly and the resulting multi-step bound becomes correspondingly less informative.

5 Experiments↩︎

This section evaluates the proposed EnCF and its likelihood-free variant EnCF-LF against classical ensemble filters and recently developed flow-based filters. The experiments are organized according to the structure of the observation model. We first consider a smooth additive-Gaussian setting as a sanity check (5.3), where classical Kalman-type filters are expected to perform well. We then examine a non-Gaussian additive noise model (5.4), which progressively moves the posterior away from a Gaussian approximation. Finally, we focus on the observation regimes targeted by the proposed formulation: a many-to-one sign-blind sensor that induces a genuinely multimodal posterior (5.5.1), and an implicit, non-additive observation relation for which the likelihood-free variant uses simulator access only (5.5.2). The computational cost of learning the controlled flow is reported separately in 5.6.

5.1 Filter configurations↩︎

The proposed filters, EnCF and EnCF-LF, are implemented following 2. The control \(a_\theta\) is represented by a convolutional neural network. The state \(z\) and the observation \(y\) are passed through separate convolutional encoders, and the resulting feature maps are concatenated and decoded by a merge convolution to produce the control. The network is initialized at zero. In EnCF, the observation energy enters the adjoint-matching regression only through the terminal condition \(\nabla J_k(Z_1;y_k)\), as described in 2. Full implementation details and hyperparameters are provided in 8.

At each assimilation step \(k\), the control network of our EnCF is trained from scratch. We therefore use more optimization in the early cycles and gradually reduce the training budget. In all experiments, the number of training epochs \(E_k\) follows a cosine schedule, inspired by [49], decreasing from \(100\) to \(1\) over the first \(20\) assimilation cycles. After this warm-up phase, each analysis step uses a single training epoch. The number of regression steps is set as \({N_\mathrm{reg}}=16\).

We compare against two classes of ensemble filters. The classical baselines are the EnKF [12], the LETKF [13], and the bootstrap particle filter [50]. For EnKF and LETKF, covariance inflation and Gaspari–Cohn localization are tuned separately for each benchmark; the tuning protocol is reported in 8. The particle filter uses the exact likelihood whenever it is available, but is included mainly as a reference method because it is known to suffer from weight degeneracy in high-dimensional problems [22]. The flow-based baselines are EnSF [16], EnFF [17], and EnSBF [40]. For EnFF, we report both the optimal-transport variant (EnFF-OT) and the filtering-to-predictive variant (EnFF-F2P). All flow-based baselines use the same number of transport-time steps \(N_\tau=100\), matching the EnCF implementation.

5.2 Metrics↩︎

For each benchmark, we run the filter for \(K\) assimilation cycles and report time-averaged quantities over a converged window \(\mathcal{C}\), defined as the last \(10\%\) of the rollout. We use three complementary metrics as follows.

5.2.0.1 Root mean square error (RMSE)

We evaluate the accuracy of the ensemble mean \(\bar x_k=N^{-1}\sum_{j=1}^N x_k^{(j)}\) by computing the root mean square error against the reference state \(x_k^\star\in\mathbb{R}^n\), defined by \[\mathrm{RMSE}_k=\sqrt{\frac{1}{n}\sum_{i=1}^n\bigl(\bar x_{k,i}-x_{k,i}^\star\bigr)^2}, \qquad \mathrm{RMSE}=\frac{1}{|\mathcal{C}|}\sum_{k\in\mathcal{C}}\mathrm{RMSE}_k, \label{eq:rmse}\tag{27}\] where \(i\) indexes coordinates. Lower RMSE indicates a more accurate ensemble mean.

5.2.0.2 Continuous ranked probability score (CRPS)

Since the analysis output is an ensemble distribution rather than only a point estimate, we also report the multivariate energy-score form of the continuous ranked probability score \[\mathrm{CRPS}_k=\frac{1}{N}\sum_{j=1}^N\left\| x_k^{(j)}-x_k^\star \right\|_2 -\frac{1}{2N^2}\sum_{j=1}^N\sum_{j'=1}^N\left\| x_k^{(j)}-x_k^{(j')} \right\|_2. \label{eq:crps}\tag{28}\] The reported CRPS is averaged over \(k\in\mathcal{C}\). The first term penalizes distance from the reference state, whereas the second term rewards ensemble spread. Thus CRPS evaluates the ensemble as a distribution, with lower values indicating a sharper and better calibrated ensemble.

5.2.0.3 Spread-to-RMSE ratio (SRR)

We use \[\mathrm{Spread}_k=\sqrt{\frac{1}{nN}\sum_{j=1}^N\sum_{i=1}^n\bigl(x_{k,i}^{(j)}-\bar x_{k,i}\bigr)^2}, \qquad\mathrm{SRR}_k=\frac{\mathrm{Spread}_k}{\mathrm{RMSE}_k},\] as a calibration diagnostic. The final SRR is again averaged over \(k\in\mathcal{C}\). Values near one indicate that the ensemble spread is consistent with the realized error. Values below one indicate under-dispersion or ensemble collapse, whereas values above one indicate over-dispersion.

All configurations are run with five independent random seeds. Tables report the mean and standard deviation over seeds, and figures show the same variability using error bars or shaded bands.

5.3 A conventional assimilation task↩︎

We begin with a standard setting in which the observation model is smooth, dense, and corrupted by additive Gaussian noise. This experiment is intended as a sanity check, and the purpose is to verify that our controlled-flow formulation remains stable and competitive before moving to the non-Gaussian, multimodal, and implicit observation models considered afterward.

We consider the Kuramoto–Sivashinsky (KS) dynamics \[u_t+uu_x+u_{xx}+u_{xxxx}=0\] on a periodic domain discretized into \(1024\) grid points. The state is observed through the componentwise nonlinear sensor \[\label{eq:ks-arctan} y=\mathcal{H}(x)+\varepsilon\in\mathbb{R}^{1024},\qquad\mathcal{H}(x)=\arctan(x),\quad \varepsilon\sim\mathcal{N}(0,0.1^2 I_{1024}).\tag{29}\] We run \(1000\) assimilation cycles with ensemble size \(N=20\).

The results in 2 are consistent with the expected behavior in this conventional regime. The tuned EnKF and LETKF achieve the lowest RMSE and CRPS, confirming that classical covariance-based updates remain preferable when the observation model is smooth, dense, and additive Gaussian. Among the flow-based methods, EnCF gives the best RMSE and CRPS, while EnCF-LF remains stable despite using a learned surrogate observation energy. Several existing flow-based baselines, especially EnFF-F2P and EnSBF, show very small SRR values, indicating strong under-dispersion or ensemble collapse. This sanity check provides a baseline for the more challenging non-Gaussian, multimodal, and implicit observation regimes considered in the following subsections.

5.4 Non-Gaussian additive noise↩︎

We next keep the KS dynamics and the arctangent observation operator in 29 , but replace the Gaussian observation noise by a symmetric two-component mixture. Coordinate-wise, the observation model is \[y=\mathcal{H}(x)+\varepsilon,\quad \mathcal{H}(x)=\arctan(x),\qquad \varepsilon_i\sim\mathcal{N}(a\xi_i\sigma,\sigma^2), \quad \sigma=0.1, \label{eq:bimodal-noise}\tag{30}\] where the signs \(\xi_i\) are independent Rademacher variables. The scalar parameter \(a\ge0\) controls the separation of the two noise modes in units of \(\sigma\). Thus \(a=0\) recovers the Gaussian setting in 5.3, while larger values of \(a\) produce increasingly separated likelihood modes. 3 (a) and (b) show how the likelihood changes as \(a\) increases. The noise density develops two separated modes, and the corresponding observation energy \(J=-\log \ell(\cdot\mid y)\) changes from a single well to a double-well landscape.

Note that for a fixed \(a\), the marginal noise distribution has mean zero and variance \(\sigma^2(1+a^2)\). Therefore, in addition to the nominal EnKF using \(R_k=\sigma^2 I\) (defined in 3 ), we also report a matched-variance EnKF using \(R_k=\sigma^2(1+a^2)I\). The latter is more reasonable since it uses the correct marginal variance of the mixture, although it still replaces the bimodal likelihood by a single Gaussian approximation. 3 (c) reports the converged RMSE for various offsets \(a\). For small \(a\), the likelihood remains close to Gaussian, and the Kalman update is the most effective choice. As \(a\) increases, the two likelihood modes separate, and a single Gaussian residual model becomes less representative of the posterior. The matched-variance EnKF partly accounts for the larger noise scale, but it still cannot represent the bimodal structure of the likelihood. In contrast, EnCF directly uses the observation energy and remains comparatively insensitive to the mode separation. The RMSE curves cross around \(a\approx4\), and beyond this point, EnCF gives lower error than both EnKF variants. This experiment therefore isolates a transition from a regime where Gaussian covariance updates are sufficient to one where the controlled-flow approach becomes advantageous.

5.5 Multimodal posteriors↩︎

We now turn to observation models that induce genuinely multimodal analysis distributions. The common difficulty in this section is that the state-to-observation relation is many-to-one, so a single observation can be compatible with multiple well-separated states. In such cases, a unimodal Gaussian update can collapse distinct posterior modes, while methods relying directly on analytic likelihood gradients may only exploit the local information exposed by the observation map. The two benchmarks below isolate this issue in a controlled setting: the first uses an explicit differentiable sign-blind sensor, and the second replaces the additive observation model by an implicit, non-additive relation.

Both benchmarks in this section share the same double-well dynamics in \(\mathbb{R}^{20}\), \[\mathop{}\!\mathrm{d}x_t=-4 x_t(x_t^2-1)\mathop{}\!\mathrm{d}t+0.5\mathop{}\!\mathrm{d}W_t\] where the nonlinear operations are applied coordinate-wise and \(W_t\) is a standard Brownian motion. The drift is induced by a double-well potential \(V(x)\) with stable wells at \(x=\pm1\), as illustrated in 4. We use \(N=40\) ensemble members and assimilate every \(20\) model steps for \(100\) assimilation cycles.

5.5.1 Sign-blind magnitude observation↩︎

To start with, we consider an explicit many-to-one observation model. The observation consists of a dominant magnitude channel and a weak linear sign channel \[y=\mathcal{H}(x)+\varepsilon\in\mathbb{R}^{40},\qquad\mathcal{H}(x)=(x^2,cx)^{\mathsf{T}},\quad \varepsilon\sim\mathcal{N}(0,0.1^2 I_{40}),\] where all operations are applied coordinate-wise. The first part is sign-blind, and thus we introduce a linear channel with slope \(c=0.02\) that carries the sign at effective noise \(0.1/0.02=5\) times the state scale, making the problem hard but identifiable over time. By this construction, each analysis step remains strongly multimodal.

The RMSE evolution in 4 shows that the Kalman-type filters struggle in this regime even after tuning inflation and localization. The existing flow-based filters improve over the Kalman baselines but remain limited by the artificial controls used to guide the transport. In contrast, EnCF achieves the lowest error among the compared methods. This supports the intended role of the controlled-flow formulation: rather than imposing a Gaussian update or relying on a prescribed likelihood-score control along the path, it learns a transport induced by the terminal observation energy and can better accommodate the two-well posterior structure.

5.5.2 Implicit and non-additive observation↩︎

The sign-blind sensor in 5.5.1 is many-to-one, but it is still given by an explicit observation map with additive Gaussian noise. We now replace this structure by an implicit, non-additive observation relation. The observation consists of a circle-constraint channel and a weak sign channel, given by \[\label{eq:obs-implicit} \mathcal{R}(x,y_{\mathrm c})=x^2+y_{\mathrm c}^2-2\sim\mathcal{N}(0,0.3^2I_{20});\qquad y_s=cx+\eta,\quad\eta\sim\mathcal{N}(0,0.1^2I_{20}).\tag{31}\] All operations are applied coordinate-wise, and we have \(y=(y_{\mathrm c},y_{\mathrm s})\in\mathbb{R}^{40}\) as observation data. In the circle channel, the noise is imposed on the constraint \(x^2+y_{\mathrm c}^2\approx2\), rather than being added to a forward map \(y_{\mathrm c}=H_{\mathrm c}(x)+\varepsilon\). Thus the observation model is not naturally represented as an additive residual in the observation variable. Moreover, the constraint is even in \(x\), so it is again sign-blind; the weak linear channel \(y_{\mathrm s}\) is included only to make the state identifiable over time.

This benchmark separates two levels of access to the observation model. We provide the exact observation relation in 31 to the flow-based methods, while the Kalman-type baselines use the natural residual representation with tuned inflation and localization. In contrast, EnCF-LF is not given the analytic form of 31 . It only receives simulator samples from the conditional observation law at forecast states and learns a surrogate energy online. This setting is intended to mimic an observation mechanism for which forward simulation is available but an explicit likelihood or differentiable observation formula is not.

The results in 3 show that the controlled-flow methods are the most effective in this implicit observation regime. With analytic access to the observation energy, EnCF obtains substantially lower RMSE and CRPS than the classical and existing flow-based baselines. The likelihood-free variant performs comparably, matching the exact-energy EnCF to within a fraction of the seed spread, despite using only simulator samples to learn the surrogate energy. This indicates that the energy-learning interface can recover the relevant observation information without requiring an explicit additive observation map.

5 explains the difference in calibration behavior. The Kalman-type filters and the EnFF variants have small SRR values in 3, indicating under-dispersion in this multimodal implicit problem. For EnKF, the ensemble visualization shows that this under-dispersion corresponds to a collapse onto a single well for each coordinate. Once the Gaussian update selects the wrong well, the small ensemble spread makes subsequent recovery difficult. By contrast, EnCF and EnCF-LF are somewhat over-dispersed according to SRR, but their lower CRPS suggests that the retained ensemble diversity is useful rather than spurious. The ensemble keeps mass near both wells while still concentrating around the sign supported by the observations. This behavior is consistent with the two-valued structure induced by the circle constraint and illustrates why controlled posterior transport is beneficial for implicit many-to-one observations.

5.6 Computational cost↩︎

The improved flexibility of EnCF comes with the additional cost of learning a control during each analysis step. We therefore report wall-clock time separately from the accuracy metrics above. The comparison is meant to quantify the practical overhead of the learned controlled flow, rather than to claim that EnCF is the cheapest filter. All timings are measured under the same implementation and hardware environment. Since absolute runtime depends on implementation details, we focus on the relative cost across methods.

6 separates the transient and steady-state costs of the learned controlled flow. In the early cycles, EnCF and EnCF-LF are substantially more expensive because the control is trained with a larger epoch budget. After the cosine schedule reaches \(E_k=1\), the per-cycle cost becomes stable. EnCF-LF remains more expensive than EnCF in this steady regime because it also learns a surrogate observation energy, while EnKF is the cheapest method and the existing flow-based filters have lower per-cycle cost than the controlled-flow methods.

The right panel relates this steady-state cost to the converged RMSE on the implicit circle benchmark. The Kalman-type filters and existing flow-based filters are faster, but they have larger errors in this implicit many-to-one setting. EnCF and EnCF-LF require higher wall-clock time per assimilation cycle, but they occupy the low-error region of the plot. Thus the learned control introduces a measurable computational overhead, while providing a favorable accuracy-cost trade-off in the observation regime for which the method is designed.

6 Discussion↩︎

This work introduces implicit data assimilation as a formulation for incorporating observations that are not naturally represented by an additive residual or an explicit forward operator. By expressing the analysis distribution as an energy tilt of the forecast law, the observation model is separated from the numerical mechanism used to realize the update. EnCF implements this update through a stochastic controlled flow, while EnCF-LF replaces an unavailable analytic energy with a surrogate learned from simulator samples. The accompanying analysis establishes ideal exactness and organizes the practical filtering error into contributions from the base flow, the learned control, the surrogate energy, and temporal discretization. Together, these results provide a numerical framework for constructing ensemble analyses from energy-defined observational information.

A broader implication is that the interface between a scientific measurement model and an ensemble filter need not be restricted to residuals, moments, or pointwise likelihood weights. These representations are highly effective when their underlying assumptions are appropriate, and the results confirm that tuned Kalman-type methods remain preferable for smooth additive-Gaussian observations. The role of EnCF is therefore complementary. It is intended for problems in which the principal difficulty lies in the geometry or accessibility of the observation model, including many-to-one, multimodal, implicit, and simulator-defined relations. In these regimes, controlled transport offers a way to retain distinct observation-compatible states rather than compressing the analysis into a local Gaussian summary or requiring a separately derived guidance rule for each sensing mechanism.

The flexibility of this formulation also creates several opportunities for systematic algorithm design. Because the forecast representation, observation energy, learned control, and transport discretization enter the analysis through distinct components, each can be improved without changing the overall filtering principle. The error decomposition developed in this work makes this modularity explicit and identifies how ensemble quality, control optimization, energy approximation, and numerical resolution affect the resulting analysis. It therefore provides a useful basis for allocating computational effort across assimilation cycles and for adapting the method to sharply concentrated observations, evolving forecast distributions, and spatially coupled systems. For EnCF-LF in particular, the separation between energy learning and controlled transport allows increasingly expressive observation models to be incorporated while preserving the same analysis solver.

These observations point toward several directions for further development. Adaptive schemes could select the base flow and allocate control optimization and transport resolution according to the difficulty of each analysis step. Localized, reduced-order, or multiscale parameterizations of the control may extend the approach to substantially larger systems. For likelihood-free assimilation, richer conditional energy models could accommodate non-Gaussian, spatially dependent, and temporally correlated observation processes while retaining forward-only simulator access. More generally, the decomposition into a forecast-generating base flow, an observation energy, and a learned control provides a modular foundation on which alternative transport dynamics, energy models, and approximation strategies can be developed. This modularity, rather than the replacement of established filters in their natural regime, is the main opportunity offered by the controlled-flow perspective.

Reproducibility↩︎

All the codes for data generation, network training, and visualization are publicly available on the GitHub repository https://github.com/zylipku/EnCF.

Declarations↩︎

During the preparation of this work, the authors used large language models (LLMs) to assist with the code writing and to polish the written text for spelling and grammar. The authors thoroughly reviewed and edited the content and take full responsibility for the final publication.

Acknowledgements↩︎

Z.L. is supported by the Ministry of Education, Singapore, under its Research Centre of Excellence award to the Institute for Functional Intelligent Materials. (Project No. EDUNC-33-18-279-V12). Y.Z. is supported by the National University of Singapore and the AI for Science Institute, Beijing, through the AISI–NUS Joint Research Initiative Fund 2025 Award.

7 Deferred proofs↩︎

7.1 Base-drift Monte-Carlo error↩︎

7.2 Proof of the energy-stability estimate↩︎

7.3 Proof of the multi-step filtering theorem↩︎

8 Experimental details↩︎

8.0.0.1 EnCF configuration

For the base flow in 14 , we use \[\alpha_\tau=\epsilon_\alpha+(1-\epsilon_\alpha)\tau, \qquad \beta_\tau=\sqrt{1-(1-\epsilon_\beta)\tau}, \qquad (\epsilon_\alpha,\epsilon_\beta)=(0.5,0.025),\] following [16]. This regularization keeps the noise schedule in 13 bounded near \(\tau=1\). All flow-based methods use \(N_\tau=100\) transport steps.

The number of control-training epochs at cycle \(k\) is \[E_k= \max\left\{ 1,\, \operatorname{round}\left[ 1+\frac{99}{2} \left( 1+\cos\frac{\pi\min(k,20)}{20} \right) \right] \right\}.\] Thus, \(E_k\) decreases from \(100\) to \(1\) over the first \(20\) cycles and remains \(1\) thereafter. The \({N_\mathrm{reg}}\) regression updates associated with one rollout are evaluated in a single batched forward pass.

8.0.0.2 Circle-sensor energy

The exact conditional energy includes the state-dependent normalizer \[Z(x)= \int \exp\left[ -\frac{(x^2+y_{\mathrm c}^2-R^2)^2}{2\sigma_{\mathrm o}^2} \right] \mathop{}\!\mathrm{d}y_{\mathrm c}.\] We evaluate \(Z(x)\) independently for each coordinate using a 512-point grid over \[|y_{\mathrm c}| \le R+5\max(\sigma_{\mathrm o},10^{-3}),\] with logsumexp stabilization. Gradients of \(\log Z(x)\) are computed by automatic differentiation through the quadrature.

8.0.0.3 Baseline tuning

For the perturbed-observation EnKF and LETKF, the inflation factor is selected from \(\{1.0,1.05,1.1,1.2\}\), while the Gaspari–Cohn localization radius is fixed at \(4\). Each configuration is evaluated over five random seeds.

8.0.0.4 Numerical integration and initialization

The Kuramoto–Sivashinsky equation is integrated spectrally with \(\Delta t=0.25\). Observations are assimilated every four model steps after a 2000-step spin-up, and the initial ensemble is obtained by adding independent unit Gaussian perturbations to the reference state. The mode-separation sweep in 5.4 uses \(100\) assimilation cycles.

The double-well dynamics are integrated by Euler–Maruyama with \(\Delta t=0.05\) and ten drift substeps per model step. The ensemble is initialized across both wells using Gaussian perturbations with standard deviation \(1\), following a 1000-step spin-up.

References↩︎

---

1. Equal contributions.↩︎

2. Corresponding to Z.L. (zy.li@nus.edu.sg).↩︎