AbstractMany-body perturbation theory methods, such as the G₀W₀ approximation, are able to accurately predict quasiparticle (QP) properties of several classes of materials. However, the calculation of the QP band structure of two-dimensional (2D) semiconductors is known to require a very dense BZ sampling, due to the sharp q-dependence of the dielectric matrix in the long-wavelength limit (q → 0). In this work, we show how the convergence of the QP corrections of 2D semiconductors with respect to the BZ sampling can be drastically improved, by combining a Monte Carlo integration with an interpolation scheme able to represent the screened potential between the calculated grid points. The method has been validated by computing the band gap of three different prototype monolayer materials: a transition metal dichalcogenide (MoS₂), a wide band gap insulator (hBN) and an anisotropic semiconductor (phosphorene). The proposed scheme shows that the convergence of the gap for these three materials up to 50meV is achieved by using k-point grids comparable to those needed by DFT calculations, while keeping the grid uniform.