  Gaussian processes are the gold standard for many real-world modeling
problems, especially in cases where a model's success hinges upon its ability
to faithfully represent predictive uncertainty. These problems typically exist
as parts of larger frameworks, wherein quantities of interest are ultimately
defined by integrating over posterior distributions. These quantities are
frequently intractable, motivating the use of Monte Carlo methods. Despite
substantial progress in scaling up Gaussian processes to large training sets,
methods for accurately generating draws from their posterior distributions
still scale cubically in the number of test locations. We identify a
decomposition of Gaussian processes that naturally lends itself to scalable
sampling by separating out the prior from the data. Building off of this
factorization, we propose an easy-to-use and general-purpose approach for fast
posterior sampling, which seamlessly pairs with sparse approximations to afford
scalability both during training and at test time. In a series of experiments
designed to test competing sampling schemes' statistical properties and
practical ramifications, we demonstrate how decoupled sample paths accurately
represent Gaussian process posteriors at a fraction of the usual cost.
