Source-linked AI summary
Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization
Marc Bocquet, Julien Brajard, Alberto Carrassi, Laurent Bertino
TL;DR
Learning chaotic dynamics from partial and noisy observations requires methods suited to long time series and unstable systems. This paper unifies data assimilation and machine learning in a Bayesian framework, using expectation-maximization and coordinate descent to estimate the model, trajectory, and model-error statistics. The methods are tested on two low-order chaotic models with distinct identifiability.
Problem
Learning high-dimensional chaotic dynamics is challenged by partial and noisy observations, long time series, and unstable dynamics.
Source: §1.2 · p. 2Method
A Bayesian framework combines data assimilation and machine learning, using expectation-maximization and coordinate descent to estimate the model, trajectory, and model-error statistics.
Source: §1.4 · p. 3, §4.2 · p. 9, §7 · p. 23Results
The full and approximate expectation-maximization schemes were successfully tested on two low-order chaotic models with distinct identifiability.
Source: §7 · p. 23Takeaways & Limitations
The framework unifies approaches for learning dynamics from partial and noisy data.
Source: §1.4 · p. 3The slow-variable-only L05III reference model is not identifiable within the surrogate-model space, precluding exact recovery of its governing equations.
Source: §6.1.2 · p. 16
Abstract
from arXiv · showhide
The reconstruction from observations of high-dimensional chaotic dynamics such as geophysical flows is hampered by (i) the partial and noisy observations that can realistically be obtained, (ii) the need to learn from long time series of data, and (iii) the unstable nature of the dynamics. To achieve such inference from the observations over long time series, it has been suggested to combine data assimilation and machine learning in several ways. We show how to unify these approaches from a Bayesian perspective using expectation-maximization and coordinate descents. In doing so, the model, the state trajectory and model error statistics are estimated all together. Implementations and approximations of these methods are discussed. Finally, we numerically and successfully test the approach on two relevant low-order chaotic models with distinct identifiability.
1. Introduction.
The introduction frames learning chaotic, potentially high-dimensional dynamics from partial and noisy observations as a problem for both data assimilation and machine learning. It proposes a Bayesian framework that includes model error statistics and uses expectation-maximization and coordinate descent to estimate the model, state trajectory, and error statistics together.
1.1. Data assimilation.
- Data assimilation infers physical-system states from observations and theoretical knowledge, including governing equations, and supports forecasting unstable geophysical flows with imperfect models.The introduction identifies the atmosphere, ocean, and climate as examples.Source: §1.1 · pp. 1–2
1.2. The model as a control variable.
- The paper identifies partial and noisy observations, data requirements, and short-term prediction versus asymptotic behavior as questions for learning dynamics.The questions are posed for inference of dynamics, not answered in the introduction.Source: §1.2 · p. 2
1.3. Surrogate model representations from machine learning.
- Machine-learning representations of dynamics can be incorporated into a data-assimilation formalism to address partial and noisy observations, a setting motivating the use of both approaches.The introduction also notes that learned surrogate models can represent dynamics through a flow rate or resolvent.Source: §1.4 · p. 3, §1.3 · p. 3, §1.2 · p. 3
1.4. Problem and objectives.
- The framework jointly estimates a surrogate model, state trajectory, and model-error statistics while combining data assimilation with machine-learning representations for partial and noisy data.The authors describe this as a generalization of the data-assimilation framework.Source: §1.4 · p. 3, §1.2 · p. 3
- The study assumes autonomous governing equations, possibly chaotic and high-dimensional systems, and partial, noisy observations; it treats observation-error statistics as known.Model error is injected between consecutive observation updates in the discretized data-assimilation system.Source: §1.4 · pp. 3–4
1.5. Outline.
- Arbitrary dynamics may require too many degrees of freedom at high dimension, while geophysical-flow properties and symmetries may yield sparse, practical representations.The paper outlines model representations, Bayesian inference, optimization, and numerical tests.Source: §2 · p. 4, §1.5 · p. 4
2. Model representation.
The paper represents surrogate dynamics through a flow rate, using locality and optional spatial homogeneity to make high-dimensional models tractable. It also describes a neural-network implementation whose integrated and composed flow-rate maps form the resolvent used to propagate states between observations.
2.1. Symmetries and constraints.
- Locality restricts each variable’s flow rate to a spatial stencil, while homogeneity can further reduce the number of model parameters.The authors describe homogeneity as less fundamental because many geophysical models are inhomogeneous.Source: §2.1 · pp. 4–5
2.2. Sleek representation of the dynamics.
- The sleek flow-rate representation uses constant, linear, and bilinear regressors; imposing locality reduces parameter scaling from quadratic to linear in state dimension.With homogeneity, the coefficient matrix becomes a vector; its size depends on stencil size and spatial dimension.Source: §2.2 · p. 5
2.3. Neural network representation of the dynamics.
- A neural network can represent the flow rate by combining learned linear and bilinear features, with biases supplying constant terms.The described dense architecture can use linear activation to mimic the sleek representation, though other choices are possible.Source: §2.3 · p. 5
- Convolutional layers encode locality and homogeneity, whereas locally connected layers encode locality alone; the paper uses convolutions to support high-dimensional systems.The authors contrast their flow-rate network with prior work that learned the resolvent directly.Source: §2.3 · p. 5
2.4. Time integration and resolvent.
- The flow-rate network is integrated into a one-step resolvent and composed Nc times to propagate the surrogate state between observations.The interval between observations is assumed to be a multiple of the integration step; gradients are generated by a machine-learning library.Source: §2.4 · p. 6, Figure 1 · p. 6, §3 · p. 6
3. Bayesian analysis.
The Bayesian analysis formulates inference over the model, state trajectory, and model-error statistics under Gaussian error assumptions. It distinguishes joint state-and-model estimation from marginal model estimation, and relates its cost function to data assimilation and, under restrictive assumptions, the classical machine-learning loss.
3.1. Prior error statistics.
- The analysis assumes Gaussian observation and model errors, with known observation-error covariance and model-error covariance that may be unknown, white in time, and independent of observation errors.Other error distributions are noted as possible.Source: §3.1 · p. 7
3.2. Joint estimation of the model, its error statistics and the state trajectory.
- Joint estimation seeks a maximum a posteriori solution for the model, trajectory, and model-error statistics conditioned on observations.The posterior hierarchy combines observation likelihood, trajectory prior, and a joint prior on the model and its error statistics.Source: §3.2 · pp. 7–8
- The resulting Bayesian formulation treats partial and noisy observations rigorously and reduces to the classical machine-learning loss when observations are complete and noiseless and model error is known.The cost function resembles weak-constraint 4D-Var.Source: §3.2 · p. 8
3.3. Joint estimation of the model and the error statistics.
- The marginal approach estimates the model and error statistics without selecting a state realization; the paper’s main optimization method targets this full marginal solution rather than a Laplace approximation.The authors also unify several optimization approaches as approximate solutions to the same Bayesian problem when model-error statistics are known.Source: §3.3 · p. 8
4. Optimization with known model error statistics.
With model-error statistics known, the paper compares joint variational optimization with coordinate-descent strategies that separate model and trajectory updates. It discusses variational and ensemble data-assimilation options for estimating the trajectory, including their practical trade-offs for long windows and model error.
4. Optimization with known model error statistics.
- Joint L-BFGS minimization is efficient for shorter training windows, but its growing inverse-Hessian approximation makes it impractical as the window length increases.The cited approach does not exploit the possible independence of distant trajectory segments under an ergodicity assumption.Source: §4.1 · p. 9
4.1. Joint variational solution.
- Coordinate descent alternates model and trajectory minimization, a strategy the authors expect to become more efficient as the training window grows.The approach distinguishes the variables’ different roles.Source: §4.2 · p. 9
4.2. Coordinate descent.
- For the model update, L-BFGS using neural-network-library gradients converged faster and systematically achieved better accuracy than stochastic gradient descent in the authors’ tests.The authors use this L-BFGS scheme for all numerical experiments.Source: §4.3 · p. 9, §4.2 · p. 9
4.3. Minimization over A.
- With the surrogate model fixed, trajectory estimation is a weak-constraint data-assimilation problem that can use weak-constraint 4D-Var or an ensemble Kalman smoother.The coordinate-descent setup assumes known model- and observation-error statistics.Source: §4.4 · p. 9, §4.4.2 · p. 10
4.4.1. With variational subproblems.
- Parallel L-BFGS over possibly overlapping time subwindows enables long-window training, but does not handle uncertain model-error statistics without modification.The authors do not emphasize this approach because ensemble methods can account for such uncertainty.Source: §4.4.1 · pp. 9–10
4.4.2. With ensemble data assimilation.
- An ensemble Kalman filter or smoother can estimate the trajectory while accounting for sampling and model error; smoothing can improve trajectory accuracy over filtering.In the cited earlier approach, only the posterior ensemble mean is used, with model error represented by adaptive inflation and additive stochastic noise.Source: §4.4.2 · p. 10
5. Full joint optimization with expectation-maximization.
The section formulates joint inference of model parameters and model-error statistics as a Bayesian marginal optimization problem, then develops EM and coordinate-descent schemes to solve it. Its approximate scheme reduces ensemble storage but fixes the Q update within each EM expectation step.
5.1. Principle of expectation-maximization.
- EM alternates an expectation over trajectories conditioned on observations and current parameters with a maximization that updates those parameters toward a local maximum.Here the parameters are A and Q, the observations are y0:K, and the marginalization is over trajectories x0:K.Source: §5.1 · p. 11, §5.1.2 · p. 11, §5.2 · p. 11
5.2. Full scheme.
- The full scheme uses an EnKS with SQRT-CORE to produce a smoothed trajectory ensemble, then minimizes the EM objective over A and Q.Joint L-BFGS optimization requires storing the full ensemble of trajectories over the training window.Source: §5.2.1 · pp. 11–12, §5.2.3 · p. 12
- A coordinate descent can alternate an A minimization with an analytic Q update, with a single A minimization sufficient when Q=qI.In that case Q only rescales the loss, so it does not change the A minimizer.Source: §5.2.3 · p. 13
5.3. Approximate scheme.
- The approximate scheme accumulates objective terms online and stores only the smoothed ensemble mean, avoiding storage of the full ensemble.It computes Q once during the expectation step and does not iterate the A and Q minimizations; the authors expect this approximation to be mild when Q=qI.Source: §5.3 · pp. 13–14
5.4. Hyperpriors.
- The Q hyperprior can matter more when training windows are short, while the scalar Jeffreys-prior correction is minor for large systems or long windows.The section also discusses regularizing A, but avoids it here because the windows are mostly very long and A is constrained by locality and/or homogeneity.Source: §5.4.1 · p. 15, §5.4.2 · p. 15
6. Numerical illustrations.
The numerical illustrations test the learning methods on two low-order chaotic reference models. The L05III experiment observes and reconstructs only slow variables, while the models may not be exactly identifiable by the neural-network surrogate.
6. Numerical illustrations.
- The experiments test the methods on the Lorenz-96 and two-scale Lorenz-05III models, using an RK4-based neural-network surrogate for both.The models are spatially distributed one-dimensional systems defined over circles.Source: §6.1 · p. 15, §6.2 · p. 16
6.1.1. Lorenz-96.
- Lorenz-96 uses 40 state variables with forcing F=8, and its reference trajectory is integrated with RK4 at time step 0.05.The classical data-assimilation setup observes every variable every 0.05 time units.Source: §6.1.1 · pp. 15–16
6.1.2. Lorenz-05III.
- Lorenz-05III has 36 slow and 360 fast variables, but only the slow variables are observed and represented in the surrogate.Because only its slow scale is reconstructed, L05III is not identifiable in the surrogate-model space; the authors also note that L96 may not be identifiable by the neural-network representation.Source: §6.1.2 · p. 16
6.2. Surrogate model.
The paper evaluates surrogate models with a compact neural-network architecture and assesses both short-term forecast behavior and asymptotic dynamics. Its metrics include forecast skill, Lyapunov spectra, and power spectra.
6.2. Surrogate model.
- The minimal neural network uses a width-5 convolution with six filters and a point-wise output layer, totaling 49 weights under locality and homogeneity assumptions.Larger widths or more filters did not noticeably improve results, and tested ReLU or tanh activations had no significant benefit or could degrade accuracy.Source: §6.2 · pp. 16–17
6.3. Evaluation metrics.
- Evaluation combines normalized forecast error with Lyapunov spectra and power spectral density to assess short-term and asymptotic behavior.Forecast skill is averaged over 5 × 10^3 runs; an NRMSE above 1 denotes a useless surrogate forecast.Source: §6.3 · p. 17
6.4. First results.
The nominal configuration gives reliable asymptotic properties for L96, while L05III shows greater variability and weaker short-term predictability. The two models differ in identifiability, and the screening results also show that forecast skill alone does not ensure optimal asymptotic indicators.
6.4.1. Preliminary screening of hyperparameters.
- The nominal setting L = 4, K = 5000 was selected for satisfactory forecast skill, but that criterion does not guarantee optimal asymptotic indicators.The screening varied EnKS lag and training-window length; the authors ran 10 experiments with new trajectories and observations to assess variability.Source: §6.4.1 · pp. 17–18
- Table 1 reports scalar indicators for nominal experiments based on L96 and L05III, with indicator statistics computed over 10 samples; it also recalls the key hyperparameters.Source: Table 1 · p. 18
- Using full Q gives a slight improvement over Q = qIx, while the inferred full covariance matrices are systematically close to diagonal.The added computational cost was mild in these experiments, so the authors used full Q in subsequent reported experiments.Source: §6.4.1 · p. 18
6.4.2. Results for the L96 reference model.
- For L96, surrogate Lyapunov exponents and power spectra closely match the reference, except for high-frequency PSD saturation from the inferred stochastic noise.Forecast skill varies modestly across samples; without the noise, high-frequency PSD agreement would be closer but misleading.Source: §6.4.2 · p. 18, Figure 2 · p. 19
6.4.3. Results for the L05III reference model.
- For L05III, all three indicators vary more than for L96, and forecast skill is more sample-sensitive, with steeper short-term error growth attributed to weaker identifiability.Forecast skill relaxes to climatology more slowly; the text relates this to heterogeneous predictability as fast-variable bursts affect the slow compartment.Source: §6.4.3 · pp. 19–20
- L05III surrogate Lyapunov exponents track the reference but are systematically smaller, possibly because the surrogate learns only the average impact of fast variables.The PSDs are also close except at high frequencies, where inferred additive noise yields saturation and acts as a stochastic subgrid parameterization.Source: §6.4.3 · p. 20
6.5. Comparison of full and approximate scheme.
The approximate scheme is less computationally demanding than the full scheme, with similar L96 results and a distinct forecast-skill advantage for L05III. Its MAP-based weight estimate may favor short-term deterministic forecast skill at the expense of uncertainty representation.
6.5. Comparison of full and approximate scheme.
- For L96, full and approximate schemes give very close results: approximate forecast skill is slightly better, while Lyapunov spectra and power spectra are almost identical.The full scheme requires storage at least Ne times larger and its cost function is Ne times more demanding than the approximate scheme's.Source: §6.5 · p. 20
- For L05III, the approximate scheme's forecast skill is distinctively better than the full scheme's, while Lyapunov spectra and power spectra remain very close.The authors attribute this advantage to the approximate scheme seeking A as the MAP rather than the average, potentially improving short-term deterministic skill while risking overfitting and less reliable uncertainty representation.Source: §6.5 · p. 20
- The proposed schemes are presented as able to handle long training windows.The supplied passage states this advantage but does not include the accompanying results.Source: §6.6 · p. 20
6.6. Sensitivity on the length of the training window.
The experiments vary training-window length K and average results over 10 observation sets. Longer windows improve short-term forecasting, while the L05III forecast gains become less clear at large K.
6.6. Sensitivity on the length of the training window.
- Figure 3 plots mean curves across K from 10 experiments with different observation sets, alongside reference-model Lyapunov and power spectra.The figure follows Figure 2's structure.Source: Figure 3 · p. 21, §6.6 · p. 21
- For L96, forecast skill improves at all lead times as K increases; L05III shows the same tendency less clearly at large K, likely reflecting performance saturation.The paper attributes the saturation to L05III's weaker identifiability than L96.Source: §6.6 · p. 21
- For both models, Lyapunov spectra converge accurately for K ≥ 400 and power spectra for K ≥ 200, except at high frequencies.High-frequency power spectra are recovered better as K increases.Source: §6.6 · p. 21
6.7. Sensitivity on the lag of the smoother.
Smoother lag, observation noise, and monitoring density affect the surrogate models differently across L96 and L05III. The paper also frames the inference as jointly estimating the trajectory, dynamics, and model-error statistics.
6.7. Sensitivity on the lag of the smoother.
- L05III forecast skill is best at L = 4, while L96 is best at L = 2 or 4; filtering with L = 0 significantly degrades forecast skill in both.Lyapunov and power spectra improve with L ≥ 2, but only marginally.Source: §6.7 · p. 22
6.8. Sensitivity on the observation noise.
- As observation noise decreases, L96 forecast skill keeps improving, whereas L05III saturates at σy = 1; high-frequency power spectra improve for both models.The paper links the differing forecast responses to the models' identifiability under the surrogate architecture.Source: §6.8 · p. 22
6.9. Sensitivity on the density of the observations.
- With denser monitoring, L96 forecast skill improves from Ny = 36 to Ny = 40, while L05III forecast skill saturates and its Lyapunov and power spectra remain density-sensitive.The dependence of Lyapunov and power spectra on monitoring density is especially marked for L05III.Source: §6.9 · p. 22
- The proposed Bayesian view treats the state trajectory, dynamical model, and model-error statistics as joint control variables.The approach is intended to learn these quantities from partial and noisy observations.Source: §7 · p. 22
7. Conclusion.
The approach combines data assimilation and machine learning to estimate a surrogate model, state trajectory, and model-error statistics. It is tested on two low-order chaotic models, while applying the techniques to more complex, higher-dimensional models remains a technical challenge.
7. Conclusion.
- The experiments include full noisy observations for L96 and L05III and partial noisy observations with varied monitoring density.Figures 4 and 5 report forecast skill and power spectra for full observations, and additionally Lyapunov spectra for partial observations.Source: Figure 4 · p. 23, Figure 5 · p. 24
- The algorithm alternates an EnKS data-assimilation step that produces a posterior ensemble with an optimization step that updates the model and model-error statistics.The resulting representation includes a deterministic surrogate and a stochastic correction for residual model errors.Source: §7 · p. 23
- A full EM scheme that stores the ensemble is compared with an approximate scheme that computes model-error statistics online during the EnKS step.Both schemes are tested on two low-order chaotic models with distinct identifiability and across training-window, noise, monitoring-density, and smoother-lag settings.Source: §7 · p. 23
- Applying the techniques to more complex, higher-dimensional models remains a technical challenge; incorporating prior information and building online learning systems are open questions.These are among the remaining questions identified by the study.Source: §7 · pp. 23–24