## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Beating the curse of dimension with accurate statistics for the Fokker–Planck equation in complex turbulent systems

Contributed by Andrew J. Majda, October 13, 2017 (sent for review September 28, 2017; reviewed by Joseph Tribbia and Xiaoming Wang)

## Significance

Solving the Fokker–Planck equation for high-dimensional complex dynamical systems is an important issue. Effective strategies are developed and incorporated into efficient statistically accurate algorithms for solving the Fokker–Planck equations associated with a rich class of high-dimensional nonlinear turbulent dynamical systems with strong non-Gaussian features. These effective strategies exploit a judicious block decomposition of high-dimensional conditional covariance matrices and statistical symmetry to facilitate an extremely efficient parallel computation and a significant reduction of sample numbers. The resulting algorithms can efficiently solve the Fokker–Planck equation in much higher dimensions even with orders in the millions and thus beat the curse of dimension. Skillful behavior of the algorithms is illustrated for highly non-Gaussian systems in excitable media and geophysical turbulence.

## Abstract

Solving the Fokker–Planck equation for high-dimensional complex dynamical systems is an important issue. Recently, the authors developed efficient statistically accurate algorithms for solving the Fokker–Planck equations associated with high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures, which contain many strong non-Gaussian features such as intermittency and fat-tailed probability density functions (PDFs). The algorithms involve a hybrid strategy with a small number of samples L, where a conditional Gaussian mixture in a high-dimensional subspace via an extremely efficient parametric method is combined with a judicious Gaussian kernel density estimation in the remaining low-dimensional subspace. In this article, two effective strategies are developed and incorporated into these algorithms. The first strategy involves a judicious block decomposition of the conditional covariance matrix such that the evolutions of different blocks have no interactions, which allows an extremely efficient parallel computation due to the small size of each individual block. The second strategy exploits statistical symmetry for a further reduction of L. The resulting algorithms can efficiently solve the Fokker–Planck equation with strongly non-Gaussian PDFs in much higher dimensions even with orders in the millions and thus beat the curse of dimension. The algorithms are applied to a 1,000-dimensional stochastic coupled FitzHugh–Nagumo model for excitable media. An accurate recovery of both the transient and equilibrium non-Gaussian PDFs requires only

- high-dimensional non-Gaussian PDFs
- hybrid strategy
- block decomposition
- statistical symmetry
- small sample size

The Fokker–Planck equation is a partial differential equation (PDE) that governs the time evolution of the probability density function (PDF) of a complex system with noise (1, 2). For a general nonlinear dynamical system,**2** involves strong non-Gaussian features with intermittency and extreme events (3⇓–5). In addition, the dimension of 𝐮 in these complex systems is typically very large, representing a variety of variability in different temporal and spatial scales (3, 6). Therefore, solving the high-dimensional Fokker–Planck equation for both the steady-state and transient phases with non-Gaussian features is an important issue. However, traditional numerical methods such as finite element and finite difference as well as the direct Monte Carlo simulations of Eq. **1** all suffer from the curse of dimension (7, 8).

Recently, the authors developed efficient statistically accurate algorithms for solving the Fokker–Planck equation associated with high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures (9). These conditional Gaussian nonlinear dynamical systems capture many strong non-Gaussian features such as intermittency and fat-tailed PDFs (10). Applications of the conditional Gaussian framework include modeling and predicting the highly intermittent time series of the Madden–Julian oscillation and monsoon (11⇓–13), state estimation of the turbulent ocean flows from noisy Lagrangian tracers (14⇓–16), dynamic stochastic superresolution of sparsely observed turbulent systems (17), and stochastic superparameterization for geophysical turbulent flows (18), etc. The efficient statistically accurate algorithms in ref. 9 involve a hybrid strategy that requires only a small number of samples. In these algorithms, a conditional Gaussian mixture in the high-dimensional subspace of

In this article, two effective strategies are developed and incorporated into the algorithms in ref. 9 (hereafter, basic algorithms) that enable the expanded algorithms to efficiently solve the Fokker–Planck equation in much higher dimensions, even with orders in the millions. In fact, the major computational cost in the basic algorithms for systems with a large dimension of state variables **1** consist of any nonlinear interactions between different components of

The second effective strategy exploits the statistical symmetry if the dynamical system in Eq. **1** is statistically homogeneous; namely, the statistics of different **1** represent a discrete approximation of some PDEs in a periodic domain with nonlinear advection, diffusion, and homogeneous external forcing. Examples include a rich class of models in geophysical turbulence and excitable media (4, 21). In light of statistical symmetry, the number of samples L in the algorithms can be greatly reduced. In fact, the effective sample size of each

These effective strategies are incorporated into the basic algorithms and then applied to two highly non-Gaussian dynamical systems. The first model is a 1,000-dimensional stochastic coupled FitzHugh–Nagumo (FHN) model for excitable media with extreme events (4, 22). It describes activation and deactivation dynamics of spiking neurons and has scale-invariant features. The block decomposition leads to solving the evolution of 500 individual covariance matrices and the statistical symmetry allows an accurate estimation of both the transient and the steady-state PDFs using only

The remainder of this article includes a detailed description of these effective strategies and their application to the stochastic coupled FHN and two-layer L96 models in solving the highly non-Gaussian PDFs at both the transient and statistical equilibrium phases, using a small number of samples.

## Algorithms

### Conditional Gaussian Framework.

The general framework of conditional Gaussian models is given as (10, 25)**3**, **3** are named as conditional Gaussian systems due to the fact that once **3** remains highly nonlinear and is able to capture many non-Gaussian features as observed in nature (10). One of the desirable features of Eq. **3** is that the conditional distribution in Eq. **4** has the following closed analytical form (25):

### Basic Algorithms.

Here, we summarize the basic efficient statistically accurate algorithms developed in ref. 9. First, we generate L independent trajectories of the variables **5**,**6** (as well as Eqs. **7** and **8** below) is taken to illustrate the statistical intuition, while the estimator is the nonasymptotic version. On the other hand, a Gaussian kernel density estimation method is used for solving the PDF of the observed variables **6** and **7**, a hybrid method is applied to solve the joint PDF of **3**], which is computationally affordable. In addition, the closed form of the L conditional distributions in Eq. **6** can be solved in a parallel way due to their independence (9), which further reduces the computational cost.

### Beating the Curse of Dimension with Block Decomposition.

The basic algorithms succeed in solving the Fokker–Planck equation with

Consider the following decomposition of state variables**3** are also decomposed into K groups, where the variables on the left-hand side of the kth group are

To develop efficient statistically accurate algorithms that beat the curse of dimension, the following two conditions are imposed on the coupled system:

#### Condition 1:

In the dynamics of each **3**, the terms

#### Condition 2:

The initial values of

*Conditions 1* and *2* are not artificial and they are actually the salient features of many complex systems with multiscale structures (18), multilevel dynamics (20), or state-dependent parameterizations (17). Under these two conditions, the conditional covariance matrix becomes block diagonal, which can be easily verified according to Eq. **5b**. The evolution of the conditional covariance of

Next, the structures of **9** allow the coupling among all of the K groups of variables in the conditional mean according to Eq. **5a**. The evolution of

### Statistical Symmetry.

The computational cost in the algorithms developed above can be further reduced if the coupled system Eq. **3** has statistical symmetry; namely,

With the statistical symmetry, collecting the conditional Gaussian ensembles *SI Appendix* provides the mathematical details of reconstructing the joint PDFs, using statistical symmetry.

## A Stochastic Coupled FHN Model

The efficient statistically accurate algorithms developed in this article work for a wide class of models in excitable media (4), including different versions of the famous FHN model that describes activation and deactivation dynamics of spiking neurons. Here, the algorithms are applied to a stochastic coupled FHN model with N elements (4, 22),**10**, the timescale ratio **10** exhibits rich dynamical behaviors. Below, we adopt constant parameters and the initial values are

### Model Behavior in Different Dynamical Regimes.

Fig. 1 shows the model behavior in three different regimes with *A*, the spatial–temporal patterns are highly coherent due to the choice of a weak noise *B*), the coherent patterns becomes much weaker and only quasi-regular periods are found in the time series of *C*), the spatial–temporal pattern becomes strongly mixed and the time series is more irregular. Correspondingly, the PDF of

It is shown in *SI Appendix* that the FHN model Eq. **10** is scale invariant in all three regimes. The scale-invariant structure means that the spatial–temporal structures in any given scale change little as the number of spatial grid points N increases. Mutual information (30) is used to quantify the dependence of different variables with strongly non-Gaussian features, the advantage of which over pattern correlation is also clearly illustrated in *SI Appendix*.

### Recovering the PDFs at Both Transient and Equilibrium Phases Exploiting Statistical Symmetry.

Now we apply the efficient statistically accurate algorithms to solve the PDFs associated with the stochastic coupled FHN system in Eq. **10**. Here we focus on the weakly coherent regime (Fig. 1*B*). Due to the statistical symmetry, the effective sample size is

The time evolutions of the first four moments associated with

Fig. 4 illustrates the recovered two-point statistics at *A* and *B* shows the PDFs

## Two-Layer Inhomogeneous L96 Model

The two-layer L96 model is a conceptual model in geophysical turbulence that is widely used as a testbed for data assimilation and parameterization in numerical weather forecasting (20, 23, 24). The model can be regarded as a coarse discretization of atmospheric flow on a latitude circle with complicated wave-like and chaotic behavior. It schematically describes the interaction between small-scale fluctuations with larger-scale motions. In the model presented here, large-scale motions are denoted by variables **11** is that the nonlinear interaction between **3** with **11a** while both the damping **11** has many desirable properties as in more complicated turbulent systems. Particularly, the smaller scales are more intermittent with stronger fat tails in PDFs. See *SI Appendix* for more details.

Although the two-layer inhomogeneous L96 model in Eq. **11** has no statistical symmetry, the model structure nevertheless allows the effective block decomposition. Below, **11** to implement the efficient statistically accurate algorithms. As a comparison, a direct Monte Carlo method requires

Fig. 6 shows the one-point statistics as a function of i regarding the first four moments of *SI Appendix*). These results indicate the success of the efficient statistically accurate algorithms in capturing the inhomogeneous behavior of the model. Fig. 7 demonstrates the skill of recovering different 2D joint PDFs. Fig. 7 *A* and *B* shows the PDFs *C* shows the joint PDFs of the two smallest-scale fluctuation variables *D* shows the joint PDFs

## Concluding Discussion

Effective strategies involving block decomposition of the conditional covariance matrix and statistical symmetry are developed and incorporated into the efficient statistically accurate algorithms in ref. 9. The resulting expanded algorithms are able to efficiently solve the Fokker–Planck equation in much higher dimensions even with orders in the millions and thus beat the curse of dimension. Applications of these effective strategies to both the stochastic coupled FHN and the two-layer inhomogeneous L96 models illustrate the efficiency and accuracy.

It is worthwhile pointing out that although only the recovered one-point and two-point statistics are shown in this article for illustration purposes, the algorithms can actually provide an accurate estimation of the full joint PDF of

## Acknowledgments

The research of A.J.M. is partially supported by the Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI) Grant N0001416-1-2161 and the New York University Abu Dhabi Research Institute. N.C. is supported as a postdoctoral fellow through A.J.M.’s ONR MURI Grant.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: jonjon{at}cims.nyu.edu or chennan{at}cims.nyu.edu.

Author contributions: N.C. and A.J.M. designed research, performed research, and wrote the paper.

Reviewers: J.T., National Center for Atmospheric Research (NCAR); and X.W., Florida State University.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1717017114/-/DCSupplemental.

- Copyright © 2017 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Gardiner CW

- ↵
- Risken H

*Methods of Solution and Applications*, Springer Series in Synergetics, ed Haken H (Springer, Berlin), Vol 18. - ↵
- Majda A

*Introduction to Turbulent Dynamical Systems in Complex Systems*, Frontiers in Applied Dynamical Systems: Reviews and Tutorials 5 (Springer, Cham, Switzerland). - ↵
- ↵
- Cousins W,
- Sapsis TP

- ↵
- Vallis GK

- ↵
- Papadrakakis M,
- Stefanou G,
- Papadopoulos V

- Pichler L,
- Masud A,
- Bergman LA

- ↵
- Robert CP

- ↵
- ↵
- Chen N,
- Majda AJ

- ↵
- Chen N,
- Majda AJ,
- Giannakis D

- ↵
- Chen N,
- Majda AJ

- ↵
- Chen N,
- Majda AJ

- ↵
- Chen N,
- Majda AJ,
- Tong XT

- ↵
- Chen N,
- Majda AJ,
- Tong XT

- ↵
- Chen N,
- Majda AJ

- ↵
- Branicki M,
- Majda AJ

- ↵
- Majda AJ,
- Grooms I

- ↵
- Chen N,
- Majda AJ,
- Tong XT

- ↵
- ↵
- Majda A,
- Wang X

- ↵
- Muratov CB,
- Vanden-Eijnden E,
- Weinan E

- ↵
- Lee Y,
- Majda AJ

- ↵
- ↵
- Liptser RS,
- Shiryaev AN

- ↵
- ↵
- Majda AJ,
- Harlim J

- ↵
- Guckenheimer J,
- Holmes PJ

- ↵
- DeVille RL,
- Vanden-Eijnden E,
- Muratov CB

- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Mathematics