Deep Bayesian Filter for Bayes-Faithful Data Assimilation

## Abstract

Abstract State estimation for nonlinear state space models (SSMs) is a challenging task. Existing assimilation methodologies predominantly assume Gaussian posteriors on physical space, where true posteriors become inevitably non-Gaussian. We propose Deep Bayesian Filtering (DBF) for data assimilation on nonlinear SSMs. DBF constructs new latent variables h t subscript ℎ 𝑡 h_{t} italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT in addition to the original physical variables z t subscript 𝑧 𝑡 z_{t} italic_z start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT and assimilates observations o t subscript 𝑜 𝑡 o_{t} italic_o start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT . By (i) constraining the state transition on the new latent space to be linear and (ii) learning a Gaussian inverse observation operator r ⁢ ( h t | o t ) 𝑟 conditional subscript ℎ 𝑡 subscript 𝑜 𝑡 r(h_{t}|o_{t}) italic_r ( italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT | italic_o start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) , posteriors remain Gaussian. Notably, the structured design of test distributions enables an analytical formula for the recursive computation, eliminating the accumulation of Monte Carlo sampling errors across time steps. DBF trains the Gaussian inverse observation operators r ⁢ ( h t | o t ) 𝑟 conditional subscript ℎ 𝑡 subscript 𝑜 𝑡 r(h_{t}|o_{t}) italic_r ( italic_h start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT | italic_o start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) and other latent SSM parameters (e.g., dynamics matrix) by maximizing the evidence lower bound. Experiments demonstrate that DBF outperforms model-based approaches and latent assimilation methods in tasks where the true posterior distribution on physical space is significantly non-Gaussian.

## 1 Introduction

Data assimilation (DA) is a crucial technique across various scientific domains. Its primary objective is to estimate the trajectory and current state of a system by integrating an imperfect model with partially informative observations. Specifically, given a series of observations $T$ time steps $o_{1:T}$, the goal is to infer the posterior distribution of the system’s physical variables $z_{t}$: $p(z_{t}|o_{1:t})$. DA has been widely applied in fields such as weather forecasting (Hunt et al., ; Lorenc, ; Andrychowicz et al., ), ocean research analysis (Ohishi et al., ), sea surface temperature prediction (Larsen et al., ), seismic wave analysis (Alfonzo & Oliver, ), multi-sensor fusion localization (Bach & Ghil, ), and visual object tracking (Awal et al., ).

A key challenge in DA arises from the non-Gaussian nature of the posterior distributions $p(z_{t}|o_{1:t})$, which results from the inherent nonlinearity in both the system dynamics and observation models. Despite this, many operational DA systems, such as those used in weather forecasting, rely on methods like the ensemble Kalman Filter (EnKF) (Evensen, ; Bishop et al., ) for sequential state filtering (i.e., $p(z_{t}|o_{1:t})$) and the four-dimensional variational method (4D-Var) for retrospective state analysis (i.e., $p(z_{t}|o_{1:T}),t<T$). These approaches assume Gaussianity in their test distributions $q(z_{t}|o_{1:t})$ or $q(z_{t}|o_{1:T})$, a simplification driven by computational constraints.
While exact methods such as bootstrap Particle Filters (PF) or sequential Monte Carlo (SMC) (Chopin & Papaspiliopoulos, ; Daum & Huang, ; Hu & van Leeuwen, ) could compute the true posterior, their performance degrades significantly when the number of particles is insufficient (Beskos et al., ). This issue is exacerbated in high-dimensional systems, making SMC approaches impractical for many physical problems.

To address these limitations, we propose a novel variational inference approach called Deep Bayesian Filtering (DBF) for posterior estimation. Our strategy consists of two main components: (i) constraining the test distribution to remain Gaussian to ensure computational tractability, and, in cases where the original dynamics are nonlinear, (ii) leveraging a nonlinear mapping to enhance the expressive capability of the test distribution. The DBF methodology diverges into two paths depending on the nature of the system dynamics, whether linear Gaussian or nonlinear:

##### Linear dynamics

When the system’s dynamics $p(z_{t+1}|z_{t})$ are linear, DBF assumes Gaussianity in the original space, similar to traditional methods. However, DBF introduces the concept of the inverse observation operator (IOO; see also Frerix et al. ) to construct Gaussian test distributions $q(z_{t}|o_{1:t})$. The IOO, along with any unknown system parameters, are trained to minimize the Kullback-Leibler divergence between the test distribution $q(z_{t}|o_{1:t})$ and the true posterior $p(z_{t}|o_{1:t})$. The IOO and the system parameters are trained without teacher signals $z_{t}$.

##### Nonlinear dynamics

In the more common case of nonlinear dynamics, DBF operates in a latent space, assuming Gaussianity in the latent variables $h_{t}$. The original physical variables are recovered through a nonlinear mapping function $\phi$, implemented via neural networks (NNs). This nonlinear mapping allows for a more flexible representation of the test distribution $q(z_{t}|o_{1:t})$. The IOO and other parameters are trained in a supervised manner (i.e., $z_{t}$ is used during training).

For state space models (SSMs) with nonlinear dynamics, DBF functions as a variational autoencoder (VAE) that adheres to the Markov property. Posterior distributions of the latent variables $h_{t}$ are expressed in a Bayesian framework. This approach is closely related to dynamical VAEs (DVAEs, Girin et al. for a review), which use VAEs to model time-series data. However, DBF distinguishes itself by its posterior design. Unlike DVAEs, where Monte Carlo sampling is required for inference (see Sec. 2.6.1), DBF allows for the analytical computation of the prediction step, recursively computing posteriors through closed-form expressions.

When applied to problems with nonlinear or unknown dynamics, DBF can be interpreted as learning the Koopman operator (Koopman, ) using NNs. The discovery of such latent spaces and operators through machine learning has been extensively studied (Takeishi et al., ; Lusch et al., ; Azencot et al., ) and will be experimentally validated through the handling of nonlinear filtering tasks involving chaotic dynamics.

Key contributions of the proposed DBF methodology include:

- •
DBF is the first VAE-based model for time-series data that maintains a posterior structure faithful to the Markov property in SSMs.
- •
For systems with linear dynamics, DBF extends the Kalman Filter (KF) to handle nonlinear observations through learnable NNs. The training process enables the model to infer unknown system parameters directly from data (see Section 3.1).
- •
For nonlinear dynamics, DBF constructs a new latent space for data assimilation, allowing for the analytical integration of time steps and preventing the accumulation of Monte Carlo sampling errors. This is accomplished through the application of Koopman operator theory, which ensures that the model’s representational power is maintained, as long as the latent space is sufficiently high-dimensional (see Sections 3.2 and 3.3).
- •
As a generative model, DBF estimates the uncertainty of the physical variables $z_{t}$, in contrast to regression methods that yield only point estimates (see Section 3.2 and Fig. 3).
- •
The linear constraint on dynamics stabilizes the training process, which is known to be unstable in standard recurrent NNs (see Section 3.3 and Fig. 6).

DBF has demonstrated superior performance over classical DA algorithms and latent assimilation methods in scenarios with highly non-Gaussian posteriors, particularly in the presence of strongly nonlinear observation operators or large observation noise.

## 2 Method

### 2.1 Inference of physical variables in a state-space model

A physical system is defined by variables $z_{t}$, with its evolution described by the dynamics model $p(z_{t+1}|z_{t})=\mathcal{N}(z_{t+1};f(z_{t}),Q)$, where $\mathcal{N}(x|\mu,\Sigma)$ denotes a Gaussian whose mean and covariance are $\mu$ and $\Sigma$. The nonlinear function $f$ is the dynamics operator and $Q$ is the system covariance. The Markov property holds, as $z_{t+1}$ depends only on $z_{t}$. An observation model $p(o_{t}|z_{t})=\mathcal{N}(o_{t};h(z_{t}),R)$ relates observations to physical variables via the observation operator $h$ and covariance $R$. The panel (a) of Fig. 1 shows the system’s graphical model. The objective of sequential DA is to compute the posterior of $z_{t}$ given $o_{1:t}$.

### 2.2 KF for linear dynamics, linear observations

In the KF, the dynamics and observation models are both linear Gaussian, with the mean being linear with respect to the conditional variable and a constant covariance matrix.
Given that the dynamics and observation operators $f,h$ are linear, we can represent them using matrices $A$ and $C$, respectively. All matrices ($A,C,Q$, and $R$) are constant. The filter distribution $p(z_{t}|o_{1:t})$ remains Gaussian, provided that the initial distribution $p(z_{1})$ is Gaussian. We can recursively compute the posterior parameters (means $\mu_{t}$ and covariance matrices $\Sigma_{t}$) using the following equations:

$$ $\displaystyle\mu_{t}$ $\displaystyle=$ $\displaystyle\Sigma_{t}(A\Sigma_{t-1}A^{T}+Q)^{-1}A\mu_{t-1}+K_{t}(o_{t}-H\mu_ {t}),$ (1) $\displaystyle\Sigma_{t}^{-1}$ $\displaystyle=$ $\displaystyle(A\Sigma_{t-1}A^{T}+Q)^{-1}+HR^{-1}H^{T},$ (2) $$

where $K_{t}=(A\Sigma_{t-1}A^{T}+Q)H^{T}(H(A\Sigma_{t-1}A^{T}+Q)H^{T}+R^{-1})^{-1}$ is the Kalman Gain.

### 2.3 DBF for linear dynamics, nonlinear observations

In this scenario, Gaussianity of the test distribution is lost during the KF update step. We introduce an inverse observation operator (IOO) $r(z_{t}|o_{t})$ (see also Frerix et al. ):

$$ $\displaystyle p(z_{t}|o_{1:t})$ $\displaystyle=$ $\displaystyle\frac{p(o_{t}|z_{t})p(z_{t}|o_{1:t-1})}{p(o_{t}|o_{1:t-1})} \propto\frac{r(z_{t}|o_{t})}{\rho(z_{t})}p(z_{t}|o_{1:t-1}),$ (3) $$

where $r(z_{t}|o_{t})=\frac{p(o_{t}|z_{t})\rho(z_{t})}{\int p(o_{t}|z_{t})\rho(z_{t})
dz_{t}}$ and $\rho(z_{t})$ is a virtual prior. By approximating both the IOO and the virtual prior as Gaussians, $r(z_{t}|o_{t})=\mathcal{N}(f_{\theta}(o_{t}),G_{\theta}(o_{t}))$ and $\rho(z_{t})=\mathcal{N}(m,V)$, respectively, the posterior $q(z_{t}|o_{1:t})$ can be analytically computed as a Gaussian, where the mean $\mu_{t}$ and covariance $\Sigma_{t}$ are given as:

$$ $\displaystyle\mu_{t}$ $\displaystyle=$ $\displaystyle\Sigma_{t}(A\Sigma_{t-1}A^{T}+Q)^{-1}A\mu_{t-1}+G_{\theta}(o_{t}) ^{-1}f_{\theta}(o_{t})-V^{-1}m,$ (4) $\displaystyle\Sigma_{t}^{-1}$ $\displaystyle=$ $\displaystyle(A\Sigma_{t-1}A^{T}+Q)^{-1}+G_{\theta}(o_{t})^{-1}-V^{-1},$ (5) $$

where $f_{\theta}(o_{t})$ and $G_{\theta}(o_{t})$ are NNs with parameters $\theta$, and $m$ and $V$ are constants set to $m=0$ and $V=10^{8}I$. These values bias the NNs’ outputs without affecting performance.

The recursive formula for the exact posterior (Equation 3) requires no approximation. Thus, DBF computes the exact posterior when the true IOO $r_{\rm true}(h_{t}|o_{t})$ is Gaussian, i.e., the SSM is a LGSS. In that case, the posterior update formula agrees with the KF (see Equations 1, 2 and 4, 5).
The key difference is that nonlinear functions are applied to both the mean, $f_{\theta}(o_{t})$, and the covariance, $G_{\theta}(o_{t})$. In the KF, $f_{\theta}(o_{t})$ is linear, and $G_{\theta}(o_{t})$ is a constant covariance matrix (see Equations 1 and 2).
The dependence of $G_{\theta}(o_{t})$ on observations provides flexibility in adjusting the impact of the new observation on the state estimation. The importance of adjusting the internal state updates based on observations has also been discussed in recent SSM-based approaches (Gu & Dao, ).

Figure: Figure 1: Panel (a) shows the graphical model for the simplest SSM. If the dynamics of the original SSM is linear, DBF assimilates on that space. Panel (b) shows the graphical model for the SSM assumed for SSM with nonlinear dynamics. Panel (c) shows the inference structure of our methodology for SSM with nonlinear dynamics.
Refer to caption: x1.png

### 2.4 DBF for nonlinear dynamics, linear/nonlinear observations

In this scenario, the Gaussianity of the test distribution is lost during the predict step, making it impossible to apply the original dynamics over the physical variables $z_{t}$. Therefore, we introduce a new set of latent variables $h_{t}$ and assume a dynamics model over $h_{t}$: $p(h_{t+1}|h_{t})=\mathcal{N}(h_{t+1}|Ah_{t},Q)$ (see panel (b) in Fig. 1). The IOO maps observations into the latent variables $h_{t}$: $r(h_{t}|o_{t})$. The recursive formula follows Equations 4 and 5.
To retrieve the distribution of the original physical variables $z_{t}$, we introduce an emission model $p(z_{t}|h_{t})=\mathcal{N}(z_{t};\phi(h_{t}),R)$, where $\phi$ is represented by a NN. By marginalizing over $h_{t}$ with this emission model, a trained DBF can generate samples of $z_{t}$ that follow the test distribution $q(z_{t}|o_{1:t})$ given observations $o_{1:t}$.

Although the dynamics operator $A$ for the latent variables $h_{t}$ is linear, it is able to express any nonlinear dynamics if the latent space is sufficiently high-dimensional. The Koopman operator (Koopman, ) provides a framework for representing nonlinear systems by mapping observables—functions of the system’s state—into a higher-dimensional space where the dynamics become linear. Mathematically, for a dynamical system $z_{t+1}=f(z_{t})$, the Koopman operator $\mathcal{K}$ is a linear operator acting on a set of observables $g(z)$, such that $\mathcal{K}g(z_{t})=g(f(z_{t}))$. This allows the system to be represented as $h_{t+1}=Ah_{t}$ in the latent space $h_{t}$, where $A$ is the dynamics matrix learned by DBF. While the physical dynamics $f(z)$ are nonlinear, the Koopman operator ensures the existence of an embedding that linearizes the dynamics, enabling recursive computation of test distributions. Discovering such embeddings in finite dimensions has been widely studied (Takeishi et al., ; Lusch et al., ; Azencot et al., ).
In high-dimensional simulations, the true degrees of freedom are often far fewer than the simulated variables, making surrogate modeling with the Koopman operator a promising approach to reducing computational costs.

### 2.5 Training

When assimilating in the physical space (i.e., when the dynamics are linear), we train the IOO (i.e., $f_{\theta}$ and $G_{\theta}$) by optimizing the evidence lower bound (ELBO) without using the teacher signal $z_{t}$:

$$ $\displaystyle\log p(o_{1:T})$ $\displaystyle=$ $\displaystyle\sum_{t=1}^{T}\log p(o_{t}|o_{1:t-1})\geq-\mathcal{L}_{\rm ELBO},$ $\displaystyle\mathcal{L}_{\rm ELBO}$ $\displaystyle=$ $\displaystyle-\sum_{t=1}^{T}\int q(h_{t}|o_{1:t})\log p(o_{t}|h_{t})dh_{t}+KL[ q(h_{t}|o_{1:t})||q(h_{t}|o_{1:t-1})],$ (6) $$

where $KL[p||q]$ denotes the Kullback-Leibler divergence between distributions $p$ and $q$. If the SSM contains any unknown parameters, we can train these parameters as well.

For SSMs with nonlinear or unknown dynamics, we have two approaches:

##### Strategy 1

Pretrain the Koopman operator, which consists of the nonlinear mapping from $z_{t}$ to $h_{t}$, the linear dynamics between $h_{t}$ and $h_{t+1}$ represented by matrix $A$, and the reverse nonlinear mapping from $h_{t}$ to $z_{t}$ denoted by $\phi$. With these components ($A$ and $\phi$) of the Koopman operator, the method designed for linear dynamics can be applied. For pretraining, we require samples of $z_{t}$ or the SSM for the physical variables to generate these samples. Pairs of $z_{t}$ and $o_{t}$ are not necessary, as the training for the linear dynamics ($A$ and $\phi$) and the IOO ($r(h_{t}|o_{t})$) can be performed separately.

##### Strategy 2

Train all components (the matrix $A$, the stochastic mapping $p(z_{t}|h_{t})=\mathcal{N}(z_{t};\phi(h_{t}),{\rm diag}[\sigma^{2}])$, and the IOO) simultaneously. In this case, samples of ($z_{t}$, $o_{t}$) pairs or the SSM for both physical and observation variables to generate these sample pairs are required. The parameters are optimized by maximizing a joint ELBO, $\mathcal{L}_{\rm ELBO,{\rm joint}}$, via supervised training:

$$ $\displaystyle\log p(o_{1:T},z_{1:T})$ $\displaystyle=$ $\displaystyle\sum_{t=1}^{T}\log p(o_{t},z_{t}|o_{1:t-1},z_{1:t-1})\geq- \mathcal{L}_{\rm ELBO,{\rm joint}},\hskip 99.58464pt$ $\displaystyle\mathcal{L}_{\rm ELBO,{\rm joint}}$ $\displaystyle=$ $\displaystyle-\sum_{t}\int q(h_{t}|o_{1:t})\log p(z_{t}|h_{t})dh_{t}+KL[q(h_{t }|o_{1:t})||q(h_{t}|o_{1:t-1})].$ (7) $$

Here, we have neglected $\log p(o_{t}|z_{t})$ from $\log p(o_{t},z_{t}|h_{t})=\log p(o_{t}|h_{t},z_{t})+\log p(z_{t}|h_{t})=\log p
(o_{t}|z_{t})+\log p(z_{t}|h_{t})$. Since $\log p(o_{t}|z_{t})$ does not depend on $h_{t}$ and its distributions, it is sufficient to include $\log p(z_{t}|h_{t})$. We have replaced $q(h_{t}|o_{1:t},z_{1:t})$ with its special case $q(h_{t}|o_{1:t})$ as our objective is to give the best estimate of $z_{t}$ given observations $o_{1:t}$.

### 2.6 Related works

We discuss related works in the order of relevance.

#### 2.6.1 Dynamical variational autoencoders

DVAEs (see Girin et al. for a review) are models closely related to DBF for nonlinear dynamics. Both incorporate time-series architecture in a VAE, but there are two key differences: (i) the posterior design and realization of the dynamics step, and (ii) the loss function.

##### posterior design

Our strategy for the test distribution is to incorporate an appropriate architecture that reflects the Markov property in the time dimension of the test distribution. The IOO, $r(h_{t}|o_{t})$, and the linear dynamics model serve as key instruments in constructing the test posterior distributions. A distinguishing feature of our methodology is that each component’s role is defined with respect to the Markov property of the state-space model (SSM) and is clearly differentiated from other components involved in posterior construction. For example, the IOO influences only the update step and does not affect the prediction step. We refer to this methodology as ”Bayes-Faithful” due to its tailored design for SSMs that exhibit the Markov property.

In contrast, the test posterior distributions in DVAEs are constructed using RNNs. The complexity of the transition model prevents the analytical computation of latent variables across time steps. As a result, these values can only be estimated via Monte Carlo sampling. Consequently, during inference, successive Monte Carlo sampling (“cascade trick”; Girin et al. ) becomes unavoidable.

##### loss function

DBF takes the ELBO from factorized density $\log p(o_{t}|o_{1:t-1})$ in $\log p(o_{1:T})=\sum_{t}\log p(o_{t}|o_{1:t-1})$:

$$ $\displaystyle\log p(o_{1:T})$ $\displaystyle\geq$ $\displaystyle\sum_{t=1}^{T}(E_{q(h_{t}|o_{1:t})}[\log p(o_{t}|h_{t})]-KL[q(h_{ t}|o_{1:t})|q(h_{t}|o_{1:t-1})]).$ (8) $$

On the other hand, DVAEs take the ELBO from probability density with all the observations at once.

$$ $\displaystyle\log p(o_{1:T})$ $\displaystyle\geq$ $\displaystyle E_{q(h_{1:T}|o_{1:T})}[\log p(h_{1:T},o_{1:T})-\log q(h_{1:T}|o_ {1:T})].$ (9) $$

Therefore, DBF seeks for the filtered distributions $q(h_{t}|o_{1:t})$ whereas DVAEs model the smoother distributions $q(h_{t}|o_{1:T})$. Again, for DVAEs, to evaluate the expected values in Equation 9, we need to undergo successive Monte-Carlo sampling over $T$ variables ($h_{1:T}$).

Assuming linear Gaussian dynamics and a Gaussian IOO, DBF allows for the analytical integration of $q(h_{t}|o_{1:t-1})$, resulting in a structured encoder. This structured posterior enables the recursive computation of the filtered distribution $q(h_{t}|o_{1:t})$ without relying on Monte Carlo sampling, setting it apart from other DVAEs. By constraining the dynamics to be linear, DBF ensures exact integration without the accumulation of Monte Carlo sampling errors across time steps.

Moreover, the linear assumption helps DBF mitigate the instability issues commonly faced when training standard RNNs. SSMs are increasingly favored for modeling long-range dependencies (Gu & Dao, ). S4 (Gu et al., ) utilizes linear dynamics in the latent space and learns the dynamics, proposing an efficient computation algorithm that outperforms transformers on datasets with long-range dependencies. LS4 (Zhou et al., ) extends S4 by introducing stochasticity through a VAE-like structure.
Both LS4 and DBF employ linear SSMs and Gaussian posterior approximations, but DBF updates the mean and covariance using a recursive formula based on Bayes’ rule, while LS4 replaces recurrence with convolutions, forgoing a recursive approach.

#### 2.6.2 KF-based methods

Various approaches have been explored to address LGSS limitations, including linearizing the model via first-order approximations like the extended Kalman Filter (EKF), approximating populations with a Gaussian distribution in the ensemble Kalman Filter (EnKF; Evensen ), and using NNs to approximate the Kalman gain (Revach et al., ). The EnKF and its variants (e.g., ETKF; Bishop et al. ) are commonly used in real-time data assimilation for weather forecasting. However, these methods rely on the KF’s posterior update equations, limiting the expressivity of the distributions they can represent. Additionally, computations for covariance matrices become challenging in high-dimensional spaces, requiring specialized techniques for computational efficiency.

#### 2.6.3 Sampling-based methods

The Particle Filter is a popular method for assimilating any posterior. However, achieving adequate particle density in high-dimensional state spaces poses significant challenges. Insufficient density of particles leads to particle degeneracy, where few particles explain the observed data (Beskos et al., ). In contrast, DBF directly learns to position density through the IOO, offering advantages for high-dimensional tasks. The Particle Flow Filter (PFF; Daum & Huang ; Hu & van Leeuwen ) addresses particle degeneracy by moving particles according to gradient flow and effectively scales to nonlinear SSMs with hidden state dimensions up to 1000 (Hu & van Leeuwen, ).

#### 2.6.4 Approximate MAP estimation method

MAP estimation is used to identify the high-density point of the posterior in high-dimensional space, such as in weather forecasting Lorenc (); Frerix et al. (). Even if the computation of the posterior $p(h_{t}|o_{1:t})$ is intractable, we can optimize $\log p(h_{t}|o_{1:t})=\log p(o_{t}|h_{t})+\log p(h_{t}|o_{1:t-1})$ if we can describe $p(o_{t}|h_{t})$ and $p(h_{t}|o_{1:t-1})=\int p(h_{t}|h_{t-1})p(h_{t-1}|o_{1:t-1})dh_{t-1}$ explicitly. In practice, we cannot access $p(h_{t-1}|o_{1:t-1})$ and therefore the integral $\int p(h_{t}|h_{t-1})p(h_{t-1}|o_{1:t-1})dh_{t-1}$, so we only compute the mean. The downside is that sequential computation of the covariance matrix of $p(h_{t}|o_{1:t-1})$ is impossible.

#### 2.6.5 NN-based PDE surrogate

Recently, there have been attempts to approximate partial differential equations (PDEs) using NNs. We also tried one of the latest methods, PDE-refiner  (Lippe et al., ), in the experiments of this study. However, its performance was poor, and we decided not to include it in the experiments. We suspect this is because PDE-refiner was designed for constructing PDE surrogates and did not account for noisy observations, making it susceptible to noise. We confirmed that it produces good predictions under noiseless observations.

## 3 Experiments

We evaluate the performance of DBF on three tasks: a linear dynamics problem (moving MNIST) and two nonlinear dynamics problems (double pendulum and Lorenz96).

##### Linear dynamics: moving MNIST

In the moving MNIST task, the goal is to identify the images, positions, and velocities of two handwritten digits as they move within the observed frames. While the dynamics of these digit images and their observation processes are provided, the actual images, positions, and velocities are not available, making supervised learning impossible.
DBF assimilates directly in the physical space. It trains the parameters of the IOO, $f_{\theta}$ and $G_{\theta}$, which are represented by NNs, along with the pixel values of the embedded images. The results are compared against conventional DA methods (EnKF, ETKF, PF) that also operate in the physical space.

##### Nonlinear dynamics: double pendulum and Lorenz96

For nonlinear dynamics problems, such as the double pendulum and Lorenz96, DBF constructs a new latent space in addition to the original physical space. Here, we took Strategy 2 in Sec.2.5 for the training: we simultaneously train NNs for the IOO, nonlinear observation operator $\phi$, the dynamics matrix $A$, and the emission model’s standard deviation.
We compare the performance of DBF with the classical DA algorithms (EnKF, ETKF, PF), state-of-the-art assimilation methodologies (PFF Daum & Huang ; Hu & van Leeuwen , KalmanNet Revach et al. ), and DVAE-based approaches (deep Kalman Filter; DKF, Krishnan et al. ; , variational recurrent neural network; VRNN, Chung et al. , and stochastic recurrent neural network; SRNN, Fraccaro et al. ).
DBF and other DVAEs are trained by optimizing the evidence lower bound (ELBO), as described in Sec. 2.5.

For all the experiments, we generate random initial conditions and let them evolve via the dynamics of the problem. Synthetic observation data are generated by applying the observation operator and additive noise. Noise levels and observation operators are further explained in the section for each problem and in the appendix. Training details, including hyperparameters and architectures, are also provided in the appendix.

### 3.1 Linear dynamics: two-body moving MNIST

This experiment demonstrates DBF’s ability to handle linear dynamics where key parameters of the observation operator are unknown. The dataset consists of 2D figures containing two embedded images, each moving at a constant speed and bouncing off frame edges. The system’s physical state is described by eight variables: the positions $(x,y)$ and velocities $(v_{x},v_{y})$ of the two embedded images. The dynamics matrix is block-diagonal, composed of four (two-body times two dimensions) translation matrices, $A_{tr}$:
$A_{tr}=\begin{pmatrix}1&1\\
0&1\\
\end{pmatrix},z_{t}=\begin{pmatrix}x_{t}\\
v_{x_{t}}\\
\end{pmatrix}$.
Observations are corrupted by additive Gaussian noise with a standard deviation of $\sigma=50$ per pixel, where the original pixel values range from 0 0 to $255$ (see panel (a) of Fig. 2 for an example of the data provided).

The aim is to show that DBF can track the linear dynamics while estimating unknown system parameters. DBF learns the pixel values of the embedded images from noisy observations, while maintaining consistency with physical motion. The observation model contains $1{,}568$ unknown parameters, corresponding to the number of pixels in the images.
In classical DA algorithms, it is not possible to train unknown system parameters. However, it may be possible to infer these parameters by incorporating them as new physical dimensions. We have adopted this strategy for classical DA algorithms (EnKF, ETKF, and PF). We tested at three different noise levels ($\sigma_{sys}$ of $1$, $0.1$, and $0.01$). DVAEs are not comparable as they need to undergo supervised training. We were unable to compare with KalmanNet, as its input dimensionality, $x_{\rm dim}^{2}=(44\times 44)^{2}$, exceeded the available GPU memory.

Fig. 2 summarizes the experiment. Panel (a) shows an example from the test set, illustrating the challenges posed by strong noise and overlapping images. Panel (b) presents the DBF learning process. In the rightmost table, we compare the success rates of DBF against model-based approaches (EnKF, ETKF, PF). We define success as achieving a root-mean-square error (RMSE) of less than 1.0 for both position ($x_{1},y_{1},x_{2},y_{2}$) and velocity ($v_{x_{1}},v_{y_{1}},v_{x_{2}},v_{y_{2}}$) of the two digits over the final ten steps.
DBF successfully performs assimilation without explicit knowledge of the images, while all the other model-based approaches fail. The KF-inspired approaches (EnKF, ETKF) failed because of very strong non-Gaussianity in the observation process and the high system dimension. Similarly, PF underperformed because the number of particles ($10{,}000$) was insufficient for the problem dimension ($z_{\rm dim}=8$ and two digits images $2\times 28\times 28=1{,}568$). Figures for visualizing the assimilation results for all the algorithms are given in the appendix (Fig. 7).

Panel (b) of Fig. 2 illustrates the evolution of the estimated figures. Initially, DBF assumes two random shapes. As training progresses, it first identifies one of the numbers (“9”) and subsequently detects the second shape (“5”). By the end of the training process, DBF nearly perfectly estimates the parameters of the observation model, including the positions of the figures, which is crucial for adjusting their reflective behavior. The ability to learn system parameters through gradient descent sets DBF apart from traditional model-based approaches.

Figure: Figure 2: Figures from the two-body Moving MNIST experiments. Panel (a) displays examples of the observation data. Panel (b) illustrates the evolution of the observation model parameters (the embedded images) during training. The table compares the success rates of four methodologies.
Refer to caption: x4.png

### 3.2 Nonlinear dynamics 1: double pendulum

Figure: Figure 3: A schematic figure (panel a) and results for double pendulum experiments. Panel (b) shows the RMSE evolution during assimilation. Panels (c) and (d) show example histograms for normalized errors in DBF and ETKF samples compared against the unit Gaussian $\mathcal{N}(x;\mu=0,\sigma^{2}=1)$. The small table compares the Jeffreys divergence of normalized errors and the unit Gaussian between DBF, EnKF, and ETKF predictions.
Refer to caption: x6.png

**Table 1: RMSE at the final ten steps of assimilation in double pendulum experiments.**
|  | $\sigma=0.1$ | $\sigma=0.3$ | $\sigma=0.5$ |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | $\theta$ | $\omega$ | $\theta$ | $\omega$ | $\theta$ | $\omega$ |
| DBF | 0.03 $\pm$ 0.01 | 0.21 $\pm$ 0.04 | 0.05 $\pm$ 0.02 | 0.26 $\pm$ 0.05 | 0.06 $\pm$ 0.01 | 0.36 $\pm$ 0.04 |
| EnKF | $0.05\pm 0.00$ | $0.33\pm 0.07$ | $0.14\pm 0.01$ | $0.71\pm 0.09$ | $0.24\pm 0.01$ | $1.17\pm 0.22$ |
| ETKF | $0.05\pm 0.01$ | $0.46\pm 0.08$ | $0.22\pm 0.05$ | $1.41\pm 0.41$ | $0.36\pm 0.08$ | $2.70\pm 1.25$ |
| PF | $0.05\pm 0.00$ | $0.63\pm 0.24$ | $0.21\pm 0.14$ | $1.41\pm 1.30$ | $0.32\pm 0.08$ | $2.36\pm 2.29$ |
| PFF | $1.27\pm 0.29$ | $1.04\pm 0.15$ | NA | $5.99\pm 1.09$ | $5.88\pm 0.67$ | NA |
| KNet | NA | NA | NA | NA | NA | NA |
| VRNN | $0.04\pm 0.01$ | $0.44\pm 0.19$ | $0.06\pm 0.02$ | $0.35\pm 0.14$ | $0.08\pm 0.04$ | $0.40\pm 0.16$ |
| SRNN | $0.05\pm 0.02$ | $0.52\pm 0.18$ | $0.06\pm 0.02$ | $0.44\pm 0.08$ | $0.08\pm 0.03$ | $0.52\pm 0.22$ |
| DKF | $0.12\pm 0.02$ | $2.70\pm 0.28$ | $0.17\pm 0.03$ | $2.61\pm 0.74$ | $0.23\pm 0.04$ | $2.61\pm 0.56$ |

This section presents our experiments with a double pendulum system, selected for its nonlinear and chaotic behavior. The pendulum consists of two 1 kg masses, P1 and P2, connected by two 1 meter bars, B1 and B2. One end of the bar B1 is fixed at the origin (“O”), with the other end attached to P1. Mass P2 is connected to P1 via bar B2. A schematic of the setup is shown in panel (a) of Fig. 3.

We use the angles $\theta_{1}$ and $\theta_{2}$, and the two angular velocities, $\omega_{1}$ and $\omega_{2}$, as target physical variables. The latent dimension for DBF, VRNN, SRNN, and DKF is set to 50. Observation data consists of the two-dimensional spatial positions of masses $P_{1}$ and $P_{2}$, corrupted by Gaussian noise. The observation operator combines trigonometric functions for $\theta_{1}$ and $\theta_{2}$, creating a highly nonlinear relationship. Experiments are conducted with noise levels of $\sigma=0.1,0.3$, and $0.5$ meters, with a time step of 0.03 seconds between observations. In the emission model $p(z_{t}|h_{t})$, we assume von Mises distributions for $\theta_{1}$ and $\theta_{2}$, while $\omega_{1}$ and $\omega_{2}$ follow Gaussian distributions.

Table 1 presents the RMSE between the physical variables and the mean of the filtered distribution. Training for KalmanNet was unsuccessful under all conditions. For the DVAEs, we exclude failed initial conditions when calculating the RMSE. DBF outperforms both model-based and latent assimilation methods across all settings, showing significant improvements in estimating $\omega$, which cannot be inferred from a single observation. Fig. 3 (b) illustrates an example of RMSE evolution during assimilation, where DBF consistently outperforms the other methods. The assimilation of $\omega$ occurs within the first $\sim 20$ steps, maintaining an excellent estimation accuracy throughout the experiment.

A key feature of DBF is its ability to generate samples of $z_{t}$ and assess the uncertainty in state estimates. To evaluate this capability, we analyze the distributions of normalized errors defined as $\epsilon_{norm,t,i}=(z_{t,sample,i}-z_{t,i})/\delta_{i}$, where $z_{t,i}$ represents the true value of dimension $i$ at time $t$, and $\delta_{i}$ is the standard deviation of $z_{t,sample,i}$. We collect $\epsilon_{norm,t,i}$ across all time steps, focusing on $i=\omega_{1}$ and $i=\omega_{2}$, since $\theta_{1}$ and $\theta_{2}$ follow von Mises distributions. If the uncertainty estimates are accurate, $\epsilon_{norm,t,i}$ should approximate a Gaussian distribution with a standard deviation of one. To quantify the accuracy, we compute the symmetric KL divergence (Jeffreys divergence) $KL_{sym}[p,q]=(KL[p||q]+KL[q||p])/2$ between the histogram of $\epsilon_{norm,t,i}$ and a unit Gaussian. DBF exhibits very low $KL_{sym}$ values, indicating accurate error estimation. Panels (c) and (d) display example histograms of $\epsilon_{norm,t,i}$ for DBF and ETKF.

### 3.3 Nonlinear dynamics 2: Lorenz96

Figure: Figure 4: A Hovmöller diagram for one of data in the test set. The observation operator is nonlinear, $o_{t,j}=min(z_{t,j}^{4},10)+\epsilon$.
Refer to caption: x10.png

In the final experiment, we focus on state estimation in the Lorenz96 model (Lorenz, ), a benchmark for testing data assimilation algorithms on noisy, nonlinear observations.
The Lorenz96 model describes the evolution of a one-dimensional array of variables, each representing a physical quantity over a spatial domain, like an equilatitude circle. The dynamics are governed by the following coupled ordinary differential equations:

$$ $\displaystyle\frac{dz_{i}}{dt}=(z_{i+1}-z_{i-2})z_{i-1}-z_{i}+F,\hskip 8.53581 pti=1,\ldots,N,$ (10) $$

where $z_{i}$ is the value at grid point $i$, $N$ is the number of grid points, and $F$ is external forcing. For our experiments, we take $F=8$ and $N=40$.

We consider two observation operators. The first adds Gaussian noise to direct observations: $o_{t,j}=z_{t,j}+\epsilon$, with noise levels $\sigma=1,3,5$. The second uses a nonlinear operator: $o_{t,j}=\min(z_{t,j}^{4},10)+\epsilon$, with the same noise levels. The dynamic range of $z_{t,j}$ is around $\pm 10$, and observations are capped at 10 when $z_{t,j}$ exceeds 1.8. This makes it highly challenging for classical DA methods, as each observation offers limited information. The filter must integrate data over long timesteps, where nonlinear dynamics distort the probability distribution. Fig. 4 illustrates observations and target values. All models use 80 observation steps with a 0.03 time interval. The latent dimension for DBF, VRNN, SRNN, and DKF is set to 800. For further details for the experiment, see Sec. A.3.

**Table 2: RMSE at the final ten steps of assimilation in Lorenz96 experiments.**
|  | direct observation | nonlinear observation |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | $\sigma=1$ | $\sigma=3$ | $\sigma=5$ | $\sigma=1$ | $\sigma=3$ | $\sigma=5$ |
| DBF | $0.53\pm 0.04$ | 0.82 $\pm$ 0.03 | 1.16 $\pm$ 0.07 | 1.08$\pm$ 0.15 | 1.29 $\pm$ 0.18 | 1.65 $\pm$ 0.17 |
| EnKF | $0.31\pm 0.01$ | $0.83\pm 0.10$ | $1.73\pm 0.12$ | $4.69\pm 0.14$ | $3.93\pm 0.08$ | $3.81\pm 0.07$ |
| ETKF | 0.30 $\pm$ 0.01 | $1.06\pm 0.15$ | $2.42\pm 0.11$ | $4.57\pm 0.25$ | $4.28\pm 0.04$ | $4.23\pm 0.07$ |
| PF | $2.80\pm 0.04$ | $3.12\pm 0.06$ | $3.62\pm 0.13$ | $6.05\pm 0.16$ | $4.95\pm 0.12$ | $4.58\pm 0.14$ |
| PFF | $0.60\pm 0.02$ | $1.00\pm 0.05$ | $2.20\pm 0.09$ | $3.75\pm 0.09$ | $3.85\pm 0.04$ | $3.83\pm 0.11$ |
| KNet | $0.60\pm 0.02$ | $1.81\pm 0.05$ | $3.02\pm 0.09$ | $2.97\pm 0.21$ | $3.47\pm 0.17$ | $3.99\pm 0.25$ |
| VRNN | $3.67\pm 0.06$ | $3.67\pm 0.06$ | $3.67\pm 0.06$ | $3.69\pm 0.04$ | $2.51\pm 0.79$ | $3.67\pm 0.06$ |
| SRNN | $3.08\pm 0.56$ | $3.63\pm 0.05$ | $3.40\pm 0.29$ | $3.30\pm 0.81$ | $3.62\pm 0.41$ | $2.96\pm 0.32$ |
| DKF | 3.70 | NA | NA | NA | NA | NA |

Figure: Figure 5: RMSE results for Lorenz96 experiments. Panels (a), (b) show results for direct observation with $\sigma=1$ and $\sigma=5$. Panel (c) shows results for nonlinear observation with $\sigma=1$.
Refer to caption: x11.png

Table 2 presents the assimilation performance across different noise levels and observation settings. DBF outperforms existing methods in direct observations with $\sigma=3,5$, and across all noise levels for nonlinear observation cases. In the $\sigma=1$ setting with direct observation, traditional algorithms like EnKF and ETKF outperform DBF.

The superior performance of EnKF and ETKF with direct observations at the lowest noise level can be attributed to the minimal non-Gaussianity in the posteriors within physical space. Non-Gaussianity can originate from both the dynamics model (predict step) and the observation model (update step). In this setting, the linearity of the observation operator prevents non-Gaussianity from being introduced during the update step, provided that the prior $q(z_{t}|o_{1:t-1})$ is Gaussian. Additionally, state estimation from each observation is highly accurate due to small noise. As a result, the prior $p(z_{t+1}|o_{1:t})$ remains close to a Gaussian distribution, as the locally linear approximation of the dynamics adequately captures the time evolution of probability distributions. The poorer performance of EnKF and ETKF in the $\sigma=5$ experiment is attributed to the increased non-Gaussianity introduced during each predict step. Similarly, when the observation operator is nonlinear, each update step introduces substantial non-Gaussianity. This results in a significant drop in performance for traditional filtering methods across all noise levels. In these scenarios, DBF consistently maintains an advantage over classical DA algorithms.

Figure: Figure 6: Histogram of 800 eigenvalues of the dynamics matrix in Lorenz96. D for direct and N for nonlinear observations.
Refer to caption: x14.png

We observe that training DVAE-based methods is highly unstable, while that for DBF exhibits stability. Dynamics in DVAEs are modeled by RNNs, which often suffer from unstable training due to exploding or vanishing gradients. In contrast, DBF employs matrix multiplication for dynamics. If the eigenvalues of the matrix exceed one by a large margin, the model predictions, and consequently the loss function, would explode irrespective of inputs.
Fig. 6 shows the histogram of the absolute values of eigenvalues at the end of training, which are distributed around or below one, indicating stable training.

## 4 Limitation

DBF’s learning of IOO requires a training phase, unlike classical model-based data assimilation methods. Specifically, when dealing with nonlinear dynamics, DBF requires either: (i) a pair of $(z_{t},o_{t})$ generated from the original SSM, (ii) a pair of $(z_{t},o_{t})$ obtained via, e.g., retrospective reanalysis (ERA5; Hersbach et al. in weather forecasting), or (iii) a pretrained Koopman operator and observed data $o_{t}$.

In the Lorenz96 experiment, DBF’s performance with direct observation with $\sigma=1$ falls short compared to EnKF and ETKF. In this setting, the non-Gaussianity of posteriors is weak, resulting in minor approximation errors due to Gaussian assumptions. Consequently, a model-based approach may be more advantageous in such situations, as it leverages complete SSM knowledge without introducing training biases.

## 5 Conclusion

We propose DBF, a novel DA method. DBF is a NN-based extension of the KF designed to handle nonlinear observations. While constraining the test distributions to remain Gaussian, DBF enhances their representational capacity by leveraging nonlinear transform expressed by a NN.
DBF is the first “Bayes-Faithful” amortized variational inference methodology, constructing test distributions that mirror the inference structure of a SSM with the Markov property.
This structured inference enables analytical computation of test distributions, preventing the accumulation of Monte Carlo sampling errors over time steps. DBF exhibits superior performance over existing methods in scenarios where posterior distributions become highly non-Gaussian, such as in the presence of nonlinear observation operators or significant observation noise.

## Appendix A Settings and additional results for experiments

##### parametrization of the dynamics matrix

We have parametrized the dynamics matrix $A$ following Lusch et al. (): we consider that $h_{dim}/2$ complex eigenvalues $\lambda_{i}(0\leq i<h_{dim}/2)$ characterize $A$. Namely, $A$ is a block-diagonal matrix of $h_{dim}/2$ blocks. Each block consists of $2\times 2$ matrix, whose components are:

$$ $A_{block}=\exp(\rho_{i})\begin{pmatrix}\cos(\omega_{i})&-\sin(\omega_{i})\\ \sin(\omega_{i})&\cos(\omega_{i})\end{pmatrix},$ (11) $$

where $\rho_{i}={\rm Re}[\lambda_{i}]$ and $\omega_{i}={\rm Im}[\lambda_{i}]$. In contrast to Lusch et al. (), we apply the same dynamics matrix at any positions on the latent space. We consider that this representation is sufficiently expressive, as it can express any matrix on a complex number field that is diagonalizable.

One key advantage of DBF is that augmenting the latent dimension only results in a linear increase in computational demand. This scaling is due to the efficient parametrization of the dynamics matrix, where the block-diagonal structure allows operations to scale linearly with the latent dimension. In contrast, methods such as Sequential Monte Carlo (SMC) suffer from exponential increases in computational demand as the latent space grows, assuming that the same density of particles must be maintained to capture posterior distributions. This makes DBF particularly well-suited for high-dimensional systems where traditional methods struggle with computational complexity.

##### Computational resources

We conduct experiments on a cluster of V100 GPUs. Each GPU has memory of 32GB.

##### hyperparameters for training

For all experiments, we have used Adam optimizer with default parameters. Table 3 shows hyperparameters employed in our experiments. Trainings for moving MNIST and double pendulum are conducted with one GPU, while that for Lorenz96 is with eight GPUs.

**Table 3: Hyperparameters for training**
|  | lr | batch size | $h_{dim}$ | $N_{data,train}$ | Epochs | train time per model |
| --- | --- | --- | --- | --- | --- | --- |
| moving MNIST | $10^{-3}$ | 64 | 8 | 480,000 | 2 | 3hr$\times$ 1GPU |
| double pendulum | $10^{-3}$ | 256 | 50 | $1.0\times 10^{7}$ | 1 | 6hr$\times$ 1GPU |
| Lorenz96 | $3\times 10^{-3}$ | 64 | 800 | $2.6\times 10^{7}$ | 1 | 15hr$\times$ 8GPUs |

### A.1 Moving MNIST

##### Dataset:

The dataset consists of a series of 2D images, where each pixel has a dynamic range from 0 to 255. The training set contains $480{,}000$ initial conditions, while the test set consists of ten initial conditions, with both datasets comprising 20 time steps each. The number of training samples and epochs is sufficiently large to ensure that the training converges effectively. A Gaussian noise with a standard deviation of $\sigma=50$ is added to all pixels. The MNIST images of the digits “9” (data point 5740) and “5” (data point 5742) move at constant speeds until they reach the edges, where reflection occurs.

##### Training:

The network weights for $G_{\theta}$ are fixed during the first epoch to facilitate the learning of $f_{\theta}$ and the image tensor for the observation model. Subsequently, $G_{\theta}$ is trained during the second epoch. In total, DBF undergoes training for two epochs.

##### Dynamics model:

Constant velocity model. The exact dynamics matrix we have used is:

$$ $z_{t+1}=Fz_{t}$ (12) $$

$$ $\displaystyle F=\begin{pmatrix}1&0&0&0&dt&0&0&0\\ 0&1&0&0&0&dt&0&0\\ 0&0&1&0&0&0&dt&0\\ 0&0&0&1&0&0&0&dt\\ 0&0&0&0&1&0&0&0\\ 0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&1\\ \end{pmatrix},z_{t}=\begin{pmatrix}x_{1,t}\\ y_{1,t}\\ x_{2,t}\\ y_{2,t}\\ v_{x_{1,t}}\\ v_{y_{1,t}}\\ v_{x_{2,t}}\\ v_{y_{2,t}}\\ \end{pmatrix},$ (13) $$

and true observation model:

$$ $\displaystyle\tilde{x}_{t}$ $\displaystyle=$ $\displaystyle\begin{cases}(x_{t}\ \text{mod 16})&\text{if x//16 is even}\\ 9-(x_{t}\ \text{mod 16})&\text{if x//16 is odd}\\ \end{cases}\ \text{, same for $y$}$ (14) $\displaystyle o_{t}$ $\displaystyle=$ $\displaystyle h(z_{t}),\mathrm{dim}(o_{t})=44\times 44\ \text{, a $28\times 28 $ image is embedded at}(\tilde{x}_{t},\tilde{y}_{t}).$ (15) $$

The formulation above addresses image reflection through the observation operator, resulting in linear dynamics while permitting multiple solutions for each observed figure. This approach presents significant challenges for the EnKF, which assumes a single-peak Gaussian distribution in the assimilating space. To ensure a fair comparison, we revise the dynamics and observation models to allow for a single solution for each figure. This adjustment notably enhances the performance of the EnKF if the image is provided. However, even with this modification, the EnKF fails to accurately estimate the position, velocity, and the embedded image.

##### Network architecture:

$f_{\theta}$: Two-dimension convolutional NNs. Below is the list of layers.

- •
conv1: nn.Conv2d(1, 2, kernel_size=3, stride=2, padding=1)
- •
conv2: nn.Conv2d(2, 4, kernel_size=3, stride=2, padding=1)
- •
conv3: nn.Conv2d(4, 4, kernel_size=3, stride=1, padding=1)
- •
conv4: nn.Conv2d(4, 4, kernel_size=3, stride=1, padding=1)
- •
fc: nn.Linear($11\times 11\times 4$, 8)

The input image, sized $44\times 44$, is sequentially processed by convolutional layers (conv1, conv2, conv3, and conv4). The output is then flattened to serve as the input for the fully connected layer (fc). Ultimately, this process yields eight variables for $f_{\theta}(o_{t})$. The network $G_{\theta}$ follows the same architecture as $f_{\theta}$, but it produces only the diagonal components of $G_{\theta}(o_{t})$ through the NN.

**Table 4: List of hyperparameters for moving MNIST experiment.**
| parameter | value |
| --- | --- |
| $R$ | diag[$e^{6}$] |
| $Q$ | diag[$e^{-4}$] |

##### Example figures:

In Fig. 8, we show example images for observations and all the algorithms in image-informed setting.

Figure: Figure 7: Example figures for two-body moving MNIST experiment. This is the setting explained in the main text. For all algorithms, the two embedded images are not explicitly informed: algorithms need to deal with many unknown parameters in the observation model.
Refer to caption: x15.png

**Table 5: The success rates of different methodologies in the two-body moving MNIST problem. For the model-based approaches, we used the same dynamics and observation models that generated the data. For DBF, the model was initialized with random image tensors and trained solely on the data.**
| Method | Success rate |
| --- | --- |
| DBF | 100% (50/50) |
| EnKF | 58% (29/50) |
| ETKF | 0% (0/50) |
| PF | 0% (0/50) |

Figure: Figure 8: Example figures for two-body moving MNIST experiment. For model-based approaches (EnKF, ETKF, PF), contrary to the experiment reported in the main text, the true images are informed. In data 0, both DBF and EnKF successfully determine and follow the position of the two images. On the other hand, in data 20 and 27, EnKF estimate becomes unstable soon after the two letters overlap. Even in that situation, DBF stably follows the positions of the embedded images.
Refer to caption: x18.png

### A.2 Double pendulum

##### Dataset:

The dataset consists of 2D coordinates representing the positions of two weights. The training set includes $10,240,000$ initial conditions, while the test set contains 10 initial conditions. The number of training samples is sufficiently large to ensure that the training converges. During DVAE training, we observed that some initial conditions resulted in training failure due to instability; however, we maintained the total number of training samples since the training was successful for at least one initial condition. Both datasets comprise 80 time steps. Numerical integration is performed using the solve_ivp function in SciPy, with relative tolerance (rtol) set to $10^{-2}$ and absolute tolerance (atol) set to $10^{-2}$.

A schematic figure explaining the problem setting is presented in panel (a) of Fig. 3 in the main text.

Dynamics model is described in [https://matplotlib.org/stable/gallery/animation/double_pendulum.html](https://matplotlib.org/stable/gallery/animation/double_pendulum.html). The length of the bars is 1 [m], and the positions of the two pendulum weights are observable with Gaussian noise of $\sigma=0.1,0.3$, or $0.5$ [m]. The observation interval is $0.03$ [s]. The task is to predict the positions of the two weights in the successive ten frames.

##### Network architecture:

$f_{\theta}$: A sequence of ten “linear blocks” composed of fully connected layers, layer normalizations, and skip connections. Namely, each linear block has three components:

- •
fc: (input dimension)$\times$ (output dimension) linear layer,
- •
norm: layer normalization,
- •
skip: skip connection.

Taking four observation variables as input, the first linear block expands the dimensionality to 100. The intermediate linear blocks maintain these 100-dimensional variables. The final linear block reduces the 100-dimensional input to a 50-dimensional output, representing 50 latent space variables. The ReLU activation function is applied throughout the network. The structure of $G_{\theta}$ mirrors that of $f_{\theta}$, while $\phi_{\theta}$ serves as the inverse of $f_{\theta}$. The initial eigenvalues are randomly sampled from the range between $e^{0}$ and $e^{0.01}$.

**Table 6: List of hyperparameters for double pendulum experiment.**
| parameter | value |
| --- | --- |
| $R_{init}$ | diag[$1$] |
| $Q$ | diag[$e^{-6}$] |
| initial concentration parameter | $e^{5}$ |

##### Training:

All training variables (network weights for the IOO ($f_{\theta}$, $G_{\theta}$), the emission model operator $\phi$, eigenvalues $\lambda$ for the dynamics matrix $A$, Gaussian noise parameter $\sigma$ for angular velocity $\omega$, and the concentration parameter for Von Mises distribution used for angular coordinate $\theta$) are trained together.

##### Examples:

Here, we show examples for assimilated $\theta$ and $\omega$ in Fig. 9. Also, we give an additional figure for the RMSE of $\theta$ for various methods.

Figure: Figure 9: Results for five example data in test set. Two left columns show evolution of $\theta_{1}$ and $\theta_{2}$ (rad) (, therefore, the values are cyclic with the period of $2\pi\simeq 6.3$, and we corrected for those periodic shifts) and the two right columns show $\omega_{1}$ and $\omega_{2}$ (rad/s).
Refer to caption: x21.png

Figure: Figure 10: Assimilation results for the angle variable $\theta$. All models successfully determine the angle coordinate in spite of the strong nonlinearity in the observation (trigonometric function). Among these, performance of DBF is the best.
Refer to caption: x41.png

### A.3 Lorenz96

##### Dataset:

The dataset consists of physical and observed variables sampled at 40 grid points. The training set includes 25,600,000 initial conditions, while the test set contains 10 initial conditions. The number of training samples is sufficiently large to ensure that the training converges in most cases. The original datasets comprise 80 time steps. Numerical integration is performed using the solve_ivp function in SciPy, with a relative tolerance $\texttt{rtol}=10^{-2}$ and an absolute tolerance of
$\texttt{atol}=10^{-2}$. Gaussian noise with standard deviations of $\sigma=1,3$, or $5$ is added to all measurements.

For KalmanNet, we attempted to train with 25,600,000 and 400,000 initial conditions; however, the process was terminated due to memory limitations. Consequently, we report results using a dataset size of 120,000. For DKF, VRNN, and SRNN, we also tried training with 25,600,000 conditions, but all models encountered a RuntimeError due to instability during the backward computation. To obtain results, we reduced the number of training samples to 512,000. With this adjustment, both SRNN and VRNN successfully completed the training procedure for some initial conditions.

A physical quantity $z_{j}$ is defined at each grid point
$j(1\leq j\leq 40)$. The time evolution of this quantity is described by the following set of differential equations:

$$ $\frac{dz(t)_{j}}{dt}=(z_{j+1}-z_{j-2})z_{j-1}-z_{j}+F,(1\leq j\leq 40)$ (16) $$

In this equation, the driving term $F$ is set to $8$. The first term models the advection of the physical quantity, while the second term represents its diffusion along a fixed latitude. With these parameters, the evolution of the physical quantity exhibits chaotic behavior.

##### Network architecture:

The NN $f_{\theta}$ consists of ten convolutional blocks followed by a fully connected layer. Each convolutional block comprises a 1D convolution, layer normalization, and a skip connection:

- •
conv1d: nn.Conv1d(
$c_{\rm in}$,
$c_{\rm out}$,
kernel_size=5,
padding=2,
padding_mode=“circular”,
)
- •
norm: layer normalization,
- •
skip: skip connection.

The first convolutional block has $c_{\rm in}=1$ and
$c_{\rm out}=20$, expanding the input by a factor of 20 in the channel dimension. The subsequent eight layers maintain 20 channels. Finally, the 20 channels and 40 physical dimensions are flattened into 800-dimensional variables, which are then fed into a fully connected layer of size $800\times 800$. For all layers, the activation function used is ReLU. The function $G_{\theta}$ is structured identically to
$f_{\theta}$, while $\phi_{\theta}$ represents the inverse of $f_{\theta}$.

**Table 7: List of hyperparameters for Lorenz96 experiment.**
| parameter | value |
| --- | --- |
| $R_{init}$ | diag[$1$] |
| $Q$ | diag[$e^{-8}$] |

##### Training:

All training variables, including the network weights for the inverse observation operator $f_{\theta}$ and $G_{\theta}$, the emission model operator $\phi$, the eigenvalues $\lambda$ for the dynamics matrix $A$, and the Gaussian noise parameter $\sigma$, are trained concurrently.

##### Examples:

We show an example figure for assimilation experiment with DBF in Fig. 11.

Figure: Figure 11: An example of assimilation output in the experiment with nonlinear observation operator. The observation is not very informative due to low threshold for saturation in the observation operator ($o_{t,j}=min(z_{t,j}^{4},10)+\epsilon,$, all cells with $z_{t,j}>1.8$ are just observed as $10+\epsilon$). In the first 20 steps, the model output resembles little with the target. However, as the step proceeds, the estimated state begins to capture features of the true state. Even with such a poor observation operator, DBF finds a latent space representation that captures the evolution of the true state.
Refer to caption: x42.png

## Appendix B training stability

We observe that the training of our proposed method is stable compared to RNN-based models. Fig. 12 shows the evolution of the real parts of eigenvalues.
Although we do not impose constraints on the real parts of eigenvalues, the values only marginally exceed one. Therefore, long-time dynamics is stable during training.

Figure: Figure 12: Evolution of histograms for the real parts of 800 complex eigenvalues in Lorenz96 experiment. Initially, eigenvalues are taken as one. As the model learns the dynamics, eigenvalues lower than 1.0 appear. However, the largest eigenvalue $\lambda_{max}$ mostly remains less than 1.02.
Refer to caption: x43.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: (1) In the non-linear case, introducing auxiliary states, which fundamentally change the model and the methodology (e.g. by requiring supervised learning), seems to me a solution that is too ad hoc. Are there other ways in which nonlinearities in state dynamics could be addressed? 

 (2) A similar question concerns the Gaussianity of the IOO. Is there any way in which a slightly more flexible distribution class could be used? 

(2) If I have understood correctly, the numerical part evaluates the different filtering algorithms based on the RMSE with respect to the _true_ states? Since the goal of these algorithms is optimal filtering, i.e. to determine the posterior distributions of the states given the data and not the true, exact states, such a comparison makes it difficult, in my opinion, to claim that one algorithm is better than another. Would it not be better to compare the methods in terms of RMSE with respect to reference values calculated on the basis of a consistent algorithm, e.g. the particle filter with a massively large number (say, half a million) of particles.
soundness: 3
strengths: The idea of learning a Gaussian, NN-parameterized IOO to obtain Gaussian approximations of the filter distributions in the case of linear state dynamics is, to my knowledge, novel. Although I am more doubtful about the relevance of the method for models with nonlinear dynamics, at least in its current form, the simplicity of the proposed method appeals to me. Furthermore, the method's performance, at least in the linear example, is promising, although it is unclear how it generalizes to other models (as the assumed Gaussianity of the IOO introduces a bias that is not easily controlled; see my comment below). The article is fairly clearly written and the methodology is not too difficult to understand. So sum up, the proposed method should have some relevance, at least in the case where the latent-state dynamics is linear.
summary: This paper deals with data assimilation in general state-space models (SSM) and proposes a new method for approximating the flow of marginal state posteriors, also known as filter distributions. In the case of a linear Gaussian latent state process (but a general emission density), by formulating the filter recursion in terms of a so-called inverse observation operator (IOO)—which models the law of the current state given the corresponding observation—and approximating it by a Gaussian distribution parameterized by a neural network, Gaussianity is preserved through the filter recursion, allowing the generation of a sequence of approximate Gaussian filter distributions. The IOO is trained along with the unknown model parameters using variational inference, and in the case of linear dynamics this can be done without supervision. When the state process is nonlinear, the authors introduce another level of auxiliary latent states with linear Gaussian dynamics, generating the original latent state process through an auxiliary emission density. In this case, training cannot be performed without supervision, but requires access to training data consisting of a known history of both latent states and observations. The authors illustrate their method numerically with three examples, one where the goal is to identify images from the MINST dataset that move with linear dynamics and two where the goal is to estimate parameters and states in models with non-linear dynamics describing a double pendulum and a Lorenz system.
weaknesses: As far as I am concerned, the proposed methodology has three major shortcomings, which makes the paper boardeline. 

(i) It seems crucial to use a Gaussian IOO, which is likely to lead to poor approximations for some models, e.g. in cases where the distribution of the current state is multimodal given the corresponding observation. Moreover, the bias implied by the assumed Gaussianity seems very difficult to control.  

(ii) In the case of nonlinear dynamics, the method requires that training is done in a supervised manner, i.e. with access to the observed history of the latent state process, which is only possible in a limited number of contexts. However, in many cases, these states have no direct physical counterparts, but are only a phenomenological tool to capture patterns and features in the observed data (think for example of the unobserved volatility driving the dynamics of asset prices in a stochastic volatility model), preventing the method from being used. This significantly limits the applicability of the methodology.

(iii) To deal with the nonlinear case, the authors introduce an auxiliary state process (the h's), and a Gaussian emission density that generates the original states (the z's) given the auxiliary ones. My concern is that this construction strips the original state process of the Markov property (in the same way as the observations of an SSM are not Markov, even if the underlying states are so), leading to a model with a significantly different dynamics than the original one. Together with point (ii), this leads me to doubt even more the relevance of the method for the nonlinear case.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Can the method be used with approximate, even crude samples of z_t obtained by another method? Would it still add value there?
soundness: 3
strengths: Designing generic-purpose VAE for modeling time series is and will remain an important problem. The authors use nice illustrations to make their point. The use of the filter as opposed to the smoother is noted and that fact that it seems to be (more than) sufficient.
summary: This paper looks at state inference in non-linear state space models, that is, the problem of estimating latent "physical" variables z_t given observations o_t. The posterior of z_t given o_t is often non-Gaussian in models with a practical use, e.g., models of weather dynamics. This poses a challenge, and solutions have taken a path along using Gaussian approximations to the posterior or computationally intensive methods such as Sequential Monte Carlo. A key idea that the paper relies on is the introduction of additional latent variables (state expansion) h_t variables for which Gaussian inference is possible; these variables are later mapped to the original state variables and marginalized. In this sense, the method acts like a dynamic VAE but without the use of Monte Carlo approximations. Essentially, the assumption that the original physical variable dynamics can be linear in a sufficiently high-dimensional space when the mapping to that space is chosen correctly is taken. Two proposed training strategies are given, on the basis of a standard ELBO objective. Some standard examples with a Lorenz model and MNIST are presented to illustrate the method and good performance is shown.
weaknesses: No idea of the computation time for the proposed method is given. Presumably, it should be faster? What if pretraining is taken into account? It is also unclear to me how to choose the dimensionality of the new state space h: is it of a comparable size to z_t or must it be much bigger? As the authors say, they need to know values of at least z_t.

---

confidence: 3
contribution: 3
initial_rating: 3
presentation: 2
questions: 1. Is your approach identifiable? I.e., is there a unique latent (linear dynamics) state trajectory corresponding to a given physical trajectory? If not, what effect does the nonidentifiability have?

2. Can DBF represent multimodal posteriors?

3. line 273: “computations for covariance matrices become challenging in high-dimensional spaces” – can you explain why these are not needed in DBF?

4. One of your key motivations seems to be that particle filters don’t scale to high dimensions, and DBF would scale better to high dimensions. But then you also need to have a significantly larger latent dimension than the actual state space to capture the nonlinear state transitions through the linear embedding: does it still reduce computational cost overall?

5. Also, how do the _state_ dimensions vs the _observation_ dimensions affect the performance of DBF vs particle filters?

6. You mention Frerix et al. (2021), where the Inverse Observation Operator is also learned – could you have included their approach in the empirical comparison?

7. How did you set up & run the baselines – can you link to which implementation you used?

8. Baseline methods could not handle the amount of training data used for DBF. To what extent is better performance due to being able to ingest more data, i.e. how would DBF have performed at the same amount of training data as baselines?

9. line 353: “We define success as …” seems rather arbitrary, and this might bias results (surprisingly related: “Are Emergent Abilities of Large Language Models a Mirage?” by Schaeffer et al., NeurIPS 2023), what are the RMSE results themselves?

10. line 365: can you explain how learning system parameters sets DBF apart from traditional approaches? I am not intimately familiar with the field but surely there are other approaches that learn system parameters.

11. I would have liked to see some sample trajectories vs inferred of the different methods.

12. line 480–485: could you quantify the (non)Gaussianity?

13. I would have liked to see an ablation study; for example, how does the choice of latent dimension affect the accuracy and compute requirements of your approach?

## Clarity questions

14. line 90: What do you mean by “regression methods” here? In general, regression can be probabilistic as well, thus quantifying uncertainty and not only providing point estimates.

15. p.3, Eq. (1): $\mu_t$ shows up in the last term on the right-hand side; is this a mistake? What would be the correct equation?

16. line 145: what do you mean by “virtual prior”, can you explain this in more detail?

17. Eq. (6): what is the KL term of the first element in the sum ($t=1$), where there are no previous observations?

18. Would you call DBF part of the DVAE family or not? In line 224 you say “DVAEs are closely related” (implying no), in line 255 you write “apart from other DVAEs” (implying yes).

19. line 265: “LS4 replaces recurrence with convolutions, forgoing a recursive approach” - does this matter?

20. $\sigma$ in line 339 vs $\sigma_{sys}$ (NB- should be $\sigma_\text{sys}$) in line 346: I don’t understand the distinction, can you explain?

21. line 348: How can the input dimensionality exceed the GPU memory – should this not also depend on sequence length, batch size, etc.?
(NB- it was unclear from the main text where the $44 \times 44$ came from, only explicitly mentioned in supplementary. It is especially confusing as in line 359 you now discuss dimensionality of digits themselves, should this not be irrelevant?)

22. Fig. 2 (b): what are the iterations – training steps of NN ?

23. Fig. 3 (b): is this RMSE of $\omega$ alone (averaged over both velocities, see below)?

24. Table 1 and line 420, “RMSE between physical variables”: RMSE is given for $\theta$ and $\omega$ separately, but there are two of each: is this the average over the two angles / velocities?

25. line 422, excluding failed initial conditions: what is the fraction of initial conditions that failed?

26. line 426: should this be time steps (of dynamics, as opposed to training steps of model)?

27. lines 483, 485: what do you mean by “prior $q(z_t | o_{1:t-1})$” vs “prior $p(z_{t+1} | o_{1:t})$”? Are they the same (except for the shift by 1)? Why $p$ vs $q$?

## Further points for improving clarity

You refer to the “test distribution” $q(z_t)$; in variational inference, this would more commonly be referred to as ‘approximate posterior’ or ‘variational distribution’ – mentioning this equivalence explicitly would make the paper a bit more accessible.

line 62, “assuming Gaussianity”: also mention explicitly, already in the introduction, the assumption of linear dynamics for the latent variables $h$.

line 156: “LGSS” is used but not introduced what it stands for.

line 288: citations for Lorenc and Frerix should be in parentheses.

line 334, “bouncing off frame edges”: at first reading I was quite confused how this clearly non-linear behaviour would match up with underlying _linear_ dynamics; it became clear from reading the supplementary, but would be helpful to mention here already that this nonlinearity is part of the _observation_ model, not the state.

Fig. 4: “for one of data” sounds a bit odd
soundness: 2
strengths: Combining a neural network model to learn nonlinearities in observation and state whilst assuming an underlying _latent_ state with linear dynamics seems original (I am not very familiar with the literature). This fills a gap between hard-to-train DVAEs, Kalman filter approaches that struggle with non-Gaussianity, and particle filters that do not scale to high dimensions. If this holds up in larger-scale experiments, it will be a significant contribution to data assimilation in applications such as weather models.

The paper clearly states the approach’s limitation to require pairs of samples and observations.

In the experiments, a wide range of different methods are included as baselines. The supplementary material provides a high level of detail on experiment settings for DBF itself (though not on the baselines).
summary: The paper proposes a variational inference approach “DBF” to data assimilation, i.e. the task of obtaining the posterior over an SSM state $z_t$ given observations $o_{1:t}$. In contrast to the linear Kalman filter, DBF supports non-linear observation models, as well as non-linear SSMs. In contrast to related dynamical VAEs (DVAEs), DBF does not rely on RNNs and can be trained more easily (less training instabilities).

DBF assumes a latent state $h_t$ that evolves linearly, with a Gaussian variational distribution for the approximate posterior $q(h_t)$, ensuring that the approximate posterior remains Gaussian with analytic updates under new observations. When the true SSM dynamics are linear, $h=z$ and only the (nonlinear) mapping from observation $o_t$ to latent state $h_t$ is learned.

When true SSM dynamics are non-linear, then $h$ is a latent state with higher dimensionality than $z$, but linear dynamics $h_{t+1}=A h_t$; DBF learns these latent dynamics $A$ as well as a latent emission model from the linear state $h_t$ to the actual SSM state $z_t$ from pairs $(z_t, o_t)$ of samples of SSM state $z_t$ and corresponding observation $o_t$.

The form of the variational approximation matches the Markov properties of the underlying SSM, hence termed “Bayes-faithful”.

The paper compares DBF to classical data assimilation algorithms (e.g. ensemble Kalman filter, particle filters) and DVAE variants in three settings (1. linear dynamics on 8-dimensional state, nonlinear observation with 784 dimensions; 2. nonlinear dynamics on 4-dimensional state, nonlinear mixed observation with 2 dimensions; 3. nonlinear dynamics on 40-dimensional state, nonlinear point-wise observation).

Amongst all compared methods, DBF consistently performs best in most settings (the one exception being when the state is directly observed with low noise, thus true posteriors are almost Gaussian and ensemble Kalman methods perform best).
weaknesses: The main weakness of this paper is that the empirical evaluation insufficiently supports the claims made about the proposed DBF.

While DBF’s assumed Gaussianity and linearity in latent space leads to analytic solutions, this can be considered to do closed-form inference in an _approximate_ model. As noted in the paper, given a sufficient number of particles, a particle filter would be able to represent the true posterior. However, the paper does not directly analyse the loss in accuracy due to the approximations made in DBF’s variational inference approach. While it does compare to particle filters, the paper only considers a single, fixed number of particles, and simply states that it is not sufficient. The evaluation that would be needed here is to compare the accuracy–compute trade-off. For particle filter, how do accuracy and resource needs scale with number of particles, and correspondingly, for DBF, how do accuracy and resource needs scale with the number of latent dimensions? Ideally, for a range of dimensions of the actual SSM.

While the paper comes with code for the DBF implementation in the supplementary material, the main text itself is short on detail of the method. For example, the ELBO training objectives are simply given, without explicit derivation.

While ‘Strategy 1’ is introduced to train solely with samples of $z_t$, not requiring pairing with observations, this is not actually included in the empirical evaluation.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. A discussion (or even a comparison) between the normalizing Kalman filter (NKF) and also the model mentioned in the DVAE paper, Kalman Variational Autoencoders (KVAEs), seems necessary, as in these two models, the generative structure is the same as the one shown in Panel (b), Figure 1. 

     - de Bézenac E, Rangapuram SS, Benidis K, et al. *Normalizing Kalman Filters for Multivariate Time Series Analysis*. *NeurIPS*, 2020;33:2995-3007.

2. The loss functions in Equations 6 and 7 still require Monte Carlo sampling if $p(o_t | h_t) $ is non-Gaussian. While DVAEs must sample the latent state trajectory due to the variational distribution is indeed model the smoothing distribution, the benefit of “Bayes-faithful” data assimilation highlighted here is not entirely clear. If Monte Carlo sampling is still required at each step, and if the neural network is shared and amortized over each time step $t$ (which is the same as in DVAEs), what is the specific advantage of this “Bayes-faithful” property? Are there empirical results demonstrating the practical benefit of this property?

3. Follow the question 2, the paper claims to present the first VAE-based model for time-series data that maintains a posterior structure faithful to the Markov property in SSMs. However, recent work on Auto-EnKF models also employs filtering distributions to construct the model evidence or ELBO. For example, Equation (32) in Lin et al. (2024), shares a similar loss function structure to Equation (6) in this paper. Including a discussion or comparison of these models would enhance the comprehensiveness of the proposed method's positioning. 

     - Chen Y, Sanz-Alonso D, Willett R. "Autodifferentiable Ensemble Kalman Filters," *SIAM Journal on Mathematics of Data Science*, 2022; 4(2):801-33.

     - Chen Y, Sanz-Alonso D, Willett R. "Reduced-order Autodifferentiable Ensemble Kalman Filters," *Inverse Problems*, 2023; 39(12):124001.

     - Lin Z, Sun Y, Yin F, Thiéry AH. "Ensemble Kalman Filtering Meets Gaussian Process SSM for Non-Mean-Field and Online Inference," *IEEE Transactions on Signal Processing*, 2024 Aug 22.



4. Other questions:

 - Adding derivations for Equations (1, 2, 4, and 5) in the Appendix would improve clarity.

 - Could the authors discuss potential applications for the proposed method? 

 -  In Section 3.1, the authors mention that “DVAEs are not comparable as they need to undergo supervised training,” which may need clarification. If I’m not mistaken, DVAEs typically require only observation data and do not rely on model knowledge or latent state data. Could you clarify why DVAE training is interpreted as supervised in this context?

 - Given the low latent dimension and the known model structure, it seems model-based methods (e.g., EnKF, PF) should perform reasonably well. Could you explain why DBF and other data-driven methods outperform these model-based approaches in Table 1?
soundness: 2
strengths: The paper is well-structured, making it easy to follow the ideas presented. The proposed method can be considered to be novel and could be particularly beneficial for certain generic high-dimensional dynamical systems, especially those that can offer underlying latent states for training.
summary: The paper presents Deep Bayesian Filtering (DBF), a method for data assimilation within nonlinear state space models (SSMs). It addresses the prevalent challenge of non-Gaussian posterior distributions encountered in many scientific applications. By incorporating new latent variables, DBF ensures Gaussianity in the latent space, facilitating linear transitions and inverse Gaussian observation operators. This approach claims that one of the advantages is to avoid errors associated with successive Monte Carlo sampling, as seen in DVAEs.
weaknesses: The primary limitation, in my view, is that the model necessitates extensive training data, including the underlying system's latent states (nonlinear dynamics), to effectively learn the IOO. This requirement could render DBF less practical for real-time applications or in domains where data availability is limited, making model-based methods potentially more viable. Additional issues and questions are discussed in the Questions part.

Meta-review: The paper 'Deep Bayesian Filter for Bayes-Faithful Data Assimilation' was reviewed by 4 reviewers who gave it an average score of 5.5 (final scores 5+5+6+6). The reviewers found the paper relevant, appreciated the experiments, and easy to follow. On the other hand, the reviewers had concerns related to the practical aspects of the approach, insufficient evaluations, and lack of details. Some of these issues were addressed during the discussion phase. However, as the paper now stands it remains borderline, and the right decision appears to be to reject the paper in its current form to allow the authors to improve their work.

Decision: Reject