  This paper considers the strong error analysis of the Euler and fast Euler
methods for nonlinear overdamped generalized Langevin equations driven by the
fractional noise. The main difficulty lies in handling the interaction between
the fractional Brownian motion and the singular kernel, which is overcome by
means of the Malliavin calculus and fine estimates of several multiple singular
integrals. Consequently, these two methods are proved to be strongly convergent
with order nearly $\min\{2(H+\alpha-1), \alpha\}$, where $H \in (1/2,1)$ and
$\alpha\in(1-H,1)$ respectively characterize the singularity levels of
fractional noises and singular kernels in the underlying equation. This result
improves the existing convergence order $H+\alpha-1$ of Euler methods for the
nonlinear case, and gives a positive answer to the open problem raised in [4].
As an application of the theoretical findings, we further investigate the
complexity of the multilevel Monte Carlo simulation based on the fast Euler
method, which turns out to behave better performance than the standard Monte
Carlo simulation when computing the expectation of functionals of the
considered equation.
