Summary: 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 [D. Fang and L. Li, ESAIM, Math. Model. Numer. Anal. 54, No. 2, 431–463 (2020; Zbl 07201590)]. 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. Finally, numerical experiments are carried out to support the theoretical results.