  When comparing approximate Gaussian process (GP) models, it can be helpful to
be able to generate data from any GP. If we are interested in how approximate
methods perform at scale, we may wish to generate very large synthetic datasets
to evaluate them. Na\"{i}vely doing so would cost \(\mathcal{O}(n^3)\) flops
and \(\mathcal{O}(n^2)\) memory to generate a size \(n\) sample. We demonstrate
how to scale such data generation to large \(n\) whilst still providing
guarantees that, with high probability, the sample is indistinguishable from a
sample from the desired GP.
