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

aDepartment of Mathematics, Courant Institute of Mathematical Sciences, New York University, New York, NY 10012;bCenter for Atmosphere Ocean Science, Courant Institute of Mathematical Sciences, New York University, New York, NY 10012;

aDepartment of Mathematics, Courant Institute of Mathematical Sciences, New York University, New York, NY 10012;bCenter for Atmosphere Ocean Science, Courant Institute of Mathematical Sciences, New York University, New York, NY 10012;cCenter for Prototype Climate Modeling, New York University Abu Dhabi, Abu Dhabi, United Arab Emirates

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 L=1 samples! In addition, the block decomposition facilitates the algorithms to efficiently capture the distinct non-Gaussian features at different locations in a 240-dimensional two-layer inhomogeneous Lorenz 96 model, using only L=500 samples.

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,d??=??(??,t)dt+??(??,t)d??,[1]with state variables ??∈?N, noise matrix ??∈?N×K, and white noise ??∈?K, the associated Fokker–Planck equation is given by??tp(??,t)=????(??(??,t)p(??,t))+12???????(??(??,t)p(??,t)),pt|t=t0=p0(??),[2]with ??=????T. In many complex dynamical systems, such as geophysical and engineering turbulence, neuroscience, and excitable media, the solution of the Fokker–Planck equation in Eq. 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 ?????? via an extremely efficient parametric method is combined with a judicious Gaussian kernel density estimation in the remaining low-dimensional subspace of ????. Particularly, the parametric method provides closed analytical formulas for determining the conditional Gaussian distributions in the high-dimensional subspace of ?????? and is therefore computationally efficient and accurate. It has been shown in a stringent set of numerical tests (9) that with an order of L～O(100) samples the mixture distribution has significant skill in capturing both the statistically steady state and the transient behavior with fat tails of non-Gaussian PDFs in up to six dimensions. Rigorous analysis (19) indicates that L does not increase exponentially as the dimension of the high-dimensional subspace of ?????? to maintain a given level of accuracy, which is fundamentally different from Monte Carlo methods.

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 ?????? comes from solving the time evolution of the conditional covariance. To overcome this difficulty, an effective strategy involving a decomposition of the state variables into different groups ??=∪k=1K??k with each ??k=(????,k,??????,k)∈?N??,k+N????,k is developed. Then the high-dimensional conditional covariance matrix of ?????? conditioned on ???? becomes a block diagonal matrix under the conditions that the nonlinear terms on the right-hand side of ??k in Eq. 1 consist of any nonlinear interactions between different components of ???? but only those between ??????,k and nonlinear functions of ????,k. Note that such conditions are not artificial and they are actually salient features of many complex dynamical systems with multiscale structures (18), multilevel dynamics (20), or state-dependent parameterizations (17). One important characteristic of the resulting covariance matrix is that the evolution of the kth block, representing the conditional covariance of ??????,k conditioned on ????, has no interactions with that of ??????,k′ for all k′≠k. This allows an extremely efficient parallel computation due to the small size of each individual block. On the other hand, the conditional means of different ??????,k are all coupled in their evolutions.

The second effective strategy exploits the statistical symmetry if the dynamical system in Eq. 1 is statistically homogeneous; namely, the statistics of different ??k are identical with each other. The statistical symmetry is often satisfied when the underlying dynamics in Eq. 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 ??k becomes L′=LK and therefore a much smaller L is needed to reach the same level of accuracy as in the situation without using statistical symmetry.

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 L=1 samples! The second model is the so-called two-layer Lorenz 96 (L96) model in geophysical turbulence (20, 23, 24), which is widely used as a testbed for data assimilation and parameterization in numerical weather forecasting. Inhomogeneous damping and coupling are adopted in the model with 240 state variables that mimic the atmosphere motion along a latitude circle with different dissipation over land and sea. Despite the absence of statistical symmetry, the block decomposition facilitates the algorithms to efficiently capture the distinct non-Gaussian features at different locations, using only L=500 samples.

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)d????=[??0(t,????)+??1(t,????)??????]dt+????(t,????)d????(t),[3a]d??????=[??0(t,????)+??1(t,????)??????]dt+??????(t,????)d??????(t),[3b]where ??=(????,??????) with both ????∈RN?? and ??????∈RN???? being multidimensional state variables. In Eq. 3, ??0,??1,??0,??1,???? and ?????? are vectors and matrices that depend only on time t and the state variables ????, and ????(t) and ??????(t) are independent Wiener processes. The systems in Eq. 3 are named as conditional Gaussian systems due to the fact that once ????(s) for s≤t is given, ??????(t) conditioned on ????(s) becomes a Gaussian process with mean ??ˉ????(t) and covariance ??????(t); i.e.,p(??????(t)|????(s≤t))～N(??ˉ????(t),??????(t)).[4]Despite the conditional Gaussianity, the coupled system in Eq. 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):d??ˉ????(t)=[??0(t,????)+??1(t,????)??ˉ????]dt+(????????1?(t,????))×(?????????)?1(t,????)[d?????(??0(t,????)+??1(t,????)??ˉ????)dt],[5a]d??????(t)={??1(t,????)??????+????????1?(t,????)+(?????????????)(t,????)?????????1?(t,????)(?????????)?1(t,????)(????????1?(t,????))?}dt.[5b]

Basic Algorithms.

Here, we summarize the basic efficient statistically accurate algorithms developed in ref. 9. First, we generate L independent trajectories of the variables ????, namely ????1(s≤t),…,????L(s≤t). Then, different strategies are used to deal with ???? and ??????. The PDF of ?????? is estimated via a parametric method that exploits the closed form of the conditional Gaussian statistics in Eq. 5,p(??????(t))=limL→∞1L∑i=1Lp(??????(t)|????i(s≤t)).[6]Note that the limit L→∞ in Eq. 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 ????,p(????(t))=limL→∞1L∑i=1LK??(????(t)?????i(t)),[7]where K??(?) is a Gaussian kernel centered at each sample point ????i(t) with covariance given by the bandwidth matrix ??(t). The kernel density estimation algorithm here involves a “solve-the-equation plug-in” approach for optimizing the bandwidth (26) that works for any non-Gaussian PDFs. Finally, combining Eqs. 6 and 7, a hybrid method is applied to solve the joint PDF of ???? and ?????? through a Gaussian mixture,p(????(t),??????(t))=limL→∞1L∑i=1L(K??(????(t)?????i(t))?p(??????(t)|????i(s≤t))).[8]Practically, L～O(100) is sufficient for the hybrid method to solve the joint PDF with N??≤3 and N????～10. Since L is small, the trajectories ????1(s≤t),…,????L(s≤t) can be obtained by running a Monte Carlo simulation for the coupled system [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 O(10) state variables. Now we develop an effective strategy with block decomposition and incorporate it into the basic algorithms. The expanded algorithms can efficiently solve the Fokker–Planck equation in much higher dimensions even with orders in the millions and beat the curse of dimension.

Consider the following decomposition of state variables??k=(????,k,??????,k)with????,k∈?N??,kand??????,k∈?N????,k,where 1≤k≤K, N??=∑k=1KN??,k, and N????=∑k=1KN????,k. Correspondingly, the full dynamics in Eq. 3 are also decomposed into K groups, where the variables on the left-hand side of the kth group are ??k. In addition, we assume both ???? and ?????? are diagonal for notation simplicity.

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 ??k in Eq. 3, the terms ??0,k and ??0,k can depend on all of the components of ???? while the terms ??1,k and ??1,k are only functions of ????,k; namely,??0,k???0,k(t,????),???0,k???0,k(t,????),??1,k???1,k(t,????,k),???1,k???1,k(t,????,k).[9]In addition, only ??????,k interacts with ??1,k and ??1,k on the right-hand side of the dynamics of ??k.

Condition 2:

The initial values of (????,k,??????,k) and (????,k′,??????,k′) with k≠k′ are independent from each other.

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 ??????,k conditioned on ???? is given byd??????,k(t)={??1,k??????,k+??????,k??1,k?+(??????,k??????,k?)..?(??????,k??1,k?)(????,k????,k?)?1(??????,k??1,k?)?}dt,which has no interaction with that of ??????,k′ for all k′≠k since ??0 and ??0 do not enter into the evolution of the conditional covariance. Notably, the evolutions of different ??????,k with k=1,…,K can be solved in a parallel way and the computation is extremely efficient due to the small size of each individual block. This facilitates the algorithms to efficiently solve the Fokker–Planck equation in large dimensions.

Next, the structures of ??0,k and ??0,k in Eq. 9 allow the coupling among all of the K groups of variables in the conditional mean according to Eq. 5a. The evolution of ??ˉ????,k, namely the conditional mean of ??????,k conditioned on ????, is given byd??ˉ????,k(t)=[??0,k+??1,k??ˉ????,k]dt+??????,k??1,k?(????,k????,k?)?1×[d????,k?(??0,k(t,????)+??1,k??ˉ????,k)dt].

Statistical Symmetry.

The computational cost in the algorithms developed above can be further reduced if the coupled system Eq. 3 has statistical symmetry; namely,p(????,k(t),??????,k(t))=p(????,k′(t),??????,k′(t)),for allt,kandk′,including the initial conditions. The statistical symmetry is often satisfied when the underlying dynamical system represents a discrete approximation of some PDEs in a periodic domain with nonlinear advection, diffusion, and homogeneous external forcing (21, 27).

With the statistical symmetry, collecting the conditional Gaussian ensembles N(??ˉ????,k(t),??????,k(t)) for a specific k in K different simulations is equivalent to collecting that for all k with 1≤k≤K in a single simulation. This also applies to N(????i(t),??(t)) that are associated with ????. Therefore, the statistical symmetry implies that the effective sample size is L′=KL, where K is the number of the group variables that are statistically symmetric and L is the number of different simulations of the coupled systems via Monte Carlo. If K is large, then a much smaller L is needed to reach the same accuracy as in the situation without using statistical symmetry, which greatly reduces the computational cost. 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),?duidt=ui?13ui3?vi+?δ1W˙ui+du(ui+1+ui?1?2ui),i=1,…,N,[10a]dvidt=ui+a+δ2W˙vi,[10b]where ui and vi are activator and inhibitor variables and they belong to ???? and ?????? in the conditional Gaussian framework, respectively, such that the conditions of the block decomposition strategy are satisfied. Periodic boundary conditions are imposed on ui variables. In Eq. 10, the timescale ratio ?=0.01?1 leads to a slow–fast structure of the model. The parameter a=1.05>1 such that the system has a global attractor in the absence of noise and diffusion (28). The random noise is able to drive the system above the threshold level of global stability and triggers limit cycles intermittently. Note that with N=1, the model reduces to the classical FHN model with a single neuron and it contains the model families with both coherence resonance and self-induced stochastic resonance (29). With different choices of the noise strength δ1,δ2 and the diffusion coefficient du, the system in Eq. 10 exhibits rich dynamical behaviors. Below, we adopt constant parameters and the initial values are ui(0)=?2,vi(0)=0.5 for all i. Therefore, the model satisfies the statistical symmetry.

Model Behavior in Different Dynamical Regimes.

Fig. 1 shows the model behavior in three different regimes with N=500. Here, the noise coefficient δ1=0.2 is fixed while different values of δ2 and du are chosen for the three regimes. In Fig. 1A, the spatial–temporal patterns are highly coherent due to the choice of a weak noise δ2=0.1 and a strong diffusion du=10. The time series of both u1 and v1 have nearly regular oscillations with large bursts in ui around every 3 units. The associated statistical equilibrium PDF of u1 is bimodal while that of v1 is skewed. With an increase of the noise δ2=0.4 and a decrease of the diffusion du=0.5 (Fig. 1B), the coherent patterns becomes much weaker and only quasi-regular periods are found in the time series of u1. The associated PDF of v1 turns into symmetric and slightly sub-Gaussian. With a further increase of the noise to δ2=0.8 (Fig. 1C), the spatial–temporal pattern becomes strongly mixed and the time series is more irregular. Correspondingly, the PDF of v1 becomes nearly Gaussian with a large variance.

Stochastic coupled FHN model. Shown are three dynamical regimes of the stochastic coupled FHN model in Eq. 10. (A) Strongly coherent regime with δ2=0.1 and du=10. (B) Weakly coherent regime with δ2=0.4 and du=0.5. (C) Strongly mixed regime with δ2=0.8 and du=0.5. Each plot shows the spatial–temporal evolutions of u and v (A–C, Top) as well as the time series and equilibrium PDFs at i=1 (A–C, Bottom). Due to the homogeneous property, the PDFs at different grid points are the same.

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. 1B). Due to the statistical symmetry, the effective sample size is L′=NL=500L, where L is the number of repeated simulations of the systems. Below, we simply take L=1 in the efficient statistically accurate algorithms, which is extremely cheap. For comparison, we take LC=300 in Monte Carlo simulations and again use statistical symmetry to generate the true PDFs and therefore the effective sample size in Monte Carlo simulation is LC′=NLC = 150,000.

The time evolutions of the first four moments associated with u1 and v1 are shown in Fig. 2, where the black circles mark three transient phases and the statistical equilibrium phase with different non-Gaussian features. Fig. 3 shows the skill of solving the one-point statistics at these phases. First, the PDFs of u1 at t=1.6,1.9 and 25 with different bimodal features are all accurately recovered. Next, p(u1) at t=2.7 has a tiny second peak around u1=1 since only about 1.5% of the events at this time instant have large bursts and these rare events contribute to a significant kurtosis of u1 at t=2.7. Nevertheless, the efficient statistically accurate algorithm succeeds in capturing the appearance of this weak peak and the recovered PDF almost overlaps with the truth. On the other hand, the skewed p(v1) at t=1.6 and the sub-Gaussian p(v1) at t=1.9,2.7, and 25 are all recovered with high accuracy by the algorithm. In addition, the recovered joint PDFs p(u1,v1) at different phases all resemble the truth.

Stochastic coupled FHN model. Shown is time evolution of 1D marginal statistics associated with u1 and v1. Due to the homogeneous property, these statistics at different grid points remain the same. The black circles in A–D mark the time instants that the recovered PDFs from the efficient statistically accurate algorithms are compared with the truth in Figs. 3 and 4. At the statistical equilibrium phase, the skewness and kurtosis of u1 are 1.5 and 3.7, respectively.

Stochastic coupled FHN model. Shown is a comparison of the truth and recovered PDFs at three transient phases t=1.6,1.9, and 2.7 as well as the statistical equilibrium phase t=25. A and B show the truth and recovered 2D PDF p(u1,v1). C and D compare the recovered 1D PDFs p(u1) and p(v1) (blue) with the truth (red). E shows the logarithm plots of D, where the black dotted curves are the Gaussian fits of the truth. Inset above the subplot of p(u1) at t=2.7 in C is also in logarithm scale, ranging from 5×10?3 to 5.

Fig. 4 illustrates the recovered two-point statistics at t=1.6. Fig. 4 A and B shows the PDFs p(u1,u2) and p(u1,u20) and those for v, respectively. The variables u1 and u2 have a strong correlation while u1 and u20 are nearly uncorrelated. These features as well as the highly non-Gaussian PDFs with multiple peaks are both captured by the recovered PDFs from the efficient statistically accurate algorithms with only L=1. Likewise, the skewed and nearly Gaussian joint PDFs p(v1,v2) and p(v1,v20) are also recovered with high accuracy.

Stochastic coupled FHN model. Shown is a comparison of the 2D PDFs with components from two different grid points at t=1.6. Rows A and B show the truth and the recovered p(u1,u2),p(v1,v2) (A) and p(u1,u20),p(v1,v20) (B), respectively.

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 ui, which are coupled to small-scale variables vi,j,duidt=ui?1(ui+1?ui?2)+∑j=1Jγi,juivi,j?dˉiui+F+σuW˙ui,??i=1,…,I,[11a]dvi,jdt=?dvi,jvi,j?γjui2+σi,jW˙vi,j,?j=1,…,J,[11b]with periodic boundary conditions in ui. One important feature of Eq. 11 is that the nonlinear interaction between ui and vi,j conserves energy, as observed in nature (3). The two-layer L96 model belongs to the conditional Gaussian framework in Eq. 3 with ????={ui} and ??????={vi,j}. Below, we take I=40 as in the standard L96 model. Associated with each ui, there are J=5 variables vi,j representing different scales of fluctuations. Thus, the total number of state variables is 240. A constant forcing F=8 is adopted in Eq. 11a while both the damping dˉi and the coupling γi,j are functions in space. Therefore, the model is inhomogeneous. These mimic the situation that the damping and coupling above the ocean are weaker than those above the land since the latter usually have stronger friction or dissipation. As a result, the large-scale wave patterns over the ocean are more significant (Fig. 5). The model in Eq. 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.

Spatial–temporal evolution of the observed variables ui,i=1,…40 in the two-layer L-96 model Eq. 11 with F=8. The profiles of damping dˉi and coupling coefficients γi as a function of i are shown below the spatial–temporal evolution of ui in solid blue curves and the black dashed lines represent the spatial averaged values. Here, dˉi=1+0.7cos(2πi/J) and γi,j=γi=0.1+0.025cos(2πi/J).

Although the two-layer inhomogeneous L96 model in Eq. 11 has no statistical symmetry, the model structure nevertheless allows the effective block decomposition. Below, L=500 trajectories of each variable ui are simulated from Eq. 11 to implement the efficient statistically accurate algorithms. As a comparison, a direct Monte Carlo method requires LC = 150,000 samples for each of the 240 variables for an accurate estimation of at least the one-point statistics. This means the total number of samples is around 4×107! For an efficient calculation of the truth, we focus only on the statistical equilibrium state here but the algorithms are not restricted to the equilibrium statistics. The true PDFs are calculated using the Monte Carlo samples over a long time series in light of the ergodicity while the recovered PDFs from the efficient statistically accurate algorithms are computed at t=25.

Fig. 6 shows the one-point statistics as a function of i regarding the first four moments of ui and each vi,j. The mean and variance resulting from the efficient statistically accurate algorithms of all variables are highly consistent with the truth. For the skewness and kurtosis, despite small fluctuations in the recovered statistics, the pattern correlations between the curves associated with the truth and the recovered statistics are significant unless both are nearly flat (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 p(ui,vi,5) at i=11 and 21. These highly non-Gaussian PDFs with distinct features are recovered accurately by the algorithms. On the other hand, Fig. 7C shows the joint PDFs of the two smallest-scale fluctuation variables v11,5 and v13,5, which are highly correlated and have strong non-Gaussian features. Fig. 7D shows the joint PDFs p(u11,u13), where the two components also have a strong correlation. The efficient statistically accurate algorithms succeed in solving all these joint PDFs with only L=500 samples.

Two-layer inhomogeneous L-96 model. (A–D) Comparison of 2D joint PDFs at a fixed grid point (A and B) and between two different grid points (C and D).

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 ??????, using a small number of samples. This is because the sample size in these algorithms does not grow exponentially as the dimension of ??????, which is fundamentally different from Monte Carlo methods. See ref. 19 for a theoretical justification. The algorithms developed here are extremely useful in understanding the causality as well as improving the parameterizations and predictions of high-dimensional complex turbulent dynamical systems with non-Gaussian features.

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

?1To 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.

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