Rigorous Analysis for Efficient Statistically Accurate Algorithms for Solving Fokker-Planck Equations in Large Dimensions

This article presents a rigorous analysis for efficient statistically accurate algorithms for solving the Fokker-Planck equations associated with high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures. Despite the conditional Gaussianity, these nonlinear systems contain many strong non-Gaussian features such as intermittency and fat-tailed probability density functions (PDFs). The algorithms involve a hybrid strategy that requires only a small number of samples $L$ to capture both the transient and the equilibrium non-Gaussian PDFs with high accuracy. Here, 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. Rigorous analysis shows that the mean integrated squared error in the recovered PDFs in the high-dimensional subspace is bounded by the inverse square root of the determinant of the conditional covariance, where the conditional covariance is completely determined by the underlying dynamics and is independent of $L$. This is fundamentally different from a direct application of kernel methods to solve the full PDF, where $L$ needs to increase exponentially with the dimension of the system and the bandwidth shrinks. A detailed comparison between different methods justifies that the efficient statistically accurate algorithms are able to overcome the curse of dimensionality. It is also shown with mathematical rigour that these algorithms are robust in long time provided that the system is controllable and stochastically stable. Particularly, dynamical systems with energy-conserving quadratic nonlinearity as in many geophysical and engineering turbulence are proved to have these properties.

1. Introduction.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 [26,65].Many complex dynamical systems in geophysical and engineering turbulence, neuroscience and excitable media have large dimensions and strong nonlinearities, the associated PDFs of which are highly non-Gaussian with intermittency and extreme events [41,38].Predicting the rare and extreme events [15,19,29,63,61,20,73], quantifying the uncertainty in the presence of intermittent instabilities [47,6,30,5] and characterizing other non-Gaussian features [62,32] all require solving high-dimensional Fokker-Planck equations with strong non-Gaussian features.
Since there is no general closed-form solution for the Fokker-Planck equation, various numerical and approximate approaches have been developed to solve the evolution of the PDF p(u, t), where u consists of the state variables and t is the time.However, traditional numerical methods such as finite element and finite difference as well as the direct Monte Carlo simulations of the underlying dynamics all suffer from the curse of dimensionality [66,22,64,35,70].Furthermore, even in the low-dimensional scenarios, substantial computational cost is already required for an accurate estimation of the fat tails of the highly intermittent non-Gaussian PDFs.On the other hand, different methods for solving the partial or the approximate solutions of p(u, t) have been proposed for special dynamical systems.For example, asymptotic expansion with truncations provides good approximate PDFs associated with the slow varying variables in non-Gaussian systems with multiscale features [26,55,56,44].Splitting methods [23,24], orthogonal functions and tensor decompositions [75,71,65] are able to provide reasonably good estimations of the steady state PDFs.If the systems are weakly nonlinear with additive noise, then equivalent linearization method [69,3] is also frequently used for solving approximate solutions.
In recent work by two of the authors [14], efficient statistically accurate algorithms have been developed for solving the Fokker-Planck equation associated with high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures [11].Decomposing the state variables u into two groups u = (u I , u II ) with u I ∈ R N I and u II ∈ R N II .The conditional Gaussian systems are characterized by the fact that once a single trajectory of u I (s ≤ t) is given, u II (t) conditioned on u I (s ≤ t) becomes a Gaussian process.Despite the conditional Gaussian structure, the coupled system of u I and u II is highly nonlinear and it is able to capture many strong non-Gaussian features such as intermittency and fat-tailed PDFs that are commonly seen in nature [11].Note that in most turbulent dynamical systems, the observed variables u I represent large scale or resolved variables, which usually have only a small dimension, while the dimension of the unresolved or unobserved variables u II can be very large [53,41].Applications of the conditional Gaussian framework to highly nonlinear turbulent dynamical systems include modelling and predicting the highly intermittent and non-Gaussian times series of the Madden-Julian oscillation (MJO) and monsoon [15,10,9], filtering the stochastic skeleton model for the MJO [12], and state estimation of the turbulent ocean flows from noisy Lagrangian tracers [16,17,13].Other studies that also fit into the conditional Gaussian framework includes the dynamic stochastic superresolution of sparsely observed turbulent systems using cheap exactly solvable forecast models [7,34], stochastic superparameterization for geophysical turbulent flows [50], physics constrained nonlinear regression models [52,31], stochastic parameterized extended Kalman filter [28,27,6,8,36] and blended particle filters for high-dimensional chaotic systems [54].
The efficient statistically accurate algorithms [14] 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 u II via an extremely efficient parametric method is combined with a judicious Gaussian kernel density estimation in the low-dimensional subspace of u I .In particular, the conditional Gaussian distributions in the high-dimensional subspace are solved via closed analytical formulae and are therefore computationally efficient and accurate.The full non-Gaussian joint PDF of the system is then given by a Gaussian mixture.One remarkable feature of these efficient hybrid algorithms is that each conditional Gaussian distribution is able to cover a significant portion of the high-dimensional PDF.This guarantees the sufficiency of using only a small number of samples, which overcomes the curse of dimensionality.It has been shown in a stringent set of numerical tests [14] that with an order of O(100) samples the mixture distribution has a significant skill in capturing both the statistically steady state and the transient behavior with fat tails of the high-dimensional non-Gaussian PDFs in up to 6 dimensions while an order of O(10 6 ) samples is required in the Monte Carlo simulation to reach the same accuracy.In [14], the restriction to 6 dimension of the hybrid method is not essential but was utilized to allow comprehensive validation of the statistics in the truth model with an instructive simple model.
This article serves as a rigorous analysis for these efficient statistically accurate algorithms.The main focus here is the accuracy of the recovered PDFs in terms of the sample size L as well as its dependence on different factors, in particular the dimension of the state variables and the time span.Throughout the article, the mean integrated square error (MISE) is used to quantify the accuracy.
Our first result Theorem 3.1 reveals that the MISE in the recovered high-dimensional PDFs associated with the unresolved variables u II is bounded by E(det(R II ) −1/2 ), where R II is the conditional covariance of u II given the trajectory of u I .Notably, R II is completely determined by the underlying dynamical systems and has no dependence on the sample size L. In contrast, if a direct kernel density method is applied to recover the PDF of u II , then the bandwidth of the kernel H is scaled as the reciprocal of L to a certain power in order to minimize the MISE and the resulting MISE is proportional to L −1/N II , which means L has to increase exponentially with N II to guarantee the accuracy in the solution.This indicates the curse of dimensionality in the direct kernel density estimation and other smoothed versions of Monte Carlo methods.Such a notorious issue is overcome by the efficient statistically accurate algorithms due to the independence between R II and L in the high-dimensional subspace of u II .Another significant feature of the efficient statistically accurate algorithms is their long term persistence, which is affirmed by Theorem 3.7 in a rigorous way provided that the joint process (u I , u II ) is controllable and stochastically stable.Theorem 3.7 also supplies a lower bound of R II using the controllability condition.In addition, Proposition 3.8 demonstrates that dynamical systems with energy conserving quadratic nonlinear interactions as in most geophysical and engineering turbulence [41] automatically satisfy all the conditions for the long time persistence, which justifies the skillful performance of the efficient statistically accurate algorithms in the numerical tests reported in [14].Further validations of the controllability and other theoretical conditions in the algorithms are demonstrated in the numerical simulations at the end of this article.
The remaining of this article is organized as follows.The high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures are summarized in section 2, which is followed by a brief review of the efficient statistically accurate algorithms in [14] for solving the PDFs of such kind of systems.The main theoretical results are shown in section 3, where the proofs are included in section 4 and the appendix.In section 5, numerical tests on a nonlinear triad model and its modified versions are used to validate the theoretical results.Conclusion and discussions are given in section 6.
2. Review of the efficient statistically accurate algorithms for solving the PDFs of nonlinear dynamical systems with conditional Gaussian structures.
2.1.High-dimensional conditional Gaussian models with nonlinear and intermittent dynamical features .The general framework of high-dimensional conditional Gaussian models is given as follows [39,11]: where the state variables are u = (u I , u II ) with both u I ∈ R N I and u II ∈ R N II being multidimensional variables.In (1), A 0 , A 1 , a 0 , a 1 , Σ I and Σ II are vectors and matrices that are functions of time t and the state variables u I , and W I (t) and W II (t) are independent Wiener processes.Here the noise coefficient matrix Σ I is non-degenerated in order to guarantee the observability while there is no special requirement for Σ II .The dynamics (1) are named as conditional Gaussian systems due to the fact that once a single trajectory u I (s) for s ≤ t is given, u II (t) conditioned on u I (s) becomes a Gaussian process with mean ūII (t) and covariance R II (t), i.e., (2) p u Despite the conditional Gaussianity, the coupled system (1) remains highly nonlinear and is able to capture the strong non-Gaussian features as observed in nature [11].One of the desirable properties of the conditional Gaussian system (1) is that the conditional distribution in (2) has the following closed analytical form [39], (3) In most geophysical and engineering turbulent dynamical systems, the nonlinear terms such as the nonlinear advection have quadratic forms and these quadratic nonlinear interactions conserve energy [31,46,52,41,55,56].The nonlinear interactions allow energy transfer between different scales that induces intermittent instabilities in the turbulent dynamical systems.Such instabilities are then mitigated by energy-conserving quadratic nonlinear interactions that transfer energy back to the linearly stable modes where it is dissipated, resulting in a statistical steady state.Note that the nonlinear turbulent systems without the energyconserving nonlinear interactions may suffer from non-physical finite-time blow up of statistical solutions and pathological behavior of the related invariant measure [58].Mathematically, the turbulent dynamical systems with energy-conserving quadratic nonlinear interactions have the following abstract forms: where −Λ = L + D. Here, L is a skew-symmetric linear operator that can represent the β effect of Earth's curvature and topography, while D is a negative definite symmetric operator representing dissipative processes such as surface drag, radiative damping and viscosity, etc [67,72,45,74].The quadratic operator B(u, u) conserves energy by itself so that it satisfies the following: Notably, a rich class of turbulent models with energy-conserving quadratic nonlinear interactions in (4) belongs to the conditional Gaussian systems (1), including the noisy version of Lorenz 63 model [40], the reduced stochastic climate model [49,42], the nonlinear triad model mimicking structural features of low-frequency variability of GCMs with non-Gaussian features [48], the modified conceptual dynamical model for turbulence [53], and the two-layer Lorenz 96 model [37].See [14] and its appendix for a general framework of conditional Gaussian systems with energy-conserving nonlinear interactions as well as concrete examples.
2.2.The efficient statistically accurate algorithms for solving the PDFs of the conditional Gaussian systems.Assume the dimension N I of the observed variables is low, while the dimension N II of the unobserved variables can be high.This is the typical scenario in most turbulent dynamical systems, where the low-dimensional variables u I represent large scales or resolved variables while the high-dimensional ones u II stand for the unresolved and unobserved variables [53,41].
Below, we summarize the procedures of the efficient statistical algorithms developed in [14].First, we generate L independent trajectories from the stochastic dynamical systems (1).In fact, the only information that is required for these algorithms is L independent trajectories of the observed variables, namely u 1 I (s ≤ t), . . ., u L I (s ≤ t).Then, different strategies are used to deal with the observed variables u I and unobserved variables u II , respectively.The PDF of u II is estimated via a parametric method that exploits the closed form of the conditional Gaussian posterior statistics (3), (6) p(u Note that the limit L → ∞ in (6) (as well as (7) and (9) below) is taken to illustrate the statistical intuition, while the estimator is the non-asymptotic version.On the other hand, a Gaussian kernel density estimation method is used for solving the PDF of the observed variables u I , (7) where H = H(t) is the bandwidth matrix, and K H (•) is a Gaussian kernel centered at each sample point with covariance H(t), Below, we simply use H to represent the bandwidth at time t for the notation simplicity.
The kernel density estimation algorithm here involves a "solve-the-equation plug-in" approach for optimizing the bandwidth, the idea of which was originally proposed in [4].The solve-the-equation approach does not impose any requirement for the profile of the underlying PDF.Therefore, it works for the non-Gaussian cases and the computational cost comes from numerically solving a scalar high order algebraic equation for the optimal bandwidth in order to minimize the asymptotic mean integrated squared error (AMISE) in the estimator.Furthermore, we adopt a diagonal matrix for H.This greatly reduces the computational costs while remains the results with reasonable accuracy.Note that in the limit L → ∞, the kernel density method is simply the Monte Carlo simulation, where the bandwidth shrinks to zero.
Finally, with ( 6) and ( 7) in hand, a hybrid method is applied to solve the joint PDF of u I and u II through a Gaussian mixture, (9) p(u One important features of these algorithms is that the solutions of both the two marginal distributions in ( 6) and ( 7) and the joint distribution in (9) are consistent with those of solving the Fokker-Planck equation for p(u II (t)), p(u I (t)) and p(u I (t), u II (t)), respectively.Practically, L ∼ O(100) is sufficient for the efficient hybrid method (9) to solve the joint PDF with N I ≤ 3 and N II ∼ 10 while an order of O(10 6 ) samples is required for solving the joint PDF using classical Monte Carlo methods to reach the same accuracy for a 6 dimensional turbulent system [14].Since L is only of order O(100), the L independent trajectories u 1 I (s ≤ t), . . ., u L I (s ≤ t) can be obtained by running a Monte Carlo simulation for the coupled system (1) with L samples, which is computationally affordable.In addition, the closed form of the L conditional distributions in (6) can be computed in a parallel way due to their independence, which further reduces the computational cost.See [14] for more details.

Main theoretical results
. The rigorous analysis of the efficient statistically accurate algorithms involving the hybrid strategy (9) is studied in this section.For comparison, the theoretical results by applying the kernel density estimation method to the full system (1) is also illustrated.Note that the kernel density estimation is essentially the Monte Carlo simulation when L is large and therefore it suffers from the curse of dimensionality.Such comparison facilitates the understanding of the advantages of the efficient algorithm (9) in recovering the high-dimensional subspace of u II using only a small number of samples.Below, p t (u I , u II ) represents the true PDF while pt (u I , u II ) and pt (u I , u II ) stand for the recovered PDFs based on the pure kernel density estimation and the efficient hybrid method (9), respectively.
Kernel density estimation for the joint PDF.
Hybrid method -kernel density estimation for u I and conditional Gaussian mixture for u II .
In (11) and (13), we let H = HC as in (9).The scalar H is the scale of the bandwidth [68,76,77,4] and c 2 i are the diagonal terms of C such that c 2 i H represents the bandwidth in one direction.In the following, we mostly concern the performance of pt and pt when L is large.
One standard metric to measure the performance of a density estimator is the mean integrated squared error (MISE).The MISE of the hybrid method, for example, is the average L 2 distance to the true density: Note that pt relies on the realization of the samples and therefore it is natural to take the expectation of the distance.
Applying the Bias-Variance decomposition [25] to the MISE yields where pt := Ep t .The variance part comes from the sampling error of the method and the bias part comes from the usage of the kernel method.See (28) for a direct proof of this decomposition.
The MISE and its decomposition (14) will be used to understand the performance of the two density estimation methods in (10) and (12), where the scenarios with a large number of samples and a large dimension of the variables N II are of particular interest.Main results are presented below and the rigorous proofs of these results are shown in section 4. Note that despite quite a few studies of kernel density estimation, especially in the asymptotic limit, exist in literature [68,76,77,33,4], no analysis has been established for the hybrid method (12).Moreover, the results here are all non-asymptotic, and therefore they hold for arbitrary choice of bandwidth parameters.This is important in practice, as the bandwidth matrix H(t) may change with t.

MISE of the hybrid method.
The main result of our analysis is the following: Theorem 3.1.The two parts of MISE in (14) for the hybrid method (12) are bounded: (15) pt Here δ is any fixed strictly positive number.E is the statistical average.
) is an upper bound of the third order directional derivative of p t in the direction of u I around (u I , u II ).That is, we assume In a practical scenario, as the sample size L increases, bandwidth H can decrease, so that both the variance and bias terms decrease to zero.By taking δ close to zero and ignoring the higher order term in the bias upper bound, we recover an upper bound similar to the asymptotic MISE (AMISE) in [76,Eqn. (2.6)], except that our method also consists a random component of R II (t): where the two terms on the right hand side represents the variance and bias, respectively.It is natural to equate the order of these two terms, that is letting ).This leads to the common choice of the bandwidth [33] (18) Notably, the variance part of MISE in (17) depends on u II through E det(πR II (t)) −1 , which indicates that the hybrid method in (12) performs better with a larger R II (t).This is consistent with the intuition that a large R II (t) corresponds to a conditional distribution N (ū II (t), R II (t)) with a wide band that is able to recover a sufficient portion of the PDF.

3.2.
Comparison between the two density estimators.Theorem 3.1 already reveals the advantage of the hybrid method (12) over the the direct kernel density method (10).For a qualitative comparison of the two methods, we can view the latter as a trivial application of the hybrid method by taking u ′ I = (u I , u II ) and u ′ II = ∅, and therefore where M ≥ M is the upper bound for third order directional derivative in R N I +N II of p t .Similar results in the asymptotic setting can be found in [76].
Practically, a large L is chosen to guarantee the accuracy of the recovered PDFs, which corresponds to a small bandwidth H. Then the variance part of the direct kernel method is several magnitudes larger than that of hybrid method, especially when the dimension N II is high.
As discussed above, one would optimize the choice of H such that the two quantities in I +N II , However, This is much worse than the MISE associated with the conditional Gaussian method (18) when N II is large.Alternatively, if one wants the performance of the direct kernel method to be the same as the conditional Gaussian one (18), then the sample size needs to increase to , which can be many magnitudes larger than L.
In conclusion, direct application of the kernel method suffers from the curse of dimensionality.This is due to the fact that the variance scales with the bandwidth as , and therefore one needs to increase sample size exponentially with the dimension in order to have a small bandwidth that guarantees the accuracy of the recovered PDFs.However, when H is small, the kernel density method approximates the standard Monte Carlo simulation, which suffers from the curse of dimensionality.On the other hand, the hybrid method resolves this issue by estimating the u II part using a parametric method where the bandwidth (or the covariance) does not depend on L. Therefore, the performance of the hybrid method (12) can be much superior than the direct kernel method (10) when N II is large.

Marginal distribution of u II (t).
There are scenarios where the focus is only on estimating the density of u II (t).Again, both methods can be applied here.The direct kernel method (10) results in the estimation of the marginal density (20) pt (u On the other hand, the hybrid method ( 12) simply becomes a conditional Gaussian mixture method which contains no kernel density estimation It is straightforward to check these density estimators are the marginal PDFs of the joint distributions in (10) and (12).Since there is no kernel involved for the conditional Gaussian method in (21), the MISE has a simple bound without the bias part: Proposition 3.2.The marginal MISE of the conditional Gaussian estimator in (21) is bounded as (22) pt Following the derivation of ( 19), the MISE of the direct kernel method in ( 20) is given by With the optimal choice H ∼ O L . The hybrid method with the conditional Gaussian mixture is clearly superior for marginal density estimation, as its MISE ( 22) is essentially O(L −1 ), and the bandwidth H has no dependence on L.
3.4.Fixed subspace.In many scenarios, only a part of u II is of practical interest.To this end, we consider here u P II = Pu II , where P : R N II → R N P II maps u II onto a lower dimensional subspace.Below, we study the estimation of the density p P t (u I , u P II ) of (u I (t), u P II (t)) using the hybrid method.
It is straightforward to show the conditional distribution of u P II (t) given u I (s ≤ t) follows the Gaussian density p(u P II |u I (s ≤ t)) of the following form det(2πPR The density of (u I (t), u P II (t)) can be estimated by Following Theorem 3.1, we can show that Corollary 3.3.Under the same assumption as in Theorem 3.1, the MISE decomposition of pP t has the following two bounds where M P is a upper bound of third order derivative of p P t in u I , as in (16).Notably, the variance term depends only on E det(πPR II (t)P * ) −1 , where PR II (t)P * is a N P II × N P II matrix that is independent of the components complementary to u P II (t).In other words, the performance of the hybrid estimator on a certain part of the components is independent of the other components.This is particularly useful when N P II is small.Note that such a property also holds for the direct kernel method but in practice the kernel method works only for the case when N II is small.3.5.Controllability and a lower bound of R II .According to Theorem 3.1, R II (t) controls the sampling variance term in the MISE.Therefore, it is desirable to derive a lower bound for R II (t).Note that in the conditional Gaussian system (1), u I can be interpreted as an observation of u II , and p(u II |u I (s ≤ t)) is essentially the optimal Kalman filter with covariance R II (t).Therefore, a lower bound of R II (t) can be guaranteed by the controllability of the associated signal-observation system.In short, the controllability condition ensures the noise in the system is regular enough such that the optimal filter is not accurate to a singular degree in any component.More discussions on the controllability of Kalman filters can be found in [18,21,51].A recent work [2] has summarized some of the major results in this area.It is noteworthy that since the term a 1 depends on realization of u I , both the controllability condition and the lower bounds rely on the realization of u I .
In our context, a standard way to characterize this notion is the following assumption: Assumption 3.4.Let E s,t be the matrix flow generated by a 1 : Suppose there are constants v > 0, m ≥ 0 and D c ≥ 1 such that for any t ≥ v and s ∈ [t − v, t], Throughout this paper, for two real symmetric matrices A and B, we use A B to indicate that B − A is a positive semi-definite matrix.While A * 1 (Σ I Σ * I ) −1 A 1 actually concerns of observability, this bound is very mild.Thus, we still call Assumption 3.4 the controllability condition.
Proposition 3.5.Suppose N II ≥ 2, and the controllability condition, Assumption 3.4 holds, then for any In particular there are constants D 1 and D 2 such that The dependence of R II (t) on u I (s)| t−v≤s≤t comes from the observational term A 1 .As is seen from has a large quadratic damping, which can bring it to a very low level.
In symmetry, an upper bound can be derived if a lower bound of A * 1 (Σ I Σ * I ) −1 A 1 is assumed.Furthermore, one can show that the Riccati flow of R II (t) is contractive, so its dependence on R II (0) is diminishing.Since these results are not directly related to the performance of the hybrid estimator, we put them in the appendix along with the verification of Proposition 3.5.

3.6.
Long time performance.The simulation of (u i I (t), u i II (t)) can be maintained continuously, and the conditional Gaussian density estimator (12) can be applied for an online estimation.One important question to ask is whether the performance, and in particular the MISE, degenerates with time.If this is the case, additional samples are needed to reinforce the estimation, which is however usually difficult to carry out in practice.In this subsection, we show that the conditional Gaussian density estimator has a long time stable performance, as long as the joint process (u I , u II ) is stable and ergodic.
In stochastic analysis, the stability and ergodicity of a process can be guaranteed by energy dissipation and non-degenerate stochastic forcing.For our purpose, we can assume the energy is dissipative, while the noise is elliptic [57].
Assumption 3.6.Suppose Σ I and Σ II are full rank, and the energy is dissipative with a rate ρ > 0 and a constant D e (23) u Theorem 3.7.Under Assumption 3.6, the following hold.1) The joint density p t converges geometrically to an ergodic measure p ∞ with a rate c > 0.
In particular, there is a constant D 0 so that Here |u| 2 + 1, p 0 denotes the quantity 2) Suppose Assumption 3.4 also holds, then for any t > 0 and δ > 0, N II ≥ 2, the two parts of the MISE using the hybrid method are bounded by where D m,N II ,v is a constant independent of L and H, and M ∞ is a bound for the third order u I -directional derivative of p ∞ as in (16).
In particular, when t → ∞, we have lim sup This leads to the same bandwidth and MISE scaling with L, namely: .
The proof strategy of Theorem 3.7 is straightforward.The first part is simply corollaries of [60,59,1].To reach a bound on the variance part in 2), it suffices to have a lower bound on This can be achieved by Proposition 3.5 and an energy dissipation argument.
For the bias term, we use the Poincaré inequality (24) to approximate it with the bias term at equilibrium.

3.7.
Conditional Gaussian turbulent dynamical systems with energy-conserving quadratic nonlinearity.Recall the turbulence model u with quadratic energy conserving nonlinear interactions (4)-( 5) The linear damping part provides a uniform dissipation, so for some λ − > 0, and the nonlinearity term B is quadratic and conserves energy.
In our conditional Gaussian setup, we can decompose the dynamics into the form below ( 25) The quantities in the brackets naturally correspond to A 0 , A 1 , a 0 and a 1 respectively.For the damping term Λ, we assume there are constants 0 The energy conservation condition, u • B(u, u) = 0, requires that See the Appendix of [14] for details.
The true density can be written as p t (u , since for any test function f , the following holds This gives the following result where * denotes the convolution.The Variance-Bias decomposition of the MISE can be made: In Lemma A.2, a Taylor expansion on (p t (u ′ I , u II ) − p t (u I , u II )) leads to the following upper bound for the bias part: Moreover, in light of the relation pt (u I , u II ) = 1 L L i=1 pi (u I , u II ) and the independence of the density samples pi , we have var pt (u I , u II )du Note that each pi (x, y) is a Gaussian density with mean (u i I (t), ūII (t)) and a block diagonal covariance, where the blocks are given by HC and R II (t), respectively.In Lemma A.1, a straightforward computation of the L 2 norm of a Gaussian density shows that .
This leads to the bound of the MISE.Proof of Proposition 3.2.Denote pi (u II ) = p(u II |u i I (s ≤ t)), then following the same proof as in Theorem 3.1, we have p t (u II ) = Ep i (u II ) and pt (u II ) = 1 L L i=1 pi (u II ).Thus, .
Proof of Corollary 3.3.The proof is identical to the one of Theorem 3.1, as long as one replaces the densities involving u II to the version for u P II .Therefore it is omitted here.

Long time result.
Proof of Theorem 3.7.Part 1): The geometric ergodicity, i.e. the following L 1 convergence, is a direct result that comes from the framework of [60,59].Its equivalence to the Poincaré type of inequality ( 24) is a result by [1].We will try to verify the conditions needed in [1].
We claim that V (u I , u II ) = |u I | 2 + |u II | 2 + 1 is a Lyapunov function of Definition 1.1 in [1].Apply the generator L of the diffusion process where b = 2ρ + 2D e + tr(Σ I Σ * I + Σ II Σ * II ), and U = {V (u I , u II ) ≤ b}.The fact that U , and actually any compact subset, is a petite set can be verified by the same proof of Lemma 3.4 in [59], since we assume Σ I and Σ II are full rank.The fact the stochastic process is irreducible can also be verified using the same argument.More details on these arguments are provided in [57] for more general conditions.
Therefore, applying theorem 1.2 of [1] leads to the L 1 convergence above.Theorem 2.1 also applies with f (u I , u II ) = p 0 p∞ (u I , u II ), which gives (24).
Following the proof of Theorem 3.1, we have the variance part To provide a bound for E|u I (t)| mN II , we verify that any fixed moment |u| 2n = (|u I | 2 + |u II | 2 ) n is also dissipative.Applying the generator of the diffusion process yields where Σ = [Σ * I , Σ * II ] * and the constant D n,Σ exists because of Young's inequality.Apply Dynkin's formula for e ρnt |u(t)| 2n , and combine it with the result above, we have the following Gronswall's inequality ( 29) To continue, we let n = mN II /2 in (29) and integrate it in time range [t − v, t], Consequently, there exists a constant D m,N II ,v such that For the bias term, |p t (u I , u II ) − p t (u I , u II )| 2 du I du II , we use the Cauchy Schwartz Recall that p∞ = K H * p ∞ .Using the same proof as in Theorem 3.1, we have Then apply (24), we have Next, recall that pt (u Consequently, Combining the results finishes the proof.

Numerical examples. Below
, numerical examples are used to support the theoretical results in section 3. The test model considered here is the following triad model [52], where A 1 + A 2 + A 3 = 0 represents the energy-conserving nonlinear interactions and d 2 > 0, d 3 > 0 are the damping terms.Note that there is no damping and dissipation in (30a) but ( 30) is a hypoelliptic diffusion [57,59].Linear stability is satisfied for u 2 , u 3 while there is only neutral stability of u 1 .Define E 2 = σ 2 2 /(2d 2 ) and E 3 = σ 2 3 /(2d 3 ).It is straightforward to show that the triad system (30) has a Gaussian invariant measure [43,52] p eq (u) = C exp − 1 2 provided that the following condition is satisfied If the condition in ( 31) is violated, namely E 1 < 0, then the variance in u 1 direction will increase unboundedly and there is no invariant measure for the triad system (30).Below, two dynamical regimes of the triad model ( 30) are studied, where the corresponding parameters are listed in the Table 1.Particularly, the triad system (30) in Regime I has a Gaussian invariant measure while there is no invariant measure in Regime II due to the fact that E 1 < 0. See Figure 1 for the time evolution of the three marginal variances and one realization of each variable and [41] for dynamical introduction about such triad models.

Table 1
Parameters of two dynamical regimes of the triad model (30) Denote u I = (u 2 , u 3 ) T and u II = u 1 .The triad system (30) belongs to the conditional Gaussian family (1).Notably, the noise coefficient in u II is Σ II = 0, which implies the system has no controllability.The initial values in the tests below are all given at origin.Here only the hybrid method ( 9) is tested and the number of samples is always L = 500.
Figure 2 shows the recovered PDF at t = 1 in Regime I of the triad model.Despite an accurate estimation of the joint PDF of the observed variables p(u 2 , u 3 ) as shown in Panel (e), the recovered PDF of the unobserved variable u 1 in Panel (f) has quite a few noisy fluctuations and the recovered joint PDFs p(u 1 , u 2 ) and p(u 3 , u 1 ) in Panel (d) and (f) are non-smooth in u 1 direction as well.Such pathological behavior results from the loss of controllability of the system, which is consistent with the theoretical discussion in subsection 3.5.In fact, the term a 1 in (1) associated with the triad system (30) is zero.Therefore, according to (3), Σ II = 0 implies the posterior variance R II = 0 and the posterior mean ūII simply follows the sampled trajectory of u II .In other words, the posterior states from the algorithm are exactly the Monte Carlo samples, as is validated in Panel (h).The same performance is found in Regime II and thus we omit the figure here.
In order to make the triad system have controllability, a small noise is added to (30a) and the resulting modified triad system is given as follows, where ǫ is the noise coefficient of u 1 with Σ II = ǫ in (1).Below we set ǫ = 0.1 ≪ σ 2 = σ 3 = 1.The other parameters in (32) remain the same as those in Table 1.
This extra noise implies the triad system is controllable, which significantly improves the accuracy of the recovered PDFs.See Figure 3 for the results in Regime I at t = 1.In particular, Panel (h) of Figure 3 shows that the posterior means are quite different from the Monte Carlo samples and the posterior variances are no longer zero.It is also shown in Figure 6 that the recovered PDFs at a long time t = 20 (i.e., statistically steady state) are very close to the truth with this extra small noise.
Similarly, Figure 4 shows the recovered PDFs of Regime II with ǫ = 0.1 at t = 1, the error in which compared with the truth is negligible.Notably, although the amplitude of u 1 has an unbounded growth in this regime due to the fact that E 1 < 0, the recovered PDFs with ǫ = 0.1 at t = 20 as illustrated in Figure 7 remain quite accurate.Next, the performance of the hybrid algorithm at a very long time in this regime is studied.Figure 5 shows the recovered PDFs at t = 400.Similar to Figure 2, the noisy fluctuations are found in the recovered PDF of u 1 .In fact, direct calculations show that the posterior variance R II in (3) is bounded from above since the unbounded signal u 1 does not enter into the evolution of R II , which is also validated by the numerical simulation in Panel (h).Since the variance of u 1 increases with time, the percentage of the portion covered by each conditional Gaussian distribution decreases in time, which reduces the skill in the recovered PDFs by the conditional Gaussian mixtures.In Figure 8, we show that by further imposing a damping in the dynamics of u 1 of the modified triad model (32) in Regime II, the model then satisfies all the conditions in Proposition 3.8 and the resulting model has an invariant measure.In such a scenario, the hybrid algorithm is skillful in both short and long time as is affirmed by Proposition 3.8.
It is also worthwhile pointing out that all the test models in [14], including the noisy version of Lorenz 63 model [40], the stochastic climate model [49,42], the nonlinear triad model mimicking structural features of low-frequency variability of GCMs with non-Gaussian features [48] and the modified conceptual dynamical model for turbulence [53], all satisfy the conditions in Proposition 3.8.Therefore, the hybrid algorithm ( 9) is able to solve the PDFs of those models with high accuracy with only a small number of samples.6. Discussion and Conclusions.This article presents a rigorous analysis for the efficient statistically accurate algorithms developed in [14], which succeed in solving both the transient and the equilibrium solutions of Fokker-Planck equations associated with high-dimensional nonlinear turbulent dynamical systems with conditional Gaussian structures.Despite the conditional Gaussianity, these nonlinear systems capture many strong non-Gaussian features such as intermittency and fat-tailed PDFs.The algorithms involve a hybrid strategy that requires only a small number of samples L to capture both the transient and the equilibrium non-Gaussian PDFs with high accuracy.
Theorem 3.1 shows that the MISE in the recovered high-dimensional PDFs associated with the unresolved variables u II is bounded by E(det(R II ) −1/2 ), where R II is completely determined by the underlying dynamical systems and it has no dependence on the sample size L.This is fundamentally different from the direct application of the kernel methods to recover the PDF of u II , in which the bandwidth of the kernel H is scaled as a reciprocal of L to a certain power and the resulting MISE is proportional to L −1/N II .This implies the curse of dimensionality in the kernel density estimation and other smoothed Monte Carlo methods due to the fact that L has to increase exponentially as N II in order to guarantee the accuracy in the solution.As is shown in Theorem 3.1, many fewer samples are needed in the efficient statistically accurate algorithms in order to reach the same accuracy as using the smoothed Monte Carlo methods, especially with a large N II .Theorem 3.7 affirms the long term persistence of the efficient statistically accurate algorithms in a rigorous way under the assumption that the joint process (u I , u II ) is controllable and stochastically stable.It also provides a lower bound of R II using the controllability condition.The validations of the controllability and other theoretical conditions in the algorithms are demonstrated in the numerical simulations in section 5. Furthermore, Proposition 3.8 illustrates that the turbulent dynamical systems with quadratic energy conserving nonlinear interactions [41] automatically satisfy all the conditions for the long time persistence.This justifies the skillful performance of the efficient statistically accurate algorithms in the numerical tests reported in [14] and provides important guidelines for future applications.
Appendix.This appendix contains the following information.Section A states and proves two lemmas that are needed to prove Theorem 3. Proof.Note the Gaussian density has form Its square can be decomposed as The second part is the Gaussian density N (a, 2Σ), so its integral is one.This concludes our proof.Lemma A.2. Suppose p t (u I , u II ) has bounded third derivative as (16).Consider filtering it with kernel K H at the u I components, define Then |p t (u I , u II ) − p t (u I , u II )| 2 du I du II is bounded by Proof.Apply Taylor's expansion to p t (u I , u II ), use (16) Denote the gradient and Hessian with respect to u I as ∇ I and ∇ 2 I .Note that and Therefore, To continue, we write where the range of s is [0, 1].Note that by the symmetry of K H (u I ) in u I , Note that K H (u ′ I − u I ) is the density of N (u I , HC), where C is diagonal with diagonal entries c 2 i .If we let Z be a random variable following this distribution Therefore Finally, by Cauchy Schwartz, Note again K H is the density of Z ∼ N (0, HC), Therefore Combine these estimates with In [2], the matrix flow generated by this equation was studied.(In [2], it was termed as φ(t)(Q)).Define the controllability and observability Gramian [2]: Define also the following processes: Theorem 4.4 in [2] has shown that by applying the comparison principal between the Riccati equation of R II (t) and the bounds in (33).
For the purpose of this proposition, we focus only on the left hand side of (33).First we find bounds for the controllability Gramian.Note that Then for t − v ≤ s ≤ t, apply Assumption 3.4 and the bounds above, we have Consequentially, by taking s = t − v, we have R II (t) −1 h t,v (u I )I N II , where In particular, this leads to detR II (t) .
Then applying Holder's inequality yields Combining the estimates above finishes the proof of Proposition 3.5.
In symmetry to Proposition 3.5, we can also find an upper bound for R II (t).For that purpose, let σ 2 A,− (t) ≥ 0 be a stochastic process that satisfies: Proposition B.1.Under Assumption 3.4, for t ≥ v, we have R II (t) g v,t (u I ), where Note that σ 2 A,− (t) essentially characterizes how well the observation of u II is made at time t.It can be zero or close to zero in certain time periods.During such time periods, little information is provided by the upper bound in Proposition B.1.
Proof of Proposition B.1.We continue our discussion from the end of the proof of Proposition 3.5.For t − v ≤ s ≤ t, the observability Gramian is bounded by Likewise, we have In conclusion, by (33) we have R II (t) g v,t (u I )I N II , where Appendix C. Contraction of the Riccati flow.One practical issue in applying the hybrid method is how to initialize R II (0), as p 0 (u I , u II ) often does not have an explicit form.In fact, the value of pt has a diminishing dependence on R II (0) as t → ∞.One way to characterize this is to consider a covariance process R ′ II (t) that follows Here a 1 and A 1 are the same as the one in (3), but R ′ II (0) = R II (0).In fact, the rescaled difference between R II (t) and R II (t) ′ converges geometrically fast, although the contraction rate depends on the realization of u I : Proposition C.1.Under the same assumptions of Proposition B.1, for any t where • denotes the spectral radius, and Note that the same lower and upper bounds in Proposition 3.5 and Proposition B.1 apply to R II (t) and R ′ II (t) here, as they are both driven by u I , which is the only thing the bounds depend on.
Remark C.2.When σ A,− is bounded uniformly away from zero and A 1 has no dependence on u I , namely m = 1 in Assumption 3.4, then f v,t is bounded away from zero uniformly in time.In this case, one can also show the initialization of ūII (0) has diminishing influence on pt , since a ū′ II (t) process driven by (3) will converge to ūII (t) even if ū′ II (0) = ūII (0) [2].Although in principle the convergence of the mean process should hold in more general scenarios, the convergence rate will have a very involved dependence on the realization of u I (s).Also it is not a significant property for our estimator, and therefore we omitted here.
Proof of Proposition C.1.We continue our discussion from the end of the proof of Proposition B.1.In Proposition 3.3 of [2], it is shown that Then in the third displayed equation of the proof of theorem 4.8 in [2], we have the identity By Proposition 3.5, we have Consequently, By Gronwall's inequality, we have and therefore, The same inequality also holds for E s,t (R ′ II ) .In conclusion, we have   Appendix E. Numerical simulations of the modified triad models.Figure 6 and Figure 7 show the recovered PDFs of the modified triad model (32) at t = 20 for Regime I and II, respectively.With the extra small noise ǫ in the dynamics of u 1 , the controllability is regained and the skill of the hybrid algorithm is greatly improved.
Figure 8 shows the statistics and the recovery of the PDFs of the following system  where an extra damping −d 1 u 1 is imposed in addition to the small noise term in the u 1 dynamics.The parameter d 1 = 0.1 and other parameters are the same as those in (32).The model (34) satisfies all the conditions in Proposition 3.8 and therefore hybrid algorithm is always skillful.Note that the triad model ( 34) is linearly stable with respect to all the three variables and no unbounded growth in the u 1 direction as the versions in (32) and (34).

( 19 ) 2 4+NI 4 4+N
are of the same order, which leads to the scaling H ∼ O L − +N II , and also the overall MISE ∼ O L −

Proposition 3 . 8 .
For the stochastic flow with energy conserving quadratic nonlinearity(25), assume that(26) and (27) hold, and Σ I and Σ II are of full rank.We have the following results: 1).Assumption 3.6 holds with ρ = 1 2 λ − and D e = 1 2λ − (|F I | 2 + |F II | 2 ).2).Assumption 3.4 holds with v = 1, m = 1 and where the constants are chosen such that |B II,1 (u I )| ≤ λ B |u I | and σ 2 I,− I N I Σ I Σ * I , σ 2 II,− I N II Σ II Σ * II σ 2 II,+ I N II .The proof of Proposition 3.8 is shown in D. The energy conservation property plays an essential role in verifying the system stability, and 4. Proofs.4.1.Finite time MISE.Proof of Theorem 3.1.Denote the one sample path density function:

Figure 1 .
Figure 1.Triad model (30).(a) Marginal variance as a function of time (t ∈ [0, 100]) in the two dynamical regimes with parameters in Table 1.(b) Sample trajectories up to t = 1000 of the two dynamical regimes.Note the unbounded growth of the amplitude of u1 in Regime II.

1 .
Section B shows the proof of Proposi-tion 3.5 regarding the controllability and observability and Section C includes the discussions of the contraction of the Riccati flow.The controllability and long time behavior of the conditional Gaussian turbulent dynamical systems with energy-conserving quadratic nonlinearity in Proposition 3.8 are demonstrated in Section D. Finally, extra numerical examples of the triad model and modified version of the triad model are shown in E. Appendix A. Two lemmas for Theorem 3.1.Lemma A.1.Let p(u II ) be the PDF of N (a, Σ), its L 2 norm is p 2 (u II )du II = 1 det(πΣ) .

Figure 8 .
Figure 8. Modified triad model (34) with extra damping d1u1, Regime II at t = 50.Panels (a)-(c): time evolution of the variance of the three components.Panels (d)-(k): same captions as in Figure 2.