  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.
