Existing distribution compression methods, like Kernel Herding (KH), were originally developed for unlabelled data. However, no existing approach directly compresses the conditional distribution of labelled data. To address this gap, we first introduce the Average Maximum Conditional Mean Discrepancy (AMCMD), a natural metric for comparing conditional distributions. We then derive a consistent estimator for the AMCMD and establish its rate of convergence. Next, we make a key observation: in the context of distribution compression, the cost of constructing a compressed set targeting the AMCMD can be reduced from $\mathcal{O}(n^3)$ to $\mathcal{O}(n)$. Building on this, we extend the idea of KH to develop Average Conditional Kernel Herding (ACKH), a linear-time greedy algorithm that constructs a compressed set targeting the AMCMD. To better understand the advantages of directly compressing the conditional distribution rather than doing so via the joint distribution, we introduce Joint Kernel Herding (JKH), a straightforward adaptation of KH designed to compress the joint distribution of labelled data. While herding methods provide a simple and interpretable selection process, they rely on a greedy heuristic. To explore alternative optimisation strategies, we propose Joint Kernel Inducing Points (JKIP) and Average Conditional Kernel Inducing Points (ACKIP), which jointly optimise the compressed set while maintaining linear complexity. Experiments show that directly preserving conditional distributions with ACKIP outperforms both joint distribution compression (via JKH and JKIP) and the greedy selection used in ACKH. Moreover, we see that JKIP consistently outperforms JKH.

We show that common choices of kernel functions for a highly accurate and massively scalable nearest-neighbour based GP regression model (GPnn: \cite{GPnn}) exhibit gradual convergence to asymptotic behaviour as dataset-size $n$ increases. For isotropic kernels such as Mat\'{e}rn and squared-exponential, an upper bound on the predictive MSE can be obtained as $O(n^{-\frac{p}{d}})$ for input dimension $d$, $p$ dictated by the kernel (and $d>p$) and fixed number of nearest-neighbours $m$ with minimal assumptions on the input distribution. Similar bounds can be found under model misspecification and combined to give overall rates of convergence of both MSE and an important calibration metric. We show that lower bounds on $n$ can be given in terms of $m$, $l$, $p$, $d$, a tolerance $\varepsilon$ and a probability $\delta$. When $m$ is chosen to be $O(n^{\frac{p}{p+d}})$ minimax optimal rates of convergence are attained. Finally, we demonstrate empirical performance and show that in many cases convergence occurs faster than the upper bounds given here.

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.