  We present new empirical constraints on the evolution of $\rho_{\rm H_2}$,
the cosmological mass density of molecular hydrogen, back to $z\approx2.5$. We
employ a statistical approach measuring the average observed $850\mu{\rm m}$
flux density of near-infrared selected galaxies as a function of redshift. The
redshift range considered corresponds to a span where the $850\mu{\rm m}$ band
probes the Rayleigh-Jeans tail of thermal dust emission in the rest-frame, and
can therefore be used as an estimate of the mass of the interstellar medium
(ISM). Our sample comprises of ${\approx}150,000$ galaxies in the UKIDSS-UDS
field with near-infrared magnitudes $K_{\rm AB}\leq25$ mag and photometric
redshifts with corresponding probability distribution functions derived from
deep 12-band photometry. With a sample approximately 2 orders of magnitude
larger than in previous works we significantly reduce statistical uncertainties
on $\rho_{\rm H_2}$ to $z\approx2.5$. Our measurements are in broad agreement
with recent direct estimates from blank field molecular gas surveys, finding
that the epoch of molecular gas coincides with the peak epoch of star formation
with $\rho_{\rm H_2}\approx2\times10^7\,{\rm M_\odot}\,{\rm Mpc^{-3}}$ at
$z\approx2$. We demonstrate that $\rho_{\rm H_2}$ can be broadly modelled by
inverting the star-formation rate density with a fixed or weakly evolving
star-formation efficiency. This 'constant efficiency' model shows a similar
evolution to our statistically derived $\rho_{\rm H_2}$, indicating that the
dominant factor driving the peak star formation history at $z\approx2$ is a
larger supply of molecular gas in galaxies rather than a significant evolution
of the star-formation rate efficiency within individual galaxies.
