  We propose a novel discrete Poisson equation approach to estimate the
statistical error of a broad class of numerical integrators for the underdamped
Langevin dynamics. The statistical error refers to the mean square error of the
estimator to the exact ensemble average with a finite number of iterations.
With the proposed error analysis framework, we show that when the potential
function $U(x)$ is strongly convex in $\mathbb R^d$ and the numerical
integrator has strong order $p$, the statistical error is
$O(h^{2p}+\frac1{Nh})$, where $h$ is the time step and $N$ is the number of
iterations. Besides, this approach can be adopted to analyze integrators with
stochastic gradients, and quantitative estimates can be derived as well. Our
approach only requires the geometric ergodicity of the continuous-time
underdamped Langevin dynamics, and relaxes the constraint on the time step.
