Improving the prediction of complex nonlinear turbulent dynamical systems using nonlinear filter, smoother and backward sampling techniques

Abstract

Predicting complex nonlinear turbulent dynamical systems using partial observations is an important topic. Despite the simplicity of the forecast based on the ensemble mean time series, several critical shortcomings in the ensemble mean forecast and using path-wise measurements to quantify the prediction error are illustrated in this article. Then, a new ensemble method is developed for improving the long-range forecast. This new approach utilizes a mixture of the posterior distributions from data assimilation and is more skillful in predicting non-Gaussian statistics and extreme events than the traditional method by simply running the forecast model forward. Next, a systematic framework of improving forecast models is established, aiming at advancing the predictions at all ranges. The starting model in this new framework belongs to a rich class of nonlinear systems with conditional Gaussian structures. These models allow an efficient nonlinear smoother for state estimation using partial observations, which in turn facilitates a rapid parameter estimation based on an expectation–maximization algorithm. Conditioned on the partially observed time series, the nonlinear smoother further advances an efficient backward sampling of the hidden trajectories, the dynamical and statistical characteristics from which allow a systematic quantification of model error through information theory. The sampled trajectories then serve as the recovered observations of the hidden variables that promote the use of general nonlinear data-driven modeling techniques for a further improvement of the forecast model. A low-order model of the layered topographic equations with regime switching and rare events is used as a test example to illustrate this framework.

Appendix

1.1 An expectation–maximization algorithm for learning nonlinear models with partial observations

For the simplicity of discussion, define $\widehat{{\mathbf {u}}}_{\mathbf {I}}=\{{\mathbf {u}}_{\mathbf {I}}^0, \ldots , {\mathbf {u}}_{\mathbf {I}}^j,\ldots , {\mathbf {u}}_{\mathbf {I}}^J\}$ and $\widehat{{\mathbf {u}}}_\mathbf {II} = \{{\mathbf {u}}_\mathbf {II}^0, \ldots , {\mathbf {u}}_\mathbf {II}^j,\ldots , {\mathbf {u}}_\mathbf {II}^J\}$, where ${\mathbf {u}}_{\mathbf {I}}^j := {\mathbf {u}}_{\mathbf {I}}(t_j)$ and ${\mathbf {u}}_\mathbf {II}^j={\mathbf {u}}_\mathbf {II}(t_j)$. In other words, the values of ${\mathbf {u}}_{\mathbf {I}}$ and ${\mathbf {u}}_\mathbf {II}$ are taken at discrete points in time ${\mathbf {u}}_{\mathbf {I}}(t_j)$ and ${\mathbf {u}}_\mathbf {II}(t_j)$, for $j = 0,1,\ldots ,J$, where $T=J\varDelta {t}$ is the total length of the time series and $\varDelta {t}\ll 1$. This can be achieved by applying an Euler–Maruyama scheme [33, 54] to the continuous system and discretizing the observational data as well.

Given an ansatz of the conditional Gaussian nonlinear model (6) with observed variables ${\mathbf {u}}_{\mathbf {I}}$ and unobserved variables ${\mathbf {u}}_\mathbf {II}$, the goal here is to maximize the objective function, which is the log likelihood,

$$\begin{aligned} {\mathcal {L}}(\varvec{\theta }) = \log p(\widehat{{\mathbf {u}}}_{\mathbf {I}}|\varvec{\theta }) = \log \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} p (\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta }){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}, \end{aligned}$$

(57)

where $\varvec{\theta }$ is the collection of model parameters.

Using any distribution $Q(\widehat{{\mathbf {u}}}_\mathbf {II})$ over the hidden variables, a lower bound on the likelihood ${\mathcal {L}}$ can be obtained in the following way [38];

$$\begin{aligned} \log \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta }){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}= & {} \log \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \frac{p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta })}{Q(\widehat{{\mathbf {u}}}_\mathbf {II})}{\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}\nonumber \\\ge & {} \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \log \frac{p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta })}{Q(\widehat{{\mathbf {u}}}_\mathbf {II})}{\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}\nonumber \\= & {} \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \log p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta }){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II} \\&- \int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \log Q(\widehat{{\mathbf {u}}}_\mathbf {II}){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}\nonumber \\:= & {} {\mathcal {F}}(Q,\varvec{\theta }),\nonumber \end{aligned}$$

(58)

where the negative value of $\int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \log p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta }){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}$ is the so-called free energy, while $-\int _{\widehat{{\mathbf {u}}}_\mathbf {II}} Q(\widehat{{\mathbf {u}}}_\mathbf {II}) \log Q(\widehat{{\mathbf {u}}}_\mathbf {II}){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}$ is the entropy. Therefore, based on the fact ${\mathcal {F}}(Q,\varvec{\theta })\le {\mathcal {L}}(\varvec{\theta })$, it is clear that maximizing the log likelihood is equivalent to maximizing ${\mathcal {F}}$ alternatively with respect to the distribution Q and the parameters $\varvec{\theta }$. This can be achieved by the expectation–maximization (EM) algorithm [29, 99, 100],

$$\begin{aligned} \text{ E-Step: }&\qquad Q_{k+1}\leftarrow \arg \max _{Q}{\mathcal {F}}(Q,\varvec{\theta }_k), \end{aligned}$$

(59a)

$$\begin{aligned} \text{ M-Step: }&\qquad \varvec{\theta }_{k+1}\leftarrow \arg \max _{\varvec{\theta }}{\mathcal {F}}(Q_{k+1},\varvec{\theta }). \end{aligned}$$

(59b)

The maximization in the E-Step is reached when Q is exactly the conditional distribution of $\widehat{{\mathbf {u}}}_\mathbf {II}$ corresponding to the smoother estimates, that is,

$$\begin{aligned} Q_{k+1}(\widehat{{\mathbf {u}}}_\mathbf {II}) = p(\widehat{{\mathbf {u}}}_\mathbf {II}|\widehat{{\mathbf {u}}}_{\mathbf {I}},\varvec{\theta }_k). \end{aligned}$$

(60)

In such a situation, the bound in (58) becomes an equality ${\mathcal {F}}(Q,\varvec{\theta })={\mathcal {L}}(\varvec{\theta })$. Note that the conditional distribution in the E-Step is very difficult to solve for general nonlinear systems. Various numerical methods and approximations are often used [37, 38], which however may suffer from both the approximation errors and the curse of dimensionality. Nevertheless, for the conditional Gaussian systems, the distribution $p(\widehat{{\mathbf {u}}}_\mathbf {II}|\widehat{{\mathbf {u}}}_{\mathbf {I}},\varvec{\theta }_k)$ is given by the closed analytic formulae of the nonlinear smoother in Theorem 2, which greatly facilitates the application of the EM algorithm to many nonlinear models.

On the other hand, since the entropy (the second term on the right hand side of (58)) does not depend on $\varvec{\theta }$, the maximum in the M-Step is obtained by maximizing the negative of the free energy

$$\begin{aligned} \varvec{\theta }_{k+1}\leftarrow \arg \max _{\varvec{\theta }}\int _{\widehat{{\mathbf {u}}}_\mathbf {II}} p(\widehat{{\mathbf {u}}}_\mathbf {II}|\widehat{{\mathbf {u}}}_{\mathbf {I}},\varvec{\theta }_k) \log p(\widehat{{\mathbf {u}}}_{\mathbf {I}},\widehat{{\mathbf {u}}}_\mathbf {II}|\varvec{\theta }){\, \mathrm d}\widehat{{\mathbf {u}}}_\mathbf {II}. \end{aligned}$$

(61)

1.2 Parameter estimation of fully observed systems

This appendix provides more details of Step 4 of Sect. 5.1. Now assume a sampled trajectory of the unobserved variables ${\mathbf {u}}_\mathbf {II}$ is available. Assume after certain combination, the parameters only appear as prefactor coefficients in front of the nonlinear terms of the system. Therefore, the equation of ${\mathbf {u}}_\mathbf {II}$ can be rewritten as

$$\begin{aligned} {\, \mathrm d}{\mathbf {u}}_\mathbf {II} = {\mathbf {M}}({\mathbf {u}}_{\mathbf {I}},{\mathbf {u}}_\mathbf {II})\varvec{\theta }{\, \mathrm d}t + \text{ noise }, \end{aligned}$$

(62)

where ${\mathbf {M}}({\mathbf {u}}_{\mathbf {I}},{\mathbf {u}}_\mathbf {II})$ contains all the nonlinear functions that depend on ${\mathbf {u}}_{\mathbf {I}}$ and conditionally linear on ${\mathbf {u}}_\mathbf {II}$ as in the original equation. The noise coefficients are provided by the EM algorithm and are not updated here. The collection of the parameters $\varvec{\theta }$ can be estimated via a least squares estimator,

$$\begin{aligned} \varvec{\theta }= \left( \sum _j{\mathbf {M}}^*(t_j){\mathbf {M}}(t_j)\right) ^{-1}\left( \sum _j{\mathbf {M}}^*(t_j)\frac{{\, \mathrm d}{\mathbf {u}}_\mathbf {II}}{{\, \mathrm d}t}(t_j)\right) \end{aligned}$$

(63)

where $t_j$ represents discrete time instants.