We develop a framework that allows the use of the multi-level Monte Carlo (MLMC) methodology (Giles in Acta Numer. 24:259–328, 2015. https://doi.org/10.1017/S096249291500001X) to calculate expectations with respect to the invariant measure of an ergodic SDE. In that context, we study the (over-damped) Langevin equations with a strongly concave potential. We show that when appropriate contracting couplings for the numerical integrators are available, one can obtain a uniform-in-time estimate of the MLMC variance in contrast to the majority of the results in the MLMC literature. As a consequence, a root mean square error of O(ε) is achieved with O(ε−2) complexity on par with Markov Chain Monte Carlo (MCMC) methods, which, however, can be computationally intensive when applied to large datasets. Finally, we present a multi-level version of the recently introduced stochastic gradient Langevin dynamics method (Welling and Teh, in: Proceedings of the 28th ICML, 2011) built for large datasets applications. We show that this is the first stochastic gradient MCMC method with complexity O(ε−2|logε|3), in contrast to the complexity O(ε−3) of currently available methods. Numerical experiments confirm our theoretical findings.