Auto-differentiable data assimilation: Co-learning of states, dynamics, and filtering algorithms
Authors: Melissa Adrian, Daniel Sanz-Alonso, Rebecca Willett
Organizations: School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907 · Department of Statistics, The University of Chicago, Chicago, IL 60637 · Department of Computer Science, The University of Chicago, Chicago, IL 60637 · NSF-Simons National Institute for Theory and Mathematics in Biology, Chicago, IL 60611
Data assimilation algorithms estimate the state of a dynamical system from partial observations, where the successful performance of these algorithms hinges on costly parameter tuning and on employing an accurate model for the dynamics. This paper introduces a framework for jointly learning the state, dynamics, and parameters of filtering algorithms in data assimilation through a process we refer to as auto-differentiable filtering. The framework leverages a theoretically motivated loss function that enables learning from partial, noisy observations via gradient-based optimization using auto-differentiation. We further demonstrate how several well-known data assimilation methods can be learned or tuned within this framework. To underscore the versatility of auto-differentiable filtering, we perform experiments on dynamical systems spanning multiple scientific domains, such as the Clohessy-Wiltshire equations from aerospace engineering, the Lorenz-96 system from atmospheric science, and the generalized Lotka-Volterra equations from systems biology. Finally, we provide guidelines for practitioners to customize our framework according to their observation model, accuracy requirements, and computational budget.
Figures & tables
Alg.
Mθ(xt−1n)
Kϕ(xt−11:N)
Y(⋅;Rt)
Ens. Size
3DVar- C
Fθ1(xt−1)
CϕHt⊤(HtCϕHt⊤+Rt)−1
yt
1
Cϕ∈Rdx×dx
EnKF
Fθ1(xt−1n)+Q~θ2ξtn
C^tθ,ϕHt⊤(HtC^tθ,ϕHt⊤+Rt)−1
yt+γtn
NEnKF
ξtn∼N(0,Idx)
( C^tθ,ϕ defined in ( 8 ))
γtn∼N(0,Rt)
Ens3DVar
Fθ1(xt−1n)+Q~θ2ξtn
C~tθ,ϕHt⊤(HtC~tθ,ϕHt⊤+Rt)−1
yt+γtn
NEns3DVar
ξtn∼N(0,Idx)
( C~tθ,ϕ defined in ( 9 ))
γtn∼N(0,Rt)
Table 1: Table of examples of named special cases of the general filtering algorithm defined in ( 4 ). We note that this list is not exhaustive.
Forecast cost
Analysis cost
3DVar- C
cF
O(dy3+dx2dy+dy2dx)
EnKF
NcF
O(dy3+dy2N+dyN2+dxN2)
Ens3DVar
NcF
O(dy3+dx2dy+dy2dx)
Table 2: Table of computational costs for the forecast and analysis in one filtering step for straight-forward implementations of the algorithms 3DVar- C , EnKF, and Ens3DVar, where we assume that N≪dy for Ens3DVar.
Figure 1: Clohessy-Wiltshire parameter estimation (Section 4.1 ). Parameter estimation error of the learned forecast parameter θ^1 compared to the ground truth parameter θ1∗=0.0013 across choices of σ02 and variance on the diagonals of R . The shaded regions correspond to the 0.25 and 0.75 quantiles of errors across 10 independent simulations. We additionally plot the median error for the initialized parameters {θ1,j(0)}j=1J for j=1,…,J=10 in the dashed black line (c) representing the errors in the absence of any learning, and the gray region corresponds to the 0.25 and 0.75 quantiles of test forecast RMSEs across these 10 initialized models.
Figure 2: Clohessy-Wiltshire filtering (Section 4.1 ). Visualization of observations, true trajectories, and estimated positions for x1 and x3 , and velocities dtdx1 and dtdx3 across each of the four filtering methods using the learned forecast parameters θ^ and learned filtering parameters ϕ^ , evaluated on test data. In this example, the observation noise is fixed at R=0.1Idy×dy , and the variance of the initialization of θ1(0) is fixed at σ02=0.03 .
Figure 3: Clohessy-Wiltshire data log-likelihood (Section 4.1 ). Plots of the forecast log-likelihood evaluated at observations from ( 10 ) for assimilating test observations into filters parameterized by the learned forecast parameters θ^ and ϕ^ learned from AD-3DVar- C , AD-EnKF ( N=25 ), and AD-Ens3DVar ( N=25 ). The algorithm used for filtering corresponds to the algorithm used in the learning task (i.e., EnKF is used as the filter with the learned parameters from AD-EnKF). The log-likelihood from the optimal solution to the filtering task, the Kalman filter (detailed in Algorithm 3 ) using the true forecast parameters θ^ , is shown in each plot as a comparison.
Figure 4: Lorenz-96 learned forecast performance (Section 4.2 ). Test performance of the learned Lorenz-96 model in ( 22 ) learned from various methods: AD-3DVar- C , AD-3DVar- K , AD-EnKF with a covariance tapering radius of 5 and N=25 , and AD-Ens3DVar with a covariance tapering radius of 5 and N=25 . The shaded regions correspond to the 0.2 and 0.8 quantiles of test forecast performances across 10 independent simulations. In the plot varying dy/dx , we additionally plot the median forecast RMSE for the initialized models Fβj for j=1,…,J=10 in the dashed black line (c), and the gray region corresponds to the 0.2 and 0.8 quantiles of test forecast RMSEs across these 10 initialized models.
Figure 5: Lorenz-96 filtering (Section 4.2 ). Visualization of filtering on test observation data using the learned θ^ and ϕ^ values across the methods AD-3DVar- C , AD-3DVar- K , AD-EnKF with N=25 ensemble members and a tapering radius of 5, and AD-Ens3DVar with N=25 ensemble members and a tapering radius of 5 compared to the noisy observations and ground truth Lorenz-96 simulations. In this plot, the settings of the filtering problem are R=Idy×dx , dy/dx=0.6 , σ02=1 , and dx=200 . For ease of visualization, only the first 50 out of 200 state components are shown in the plot. Filtering is done with the base-filtering algorithm used during training (i.e., the results presented for AD-EnKF use the learned θ^ and ϕ^ in an EnKF algorithm).
Figure 6: Example interaction matrix A for 16 interacting species in the generalized Lotka-Volterra equations. The values ±a1:10 denote the shared parameter values within each block.
Figure 7: Generalized Lotka-Volterra parameter estimation (Section 4.3 ). Parameter estimation mean absolute errors (MAEs) of the learned generalized Lotka-Volterra model in ( 26 ) learned from the methods AD-3DVar- C , AD-EnKF with N=dx , and AD-Ens3DVar with N=dx . The shaded regions correspond to the 0.2 and 0.8 quantiles of test forecast performances across 10 independent simulations. In each plot, we plot the median MAEs for the initialized models Fθ1j for j=1,…,J=10 in the dashed black line, and the gray region corresponds to the 0.2 and 0.8 quantiles of test forecast RMSEs across these 10 initialized models. The top row of plots correspond to the error in the estimated interaction matrix A^ component of θ^1 , and the second row of plots correspond to the error in the estimated population change vector r^ component of θ^1 .
Figure 8: Generalized Lotka-Volterra filtering (Section 4.3 ). Visualization of filtering on test observational data using the learned θ^ and ϕ^ values across the methods AD-3DVar- C , AD-EnKF with N=500 ensemble members, and AD-Ens3DVar with N=500 ensemble members where dx=500,Rt=0.05Idyt×dx , σ02=0.1, and dyt/dx=0.5 . For ease of visualization, only 10 randomly selected species abundance trajectories out of the 500 total species are visualized. In each plot in the first row, the true trajectories are given by the dashed black lines, the sparse, noisy observation are plotted as a scatterplot, and the analyses are plotted as solid colored lines. The second row plots the errors in the analyses compared to the ground truth trajectories.
Time-
High-
Algorithm
Ht≡H
varied Ht
dimensional dx
Memory
Learning Fθ1
AD-EnKF (Modified)
Reasonably
✓
May require
High; scales
Accurate
( Chen et al.,, 2022 )
recovers
non-trivial
with N
all states
modifications
AD-Ens3DVar
Reasonably
✓
May require
High; scales
Accurate
(Ours)
recovers
modifications
with N
Table 3: Summary table of particular examples of auto-filtering algorithms in terms of their ability to learn accurate filtering parameters when Ht≡H when dy<dx , ability to filter with a time-varying Ht , relative scalability to high-dimensional systems, relative memory footprint, and relative accuracy in improving the forecast model Fθ1 .
Appendix figures & tables6 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 9: Experiment based on the Lorenz-96 system described in Section 4.2 . (Left) Median test forecast RMSE for learning a low rank update to the background covariance such that Cϕ=C0+BϕBϕ⊤ , where Bϕ∈Rdx×p and p<dx . In this experiment, dx=200,dy/dx∈[0.3,0.5,0.8] , σ02=1 , and observations at each time t are randomly sampled across state dimensions (i.e., Ht varies with time). Each training task was replicated for J=10 independent simulations. We additionally plot the median forecast RMSE for the initialized models Fβ(j) for j=1,…,M in the dashed black line, and the gray region corresponds to the 0.2 and 0.8 quantiles of test forecast RMSEs across these 10 initialized models. (Right) Plot of millions of floating point operations per second (MFLOPS) across values of p for processing one subsequence (both the forward and backward pass) for the setting where dx=200 , p∈[50,75,100,150,175,200],L=20, and dy/dx∈[0.3,0.5,0.8] .
Figure 10: Heatmap of the EnKF filter RMSE values when filtering with the true Fθ1∗ on unseen test data across different choices of tapering radius and covariance inflation for dx=40 , dy/dx=1 , and σ02=1 . The red trajectory corresponds to the convergence of the analysis parameter vector ϕ∈R2 for the covariance tapering and inflation for AD-EnKF (Modified), which was initialized at tapering = 5 and inflation = 1.2. The red star corresponds to location of the converged value for the (tapering, inflation) pair.
Figure 11: Varying the ensemble size for AD-EnKF (Modified) and AD-Ens3DVar for the Clohessy-Wiltshire (Section 4.1 ), Lorenz-96 (Section 4.2 ), and generalized Lotka-Volterra (Section 4.3 ) examples. The shaded regions correspond to the 0.2 and 0.8 quantiles of forecast performance or forecast parameter estimation error across 10 independent simulations.
Figure 12: Visualization of the test filter RMSEs across examples and methods. The lines in each plot correspond to the test filter RMSE across parameter choices for either the covariance mixing parameter ϕ1 in Ens3DVar or the covariance inflation factor ϕ in EnKF. The top row corresponds to filtering results using EnKF, the bottom row corresponds to filtering results using Ens3DVar, and each column corresponds to results from one of our experimental examples. The star in each plot corresponds to the value chosen by the learning algorithm. The EnKF evaluations use the learned forecast parameters θ^ in the filter, and the Ens3DVar evaluations use both the learned forecast parameters θ^ and the learned 3DVar background covariance parameters ϕ2 in the filter.
Figure 13: Clohessy-Wiltshire filtering (Section 4.1 ). Visualization of observations, true trajectories, and estimated positions for y and z , and velocities dtdy and dtdz across each of the four filtering methods using the learned forecast parameters θ^ and learned filtering parameters ϕ^ , evaluated on test data. In this example, the observation noise is fixed at R=0.1Idy×dy , and the variance of the initialization of θ1(0) is fixed at σ02=0.03 .
Figure 14: Visualization of each component of the forecast log-likelihood evaluated at observations, given in ( 14 ) for filtering with the learned θ^ and ϕ^ from AD-3DVar- C , AD-EnKF, and AD-Ens3DVar compared to the optimal Kalman filter. (Left) Component of the log-likelihood that computes −21logdetSt for each time t . (Right three) Residual component of the loss that computes −21∥yt−Htm^tθ∥St2 for each time t .
Data assimilation (DA) is an essential tool for prediction and understanding in the geosciences. DA combines simulation programs representing scientific knowledge with observations that constrain system dynamics, resulting in analyses and forecasts that incorporate both knowledge and data. DA tasks can be addressed with a diverse toolset, including variational, ensemble and learning-based methods. In particular, many recent works have proposed using automatic differentiation tools for variational, learning-based or hybrid methods. However, comprehensive comparisons across algorithms and dynamical systems remain challenging, due to the incompatibility of simulation and assimilation codes, inflexible handling of spatial and temporal discretizations, specialization of DA methods to specific simulations, and limited support for automatic differentiation and parallel computation in simulations. To address this challenge, we introduce Automatic Differentiation for Data Assimilation (ADDA), a software framework for defining and working with system states, simulations, observation schemes and DA methods. ADDA provides a powerful and flexible set of base classes for representing dynamical systems and observation operators, with support for collocated and staggered grids, unstructured meshes, Lagrangian state variables and irregular or continuous-time observations. Parallel processing and differentiability are first-class features, with support for batch axes and automatic differentiation throughout. ADDA is implemented in PyTorch library, but supports DA for JAX-based computation of dynamics and their gradients. To demonstrate its features, we further provide differentiable, ADDA-compatible implementations of 10 dynamical systems of various dimensionalities and scales, from which we design multiple illustrative DA examples. All of our code is publicly available at https://github.com/m-dml/ADDA.
Anthony Frion, Vien Minh Nguyen-Thanh, Ali Can Bekar +3
Model-Driven Machine Learning Institute of Coastal Systems - Analysis and Modeling Helmholtz-Zentrum Hereon Germany
Data assimilation (DA) estimates the state of an evolving dynamical system from noisy, partial observations, and is widely used in scientific simulation as well as weather and climate science. In practice, filtering methods rely on frame-to-frame transition models. However, these models are fragile when observations are non-Markovian (when they form only a partial slice of a higher-dimensional latent state as in real-world weather data): they tend to accumulate errors over long horizons. At the same time, learned DA methods typically commit to a single regime, either filtering (nowcasting, real-time forecasting) or smoothing (retrospective reanalysis), which splits what should be a shared prior across application-specific pipelines. To address both issues, we introduce ForcingDAS, a unified and robust DA framework. Built on Diffusion Forcing with an independent noise level assigned to each frame, ForcingDAS learns a joint-trajectory prior instead of frame-to-frame transitions. This allows it to capture long-horizon temporal dependencies and reduce error accumulation. In addition, the same trained model spans the full filtering to smoothing spectrum at inference time. Specifically, nowcasting, fixed-lag smoothing, and batch reanalysis are selected through the inference schedule alone, without retraining. We evaluate ForcingDAS on 2D Navier-Stokes vorticity, precipitation nowcasting, and global atmospheric state estimation. Across all settings, a single model is competitive with or outperforms both learned and classical baselines that are specialized for individual regimes, with the largest gains observed on real-world weather benchmarks.
Yixuan Jia, Siyi Chen, Yida Pan +9
University of Michigan, Ann Arbor · Michigan State University · Los Alamos National Laboratory
Bayesian filtering and smoothing are central to data assimilation in nonlinear dynamical systems. Recent advances in deep generative models provide flexible approximations of the associated non-Gaussian posterior distributions. However, several existing approaches require the score to be re-estimated or a transport map to be constructed at each assimilation step. We propose an amortized framework for filtering and smoothing that reuses trained conditional models across observation sequences and assimilation times. The framework jointly learns a shared recurrent summary network and two conditional normalizing flows from simulated state and observation trajectories. The recurrent network represents each observation history by a fixed-dimensional summary that conditions the filtering approximation and, together with the next state, the backward transition approximation. An information-theoretic analysis shows that, under the Markov assumption, a summary sufficient for filtering is also sufficient for the backward kernel. Combining the terminal filtering approximation with the backward kernels yields an approximate joint smoothing distribution. Numerical experiments on advection-diffusion, Burgers, and Lorenz systems demonstrate the accuracy of the proposed filtering and smoothing approximations and characterize the evolution of their errors beyond the training horizon.
Tiangang Cui, Xiaodong Feng, Chenlong Pei +2
School of Mathematics and Statistics, The University of Sydney, Camperdown, NSW 2006, Australia · Institute of Advanced Study, Beijing Normal-Hong Kong Baptist University, Zhuhai 519087, China · SKLMS & Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and Systems Science, Chinese Academy of Sciences +1