  Gaussian processes are the ideal tool for modelling the Galactic ISM,
combining statistical flexibility with a good match to the underlying physics.
In an earlier paper we outlined how they can be employed to construct
three-dimensional maps of dust extinction from stellar surveys. Gaussian
processes scale poorly to large datasets though, which put the analysis of
realistic catalogues out of reach. Here we show how a novel combination of the
Expectation Propagation method and certain sparse matrix approximations can be
used to accelerate the dust mapping problem. We demonstrate, using simulated
Gaia data, that the resultant algorithm is fast, accurate and precise.
Critically, it can be scaled up to map the Gaia catalogue.
