# How to Allocate your Label Budget? Choosing between
Active Learning and Learning to Reject in Anomaly Detection

Lorenzo Perini,1 Daniele Giannuzzi, Jesse Davis 1

###### Abstract

Anomaly detection attempts at finding examples that deviate from the expected
behaviour. Usually, anomaly detection is tackled from an unsupervised
perspective because anomalous labels are rare and difficult to acquire.
However, the lack of labels makes the anomaly detector have high uncertainty
in some regions, which usually results in poor predictive performance or low
user trust in the predictions. One can reduce such uncertainty by collecting
specific labels using Active Learning (AL), which targets examples close to
the detector’s decision boundary. Alternatively, one can increase the user
trust by allowing the detector to abstain from making highly uncertain
predictions, which is called Learning to Reject (LR). One way to do this is by
thresholding the detector’s uncertainty based on where its performance is low,
which requires labels to be evaluated. Although both AL and LR need labels,
they work with different types of labels: AL seeks strategic labels, which are
evidently biased, while LR requires i.i.d. labels to evaluate the detector’s
performance and set the rejection threshold. Because one usually has a unique
label budget, deciding how to optimally allocate it is challenging. In this
paper, we propose a mixed strategy that, given a budget of labels, decides in
multiple rounds whether to use the budget to collect AL labels or LR labels.
The strategy is based on a reward function that measures the expected gain
when allocating the budget to either side. We evaluate our strategy on 18
benchmark datasets and compare it to some baselines.

## Introduction

Anomaly detection is the task of automatically detecting examples that do not
follow expected patterns (Chandola, Banerjee, and Kumar 2009). These examples,
named anomalies, are usually indicative of critical events such as water leaks
in stores (Perini, Vercruyssen, and Davis 2022), breakdowns in gas turbines
(Zhao, Wen, and Li 2016), or failures in the petroleum extraction (Martí et
al. 2015). Such critical events usually come along with elevated (maintenance)
costs or with substantial natural damages (e.g., dispersion of petroleum or
gas). Thus, detecting anomalies in time is a relevant task that limits such
resource waste.

Collecting labels, especially for anomalies, is often a hard task because
anomalies are costly events (e.g., machine failures cannot be voluntarily
induced), or simply time-consuming (e.g., you may need to label 100s of
examples before getting an anomaly). Thus, anomaly detection is often tackled
from an unsupervised perspective. However, the lack of labels usually forces
the unsupervised detector to have high uncertainty on specific regions of the
example space (Perini, Vercruyssen, and Davis 2020). High uncertainty is
undesirable because it is often associated with poor predictive performance or
reduced trust in the predictions.

This uncertainty can be tackled in two complementary ways. On the one hand,
one can try to learn a more accurate detector by acquiring a limited number of
labels using Active Learning (AL) (Abe, Zadrozny, and Langford 2006). On the
other hand, it is possible to increase the user trust in the detector’s
outputs by allowing the detector to abstain from making a prediction when it
is highly uncertain, which is called Learning to Reject (LR) (Hendrickx et al.
2021; De Stefano, Sansone, and Vento 2000). One way to do this is to set a
rejection threshold on the detector’s uncertainty based on where its
performance is poor (Cortes, DeSalvo, and Mohri 2016). However, evaluating the
detector performance requires labels.

Both of these approaches rely on labeled data. However, the types of labels
needed for each approach are quite different. Many AL strategies rely on
biased sampling strategies such as explicitly targeting acquiring labels, for
example, for which the detector is highly uncertain (i.e., near the detector’s
current decision boundary) as these are known to yield better performance
(Pimentel et al. 2020; Culotta and McCallum 2005). Alas, using such labels to
evaluate the detector’s performance, as required when setting the threshold in
LR, will yield a biased performance estimate and hence a sub-optimal threshold
(Marrocco, Molinara, and Tortorella 2007). Thus, if a user has a fixed budget
for acquiring labels there is a tension between collecting (a) strategic
labels that can be used to train a better detector, or (b) i.i.d. labels that
can be used to evaluate performance and set a proper rejection threshold.
Therefore, a data scientist is confronted with the challenging question of how
they should optimally allocate their label budget between these two purposes.

In this paper, we assume that the label budget can be split and allocated in
multiple rounds. We introduce ballad (Budget allocation for Active Learning
and Learning to reject in Anomaly Detection) a novel adaptive strategy that,
in each allocation round, (1) measures the potential reward obtained by
assigning the budget to either AL or LR, and (2) chooses the highest reward
option to collect the labels.

## Preliminaries and Related Work

#### Anomaly Detection.

Let $X$ be a $d-$dimensional random variable with unknown $p(X)$. We are given
a dataset $D=\\{x_{1},\dots,x_{n}\\}$ with $n$ examples and $d$ features is
drawn i.i.d. from $p(X)$. Let $V=\\{x_{n+1},\dots,x_{m}\\}\sim_{i.i.d}p(X)$,
$m>n$, be a validation set. Let $Y$ be the label random variable, such that
$Y|X=x$ indicates the class label ($1$ if anomaly, $0$ if normal) for
$x\in\mathbb{R}^{d}$. An anomaly detection problem is the task of finding an
anomaly score function $h\colon\mathbb{R}^{d}\to\mathbb{R}$ and a threshold
$t\in\mathbb{R}$ such that $Y=h^{(t)}(X)$, where $h^{(t)}(x)=1$ if $h(x)\geq
t$, $0$ otherwise. Usually, one sets $t$ based on the contamination factor
$\gamma$, i.e. the proportion of anomalies (Perini, Buerkner, and Klami 2022).

#### Pool-based Active Learning (AL).

The goal of pool-based AL strategies is to reduce the detector’s uncertainty
by selecting the most informative training instances. The AL approaches can be
classified into $3$ categories (Monarch 2021): _uncertainty-based_ sampling
strategies aim to select the unlabeled data samples with the highest
uncertainty (Hacohen, Dekel, and Weinshall 2022), _diversity strategies_
capture the diversity among the training data (Abe, Zadrozny, and Langford
2006; Dagan and Engelson 1995), _combined strategies_ integrate the advantages
of uncertainty-based and diversity-based criteria (Ebert, Fritz, and Schiele
2012).

#### Learning to Reject (LR).

The goal of a detector’s reject option is to abstain from making a prediction
when a detector is too uncertain about predicting a test example (Hendrickx et
al. 2021; Cortes, DeSalvo, and Mohri 2016). Our goal is to develop a detector-
agnostic strategy that does _ambiguity rejection_ , as novelty rejection would
reject all anomalies. Thus, we use a _dependent_ rejector architecture (Chow
1970). We indicate by $\mathcal{C}(x)$ the detector’s confidence for
predicting $x\in V$, and with $\tau\in[0,1]$ the rejection threshold. If the
confidence is below $\tau$, the prediction is rejected
$h^{t}(x)=\text{\textregistered}$. _Note that for appropriate inference, we
need to collect validation labels randomly (i.i.d.)_.

## A strategy to allocate the label budget

This paper tackles the following problem:

Given:

initially unlabeled training set $D$ and validation set $V$, the dataset’s
contamination factor $\gamma$, an anomaly detector $h$, and a label budget
$B$;

Do:

decide whether, in each allocation round $k$, to acquire labels for $D$ (AL)
or for $V$ (LR).

Both training the detector with more labels (AL) and learning a threshold
using larger validation data (LR) improve the detector’s performance. However,
choosing the side to maximize such improvement is challenging for multiple
reasons. First, it requires measuring the reward of either side, i.e. the
expected gain in terms of the detector’s improvement. Second, the rewards need
to be on a similar scale such that neither side is privileged during the
process. Third, comparing a standard detector to one with the reject option is
challenging because the latter needs ad-hoc metrics to overcome the problem of
predicting three classes (anomaly, normal, reject) (Nadeem, Zucker, and
Hanczar 2009).

In this paper, we introduce ballad, a strategy that measures the reward of
allocating the budget for AL, i.e. collecting _strategic_ labels on the
training set, and for LR, i.e., collecting _random_ labels on the validation
set. Let $B=k\cdot b\in\mathbb{N}$ be our labelling budget. We perform $k$
rounds and the labels of $b$ examples are queried in each round. We initialize
the problem by (1) training the detector with no labels and setting a default
rejection threshold, and (2) collecting $b$ random labels for $V$ (LR) and for
$D$ (AL) for a total of $2b$ labels. This allows us to compute the initial
rewards by measuring how the detector varies from (1) to (2): for LR, we
measure the variation after re-setting the validation threshold; for AL, we
measure the variation after re-training the detector with the new labels.
Then, we start the allocation loop. In each round, we allocate the budget to
the option (LR or AL) with the highest reward, and we update the reward using
the new labels. We propose two alternative reward functions: 1) the _entropy_
reward looks at the detector’s probabilities, either for prediction (AL), or
for rejection (LR); 2) the _cosine_ reward considers the predicted class
labels, either anomaly yes/no (AL), or rejection yes/no (LR).

### Measuring the reward

Because we do not know how beneficial the next label allocation would be for
the detector, we look at the past and measure the effect of the last
allocation round. Our challenge is to design a reward function that reflects
the gain when querying the labels. We use the following methods to derive the
reward for both AL and LR, by using as detector’s probabilities either the
probability of predicting anomaly (AL), or the probability of rejecting the
example (LR). Similarly to Vercruyssen et al. (2022), we consider two
scenarios:

#### Entropy.

Adding more labels has the ability to decrease the overall uncertainty of the
anomaly detector. Thus, we measure the variation of the detector’s
probabilities as:

$~{}\mathcal{R}^{e}(k)=\mathbb{E}_{x\sim
X}\left[\left|\mathcal{H}(h_{k}(x))-\mathcal{H}(h_{k-1}(x))\right|\right],$
(1)

where $\mathcal{H}(h(x))=-p\log_{2}p$ is the entropy of the detector’s
probabilities $p$, and the subscript indicates the query round (for $k>2$). A
large difference in entropy means a large detector variation, which indicates
a large impact of the new labels and, in turn, a large reward
$\mathcal{R}^{e}$.

#### Cosine.

More directly, one can measure the impact of the labels in terms of variation
of class predictions. Given the detector’s probabilities, we threshold them at
$0.5$ and assign value $1$ to higher probabilities and $0$ to lower ones.
Thus, we measure the cosine similarity between different outputs as

$\mathcal{R}^{c}(k)=\mathbb{E}_{D\sim X}\left[1-\frac{h_{k}(D)\cdot
h_{k-1}(D)}{\|h_{k}(D)\|\cdot\|h_{k-1}(D)\|}\right],$ (2)

where $h(D)$ is a vector containing the outputs ($0$ or $1$) by the detector
$h$, and $\|\cdot\|$ is the Euclidean norm. This metric is less sensitive to
little variations in the detector and discriminates more in case the new
labels change the predicted class.

### Deriving the detector’s probabilities

Measuring the reward needs some probabilities, which are not easy to derive
due to the partially supervised setting. For both prediction and rejection, we
exploit the squashing function: given a positive real score
$s\in\mathbb{R}_{+}$ and a threshold $\lambda\in\mathbb{R}_{+}$, the squashing
function

$\mathcal{S}_{\lambda}\colon\mathbb{R}_{+}\to(0,1),\quad\mathcal{S}_{\lambda}(s)=1-2^{-\frac{s^{2}}{\lambda^{2}}}$

maps $s$ to a probability $>0.5$ if $s>\lambda$, and $\leq 0.5$ otherwise.
Roughly speaking, $\mathcal{S}_{\lambda}$ calibrates the probabilities by
centering $\lambda$ as the decision threshold.

#### Detector’s posterior probabilities.

Given the contamination factor $\gamma$, a common approach to set the
threshold $t$ is by forcing the detector to have a training positive class
rate equal to $\gamma$. Thus, one can center the probabilities to $t$ by
transforming the anomaly scores $h(x)$ through the squashing function:

$\mathbb{P}(h^{t}(x)=1)=\mathcal{S}_{t}(h(x))\ \text{ s. t. }\
t=Q_{h}(1-\gamma).$

We set $t$ as the $1-\gamma$th quantile of the score distribution such that
only a proportion of $\gamma$ scores have $\mathbb{P}(h^{t}(x)=1)\geq 0.5$.

#### Rejection probabilities.

Given a validation set with some labels, we (1) set a specific detector
confidence $\mathcal{C}(x)$, and (2) set the rejection threshold
$\tau\in[0,1]$. For the former, we use the detector’s posterior probabilities:

$\mathcal{C}(x)=2\times\left|\mathbb{P}(h^{t}(x)=1)-0.5\right|\in[0,1].$

Thus, the closer $\mathbb{P}(h^{t}(x)=1)$ is to $0.5$ (high uncertainty), the
lower the detector confidence. For the latter, we optimize the threshold
$\tau$ over the validation set (only the labeled examples) by minimizing a
cost function $\mathcal{M}(h^{t})$. Finally, we compute the rejection
probabilities by centering 1-confidence values to the rejection threshold,
i.e. by applying the squashing function

$\mathbb{P}(h^{t}(x)=\text{\textregistered})=\mathcal{S}_{\tau}(1-\mathcal{C}(x)).$

### The cost-based evaluation metric

Given a detector with a reject option and a detector without it, we cannot
compare their performance on the non-rejected examples, as they would have
different test sets. Thus, we introduce a cost-based evaluation metric.
Formally, given a rejection cost $c_{r}>0$, a false positive cost $c_{fp}>0$,
and a false negative cost $c_{fn}>0$, the detector is evaluated as:

$\begin{split}\mathcal{M}_{h}&=c_{r}\cdot\mathbb{P}(h^{t}(X)=\text{\textregistered})+c_{fp}\cdot\mathbb{P}(h^{t}(X)=1|Y=0)\\\
&+c_{fn}\cdot\mathbb{P}(h^{t}(X)=0|Y=1).\end{split}$

Note that we assume cost null for the correct predictions, while every
misprediction as well as the rejection gets penalized. Because rejecting is
assumed to be less costly than mispredicting, the rejection cost needs to
satisfy the inequality $c_{r}\leq
min\\{c_{fp}\times(1-\gamma),c_{fn}\times\gamma\\}$, otherwise one could
predict either always normal and pay an expected cost of $c_{fn}\times\gamma$,
or always anomaly and pay $c_{fp}\times(1-\gamma)$.

## Experiments

We experimentally answer the following questions:

  * Q1.

Does ballad result in lower costs when compared to using only AL or LR?

  * Q2.

Which reward metric is better?

  * Q3.

Is the reward function on a similar scale for AL and LR?

  * Q4.

How does our strategy behave when varying $c_{fp}$, $c_{fn}$ ?

### Experimental setup

#### Methods.

We compare ballad111Code available at https://github.com/Lorenzo-Perini/Ballad
to two baselines: All-in-AL allocates all the budget for active learning and
sets the rejection threshold using the (biased) training labels; on the
contrary, All-in-LR allocates all the budget for learning to reject and uses
an unlabeled training set.

#### Data.

Table 1: Properties of the $18$ datasets used. Dataset | # Examples | # Features | $\gamma$
---|---|---|---
ALOI | 12384 | 27 | 0.0304
Annthyroid | 7129 | 21 | 0.0749
Arrhythmia | 271 | 259 | 0.0996
Cardiotocography | 1734 | 21 | 0.0496
Glass | 214 | 7 | 0.0421
InternetAds | 1682 | 1555 | 0.0499
KDDCup99 | 48113 | 40 | 0.0042
PageBlocks | 5473 | 10 | 0.1023
PenDigits | 9868 | 16 | 0.0020
Pima | 526 | 8 | 0.0494
Shuttle | 1013 | 9 | 0.0128
SpamBase | 2661 | 57 | 0.0499
Stamps | 340 | 9 | 0.0912
WBC | 223 | 9 | 0.0448
WDBC | 367 | 30 | 0.0272
WPBC | 160 | 33 | 0.0562
Waveform | 3443 | 21 | 0.0290
Wilt | 4655 | 5 | 0.0199

We carry out our study on $18$ publicly available benchmark datasets, which
are widely used in the literature (Campos et al. 2016). See Table 1 for the
properties.

#### Setup.

For each of the 18 benchmark datasets, we go as follows: (i) we split the
dataset into training, validation and test sets using the proportions
$40-40-20$ (we have a large validation set to better measure the impact of
rejection); (ii) we fit the anomaly detector on the unlabeled dataset and set
the rejection threshold to the default value of $0.1$; (ii) we allocate a
budget $b$ to LR and AL by randomly selecting the initial examples; (iii) we
optimize the rejection threshold and measure the LR reward; (iv) we train the
anomaly detector on the partially labeled training set and measure the AL
reward; (v) we allocate the next round budget $b$ to the option with the
highest reward and repeat (iii) or (iv) until the whole budget $B$ is used.
During each of the steps, we measure the detector performance on the test set
using our cost function. We set $B$ to the $30\%$ of the training set’s size,
and $b$ to $2\%$ of it, such that we run $15$ allocation rounds. We repeat (i
- v) $10$ times and report the average results. In total we run $18\times
15\times 10=2700$ experiments.

#### Costs and hyperparameters.

We set $c_{fp}=c_{fn}=1$ and $c_{r}=\gamma$, following the cost inequality.
ssdo (Vercruyssen et al. 2018) with its default parameters is used as the
semi-supervised anomaly detector (Soenen et al. 2021). We use IForest (Liu,
Ting, and Zhou 2008) as its unsupervised prior. We use Uncertainty Sampling as
the active learning strategy (Zhan et al. 2021), and the entropy as default
reward. For setting the rejection threshold, we use Bayesian Optimization
(gp$\\_$minimize implemented in Skopt) with $20$ calls (Frazier 2018) and
limit the rejection rate on the validation set to $50\%$.

### Experimental results

Figure 1: Comparison between ballad and the All-in strategies on the 18
benchmarks. The x-axis reports the $15$ rounds of $2\%$ labels each. The
y-axis shows the average cost per test example. ballad obtains lower costs in
the majority of cases.

Q1. Comparing ballad to the All-in strategies. Figure 1 shows the comparison
between ballad with the entropy reward and the All-in strategies on the $18$
benchmark datasets. On $8$ datasets (Arrhythmia, Glass, KddCup99, Pima,
SpamBase, Wbc, Wdbc, Wpbc), ballad results in evident lower costs, although
sometimes the difference is small. On $5$ datasets (Cardiotocography,
InternetAds, PageBlocks, Stamps, Waveform) ballad performs similar/worse than
All-in-AL. This happens because ssdo has an overall high performance and a
contained uncertainty in the predictions. On the other hand, in $3$ cases
(Aloi, Annthyroid, Wilt), allocating all the budget for LR has a lower cost.
This is due to the detector being inaccurate and unable to learn from the
training labels, which makes learning an optimal threshold more convenient. As
support for this intuition, we analyze the plain test AUC of SSDO on the whole
test set (no rejection) for each of the three previous cases. By aggregating
over the rounds, SSDO obtains an average AUC equal to $0.86$, $0.88$, and
$0.57$ when the best strategy is, respectively, ballad, All-in-AL, and All-in-
LR. Finally, ballad obtains an overall average cost of $0.043$, which is
$\approx 20\%$ lower than the baselines’ average cost ($0.055$ for All-in-AL,
$0.054$ for All-in-LR).

Budget | Entropy $\mathcal{R}^{e}$ | Cosine $\mathcal{R}^{c}$
---|---|---
2% | 0.0536 $\pm$ 0.0401 | 0.0399 $\pm$ 0.0303
4% | 0.0465 $\pm$ 0.0330 | 0.0411 $\pm$ 0.0334
6% | 0.0443 $\pm$ 0.0284 | 0.0398 $\pm$ 0.0321
8% | 0.0436 $\pm$ 0.0303 | 0.0399 $\pm$ 0.0306
10% | 0.0420 $\pm$ 0.0299 | 0.0411 $\pm$ 0.0325
12% | 0.0416 $\pm$ 0.0303 | 0.0433 $\pm$ 0.0347
14% | 0.0413 $\pm$ 0.0306 | 0.0448 $\pm$ 0.0367
16% | 0.0408 $\pm$ 0.0288 | 0.0456 $\pm$ 0.0372
18% | 0.0403 $\pm$ 0.0290 | 0.0457 $\pm$ 0.0355
20% | 0.0412 $\pm$ 0.0297 | 0.0451 $\pm$ 0.0363
22% | 0.0407 $\pm$ 0.0301 | 0.0451 $\pm$ 0.0359
24% | 0.0421 $\pm$ 0.0325 | 0.0438 $\pm$ 0.0361
26% | 0.0417 $\pm$ 0.0345 | 0.0438 $\pm$ 0.0363
28% | 0.0416 $\pm$ 0.0332 | 0.0438 $\pm$ 0.0380
30% | 0.0418 $\pm$ 0.0345 | 0.0427 $\pm$ 0.0354
Table 2: Average ($\pm$ std) cost per test example over the datasets grouped
by allocation round for each of the two reward functions. For low budgets, the
cosine reward obtains lower costs, while not being competitive for high
budgets.

Q2. Which reward function works better? We analyze both types of reward
functions that we introduced in Eq. 1 and Eq. 2. Table 2 shows the mean and
standard deviation of the cost, divided by allocation round. Overall, using
the cosine reward builds a strategy that produces on average low costs for
little budget ($\leq 10\%$), whereas, for a higher budget, the entropy reward
obtains better average costs. This is due to the highly imbalanced choices
made by the cosine reward: the strategy opts for AL in $93\%$ of the cases,
which usually improves a lot the detector’s performance with few labels but
tends to produce little effect when enough labels are given. On the other
hand, the entropy reward is more balanced and opts for AL in $63\%$ of the
cases. This allows the detector to keep decreasing the costs while learning
during the allocation rounds and obtain more steady performance.

Q3. Is the entropy reward balanced for AL and LR?

Figure 2: Distribution of the difference between AL’s and LR’s entropy reward.
The median close to $0$ indicates the absence of a predominant strategy.

Figure 2 shows the distribution of the difference between AL and LR entropy
rewards over all the $2700$ experiments. Negative values indicate that the LR
reward is higher than AL’s one, while the opposite holds for positive values.
Overall, the median is close to $0$, which means that there is no clearly
predominant strategy. Because the left tail of the density is larger than the
right one, we conclude that LR rewards have higher variability (std $=0.07$ vs
$0.03$).

Q4. The impact of varying $c_{fp}$, and $c_{fn}$. In this experiment, we
penalize more false positives and false negatives by setting, one at a time,
$c_{fp}$ and $c_{fn}$ to $10$. We compare ballad to the two All-in baselines.
For $c_{fp}=10$, our strategy is still the best for low budgets ($<15\%$),
reducing the relative cost by between $5\%$ and $25\%$ with respect to the
runner-up All-in-LR. However, for higher budgets ($>15\%$), All-in-LR becomes
the best strategy as it reduces ballad’s cost by around $20\%$ and All-in-AL’s
cost by more than $40\%$. This happens because the anomaly detector produces
too many false positives, which, if rejected, allow us to reduce the cost. For
$c_{fn}=10$, ballad performs much better than the baselines, reducing their
cost by around $20\%$ (vs All-in-LR) and $24\%$ (vs All-in-AL).

## Conclusion

We proposed ballad, a novel strategy to decide whether to allocate the budget
for Active Learning (AL), i.e. labeling strategic training instances, or for
Learning to Reject (LR), i.e. labeling a random validation set. Our key
insight is that we can measure the expected reward when labeling either set
and allocate the label in the next round to the option with the highest
reward. We proposed two reward functions (entropy and cosine similarity
based). Experimentally, we evaluated ballad on $18$ datasets, and show that it
performs better than simply allocating all the labels to either AL or LR.

#### Acknowledgements.

This work was presented at the 1st AAAI Workshop on Uncertainty Reasoning and
Quantification in Decision Making (UDM23).
This research is supported by an FB Ph.D. fellowship by FWO-Vlaanderen (grant
1166222N) [LP], the Flemish Government under the “Onderzoeksprogramma
Artificiële Intelligentie (AI) Vlaanderen” programme [JD], and KUL Research
Fund iBOF/21/075 [JD].

## References

  * Abe, Zadrozny, and Langford (2006) Abe, N.; Zadrozny, B.; and Langford, J. 2006.  Outlier detection by active learning.  In _Proceedings of ACM SIGKDD_.
  * Campos et al. (2016) Campos, G. O.; Zimek, A.; Sander, J.; Campello, R. J.; Micenková, B.; Schubert, E.; Assent, I.; and Houle, M. E. 2016.  On the evaluation of unsupervised outlier detection: measures, datasets, and an empirical study.  _Data mining and knowledge discovery_.
  * Chandola, Banerjee, and Kumar (2009) Chandola, V.; Banerjee, A.; and Kumar, V. 2009.  Anomaly detection: A survey.  _ACM computing surveys (CSUR)_.
  * Chow (1970) Chow, C. 1970.  On optimum recognition error and reject tradeoff.  _IEEE Transactions on information theory_.
  * Cortes, DeSalvo, and Mohri (2016) Cortes, C.; DeSalvo, G.; and Mohri, M. 2016.  Learning with rejection.  In _International Conference on Algorithmic Learning Theory_. Springer.
  * Culotta and McCallum (2005) Culotta, A.; and McCallum, A. 2005.  Reducing labeling effort for structured prediction tasks.  In _AAAI_.
  * Dagan and Engelson (1995) Dagan, I.; and Engelson, S. P. 1995.  Committee-based sampling for training probabilistic classifiers.  In _Machine Learning Proceedings 1995_. Elsevier.
  * De Stefano, Sansone, and Vento (2000) De Stefano, C.; Sansone, C.; and Vento, M. 2000.  To reject or not to reject: that is the question-an answer in case of neural classifiers.  _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_.
  * Ebert, Fritz, and Schiele (2012) Ebert, S.; Fritz, M.; and Schiele, B. 2012.  Ralf: A reinforced active learning formulation for object class recognition.  In _2012 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE.
  * Frazier (2018) Frazier, P. I. 2018.  Bayesian optimization.  In _Recent advances in optimization and modeling of contemporary problems_. Informs.
  * Hacohen, Dekel, and Weinshall (2022) Hacohen, G.; Dekel, A.; and Weinshall, D. 2022.  Active learning on a budget: Opposite strategies suit high and low budgets.  _arXiv preprint arXiv:2202.02794_.
  * Hendrickx et al. (2021) Hendrickx, K.; Perini, L.; Van der Plas, D.; Meert, W.; and Davis, J. 2021.  Machine learning with a reject option: A survey.  _arXiv preprint arXiv:2107.11277_.
  * Liu, Ting, and Zhou (2008) Liu, F. T.; Ting, K. M.; and Zhou, Z.-H. 2008.  Isolation forest.  In _8th IEEE international conference on data mining_.
  * Marrocco, Molinara, and Tortorella (2007) Marrocco, C.; Molinara, M.; and Tortorella, F. 2007.  An empirical comparison of ideal and empirical ROC-based reject rules.  In _International Workshop on Machine Learning and Data Mining in Pattern Recognition_. Springer.
  * Martí et al. (2015) Martí, L.; Sanchez-Pi, N.; Molina, J. M.; and Garcia, A. C. B. 2015.  Anomaly detection based on sensor data in petroleum industry applications.  _Sensors_.
  * Monarch (2021) Monarch, R. M. 2021.  _Human-in-the-Loop Machine Learning: Active learning and annotation for human-centered AI_.  Simon and Schuster.
  * Nadeem, Zucker, and Hanczar (2009) Nadeem, M. S. A.; Zucker, J.-D.; and Hanczar, B. 2009.  Accuracy-rejection curves (ARCs) for comparing classification methods with a reject option.  In _Machine Learning in Systems Biology_. PMLR.
  * Perini, Buerkner, and Klami (2022) Perini, L.; Buerkner, P.; and Klami, A. 2022.  Estimating the Contamination Factor’s Distribution in Unsupervised Anomaly Detection.  _arXiv preprint arXiv:2210.10487_.
  * Perini, Vercruyssen, and Davis (2020) Perini, L.; Vercruyssen, V.; and Davis, J. 2020.  Quantifying the confidence of anomaly detectors in their example-wise predictions.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer.
  * Perini, Vercruyssen, and Davis (2022) Perini, L.; Vercruyssen, V.; and Davis, J. 2022.  Transferring the Contamination Factor between Anomaly Detection Domains by Shape Similarity.  In _Proceedings of the AAAI Conference on Artificial Intelligence_.
  * Pimentel et al. (2020) Pimentel, T.; Monteiro, M.; Veloso, A.; and Ziviani, N. 2020.  Deep active learning for anomaly detection.  In _2020 International Joint Conference on Neural Networks_. IEEE.
  * Soenen et al. (2021) Soenen, J.; Van Wolputte, E.; Perini, L.; Vercruyssen, V.; Meert, W.; Davis, J.; and Blockeel, H. 2021.  The effect of hyperparameter tuning on the comparative evaluation of unsupervised anomaly detection methods.  In _Proceedings of the KDD_.
  * Vercruyssen et al. (2018) Vercruyssen, V.; Meert, W.; Verbruggen, G.; Maes, K.; Baumer, R.; and Davis, J. 2018\.  Semi-supervised anomaly detection with an application to water analytics.  In _2018 IEEE international conference on data mining (icdm)_. IEEE.
  * Vercruyssen et al. (2022) Vercruyssen, V.; Perini, L.; Meert, W.; and Davis, J. 2022.  Multi-domain Active Learning for Semi-supervised Anomaly Detection.  _ECML 2022 published proceedings_.
  * Zhan et al. (2021) Zhan, X.; Liu, H.; Li, Q.; and Chan, A. B. 2021.  A Comparative Survey: Benchmarking for Pool-based Active Learning.  In _IJCAI_.
  * Zhao, Wen, and Li (2016) Zhao, N.; Wen, X.; and Li, S. 2016.  A review on gas turbine anomaly detection for implementing health management.  _Turbo Expo: Power for Land, Sea, and Air_.


# Deep Gaussian Processes: A Survey

Kalvik Jakkala
Department of Computer Science
University of North Carolina at Charlotte
Charlotte  NC  USA
<EMAIL_ADDRESS>
###### Abstract

Gaussian processes are one of the dominant approaches in Bayesian learning.
Although the approach has been applied to numerous problems with great
success, it has a few fundamental limitations. Multiple methods in literature
have addressed these limitations. However, there has not been a comprehensive
survey of the topics as of yet. Most existing surveys focus on only one
particular variant of Gaussian processes and their derivatives. This survey
details the core motivations for using Gaussian processes, their mathematical
formulations, limitations, and research themes that have flourished over the
years to address said limitations.

Furthermore, one particular research area is Deep Gaussian Processes (DGPs),
it has improved substantially in the past decade. The significant publications
that advanced the forefront of this research area are outlined in their
survey. Finally, a brief discussion on open problems and research directions
for future work is presented at the end.

##  1 Introduction

There have been numerous advances in the field of machine learning in recent
years. Most of these advances can be attributed to backpropagation, large
datasets, and computational resource improvements. However, most of the
currently popular machine learning approaches mainly, deep learning methods,
are based on frequentist approaches, which entail making any prediction
decisions by studying the correlations among features and predictions in a
dataset. The problem with such an approach is that it is easy to overfit to
the available dataset and risk learning unwanted biasing in the dataset.

Furthermore, current approaches make it difficult and unintuitive to introduce
any prior domain knowledge into a prediction model. Several real-world
problems have domain experts; incorporating their knowledge could result in
substantially better models. However, most deep learning approaches do not
accommodate such incorporations and require application-specific methods to be
developed to address such an issue.

Prediction uncertainty is an important metric that needs to be estimated by
reliable models. Most data sources contain non-negligible quantities of noise
that might hinder the performance of a prediction model. It is also not
uncommon to have test data samples that do not closely resemble the training
dataset distribution. In such cases, it is essential to know the model’s
prediction uncertainty. If the model were to be used in a mission-critical
task without accounting for its prediction uncertainty, it could lead to
catastrophic results.

Another major drawback of conventional deep learning approaches is model
comparison. Deep learning approaches are parametric and require explicit
definitions for the model architecture. Moreover, model architectures are
application specific. Often multiple model architectures need to be compared
against each other to determine which is the best model for a task. However,
it is nontrivial to factor in model size in terms of its parameter count and
accuracy for comparison.

The limitations mentioned above are addressed by Bayesian approaches with
varying degrees of ease and efficiency. We can incorporate domain knowledge
with a prior distribution, prediction uncertainty can be estimated with
prediction variance, and models can be aptly compared against each other with
the Bayes factor.

Apart from the advantages mentioned above, another interesting feature of
Bayesian methods is that they facilitate causal modeling of any system or
process. Indeed, most classification or regression problems require a chain of
sub-decisions, each of which would lead to the final prediction. However,
conventional deep learning approaches are not particularly amenable to
specifying such causal models. The Bayesian framework, along with do-calculus
[Pearl, 2000, Pearl and Mackenzie, 2018], can be used to specify such
structures in models.

The advantages of Bayesian approaches raise the question of why they dont have
widespread adaptation yet. Bayesian approaches often incur heavy computational
expense or outright intractabilities making them infeasible for several
problems. Nonetheless, these methods are steeped with history and have been
used to solve numerous problems with substantial ramifications [McGrayne,
2011]. Time and time again, the Bayesian framework has proved itself to be
worthy of further research.

This paper considers one particular type of Bayesian approach, i.e., Gaussian
processes [Rasmussen and Williams, 2006]. The method has its roots in
Stochastic processes—a field of study dedicated to modeling random processes
with Probability theory [Klebaner, 2012, Rosenthal, 2006]. Most problems of
interest are usually not deterministic processes, or even if they are, one
might not have access to all the information needed to model it as such.
Stochastic processes mathematically accommodate such uncertainties, and
Gaussian processes are one particular variant of Stochastic processes.

I start my exposition by detailing Gaussian processes, their advantages, and
their disadvantages. However, the main focus of this survey is Deep Gaussian
Processes (DGPs). I will describe some of the Gaussian Processes’ prominent
variants crucial for building DGPs and explain the critical DGP approaches.

##  2 Preliminary

A collection of random variables is called a Gaussian Process(GP) if the joint
distribution of any finite subset of its variables is a Gaussian. Gaussian
processes can be thought of as a reinterpretation or a generalization of
Gaussian distributions. Gaussian processes are essentially Gaussian
distributions with infinite dimensions. Each dimension of the distribution
corresponds to an input data sample, and the distribution mean represents each
data input’s associated label.

Gaussian processes are fascinating because of their analytical properties.
Gaussian processes being Gaussian distributions, have closed-form analytical
solutions to summation, conditioning, and marginalization under certain
conditions [Murphy, 2012]. Given a Gaussian distribution on a dataset
$\mathbf{X}$ with mean $\boldsymbol{\mu}$ and covariance $\mathbf{\Sigma}$,
the marginals and conditionals are given by the following

$\mathbf{X}=\left(\begin{array}[]{l}\mathbf{X}_{1}\\\
\mathbf{X}_{2}\end{array}\right),\quad\boldsymbol{\mu}=\left(\begin{array}[]{l}\boldsymbol{\mu}_{1}\\\
\boldsymbol{\mu}_{2}\end{array}\right),\quad\boldsymbol{\Sigma}=\left(\begin{array}[]{ll}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{array}\right)$
$\displaystyle p\left(\mathbf{X}_{1}\right)$
$\displaystyle=\mathcal{N}\left(\mathbf{X}_{1}\mid\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11}\right)$
$\displaystyle p\left(\mathbf{X}_{2}\right)$
$\displaystyle=\mathcal{N}\left(\mathbf{X}_{2}\mid\boldsymbol{\mu}_{2},\mathbf{\Sigma}_{22}\right)$
$\displaystyle p\left(\mathbf{X}_{1}\mid\mathbf{X}_{2}\right)$
$\displaystyle=\mathcal{N}\left(\mathbf{X}_{1}\mid\boldsymbol{\mu}_{1\mid
2},\boldsymbol{\Sigma}_{1\mid 2}\right)$ $\displaystyle\boldsymbol{\mu}_{1\mid
2}$
$\displaystyle=\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\left(\mathbf{X}_{2}-\boldsymbol{\mu}_{2}\right)$
$\displaystyle\boldsymbol{\Sigma}_{1\mid 2}$
$\displaystyle=\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21}$

Moreover, as mentioned before, deep neural networks require the specification
of model architectures which are usually application-specific. However, GPs do
not suffer from such issues as they are non-parametric. GPs consider the
entire function space and marginalize it to obtain their predictions. The
approach has two main advantages; first, one need not decide the model’s
complexity. The second, the marginalization process, induces Occum’s Razor
[Rasmussen and Williams, 2006], thereby overcoming any overfitting issues.

One might wonder if modeling any arbitrary process as a Gaussian would result
in poor performance. However, Gaussians are omnipresent, and the Central Limit
Theorem applies to them [Rasmussen and Williams, 2006]. The theorem states
that any dataset’s mean follows a Gaussian distribution as the number of
samples from the distribution increases. So, assuming enough data is
available, usually, a Gaussian approximation results in reasonable estimates.

Another attractive property of Gaussians is that the distribution has the
maximum entropy given a dataset’s mean and variance [Murphy, 2012]. The first
two moments, mean and variance, are usually the only metrics that we can
accurately compute from data. As a result, using Gaussians to estimates
distributions results in the most informative data distribution estimates.

Finally, another essential property of Gaussian processes was established by
Neil [Neal, 1996] where he showed that GPs are equivalent to a single-layered
Perceptron with infinite hidden units. Perceptrons were, in turn, known to be
universal approximators that can approximate any function given enough data
[Hornik et al., 1989].

##  3 Gaussian Processes

I have detailed the critical advantages of Bayesian approaches and why
researchers are interested in Gaussian processes specifically. This section
further elaborates on GPs. I give an intuition of GPs; their mathematical
formulation [Rasmussen and Williams, 2006, Murphy, 2012], and an intuitive
interpretation of the terms in its formulation. Furthermore, I will explain
Kernel functions and list some limitations of GPs.

Figure 1: Illustration of a Gaussian Process on the MNIST dataset.

GPs, as mentioned before, are Gaussian distributions interpreted as
probabilistic mappings from an input data space to its output label space.
This mapping has a convenient interpretation in function space. Formally, we
assume a function $f$, maps an input space $\mathbf{X}$ to an output space
$\mathbf{y}$. This mapping is characterized by a noise term $\epsilon$.

$\mathbf{y}=f(\mathbf{X})+\boldsymbol{\epsilon},\boldsymbol{\epsilon}\sim\mathcal{N}\left(\mathbf{0},\beta^{-1}\mathbf{I}\right)$

In a non-probabilistic approach, a parametric form for the function $f$ is
considered, and Maximum A Posteriori (MAP) is used to estimate the function
parameters. GPs, on the other hand, assume a distribution over the functions
$f$. The function values are also called latent variables, and any predictions
in GPs involve marginalizing over the entire function space. Fig. 1
illustrates this process. Furthermore, the distribution over the functions is
parametrized by a mean vector and a covariance matrix. The covariance is
computed with a similarity metric known as a kernel function $K$.

$\displaystyle\mathbf{f}|\mathbf{X}$
$\displaystyle\sim\mathcal{G}\mathcal{P}\left(\mu(\mathbf{X}),K\left(\mathbf{X},\mathbf{X}\right)\right)$
$\displaystyle\mu(\mathbf{X})$
$\displaystyle=\mathbb{E}[\mathbf{f}(\mathbf{X})]$ $\displaystyle
K\left(\mathbf{X},\mathbf{X}\right)$
$\displaystyle=\mathbb{E}\left[(\mathbf{f}(\mathbf{X})-\mu(\mathbf{X}))\left(\mathbf{f}\left(\mathbf{X}\right)-\mu\left(\mathbf{X}\right)\right)\right]$

Moreover, by treating the function outputs as noisy versions of the labels,
the distribution over $\mathbf{y}$ conditioned on the input $\mathbf{X}$ can
be achieved as follows.

$\displaystyle p(\mathbf{y}\mid\mathbf{f})$
$\displaystyle=\mathcal{N}\left(\mathbf{f},\beta^{-1}\mathbf{I}\right)$
$\displaystyle p(\mathbf{y}\mid\mathbf{X})$ $\displaystyle=\int
p(\mathbf{y}\mid\mathbf{f})p(\mathbf{f}\mid\mathbf{X})d\mathbf{f}$

The above marginalization is the likelihood of the dataset. It is used to fit
hyperparameters involved in the GP formulation by maximizing the training
dataset’s likelihood. However, we are also interested in predicting labels for
new data samples. We can compute the predictions by considering a joint
distribution over the function values of data $\mathbf{X}$ and
$\mathbf{X}_{*}$ from the training and testing datasets denoted by
$\mathbf{f}$ and $\mathbf{f}_{*}$ respectively.

$\left[\begin{array}[]{l}\mathbf{f}\\\
\mathbf{f}_{*}\end{array}\right]\sim\mathcal{N}\left(\mathbf{0},\left[\begin{array}[]{ll}K(\mathbf{X},\mathbf{X})&K\left(\mathbf{X},\mathbf{X}_{*}\right)\\\
K\left(\mathbf{X}_{*},\mathbf{X}\right)&K\left(\mathbf{X}_{*},\mathbf{X}_{*}\right)\end{array}\right]\right)$

We are interested in $\mathbf{f}_{*}$ as it is the noise-free realization of
$\mathbf{y}_{*}$. The posterior distribution of $\mathbf{f}_{*}$ can be
obtained by conditioning it with $\mathbf{X}_{*}$, $\mathbf{X}$, and
$\mathbf{y}$ as follows.

$p(\mathbf{f_{*}}\mid\mathbf{X}_{*},\mathbf{X},\mathbf{y})=\frac{1}{p({\mathbf{y}})}\int
p(\mathbf{y}\mid\mathbf{f})p(\mathbf{f},\mathbf{f_{*}})d\mathbf{f}$
$\displaystyle
p(\mathbf{f_{*}}\mid\mathbf{X}_{*},\mathbf{X},\mathbf{y})\sim\mathcal{N}($
$\displaystyle
K\left(\mathbf{X}_{*},\mathbf{X}\right)K(\mathbf{X},\mathbf{X})^{-1}\mathbf{f},$
$\displaystyle\left.K\left(\mathbf{X}_{*},\mathbf{X}_{*}\right)-K\left(\mathbf{X}_{*},\mathbf{X}\right)K(\mathbf{X},\mathbf{X})^{-1}K\left(\mathbf{X},\mathbf{X}_{*}\right)\right)$

The above conditional prediction distribution’s mean is treated as the model
prediction since a Gaussian distribution’s mode coincides with its mean. The
mean term of the predictive distribution can be interpreted as a weighted
average of the training set labels. The weights are defined by the kernel
matrices, which assign higher weights to data points close to the test samples
in feature space.

The variance term has two parts; the first term on the left can be assumed to
be the prior variance from the testing data. And the second term can be
thought of as the variance induced by the training data. By subtracting the
two terms, the initial test data variance is reduced by the training data’s
information.

###  3.1 Kernel Function

Unlike a conventional Gaussian distribution, a GP has infinite dimensions;
therefore, explicitly defining a covariance matrix would be infeasible. GPs
resort to using kernel functions to determine the covariance matrix, allowing
for a compact representation of the kernel matrix.

Kernel functions measure data similarity, and GPs use it to produce a positive
semi-definite (PSD) matrix. Kernel functions can be interpreted as dot
products between data points in a high dimensional manifold without explicitly
mapping them to such a manifold. Hence, these functions allow us to define
meaningful similarity metrics relevant to a given task in a computationally
efficient manner.

Numerous kernel functions have been used in literature, each with its own set
of properties that make it useful for different tasks. Enumerating and
explaining all such functions is out of this paper’s scope, and the readers
are referred to [Rasmussen and Williams, 2006] for a detailed overview.
However, two kernels are often used in literature, particularly the Squared
Exponential kernel and the Automatic Relevance Detection (ARD) kernel
functions detailed below.

The Squared Exponential kernel, also known as the Radial Basis Function and
Gaussian kernel, has an exponential form. It has two parameters, the length
scale $\ell$, which determines how close a pair of data points should be in
order to be correlated. And it has a variance scale $\sigma$, which quantifies
the variance induced from data points. The function is parameter efficient as
it has only two parameters, and its exponential form allows for analytical
computations of integrals in some cases. For a pair of inputs $\mathbf{x}$,
$\mathbf{x}^{{}^{\prime}}\in\mathbb{R}^{q}$, the Gaussian kernel is defined as
follows

$K\left(\mathbf{x},\mathbf{x}^{{}^{\prime}}\right)=\sigma^{2}\exp\left(-\frac{1}{2\ell^{2}}\sum_{j=1}^{q}\left(\mathbf{x}_{j}-\mathbf{x}^{{}^{\prime}}_{j}\right)^{2}\right)$

Another important kernel function is the Automatic Relevance Detection (ARD)
kernel. It has a similar form to the Gaussian kernel. However, it has weight
parameters for each dimension of the input space. We can use the weights to
scale down irrelevant features and increase the importance of the relevant
ones.

$K\left(\mathbf{x},\mathbf{x}^{{}^{\prime}}\right)=\sigma^{2}\exp\left(-\frac{1}{2}\sum_{j=1}^{q}w_{j}\left(\mathbf{x}_{j}-\mathbf{x}^{{}^{\prime}}_{j}\right)^{2}\right)$

The weight parameter in the ARD kernel function can be defined manually.
However, such an approach is infeasible as we usually don’t know in advance if
a feature is essential or not. In practice, the weights are determined from
the training dataset, which allows us to obtain kernel function
parametrizations that are most suited for the task. Indeed, there are two
primary approaches to determine the function parameters, Maximum A Posteriori
(MAP) and the Bayesian approach.

The MAP approach considers the log probability of the dataset $\log
p(\mathbf{y}\mid\mathbf{X})$. The optimal parameters are computed by taking
the log probability derivative with respect to the kernel parameters.
Depending on the dataset size, the parameters can be calculated analytically
but, in some cases, one needs to resort to approaches like gradient descent to
determine the parameters.

$\log
p(\mathbf{f}\mid\mathbf{X})=-\frac{1}{2}\mathbf{f}^{\top}K^{-1}\mathbf{f}-\frac{1}{2}\log|K|-\frac{n}{2}\log
2\pi\\\ $ $\log
p(\mathbf{y}\mid\mathbf{X})=-\frac{1}{2}\mathbf{y}^{\top}\left(K+\sigma_{n}^{2}I\right)^{-1}\mathbf{y}-\frac{1}{2}\log\left|K+\sigma_{n}^{2}I\right|-\frac{n}{2}\log
2\pi$

A limitation of using MAP to determine the kernel parameters is that it could
lead to overfitting. The problem’s severity varies depending on the dataset
size, the GP’s variant, and the kernel function being used. Nonetheless, we
can mitigate the problem by using a Bayesian approach.

A Bayesian approach defines a distribution over the kernel parameters and
computes the posterior parameter distribution conditioned on the training
dataset. Although the Bayesian approach is ideal, it is usually intractable to
analytical computations. As such, one has to resort to sampling approaches
such as Hamiltonian Monte Carlo (HMC) methods [Bishop, 2016] to determine the
parameter distributions. However, HMC methods have their limitations and
introduce additional computational costs that might be infeasible.

###  3.2 Limitations

Although GPs have several advantages, they also have a few key limitations
which hinder their use in most machine learning problems. Particularly, there
are three main problems:

  1. 1.

Computation cost

  2. 2.

Storage cost

  3. 3.

Hierarchical feature extraction

Figure 2: Overview of Gaussian Processes and their variants.

The computational cost of GPs can be rather substantial; one needs to invert
the kernel matrix to get the predictive distribution of a GP. The kernel
matrix’s size is $n\times n$ where $n$ is the number of data points in the
training dataset. Inverting such a matrix takes $O(n^{3})$ computation time.
Furthermore, once the kernel matrix inverse is available, it takes $O(n)$ and
$O(n^{2})$ time to determine the predictive distribution’s mean and variance
of a new data point.

Additionally, since GPs require the entire training dataset’s storage, the
storage cost is $O(n^{2})$. Depending on the size of the dataset, the storage
costs substantially limit the scalability of the approach. Moreover, if the GP
were to be used in an environment where the training dataset size keeps
increasing, the computation and storage costs could overwhelm the entire
process, rendering the benefits of GP far too expensive. As such, GPs are
usually only feasible for datasets with about $1000-3000$ data points.

Another major drawback of GPs is the lack of kernel functions capable of
handling structured data wherein one needs to consider hierarchical feature
extraction to determine the similarity of a pair of data points properly. Such
an issue often arises in data such as images but, it is also prevalent in
simpler vector datasets. Conventional kernel functions are ill-equipped to
handle such correlations, and one needs to resort to deep feature extraction
like the ones used in Deep learning models. However, such feature extraction
still needs to be confined to a Bayesian framework to retain a GP’s
advantages.

Sparse Gaussian Processes address the computational and storage costs. And the
feature extraction issue is addressed by Deep Gaussian Processes. I explain
some of the prominent methods for Sparse and Deep GPs that have been developed
over the last two decades in the following sections. Fig: 2 shows a flow chart
that relates the limitations to the GP varients that address them.

##  4 Sparse Gaussian Processes

Figure 3: Overview of Sparse Gaussian Processes.

Given the computational and storage requirements that hinder the widespread
use of GPs, a substantial number of papers have tried to address the problem
and are collectively referred to as Sparse Gaussian Processes (SGPs), Fig. 3
depicts the prominent methods coverd in this section. The terminology stems
from the way most of these approaches address the issue. Because the primary
problem is the inversion of a covariance matrix, most methods try to introduce
sparsity and reduce the matrix size that needs to be inverted while retaining
the original matrix’s performance.

This section focuses on some of the well-known methods crucial for developing
some of the Deep Gaussian Process methods that will be detailed in the
upcoming section. A complete overview of all SGPs is outside this survey’s
scope; the readers are referred to [Liu et al., 2020] for a thorough summary.

The Nyström approximation of [Williams and Seeger, 2001] is a well know
approach to reducing the inversion cost of the covariance matrix in GPs. The
Nyström approximation allows one to generate a low rank approximation of any
kernel matrix. The approach is applied to GPs by selecting $m$ data points
with $m<<n$ from the training set. Then a low rank approximation $\tilde{K}$
to the kernel matrix is computed as shown below

$\tilde{K}=K_{n,m}K_{m,m}^{-1}K_{m,n}$

Here, $K_{n,m}$ represents the kernel matrix computed from the $n$ and $m$
data points in the training dataset and the selected subset, respectively. The
same notation is used for the other kernel matrices. The approximation needs
the inversion of only an $m\times m$ matrix, thereby reducing the computation
cost from $O(n^{3})$ to $O(m^{3})$.

However, the approximation assumes that the data comes from a low rank
manifold, which would be the case if the data dimension $d<n$. In such a case,
the low rank approximation would be exact and result in no information loss.
But, the selected $m$ data points also influence the approximation. It is
possible to have data points that result in a poor approximation even if the
data is from a low dimensional manifold.

In practice, most datasets have more data points than the number of features;
thus, the method is applicable in most cases. However, selecting the data
points is crucial to the approximation’s performance. Williams and Seeger used
a random subset of $m$ data points in their approach [Williams and Seeger,
2001]. Although the method works, the unsophisticated data selection process
limited the method’s performance.

Snelson and Ghahramani [Snelson and Ghahramani, 2006] addressed the subset
selection for Nyström approximation by treating the subset as model parameters
and termed them pseudo data. The pseudo data is assumed to be synthetic and
does not necessarily correspond to any data point available in the training
dataset. Indeed, they could take values that are some combination of the
training dataset.

The pseudo points were computed using maximum likelihood. However, to use
maximum likelihood, one needs to parameterize the GP with pseudo points
appropriately. Snelson and Ghahramani [Snelson and Ghahramani, 2006]
introduced a distribution over the pseudo points and considered a joint
distribution over latent representations of data from the training, testing
and, pseudo data points given by $\mathbf{f},\mathbf{f_{*}}$ and $\mathbf{m}$.
The authors then marginalized the pseudo points to obtain the posterior
distribution as shown below

$\displaystyle p(\mathbf{y}\mid\mathbf{X},\mathbf{m})$
$\displaystyle=\mathcal{N}\left(\mathbf{y}\mid
K_{nm}K_{mm}^{-1}\mathbf{m},\mathbf{\Lambda}+\sigma^{2}\mathbf{I}\right)$
$\displaystyle\mathbf{\Lambda}$
$\displaystyle=\operatorname{diag}(\boldsymbol{\lambda})$
$\displaystyle\lambda_{n}$
$\displaystyle=K_{nn}-K_{mn}^{\top}K_{mm}^{-1}K_{mn}$ $\displaystyle
p\left(\mathbf{f}_{*},\mathbf{f}\right)$ $\displaystyle=\int
p\left(\mathbf{f}_{*},\mathbf{f},\mathbf{m}\right)\mathrm{d}\mathbf{m}=\int
p\left(\mathbf{f}_{*},\mathbf{f}\mid\mathbf{m}\right)p(\mathbf{m})\mathrm{d}\mathbf{m}$
$\displaystyle p(\mathbf{m})$
$\displaystyle=\mathcal{N}\left(\mathbf{0},K_{\mathbf{m},\mathbf{m}}\right)$

Although using maximum likelihood to determine the pseudo points’ distribution
works in practice, using maximum likelihood runs the risk of overfitting. It
would be ideal to use a Bayesian approach to compute the pseudo points
distribution conditioned with the training set. Unfortunately, such an
approach is infeasible as it becomes intractable to analytical solutions of
the pseudo points. Moreover, the approach works by assuming that the joint
distribution $p\left(\mathbf{f}_{*},\mathbf{f}\right)$ can be partitioned as
follows

$p\left(\mathbf{f}_{*},\mathbf{f}\right)=\int
p\left(\mathbf{f}_{*},\mid\mathbf{m}\right)p\left(\mathbf{f}\mid\mathbf{m}\right)p(\mathbf{m})\mathrm{d}\mathbf{m}\\\
$

The assumption limits the information a GP obtains from the training set to be
induced only through the pseudo set. Hence, the pseudo points are also
referred to as inducing inputs. The factorization assumption limits the
capacity of the model and affects the accuracy of the model. Notably, the
prior distribution that is assumed for the pseudo set substantially influences
the results.

Snelson and Ghahramani [Snelson and Ghahramani, 2006] treated the pseudo
points as hyperparameters and introduced a prior over them, resulting in an
inaccurate posterior compared to a vanilla GP. The inaccuracy was a
consequence of the formulation of the kernel approximation. Titsias [Titsias,
2009] addressed the outfitting and inexact posterior issue by considering a
variational approach. The approach introduces a lower bound that can be
optimized to determine the inducing inputs and the kernel hyperparameters. The
bound shown below, can be used to solve for the inducing points and kernel
hyperparameters. We can then use the inducing points for computing the
predictive distribution.

$\displaystyle\mathcal{F}$
$\displaystyle=\log\left[\mathcal{N}\left(\mathbf{y}\mid\mathbf{0},\sigma^{2}I+Q_{nn}\right)\right]-\frac{1}{2\sigma^{2}}\operatorname{Tr}(\widetilde{K})$
$\displaystyle Q_{nn}$ $\displaystyle=K_{nm}K^{-1}_{mm}K_{mn}$
$\displaystyle\widetilde{K}$ $\displaystyle=K_{nn}-K_{nm}K^{-1}_{mm}K_{mn}$

However, the marginal in [Titsias, 2009] didn’t have the factorization
required to apply stochastic gradient descent. Hensman et al. [Hensman et al.,
2013] improved on the work of Titsias [Titsias, 2009] by developing a new
bound that could be optimized with stochastic gradient descent. Unlike
Titsias’s approach, the method in [Hensman et al., 2013] did not require the
entire dataset at once to compute the variational parameters. It used the
following bound that could be optimized with stochastic gradient descent.

$\displaystyle\mathcal{F}=\sum_{i=1}^{n}\\{$
$\displaystyle\log\mathcal{N}\left(\mathbf{y}_{i}\mid
K_{i}^{\top}K_{mm}^{-1}\mathbf{v},\beta^{-1}\right)$
$\displaystyle\left.-\frac{1}{2}\beta\widetilde{k}_{i,i}-\frac{1}{2}\operatorname{tr}\left(\mathbf{S}\mathbf{\Lambda}_{i}\right)\right\\}$
$\displaystyle-\mathrm{KL}(q(\mathbf{u})\|p(\mathbf{u}))$

$\displaystyle q(\mathbf{u})$
$\displaystyle=\mathcal{N}(\mathbf{u}\mid\mathbf{v},\mathbf{S})$
$\displaystyle\mathbf{\Lambda}$ $\displaystyle=\beta
K^{-1}_{mm}k_{i}k_{i}^{\top}K^{-1}_{mm}$

Here, $\mathbf{u}$ is the set of feature space representations of the inducing
points, an0d $k_{i}$ is the $i^{th}$ column of $K_{mn}$. Hensman et al. showed
the approach to scale well to large datasets while retaining the reduced model
complexity of $O(m^{3})$.

##  5 Gaussian Process Latent Variable Model

The methods discussed so far primarily addressed the computation and storage
cost issues. This section introduces a variant of GPs that can be trained in
an unsupervised manner. Gaussian Process Latent Variable Models (GPLVMs)
[Lawrence, 2004, Lawrence, 2005] assume the feature space is a latent space
with unknown data distribution. The latent space distribution in then learned
during the training phase. Although the method seems irrelevant, it plays a
significant role in some of the Deep Gaussian Processes that I will cover in
the next section.

Lawrence [Lawrence, 2004, Lawrence, 2005] showed that if the function space of
a GP is constrained to linear function space, the GP can be interpreted as a
probabilistic variant of Principal Component Analysis (PCA). Moreover, if the
function space is relaxed to a non-linear space defined by a kernel function,
it can be interpreted as probabilistic non-linear PCA. The approach assumes a
standard Gaussian prior over the input space and maximizes the log probability
of the dataset $p(\mathbf{y}\mid\mathbf{X},\beta)$ with respect to the input
data likelihood $\mathbf{X}$.

The input data or latent space distribution can not be computed analytically
because of the non-linearities introduced by the kernel function. However,
Lawrence [Lawrence, 2004, Lawrence, 2005] showed that the distribution could
be estimated using the Expectation-Maximization algorithm. But, the approach
returns only the mode of the data distribution as a consequence.

Additionally, GPLVMs were shown to be useful for reconstructing inputs with
partially observed features. Such a scenario often occurs in image
reconstruction or denoising tasks.

Although GPLVMs are very useful for unsupervised tasks, the original method
assumed access to the full kernel matrix, which entails storing the entire
training dataset and inverting an $n\times n$ kernel matrix. Lawrence
[Lawrence, 2007] addressed this issue by showing that most sparse Gaussian
process approaches can be turned into GPLVMs. However, the approach still gave
a MAP estimate of the latent space and risked overfitting to the training
dataset.

Titsias and Lawrence [Titsias and Lawrence, 2010] addressed the overfitting
issue by proposing a Bayesian approach. Instead of finding a MAP solution to
the latent space, they proposed a variational approach. However, naive use of
the variational approach to finding the data distribution introduces
intractabilities. Titsias and Lawrence addressed the problem by incorporating
Titsias’s variational approach to SGPs [Titsias, 2009]. The introduction of
pseudo points [Titsias, 2009] canceled out the intractable terms in the
variational bound for GPLVMs and resulted in a feasible optimization bound
shown below

$\displaystyle\mathcal{F}(q)\geq\log\left(\int
e^{\left\langle\log\mathcal{N}\left(\mathbf{y}_{d}\mid\boldsymbol{\alpha}_{d},\beta^{-1}I\right)\right\rangle_{q(\mathbf{X})}}p\left(\mathbf{m}_{d}\right)d\mathbf{m}_{d}\right)$
$\displaystyle-\frac{\beta}{2}\operatorname{Tr}\left(\left\langle
K_{nn}\right\rangle_{q(\mathbf{X})}\right)+\frac{\beta}{2}\operatorname{Tr}\left(K_{mm}^{-1}\left\langle
K_{mn}K_{nm}\right\rangle_{q(\mathbf{X})}\right)$

Here, $q$ is the variational distribution over pseudo points $\mathbf{m}$ and
$\mathbf{\alpha}=K_{nm}K^{-1}_{mm}\mathbf{m}$. The subscript $d$ is used to
denote each dimension of the features.

GPLVMs have been applied to numerous applications and their variants that we
did not cover here. The readers are referred to [Li and Chen, 2016] for an in-
depth survey of GPLVMs.

##  6 Deep Gaussian Processes

Although SGPs addressed the computation cost issue, GPs remain inapplicable to
a lot of applications. The reason being the kernel function. The most commonly
used kernel functions have relatively simple similarity metrics. However, in
specific datasets, one might have to use a different similarity metric in
different regions of the input space. A similarity metric that can extract
such features would have to utilize a hierarchical structure for feature
extraction.

One strategy to the problem would be to stack GPs similar to how Perceptrons
are stacked in an MLP. But, stacking GPs such that one layer’s output becomes
the following layer’s input makes them highly non-linear and intractable to
analytical solutions. Moreover, a stacked GP would not even correspond to a GP
anymore as the posterior distribution can take on any arbitrary distribution.
However, such methods are usually referred to as Deep Gaussian Processes
(DGPs). Several authors have tried to model and fit such models; this section
explains such methods’ development.

One of the earliest DGP approaches is that of Lawrence and Moore [Lawrence and
Moore, 2007]. They considered a GPLVM model, but a GP was assumed for the
input space prior distribution, making it a two layered DGP. The DGP resulted
in the following likelihood function, which cannot be marginalized
analytically.

$p\left(\mathbf{y}\mid\mathbf{t}\right)=\int
p\left(\mathbf{y}\mid\mathbf{X}\right)p\left(\mathbf{X}\mid\mathbf{t}\right)\mathrm{d}\mathbf{X}\\\
$

Here, $\mathbf{t}$ is the input GP at the input layer, and $\mathbf{X}$ are
the intermediate representations passed to the second layer GP. Lawrence and
Moore considered a MAP solution to the above problem. This was achieved by
maximizing the following.

$p\left(\mathbf{X}\mid\mathbf{y,t}\right)=\log
p\left(\mathbf{y}\mid\mathbf{X}\right)+\log
p\left(\mathbf{X}\mid\mathbf{t}\right)+\text{const}\\\ $

The authors also showed that deeper hierarchies could be modeled with such an
approach. However, the model was limited to a MAP solution which is highly
susceptible to overfitting.

Damianou et al. [Damianou et al., 2011, Damianou, 2015] proposed a variational
approach to the overfitting problem. They also considered a 2-layered stacked
GP, but a naive variational bound for such a model introduces intractabilities
similar to a GPLVM. However, the authors showed that the variational approach
used for Bayesian GPLMVs in [Titsias and Lawrence, 2010] could also be
utilized to formulate a variational bound for a 2-layered GP. The final bound
is shown below, with $q(\mathbf{X})$ as the variational distribution

$\mathcal{F}=\hat{\mathcal{F}}-\operatorname{KL}(q(X)\|p(X\mid\mathbf{t}))\\\
\hat{\mathcal{F}}=\int q(X)\log p(y\mid\mathbf{f})p(\mathbf{f}\mid
X)\mathrm{d}X\mathrm{~{}d}\mathbf{f}$

Furthermore, Damianou and Lawrence [Damianou and Lawrence, 2013] improved the
bound shown above by generalizing the variational bound to a DGP with an
arbitrary number of layers. The bound shown below can be used on a DGP with
two or more layers.

$\displaystyle\mathcal{F}$
$\displaystyle=\mathbf{g}_{Y}+\mathbf{r}_{X}+\mathcal{H}_{q(\mathbf{X})}-\operatorname{KL}(q(\mathbf{Z})\|p(\mathbf{Z}))$
$\displaystyle\mathbf{g}_{Y}$
$\displaystyle=g\left(\mathbf{Y},\mathbf{F}^{Y},\mathbf{U}^{Y},\mathbf{X}\right)$
$\displaystyle=\left\langle\log
p\left(\mathbf{Y}\mid\mathbf{F}^{Y}\right)+\log\frac{p\left(\mathbf{U}^{Y}\right)}{q\left(\mathbf{U}^{Y}\right)}\right\rangle_{p\left(\mathbf{F}^{Y}\mid\mathbf{U}^{Y},\mathbf{X}\right)q\left(\mathbf{U}^{Y}\right)q(\mathbf{X})}$
$\displaystyle\mathbf{r}_{X}$
$\displaystyle=r\left(\mathbf{X},\mathbf{F}^{X},\mathbf{U}^{X},\mathbf{Z}\right)$
$\displaystyle=\left\langle\log
p\left(\mathbf{X}\mid\mathbf{F}^{X}\right)+\log\frac{p\left(\mathbf{U}^{X}\right)}{q\left(\mathbf{U}^{X}\right)}\right\rangle_{p\left(\mathbf{F}^{X}\mid\mathbf{U}^{X},\mathbf{Z}\right)q\left(\mathbf{U}^{X}\right)q(\mathbf{X})q(\mathbf{Z})}$

Here, $\mathbf{Y}$ represents the output multidimensional label space,
$\mathbf{Z}$ represents the latent variables in the input layer, and
$\mathbf{X}$ represents the latent inputs in intermediate layers. $\mathbf{U}$
and $\mathbf{F}$ are the values of the latent function corresponding to
inducing points and latent inputs, respectively; their superscript denotes the
layer to which they belong. Additionally, $\mathcal{H}$ represents the entropy
of the distribution shown in its subscript, and $\operatorname{KL}$ is the
standard KL-divergence. Fig. 4 shows the DGP model architecture from [Damianou
and Lawrence, 2013].

Figure 4: Deep Gaussian Process model overview [Damianou and Lawrence, 2013].

Again, the crux of the approach relied on the variational trick of introducing
inducing points as presented in [Titsias and Lawrence, 2010]. Damianou and
Lawrence conducted experiments on the MNIST dataset wherein they showed a
5-layer DGP could be used for the image classification task.

A limitation of the approach in [Damianou and Lawrence, 2013] is that the
number of variational parameters that need to be learned increases linearly
with the number of data points in the training set. And it involved inverting
matrices which is a computationally expensive operation, thereby limiting its
scalability. Dai et al. [Dai et al., 2015] addressed this by introducing a
back-constraint. The constraint allowed them to define the latent variables’
mean terms as deterministic functions of the latent variables themselves via
an MLP. The approach reduced the number of variational parameters. Moreover,
Dai et al. [Dai et al., 2015] also showed that their approach could be trained
in a distributed manner, allowing the model to be scaled to large datasets.

Salimbeni and Deisenroth [Salimbeni and Deisenroth, 2017] recently proposed an
approach that addressed the layer independence issue of prior DGP approaches.
The DGP in [Damianou and Lawrence, 2013] assumed independence of GPs across
layers and only considered the correlations within a layer. However, Salimbeni
and Deisenroth argue that such an approach is equivalent to a single GP with
each GP’s input coming from a GP itself. The authors also stated that they
found that some layers get turned off when using the DGP in [Damianou and
Lawrence, 2013].

Salimbeni and Deisenroth [Salimbeni and Deisenroth, 2017] presented a new
variational bound that retains the exact model posterior similar to [Damianou
and Lawrence, 2013] while maintaining the correlations both within and across
adjacent layers. However, Salimbeni and Deisenroth showed such an approach is
infeasible to analytic computations, but the bound could still be optimized
with MCMC sampling techniques. Such an approach is computationally expensive.
But, it can be parallelized by exploiting the factorization of the DGP across
output dimensions. Furthermore, the method also required sampling methods
during inference, but its performance is substantially better than prior
works.

$\mathcal{F}=\sum_{i=1}^{N}\mathbb{E}_{q\left(\mathbf{f}_{i}^{L}\right)}\left[\log
p\left(\mathbf{y}_{n}\mid\mathbf{f}_{n}^{L}\right)\right]-\sum_{l=1}^{L}\mathrm{KL}\left[q\left(\mathbf{U}^{l}\right)\|p\left(\mathbf{U}^{l};\mathbf{Z}^{l-1}\right)\right]$

In the optimization bound shown above from [Salimbeni and Deisenroth, 2017],
the subscripts are used to denote each data sample in the dataset, and
superscripts are used to indicate a layer in the DGP. The rest of the terms
follow the same convention as the one used by Damianou et al. [Damianou and
Lawrence, 2013].

I briefly mentioned that DGPs do not necessarily correspond to Gaussian
processes. Still, the methods discussed so far do model the posterior
distribution as a Gaussian, each with its assumptions. Havasi et al. [Havasi
et al., 2018] presented a technique that departs even more from a conventional
GP. The authors show that since a Gaussian is uni-modal, using it to model the
posterior will result in poor results. Instead, they suggest using a multi-
modal distribution that can better capture the true posterior distribution.

However, it is not possible to formulate an analytical solution to a multi-
modal posterior. We can use variational inference to learn a multi-modal
posterior. Still, one needs to determine the exact form of the variational
distribution, which is difficult as we usually don’t know the posterior
distribution in advance. Havasi et al. [Havasi et al., 2018] circumvent this
issue by using Stochastic Gradient Hamiltonian Monte Carlo (SGMCMC) [Chen et
al., 2014] method to estimate the posterior. The approach can determine the
inducing points by sampling from the true posterior instead of using a
variational distribution.

Although the approach far exceeds the performance of prior DGPs and is the
current state-of-the-art, it still has its limitations. Remarkably, the SGMCMC
method is difficult to tune as it introduces its own parameters in addition to
the ones that already have to be estimated for DGP. Several MCMC method
variants attempt to improve upon SGMCMC, but none of those approaches have
been applied to DGPs.

The DGPs we discussed so far attempted to develop variants of GPs that can
model hierarchical features in data, which was done by assuming a feed forward
network with each node of the network being modeled as a GP. It is the most
popular method to address the problem and has resulted in approaches that can
get reasonably promising results. However, there are other approaches as well
which do not consider such an explicit feed forward network.

Wilson et al. [Wilson et al., 2016] proposed a method that uses deep neural
networks as the kernel function, termed deep kernel. Unlike a Gaussian kernel,
the deep kernel produced a vector output, and a GP was assigned to each of the
vector elements. Wilson et al. further combined the GPs with an additive
structure that facilitated its training with an analytical bound.

Wilson et al. [Wilson et al., 2016] showed that their approach was good at
several tasks. However, the deep neural network architecture needs to be task
specific, and its parameters are susceptible to overfitting given its large
parameter count.

Another interesting perspective on DGPs was presented by Lee et al. [Lee et
al., 2018]. Until now, all the discussed GPs with linear latent functions are
combined in different ways to achieve an aggregate non-linear latent function
space. Lee et al. developed a method that considers an entire function space
that consists of non-linear functions. Unlike prior approaches, the function
space was not constrained to a specific subspace as a consequence of a
particular kernel function being used. The approach can be regarded as a
generalization of Neil [Neal, 1996] who showed the equivalence of an
infinitely wide single layer neural network to a GP. Lee et al. showed the
equivalence of a GP to an infinitely wide deep neural.

Lee et al. [Lee et al., 2018] showed the approach was on par with some neural
networks trained with gradient descent while retaining its uncertainty
estimation. Furthermore, the uncertainty estimates were proportional to the
model accuracy. But, the approach has a polynomial increasing kernel matrix
making it infeasible for some problems. Moreover, the approach only considered
deep neural networks with fully connected layers and Relu activation
functions.

Garnelo et al. [Garnelo et al., 2018] presented an approach with similar ethos
and introduced Neural Processes (NPs). However, instead of considering a
neural network with asymptotically increasing depth and width, a deep neural
network was used in place of a Gaussian distribution parametrized by a kernel
function to define $p(\mathbf{f}\mid\mathbf{X})$.

$\displaystyle p(\mathbf{y}\mid\mathbf{f})$
$\displaystyle=\mathcal{N}\left(\mathbf{f},\beta^{-1}\mathbf{I}\right)$
$\displaystyle p(\mathbf{y}\mid\mathbf{X})$ $\displaystyle=\int
p(\mathbf{y}\mid\mathbf{f})\underbrace{p(\mathbf{f}\mid\mathbf{X})}_{\text{Deep
neural network}}d\mathbf{f}$

The deep neural network was trained using amortized variational inference. The
consequence of such an approach is that the function space defined by deep
neural networks allows us to extract hierarchical features and retain a
probabilistic interpretation. The model, however, needs to be trained with
Meta-learning, an approach wherein multiple diverse datasets or tasks are used
to train the same model. Meta-learning is used because each function in the
function space corresponds to a sequence of inputs or a task. Considering
multiple tasks allows the DNN to approximate the variability of the function
space. While training, a context vector $\mathbf{r}_{c}$ is passed to the DNN
to indicate the task currently being considered as shown in Fig: 5.

Figure 5: Neural Processes model architecture [Kim et al., 2019].

Moreover, to retain the probabilistic interpretation, a latent variable $z$
was introduced, which captured uncertainty in the context data. The
implication being that, unlike vanilla GPs whose uncertainty comes from the
kernel function and its function space, NPs do this with data. As such, the
context being provided could significantly influence the model’s performance
and can be considered to be similar to inducing points in SGPs.

Additionally, the model does not assume a Gaussian prior or posterior,
allowing one to fit any data distribution. Garnelo et al. showed that their
approach produces good predictive distributions while being parameter
efficient and fast compared to a vanilla Gaussian process.

Nonetheless, the method assumed a predefined DNN model architecture which
needs to be task specific. Also, the model is only an approximation to some
stochastic process using a DNN. But, it isn’t possible to guarantee the
approximation quality of the DNN. Furthermore, the meta-learning requirement
imposes a large training computation requirement, and the datasets considered
have to be similar to the primary dataset of interest.

Finally, Yang et al.[Yang et al., 2020] recently proposed an Energy-based
Process (EBP). EBPs are a generalization of Neural Processes as they utilize
Energy-based models [LeCun et al., 2006] to approximate
$p(\mathbf{f}\mid\mathbf{X})$ instead of a MAP trained DNN as shown below,
where $f_{w}$ is the Energy model and $Z$ is the partition function.

$p(\mathbf{f}\mid\mathbf{X})=\frac{\exp(f_{w}(\mathbf{X}))}{Z(f_{w})}$

However, by utilizing an Energy based model, the authors were able to show
that vanilla GPs and NPs can be recovered from EBPs as special cases. The
Energy based formulation also allows one to approximate the conditional
$p(\mathbf{f}\mid\mathbf{X})$ with arbitrary distributions, unlike GP and NPs,
which are limited to Gaussian distributions and the distribution defined by
the DNN, respectively.

Unlike feedforward networks which are trained to predict a label $y$ given an
input $\mathbf{X}$, Energy based models predict the Energy between a pair of
$(\mathbf{X},y)$. A well trained Energy-based will output low Energy for well
matched $(\mathbf{X},y)$ pairs and high Energy for mismatched pairs. As such,
the prediction task in these models becomes a minimization task, wherein one
needs to find the label $y$, which has low Energy for the given data
$\mathbf{X}$.

The consequence of training such a model to approximate our conditional in the
stochastic process is that the function space is not constrained to any
predefined subspace. However, Energy based models are challenging to train and
require several tricks to stabilize the training process. Furthermore, it
takes longer to train such models similar to models trained with Meta-
learning.

##  7 Discussion

Gaussian processes have come a long way from their origins. Although many of
the limitations have been addressed, there remain open problems and research
directions that have not been thoroughly explored.

One such problem is the assumption of factored output dimensions. The
assumption has been made in all the methods mentioned in this paper. It
dictates that each output dimension is independent of the other. The
assumption allowed factorizations that simplify some derivations, and in some
cases, the assumption is required for the method to be tractable. However, the
assumption might not hold in some datasets. Addressing this factorization
assumption would be an interesting line of research.

Another issue is that most SGP and DGP methods require careful model
initializations and hyperparameter tuning without which, the model does not
converge. Nevertheless, there aren’t any formalized rules for determining
model initializations and hyperparameters for most methods that could
guarantee good model convergence. The tuning issue is particularly a problem
when using MCMC approaches and remains to be solved.

Additionally, MCMC methods have been shown to be successful at training DGP.
But, the method need not be limited to DGP. Indeed, the technique might result
in good results even for SGP. The primary motivation of using MCMC methods was
to address the non-Gaussian posterior. Although vanilla GPs might not have a
non-Gaussian posterior, the assumptions made in inducing methods often change
this. As such, it might be worth exploring the feasibility of MCMC for
training SGPs.

Similarly, there are several variants of SG-MCMC methods that have not been
benchmarked against DGPs. Havasi et al. [Havasi et al., 2018] only considered
the original SGMCMC [Chen et al., 2014] approach. However, numerous variants
of the method have been introduced that improve upon SGMCMC, some of which
might result in stable training dynamics.

Deep kernels in themselves were found to be highly susceptible to overfitting.
However, Wilson et al. [Wilson et al., 2016] only considered vanilla DNNs.
But, DNNs could be used as Bayesian approximators as shown by Gal and
Ghahramani [Gal and Ghahramani, 2016]. Such an approach might mitigate some of
the overfitting issues. ’Furthermore, one might also consider methods such as
Bayes by Backprob [Blundell et al., 2015] to train deep kernels. It would be
interesting to find out the ramifications of such an approach to deep kernels.

Garnelo et al. [Garnelo et al., 2018] considered a similar approach but, they
consider a DNN approximation to a distribution over the function space itself.
It required an explicit definition of a DNN, which needs to be task specific.
Also, the model performance is dependent on the context vector to estimate
prediction uncertainty which isn’t consistent with a proper stochastic
process. Kim et al. [Kim et al., 2019] introduced a variant of Neural
Processes which uses attention to improve the context vectors. Perhaps we
could modify the attention mechanism to take the test data into account and
generate an uncertainty estimate that incorporates test data.

Moreover, most methods assume relatively constrained function space, either
with the formulation of a kernel function or a DNN. However, that need not be
the case; perhaps we can consider multiple function spaces by utilizing models
such as chunked hypernetworks [von Oswald et al., 2020] to generate both the
model parameters and the model architecture. Thus greatly expanding the
stochastic process’s modeling capacity.

Energy-based models appear to be another viable approach to expanding the
function space but, the method is challenging to train and incurs substantial
computational costs. Moreover, even model inference is an expensive operation
requiring Hamiltonian Markov Chain methods for sampling.

Finally, there is the issue of scalability. Although some DGP methods have
been shown to scale well to large datasets, they have not been thoroughly
benchmarked on highly structured datasets such as Imagenet [Deng et al.,
2009]. The problem lies in the model depth required to achieve good
performance on such a dataset. Unlike MNIST, Imagenet requires DNNs that are
substantially deeper. DGPs, however, are usually only tested on models with up
to ten layers. It would be indispensable to study and understand how DGPs
scale to such a dataset.

##  8 Conclusion

Gaussian Processes in it themselves are fascinating. Their non-parametric
form, analytical properties, and ability to model uncertainty are coveted in
machine learning. However, they are plagued with limitations, particularly
their significant computation and storage cost. Also, conventional kernel
functions restrict the family of functions that a GP could model.

Sparse Gaussian Processes attempt to address the storage and computation cost.
One dominant approach to SGPs is to use the Nyström approximation. The
approach entails using variational methods to model the distribution of pseudo
points for a fully Bayesian treatment. Several methods have been proposed
along this line of research, each with its advantages and limitations.

Furthermore, GPLVM was a step towards DGP. However, hierarchical feature
representation was not the intended use case. It was proposed as an approach
for probabilistic PCA and unsupervised learning. The Bayesian GPLVM on
improved upon the original method by introducing a purely Bayesian training
approach. BGPLVMs facilitated the propagation of latent space uncertainty to
the posterior, thereby establishing a technique to propagate uncertainties
through non-linearities in GPs.

Most of the DGP approaches considered SGPs and GPLVMs to address the issue of
hierarchical feature representation. The primary trend in DGPs is to stack GPs
in a feedforward manner and train them with approaches used to train SGPs and
GPLVMs. However, such an approach has its limitations. The optimization bounds
that were developed were not always tight, and some methods were limited to
analytical solutions, which imposed scalability limitations on such
techniques.

Moreover, stacking GPs makes the model parametric as it requires a predefined
model depth and layer width. Lee et al. [Lee et al., 2018] considered these
issues and attempted to solve them by modeling the latent function space as
the space of deep neural networks. But, the approach is not feasible for real
world applications as of yet and requires more work to be done to get there.

Garnelo et al. [Garnelo et al., 2018] consider a stochastic process
parametrized with a DNN instead of a Gaussian distribution with a kernel
function to define the latent function space. Still, the approach requires
modeling a task specific neural network and is only an approximation to an
unknown stochastic process. Energy-based Processes address this limitation,
but the method isn’t mature enough as of yet.

In conclusion, GPs are an excellent approach to model datasets. The overall
trend of the field seems to be transitioning away from the Gaussian assumption
and consider general Stochastic processes. The method has come a long way from
its infancy, but there are still open problems that need to be addressed for
it to be elevated to the prominence that it deserves.

## References

  * [Bishop, 2016] Bishop, C. (2016).  Pattern Recognition and Machine Learning.  Information Science and Statistics. Springer New York.
  * [Blundell et al., 2015] Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. (2015).  Weight Uncertainty in Neural Network.  In Bach, F. and Blei, D., editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1613–1622, Lille, France. PMLR.
  * [Chen et al., 2014] Chen, T., Fox, E., and Guestrin, C. (2014).  Stochastic Gradient Hamiltonian Monte Carlo.  In Xing, E. P. and Jebara, T., editors, Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 1683–1691, Bejing, China. PMLR.
  * [Dai et al., 2015] Dai, Z., Damianou, A., González, J., and Lawrence, N. (2015).  Variational Auto-encoded Deep Gaussian Processes.  arXiv preprint arXiv:1511.06455.
  * [Damianou, 2015] Damianou, A. (2015).  Deep Gaussian Processes and Variational Propagation of Pncertainty.  PhD thesis, University of Sheffield.
  * [Damianou and Lawrence, 2013] Damianou, A. and Lawrence, N. D. (2013).  Deep Gaussian Processes.  In Carvalho, C. M. and Ravikumar, P., editors, Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics, volume 31 of Proceedings of Machine Learning Research, pages 207–215, Scottsdale, Arizona, USA. PMLR.
  * [Damianou et al., 2011] Damianou, A. C., Titsias, M. K., and Lawrence, N. D. (2011).  Variational Gaussian Process Dynamical Systems.  In Proceedings of the 24th International Conference on Neural Information Processing Systems, NIPS’11, page 2510–2518, Red Hook, NY, USA.
  * [Deng et al., 2009] Deng, J., Dong, W., Socher, R., Li, L., Kai Li, and Li Fei-Fei (2009).  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255.
  * [Gal and Ghahramani, 2016] Gal, Y. and Ghahramani, Z. (2016).  Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.  In Balcan, M. F. and Weinberger, K. Q., editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059, New York, New York, USA. PMLR.
  * [Garnelo et al., 2018] Garnelo, M., Schwarz, J., Rosenbaum, D., Viola, F., Rezende, D. J., Eslami, S., and Teh, Y. W. (2018).  Neural Processes.  arXiv preprint arXiv:1807.01622.
  * [Havasi et al., 2018] Havasi, M., Hernández-Lobato, J. M., and Murillo-Fuentes, J. J. (2018).  Inference in Deep Gaussian Processes Using Stochastic Gradient Hamiltonian Monte Carlo.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, page 7517–7527, Red Hook, NY, USA. Curran Associates Inc.
  * [Hensman et al., 2013] Hensman, J., Fusi, N., and Lawrence, N. D. (2013).  Gaussian Processes for Big Data.  In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, UAI’13, page 282–290, Arlington, Virginia, USA. AUAI Press.
  * [Hornik et al., 1989] Hornik, K., Stinchcombe, M., and White, H. (1989).  Multilayer Feedforward Networks Are Universal Approximators.  Neural Netw., 2(5):359–366.
  * [Kim et al., 2019] Kim, H., Mnih, A., Schwarz, J., Garnelo, M., Eslami, A., Rosenbaum, D., Vinyals, O., and Teh, Y. W. (2019).  Attentive Neural Processes.  In International Conference on Learning Representations.
  * [Klebaner, 2012] Klebaner, F. (2012).  Introduction to Stochastic Calculus with Applications.  Introduction to Stochastic Calculus with Applications. Imperial College Press.
  * [Lawrence, 2004] Lawrence, N. (2004).  Gaussian Process Latent Variable Models for Visualisation of High Dimensional Data.  In Thrun, S., Saul, L., and Schölkopf, B., editors, Advances in Neural Information Processing Systems, volume 16. MIT Press.
  * [Lawrence, 2005] Lawrence, N. (2005).  Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models.  Journal of Machine Learning Research, 6(60):1783–1816.
  * [Lawrence, 2007] Lawrence, N. (2007).  Learning for Larger Datasets with the Gaussian Process Latent Variable Model.  In Meila, M. and Shen, X., editors, Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, volume 2 of Proceedings of Machine Learning Research, pages 243–250, San Juan, Puerto Rico. PMLR.
  * [Lawrence and Moore, 2007] Lawrence, N. D. and Moore, A. J. (2007).  Hierarchical Gaussian Process Latent Variable Models.  In Proceedings of the 24th International Conference on Machine Learning, ICML ’07, page 481–488, New York, NY, USA. Association for Computing Machinery.
  * [LeCun et al., 2006] LeCun, Y., Chopra, S., Hadsell, R., Ranzato, M., and Huang, F. (2006).  A tutorial on energy-based learning.  Predicting structured data, 1(0).
  * [Lee et al., 2018] Lee, J., Sohl-dickstein, J., Pennington, J., Novak, R., Schoenholz, S., and Bahri, Y. (2018).  Deep Neural Networks as Gaussian Processes.  In International Conference on Learning Representations.
  * [Li and Chen, 2016] Li, P. and Chen, S. (2016).  A Review on Gaussian Process Latent Variable Models.  CAAI Transactions on Intelligence Technology, 1(4):366–376.
  * [Liu et al., 2020] Liu, H., Ong, Y.-S., Shen, X., and Cai, J. (2020).  When Gaussian Process Meets Big Data: A Review of Scalable GPs.  IEEE transactions on neural networks and learning systems, 31(11):4405—4423.
  * [McGrayne, 2011] McGrayne, S. (2011).  The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, & Emerged Triumphant from Two Centuries of C.  Matematicas (E-libro). Yale University Press (Ignition).
  * [Murphy, 2012] Murphy, K. (2012).  Machine Learning: A Probabilistic Perspective.  Adaptive Computation and Machine Learning series. MIT Press.
  * [Neal, 1996] Neal, R. (1996).  Bayesian Learning for Neural Networks.  Lecture Notes in Statistics. Springer New York.
  * [Pearl, 2000] Pearl, J. (2000).  Causality: Models, Reasoning, and Inference.  Cambridge University Press.
  * [Pearl and Mackenzie, 2018] Pearl, J. and Mackenzie, D. (2018).  The Book of Why: The New Science of Cause and Effect.  Basic Books.
  * [Rasmussen and Williams, 2006] Rasmussen, C. and Williams (2006).  Gaussian Processes for Machine Learning.  Adaptive computation and machine learning. MIT Press.
  * [Rosenthal, 2006] Rosenthal, J. (2006).  A First Look at Rigorous Probability Theory.  World Scientific.
  * [Salimbeni and Deisenroth, 2017] Salimbeni, H. and Deisenroth, M. (2017).  Doubly Stochastic Variational Inference for Deep Gaussian Processes.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * [Snelson and Ghahramani, 2006] Snelson, E. and Ghahramani, Z. (2006).  Sparse Gaussian Processes using Pseudo-inputs.  In Weiss, Y., Schölkopf, B., and Platt, J., editors, Advances in Neural Information Processing Systems, volume 18. MIT Press.
  * [Titsias, 2009] Titsias, M. (2009).  Variational Learning of Inducing Variables in Sparse Gaussian Processes.  In van Dyk, D. and Welling, M., editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA. PMLR.
  * [Titsias and Lawrence, 2010] Titsias, M. and Lawrence, N. D. (2010).  Bayesian Gaussian Process Latent Variable Model.  In Teh, Y. W. and Titterington, M., editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 844–851, Chia Laguna Resort, Sardinia, Italy. JMLR Workshop and Conference Proceedings.
  * [von Oswald et al., 2020] von Oswald, J., Henning, C., Sacramento, J., and Grewe, B. F. (2020).  Continual Learning with Hypernetworks.  In International Conference on Learning Representations.
  * [Williams and Seeger, 2001] Williams, C. and Seeger, M. (2001).  Using the Nyström Method to Speed Up Kernel Machines.  In Leen, T., Dietterich, T., and Tresp, V., editors, Advances in Neural Information Processing Systems, volume 13. MIT Press.
  * [Wilson et al., 2016] Wilson, A., Hu, Z., Salakhutdinov, R., and Xing, E. (2016).  Stochastic Variational Deep Kernel Learning.  Advances in Neural Information Processing Systems, pages 2594–2602.
  * [Yang et al., 2020] Yang, M., Dai, B., Dai, H., and Schuurmans, D. (2020).  Energy-based processes for exchangeable data.  In III, H. D. and Singh, A., editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 10681–10692. PMLR.