In non-abelian gauge theories without matter fields, expectation values of large Wilson loops and loop correlation functions are difficult to compute through numerical simulation, because the signal-to-noise ratio is very rapidly decaying for increasing loop sizes. Using a multilevel scheme that exploits the locality of the theory, we show that the statistical errors in such calculations can be exponentially reduced. We explicitly demonstrate this in the SU(3) theory, for the case of the Polyakov loop correlation function, where the efficiency of the simulation is improved by many orders of magnitude when the area bounded by the loops exceeds 1 fm 2 .