0

# Self-Distillation for Gaussian Process Regression and Classification

Kenneth Borup
Department of Mathematics
Aarhus University
<EMAIL_ADDRESS>
&Lars Nørvang Andersen
Department of Mathematics
Aarhus University
<EMAIL_ADDRESS>
###### Abstract

We propose two approaches to extend the notion of knowledge distillation to
Gaussian Process Regression (GPR) and Gaussian Process Classification (GPC);
data-centric and distribution-centric. The data-centric approach resembles
most current distillation techniques for machine learning, and refits a model
on deterministic predictions from the teacher, while the distribution-centric
approach, re-uses the full probabilistic posterior for the next iteration. By
analyzing the properties of these approaches, we show that the data-centric
approach for GPR closely relates to known results for self-distillation of
kernel ridge regression and that the distribution-centric approach for GPR
corresponds to ordinary GPR with a very particular choice of hyperparameters.
Furthermore, we demonstrate that the distribution-centric approach for GPC
approximately corresponds to data duplication and a particular scaling of the
covariance and that the data-centric approach for GPC requires redefining the
model from a Binomial likelihood to a continuous Bernoulli likelihood to be
well-specified. To the best of our knowledge, our proposed approaches are the
first to formulate knowledge distillation specifically for Gaussian Process
models.

Code: github.com/Kennethborup/gaussian_process_self_distillation

##  1 Introduction and Problem Setting

In this paper, we propose a notion of knowledge distillation (KD) in a
Gaussian process (GP) context. Despite numerous research and applications of
(extensions of) KD in the context of deep learning, research on knowledge
distillation for other methods of machine learning than deep learning is
lacking. Recent work by Borup & Andersen (2021); Mobahi et al. (2020); Phuong
& Lampert (2019); Frosst & Hinton (2018) amongst others, provide results on KD
for simplified settings, such as for (kernel ridge) regression. A natural
extension to these results is to investigate KD for GPs, but since no prior
work has been performed in this context, we initially propose multiple
approaches to how one could define knowledge distillation for GPs. In this
paper, we restrict our analysis to self-distillation, i.e. where the student
and teacher models are of an identical class, which in practice corresponds to
the kernel function being identical for both models.

Unlike e.g. deep learning, where samples and models are deterministic during
training, GPs include stochasticity in the form of a prior and thus posterior
predictive distribution. In ordinary GP regression, one defines a prior, and
based on a set of observations $\mathcal{D}$, a posterior model is obtained
which can be used for inference on previously unseen samples. We propose two
different and straightforward approaches to self-distillation of GPs: _data-
centric_ and _distribution-centric_ (see Figure 2).

Figure 1: Conceptual illustration of our two proposed methods of performing
distillation for Gaussian processes. In data-centric self-distillation, we re-
use the mean function of the previous solution in the next step, and in
distribution-centric self-distillation, we re-use the full distribution.

#### Data-centric self-distillation.

In data-centric distillation, we treat the output predictions of a model at
step $t$ as the input to the model at step $t+1$. This procedure aligns well
with most current distillation methods known from deep learning (Hinton et
al., 2015). However, unlike in deep learning, the output of our models is not
merely scalars, but distributions. We propose to discard the stochasticity of
our predictions at step $t$ and merely keep the mean predictions for training
of the model at step $t+1$. We show that for GP regression the mean after
self-distillation is closely related to distillation with kernel ridge
regression as in Borup & Andersen (2021); Mobahi et al. (2020). Furthermore,
we argue that for distillation in a classification setting, one needs to adapt
the likelihood function to support continuous targets rather than binary
targets. We address this by utilization of the continuous Bernoulli
distribution (Loaiza-Ganem & Cunningham, 2019). We investigate data-centric
distillation for regression in Section 3.1 and for classification in Section
4.1.

#### Distribution-centric self-distillation.

In distribution-centric distillation, we replace the prior at a given step,
$t$, with the posterior of the previous step, $t-1$. In GP regression the
posterior distribution is a GP itself, and we can directly consider it as a
prior for the succeeding step. This yields an iteratively refined and data-
dependent prior with each distillation step. For classification, more care has
to be taken as the posterior is not known analytically. However, by utilizing
the Laplace approximation of the posterior as the prior of the latent GP
model, we are able to define a meaningful notion of distribution-centric
distillation in the classification setup as well. We investigate distribution-
centric for regression in Section 3.2 and for classification in Section 4.2.

#### Combined and extended directions.

As is commonly practiced in the literature, one can consider the extension of
the data-centric distillation approach, where we not only consider the
posterior predictions of the previous step but rather consider a convex
combination (weighted by $\alpha\in(0,1)$) of these posterior predictions and
the original observations. Furthermore, since the data-centric and
distribution-centric approaches are not mutually exclusive, we could consider
a combined method, with an inner and outer distillation loop. In the outer
loop the prior is replaced by the posterior of the preceding step as in
distribution-centric distillation, but in the inner loop, the posterior is
distilled a set number of steps using data-centric distillation, with the
prior fixed in the outer loop. We leave such intricate, but intriguing setups
for future research.

We summarize our contributions as:

  * •

We propose two different approaches to self-distillation for both GP
regression and GP classification, reusing either the mean predictions (data-
centric) or full predictive distribution (distribution-centric) of the
previous iteration of the model.

  * •

We show a connection between data-centric self-distillation for GP regression
and self-distillation of kernel ridge regression. Furthermore, we prove that
distribution-centric self-distillation for GP regression is equivalent to
ordinary GP regression with a particular choice of hyperparameter.

  * •

We find that a naïve approach to data-centric self-distillation for GP
classification yields a misspecified model due to the continuous form of the
predictions, and we propose to alleviate this by utilization of the continuous
Bernoulli distribution. Furthermore, we find that distribution-centric self-
distillation can be efficiently approximated by a particular scaling of the
covariance function for any number of distillation steps.

The structure of the paper is as follows. In Section 2 we present preliminary
results on GP regression and GP classification. Next, we focus on self-
distillation for GP regression in Section 3, and in particular on data-centric
distillation in Section 3.1 and distribution-centric distillation in Section
3.2. Afterward, in Section 4 we investigate the classification setup and in
particular the data-centric setting in Section 4.1 and distribution-centric
setting in Section 4.2. Finally, we relate our work to existing research in
Section 5. Additionally, in the appendix we provide additional mathematical
details, implementation details, experiment ablations, as well as all proofs.

(a) Classical.

(b) Distribution-centric.

(c) Data-centric.

Figure 2: In (2(a)) we illustrate a typical procedure to fit a GP to an
observed dataset, $\mathcal{D}=\\{(x_{i},y_{i})\\}_{i=0}^{n}$, and compute the
output distribution for some unseen sample $x_{\ast}$. In (2(b)) we illustrate
the distribution-centric self-distillation procedure, where we reuse the
posterior as a prior. In (2(c)), we illustrate the data-centric self-
distillation procedure, where we iteratively replace the training observations
$\mathcal{D}_{t}=\\{(x_{i},y_{t,i})\\}_{i=0}^{n}$ with the posterior
predictions of the current model,
$\mathcal{D}_{t+1}=\\{(x_{i},y_{t+1,i})\\}_{i=0}^{n}$, and refit the model to
these samples.

##  2 Preliminaries

Our starting point is the standard set-up for Gaussian Processes regression.
These results are well-known (see e.g. Rasmussen & Williams (2006)).

###  2.1 Gaussian Process regression

We assume a generic input-output pair
$(\mathbf{x},y)\in\mathbb{R}^{d}\times\mathbb{R}$ is related by
$y=f(\mathbf{x})+\varepsilon$, where $f(\cdot)$ is a function, which we wish
to infer, based on the assumption that $f$ is a random function whose prior
distribution is $f\sim\mathcal{GP}(m,k)$ for mean function
$m:\mathbb{R}^{d}\to\mathbb{R}$ and positive semi-definite covariance function
(often called kernel) $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$, and
$\varepsilon$ is an independent noise term with
$\varepsilon\sim\mathcal{N}(0,\gamma)$. For notational convenience we will
assume that our input is univariate, i.e. $d=1$. For a training set
$\mathcal{D}=\\{(x_{i},y_{i})\mid i=1,\dots,N\\}$ and test set $(x_{\ast
1},\dots,x_{\ast M})^{\intercal}$, the first task is to describe the posterior
distribution of $\mathbf{f}_{\ast}=(f(x_{\ast 1}),\dots,f(x_{\ast
M}))^{\intercal}$ given $\mathcal{D}$. The posterior predictive distribution
of $\mathbf{f}_{\ast}$ given $\mathcal{D}$ is

$\displaystyle\mathbf{f}_{\ast}$
$\displaystyle\mid\mathbf{x},\mathbf{y},\mathbf{x}_{\ast}\sim\mathcal{N}(\bm{\mu}_{\ast},\mathbf{\Sigma}_{\ast}),\quad\text{where}$
$\displaystyle\bm{\mu}_{\ast}$
$\displaystyle=m(\mathbf{x}_{\ast})+k(\mathbf{x}_{\ast},\mathbf{x}^{\intercal})(\mathbf{K}+\gamma\mathbf{I}_{N})^{-1}\left(\mathbf{y}-m(\mathbf{x})\right),$
$\displaystyle\mathbf{\Sigma}_{\ast}$
$\displaystyle=k(\mathbf{x}_{\ast},\mathbf{x}_{\ast}^{\intercal})-k(\mathbf{x}_{\ast},\mathbf{x}^{\intercal})(\mathbf{K}+\gamma\mathbf{I}_{N})^{-1}k(\mathbf{x},\mathbf{x}_{\ast}^{\intercal}).$

and where $\mathbf{x}$ and $\mathbf{x}_{\ast}$ are stacked training and test
input, respectively. In the expressions for $\bm{\mu}_{\ast}$ and
$\mathbf{\Sigma}_{\ast}$, we have used the notation

$\displaystyle g(\mathbf{x},\mathbf{y}^{\intercal})$
$\displaystyle=\left\\{g(x_{i},y_{j})\right\\}_{i,j=1}^{N,M}\in\mathbb{R}^{N\times
M}$

where $g:\mathbb{R}^{2}\ni(x,y)\mapsto g(x,y)\in\mathbb{R}$, and
$\mathbf{x}=(x_{1},\dots,x_{N})^{\intercal}\in\mathbb{R}^{N}$,
$\mathbf{y}=(y_{1},\dots,y_{M})^{\intercal}\in\mathbb{R}^{M}$ are column
vectors. Notice, that if either argument is a scalar we have

$\displaystyle g(x,\mathbf{x}^{\intercal})$
$\displaystyle=\left(g(x,x_{1}),\dots,g(x,x_{N})\right)\in\mathbb{R}^{1\times
N}$ $\displaystyle g(\mathbf{x},y)$
$\displaystyle=\left(g(x_{1},y),\dots,g(x_{N},y)\right)^{\intercal}\in\mathbb{R}^{N\times
1}$

where $g(x,\mathbf{x}^{\intercal})$ is a row-vector and $g(\mathbf{x},y)$ a
column vector.

Throughout the paper, we will write $\mathbf{K}$ for the matrix
$\mathbf{K}=k(\mathbf{x},\mathbf{x}^{\intercal})$ which we assume is
invertible. In practice, if $\mathbf{K}$ is not invertible, one can add a
small constant to the diagonal. Furthermore, we will typically omit the
conditioning on $\mathbf{x}$ and $\mathbf{x}_{\ast}$ for ease of exposition.

###  2.2 Gaussian Process classification

We follow the usual setup for classification with Gaussian Processes in which
each input-output pair $(\mathbf{x},y)$ is related by the assumption that the
conditional distribution $y\mid f(\mathbf{x})$ is a Bernoulli distribution
with probability $\sigma(f(\mathbf{x}))$ where $\sigma(\cdot)$ is the logistic
function and the prior distribution of $f$ is $\mathcal{GP}(m,k)$. For a
training set $\mathcal{D}=\\{(x_{i},y_{i})\mid i=1,\dots,N\\}$ we denote
$\mathbf{f}=(f_{n})_{n=1}^{N}=(f(\mathbf{x}_{1}),\dots,f(\mathbf{x}_{N}))^{\intercal}$
and $\mathbf{y}=(y_{1},\dots,y_{n}){\intercal}$ so that

$\displaystyle\mathbf{y}\mid\mathbf{f}\sim\prod_{n=1}^{N}\sigma(f_{n})^{y_{n}}(1-\sigma(y_{n}))^{1-y_{n}}\,.$
(1)

The log posterior has the form

$\displaystyle\psi(\mathbf{f})\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\log p(\mathbf{f}\mid\mathbf{y})\underset{\mathbf{f}}{\propto}\log
p(\mathbf{f})+\log p(\mathbf{y}\mid\mathbf{f}),$

where $p(\mathbf{f})$ is the pdf of the multivariate Gaussian distribution and
it follows from (1) that

$\displaystyle\log
p(\mathbf{y}\mid\mathbf{f})=\mathbf{y}^{\intercal}\mathbf{f}-\sum_{n=1}^{N}\log(1+\exp(f_{n}))\
.$

Unlike the classification case, the posterior distribution
$\mathbf{f}\mid\mathbf{y}$ is not analytically tractable and needs to be
approximated. This issue is not specific to distillation but arises in any
practical application involving Gaussian process classification, and a number
of approaches have been suggested for obtaining an approximation of the
posterior distribution. In this paper, we will focus on the Laplace
approximation, where the posterior distribution is approximated by

$\displaystyle
q(\mathbf{f}\mid\mathbf{y})=\mathcal{N}\left(\mathbf{f}\mid\hat{\mathbf{f}},(\mathbf{K}^{-1}+\mathbf{W})^{-1}\right)$

where $\hat{\mathbf{f}}=(\hat{f}_{n})_{n=1}^{N}$ is the mode of
$\psi(\mathbf{f})$ and $\mathbf{W}$ is the $N\times N$ matrix
$\operatorname{diag}_{n}\\{\sigma(\hat{f}_{n})(1-\sigma(\hat{f}_{n}))\\}$. The
approximate posterior predictive distribution is found in Rasmussen & Williams
(2006, See Section 3.4.2) in the case where $m=0$ and for a general prior mean
$m$, we have
$p(f_{\ast}\mid\mathbf{x},\mathbf{y})\approx\mathcal{N}(\bm{\mu}_{\ast},\mathbf{\Sigma}_{\ast})$
where

$\displaystyle\bm{\mu}_{\ast}$
$\displaystyle=\mathbb{E}_{q}[\mathbf{f}_{\ast}\mid\mathbf{y},\mathbf{x}_{\ast}]$
$\displaystyle=m(\mathbf{x}_{\ast})+k(\mathbf{x}_{\ast},\mathbf{x}^{\intercal})\mathbf{K}^{-1}(\hat{\mathbf{f}}-m(\mathbf{x}))$
(2) $\displaystyle\mathbf{\Sigma}_{\ast}$
$\displaystyle=\mathrm{Cov}_{q}(\mathbf{f}_{\ast}\mid\mathbf{y},\mathbf{x}_{\ast})$
$\displaystyle=k(\mathbf{x}_{\ast},\mathbf{x}_{\ast}^{\intercal})-k(\mathbf{x}_{\ast},\mathbf{x}^{\intercal})(\mathbf{K}+\mathbf{W}^{-1})^{-1}k(\mathbf{x},\mathbf{x}_{\ast}^{\intercal}).$
(3)

See Section A.1 for additional details on Gaussian Process classification. We
return to the classification setting in Section 4.

##  3 Self-distillation for GP Regression

In the following we consider data- and distribution-centric self-distillation
for GP regression (GPR) separately.

###  3.1 Data-centric Self-Distillation

(a) Data-centric distillation for ten steps.

(b) Distribution-centric distillation for ten steps.

Figure 3: Ten steps of self-distillation with fixed hyperparameters and
$(\gamma_{1},\dots,\gamma_{10})=(0.1,\dots,1)$ for both the data-centric and
distribution-centric approach. We plot both the mean function and the 2.5 and
97.5 percentiles of the uncertainty estimate at each step.

For the general setup of data-centric distillation, we consider a mean-zero
prior distribution, $f\sim\mathcal{GP}(0,k)$, and the posterior distribution
on $\mathbf{x}=(x_{1},\dots,x_{N})^{\intercal}$ at step $t\geq 1$ as
$\mathbf{f}_{t}\mid\mathbf{y}_{t-1}\sim\mathcal{N}(\bm{\mu}_{t},\mathbf{\Sigma}_{t})$,
where

$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\mathbf{K}(\mathbf{K}+\gamma_{t}\mathbf{I}_{N})^{-1}\mathbf{y}_{t-1},$
$\displaystyle\mathbf{\Sigma}_{t}$
$\displaystyle=\mathbf{K}-\mathbf{K}(\mathbf{K}+\gamma_{t}\mathbf{I}_{N})^{-1}\mathbf{K},$

and $\gamma_{t}>0$ for $t\geq 1$. Data-centric self-distillation corresponds
to using the mean predictions $\bm{\mu}_{t}$ in place of $\mathbf{y}_{t}$ in
fitting of the succeeding step, thereby fitting the subsequent model to the
mean predictions of the current model. In this setting the behavior of the
mean function becomes equivalent to that of Mobahi et al. (2020) and Borup &
Andersen (2021) as evident from the following result.

######  Proposition 3.1.

Define $\mathbf{y}_{t}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\bm{\mu}_{t}$ and
$\mathbf{y}_{0}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\mathbf{y}$. Then it holds for any $t\geq 1$ that
$\mathbf{f}_{t}\mid\bm{\mu}_{t-1}\sim\mathcal{N}\left(\bm{\mu}_{t},\mathbf{\Sigma}_{t}\right)$,
where

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\left(\prod_{s=1}^{t}\mathbf{K}(\mathbf{K}+\gamma_{s}\mathbf{I}_{N})^{-1}\right)\mathbf{y},$
(4) $\displaystyle\mathbb{E}[\mathbf{f}_{t,\ast}\mid\bm{\mu}_{t-1}]$
$\displaystyle=k(\mathbf{x}_{\ast},\mathbf{x})(\mathbf{K}+\gamma_{t}\mathbf{I}_{N})^{-1}\mathbf{y}_{t-1},$
(5)

and
$\mathbf{\Sigma}_{t}=\mathbf{K}-\mathbf{K}(\mathbf{K}+\gamma_{t}\mathbf{I}_{N})^{-1}\mathbf{K}$
as well as $\mathbf{y}_{t}$ is only affected by the choice of $\gamma_{t}$.
(Proof in appendix.)

It can be shown that the behavior of $\mathbf{y}_{t}$ as $t$ increases
corresponds to a progressively sparsifying the underlying basis functions by
shrinking the basis functions with smallest eigenvalues the most, thereby
imposing a regularization effect on the predictions (Mobahi et al., 2020;
Borup & Andersen, 2021). Furthermore, since $\mathbf{y}_{t}$ only depend on
$t$ through $\\{\gamma_{s}\\}_{s=1}^{t}$, once can efficiently compute (4) for
any choice of $t$, by using a spectral decomposition of $\mathbf{K}$.

Furthermore, following Borup & Andersen (2021) the results can be extended to
include a weighting (by $\alpha$) of the original training observations and
(by $1-\alpha$) of the predictions from the previous iteration of the model.

###  3.2 Distribution-centric Self-Distillation

The central observation of the distribution-centric approach is that the
posterior distribution is itself a Gaussian Process, namely
$\mathcal{GP}(m_{\ast},k_{\ast})$, where

$\displaystyle m_{\ast}(x)$
$\displaystyle=m(x)+k(x,\mathbf{x}^{\intercal})\left(\mathbf{K}+\gamma\mathbf{I}_{N}\right)^{-1}\left(\mathbf{y}-m(\mathbf{x})\right)$
$\displaystyle k_{\ast}(x,y)$
$\displaystyle=k(x,y)-k(x,\mathbf{x}^{\intercal})\left(\mathbf{K}+\gamma\mathbf{I}_{N}\right)^{-1}k(\mathbf{x},y)\,.$

Since the posterior distribution is a Gaussian process, we may ”distill” the
posterior distribution, by using it as a prior distribution in our initial
set-up, and we may iterate this procedure for any finite number of steps, see
Figure 1 (2(b)). Furthermore, we may select the noise terms $\gamma_{t}$ to be
different in each step, and we see that the distillation procedure leads to a
series of Gaussian Processes $\mathcal{GP}(m_{t},k_{t})$ where the mean- and
covariance function satisfy the recursions

$\displaystyle m_{t+1}(x)$
$\displaystyle=m_{t}(x)+k_{t}(x,\mathbf{x}^{\intercal})[\mathbf{K}_{t}+\gamma_{t}\mathbf{I}_{N}]^{-1}(\mathbf{y}-m_{t}(\mathbf{x}))$
(6) $\displaystyle k_{t+1}(x,y)$
$\displaystyle=k_{t}(x,y)-k_{t}(x,\mathbf{x}^{\intercal})[\mathbf{K}_{t}+\gamma_{t}\mathbf{I}_{N}]^{-1}k_{t}(\mathbf{x},y)$
(7)

where $\mathbf{K}_{t}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}k_{t}(\mathbf{x},\mathbf{x}^{\intercal})$,
$m_{0}(x)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}0$,
$k_{0}(x,y)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}k(x,y)$ and
$\gamma_{t}>0$, $t=0,1,\dots$.

The solutions to the above recursions (6) and (7) are given in Theorem 3.3
below, where we show that iterating the recursions $t$ times with penalty
parameters $\gamma_{0},\gamma_{1},\dots,\gamma_{t-1}$ yields the same result
as performing a single iteration with penalty parameter $1/\gamma_{t-1}^{-}$
where $\gamma_{t}^{-}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\sum_{s=0}^{t}1/\gamma_{s}$ and
$\gamma_{-1}^{-}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}0$. To
show this, we first need a lemma.

######  Lemma 3.2.

Let $\lambda_{i}^{(t)}$, for $i=1,\dots,N$, denote the eigenvalues of
$\mathbf{K}_{t}$. Write
$\mathbf{K}_{0}=\mathbf{O}\operatorname{diag}_{i}(\lambda_{i}^{(0)})\mathbf{O}^{\intercal}$
for the spectral decomposition of $\mathbf{K}_{0}$. Then we have for
$t=0,1,\dots$ that

$\displaystyle\mathbf{K}_{t}=\mathbf{O}\underset{i}{\operatorname{diag}}\left(\lambda_{i}^{(0)}\zeta_{i}^{(t)}\right)\mathbf{O}^{\intercal}$
(8)

where $\zeta_{i}^{(t)}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\frac{1}{\lambda_{i}^{(0)}\gamma_{t-1}^{-}+1}$. For
$k_{t}(x,\mathbf{x}^{\intercal})$ and $m_{t}(\mathbf{x})$ we have

$\displaystyle k_{t}(x,\mathbf{x}^{\intercal})$
$\displaystyle=k_{0}(x,\mathbf{x}^{\intercal})\mathbf{O}\underset{i}{\operatorname{diag}}\left(\zeta_{i}^{(t)}\right)\mathbf{O}^{\intercal}$
(9) $\displaystyle m_{t}(\mathbf{x})$
$\displaystyle=\mathbf{O}\underset{i}{\operatorname{diag}}\left(\lambda_{i}^{(0)}\gamma_{t-1}^{-}\zeta_{i}^{(t)}\right)\mathbf{O}^{\intercal}\mathbf{y}.$
(10)

(Proof in appendix.)

We may now prove the main result of this section.

######  Theorem 3.3.

The solution to the recursions (6) and (7) is given by

$\displaystyle m_{t}(x)$
$\displaystyle=k_{0}(x,\mathbf{x}^{\intercal})(\mathbf{K}+\mathbf{I}_{N}/\gamma_{t-1}^{-})^{-1}\mathbf{y}$
(11) $\displaystyle k_{t}(x,y)$
$\displaystyle=k_{0}(x,y)-k_{0}(x,\mathbf{x}^{\intercal})(\mathbf{K}+\mathbf{I}_{N}/\gamma_{t-1}^{-})^{-1}k_{0}(\mathbf{x},y)$

for $t=1,2,\dots$. (Proof in appendix.)

It follows from Theorem 3.3, that performing multiple steps of distribution-
centric self-distillation is equivalent to fitting an ordinary GP with a
particular choice of noise parameter.

Next, we note that if we assume that the noise parameter is identical in all
distillation steps, i.e. $\gamma_{t}=\gamma$ for $t=0,1,\dots$, then
$\gamma_{t-1}^{-}=t/\gamma$ which corresponds to fitting a GP to a dataset
$\mathcal{D}_{t}$, which consists of $t$ replications of $\mathcal{D}$ with
noise parameter $\gamma$, in a sense made precise by the following corollary.

######  Corollary 3.4.

If we let

$\mathcal{D}_{t}=\\{(x_{ij},y_{ij})\mid i=1,\dots,N\,,\quad j=1,\dots,t\\}$

where $x_{ij}=x_{i}$ for all $i$ and $j$, then the posterior distribution of
$\mathbf{f}$ is $\mathcal{N}(\bm{\mu}_{t},\mathbf{\Sigma}_{t})$ where

$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\begin{pmatrix}&\mathbf{K}(\mathbf{K}+\gamma/t\mathbf{I}_{N})^{-1}\mathbf{y}\\\
&\vdots\\\
&\mathbf{K}(\mathbf{K}+\gamma/t\mathbf{I}_{N})^{-1}\mathbf{y}\end{pmatrix}$
(12) $\displaystyle\mathbf{\Sigma}_{t}$
$\displaystyle=\mathbf{1}\mathbf{1}^{\intercal}\otimes(\mathbf{K}-\mathbf{K}(\mathbf{K}+\gamma/t\mathbf{I}_{N})^{-1}\mathbf{K})$
(13)

where $\otimes$ is the Kronecker product. (Proof in appendix.)

###  3.3 Illustratory Example

We now consider an illustratory example, where we assume a true function
$g(z)=z\sin{(z)}$, 10 training samples $\mathbf{x}$ equidistantly distributed
between 0 and 10, and observed values $y_{i}=g(x_{i})+\varepsilon_{i}$, where
$\varepsilon_{i}$ is standard normally distributed. We use a scaled radial
basis kernel with two hyperparameters $\sigma_{f}$ and $l$ as well as noise
parameters $\gamma_{t}$. See Section B for further details on the setup.

We now consider $10$ steps of both data-centric and distribution-centric self-
distillation with fixed hyperparameters for all steps. In Figure 3(a) and
Figure 3(b) we plot the results for data-centric distillation with
$(\gamma_{1},\dots,\gamma_{10})=(0.1,\dots,1)$. Here, we observe the expected
increase in regularization of the mean, and that the uncertainty increases
with each step as expected from the choice of
$(\gamma_{1},\dots,\gamma_{10})$. This aligns well with our theory. However,
for the distribution-centric case, the mean function is nearly unaffected by
the distillation procedure, and the uncertainty estimate does not change much
either. This is due to the relatively small change in $1/\gamma_{t}^{-}$ when
$t$ goes from $1$ to $10$ as $1/\gamma_{1}^{-}=0.1$ and
$1/\gamma_{10}^{-}=0.034$, which yield a small change for each distillation
step. See Appendix B for additional examples with different choices of
parameters.

##  4 Self-Distillation for GP Classification

In the following we separately investigate data-centric and distribution-
centric self-distillation for GP classification.

###  4.1 Data-centric Self-Distillation

Figure 4: Ordinary GPC model and one step of data-centric distillation, both
with continuous Bernoulli and ordinary Bernoulli likelihood. We note that the
Bernoulli likelihood, although, well-performing does not constitute a well-
specified model.

Figure 5: One step of data-centric distillation with continuous Bernoulli
likelihood as in Figure 5, both with and without regularization imposed by a
noise assumption on the training observations. Adjusting the noise parameter
$\gamma_{1}$ adjust the regularization.

Similarly to the data-centric approach for GPR, we will in this section
consider the posterior mean predictions as deterministic targets for the next
iteration of the model. In particular, at each step, $t$, we assume the prior
of $f_{t}$ to be $\mathcal{GP}(0,k)$, and for $t=1$ we will obtain
predictions, $\hat{\mathbf{y}}_{0}$ as the average posterior prediction based
on the Laplace approximation e.g.
$\hat{\mathbf{y}}_{0}=\mathbb{E}_{q}[\sigma(\mathbf{f}_{1})\mid\mathbf{y}_{0}]$.
These predictions will then be used in place of the original observations for
step $t=2$, which will lead to new predictions $\hat{\mathbf{y}}_{2}$, that we
can use in the succeeding step, and so on.

Since the original observations, $\mathbf{y}_{0}$, are in $\\{0,1\\}^{N}$ and
the predictions, $\hat{\mathbf{y}}_{1}$, are in $[0,1]^{N}$, directly reusing
the predictions, $\hat{\mathbf{y}}_{1}$, in place of the observations,
$\mathbf{y}_{0}$, in the subsequent step yields a misspecified model as the
Bernoulli distribution only has support on $\\{0,1\\}$. To alleviate this
misspecification, we redefine our model for $t\geq 2$ to assume a conditional
continuous Bernoulli distribution rather than the usual conditional Bernoulli
(Loaiza-Ganem & Cunningham, 2019). We will occasionally refer to this method
as C-GPC, and in practice, this requires that we multiply each element of (1)
by the appropriate normalizing constant, $C$ (see (20) and Appendix C). That
is, for $t\geq 2$ we assume

$\displaystyle\mathbf{y}_{t-1}\mid\mathbf{f}_{t}\sim\prod_{n=1}^{N}C(\sigma(f_{t,n}))\sigma(f_{t,n})^{y_{t-1,n}}(1-\sigma(f_{t,n}))^{1-y_{t-1,n}},$

which in turn affects the log posterior
$\log\tilde{p}(\mathbf{f}_{t}\mid\mathbf{y}_{t-1})$ through
$\log\tilde{p}(\mathbf{y}_{t-1}\mid\mathbf{f}_{t})$ which now becomes

$\displaystyle\log\tilde{p}(\mathbf{y}_{t-1}\mid\mathbf{f}_{t})$
$\displaystyle=\log p(\mathbf{y}_{t-1}\mid\mathbf{f}_{t})+\sum_{n=1}^{N}\log
C(\sigma(f_{t,n})).$

Furthermore, this also affects the gradient and hessian of
$\log\tilde{p}(\mathbf{f}_{t}\mid\mathbf{y}_{t-1})$ used to obtain the Laplace
approximation, and although $\log C(\lambda)$ nor $\frac{d}{d\lambda}\log
C(\lambda)$ attains no simple expressions, we show in Proposition 4.1 that
both $C(\sigma(a))$ and the derivatives $\frac{d}{da}\log C(\sigma(a))$ and
$\frac{d^{2}}{d^{2}a}\log C(\sigma(a))$ yield surprisingly simple expressions
of standard hyperbolic functions (see also Figure 14 for plots of the
functions).

######  Proposition 4.1.

Let $C(\sigma(a))$ be defined as in (20) with $\lambda=\sigma(a)$, then we
have that

$\displaystyle C(\sigma(a))$ $\displaystyle=\begin{cases}2&\text{if }a=0\\\
a\mathrm{coth}\left(\frac{a}{2}\right)&\text{otherwise},\end{cases}$
$\displaystyle\frac{d}{da}\log(C(\sigma(a)))$
$\displaystyle=\begin{cases}0&\text{if }a=0\\\
\frac{1}{a}-\frac{1}{\mathrm{sinh}(a)}&\text{otherwise},\end{cases}$
$\displaystyle\frac{d^{2}}{d^{2}a}\log(C(\sigma(a)))$
$\displaystyle=\begin{cases}\frac{1}{6}&\text{if }a=0\\\
-\frac{1}{a^{2}}+\frac{\mathrm{coth}(a)}{\mathrm{sinh}(a)}&\text{otherwise},\end{cases}$

where $\mathrm{coth}$, and $\mathrm{sinh}$ are the hyperbolic cotangent, and
sine functions, respectively. (Proof in appendix.)

Thus, if we denote
$\mathbf{C}_{t}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\sum_{n=1}^{N}\log(C(\sigma(f_{t,n})))$, we get that the gradient
and hessian of the log posterior becomes

$\displaystyle\nabla\tilde{\psi}(\mathbf{f}_{t})$
$\displaystyle=\mathbf{y}_{t-1}-\bm{\sigma}_{t-1}-\mathbf{K}^{-1}\mathbf{f}_{t}+\nabla\mathbf{C}_{t},$
$\displaystyle-\nabla\nabla\tilde{\psi}(\mathbf{f}_{t})$
$\displaystyle=\mathbf{W}_{t}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{t},$

where
$\mathbf{W}_{t}=\operatorname{diag}_{n}\\{\sigma(f_{t,n})(1-\sigma(f_{t,n}))\\}$
and
$\nabla\nabla\mathbf{C}_{t}=\operatorname{diag}_{n}\\{\frac{d^{2}}{d^{2}f_{t,n}}\log
C(\sigma(f_{t,n}))\\}$, respectively. These directly affect the Laplace
approximation of the posterior as is also evident from e.g. Figure 5.

Unfortunately, unlike data-centric self-distillation for GPR, we are unable to
obtain a closed-form solution for the distilled models due to the numerical
approximations necessary to obtain the posterior (predictions). However, we do
empirically investigate the effect of using the continuous Bernoulli rather
than the discrete Bernoulli distribution in Figure 5. In particular, we plot
the solutions for a single step of distillation with C-GPC and with the
discrete Bernoulli where we use continuous observations despite the support of
the Bernoulli distribution being discrete. We expand on this example in
Section 4.3 below.

###  4.2 Distribution-centric Self-Distillation

(a) Solutions for 10 steps of distillation.

(b) Approximation error between iterated and scaled solutions.

Figure 6: Ten steps of distribution-centric self-distillation for GP
classification. In (6(a)) we plot the solutions using the scaled covariance
approximation, and in (6(b)) we plot the approximation error between the
iterated distillation procedure and scaled approximation over the test set of
90 equidistant point on the interval $[-2,7]$.

In the distribution-centric approach, we proceed in a similar manner as in
Section 3.2. Firstly, we observe from (2) and (3) that our approximate
posterior distribution is a Gaussian process $\mathcal{GP}(m_{\ast},k_{\ast})$
with

$\displaystyle m_{\ast}(x)$
$\displaystyle=m(\mathbf{x})+k(\mathbf{x},\mathbf{x}^{\intercal})\mathbf{K}^{-1}(\hat{\mathbf{f}}-m(\mathbf{x}))$
$\displaystyle k_{\ast}(x,y)$
$\displaystyle=k(x,y)-k(x,\mathbf{x}^{\intercal})(\mathbf{K}+\mathbf{W}^{-1})^{-1}k(\mathbf{x},y)\,.$

This may then be used as a prior distribution, which again implies a series of
Gaussian processes $\mathcal{GP}(m_{t},k_{t})$, $t=0,1,2,\dots$ where analogue
to equations (6) and (7) in this case becomes

$\displaystyle m_{t+1}(x)$
$\displaystyle=m_{t}(x)+k_{t}(x,\mathbf{x}^{\intercal})\mathbf{K}_{t}^{-1}(\hat{\mathbf{f}}_{t}-m_{t}(\mathbf{x}))$
(14) $\displaystyle k_{t+1}(x,y)$
$\displaystyle=k_{t}(x,y)-k_{t}(x,\mathbf{x}^{\intercal})(\mathbf{K}_{t}+\mathbf{W}^{-1}_{t})^{-1}k_{t}(\mathbf{x},y).$
(15)

Here, $\hat{\mathbf{f}}_{t}$ is the mode of $\psi_{t}(\mathbf{f})$, the log
posterior for the $t$’th iteration, and similarly
$\mathbf{W}_{t}=\operatorname{diag}_{n}\\{\sigma(\hat{f}_{t,n})(1-\sigma(\hat{f}_{t,n}))\\}$.

Since $\hat{\mathbf{f}}_{t}$ is not given analytically, we cannot solve the
above recursions in the classification case, as we could in the regression
setting and must resort to numerical results. However, we show an analogous
result to the data duplication result from Corollary 3.4 and it follows that
iterated distillation in the distribution-centric setup may be approximated by
scaling the covariance in the first distillation step. These observations are
summarized in the following proposition, where $\mathcal{D}_{t}$ as in
Corollary 3.4 consists of the dataset $\mathcal{D}$ replicated $t$ times, and
an example is illustrated in Figure 6(a) and 6(b).

######  Proposition 4.2.

$(\mathbf{A})$ A single step of distribution-centric distillation using
$\mathcal{D}_{t}$ and a Gaussian prior $\mathcal{GP}(0,k)$ yields the same
posterior distribution as a single step of distribution-centric distillation
using a Gaussian prior $\mathcal{GP}(0,tk)$. $(\mathbf{B})$ Performing $t$
iterations of the distribution-centric distillation using $\mathcal{D}$ and
starting with an inital Gaussian prior $\mathcal{GP}(0,k)$ yields
approximately the same posterior as a single step of distribution-centric
distillation using a Gaussian prior $\mathcal{GP}(0,tk)$. (Proof in appendix.)

###  4.3 Illustratory Example

We now consider an illustratory example, where we consider the true underlying
function $g(x)=2\mathrm{sin}\left(\frac{x\pi}{2}\right)$. We sample training
observations, $\mathbf{x}$, from a $\mathcal{U}(0,5)$ distribution and binary
targets from a Bernoulli distribution with probability parameterized by
$g(x)$.

#### Data-centric GPC distillation

For data-centric distillation, we initially fit an ordinary GP classification
model with Bernoulli distribution to the binary observations and fit a
secondary model to the continuous predictions of this first model. In Figure 5
we plot the predictions of the initial ordinary GPC, along with the
predictions of the secondary model, when we use either a continuous Bernoulli
distribution (as argued in Section 4) or a discrete Bernoulli distribution
where we intentionally violate the discrete support of the distribution. The
hyperparameters of the models are chosen to minimize the negative log-
likelihood (NLL) over the input data, and while the distilled Bernoulli model
is slightly flattened, yet well-performing, compared to the initial fit, the
C-GPC model yields more ”extreme” predictions. We conjecture that the implicit
regularization imposed by the misspecification of the distilled Bernoulli
model is absent in the C-GPC model. Furthermore, by Section C we expect the
continuous Bernoulli to favor more extreme values compared to the ordinary
Bernoulli, as is also evident in Figure 5. To dampen the degree of overfitting
in the C-GPC model we reintroduce the noise parameter $\gamma_{t}>0$ along the
diagonal of the kernel matrix of the training data. We plot both the
regularized and non-regularized solutions in Figure 5, and observe that with
the noise parameter, we get not only a well-specified model, but also get more
control of the degree of regularization. See Appendix B for more experiments.

#### Distribution-centric GPC distillation

Similarly to the data-centric case, we initially fit an ordinary GPC, but
unlike the data-centric case, for distribution-centric distillation, we do not
change the model assumptions for the distillation step. In Figure 6(a) we plot
the solutions for 10 steps of distillation and observe a significantly
different behavior compared to data-centric distillation; the effect of
distillation is much more subtle and more closely resembles that of
distillation in the regression case and the solutions progressively fit closer
and closer to the original data. Furthermore, in Figure 6(b) we plot the mean
squared approximation error between the solutions obtained from (6)-(7) and
those obtained by scaling the kernel matrix as argued in Proposition 4.2.
Although the error is almost linearly increasing, it does so by less than
$0.5\times 10^{-5}$ per step, and the approximation error is on the order of
$10^{-5}$ for 10 steps. Thus, since the scaling approximation is much faster
(especially for many distillation steps) it is a great alternative for naïve
computations.

##  5 Related Works

To the best of our knowledge, no literature on knowledge distillation for
Gaussian process models exists. However, for good measure in the following we
relate our proposed methods to the existing related literature on knowledge
distillation and Gaussian process regression/classification.

#### Knowledge Distillation

The idea of knowledge distillation originates back to Ba & Caruana (2014);
Bucil & Caruana (2015), but is typically most known from Hinton et al. (2015).
It was originally aimed at model compression, by training a small (in an
appropriate measure of model size) student model on the outputs of a larger
teacher model.

In recent years, the idea of distillation has expanded in numerous directions,
and for various types of applications. Of particular relevance is the
direction of self-distillation, where the student model is considered to be of
the same model class/architecture as the teacher model (Furlanello et al.,
2018; Zhang et al., 2019; Allen-Zhu & Li, 2020). Other interesting directions
of research on distillation are on the impact of the size-gap between teacher
and student (Mirzadeh et al., 2019), matching on other information than the
output (Park et al., 2019; Romero et al., 2015; Srinivas & Fleuret, 2018), as
well as multi-teacher and multi-source distillation (Li et al., 2021; Liu &
Zenke, 2020; You et al., 2017; Anonymous Authors, 2023). Finally, due to the
domain-agnostic nature of most distillation methods, they are applied in a
wide range of domains (Krishnan Parthasarathi & Strom, 2019; Borup et al.,
2023).

Despite much empirical success of a wide array of distillation techniques, a
rigorous understanding of distillation is still lacking behind. However,
recent work has provided insights into simplified settings such as linear
models, kernel ridge regression, and decision trees (Phuong & Lampert, 2019;
Borup & Andersen, 2021; Frosst & Hinton, 2018).

#### Gaussian Processes

Many results on Gaussian Process models have been well-known for years and are
collected in foundational books such as Rasmussen & Williams (2006); Bishop
(2007); Murphy (2013). However, in recent years, an increasing amount of
research on extending the capabilities of GP models and adapting them to
modern compute availability and dataset sizes has been made. This includes
theoretical connections between neural networks and Gaussian processes (Lee et
al., 2018; Yang, 2019; Garriga-Alonso et al., 2019), stochastic and sparse GPs
(Titsias, ; Hensman et al., 2013, 2015), and the possibility to exploit the
successes of neural networks in combination with Gaussian processes e.g. deep
kernel learning (Wilson et al., 2016a, b).

##  6 Conclusion

In conclusion, we propose two approaches to extend knowledge distillation to
Gaussian process regression and classification; data-centric and distribution-
centric. We show that some of these approaches are closely related to known
results or very particular choices of hyperparameters, but also that some
methods require careful redefinition of our model assumptions to be
appropriate. To the best of our knowledge, this paper is the first attempt to
propose knowledge distillation specifically for Gaussian Process models,
opening new avenues for further research with potential implications in
machine learning. Interesting directions of future research include, amongst
others, utilization of a convex combination of both the original and predicted
targets in the distillation steps, as well as combinations of data-centric and
distribution-centric distillation approaches.

## References

  * Allen-Zhu & Li (2020) Allen-Zhu, Z. and Li, Y.  Towards Understanding Ensemble, Knowledge Distillation and Self-Distillation in Deep Learning.  2020\.  URL http://arxiv.org/abs/2012.09816.
  * Anonymous Authors (2023) Anonymous Authors.  Learning Efficient Models From Few Labels By Distillation From Multiple Tasks.  _Unpublished/Under review_ , 2023.
  * Ba & Caruana (2014) Ba, L. J. and Caruana, R.  Do Deep Nets Really Need to be Deep?  _Advances in Neural Information Processing Systems_ , 3(January):2654–2662, 2014.  ISSN 10495258.
  * Bishop (2007) Bishop, C. M.  _Pattern Recognition and Machine Learning (Information Science and Statistics)_.  Springer, 2007.  ISBN 0387310738.
  * Borup & Andersen (2021) Borup, K. and Andersen, L. N.  Even your Teacher Needs Guidance: Ground-Truth Targets Dampen Regularization Imposed by Self-Distillation.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P. S., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 5316–5327. Curran Associates, Inc., 2021.  URL http://arxiv.org/abs/2102.13088https://proceedings.neurips.cc/paper/2021/file/2adcefe38fbcd3dcd45908fbab1bf628-Paper.pdf.
  * Borup et al. (2023) Borup, K., Kidmose, P., Phan, H., and Mikkelsen, K.  Automatic sleep scoring using patient-specific ensemble models and knowledge distillation for ear-EEG data.  _Biomedical Signal Processing and Control_ , 81(May 2022):104496, 2023.  ISSN 17468108.  doi: 10.1016/j.bspc.2022.104496.  URL https://doi.org/10.1016/j.bspc.2022.104496.
  * Bucil & Caruana (2015) Bucil, C. and Caruana, R.  Model Compression Cristian.  _IEEE Communications Magazine_ , 54(1):1–9, 2015\.  ISSN 0218-0014.  doi: 10.1111/j.1467-985x.2009.00624{\\_}10.x.  URL http://users.iems.northwestern.edu/~nocedal/PDFfiles/dss.pdf%0Ahttp://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.67.7952&rep=rep1&type=pdf%255Cnhttp://books.nips.cc/papers/files/nips15/AA03.pdf%250Ahttps://ieeexplore.ieee.org/document/85783.
  * Frosst & Hinton (2018) Frosst, N. and Hinton, G. r.  Distilling a neural network into a soft decision tree, 2018.  ISSN 16130073.
  * Furlanello et al. (2018) Furlanello, T., Lipton, Z. C., Tschannen, M., Itti, L., and Anandkumar, A.  Born Again Neural Networks.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 1607–1616, Stockholmsmässan, Stockholm Sweden, 2018. PMLR.
  * Garriga-Alonso et al. (2019) Garriga-Alonso, A., Aitchison, L., and Rasmussen, C. E.  Deep convolutional networks as shallow Gaussian processes.  _7th International Conference on Learning Representations, ICLR 2019_ , pp. 1–16, 2019.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  In _Uncertainty in Artificial Intelligence - Proceedings of the 29th Conference, UAI 2013_ , pp. 282–290, 2013.
  * Hensman et al. (2015) Hensman, J., Matthews, A. G., and Ghahramani, Z.  Scalable variational Gaussian process classification.  _Journal of Machine Learning Research_ , 38:351–360, 2015\.  ISSN 15337928.
  * Hinton et al. (2015) Hinton, G., Vinyals, O., and Dean, J.  Distilling the Knowledge in a Neural Network.  _arXiv preprint arXiv:1503.02531_ , pp. 1–9, 2015.  URL http://arxiv.org/abs/1503.02531.
  * Krishnan Parthasarathi & Strom (2019) Krishnan Parthasarathi, S. H. and Strom, N.  Lessons from Building Acoustic Models with a Million Hours of Speech.  _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings_ , 2019-May:6670–6674, 2019.  ISSN 15206149.  doi: 10.1109/ICASSP.2019.8683690.
  * Lee et al. (2018) Lee, J., Bahri, Y., Novak, R., Schoenholz, S. S., Pennington, J., and Sohl-Dickstein, J.  Deep neural networks as Gaussian processes.  _6th International Conference on Learning Representations, ICLR 2018 - Conference Track Proceedings_ , 2018.
  * Li et al. (2021) Li, Z., Ravichandran, A., Fowlkes, C., Polito, M., Bhotika, R., and Soatto, S.  Representation Consolidation for Training Expert Students, 2021.  URL http://arxiv.org/abs/2107.08039.
  * Liu & Zenke (2020) Liu, T. and Zenke, F.  Finding trainable sparse networks through Neural Tangent Transfer, 2020\.
  * Loaiza-Ganem & Cunningham (2019) Loaiza-Ganem, G. and Cunningham, J. P.  The continuous bernoulli: Fixing a pervasive error in variational autoencoders.  _Advances in Neural Information Processing Systems_ , 32(NeurIPS):1–11, 2019.  ISSN 10495258.
  * Mirzadeh et al. (2019) Mirzadeh, S.-I., Farajtabar, M., Li, A., Ghasemzadeh, H., Levine, N., Matsukawa, A., and Ghasemzadeh, H.  Improved Knowledge Distillation via Teacher Assistant: Bridging the Gap Between Student and Teacher.  _arXiv preprint arXiv:1902.03393_ , 2019.  URL http://arxiv.org/abs/1902.03393.
  * Mobahi et al. (2020) Mobahi, H., Farajtabar, M., and Bartlett, P. L.  Self-Distillation Amplifies Regularization in Hilbert Space.  _Advances in Neural Information Processing Systems_ , 2020.
  * Murphy (2013) Murphy, K. P.  _Machine learning: A probabilistic perspective_.  MIT Press, 2013.  ISBN 9780262018029 0262018020.
  * Park et al. (2019) Park, W., Kim, D., Lu, Y., and Cho, M.  Relational Knowledge Distillation.  In _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , 2019.  ISBN 9781728132938.  doi: 10.1109/CVPR.2019.00409.
  * Phuong & Lampert (2019) Phuong, M. and Lampert, C. H.  Towards understanding knowledge distillation.  _36th International Conference on Machine Learning, ICML 2019_ , 2019-June(2014):8993–9007, 2019.
  * Rasmussen & Williams (2006) Rasmussen, C. E. and Williams, K.  _Gaussian Processes for Machine Learning_.  2006\.
  * Romero et al. (2015) Romero, A., Ballas, N., Kahou, S. E., Chassang, A., Gatta, C., and Bengio, Y.  FitNets: Hints for thin deep nets.  _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ , pp. 1–13, 2015.
  * Srinivas & Fleuret (2018) Srinivas, S. and Fleuret, F. F.  Knowledge transfer with jacobian matching.  _35th International Conference on Machine Learning, ICML 2018_ , 11:7515–7523, 2018.
  * (27) Titsias, M. K.  Variational Learning of Inducing Variables in Sparse Gaussian Processes Michalis.
  * Wilson et al. (2016a) Wilson, A. G., Hu, Z., Salakhutdinov, R., and Xing, E. P.  Deep kernel learning.  _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016_ , (1998):370–378, 2016a.
  * Wilson et al. (2016b) Wilson, A. G., Hu, Z., Salakhutdinov, R., and Xing, E. P.  Stochastic variational deep kernel learning.  _Advances in Neural Information Processing Systems_ , (Nips):2594–2602, 2016b.  ISSN 10495258.
  * Yang (2019) Yang, G.  Scaling Limits of Wide Neural Networks with Weight Sharing: Gaussian Process Behavior, Gradient Independence, and Neural Tangent Kernel Derivation, 2019.  URL http://arxiv.org/abs/1902.04760.
  * You et al. (2017) You, S., Xu, C., Xu, C., and Tao, D.  Learning from Multiple Teacher Networks.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’17, pp. 1285–1294, New York, NY, USA, 2017. Association for Computing Machinery.  ISBN 9781450348874.  doi: 10.1145/3097983.3098135.  URL https://doi.org/10.1145/3097983.3098135.
  * Zhang et al. (2019) Zhang, L., Song, J., Gao, A., Chen, J., Bao, C., and Ma, K.  Be your own teacher: Improve the performance of convolutional neural networks via self distillation.  _Proceedings of the IEEE International Conference on Computer Vision_ , 2019-Octob:3712–3721, 2019.  ISSN 15505499.  doi: 10.1109/ICCV.2019.00381.

##  Appendix A Additional details on Self-distillation for GP classification

In the following, we restate some known results for GP classification with
more details than in the main paper. We also provide some additional details
on our self-distillation results for GP classification.

###  A.1 Gaussian Process Classification

Following the usual setup for classification with Gaussian Processes, we
assume $f\sim\mathcal{GP}(m,k)$, and

$\displaystyle\mathbf{y}\mid\mathbf{f}\sim\prod_{n=1}^{N}\sigma(f_{n})^{t_{n}}(1-\sigma(f_{n}))^{1-t_{n}},$

and denote $\mathbf{y}=(y_{1},\dots,y_{n})^{\intercal}$ with
$y_{i}\in\\{0,1\\}$ and similarly $\mathbf{f}$ the vector of realizations of
the GP-prior at the respective $x_{i}$ values. Then we have the log density

$\displaystyle\psi(\mathbf{f})\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\log p(\mathbf{f}\mid\mathbf{y})\propto\log p(\mathbf{f})+\log
p(\mathbf{y}\mid\mathbf{f}),$

where $\log p(\mathbf{f})$ is a multidimensional Gaussian, and one can show
that

$\displaystyle\log
p(\mathbf{y}\mid\mathbf{f})=\mathbf{y}^{\intercal}\mathbf{f}-\sum_{n=1}^{N}\log(1+\exp(f_{n}))\,,$

so that (up to the proportionality in $\mathbf{f}$):

$\displaystyle\psi(\mathbf{f})=-\frac{1}{2}(\mathbf{f}_{N}-\mathbf{m})^{\intercal}\mathbf{K}^{-1}(\mathbf{f}-\mathbf{m})+\mathbf{y}^{\intercal}\mathbf{f}-\sum_{n=1}^{N}\log(1+\exp(f_{n}))\,.$
(16)

Since the posterior is non-Gaussian, we apply Laplace approximation to get a
Gaussian approximation.

Thus, we need the gradient and Hessian to obtain $\hat{\mathbf{f}}$ (the MAP
estimate) and the covariance matrix. It follows from (16) that

$\displaystyle\nabla\psi(\mathbf{f})$
$\displaystyle=\mathbf{y}-\bm{\sigma}-\mathbf{K}^{-1}\mathbf{f}+\mathbf{K}^{-1}\mathbf{m}$
$\displaystyle-\nabla\nabla\psi(\mathbf{f})$
$\displaystyle=\mathbf{W}+\mathbf{K}^{-1}$

where $(\bm{\sigma})_{n}=\sigma(f_{n})$ (this term comes from
$\sum_{n=1}^{N}\log(1+{\mathrm{e}}^{f_{n}})$) and
$\mathbf{W}=\operatorname{diag}_{n}\\{\sigma(f_{n})(1-\sigma(f_{n}))\\}$.

Using these can use the Newton-Raphson procedure to find $\hat{\mathbf{f}}$,
and it follows from the expression for the gradient, that $\hat{\mathbf{f}}$
will be the solution to

$\displaystyle\hat{\mathbf{f}}=\mathbf{m}+\mathbf{K}(\mathbf{y}-\bm{\sigma})\,.$
(17)

The Gaussian approximation of $p(\mathbf{f}\mid\mathbf{y})$ comes out to

$q(\mathbf{f})=\mathcal{N}(\mathbf{f}\mid\hat{\mathbf{f}},(\mathbf{W}+\mathbf{K}^{-1})^{-1})$

The posterior predictive distribution is computed as

$\displaystyle p(f_{\ast}\mid\mathbf{x},\mathbf{y},\mathbf{x}_{\ast})=\int
p(f_{\ast}\mid\mathbf{f})p(\mathbf{f}\mid\mathbf{y})\mathrm{d}\mathbf{f}$ (18)

(see Bishop (2007)) where the $\mathcal{GP}(m,k)$ prior gives the conditional
distribution

$\displaystyle p(f_{\ast}\mid\mathbf{f})=\mathcal{N}(f_{\ast}\mid
m(x)+\mathbf{k}^{T}\mathbf{K}^{-1}(\mathbf{f}-\mathbf{m}),k-\mathbf{k}^{T}\mathbf{K}^{-1}\mathbf{k})\,.$
(19)

We can obtain predictions in two different ways: Either by substituting
$q(\mathbf{f})$ for $p(\mathbf{f}\mid\mathbf{y})$ in (19) (this leads to (2)
and (3)), or simply as the predictions of the mean

$\displaystyle\sigma(\mathbb{E}_{q}[\mathbf{f}\mid\mathbf{y}])=\sigma(\hat{\mathbf{f}}).$

We note that while the two are different in general, the decision boundary is
identical for binary classification (Bishop, 2007, Sec. 10.3).

###  A.2 Data-centric self-distillation for GP classification

For a single step of data-centric distillation, we can choose either of the
above predictions as our new targets and will denote them as $\mathbf{y}_{1}$
irrespective of the choice. However, while $\mathbf{y}\in\\{0,1\\}^{N}$ we now
have $\mathbf{y}_{1}\in[0,1]^{N}$, and we can no longer assume a conditional
Bernoulli distribution over $\mathbf{y}_{1}$. However, by extending to the
continuous Bernoulli distribution (Loaiza-Ganem & Cunningham, 2019), we can
avoid this problem, by introducing the appropriate normalizing constant

$\displaystyle C(\lambda)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\begin{cases}2&\text{if }\lambda=\frac{1}{2},\\\
\frac{2\mathrm{tanh}^{-1}(1-2\lambda)}{1-2\lambda}&\text{otherwise}.\end{cases}$
(20)

See Figure 14 for plots of the normalizing constant, its gradient, and
Hessian. We now assume $f_{1}\sim\mathcal{GP}(0,k)$ and

$\displaystyle\mathbf{y}_{1}\mid\mathbf{f}_{1}\sim\prod_{n=1}^{N}C(\sigma(f_{1,n}))\sigma(f_{1,n})^{y_{1,n}}(1-\sigma(f_{1,n}))^{1-y_{1,n}},$

which in turn yields that

$\displaystyle\log
p(\mathbf{y}_{1}\mid\mathbf{f}_{1})=\mathbf{y}_{1}^{\intercal}\mathbf{f}_{1}-\sum_{n=1}^{N}\log(1+\exp(f_{1,n}))+\sum_{n=1}^{N}\log
C(\sigma(f_{1,n})),$

where we note that $\sum_{n=1}^{N}\log C(\sigma(f_{1,n}))>0$. In Proposition
4.1 we find the derivative and second derivative of $\log C(\sigma(a))$. Thus,
if we denote $\mathbf{C}_{1}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\sum_{n=1}^{N}\log(C(\sigma(f_{1,n})))$, we get that the gradient of
the conditional density is

$\displaystyle\nabla\psi_{1}(\mathbf{f}_{1})$
$\displaystyle=\mathbf{y}_{1}-\bm{\sigma}_{1}-\mathbf{K}^{-1}\mathbf{f}_{1}+\nabla\mathbf{C}_{1},$

where we can compute the last term by use of the derivative above. Similarly,
it follows that the Hessian can be computed as

$\displaystyle-\nabla\nabla\psi_{1}(\mathbf{f}_{1})$
$\displaystyle=\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1},$

where we note that $\nabla\nabla\mathbf{C}_{1}$ is a diagonal matrix, and that
$\bm{\sigma}_{1}$ and $\mathbf{W}_{1}$ are defined analogously to the
classical case. With the gradient and Hessian, we can follow the usual setup
and compute the MAP estimate for the mean of the Laplace approximation, and
use the inverse Hessian as covariance. That is, we use

$\displaystyle\mathbf{f}_{1}^{\text{new}}$
$\displaystyle=\mathbf{f}_{1}-\left(-\mathbf{W}_{1}-\mathbf{K}^{-1}+\nabla\nabla\mathbf{C}_{1}\right)^{-1}\left(\mathbf{y}_{1}-\bm{\sigma}_{1}-\mathbf{K}^{-1}\mathbf{f}_{1}+\nabla\mathbf{C}_{1}\right)$
$\displaystyle=\left(\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}\right)^{-1}\left(\left(\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}\right)\mathbf{f}_{1}+\mathbf{y}_{1}-\bm{\sigma}_{1}-\mathbf{K}^{-1}\mathbf{f}_{1}+\nabla\mathbf{C}_{1}\right)$
$\displaystyle=\mathbf{K}\mathbf{K}^{-1}\left(\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}\right)^{-1}\left(\left(\mathbf{W}_{1}-\nabla\nabla\mathbf{C}_{1}\right)\mathbf{f}_{1}+\mathbf{y}_{1}-\bm{\sigma}_{1}+\nabla\mathbf{C}_{1}\right)$
$\displaystyle=\mathbf{K}\left(\left(\mathbf{W}_{1}-\nabla\nabla\mathbf{C}_{1}\right)\mathbf{K}+\mathbf{I}\right)^{-1}\left(\left(\mathbf{W}_{1}-\nabla\nabla\mathbf{C}_{1}\right)\mathbf{f}_{1}+\mathbf{y}_{1}-\bm{\sigma}_{1}+\nabla\mathbf{C}_{1}\right)$

That is, we now have the Laplace approximation

$\displaystyle p(\mathbf{f}_{1}\mid\mathbf{y})\approx
q(\mathbf{f}_{1}\mid\mathbf{y})=\mathcal{N}\left(\mathbf{f}_{1}\mid\hat{\mathbf{f}}_{1},\left(\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}\right)^{-1}\right),$

where $\hat{\mathbf{f}}_{1}$ is the MAP estimate from above, and we can also
use the Laplace approximation to get the marginal log-likelihood. First by a
Taylor approximation of $\psi_{1}(\mathbf{f}_{1})$ at $\hat{\mathbf{f}}_{1}$,
we get that
$\psi_{1}(\mathbf{f}_{1})\approx\psi_{1}(\hat{\mathbf{f}}_{1})-\frac{1}{2}(\mathbf{f}_{1}-\hat{\mathbf{f}}_{1})^{\intercal}\mathbf{H}^{-1}(\mathbf{f}_{1}-\hat{\mathbf{f}}_{1})$,
where $\mathbf{H}^{-1}$ is the covariance matrix of $q$, i.e.
$\mathbf{H}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}$. Thus, we
have that

$\displaystyle\log p(\mathbf{y}_{1})$
$\displaystyle=\int\psi(\mathbf{f}_{1})d\mathbf{f}_{1}$
$\displaystyle\approx\psi_{1}(\hat{\mathbf{f}}_{1})+\int-\frac{1}{2}(\mathbf{f}_{1}-\hat{\mathbf{f}}_{1})^{\intercal}\mathbf{H}^{-1}(\mathbf{f}_{1}-\hat{\mathbf{f}}_{1})d\mathbf{f}_{1}$
$\displaystyle=\psi_{1}(\hat{\mathbf{f}}_{1})+\frac{N}{2}\log(2\pi)-\frac{1}{2}\log|\mathbf{H}|$
$\displaystyle=\mathbf{y}_{1}^{\intercal}\hat{\mathbf{f}}_{1}-\sum_{n=1}^{N}\log(1+\exp(\hat{f}_{1,n}))+\sum_{n=1}^{N}\log
C(\sigma(\hat{f}_{1,n}))-\frac{1}{2}\hat{\mathbf{f}}_{1}^{\intercal}\mathbf{K}^{-1}\hat{\mathbf{f}}_{1}$
$\displaystyle\qquad-\frac{1}{2}\log|\mathbf{K}|-\frac{N}{2}\log(2\pi)+\frac{N}{2}\log(2\pi)-\frac{1}{2}\log|\mathbf{H}|$
$\displaystyle=\mathbf{y}_{1}^{\intercal}\hat{\mathbf{f}}_{1}-\sum_{n=1}^{N}\log(1+\exp(\hat{f}_{1,n}))+\sum_{n=1}^{N}\log
C(\sigma(\hat{f}_{1,n}))$
$\displaystyle\qquad-\frac{1}{2}\hat{\mathbf{f}}_{1}^{\intercal}\mathbf{K}^{-1}\hat{\mathbf{f}}_{1}-\frac{1}{2}\log(|\mathbf{K}|)-\frac{1}{2}\log(|\mathbf{W}_{1}+\mathbf{K}^{-1}-\nabla\nabla\mathbf{C}_{1}|)$

The marginal log-likelihood is useful for hyperparameter optimization.

We can then iterate this distillation procedure any number of times. Usually,
we consider the inversion of $\mathbf{K}$ the most computationally expensive
step of fitting a GP, and since we can re-use $\mathbf{K}^{-1}$ for all our
steps, the addition of distillation steps is, computationally, not very
demanding. In the above MAP estimation, we even rewrite the optimization to
avoid the use of the inverse of $\mathbf{K}$.

##  Appendix B Setup for toy-examples and additional results

Throughout the main paper we have included illustrations and plots based on
various toy-examples. However, all these illustratory examples have been
computed with a radial basis kernel function

$\displaystyle
k(\mathbf{x}_{1},\mathbf{x}_{2})=\sigma_{f}^{2}\exp\left(-\frac{\left\lVert\mathbf{x}_{1}-\mathbf{x}_{2}\right\rVert^{2}}{2l}\right),$

where $\sigma_{f}^{2}>0$ and $l>0$ are both hyperparameters tunable for the
particular example. We also (mainly for the regression setup) consider a noise
parameter $\gamma\geq 0$ which is added to the diagonal of $\mathbf{K}$, i.e.
assume noise on the training observations. Even without this assumptions we
typically use a small $\gamma\approx 10^{-8}$ to avoid numerical instability
both in regression and classification examples.

In Figure 7 we plot the negative log-likelihood over a grid of values for
$\sigma_{f}$ and $l$ both in the case of using the ordinary GP classification
with Bernoulli likelihood and the adjusted setting with continuous Bernoulli
likelihood, when the input are the true continuous underlying function
$g(x)=2\sin(x\pi/2)$.

(a) Ordinary GPC

(b) C-GPC

Figure 7: Grid search of hyperparameters of two different types of GP
classification. We indicate the minimum of each grid with a white dot. Note,
the absolute values can not be compared between the two models due to
different likelihood functions.

###  B.1 Additional experiments and results

In the following, we present additional results on the same illustratory
experiments as in the main paper, but with changes to e.g. the model
parameters or data-assumptions. In particular, in Figure 8 we present
different solutions to the data-centric GPR example in the main text, when
varying the choices of $(\gamma_{1},\dots,\gamma_{10})$. Similarly, in Figure
9 we repeat the same experiments for distribution-centric GPR. Finally, in
Figure 10 we plot the solutions of data-centric GPC, when we fit the distilled
model, not to the continuous predictions of the initial model, but to the
$(0,1)$-encoded targets obtained by thresholding the initial model at $0.5$.
In this setting, using a ordinary Bernoulli assumption does not result in a
misspecified model, but the thresholded neglects some of the information
contained in the continuous targets, and as evident from Figure10, the
distilled models in this case is largely similar to that of a C-GPC on the
continuous targets.

(a) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(0.2,\dots,0.2)$

(b) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(1,\dots,0.1)$

(c) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(0.1,\dots,3)$

(d) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(3,\dots,0.1)$

Figure 8: Ten steps of data-centric self-distillation for GPR with fixed
hyperparameters, but with varying choices of noise parameters. We note the
progressively increasing regularization present in all examples.

(a) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(0.2,\dots,0.2)$

(b) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(1,\dots,0.1)$

(c) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(0.1,\dots,3)$

(d) Noise parameters: $(\gamma_{1},\dots,\gamma_{10})=(3,\dots,0.1)$

Figure 9: Ten steps of distribution-centric self-distillation for GPR with
fixed hyperparameters, but with varying choices of noise parameters. We note
that while (c) and (d) yield the same solution at step 10, the path to that
solution varies significantly. Figure 10: Same setup as in Figure 5, but
instead of fitting the distilled model to the soft output of the first, we
also include (in dashed) the ordinary GPC and C-GPC when fitted to the binary
predictions (thresholded at 0.5). We see that the overfitting of the C-GPC is
amplified, while the ordinary GPC closely resembles the initial C-GPC on soft
targets.

##  Appendix C Some details on the continuous Bernoulli distribution

The continuous Bernoulli distribution (denoted by $\mathcal{CB}(\lambda)$) was
introduced in Loaiza-Ganem & Cunningham (2019) and dealt with the relatively
common case, where data is continuous values in $[0,1]$ rather than discrete
values in $\\{0,1\\}$. When performing distillation this is indeed the case,
and merely using a (implicit or explicit) Bernoulli assumption is not
suitable. For good measure, we repeat some key results of $\mathcal{CB}$ from
Loaiza-Ganem & Cunningham (2019) in the following section. Note, in the
following, we will treat Bernoulli variables as if they could attain
continuous values in $[0,1]$, which is the incorrect approach currently often
used in practice.

First, we recall that the density of the Bernoulli distribution with parameter
$\lambda\in(0,1)$ is given by

$\displaystyle
X\sim\mathcal{B}(\lambda)\iff\tilde{p}(x\mid\lambda)=\lambda^{x}(1-\lambda)^{1-x}.$

Second, by Loaiza-Ganem & Cunningham (2019) the density of $\mathcal{CB}$
differs from the density of $\mathcal{B}$ merely by a multiplicative
normalizing constant (dependent on the parameter $\lambda$), denoted by
$C(\lambda)$. See Figure 11 for a plot of both densities. In particular, we
have that for $\lambda\in(0,1)$ then

$\displaystyle X\sim\mathcal{CB}(\lambda)\iff
p(x\mid\lambda)=C(\lambda)\lambda^{x}(1-\lambda)^{1-x}$

where the normalizing constant can be shown to be

$\displaystyle C(\lambda)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\begin{cases}2&\text{if }\lambda=\frac{1}{2},\\\
\frac{2\mathrm{tanh}^{-1}(1-2\lambda)}{1-2\lambda}&\text{otherwise},\end{cases}$

and $\mathrm{tanh}^{-1}$ is the inverse hyperbolic tangent function.

Figure 11: Density of Bernoulli and continuous Bernoulli distributions with
different choices of parameter $\lambda$. Figure 12: Negative log-likelihood
of Bernoulli and continuous Bernoulli distributions with different choices of
parameter $\lambda$. We illustrate the optimizing $\lambda$ for each $x$ by a
cross $\times$, and note that for a given observation $x$ the optimal choice
of $\lambda$ for the Bernoulli is $x$, while it for the continuous Bernoulli
is a non-linear function of $x$, which is generally biased towards more
extreme values of $\lambda$. See also Figure 13.

In Figure 13 we plot the normalizing constant for both the ordinary (equal to
$1$ for all $\lambda$) and continuous Bernoulli distribution, and we also
observe that since the mean of a $\mathcal{CB}(\lambda)$ distributed variable
is not merely $\lambda$, the optimal $\lambda$ (i.e. maximizing the density)
for a single observation $x$ is different from $\lambda$.

Figure 13: Normalizing constant and optimal $\lambda$ values (i.e. maximizing
likelihood for a single observation of $x$) for both ordinary and continuous
Bernoulli distributions.

Finally, in Figure 14 we both plot the log-normalizing constant, as well as
the first and second derivatives wrt. $a$ versus either $\lambda=\sigma(a)$ or
$a$. See Proposition 4.1 for the derivations of these.

Figure 14: Log-normalizing constant for the continuous Bernoulli distribution
$C(\lambda)$, along with first and second derivatives wrt. $a$.

##  Appendix D Time and complexity analysis

Despite the main aim of our proposed self-distillation methods being of a
theoretical nature, in the following, we investigate the compute requirements
of our proposed methods.

###  D.1 (Efficient) Implementation

We provide a public implementation of our methods in python on
github.com/Kennethborup/gaussian_process_self_distillation. We provide four
classes; data-centric and distribution-centric self-distillation for both GP
regression and GP classification. The classes are named straightforwardly as
DataCentricGPR, DistributionCentricGPR, DataCentricGPC, and
DistributionCentricGPC. Each method has .fit, and .predict methods that are
similar to the scikit-learn API. For DataCentricGPR we provide both a naïve
and efficient implementation, which significantly improves the training speed
for multiple steps of distillation. In particular, inspired by Borup &
Andersen (2021) we utilize the singular values decomposition of $\mathbf{K}$
as $\mathbf{V}\mathbf{D}\mathbf{V}^{\intercal}$ (here
$\mathbf{V}\in\mathbb{R}^{N\times N}$ is an orthogonal matrix with the
eigenvectors of $\mathbf{K}$ as rows and $\mathbf{D}\in\mathbb{R}^{N\times N}$
is a non-negative diagonal matrix with the associated eigenvalues in the
diagonal), to rewrite (4) for $t=1$ into

$\displaystyle\mathbf{y}_{1}$
$\displaystyle=\mathbf{V}\mathbf{D}\left(\mathbf{D}+\gamma_{1}\mathbf{I}_{N}\right)^{-1}\mathbf{V}^{\intercal}\mathbf{y}$

which by iterative inserting yields the more general expression

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\left(\prod_{s=1}^{t}\mathbf{K}(\mathbf{K}+\gamma_{s}\mathbf{I}_{N})^{-1}\right)\mathbf{y},$
$\displaystyle=\mathbf{V}\left(\prod_{s=1}^{t}\mathbf{A}_{s}\right)\mathbf{V}^{\intercal}\mathbf{y},$

where
$\mathbf{A}_{s}=\mathbf{D}(\mathbf{D}+\gamma_{s}\mathbf{I}_{N})^{-1}=\operatorname{diag}_{n=1,\dots,.N}\left(\frac{d_{n}}{d_{n}+\gamma_{s}}\right)$
for $s=1,\dots t$ are diagonal matrices. Thus, when $\mathbf{V}$ and
$\mathbf{D}$ are computed, the subsequent distillation steps are
computationally cheap due to the diagonal structure of the $\mathbf{A}_{s}$’s.
One can similarly rewrite the prediction formula (5) - see the source code for
more implementation details.

###  D.2 Time requirements

To evaluate the practical implications of our proposed methods, we compared
the time it takes to fit our methods (for a varying number of distillation
steps) with an ordinary fit using the standard implementation of GP regression
and GP classification in sklearn; namely the GaussianProcessRegressor and
GaussianProcessClassifier. In particular, we fit a model with each of our
methods for a particular number of steps and divide the fitting time with the
fitting time of the corresponding ordinary sklearn method. We repeat our
experiments 30 times and report the mean and 10 and 90 empirical quantiles in
Figure 15. All experiments are performed on a Mac M1 Pro CPU.

For regression, we observe that distribution-centric self-distillation (as
expected from Theorem 3.3) requires no more computation than ordinary GPR with
the sklearn implementation of GP regression. This is simply due to
distribution-centric self-distillation for regression stays within the GP
setting. For data-centric self-distillation the story is different; the time
complexity of a naïve iterative implementation where each successive model is
fitted to the predictions of the preceding model scales linearly with the
number of distillation steps (at a rate of $\approx 0.11$). However, by
utilization of an eigendecomposition of $\mathbf{K}$ in (4), one can rewrite
the solution to optimize the fitting speed significantly. In particular, each
additional step of distillation merely requires updating a diagonal matrix and
a few matrix products (of re-used matrices), which is much cheaper. Thus, with
this efficient implementation, we are able to fit data-centric self-
distillation to any number of distillation steps at a constant time
complexity. Thus, both data-centric and distribution-centric self-distillation
for regression does not take longer to fit than an ordinary GP regression from
sklearn.

For classification, we observe that data-centric self-distillation scales
linearly with the number of distillation steps at a rate of $\approx 0.26$,
and multiple steps of self-distillation in this setting are quickly
computationally costly. However, the time to fit with distribution-centric
self-distillation is (as expected) constant with the number of distillation
steps. This clearly expectable from Proposition 4.2, and fits as quickly as
the sklearn implementation.

Figure 15: Relative time to fit various self-distilled GP models compared to a
single fit with the standard sklearn implementation of the regression and
classification model respectively. On the (left) we plot the relative time to
fit regression models. For data-centric self-distillation, we use both a naïve
and efficient approach. While the naïve approach scales linearly with the
distillation steps, the efficient approach is constant in relative time for
all distillation steps. On the (right) we plot the relative time for the
classification models compared to the sklearn implementation and observe that
the data-centric approach scales linearly with distillation steps, while (as
expected) the distribution-centric approach is constant with distillation
steps. The shaded regions represent the 10 and 90 empirical quantiles of 30
repetitions for each method and the number of distillation steps.

###  D.3 Computational requirements

In the following, we consider the computational requirements of our methods.
We will consider the inversion of matrices, matrix products, Cholesky
decomposition, and SVD to all be $\mathcal{O}(N^{3})$ ignoring constants,
although some algorithms exist to reduce the order of computation for
different methods to $\mathcal{O}(N^{\omega})$ where typically $2<\omega\leq
3$. Our implementation is based on Algorithm 3.1, 3.2, and 5.1 in Rasmussen &
Williams (2006), which largely depends on using the Cholesky decomposition for
back substitution.

We note that while the framing of our setup differs notably from the ordinary
GP regression and GP classification setups, our proposed methods still stay
close to or within the ordinary GP framework simplifying the analysis. In
particular, for distribution-centric GP regression, the solution is an
ordinary GPR model with a particularly scaled noise parameter (see Theorem
3.3), and any existing GPR implementation can be used. Thus, the complexity of
this method is $\mathcal{O}(N^{3})$ as usual. For distribution-centric GP
classification, we consider the approximate method from Proposition 4.2, and
note, that this too is equivalent to an ordinary GPC model, and is in
$\mathcal{O}(N^{3})$. For data-centric GP regression, a naive implementation
would require $t$ iterative fits of an ordinary GPR model, but utilizing the
efficient implementation above, we replace the Cholesky decomposition with the
SVD, and can thus reuse $\mathbf{V}$, and $\mathbf{D}$ for multiple steps of
distillation. Thus, each additional distillation step merely requires matrix
products of diagonal matrices and with the $N\times N$ matrices, which makes
this method $\mathcal{O}(N^{3})$ as well. Finally, due to the iterative nature
of data-centric GP classification, which does not allow a simplified solution,
this method, unfortunately, scales linearly with the number of distillation
steps (although still in $\mathcal{O}(N^{3})$). In particular, for $t>1$ we
change the likelihood function to the continuous Bernoulli, which is a
negligent change due to the simple and closed-form solutions provided in
Proposition 4.1, but since we do not have a closed form solution for our
model, we need to fit each model in an ordinary fashion, thus requiring
$\mathcal{O}^{(}tN^{3})$ with $t$ the number of distillation steps. This can
be computationally demanding for large $t$. Generally, we observe that,
despite data-centric GP classification scaling linearly with $t$, all our
proposed methods are still primarily limited by the scaling of $N$. Finally,
due to the minor changes in our methods, methods to e.g. efficiently
compute/approximate $\mathbf{K}$ should largely still be compatible with our
methods.