Source-linked AI summary

Evaluating Data Assimilation Algorithms

K. J. H. Law, A. M. Stuart

arXiv:1107.4118v4physics.data-ancs.SYeess.SYmath.OCmath.PRphysics.ao-ph

TL;DR

Data-assimilation algorithms can be difficult to evaluate, and the Bayesian posterior is not computable in practical weather-forecasting scenarios. This paper uses the posterior as a gold standard and finds that approximate filters can reproduce the mean but typically perform poorly on covariance.

  • Problem

    The effects of data-assimilation algorithms can be hard to evaluate, while the Bayesian posterior is not computable in practical weather-forecasting scenarios.

    Source: 241. Introduction · p. 3
  • Method

    The paper evaluates data-assimilation algorithms against the Bayesian posterior probability distribution on the system state as a gold standard.

    Source: 241. Introduction · p. 3
  • Results

    Approximate filters can reproduce the desired distribution’s mean with appropriate parameters, but typically perform poorly at reproducing its covariance.

    Source: 241. Introduction · pp. 3, 7
  • Takeaways & Limitations

    Filters can be useful for predicting mean behavior but should be viewed with caution as predictors of uncertainty.

    Source: 241. Introduction · p. 3

    Because of the cost of computing the posterior, the study considers only one initial interval of observations.

    Source: 241. Introduction · p. 8

Abstract

from arXiv · show

Data assimilation leads naturally to a Bayesian formulation in which the posterior probability distribution of the system state, given the observations, plays a central conceptual role. The aim of this paper is to use this Bayesian posterior probability distribution as a gold standard against which to evaluate various commonly used data assimilation algorithms. A key aspect of geophysical data assimilation is the high dimensionality and low predictability of the computational model. With this in mind, yet with the goal of allowing an explicit and accurate computation of the posterior distribution, we study the 2D Navier-Stokes equations in a periodic geometry. We compute the posterior probability distribution by state-of-the-art statistical sampling techniques. The commonly used algorithms that we evaluate against this accurate gold standard, as quantified by comparing the relative error in reproducing its moments, are 4DVAR and a variety of sequential filtering approximations based on 3DVAR and on extended and ensemble Kalman filters. The primary conclusions are that: (i) with appropriate parameter choices, approximate filters can perform well in reproducing the mean of the desired probability distribution; (ii) however they typically perform poorly when attempting to reproduce the covariance; (iii) this poor performance is compounded by the need to modify the covariance, in order to induce stability. Thus, whilst filters can be a useful tool in predicting mean behavior, they should be viewed with caution as predictors of uncertainty. These conclusions are intrinsic to the algorithms and will not change if the model complexity is increased, for example by employing a smaller viscosity, or by using a detailed NWP model.

24 1. Introduction

The study evaluates practical data-assimilation algorithms against the Bayesian posterior as a gold standard in a 2D Navier–Stokes model, focusing on their reproduction of posterior means and covariances.

24 1. Introduction

  • Carefully chosen filters can predict mean behavior well, but typically perform poorly when reproducing uncertainty, including covariance.The introduction identifies this contrast as a central finding.Source: 241. Introduction · p. 3
  • The paper evaluates practical filters and variational methods against the Bayesian posterior distribution, which serves as an accurate but expensive gold standard.The posterior cannot be computed for practical weather-forecasting scenarios, but can assess approximate methods used in practice.Source: 241. Introduction · pp. 3, 6
  • The benchmark uses periodic 2D Navier–Stokes dynamics: a high-dimensional, complex system small enough for accurate offline posterior computation by statistical sampling.This allows the posterior to be computed accurately for comparison, despite its infeasibility for operational NWP.Source: 241. Introduction · p. 3
  • Particle filters approximate the desired distribution as particle counts increase, but high-dimensional degeneracy makes them impractical for geophysical data assimilation.The introduction notes ongoing research to address this limitation.Source: 241. Introduction · p. 4
  • In the Bayesian interpretation, 4DVAR corresponds to maximizing the posterior distribution, or finding a maximum a posteriori estimator.Variational methods choose an initial state to best match observations over a time window.Source: 241. Introduction · pp. 5–6
  • With appropriate parameter choices, approximate filters reproduce the posterior mean well, but typically perform poorly at reproducing covariance and should be treated cautiously as predictors of uncertainty.The authors characterize filters as useful for predicting mean behavior, not uncertainty.Source: 241. Introduction · p. 7
  • Filter covariance modifications needed for stability compound poor covariance performance; the authors say these limitations are intrinsic to the algorithms and persist with more complex models or smaller viscosity.The modifications are intended to induce filter stability and avoid divergence.Source: 241. Introduction · pp. 7–8
  • Because posterior computation is costly, the study evaluates one initial observation interval, does not address sensitivity to prior choice, and restricts comparison to the perfect-model scenario.Its conclusions assume a well-defined initial prior, observational error, and Bayesian posterior distribution.Source: 241. Introduction · p. 8
  • 4DVAR could perform better when using its approximate filtering distribution with the propagated Hessian, but that approach is too expensive in practice and constant covariance can limit posterior-variation reproduction.The limitation concerns reproducing the statistical variation of the posterior filtering distribution used as the next cycle’s prior.Source: 241. Introduction · p. 8

204 2. Statement of the Model

This section introduces the dynamical model and the filtering and smoothing problems posed by assimilating data into it. The discussion precedes discretization, with numerical implementation details deferred to the Appendix.

a. Dynamical Model: Navier-Stokes Equation

The paper uses a periodic two-dimensional incompressible Navier–Stokes model whose viscosity can span regimes with different predictability. It formulates assimilation as estimating a Bayesian posterior from noisy observations and explains why nonlinear dynamics can make Gaussian approximations misleading.

a. Dynamical Model: Navier-Stokes Equation

  • The model is the two-dimensional incompressible Navier–Stokes equation in a periodic box with side length two.The velocity field is projected onto divergence-free fields.Source: a. Dynamical Model: Navier-Stokes… · p. 10
  • Changing viscosity spans regimes from a stationary point to periodic, quasi-periodic, chaotic, and strongly chaotic behavior.Increasingly many positive Lyapunov exponents accompany increasing unstable-manifold dimension and decreasing predictability.Source: a. Dynamical Model: Navier-Stokes… · p. 10
  • This tunable predictability and the possibility of high-dimensional effective dynamics make the model useful for studying atmospheric data-assimilation issues.Source: a. Dynamical Model: Navier-Stokes… · p. 10

b. Inverse Problem

  • The inverse problem estimates the initial condition, or the state at time t, from noisy observations under the model dynamics.It is ill-posed when estimating the initial condition from later observations.Source: b. Inverse Problem · pp. 11–12
  • A Bayesian prior places a Gaussian random-field assumption on the initial condition, with prior mean and covariance encoding background state and variability across length scales.Source: b. Inverse Problem · p. 11
  • The observations measure velocity but not pressure, are equally spaced for simplicity, and have independent observational noise.Source: b. Inverse Problem · p. 11
  • With deterministic dynamics, the smoothing distribution at t = 0 maps forward to the exact filtering distribution.The corresponding forward map gives the predicting distribution; an inverse map can recover smoothing from filtering when it exists.Source: b. Inverse Problem · p. 12
  • For nonlinear dynamics, the posterior may require infinitely many samples to represent, and Gaussian approximations can become increasingly misleading as nonlinearity and attractor dimension increase.Sufficient nonlinearity may arise through the dynamics or the observation frequency.Source: b. Inverse Problem · p. 12

1) Smoothing

The paper formulates smoothing and filtering posteriors using Bayes’ rule, then uses MCMC sampling to approximate the posterior and compares it with variational and Gaussian-filter methods.

1) Smoothing

  • The initial-time smoothing posterior determines later-time posteriors by mapping it forward under the dynamics.Bayes’ rule characterizes the initial posterior relative to the prior.Source: 1) Smoothing · p. 13

2) Filtering

  • Filtering applies Bayes’ rule sequentially to obtain the distribution conditioned on observations through the current time.The prediction and update steps together yield the final filtering distribution.Source: 2) Filtering · p. 14

a. MCMC Sampling of the Posterior

318 3. Overview of Methods

  • The paper compares 3DVAR, Fourier Diagonal, Extended Kalman, and Ensemble Kalman filters, which use Gaussian approximations to represent state distributions.MCMC sampling approximates the Bayesian posterior, while 4DVAR computes its MAP estimator.Source: 3183. Overview of Methods · p. 15

b. Variational Methods: 4DVAR

The section describes 4DVAR and approximate Gaussian filters as computational approaches to posterior estimation. Their covariance handling ranges from a fixed 3DVAR choice to dynamics-based FDF updates.

b. Variational Methods: 4DVAR

  • 4DVAR uses the minimizer of the posterior functional as its approximation, and the paper solves this optimization with Newton’s method initialized by homotopy.When posterior variance is small, the MAP estimate can anchor a Gaussian approximation that is propagated forward, with linearized covariance.Source: b. Variational Methods: 4DVAR · p. 16

c. Approximate Gaussian Filters

  • High-dimensional PDE problems make exact distributional calculations extremely costly, motivating Gaussian approximations that make update formulas tractable.For nonlinear dynamics, a Gaussian filtering distribution need not remain Gaussian after prediction, although a Gaussian predicted distribution remains Gaussian after the update.Source: c. Approximate Gaussian Filters · p. 17

1) Constant Gaussian filter (3DVAR)

  • 3DVAR uses a constant covariance in its updates, with the initial-condition prior covariance as a natural theoretical choice.The authors note that other issues may require different choices.Source: 1) Constant Gaussian filter (3DVAR) · p. 18

2) Fourier Diagonal Filter (FDF)

  • The Fourier Diagonal Filter replaces nonlinear dynamics with a linear stochastic model to derive tractable covariances, including from climatological statistics.Its covariance converges to an explicitly computable limit, so asymptotically the method behaves like 3DVAR with a systematic covariance choice.Source: 2) Fourier Diagonal Filter (FDF) · pp. 18–19

3) Low Rank Extended Kalman Filter (LRExKF)

The LRExKF reduces the cost of extended-Kalman covariance propagation with a low-rank approximation, while the EnKF uses an ensemble to estimate filtering statistics. The accuracy assessment compares these methods with posterior moments and excludes model error.

3) Low Rank Extended Kalman Filter (LRExKF)

  • Extended-Kalman covariance propagation requires tangent-linear and adjoint solves for every system dimension, so the paper uses a low-rank approximation to reduce this high-dimensional cost.The method retains dominant covariance eigenpairs; when eigenvalues are well separated, a small number captures most of the action.Source: 3) Low Rank Extended Kalman Filter… · p. 21

4) Ensemble Kalman Filter (EnKF)

  • The EnKF updates an ensemble of particles and uses its empirical mean and covariance; the comparisons instead use the analytical covariance update.For small ensembles, the analytical and sample covariances can differ substantially, with either positive or negative effects on filter divergence.Source: 4) Ensemble Kalman Filter (EnKF) · p. 22

464 4. Filter Accuracy

  • The accuracy assessment compares algorithms with posterior distributions and truth, characterizing the posterior by its mean, variance, and MAP estimate.Model error is excluded: uncertainty in the dynamical model comes only through the initial condition.Source: 4644. Filter Accuracy · p. 23

a. Nature of Approximations

The paper identifies Gaussian approximation and, for EnKF, sampling as sources of error, then assesses filters across stationary and strongly chaotic regimes. Most methods reproduce posterior means better than covariance, with accuracy depending on method and regime.

a. Nature of Approximations

  • Gaussian approximations become biased when dynamics are sufficiently nonlinear or covariance becomes large, while EnKF also incurs sampling error when its ensemble is too small.The paper expects these approximation-breaking factors to grow with Reynolds number and time between observations; EnKF sample requirements also tend to increase with dimension.Source: a. Nature of Approximations · p. 24
  • The ExKF performs well in the stationary regime and under high-frequency observations in the strongly chaotic regime, but fails to recover accurate mean and covariance when chaotic observations are sparse.In the latter setting, the Gaussian approximation is no longer valid.Source: a. Nature of Approximations · p. 25

b. Illustration via Two Regimes

  • Most algorithms reproduce posterior means reasonably well, but algorithms other than LRExKF and 4DVAR perform poorly on posterior variance.LRExKF and 4DVAR reproduce variance reasonably well when the Reynolds number is sufficiently small and/or observations are frequent; otherwise ad hoc filter approximations can fail to recover covariance accurately.Source: b. Illustration via Two Regimes · pp. 26–27
  • LRExKF usually has the smallest error relative to the posterior mean and sometimes accurately recovers variance, while variance errors can overestimate some wavenumber modes and underestimate others.The posterior smoothing distribution can become noticeably non-Gaussian while remaining unimodal, even as the filtering distribution stays close to Gaussian.Source: b. Illustration via Two Regimes · p. 27

c. Stationary Regime

In the stationary regime, the posterior mean and MAP estimate on the initial state are rougher than the truth, but the solution dynamics quickly smooth this structure and contract errors by the final time. As observation intervals increase, 4DVAR and 3DVAR show distinct trends in mean and variance error.

c. Stationary Regime

  • The initial posterior mean and MAP estimate are rougher than the truth, but the MAP is indistinguishable from the posterior mean in the illustrated case.The roughness is attributed to the prior mean; the solution operator damps high wavenumbers and removes it.Source: c. Stationary Regime · p. 28
  • The truth, posterior distribution, and MAP estimator are mutually closer at t = T than at t = 0, as errors contract under dynamics attracted to a unique stationary state.The filtering-distribution mean at t = T likewise appears indistinguishable from the truth and MAP estimator.Source: c. Stationary Regime · p. 28
  • With T = 2 and h increasing through 0.2, 1, and 2, 4DVAR smoothing mean error is much smaller at t = T than t = 0, increasing at T and decreasing at t = 0.The 4DVAR mean and variance errors at t = 0 and t = T are close to or below the MCMC accuracy threshold.Source: c. Stationary Regime · pp. 28–29
  • Comparisons use the stationary regime with ν = 0.1, filtering time T = 2, and observation intervals h = 0.2, 1, and 2.The basic time step is dt = 0.05.Source: c. Stationary Regime · p. 28

d. Strongly Chaotic Regime

In the strongly chaotic regime, uncertainty can persist or increase from initial to final time, and filter performance depends on observation interval and method. FDF recovers the mean consistently as h increases, while LRExKF degrades at longer intervals, especially when nonlinear effects become prominent.

d. Strongly Chaotic Regime

  • The strongly chaotic regime uses ν = 0.01, with filtering times T = 0.2 or 1 and observation intervals h = 0.02, 0.1, 0.2, or 0.5.The dynamics are described as more nonlinear and less predictable, with a high-dimensional attractor spanning many scales.Source: d. Strongly Chaotic Regime · pp. 29–30
  • In this regime, truth and MAP profiles are indistinguishable from the mean throughout the time interval, and the prior mean's influence is not noticeable initially.Solutions are described as rougher than in the stationary regime.Source: d. Strongly Chaotic Regime · p. 29
  • Uncertainty in u11 barely decreases from t = 0 to t = T, while uncertainty in the moderately high mode u55 increases.The passage attributes the regime's contrast with near-stationary behavior primarily to a large-dimensional attractor and unpredictable dynamics.Source: d. Strongly Chaotic Regime · pp. 29–30
  • When computable, 4DVAR smoothing variance error is smaller at t = 0 than at t = T, while its mean error at t = 0 for h ≤ 0.1 is below MCMC accuracy threshold.Variance error at T = 1 cannot be computed accurately because long-time accumulated error affects the discretized analytical adjoint approximation.Source: d. Strongly Chaotic Regime · p. 30
  • At h = 0.5 and T = 1, FDF performs best; its mean recovery is consistent as h increases while other algorithms deteriorate, and its variance error decreases with h.The FDF variance error is very large when Q is an order of magnitude larger than Γ; the reported favorable properties are associated with stability.Source: d. Strongly Chaotic Regime · p. 30
  • LRExKF has substantial error at T = 1 and majorly fails at h = 0.5, where nonlinear effects become prominent and the linear approximation is not expected to remain valid.At h = 0.1 and T = 0.2, its mean error first increases and approaches its variance error, attributed to bias from neglecting the next-order dynamics expansion.Source: d. Strongly Chaotic Regime · p. 30

650 5. Filter Stability

This section examines how variance inflation and other stabilization choices affect filter stability and mean and variance accuracy across time horizons, filters, and observation times.

Filter Stability

  • Variance inflation lets filters run longer before diverging, but may worsen mean accuracy before divergence and variance predictions.The section describes the mechanism and general effects of stabilization rather than claiming an optimal inflation implementation.Source: 6505. Filter Stability · p. 31
  • Stability balances observational correction against unstable model dynamics; increasing forecast covariance moves the Kalman gain toward the identity and strengthens observation weighting.Growing directions are typically associated with overestimated variance during stable assimilation, while unstable cases include times when variance is underestimated.Source: 6505. Filter Stability · p. 32
  • Filter instability begins when the estimated mean’s distance from truth persistently exceeds the estimated standard deviation.The section also states that observations must be sufficient to control instabilities inherent in the model dynamics for this distance to remain small.Source: 6505. Filter Stability · pp. 31, 33
  • The average filter error is set by observation-error scale, while covariance updates—and thus the Kalman gain—affect its exact size.The section treats the observation-only filter’s error as an upper bound for a reasonable filter and uses it to define instability.Source: 6505. Filter Stability · p. 33
  • In the strongly chaotic 3DVAR example, stabilization keeps error below the desired bound and runs longer; stabilized LRExKF shows similar behavior.Both simulations exceed the lower error bound, but only the unstable algorithm exceeds the desired bound; the stabilized model also tracks a high-wavenumber trajectory better than observations alone.Source: 6505. Filter Stability · pp. 34–35
  • At T = 0.2, stabilized algorithms are less accurate for the mean at both observation times, but at T = 1.0 inflating variance reverses this behavior by enhancing stability.The passage attributes the short-horizon loss of accuracy to variance inflation.Source: 6505. Filter Stability · p. 35
  • At T = 1.0 and h = 0.2, stabilized 3DVAR outperforms its original version, whereas stabilized EnKF and LRExKF are less accurate than their original versions.These outcomes show that the effect of stabilization differs among filters at this observation time.Source: 6505. Filter Stability · p. 35
  • At h = 0.5 and T = 1.0, the stabilized versions improve upon the original algorithms in all three cases; without stabilization, FDF outperforms the other three filters and their stabilized counterparts.Source: 6505. Filter Stability · p. 35

770 6. Conclusion

The paper evaluates assimilation algorithms against a Bayesian posterior computed for 2D Navier–Stokes, comparing their mean and covariance. Most methods reproduce the mean reasonably well, but covariance estimates are often inaccurate; stabilization can extend useful mean prediction while worsening other accuracy.

770 6. Conclusion

  • The evaluation compares algorithm means and covariance structures with the Bayesian posterior for strongly constrained 2D Navier–Stokes regimes, including stationary and strongly chaotic dynamics.The setup makes the posterior unimodal and approximately Gaussian, allowing its distribution to be sampled using MCMC.Source: 7706. Conclusion · pp. 36–37
  • Most studied filtering and variational algorithms reproduce the mean reasonably well, but often fail to reproduce covariance accurately when their approximations are unjustified.The paper also reports that most filtering algorithms can become unstable over longer intervals, losing accuracy even in mean prediction.Source: 7706. Conclusion · p. 37
  • Variance inflation can stabilize filters over longer windows and improve long-term mean accuracy, but may reduce short-window accuracy and makes covariance prediction impossible.The conclusions are attributed to approximations used to make online algorithms tractable and are not expected to change with different dynamical models or algorithm parameters.Source: 7706. Conclusion · p. 37
  • The authors recommend caution when using practical assimilation algorithms to make uncertainty claims, while noting that properly tuned filters often track the signal mean over fairly long intervals.They attribute the limits to tractability-driven algorithmic approximations, not to the chosen dynamics or parameter values.Source: 7706. Conclusion · p. 37
  • Choosing the prior distribution on the initial condition remains an important issue not addressed by the paper, and initializing statistical forecasting algorithms remains challenging.The authors note that further investigation has conceptual and practical challenges.Source: 7706. Conclusion · pp. 37–38
  • Uniform component-wise rank histograms prevent consistent variance underestimation in any direction and are necessary, but not sufficient, for accurate posterior approximation.The paper illustrates the gap with a 3DVAR estimator constructed to have uniform rank histograms even though its constant prior does not match the actual prior.Source: 7706. Conclusion · p. 38
  • The FDF variance recovers posterior-variance decay but is about an order of magnitude larger, and it overestimates high-variance directions while underestimating low-variance ones.The high- and low-variance directions correspond, respectively, to unstable and stable directions in this setup.Source: 7706. Conclusion · p. 38

852 7. Appendix: Some numerical details

The appendix documents the numerical solvers, optimization, filtering, and sampling procedures used to compute the posterior and compare algorithms. It also describes the computational costs and accuracy limits that constrain those computations.

Model setup

  • The prior is a smoothness prior, and the observational covariance is chosen to emphasize low wavenumbers where the unstable manifold has support.The prior-to-observation covariance ratio tends to zero as |k| increases.Source: 8527. Appendix: Some numerical details · p. 39

Dynamical model

  • The dynamical solver uses a Galerkin approximation in a divergence-free Fourier basis, with ETD4RK and FFT-based nonlinear computations.The heat semigroup and Duhamel’s principle handle diffusion, while mode truncation prevents aliasing.Source: 8527. Appendix: Some numerical details · p. 40
  • The numerical objective is to resolve unstable, low-wavenumber modes and achieve reproducible high temporal resolution, rather than high spatial resolution.The authors seek to ensure that trajectory divergence reflects model instability rather than the numerical time-stepping scheme.Source: 8527. Appendix: Some numerical details · p. 40

Adjoint solver

  • Adjoint accuracy depends on the forward solver; over longer intervals, the linearization becomes less accurate and accurate gradients and Hessians cannot be computed.The adjoint uses the same time-stepper with twice the forward time-step, and the implementation follows an “optimize then discretize” approach.Source: 8527. Appendix: Some numerical details · pp. 40–41

Linear algebra and optimization

  • LRExKF’s low-rank covariance approximation captures more than 99% of the full-rank version in Frobenius norm.GMRES solves the Newton-method linear systems for 4DVAR, while Arnoldi supplies low-rank covariance approximations for LRExKF and 4DVAR.Source: 8527. Appendix: Some numerical details · p. 41

MCMC computations

  • About 1 s per sample implies roughly 10 days for millions of serial MCMC samples, making posterior sampling the most costly algorithm discussed.The authors accelerate burn-in with the 4DVAR MAP estimate and run independent batches in parallel; early rejection can stop rejected proposals before the assimilation window ends.Source: 8527. Appendix: Some numerical details · pp. 41–42

MCMC convergence

  • With 3.2 × 10^5 samples per chain, the eight-chain variance discrepancy is O(10^-2) and the mean discrepancy is O(10^-3).The appendix cautions that sampling error and cost limit reliable precision: relative norms below 10^-2 for variance or 10^-3 for mean cannot be relied on, and more than 10^7 samples is not tractable.Source: 8527. Appendix: Some numerical details · p. 42
Loading 1107.4118v4…