# Literature Review: Graph Kernels in Chemoinformatics

James Young

(Summer 2022)

## Preface

The purpose of this review is to introduce the reader to _graph kernels_ and
the corresponding literature, with an emphasis on those with direct
application to _chemoinformatics_. Graph kernels are functions that allow for
the inference of properties of molecules and compounds, which can help with
tasks such as finding suitable compounds in drug design. The use of kernel
methods is but one particular way two quantify similarity between graphs. We
restrict our discussion to this one method, although popular alternatives have
emerged in recent years, most notably graph neural networks.

The first two chapters provide the necessary background on statistical
learning theory and common techniques. Chapter 1 is a crash course on basic
_statistical learning_ concepts that provides motivations for the methods seen
throughout the review. In chapter 2, a classification model based on geometric
separation of data, called _support vector machines_ , will be discussed.
Support vector machines make use of a class of function known as _kernels_ ,
which act as a similarity measure between data points. Different kernels
emphasize different properties of the underlying data, thus the task of
choosing the right kernel for a given problem is of great interest. Chapter 3
is the main body of the review. Here, we give a broad overview of the
literature of kernels defined on graphs, with an emphasis on kernels that may
be useful in chemoinformatics. Our goal is to provide a list of popular graph
kernels along with their properties, computational complexity, and test
results so that the reader can become acquainted with many of the important
terms, techniques, and papers in the literature so that they can determine if
graph kernels can be helpful in their research.

###### Contents

  1. Preface
  2. 1 Supervised Statistical Learning
    1. 1.1 Regression Problems
      1. 1.1.1 The Setup
      2. 1.1.2 General Modelling Approaches and Measuring Error
      3. 1.1.3 Resampling
    2. 1.2 Classification Problems
  3. 2 Support Vector Machines
    1. 2.1 Hyperplane Classification
      1. 2.1.1 Separable Observations and the Maximal Margin Classifier
      2. 2.1.2 Non-Separable Observations and the Support Vector Classifier
    2. 2.2 Nonlinear Classification
      1. 2.2.1 Reproducing Kernel Hilbert Spaces
      2. 2.2.2 Using Kernels for Nonlinear Classification
      3. 2.2.3 SVMs for Multiclass Classification
  4. 3 Graph Kernels
    1. 3.1 Introductory Concepts
      1. 3.1.1 Essential Graph Theory
      2. 3.1.2 Convolution kernels
    2. 3.2 Baseline Graph Kernels
      1. 3.2.1 Simple Examples
      2. 3.2.2 Direct Product Kernel
      3. 3.2.3 Shortest-Path Graph Kernel
    3. 3.3 Specialized Examples
      1. 3.3.1 Marginalized Kernels
      2. 3.3.2 Tree-Pattern Graph Kernel
      3. 3.3.3 Tanimoto Kernels
    4. 3.4 Boosting Kernel Performance
      1. 3.4.1 Graph RBF Kernel
      2. 3.4.2 Weisfeiler-Lehman Kernels
      3. 3.4.3 Optimal Assignment Kernels
      4. 3.4.4 Notes on Explicit and Implicit Computation
    5. 3.5 Software and Data sets
      1. 3.5.1 Common Data Sets
      2. 3.5.2 Computing Graph Kernels
    6. 3.6 Further Reading

##  Chapter 1 Supervised Statistical Learning

In this first chapter, we briefly touch on what statistical learning is and
what we are trying to accomplish with it. While mostly self-contained, it is
not intended to be a comprehensive source on the subject’s foundations. As
such, we are only going to cover the most essential theoretical basics, as
well as outline the general notation and terminology used throughout the
review. For a more detailed and example driven treatment of this material, see
the book (James et al., 2021) in which this chapter is based upon, or (Hastie
et al., 2009) for a more advanced and mathematically rigorous approach.

The general setting is the following: We receive a sample containing $n$
observations $y_{i}$ ($i\in\\{1,\,\cdots,n\\}$), each with a corresponding
vector containing $p$ predictors (or _inputs_)
$\textbf{x}_{i}=(x_{1},\,\cdots,x_{p})^{T}$. Often, it will be useful to think
of the observations and predictors as a set of ordered pairs
$(\textbf{x}_{i},\,y_{i})_{i=1}^{n}$ inhabiting a higher dimensional space,
with the $y_{i}$’s being the response to the inputs $\textbf{x}_{i}$. This
notation will appear throughout the review. The $p$-dimensional space of
predictors will be referred to as _input space_. Now, given this sample, what
can we say about other members in the population? In reality, it is not
feasible to measure every single input of every member of an entire
population. The basic idea of supervised statistical learning is to build a
statistical model from our sample that allows us to predict observations based
on inputs not found in our sample. We typically won’t fit a model to _all_ of
our sample. The subset in which we do is called the _training data_ , which is
used to _learn_ how the predictors relate to the observations. The remaining
members in the sample are apart of the _testing data_. Inputs from this subset
can be plugged into the model, allowing us to observe the predictive power of
the model on new data. The term _supervised_ is in reference to the
observations being present in our samples, which means we can measure the
predictive power of our model by quantifying how close the predicted
observations are to the true observations. There is also what is known as
unsupervised learning, where no response is observed from any predictor in the
sample. In this situation, we have nothing to base predictions off of, so
instead the analysis is based in understanding how particular observations
relate to one another. Our focus for the duration of the review will be on
supervised learning.

There are two types of observations we can encounter. Either the observations
have quantitative values, meaning they are assigned a numerical quantity in
which typical calculations make sense (eg. the real numbers $\mathbb{R}$), or
they have qualitative values, meaning the value (or often label) assigned to
them is meant to denote a classification, rather than a continuous quantity.
Under this assumption, there is no meaningful way to combine the class labels
as we do for quantitative variables.

###  1.1 Regression Problems

Here, we begin our journey with the foundational elements of supervised
statistical learning. The theory has been framed in terms of regression
problems, although much of it applies to classification problems as well. The
differences in the material will be highlighted in section 1.2 when we look at
equivalent tools for classification problems.

####  1.1.1 The Setup

The question now is how do we build such a model described in the opening
paragraph of this chapter? As per the naming convention, the predictors are
thought to have a level of influence over their respective responses. We can
model this relationship by writing

$y=f(\textbf{x})+\varepsilon$ (1.1)

where $f$ is an unknown function representing the relationship between x and
$y$ (or the information about $y$ provided by x), and $\varepsilon$ represents
a random error term that is independent of x. This error term essentially
accounts for unmeasured variables that, in theory, would be useful in
predicting $y$. We assume that $\varepsilon$ has mean zero, as for larger and
larger samples we expect the responses predicted by $f$ to over and
underestimate $y$ equally. Of course in practice, $f$ is not known, so we must
resort to estimating it. We represent this estimation by

$\hat{y}=\hat{f}(\textbf{x})$ (1.2)

where $\hat{f}$ is the estimate of $f$ and $\hat{y}$ is its prediction for
$y$. We can gain a theoretical understanding of the accuracy of our estimate
by understanding the types of errors it introduces. The _reducible error_ is
the error generated by $\hat{f}$ not being a perfect estimate of $f$. We call
it reducible because if we came up with a better estimate, the error would
decrease. The _irreducible error_ is the error generated by the variance of
our random error term $\varepsilon$. It is irreducible because unlike
$\hat{f}$, we have no control over it, and thus cannot make it smaller. Even
if we could estimate $f$ perfectly so that $\hat{f}\equiv f$ and consequently
$\hat{y}=f(\textbf{x})$, our prediction $\hat{y}$ would still have error built
into it because $y$ also depends on the random variable $\varepsilon$, which
is independent of x. Thus, it cannot be predicted by $f$. Symbolically, let us
recall the expectation operator $E$. The expectation of a random variable $X$,
written $E[X]$, is (informally) a measure of the average value of independent
outcomes of $X$. Depending on context, this can involve a sum over all
outcomes of $X$, each multiplied by their corresponding probability of
occurring, or as an integral involving a probability density function of $X$,
for example. These formulations imply that it is a linear operator that fixes
constants. Returning to our discussion on the theoretical accuracy of
$\hat{y}$ as a predictor of $y$, with a little algebra it can be shown that if
$\hat{f}$ and x are fixed, then

$E[(y-\hat{y})^{2}]=[f(\textbf{x})-\hat{f}(\textbf{x})]^{2}+\text{Var}(\varepsilon)$
(1.3)

where $\text{Var}(X)$ is the variance of a random variable $X$ defined as

$\text{Var}(X)\stackrel{{\scriptstyle\textup{def}}}{{=}}E[(x-\mu)^{2}]$ (1.4)

where $\mu=E[X]$ is the mean of $X$. Crucially for equation 1.3 to hold,
$\varepsilon$ has mean zero. The first term on the left-hand side of equation
1.3 is the reducible error and the variance term is the irreducible error.

####  1.1.2 General Modelling Approaches and Measuring Error

The pressing question is: how do we produce a satisfactory estimate $\hat{f}$
for $f$? Furthermore, if we don’t know what $f$ _is_ , how do we quantify the
error introduced by our estimation? The unsatisfying answer for both is: it
depends. While theory certainly is a factor when trying to estimate $f$ for a
particular problem, often one must explore many different approaches to find
one that works well for the problem and the sample available. There are two
general classes of methods used to estimate $f$. The first are known as
_parametric methods_. Here, a common form for $f$ is assumed, such as a
polynomial or a linear combination of elementary functions. For example,
suppose that we think the relationship between predictors and responses is
linear. We would then guess that $f$ has the form

$\displaystyle f(\textbf{x})$ $\displaystyle=\textbf{w}^{T}\textbf{x}+b$
$\displaystyle=w_{p}x_{p}+\cdots+w_{1}x_{1}+b.$ (1.5)

Fitting the model with a method like linear least squares would yield the
coefficients that result in the best fit for the training data, and thus
produce an estimate for the true $f$. An advantage of this approach is that
the problem of finding a function has been reduced to finding $p+1$
parameters. The downside is that the form is often too simple, and doesn’t
match the true $f$. One could trade in some of the simplicity by increasing
the number of parameters, but the cost is that this can lead to _overfitting_
the training data. This phenomenon happens when the estimate is attempting to
match the training data so closely that it picks up on patterns that only
exist due to randomness in the training data (sometime called _noise_). This
issue often becomes apparent once one measures the prediction error of the
testing data, as an overfit model tends to not produce good predictions.

The second option is to use a non-parametric method, which makes no assumption
about the form of $f$ and simply attempts to balance matching the training
data while remaining “smooth”, as in some differentiability condition imposed
at training data points. With no restriction put on the shape of $f$, these
methods can potentially match the true function $f$ more accurately than a
parametric method. However, without imposing more restrictions, we require far
more sample data with this method to be confident in an accurate prediction.

We won’t go into any specific examples of the methods outlined above in this
chapter, rather we will continue to examine the principles of statistical
learning in generality. These methods are a prominent feature in many
statistics textbooks. Chapter 2 is all about a hybrid method in the context of
classification. Instead, we will begin to answer our second question of this
section regarding the performance of our model on the test data. The function
we will use for regression problems is the _mean squared error_ (MSE), which
acts on a sequence of data points as follows:

$(\textbf{x}_{i},\,y_{i})_{i=1}^{m}\mapsto\frac{1}{m}\sum_{i=1}^{m}[y_{i}-\hat{f}(\textbf{x}_{i})]^{2}.$
(1.6)

Why do we chose this particular measurement, say over any other measurement
such as a simple average of the absolute difference? There are a few reasons
that may not be immediately apparent. For one, squaring the difference has the
effect of minimizing the contribution of small errors, while amplifying the
effect of large errors. The benefit of this is that it provides a simple way
to increase the penalty associated with greater prediction errors, while being
more lenient with small errors, since it is expected that our model will not
be perfect in reality. Furthermore, squaring remains a relatively easy
computation while being differentiable everywhere with respect to
$\hat{y}=\hat{f}(\textbf{x})$. To be more precise, it is highly desirable to
be able to differentiate an error function with respect to some parameter of
the model $\hat{y}=\hat{f}(\textbf{x})$. For example, the linear regression
model takes the form $\hat{y}=\textbf{w}^{T}\textbf{x}+b$. Thus
differentiating with respect to w can done to find the global minimum (as it
is a convex optimization problem) of the MSE with respect to w, which is how w
is typically estimated in practice. If the data points in the input of the
mapping 1.6 consist of training data, then we call this measurement the
_training MSE_. When they consist of the test data points, it is then called
the _test MSE_. Ultimately, we do not care as much about the training MSE. We
want to make predictions, so our goal is to minimize the test MSE. Moreover,
it is not generally true that minimizing training MSE will result in a minimum
test MSE. A fundamental truth in statistical learning is that as model
flexibility increases (i.e. how much the model is influenced by the given
data), we tend to observe a monotonic decrease in the training MSE, but the
test MSE traces out a characteristic _U_ -shape. Starting small, as
flexibility begins to increase, the test MSE decreases towards a minimum
before increasing again. With our new terminology, we can restate the
phenomenon of overfitting data as obtaining a small training MSE, but a large
test MSE, which tends to happen for a model that is very flexible. Why does
the _U_ -shape appear in the test MSE? It comes from competing properties
known as _bias_ and _variance_. Let $\hat{\theta}$ be an estimator of a fixed
parameter $\theta$. If we define the MSE of $\hat{\theta}$ as

$\text{MSE}(\hat{\theta})\stackrel{{\scriptstyle\textup{def}}}{{=}}E[(\hat{\theta}-\theta)^{2}]$
(1.7)

and define the bias of the estimator $\hat{\theta}$ as

$\text{Bias}(\hat{\theta})\stackrel{{\scriptstyle\textup{def}}}{{=}}E[\hat{\theta}]-\theta$
(1.8)

then it can be shown that

$\text{MSE}(\hat{\theta})=\text{Bias}(\hat{\theta})^{2}+\text{Var}(\hat{\theta})$
(1.9)

which shows that a decrease in bias is offset by an increase is variance, and
vice-versa. In our situation where we are measuring the expected test MSE of a
learning method for a fixed observation pair $(\textbf{x}_{0},\,y_{0})$, this
Bias-Variance trade-off is found in the equation

Expected Test MSE
$\displaystyle\stackrel{{\scriptstyle\textup{def}}}{{=}}E[(y_{0}-\hat{f}(\textbf{x}_{0}))^{2}]$
$\displaystyle=\text{Bias}(\hat{f}(\textbf{x}_{0}))^{2}+\text{Var}(\hat{f}(\textbf{x}_{0}))+\text{Var}(\varepsilon)$
(1.10)

which can be derived in a similar way as equation 1.3. The squared bias and
the variance of the estimate in this equation are always non-negative, so the
lower bound on the expected test MSE is the variance of the random noise term
$\text{Var}(\varepsilon)$. Intuitively, the expected test MSE at
$\textbf{x}_{0}$ is the value that the average test MSE at $\textbf{x}_{0}$
would approach if we repeatedly estimated $f$ on more and more training sets.
The bias of a statistical learning method represents the error introduced by
approximating a (possibly complex) real-life problem by a (possibly
simplistic) model. An example where high bias can occur is if we use linear
regression to model a relationship between variables that is truly not linear.
Think of this in terms of the colloquial use of the term bias. The model is
essentially ignoring certain information and only looking to validate the
assumption that a linear relationship holds. Variance on the other hand is a
measure of how much our $\hat{f}$ would change if we followed the same
procedure for finding it, but instead used a different training set for
learning.

We can now explain the general behaviour of the test MSE. The relative rate of
change of bias and variance dictates whether test MSE will rise or fall.
Starting with low flexibility, typically what we will see as the flexibility
of a model increases is an initial rapid decrease in bias in comparison to the
increase in variance, corresponding to a decrease in expected test MSE. At
some point however, the increasing flexibility of the model will yield
diminishing returns on lower bias, but will cause variance to rise
dramatically as the model begins overfitting the data, corresponding to an
increase in expected test MSE. This is what gives the plot a characteristic
U-shape. In reality, the true $f$ is not known, so we have no way to compute
the expected test MSE, the bias, or the variance. It is still important to be
aware of them as they form the theoretical underpinning of certain model
behaviour. We do however have methods to estimate the test MSE, which we will
see shortly.

####  1.1.3 Resampling

Continuing with the theme of exploring model error, we now turn our attention
to a fundamental problem in statistics: sample size. The situation of not
having “enough” members in a sample in very prominent, and attempts to rectify
some of the issues it presents is of great interest. One occurrence of this
issue is known as the _curse of dimensionality_. While this ominous term very
generally refers to the differences between the properties of high and low-
dimensional data, a particular instance as it relates to sampling is where
increasing the number of inputs (dimension) causes an exponential increase in
the “size” or “volume” of the input space. What ends up happening is that a
correspondingly exponential increase in sample size is needed to accommodate
this increase in volume, as we would like for our model to have seen many
different possible combinations of inputs (otherwise we risk having high
variance).

Consider as before a partition of the samples into a training and a testing
set. The way that the partition is formed is random in an effort to mitigate
any implicit biases such as the way the data was collected. But what if our
random choices ends up selecting a subset that is not very representative of
the sample as a whole? The techniques we are going to describe are known as
_resampling methods_. These involve drawing many different partitions of our
sample, and fitting the model each time on the new training sets in order to
gain more accurate information about the performance of the model.
Essentially, this is a sort of smoothing process that seeks to lessen the
obscuring effect that a random partition can have on our measurement of the
test error. We are going to look at two very common methods for resampling in
statistical learning. The first method is called _cross-validation_ (CV). It
works by holding out on a subset of the observations during each iteration to
use as the testing set, while the remaining observations are used for
training. There are a few approaches to implementing cross-validation we will
discuss, each of which is defined by the size of the holdout set.

The first and simplest approach is the _validation set_ approach. During each
resampling, we simply take half of the observations for the training set, and
the other half for the validation (testing) set. As is typical, the selections
for each set is random after deciding on the proportions. Once this is
completed, the model can be fit to the training data and the _validation error
rate_ can be computed as an estimate for the test error rate. This is
typically done using the MSE function. As it stands, the validation set
approach performs poorly compared to its modified counterparts, as the
validation error rate can be highly variable. In practice, statistical models
tend to perform worse when trained with fewer observations. What results is a
tendency to overestimate the test error rate of the data.

On the other side of the spectrum, we have _leave-one-out cross-validation_
(LOOCV). During each iteration of the resampling, we again randomly partition
the sample. This time however, only a single observation-predictor pair is
selected for the validation set, say $(\textbf{x}_{1},\,y_{1})$, which leaves
the rest for training. We fit the model on the $n-1$ observations in the
training set and then compute the test MSE for this iteration, which is given
by $\text{MSE}_{1}=(y_{1}-\hat{y}_{1})^{2}$. This process is repeated, each
time selecting a different observation for the validation set and obtaining
another estimate for the test error. After exhausting all $n$ unique
validation sets, we take the average of each of the $n$ MSE terms as the LOOCV
estimate for the test error:

$\text{CV}_{n}\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{1}{n}\sum_{i=1}^{n}\text{MSE}_{i}.$
(1.11)

As each $\text{MSE}_{i}$ term only makes its measurement on a single
observation and the model is trained on nearly the entire sample set, it
follows that each MSE calculation is an approximately unbiased estimate for
the true test error. Mathematically, this can be shown by first assuming our
model is correct (a big assumption typically made for theoretical
calculations) and fixing an observation-predictor pair
$(\textbf{x}_{0},y_{0})$. Then we have that $E[y_{0}]=\hat{y}_{0}$, and using
the fact that $E[\varepsilon]=0$, it can be shown that

$\displaystyle E[\text{MSE}_{0}]$
$\displaystyle=E[(f(\textbf{x}_{0})+\varepsilon-\hat{f}(\textbf{x}_{0}))^{2}]$
$\displaystyle=E[\varepsilon^{2}]$ $\displaystyle=\text{Var}(\varepsilon)$
(1.12)

and thus

$\text{Bias}(\text{MSE}_{0})=E[\text{MSE}_{0}]-(\text{True Test Error})=0.$
(1.13)

However, there is a large variance associated with this estimation. It can be
reduced slightly when we compute the overall estimate $\text{CV}_{n}$, but
this comes at the cost of added computational time for fitting a model $n$
times. The result of the decreased bias is that the LOOCV estimate doesn’t
overestimate the test error as much as the validation set approach. Moreover,
as the validation sets are picked systematically, the LOOCV method has no
randomness in its partitioning, versus the validation set approach which
yields different results each time it is applied. The LOOCV method is very
general, and can be paired with nearly any kind of predictive model.

As mentioned, one drawback to LOOCV is that fitting a model $n$ times can be
computationally expensive —especially when $n$ is large or the model is
complex. The next method attempts to deal with this drawback, along with
finding a more desirable balance of bias and variance, while still attempting
to capture the essence of LOOCV.

######  Remark.

In the case that one uses linear or polynomial regression, a remarkable
formula found in (James et al., 2021, ch. 5) computes the LOOCV estimate for
the test error $\text{CV}_{n}$ with the same computational cost as one model
fit, as opposed to $n$ model fits.

The final approach we will look at, called _k-fold cross-validation_ , is one
generalization of LOOCV that performance-wise inhabits the space between the
validation set approach and LOOCV. For a natural number $k$, we randomly
divide the set of $n$ observations into $k$ groups (called _folds)_ of roughly
equal sizes. Choosing one fold at a time for the validation set, we fit our
model to the $k-1$ remaining folds, and then compute the MSE using the
validation fold. This yields $k$ estimates for the test error, and we use the
average of these quantities as the $k$-fold CV estimate:

$\text{CV}_{k}\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{1}{k}\sum_{i=1}^{k}\text{MSE}_{i}.$
(1.14)

Notice that LOOCV is a special case of this approach, corresponding to
$n$-fold CV. In practice, it is common to see 5 or 10-fold implementation, 5
for smaller data sets (say $n<100$) and 10 for large data sets ($n>1000$).
These values tend to strike a good balance between training with lots of
observations, while leaving enough for testing to get an adequate estimation
of performance. As suggested previously, the $k$-fold method performs
somewhere between the other, more extreme methods. While LOOCV tends to have
lower bias due to the number of training observations being maximal, the
viewpoint of modern theory is that in practice, it can be better to have a
biased estimate with smaller variance. Recalling the bias-variance trade-off,
LOOCV will indeed have higher variance that $k$-fold CV. This can be further
explained due in part to the high positive correlation between the outputs of
LOOCV. Notice that the training sets in LOOCV are nearly identical from
iteration to iteration, and hence knowing the performance of a model trained
on one gives us a good idea of the performance of the same model trained on
the others. In contrast, outputs in $k$-fold CV are less positively correlated
due to having more varied training sets, and hence a comparatively smaller
variance. This can be expressed mathematically using the identity for the
variance of the sum of random variables. Let $F_{i}$ and $F_{j}$ represent two
different folds with randomly chosen members. Then

$\text{Var}(F_{i}+F_{j})=\text{Var}(F_{i})+\text{Var}(F_{j})+2\text{Cov}(F_{i},F_{j}).$
(1.15)

We see that the variance becomes smaller when the covariance term is closer to
0. As a final note, it is clear to see that $k$-fold CV requires far less
computational time to execute when compared to LOOCV. As important as test
error is to assessing models, its actual value may not always be of interest.
If we want to identify a method or level of flexibility that results in lowest
test error, it really doesn’t matter what this value is. Rather, we want to
know the _location_ of the minimum test error in terms of level of
flexibility, either across multiple curves (in the flexibility-MSE plane)
representing different methods, or along the same curve representing one
method.

The second, very general method for resampling we will briefly mention is the
_bootstrap_ method. It is a widely employed tool for quantifying the accuracy
of estimators and statistical methods. It finds many use cases, especially
when these measurements are quite difficult to compute or statistical software
doesn’t compute them automatically. The outline for how it is performed is as
follows: Choose $n$ observations from the sample set of size $n$ with
replacement, meaning the same observation may be chosen multiple times, or not
at all. This is called the bootstrap data set. With this new data set, we make
the measurement or calculate the statistic that we ultimately want to
estimate. This process is repeated many times, and the mean of theses values
is used as the bootstrap estimate for its true value. The effect of resampling
using the bootstrap method is that it simulates the process of gathering new
samples from the population.

###  1.2 Classification Problems

Next, we will touch on some basic supervised statistical learning practices as
it pertains to classification problems. Recall that the distinction with
classification problems is the observed responses are placed into categories
rather than given a continuous numerical value. Most of the concepts from
section 1.1 on regression problems remain largely the same, with differences
having a natural “discrete” analogue. We will mainly touch on how some of the
tools change in this setting, as the core statistical principles are
essentially invariant to the domain of the observations.

Suppose that we wish to estimate $f$ on the basis of a subset of our
observations-predictor pairs as we did in the previous section, except now the
observations $y_{i}$ are _qualitative_. We call the estimator $\hat{f}$ for
$f$ a _classifier_ , as it takes a vector containing $p$ features as input,
and outputs a class label prediction (typically represented with an integer).
The function we will use to quantify the empirical error introduced by
$\hat{f}$ is defined by the following map:

$(\textbf{x}_{i},\,y_{i})_{i=1}^{m}\mapsto\frac{1}{m}\sum_{i=1}^{m}\textbf{I}(y_{i}\neq\hat{y}_{i}).$
(1.16)

The function $\textbf{I}(y_{i}\neq\hat{y}_{i})$ is called the _$0\text{-}1$
loss function_, where I denotes the indicator variable that outputs 1 if the
input is true and 0 otherwise. Thus, the function defined by the mapping in
equation 1.16 can be used to compute the ratio of incorrect classifications
made by $\hat{f}$. The _training error rate_ replaces the training MSE from
section 1.1, and is computed using the mapping in equation 1.16 with the
training set as input. Similarly, the _test error rate_ replaces the test MSE
and is found by computing the mapping in equation 1.16 on the test set. Once
again, we are most interested in reducing the test error rate rather than the
training error rate, as our goal is to make predictions for new data.

How can we use the training data to find a classifier? In theory, what is
known as the _Bayes_ classifier would be optimal in the sense that on average,
it is the classifier that minimizes the test error rate. It works simply by
assigning each observation to the class that it most likely belongs to given
its inputs. That is, a predictor $\textbf{x}_{0}$ is assigned to the class $j$
that maximizes the conditional probability
$P(y=j\mid\textbf{x}=\textbf{x}_{0})$. If we have $r$ possible
classifications, we could define the Bayes classifier by the mapping

$\textbf{x}_{0}\mapsto\operatorname*{arg\,max}_{j\in\\{1,\cdots,r\\}}P(y=j\mid\textbf{x}=\textbf{x}_{0}).$
(1.17)

A derivation of the Bayes classifier can be found in (Hastie et al., 2009).
When there are only two possible classifications, the set of points in the
input space for which $P(y=j\mid\textbf{x}=\textbf{x}_{0})=\frac{1}{2}$ is
called the _Bayes decision boundary_. This boundary splits the input space
between the two classes, although the set of points representing a class need
not be connected. For the points satisfying
$P(y=j\mid\textbf{x}=\textbf{x}_{0})>\frac{1}{2}$, the classifier will put
them into class _A_ , and for points satisfying
$P(y=j\mid\textbf{x}=\textbf{x}_{0})<\frac{1}{2}$, it will put them into class
_B_. One can visualize this by imagining a plane where each axis represents a
particular input (i.e. the input space is two dimensional), then splitting the
plane with a curve (which will be the decision boundary). We can imagine the
members of class _A_ inhabiting one side of the curve, and the members of
class _B_ inhabiting the other side. Do note that the Bayes classifier is the
best at minimizing test error on average, and is not perfect. With real data,
some of the observations may lie on the “wrong side” of the decision boundary.
Moreover, in general we won’t always be able to cleanly divide the
observations into distinct groups. The Bayes error rate at a point
$\textbf{x}_{0}$ is given by $1-\max_{j}P(y=j\mid\textbf{x}=\textbf{x}_{0})$,
and the overall Bayes error rate is $1-E[\max_{j}P(y=j\mid\textbf{x})]$, where
the expectation is taken with respect to x. The Bayes error rate is a form of
irreducible error, just as $\textup{Var}(\varepsilon)$ was in section 1.1. As
mentioned previously, the Bayes classifier remains an unobtainable theoretical
best. This is due to the conditional probability found in its definition,
which cannot be computed as we do not know the conditional distribution of $y$
given x. Next, we give an example of a simple method that estimates this
conditional probability.

Given a test observation $\textbf{x}_{0}$, it is quite natural to base the
prediction for its classification on nearby training observations. This is the
core idea behind the _K_ -Nearest Neighbours (KNN) classifier. First, we
decide on how many neighbours $K$ can influence classification. Next, it finds
the closest _K_ points to $\textbf{x}_{0}$. This requires a notion of
distance, which in this context we simply take as the regular Euclidean
distance. We label this set of neighbours with the notation
$N_{K}(\textbf{x}_{0})$. It then estimates the probability that
$\textbf{x}_{0}$ is in the _j_ th class using the following as an
approximation:

$P(y=j\mid\textbf{x}=\textbf{x}_{0})\approx\frac{1}{K}\sum_{\textbf{x}_{i}\in
N_{K}(\textbf{x}_{0})}\textbf{I}(y_{i}=j).$ (1.18)

Now, Bayes rule of classification can be applied, again picking the class for
which the estimated probability is highest. KNN classification tends to work
quite well in practice. One can adjust the flexibility of the KNN classifier
by simply choosing different values for $K$. Remember that we are most
interested in minimizing the test error and not so much the training error.
Thus to achieve the best results, we can apply the method numerous times while
varying $K$ to find the one in which test MSE is minimized.

Finally, we may want to employ resampling methods to gather more information
from our samples. Cross validation, as explored in section 1.1, works in much
the same way as we saw before. Only this time, our estimate for the test error
uses the training error rate found in equation 1.16 instead of MSE.

##  Chapter 2 Support Vector Machines

Given some characteristics of an object, what from those characteristics
differentiates it from other, similar objects? For example, how does your
email provider decide whether a given message belongs in your inbox, or in the
spam folder? How does a bank determine which transactions are fraudulent?
Based on patient symptoms and data collected from samples, can we say with any
level of confidence if a tumor is cancerous or benign? Classifications
problems are ubiquitous, and having a general, easy to implement, and
efficient way of solving these problems is crucial for many aspects of modern-
day life. At the end of chapter 1, we briefly touched on an approach to
classification known as $K$-nearest neighbours that classified observations
based on the prominent class amongst similar observations. This was in an
attempt to approximate the optimal Bayes classifier, which itself is pulling
from the true conditional distribution of the predictor and response
variables. We will now explore an alternative approach to classification using
what are known as _support vector machines_. In essence, support vector
machines make one assumption of the given data, which is that there is level
of separability of the classes. This differentiates them from the $K$-nearest
neighbours classifier as no assumption about a distribution is made. Rather,
they seek to partition classes from a geometric point of view. Picture it as
drawing a line or curve through the data, keeping an equal distance from both
classes, instead of focusing on particular points and factoring in the local
behaviour of other points. We mainly concern ourselves with the case of binary
classification problems, which is the setting in which support vector machine
were created and often employed. There have been attempts to extend the theory
to a greater number of classes, which will be briefly mentioned at the end of
the chapter.

The chapter will proceed roughly as follows. First, we will start with a light
introduction to the concept of classifying observations based on “separation”
by assuming that our observations can be perfectly separated into two distinct
classes by a linear function called a _hyperplane_ (section 2.1.1). In
practice, this method by itself is not very useful since classes tend to
overlap, so we will then move to the _support vector classifier_ that attempts
to deal with observations that seem to cross what we would perceive to be the
decision boundary (section 2.1.2). This support vector classifier is but a
special case of the general _support vector machine_ , or _SVM_ (section 2.2).
It is not always true that a linear decision boundary accurately describes
where one class should end and another should begin at every point in space.
In general, SVMs replace such a “dividing line” with one or more curves that
can bend around data points in the space. It will be here that we discuss a
very special class of functions known as “kernels”, which measure how similar
two observations are in some very specific way. We will assume that our input
space is the real vector space $\mathbb{R}^{p}$. This presentations is adapted
from James et al. (2021, ch. 9), Theodoridis (2015, ch. 11), and Moguerza and
Muñoz (2006).

###  2.1 Hyperplane Classification

For a $p$-dimensional space (eg. the real vector space $\mathbb{R}^{p}$), a
_hyperplane_ is a $p-1$ dimensional affine subspace. The term “affine” means
that the space need not pass through the origin, as a typical vector subspace
would, but must still retains its closure properties. For example, a
hyperplane $H$ in $V$ could be written as $H=\textbf{v}_{0}+W$ for some vector
$\textbf{v}_{0}$ in $V$, and subspace $W$ of $V$. In the regular $xy$-plane, a
hyperplane would be any straight line dividing the plane. Similarly in
$\mathbb{R}^{3}$, a hyperplane would look like a two dimensional flat sheet
dividing the space, or in other words a “copy” of $\mathbb{R}^{2}$ inside
$\mathbb{R}^{3}$. More precisely, we can parameterize a hyperplane so that it
can be written as the set of points $\textbf{x}=(x_{1},\cdots,x_{p})^{T}$ in
the space satisfying

$\textbf{w}^{T}\textbf{x}+b=w_{p}x_{p}+\cdots+w_{1}x_{1}+b=0$ (2.1)

for some fixed vector $\textbf{w}=(w_{1},\,\cdots,w_{p})^{T}$ and constant
$b$. The intuitive picture of the dividing nature of the hyperplane can be
formulated mathematically by saying the set of points satisfying
$\textbf{w}^{T}\textbf{x}+b>0$ are on one side of hyperplane, while the set of
points satisfying $\textbf{w}^{T}\textbf{x}+b<0$ are on the other side. When
such a hyperplane exists, a test data point $\textbf{x}_{i}$ can be classified
based on the sign of $\textbf{w}^{T}\textbf{x}_{i}+b$. Thus, we represent the
two classes with $y\in\\{1,-1\\}$. The quantity
$\displaystyle\left\lvert\textbf{w}^{T}\textbf{x}_{i}+b\right\rvert$ also
carries with it useful information. The larger its value, the further away
from the hyperplane the observation will lie, and hence the more confident we
are in the correctness of the classification.

Next we establish some terminology and common assumptions regarding hyperplane
classification. Recall that the perpendicular distance from a point
$\textbf{x}_{0}$ to a hyperplane is given by the quantity

$\frac{|\textbf{w}^{\,T}\textbf{x}_{0}+b|}{\|\textbf{w}\|}.$ (2.2)

As the hyperplane is invariant to scaling, if $\textbf{x}_{0}$ is the closest
training data point to the hyperplane, we can arrange for
$|\textbf{w}^{T}\textbf{x}_{0}+b|=1$ so that the perpendicular distance is
$1/\|\textbf{w}\|$. This “normalizing” of the hyperplane will be assumed
henceforth. The goal is to build an _optimal margin classifier_ , meaning the
classifier doesn’t favour one class over the other. The optimal hyperplane
will have at least two points that are “closest” to it, that is, two points
with perpendicular distance $1/\|\textbf{w}\|$ to the hyperplane. These points
lie on either side of the hyperplane and represent the closest points in each
training data class to the hyperplane. Thus, the distance between the classes,
called the _margin_ , is $2/\|\textbf{w}\|$. The two hyperplanes defined by
$\textbf{w}^{T}\textbf{x}+b=\pm 1$ are together called the _boundary of the
margin_. The goal can be stated more concretely as maximizing the margin, or
equivalently, minimizing $\|\textbf{w}\|$.

####  2.1.1 Separable Observations and the Maximal Margin Classifier

With the general idea of hyperplane classification in mind, we now begin with
the simplest approach. This occurs when the data is _linearly separable_ ,
meaning there exists a hyperplane such that one side of it contains only one
class and the other side contains only the other. Our goal is to construct an
optimization problem that can find a vector w and constant $b$ parameterizing
a hyperplane that sits perfectly between the two classes, meaning that no
training point is misclassified and that the plane is equidistant from the
nearest points in each class. This process will serve as a template for each
successive, more general classification procedure. For the reader looking for
a more thorough exploration of the theory, a standard source for (convex)
optimization theory is Boyd and Vandenberghe (2004). Outlines of the solutions
to the optimization presented here can be found in Theodoridis (2015, ch. 11).

For a fixed training set $(\textbf{x}_{i},y_{i})_{i=1}^{m}$, we want to train
a model $\hat{f}$ (also called a decision rule or discriminant function in the
context of classification problems) that minimizes the empirical loss from
misclassification

$\hat{f}\mapsto\sum_{i=1}^{m}L(y_{i},\hat{f}(\textbf{x}_{i}))$ (2.3)

where $L$ is a _loss function_. We saw an example of a loss function in
section 1.2, which was the 0-1 loss function. A different loss function will
be used to construct SVMs for a few reasons. While the 0-1 loss function is
very simple and quite natural, the presence of jump discontinuities make it
incompatible with many optimization methods that require levels of
differentiability. It is possible to find a work around if a loss function has
points where it is only continuous by using a smooth function to approximate
the loss function, but this isn’t viable for the 0-1 loss function since it is
discrete. Instead, the _hinge loss function_ is used:

$L_{H}(y_{i},\hat{f}(\textbf{x}_{i}))\stackrel{{\scriptstyle\textup{def}}}{{=}}\max(0,\,1-y_{i}\hat{f}(\textbf{x}_{i})).$
(2.4)

Since we know that the sign of a normalized hyperplane will be used for
$\hat{f}$, it is immediately apparent how $L_{H}$ will penalize
classifications. A penalty is incurred whenever
$y_{i}\hat{f}(\textbf{x}_{i})\leq 1$, corresponding to a violation of the
margin. Otherwise, the point has been correctly classified, so no penalty is
added. What this loss function does is give us a _sparse solution_ , since the
only points that are going to influence the hyperplane are the points lying on
the boundary of the margin, which generally consists of a very small subset of
the training data. This is desirable as it makes the classifier resilient to
outliers and it gives them an advantage when dealing with large amounts of
data. Thus, the problem of designing an optimal classifier with the desired
properties can be recast as finding the minimum of the following cost
function:

$(\textbf{w},b)\mapsto\frac{1}{2}\|\textbf{w}\|^{2}+C\sum_{i=1}^{n}L_{H}(y_{i},\,\textbf{w}^{T}\textbf{x}_{i}+b).$
(2.5)

Let us dissect the meaning of this cost function. First, as mentioned at the
beginning of section 2.1, we want to maximize the margin subject to some
constraints. This is equivalent to minimizing $\|\textbf{w}\|^{2}$, as
squaring the norm doesn’t change the location of the minimum. The squaring is
done to simplify the computation. The factor of $\frac{1}{2}$ is to cancel the
factor of 2 that results from differentiation. This scaling has no effect on
the second term, as it itself is scaled by a user-controlled hyperparameter
$C>0$ that controls the weight of the accumulated losses. In the case of
linearly separable data, minimizing this cost function is equivalent to the
following optimization problem

minimize $\displaystyle\frac{1}{2}\|\textbf{w}\|^{2}$ (2.6) subject to
$\displaystyle y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)\geq
1,\,i\in\\{1,2,\,\cdots,n\\}.$

The solution is given by the finite linear combination

$\textbf{w}=\sum_{i\in\mathcal{S}}\lambda_{i}y_{i}\textbf{x}_{i}$ (2.7)

where $\mathcal{S}$ is the set of indices of the nonzero Lagrange multipliers
$\lambda_{i}$ obtained from the dual representation of the problem. Each
Lagrange multiplier is associated to a constraint
$y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)\geq 1$, and remarkably, the nonzero
Lagrange multipliers correspond to precisely those constraints for which
$y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)=1$, which as described earlier are the
closest points to the hyperplane lying on the boundary of the margin. These
training points are called _support vectors_ , as they completely determine
the hyperplane. In effect, the solution is able to ignore all but the most
essential observations. Now, _b_ can be found by taking one of the constraints
for which $\lambda_{i}\neq 0$ ($i\in\mathcal{S}$) and solving

$y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)=1$

for $b$. In practice, different constraints may have slightly different values
of $b$ due to rounding errors, so for the sake of numerical stability, the
average over the set $\mathcal{S}$ is taken. Therefore, our classifier
$\hat{f}$ predicts the class $\hat{y}\in\\{1,-1\\}$ for a test observation x
using the formula

$\hat{f}(\textbf{x})=\operatorname{sgn}\left(\sum_{i\in\mathcal{S}}\lambda_{i}y_{i}\left(\textbf{x}^{T}\textbf{x}_{i}\right)+b\right)$
(2.8)

where $\operatorname{sgn}$ is the sign function that outputs 1 if the input is
positive, and -1 if the input is negative. The function in equation 2.8 is
what is known as the _maximal margin classifier_.

####  2.1.2 Non-Separable Observations and the Support Vector Classifier

Next, we describe the procedure for dealing with two classes that are not
linearly separable. We can imagine a scenario where their are two
concentrations of observation in our data, but the border between one and the
other is not entirely clear due to noise. Some observation appear to be on the
“wrong side” of the plot. For now, we are still only considering the case
where a linear classifier is most appropriate. To codify this, we need an
optimization problem with a built-in level of tolerance for how many training
observations can be misclassified, and by how much. This is accomplished
through the following optimization problem:

minimize
$\displaystyle\frac{1}{2}\|\textbf{w}\|^{2}+C\sum_{i=1}^{n}\varepsilon_{i}$
(2.9) subject to $\displaystyle y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)\geq
1-\varepsilon_{i},$ $\displaystyle\varepsilon_{i}\geq
0,\,i\in\\{1,2,\,\cdots,n\\}.$

This problem will maximize the margin around the decision boundary, while
ensuring a minimum number of misclassifications. In this case, we say that
there is a soft-margin, meaning that observations may cross it, but at a cost.
The user-defined variable $C\geq 0$ is a hyperparameter that acts like a
budget. It give an upper bound on how many violations we are willing to
accept, and how severely. When $C=0$, the problem is reduced to equation 2.6
corresponding to the separable case. When $C>0$, then at most $\lfloor
C\rfloor$ observations are allowed to cross the hyperplane, and consequently
be misclassified. An increase in $C$ will widen the margin, while a decrease
in $C$ will shrink the margin. To tie this idea back to our discussions in
chapter 1, $C$ is controlling the bias-variance trade-off of the method. A
larger $C$ allows for more violations of the margin, which generally results
in a high bias but a low variance. Conversely, a smaller $C$ allows for fewer
violations, so the classifier will be highly fit to the training data, which
suggests a low bias but a high variance.

The solution to this optimization problem takes the same form as in the
separable case. The only difference is which constraints have a nonzero
Lagrange multiplier. Let us briefly describe the three types of observations
and the values of their associated Lagrange multipliers. A predictor
$\textbf{x}_{i}$ whose constraint satisfies
$y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)\geq 1$ has been classified correctly
and lies on the boundary or outside of the margin. These points incur no
penalty, so $\varepsilon_{i}=0$ and $\lambda_{i}=0$. Now, if
$0<y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)<1$, then $\textbf{x}_{i}$ has been
classified correctly, but lies inside the margin. These points are said to
violate the margin, and incur a penalty $0<\varepsilon_{i}<1$. Their
associated Lagrange multiplier in this case are nonzero. Finally, if
$y_{i}(\textbf{w}^{T}\textbf{x}_{i}+b)\leq 0$, then $\textbf{X}_{i}$ has been
incorrectly classified. These points incur a penalty $\varepsilon_{i}>1$, and
again have a nonzero Lagrange multiplier. In this context, not only are the
points on the boundary of the margin support vectors, but so are all points
within the margin and outside of the margin on the wrong side of the
classifier. These classifiers are commonly called _support vector classifiers_
, although the terminology in the literature may simply call them support
vector machines for simplicity. We reserve this title for the more general
classifier in the next section.

###  2.2 Nonlinear Classification

In the previous section, we saw how a classifier can be constructed using a
linear combination of dot products with support vectors. We typically think of
the dot product as a similarity measure, taking into account the angle between
two points through the identity
$\textbf{x}^{\,T}\textbf{y}=\|\textbf{x}\|\|\textbf{y}\|\cos{\theta}$, with
$\theta$ being the angle between the points. A dot product is a particular
example of an _inner product_ , which is a function used to quantify this same
type of similarity in general vector spaces. This fact is used to construct
support vector machines that solve the task of classifying data that is not
only nonseparable, but also where a nonlinear decision boundary will lead to a
lower test error.

####  2.2.1 Reproducing Kernel Hilbert Spaces

Hilbert space is a very important generalization of Euclidean space found
throughout Mathematics, Physics, and Statistics. It is a special type of
vector space endowed with many of the familiar notions of distance, angle,
projection, and “continuity”. More formally, a Hilbert space $\mathbb{H}$ is a
_complete inner product space_. This means that it is a complete vector space
(possibly infinite dimensional) over a field of scalars (typically the real or
complex numbers) equipped with an inner product. An inner product is a scalar-
valued function that gives a way to quantify similarities between vectors.
Now, inner products induce a norm (given by the square root of the inner
product), which is interpreted as a measure of length of vectors, and can be
used to define distance between vectors. The final piece of the definition is
_completeness_ , meaning that any sequence of points in the space that get
arbitrarily close together with respect to the norm must also get arbitrarily
close (or converge) to an element in the space. Sequences with this property
are known as _Cauchy sequence_. This can intuitively be thought of as the
space having “no holes”. A familiar example of a complete space is the real
numbers with the Euclidean metric. This space can be viewed as the completion
of the space of rational numbers, which is famously incomplete, not containing
fundamental constants such as $\sqrt{2},\,\pi,\,\textup{and }e$. Here we are
only interested in Hilbert spaces consisting of real-valued functions defined
on a subset of $\mathbb{R}^{n}$, plus an extra bit of structure described in
the following definition.

######  Definition.

A Hilbert space of real-valued functions $\mathbb{H}$ on a set
$U\subset\mathbb{R}^{n}$ is called a Reproducing Kernel Hilbert Space (RKHS)
if there is a real-valued function $k:U\times U\to\mathbb{R}$ with the
following properties:

  1. i

. For every point $x$ in $U$, the function
$k_{x}(\cdot)\stackrel{{\scriptstyle\textup{def}}}{{=}}k(\cdot,x)$ is in
$\mathbb{H}$.

  2. ii

. For every function $f$ in $\mathbb{H}$ and point $x$ in $U$, we have that
$f(x)=\langle f,\,k_{x}(\cdot)\rangle$.

Property ii. is called the reproducing property of k, and k is called the
reproducing kernel of $\mathbb{H}$.

Notice that the function $k$ and the inner product can be used to evaluate any
function in the space at any point. A consequence of the reproducing property
of $k$ is that for any points $x,\,y\in U$,

$k(x,y)=k(y,x)=\langle k_{x}(\cdot),k_{y}(\cdot)\rangle.$ (2.10)

The mapping $x\mapsto k_{x}(\cdot)$ is called the _feature map_ , and in
combination with the above inner product, will ultimately allow us to access
the rich structure of the RKHS _implicitly_. The procedure for this is as
follows: Starting with a set containing training data
$(\textbf{x}_{i},\,y_{i})_{i=1}^{m}$, the feature map is employed implicitly
through equation 2.10 by evaluating the kernel function
$k(\textbf{x}_{i},\textbf{x}_{j})$ on the left-hand side. This has the effect
of measuring the similarity of the two functions in the RKHS associated to the
training points, which can be leveraged in classification problems as the
training points might not be linearly separable in the input space, but the
associated functions might be in the RKHS. We know that this approach is
potentially viable due to Cover’s theorem, which proves that data becomes
arbitrarily separable as the number of features becomes arbitrarily large
(Cover (1965), Theodoridis (2015)). In the context of support vector machines,
this RKHS is called the _feature space_. The feature space typically has a
much higher dimension than the input space (sometimes an infinite number of
dimensions), and the feature map arrives at these extra dimensions by
effectively building new predictors out of the ones in the input space. Next,
we will highlight some relevant technical details of RKHS’s.

######  Proposition.

Let $\mathbb{H}$ be a RKHS on a set U with kernel k. Then

$\mathbb{H}=\overline{\operatorname{span}\\{k_{x}(\cdot):x\in U\\}}.$ (2.11)

What this result is saying is that any function in the RKHS can be generated
from the kernel function, either by a finite linear combination of
$k_{x_{i}}$’s, or by an infinite series. The idea behind why this is true is
that the only function orthogonal to all elements in the span of the $k_{x}$’s
is the zero function, due to the reproducing property. Then, one can use the
well-known result that any Hilbert space can be decomposed into the direct sum
of any closed subspace and its orthogonal complement (see Rudin, 1991, Theorem
12.4).

It is important to note that RKHS’s have a few different characterizations,
and one may encounter them from sources using an alternative definition as a
starting point. (Schölkopf et al., 2002) contains complete presentation of
these characterizations. We will briefly cover some of these in the remaining
paragraphs of this section. An equivalent (and more intuitive) formulation of
an RKHS is that every evaluation functional defined on the entire space is
continuous (Aronszajn, 1950). In other words, two functions that are close
with respect to the norm are also point-wise close throughout the underlying
set. This is an important property because functions in an RKHS are vectors in
and of themselves, independent of the numerical value they take at points in
their domain. Next, we present the definition of a very important class of
functions.

######  Definition.

Let U be any set. A symmetric real-valued function $k:U\times
U\rightarrow\mathbb{R}$ is said to be a positive definite kernel if for every
$n\in\mathbb{N}$, we have that

$\sum_{i=1}^{n}\sum_{j=1}^{n}a_{i}a_{j}k(x_{i},x_{j})\geq 0$ (2.12)

for any $n$ points $x_{1},\,\cdots,x_{n}$ in $U$ and $n$ real numbers
$a_{1},\,\cdots,a_{n}$.

The kernel matrix (also known as the _Gram matrix_) $\mathcal{K}$ of a kernel
$k$ with respect to a set of data points
$\textbf{x}_{1},\cdots,\textbf{x}_{p}$ in input space is a $p\times p$ matrix
whose entries consists of the kernel evaluated at the associated data points,
that is, the $ij$th entry contains the real number
$k(\textbf{x}_{i},\textbf{x}_{j})$. This is a common way of storing the values
of a kernel acting on a data set. The positive definite property of $k$ is
equivalent to the Gram matrix $\mathcal{K}$ being _positive semidefinite_ ,
meaning that $A^{T}\mathcal{K}A\geq 0$ for every set of data points in the
domain of $k$, and every $A\in\mathbb{R}^{p}$.

Another alternate characterization is that the reproducing kernel of an RKHS
is a symmetric positive definite (p.d.) kernel. In the case that _k_ is the
kernel of an RKHS, this property follows from the reproducing property and the
fact that inner products are p.d. kernels. Conversely, every symmetric p.d.
kernel induces a unique RKHS on its underlying set (see Aronszajn, 1950). This
can be particularly useful for building new kernels out of old ones, as one
can apply easier-to-prove properties of p.d. kernels to find new kernels. Some
properties include: the sum (and product) of two p.dk̇ernels is again a p.d.
kernel, and a composition of a p.d. kernel with any function is again a p.d.
kernel. A more extensive list can be found in Theodoridis (2015), and the book
Schölkopf et al. (2002) contains an entire chapter dedicated to the design of
kernels. The latter book also shows that new positive definite kernels can be
constructed using any inner product space: If $\varphi$ is a function mapping
a set $U$ into an inner product space, then
$k(x,y)=\langle\varphi(x),\varphi(y)\rangle$ is a p.d. kernel. Finally, due to
Mercer, we know that a continuous symmetric p.d. kernel has a series expansion
of the form

$k(x,y)=\sum_{i=1}^{\infty}\alpha_{i}^{2}\varphi_{i}(x)\varphi_{i}(y)$ (2.13)

where each $\alpha_{i}$ is a nonnegative real number, and each $\varphi_{i}$
is a real-valued function with some special properties. The details can be
found in Shawe-Taylor et al. (2004) or Mercer (1909), and the converse
statement can be found in Schölkopf et al. (2002). This gives yet another way
to check whether a given function is the kernel of an RKHS. We can write
equation 2.13 more compactly as $k(x,y)=\varphi(x)^{T}\varphi(y)$, where
$\varphi(x)=(\alpha_{1}\varphi_{1}(x),\,\alpha_{2}\varphi_{2}(x),\,\cdots)$.
This convenient notation shows us that Mercer kernels act as a “dot product”
in the Hilbert space.

How does this connect back to statistical learning and more specifically to
support vector machines? As a partial answer to the first question, there is a
famous theorem known as the _representer theorem_ that very roughly states
that any minimizing function of a minimization task can be represented by a
finite linear combinations of functions $k_{\textbf{x}_{i}}$ where $k$ is the
reproducing kernel belonging to the space of the minimizing function, and the
$\textbf{x}_{i}$’s are training observations (Moguerza and Muñoz, 2006). This
concept will come in handy later as we answer the second question.

####  2.2.2 Using Kernels for Nonlinear Classification

In more general applications, linear classification via hyperplane simply
doesn’t make for effective classification. If the decision boundary is highly
nonlinear, it may not even be possible to find a hyperplane that divides the
data in any meaningful way. When dealing with two dimensional data, a plot may
reveal that a linear classifier cannot accurately capture where one class
should end and the other should begin. Perhaps a parabola would result in a
better fit if they are still separated from each other. What if one class is
completely surrounded by another class? How do we come up with a classifier
that can capture this geometry? First, we fix a kernel $k$ to build our
support vector machine with. Selecting a kernel that minimizes the test error
is a nontrivial task, and often many kernels will be used in the fitting
process to find one that works best for the particular problem. We can then
solve the same optimization problem as in the nonseparable case, but instead
we replace the inner product in the dual representation of the problem with
our choice of kernel $k$ (this is the so-called _kernel trick_) (see
Theodoridis, 2015, ch. 11). Once again, the solution will yield a set of
nonzero Lagrange multipliers associated to the constraints, and we can solve
for $b$ in the same way as before. Thus for a choice of kernel $k$, the
corresponding decision function $\hat{f}$ can be written as

$\hat{f}(\textbf{x})=\operatorname{sgn}\left(\sum_{i\in\mathcal{S}}\lambda_{i}y_{i}k(\textbf{x},\textbf{x}_{i})+b\right).$
(2.14)

Throughout the process of solving the dual problem, as well as evaluating the
decision function, we never actually need to perform any computations in the
RKHS associated to $k$. All of the benefit brought by employing an RKHS can be
accessed simply through evaluating $k$, which in many cases is a function of a
familiar form. Some of the basic kernels used in practice are as follows. The
first example is called the _linear kernel_ , and is given by the formula
$k(\textbf{x},\textbf{y})=\textbf{x}^{\,T}\textbf{y}$. This is simply the dot
product in $\mathbb{R}^{p}$, and was the kernel used in section 2.1 to solve
the problem of classifying observations using a linear decision boundary. It
can readily generalize to the _polynomial kernel_ , given by
$k(\textbf{x},\textbf{y})=(\textbf{x}^{\,T}\textbf{y}+c)^{d}$, where _c_ and
_d_ are hyperparameters that control the weight of lower-order terms and the
degree of the polynomial, respectively. Naturally, the decision boundary that
results from a polynomial kernel will take a very similar shape to the graph
of a polynomial of the same degree. The Gaussian kernel, or _radial basis
function_ (RBF) kernel is a very popular choice as it incorporates Euclidean
distance, which in effect allows it to create highly nonlinear, “oval-shaped”
decision boundaries that bunch clusters of observations together. It is given
by the formula
$k(\textbf{x},\textbf{y})=\exp{(-\lambda\|\textbf{x}-\textbf{y}\|^{2})}$,
where $\lambda>0$ is a hyperparameter that controls the range of influence of
the support vectors. It is common to see $\lambda$ parameterized with
$1/(2\sigma^{2})$, where $\sigma$ is taken as the tuning parameter.
Interestingly, the feature space of the RBF kernel has an infinite number of
dimensions. Quite often, the hyperparameters found in kernels are chosen via
cross-validation (discussed in section 1.1.3).

Through rigorous testing, SVMs have been found to offer good generalization
performance (Theodoridis, 2015, ch. 11.10). What does this mean exactly? After
fitting a model to a training set and tuning hyperparameters to reduce test
error, one can bring in an entirely new data set not yet seen by the model,
and find the _generalization error_ , which measures the accuracy of
predictions. Empirical observations show that SVMs perform quite well in terms
of having a low generalization error. While sparsification helps to reduce the
computational requirement in practice, large data sets do pose a problem.
There can still be many support vectors, and thus the computational load can
still be high, but more importantly the training time scales poorly with
larger data sets. For example, during the optimization process, an $n\times n$
matrix (called a _Gram matrix_) is used for each possible comparison of
predictor variables. Its size scales quadratically, and along with more
technical details regarding the optimization problem, a major slow down in
training time occurs for large data sets. Hsu et al. (2003) provide a guide on
good practices when performing support vector classification. One of the
suggestions they make is to scale the entries of a Gram matrix before passing
it to a function that fits an SVM. This can greatly improve runtime
performance of the model fitting, as well as lead to accuracy improvements.
Typically, one would linearly scale the training data to the unit interval,
and then scale subsequent testing data using the exact same function. This
technique is very important when dealing with _graph kernels_ (chapter 3).

####  2.2.3 SVMs for Multiclass Classification

Although SVMs were not designed with multiclass classification in mind, there
are ways that they may be implemented in this setting. We will cover the two
most common approaches. The first approach is known as _one-versus-one
classification_. If we have $K>2$ possible classifications, then for any two
classes, the methods discussed in this chapter are used to construct an SVM
solely based on these classes. This requires $\binom{K}{2}$ separate SVMs to
cover all possible pairings of classes. Points are then classified according
to which class they appear in the most across each SVM. The second approach is
similar in spirit and is known as _one-versus-all classification_. This time,
for each fixed class, we treat the remaining $K-1$ classes as one large class,
and fit an SVM between the two “classes” as before. Points can then be
classified based on which of the $K$ decision functions has the largest
absolute value (in terms of equation 2.14, we would simply ignore the sign
function and take the absolute value of the inside function). As mentioned
previously in the context of hyperplanes, the larger the magnitude of this
inner function, the more confident we can be that the point has been
classified correctly.

##  Chapter 3 Graph Kernels

The success of the kernel trick discussed in chapter 2 as a solution to
nonlinear binary classification has sparked the investigation of more general
applications of kernels in machine learning. In this chapter, we discuss the
use of kernels to measure the similarity between _graphs_ , which are network-
like structures consisting of vertices and edges that are used to model
discrete data, such as molecules and compounds. We will mostly concern
ourselves with various types of graph kernels and experimental results, and
not necessarily be discussing them in conjunction with a learning method,
although the experimental results found in the literature are usually obtained
with a support vector machine.

In chemoinformatics, we are interested in the properties that chemical
compounds have. As there are an enormous number of compounds, it is impossible
to measure every property of every compound. Hence, it is useful to be able to
predict properties that a chemical compound has based on knowledge of other
compounds. This is done on the principal that similar compounds tend to have
similar properties (Brown, 2009). Similarity of course has no clear cut
definition, which forms the basis of research into graph kernels. As molecules
are naturally represented as undirected labeled graphs, researchers come up
with different measures that attempt to capture the similarities of these
graphs as a way of quantifying how similar two molecules are to each other.
Formulating these measure as positive definite kernels opens up the powerful
_kernel trick_ for classification, which allows support vector machines to be
nonlinear binary classifiers. This classification is how we infer similarity,
and hence how properties may be predicted.

The way Graph Kernels are usually constructed is based on the framework
introduced by Haussler (1999). Given two graphs $G_{1}$ and $G_{2}$, one can
view them as the composition of subgraphs with certain properties, such as
directed and undirected subgraphs, or paths, cycles, etc. (Ralaivola et al.,
2005). A graph kernel is then created from applying a positive-definite kernel
on each pairing of subgraphs from $G_{1}$ and $G_{2}$, and taking the sum.
Graph kernel designs are typically centred around one specific structural
characteristic, dependent on what the application calls for. This is due to
the fact that determining many graph properties without restriction are either
NP-hard, NP-complete, or thought to be either NP-intermediate or NP-hard. For
example, see Gärtner et al. (2003) for a proof that capturing all the
information from a graph via subgraph isomorphism is NP-Hard. Borgwardt and
Kriegel (2005) also showed a similar result for the “all-paths” kernel. In
practice, this places a limit on the expressivity of graph kernels as
similarity measures.

Section 3.1 is dedicated to establishing the basic definitions and terminology
regarding graphs, as well as an explanation of how kernels can use graph data
to measure similarity. Section 3.2 details the initial graph kernels that
began showing up in the early to mid 2000’s that were popular or influential
in the literature. Section 3.3 contains descriptions of kernels that are
slightly more complex and less general, but perform better than those of the
previous section. Section 3.4 outlines techniques that have been developed to
improve run time or classification accuracy of classes of kernels. Section 3.5
goes over common data sets and software implementations of graph kernels and
support vector machines found in experiments in the literature. Finally,
section 3.6 contains references to further readings, including more in depth
treatments of the kernels discussed here, as well as kernels relevant to
chemoinformatics that were not discussed here.

###  3.1 Introductory Concepts

####  3.1.1 Essential Graph Theory

To begin, we review some basic definitions from graph theory. For more on the
subject, see the standard textbook by Diestel (2017).

######  Definition 3.1.1.

A graph $G=(V,E)$ consists of two sets; a vertex set
$V=\\{v_{1},\cdots,v_{n}\\}$ containing $n$ vertices (or nodes) and an edge
set $E\subset\\{(v_{i},v_{j})\in V\times V:v_{i}\neq v_{j}\\}$ containing
pairs of distinct vertices.

Vertices are often thought of as representing the individual objects of
interest, while each edge is thought of as a one-way connection between two
objects. For example, we could represent a road network with a graph, where
the vertices are intersections and the edges are roads. We could also model a
molecule using a graph, taking the atoms to be vertices and the bonds between
the atoms as edges. Unless otherwise stated, we shall assume that an arbitrary
graph $G$ has a vertex set $V$ with $n$ vertices and an edge set $E$ with $m$
edges.

While intuitive, this representation cannot be readily implemented with
typical operations and algorithms. Throughout subsequent sections, specific
representations of graphs will be introduced as need. First, we look at a very
important representation (even outside our scope of machine learning) known as
the _adjacency matrix_.

######  Definition 3.1.2.

The adjacency matrix A of a graph $G$ is a $|V|\times|V|$ matrix where
$[A]_{ij}=1$ if $(v_{i},v_{j})$ is an edge in $G$ and 0 otherwise.

This representation is particularly useful as the underlying graph may be
examined and manipulated using well-understood techniques.

######  Definition 3.1.3.

A graph $G$ is said to be undirected if for each edge $(v_{i},v_{j})\in E$, we
have that $(v_{j},v_{i})\in E$.

If the above property does not hold for every edge, we have what is known as a
_directed graph_. All chemical graphs will be undirected, however directed
graphs are important for a specific class of graph known as _trees_ (section
3.3.2). If a graph is undirected, then its adjacency matrix will be symmetric.

######  Definition 3.1.4.

A walk of length m in a graph $G$ is a sequence of nodes $(v_{i})_{i=1}^{m+1}$
where $(v_{i},v_{i+1})$ is an edge in $G$ for every $i\in\\{1,\,\cdots,m\\}$.
If $i\neq j$ implies that $v_{i}\neq v_{j}$, then we say that
$(v_{i})_{i=1}^{m+1}$ is a path. A walk with distinct vertices except for the
start and end points is called a cycle.

The existence of a walk between vertices implies a level of connection between
them, depending on the context. For an undirected graph $G$, we say that $G$
is _connected_ if there exists a walk between any two vertices in $G$.
Chemical graphs will always be connected. If $v$ is a vertex in $G$, the
_degree_ of the vertex $v$, denoted by $\textup{deg}(v)$, is the number of
neighbours of $v$, that is, the number of distinct walks of length one
beginning at $v$. Finally, we have the notion of a _labeled graphs_ (for
example, see Mahé et al. (2004)).

######  Definition 3.1.5.

A labeled graph $G=(V,\,E,\,L)$ is a graph equipped with a labeling function
$l:V\cup E\rightarrow A$, where $A$ is a set of labels.

A graph’s labeling function provides a way of differentiating both vertices
and edges, and assigning properties. Graphs of molecules will always be
labeled: the vertices represent atoms, and are labeled by their chemical
symbol. The edges represent covalent bonds, and are sometimes labeled to
indicate a type of covalent bond (single, double, etc.). In general, when
referring to a graph $G$, we will assume it to be labeled, undirected, and
connected, unless otherwise stated. The _label sequence_
$(l(v_{1}),l(v_{1},v_{2}),\cdots,l(v_{m},v_{m+1}),l(v_{m+1}))$ associated with
a walk $(v_{i})_{i=1}^{m+1}$ is a sequence containing the labels of every
vertex and edge in the walk, in respective order (Mahé et al., 2004). For
simplicity, we will sometimes refer to a label sequence associated with a walk
as simply a labeled walk.

####  3.1.2 Convolution kernels

A very general construction of kernels on graphs is based on the framework of
the _convolution kernel_ , introduced by Haussler (1999). We present the
definition of the convolution kernel based on the form found in Vishwanathan
et al. (2010).

######  Definition 3.1.6.

Let $\mathcal{X}$ be a set of discrete objects, and let $x\in\mathcal{X}$ be
an object with decomposition into $N$ components
$\textup{{x}}_{d}=(x_{1},\,\cdots,x_{N})$. Define a function
$R(x,\textup{{x}}_{d})$ that outputs true if $\textup{{x}}_{d}$ is a valid
decomposition of $x$, and false otherwise. If $R^{-1}(x)$ denotes the set of
all valid decomposition’s of a discrete object $x$, and $k_{i}$ is a kernel
measuring the similarity between the ith components of $x$ for each
$i\in\\{1,\,\cdots,N\\}$, then the convolution of the kernels
$k_{1},\,\cdots,k_{N}$ is defined as

$k(x,x^{\prime})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{\textup{{x}}_{d}\in
R^{-1}(x)}\sum_{\textup{{x}}_{d}^{\prime}\in
R^{-1}(x^{\prime})}\prod_{i=1}^{N}k_{i}(x_{i},\,x^{\prime}_{i}).$ (3.1)

######  Remark.

The feature map representation of the convolution kernel is given in Kriege et
al. (2014).

Haussler (1999) showed that this convolution is itself a positive definite
kernel on $\mathcal{X}$, simply called a _convolution kernel_. Many of the
graph kernels encountered in the literature are instances of convolution
kernels. For us, the sets $R^{-1}$ will contain decompositions of a graph into
subgraphs (Rupp and Schneider, 2010). One may wish to consult the original
paper of Haussler (1999) for more technical and general details.

The kernels $k_{i}$ in the context of graph kernels often are used to compare
vertices or edges of graphs, or their associated labels. The _Dirac_ (or
Kronecker) kernel defined as $k(x_{i},y_{i})=I(x_{i}=y_{i})$ (Rupp and
Schneider, 2010) is the most popular choice of base kernel in the literature,
due to its efficiency and interpretability. See Kriege et al. (2019) for a
description of its feature map.

###  3.2 Baseline Graph Kernels

Now that we have a general understanding of how graph similarity can be
quantified, we can discuss the first real attempts at constructing graph
kernels.

####  3.2.1 Simple Examples

The convolution kernel is perhaps the most natural extension for positive
definite kernel to discrete objects. In this short section, we present first-
examples of the convolution kernel on graphs, as found in Kriege et al.
(2020). They operate by comparing the labels between two graphs, while largely
ignoring structure.

######  Definition 3.2.1.

Let $G_{1}$ and $G_{2}$ be labeled graphs. The vertex label histogram (VH)
kernel is defined as

$K_{\textup{VH}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{v_{1}\in
V_{1}}\sum_{v_{2}\in V_{2}}\textup{{I}}(l(v_{1})=l(v_{2})).$ (3.2)

There is always a trade-off between _expressivity_ of a graph kernel (Kriege
et al., 2020), that is, the amount of structure and nuance it can take into
account, versus the computational complexity. This kernel is one of the
easiest to compute (having complexity $O(n^{2})$), but treats graphs as a _bag
of components_ , rather than a structured object. For example, two molecules
are said to be isomers if they share the same number of atoms of each type,
while possibly having different arrangements. This kernel would be unable to
distinguish the graphs of a family of isomers. This is but one example
illustrating the “low-resolution” of the VH kernel. While not a focus of this
review, the computational limits of graph kernels generally are discussed in
some papers such as Gärtner et al. (2003), Borgwardt and Kriegel (2005), and
Kriege et al. (2020).

######  Definition 3.2.2.

Let $G_{1}$ and $G_{2}$ be labeled graphs. The edge label histogram (EH)
kernel is defined as

$k_{\textup{EH}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{(v_{1},v_{1}^{\prime})\in
E_{1}}\sum_{(v_{2},v_{2}^{\prime})\in
E_{2}}\textup{{I}}[l(v_{1})=l(v_{2})]\cdot\textup{{I}}[l(v_{1},v_{1}^{\prime})=l(v_{2},v_{2}^{\prime})]\cdot\textup{{I}}[l(v_{1}^{\prime})=l(v_{2}^{\prime})].$
(3.3)

It is clear that the EH kernel has complexity $O(m^{2})$. The corresponding
feature maps for both of these kernels can be found in Sugiyama and Borgwardt
(2015). These kernels are not very accurate as is. The experiments in the
literature that include them typically do so for a lower-bound on performance.
However, as noted by Sugiyama and Borgwardt (2015), they are especially
important to measure against random walk kernels such as the geometric random
walk kernel, as a process known as halting can make those much more complex
kernels perform similarly to these “trivial” kernels. See section 3.2.2 for
more details. These simple kernels can be surprisingly effective with slight
modification; see section 3.4.1 for a brief discussion on composing linear
kernels with the RBF kernel.

####  3.2.2 Direct Product Kernel

The first graph kernels began appearing in the literature in 2003, with papers
by Gärtner et al. (2003) and Kashima et al. (2003) (Kriege et al., 2020).
These kernels measure graph similarity by comparing all possible pairings of
labeled walks between two graphs. This section is focused on the _direct
product kernel_ of Gärtner et al. (2003) and its variations. We postpone the
discussion of the marginalized kernel of Kashima et al. (2003) until section
3.3.1. Let us begin with the definition of the _direct product graph_ (see
Gärtner et al. (2003)). We use the notation from Vishwanathan et al. (2010).

######  Definition 3.2.3.

Let $G_{1},\,G_{2}$ be two labeled graphs. The direct product graph
$G_{\times}$ of $G_{1}$ and $G_{2}$ is a graph with vertex set

$V_{\times}\stackrel{{\scriptstyle\textup{def}}}{{=}}\\{(v_{1},v_{2})\in
V_{1}\times V_{2}:l(v_{1})=l(v_{2})\\}$

and edge set

$E_{\times}\stackrel{{\scriptstyle\textup{def}}}{{=}}\\{\left((v_{1},v_{2}),(v_{1}^{\prime},v_{2}^{\prime})\right)\in
V_{\times}\times V_{\times}:(v_{1},v_{2})\in
E_{1},\,(v_{1}^{\prime},v_{2}^{\prime})\in
E_{2},\,l((v_{1},v_{2}))=l(v_{1}^{\prime},v_{2}^{\prime})\\}.$

The graphs $G_{1}$ and $G_{2}$ are called the factor graphs of $G_{\times}$.

######  Remark.

Note that this formulation of the direct product graph requires vertices and
edges from the the factor graphs to have matching labels. Graph-theoretic
applications involving direct product graphs may omit this restriction from
the definition.

The direct product graph is useful as it allows for simultaneous labeled-walk
comparison on the factor graphs. It can be shown that the number of walks with
a particular label sequence on $G_{\times}$ is exactly the product of the
number of walks with that label sequence in each factor graph $G_{1}$ and
$G_{2}$ (Gärtner et al., 2003, Prop. 3).

For simplicity, we first present a special case of the direct product kernel
called the _vertex edge label histogram kernel_ (VEH), which was introduced by
Sugiyama and Borgwardt (2015). It combines both the vertex (definition 3.2.1)
and edge (definition 3.2.2) label histogram kernels; however, the direct
product adjacency matrix representation (Sugiyama and Borgwardt, 2015, Lemma
1) makes its function rule more compact and helps illustrate its connection
with the other direct product kernels. The feature map representation can be
found in the original paper.

######  Definition 3.2.4.

Let $G_{1}$ and $G_{2}$ be two labeled graphs, with $G_{\times}$ denoting
their direct product. The vertex edge label histogram kernel is defined as

$k_{\textup{VEH}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i=1}^{|V_{\times}|}\sum_{j=1}^{|V_{\times}|}\left[A_{\times}\right]_{ij}$
(3.4)

where $A_{\times}$ is the adjacency matrix of $G_{\times}$.

This kernel, counting the entries of $A_{\times}$, in effect counts the number
of edges in both graphs with matching edge labels and associated vertex
labels, or equivalently, counting random labeled walks of length one.

Next, we present the general direct product graph kernel of Gärtner et al.
(2003). It takes into account _every_ labeled walk of two given graphs.

######  Definition 3.2.5.

Let $G_{1}$ and $G_{2}$ be two graphs, and $G_{\times}$ their direct product
graph with adjacency matrix $A_{\times}$. Let $(\lambda_{n})_{n=0}^{\infty}$
be a sequence of weights consisting of nonnegative real numbers with the
property that the matrix power series

$\sum_{n=0}^{\infty}\lambda_{n}A_{\times}^{n}$ (3.5)

converges to a matrix $M$. Then the direct product kernel $k_{\times}$ is
defined as the entry-wise sum of the matrix $M$:

$k_{\times}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\,\sum_{i=1}^{|V_{\times}|}\sum_{j=1}^{|V_{\times}|}\left[M\right]_{ij}.$
(3.6)

######  Remark.

We avoid discussing the precise notion of convergence of sequences and series
of matrices, instead mentioning a couple of common choices for the weight
sequence that will ensure the matrix power series converges and has a familiar
closed-form. See Horn and Johnson (2012) for details regarding power series of
matrices.

The kernel uses the adjacency matrix representation of graphs, as it can
generate the number of walks between vertices in a graph. For a graph $G$ with
adjacency matrix $A$, the $ij$th element of $A^{k}$ contains the number of
walks of length $k$ from vertex $v_{i}$ to $v_{j}$ (see Duncan (2004), for
example). Thus for a direct product graph $G_{\times}$ and
$\textbf{v}_{i}=(v_{1},v_{2}),\,\textbf{v}_{j}=(v_{1}^{\prime},v_{2}^{\prime})\in
V_{\times}$, the $ij$th entry of $A_{\times}^{k}$ contains the number of walks
of length $k$ between $v_{1}$ and $v_{1}^{\prime}$, and $v_{2}$ and
$v_{2}^{\prime}$ that have the same label sequence.

Gärtner et al. (2003) present two examples of weight sequences with the
desired convergence property. The first corresponds to the coefficients of the
power series of the exponential function, which are given by
$\lambda_{n}=\frac{1}{n!}$. The advantage is that the resulting matrix power
series converges for any square matrix with real entries (Horn and Johnson,
2012). However, due to the nature of matrix multiplication, the exact limit
may be difficult to compute for general matrices. For diagonalizable matrices,
there is a known closed-form (Gärtner et al., 2003). Otherwise, one can sum
the first $N$ terms to estimate the kernel value. The second sequence—more
widely seen in the literature—chooses the weights so that the series converges
like a geometric series. More specifically, let $\gamma\in(0,\frac{1}{a})$,
where $a\geq\max_{v\in V_{\times}}\text{deg}(v)$, and set
$\lambda_{n}=\gamma^{n}$. Then it can be shown that the matrix power series
converges to

$\sum_{n=0}^{\infty}\left(\gamma A_{\times}\right)^{n}=(I-\gamma
A_{\times})^{-1}$ (3.7)

where $I$ is the identity matrix. The corresponding kernel is called the
_geometric random walk kernel_ (GRW) (Sugiyama and Borgwardt, 2015):

$k_{\textup{GRW}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\,\sum_{i=1}^{|V_{\times}|}\sum_{j=1}^{|V_{\times}|}\left[(I-\gamma
A_{\times})^{-1}\right]_{ij}.$ (3.8)

In this case, the computational complexity is dictated by the complexity of
the matrix inversion algorithm, which are approximately $O(n^{3})$ for an
$n\times n$ matrix. In the worst-case, the dimension of $G_{\times}$ is
$n^{2}\times n^{2}$ where $n$ is the number of nodes in both factor graphs,
hence the complexity of computing the direct product kernel is $O(n^{6})$
(Borgwardt et al., 2020). This is quite slow, but has been drastically
improved by Vishwanathan et al. (2010), who were able to achieve a complexity
of $O(dn^{3})$ for this kernel, where $d$ is the total number of labels, using
fixed-point iterations and conjugate gradient methods. This was used in an
experiment by Sugiyama and Borgwardt (2015) with labeled graphs, which we
discuss below. Furthermore, Kang et al. (2012) introduced methods to
approximate the kernel that allows for further reduction in complexity.

While Gärtner et al. (2003) provide no experimental results alongside their
direct product kernel, there are other instances in the literature where it
has been tested. For example, Sugiyama and Borgwardt (2015) studied the
problem of _halting_ found in random walk kernels, where the weights
associated with longer walks are so small that walks of length one dominate
the kernel value. This issue stems from the choice of coefficients that ensure
convergence, which approach zero rather quickly for the geometric and
exponential variants. To test halting, the authors introduce the _N-step
random walk kernel_ , which modifies the direct product kernel by truncating
the infinite series after a finite number of terms, meaning only finite walks
are incorporated.

######  Definition 3.2.6.

Let $G_{1}$ and $G_{2}$ be two labeled graphs, with $G_{\times}$ denoting
their direct product, and $N\in\mathbb{N}$. The N-step random walk kernel is
defined as

$k_{\times}^{N}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i=1}^{|V_{\times}|}\sum_{j=1}^{|V_{\times}|}\sum_{n=0}^{N}\lambda_{n}\left[A_{\times}^{n}\right]_{ij}$
(3.9)

for some real coefficients $\lambda_{n},\,n\in\\{0,\,\cdots,N\\}$.

With convergence no longer an issue, a more natural weighting scheme could be
applied, making the kernel more flexible. For their experiments, the weights
of this kernel were all set to 1 and the number of steps $N$ range from 1 to
10. Other kernels were used in the experiments as well. The VEH kernel was
used as a point of comparison against the GRW kernel to test their theoretical
results on halting, and the _Weifeiler-Lehman subtree_ (WLS) kernel (see
section 3.4.2) was used, as it is known to be one of the most consistently
accurate graph kernels. The data sets used are the standard ones found in the
literature: ENZYMES, NCI1, NCI109, MUTAG, and D&D, which all consist of
labeled nodes (Shervashidze et al., 2011). While the WLS kernel came out on
top in terms of accuracy of predictions, the $N$-step random walk kernel with
an optimal choice of steps (4-6 steps in this experiment) outperformed the GRW
kernel, sometimes by up to 5%. The effect of halting was amplified when the
parameter in the GRW kernel was set much smaller than its theoretical maximum.
The higher the parameter is set, the more likely convergence will become a
problem should a previously unseen pair of graphs create a direct product
graph with a larger maximum degree. If one wants to implement the GRW kernel,
Sugiyama and Borgwardt (2015) recommends including the simple vertex and edge
label histogram kernels (see section 3.2.1) in the initial testing as a
baseline to measure performance against.

####  3.2.3 Shortest-Path Graph Kernel

The next big advancement in the graph kernel literature was the _shortest-path
kernel_ , introduced by Borgwardt and Kriegel (2005). One of the main goals
for the authors was to improve runtime over the other kernels available at the
time that suffered from long computation times, such as the direct product
kernel discussed in the previous section.

First, we provide the definition of an edge walk from Borgwardt and Kriegel
(2005) that is very similar to that of a walk, but instead emphasises the
edges rather that the vertices. For a graph $G$, an _edge walk_ is a sequence
of edges $(e_{i})_{i=1}^{m}$ in $E$ with the property that for each $1<i\leq
m$, $e_{i-1}=(v_{i-1},v_{i-1}^{\prime})$ and $e_{i}=(v_{i},v_{i}^{\prime})$,
we have $v_{i-1}^{\prime}=v_{i}$. The shortest-path graph kernel uses the
adjacency matrix representation of graphs. However, it does so after what is
known as the _Floyd-transformation_ is applied to the graphs. This
transformation can be done using the _Floyd-Warshall_ algorithm (Floyd, 1962),
which takes as input the adjacency matrix of a graph and a matrix containing
the weights of edges (or distance between vertices). The output is a new
matrix, where the $ij$th element contains the length of the shortest path
between the $i$th and $j$th vertex. This new matrix is an adjacency matrix for
the transformed graph. If $n$ denotes the number of vertices in the original
graph, then the runtime of this algorithm is $O(n^{3})$. For chemical graphs,
one may use the adjacency matrix as the weight matrix. Pseudo-code for the
algorithm can be found in the paper by Borgwardt and Kriegel (2005). We are
now able to present the general form of the _shortest path kernel_ , which is
yet another instance of the convolution kernel framework.

######  Definition 3.2.7.

Let $G_{1}$ and $G_{2}$ be two graphs, with corresponding Floyd-transformed
graphs $F_{1}$ and $F_{2}$. The shortest-path graph kernel is defined as

$k_{\textup{SP}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{e_{1}\in
E_{F_{1}}}\sum_{e_{2}\in E_{F_{2}}}k_{walk}^{1}(e_{1},e_{2})$ (3.10)

where $k_{\textup{walk}}^{1}$ is a positive definite kernel on edge walks of
length 1.

######  Remark.

The feature map of $k_{SP}$ can be found in Shervashidze et al. (2011).

As the number of edges in the Floyd-transformed graph is $n^{2}$, where $n$ is
the number of vertices in the original graph, it follows from equation 3.10
that the shortest-path kernel has a computational complexity of $O(n^{4})$. An
immediate improvement over many standard forms of random walk kernels is the
built-in prevention of _tottering_ , which helps eliminate noise. Tottering
occurs in a walk when a vertex is immediately returned to after stepping away
from it, and tends to record uninformative features. Typically, the base
kernel $k_{walk}^{1}$ is the product of three kernels. If
$e_{1}=(v_{1},v_{1}^{\prime})$, $e_{2}=(v_{2},v_{2}^{\prime})$, then

$k_{\textup{walk}}^{1}(e_{1},e_{2})=k_{l}(l(v_{1}),l(v_{2}))\cdot
k_{d}(d(v_{1},v_{1}^{\prime}),d(v_{1},v_{1}^{\prime}))\cdot
k_{l}(l(v_{1}),l(v_{2}))$ (3.11)

where $k_{l}$ is a kernel on vertex labels, such as the Dirac kernel, and
$k_{d}$ is a kernel defined on the shortest-path distance. When analyzing
chemical graphs, a default distance function is used, where the distance
between two vertices is the number of edges in the shortest path between them.
Moreover, chemical graphs are usually connected, thus every pair of vertices
can be assigned a distance. In experiments, $k_{d}$ is also taken to be the
Dirac kernel.

For testing, Borgwardt and Kriegel (2005) used 10-fold one-versus-all SVM
classification on a bioinformatic data set containing 540 proteins. They
elected to use the Dirac kernel on vertex labels, but used the _Brownian
bridge kernel_ on both edge and vertex length (adding another kernel to the
product in equation 3.11). Its definition has been omitted in this review as
weighted graphs are rarely found in chemoinformatics. The test results showed
a clear improvement over random walk-kernels, both in terms of runtime and
classification accuracy. The shortest-path kernel has been used in many other
tests as well. Shervashidze et al. (2011) compared the shortest-walk kernel to
variants of the _Weisfeiler-Lehman_ graph kernels (section 3.4.2). Their
experiments showed that the shortest-path kernel outperformed many of the
kernels that came before it by a large margin, including variations of the
random walk kernel. Now, the Weisfeiler-Lehman techniques introduced can be
applied to many different kernels, as it functions as a label refinement
algorithm. Applying this to the shortest-path kernel, they found that it
resulted in the most accurate classifier, while remaining quick to compute on
most data sets.

###  3.3 Specialized Examples

The kernels in the last section are standard in the wider field of graph
kernels. In this section, we will look at graph kernels the were designed with
chemical graphs and chemoinformatic applications in mind. It is quite common
to find the kernels of section 3.2 in software packages for graph kernels.
More specialized kernels are less readily-available in general, however each
kernel in this section can be computed using the freely-available ChemCPP
package.

####  3.3.1 Marginalized Kernels

Many iterations of marginalized kernels have appeared throughout the
literature. First introduced by Kashima et al. (2003), these random walk
kernels offer more control over the direct product kernel by having user-set
probabilities for each labeled walk. We follow both Mahé et al. (2004) and
Mahé et al. (2005), which extend the original formulation by Kashima et al.
(2003) with the goal of achieving better results on chemical graphs.
Vishwanathan et al. (2010) provided a generalization of this kernel, along
with algorithms that improve computational complexity. Unfortunately, vertex
labels are not addressed, which is not particularly useful in traditional
chemoinformatic applications. The computational improvements have still been
applied to labeled graphs, as discussed in section 3.2.2.

Let $G$ be a labeled graph. The set $V^{*}=\cup_{n\in\mathbb{N}}V^{n}$ is used
to denote the set of finite-length sequences of vertices. If $\textbf{v}\in
V^{*}$, then $l(\textbf{v})$ is used to denote the associated label sequence,
which was defined in section 3.1.1. The general form of the marginalized graph
kernel of Kashima et al. (2003) is given in the following definition from Mahé
et al. (2004).

######  Definition 3.3.1.

Let $G_{1}$ and $G_{2}$ be two labeled graphs. The marginalized graph kernel
is defined as

$k_{\textup{M}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{h_{1}\in
V_{1}^{*}}\sum_{h_{2}\in
V_{2}^{*}}p_{1}(h_{1})p_{2}(h_{2})k_{L}(l(h_{1}),l(h_{2})).$ (3.12)

where $p_{1}$ and $p_{2}$ are probability distribution on $V_{1}^{*}$ and
$V_{2}^{*}$, respectively, and $k_{L}$ is a kernel defined on label sequences.

One could simply use the Dirac kernel as a means of comparing label sequences.
Mahé et al. (2004) provide a brief explanation on how to compute this kernel,
which makes use of product graphs and geometric series, as well as how to set
the probability distributions. The authors then present two ways of modifying
the kernel in equation 3.12. First is an iterative vertex label transformation
called the _Morgan index_ (Morgan, 1965) process that seeks to highlight paths
containing more relevant information, improving performance and cutting down
computational time. The process is as follows: to start, each vertex is
labeled with “1”. For each iteration $i$ in some predefined range and each
atom $A$ in the molecule, the label of $A$ is increased by the sum of the
label values of its direct neighbours from the previous iteration. To compute
this in practice, let $M_{0}$ be the vector of dimension equal to the number
of atoms, and populate $M_{0}$ with ones to represent the labels. On the $i$th
iteration, the vector defined by $M_{i}=AM_{i-1}=A^{i}M_{0}$, where $A$ is the
adjacency matrix of the graph, contains the desired labels. The Morgan index
has also been used with the _treelet kernel_ (Gaüzere et al., 2012).

The second way to modify the marginalized graph kernel is to change its
underlying probabilistic model. This is done to prevent the phenomenon of
_tottering_ , which occurs when a random walk revisits a vertex immediately
after stepping away from it. Walks containing totters are thought to introduce
unwanted noise when attempting to learn a model. This is theoretically
implemented using a 2nd-order Markov model; however, in practice a graph
transformation is applied to the two graphs which eliminates the possibility
of tottering walks, meaning that the kernel ends up using the same 1st-order
Markov process as the original formulation, except on a larger transformed
graph with an increased complexity (from $O(n^{2})$ to $O((n+m)^{2})$). Given
two graphs $G_{1}$ and $G_{2}$ and their corresponding transformed graphs
$G_{1}^{\prime}$ and $G_{2}^{\prime}$, respectively, the kernel can be written
as

$k_{\textup{MNT}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{h_{1}^{\prime}\in(V_{1}^{\prime})^{*}}\sum_{h_{2}^{\prime}\in(V_{2}^{\prime})^{*}}p_{1}^{\prime}(h_{1}^{\prime})p_{2}^{\prime}(h_{2}^{\prime})k_{L}(l^{\prime}(h_{1}^{\prime}),l^{\prime}(h_{2}^{\prime})).$
(3.13)

The paper provides experimental results for these kernels using SVMs, which
were implemented using GIST. The first data set used was the MUTAG data set.
The use of 1st and 2nd order Markov random walk models had a very small impact
on classification ability. On the high end, less than 2% in performance was
gained when the 2nd order model was used to eliminate tottered walks. Using
the Morgan index iteration on the graph vertices, a small increase in
classification was generally observed when between one and three iterations
were performed. Applying these iterations did lead to a dramatic decrease in
computational time. Two iterations alone reduced computational time by a
factor of approximately 240. The second data set used was the PTC data set.
The best results were obtained after between 8 and 10 Morgan index iterations,
peaking at $\sim\\!63$ ROC area (for context, we would expect about 50 ROC
area if we randomly classified data).

The ChemCPP toolbox contains functions to compute both the marginalized graph
kernels of Kashima et al. (2003) and the extended marginalized kernel of Mahé
et al. (2005).

####  3.3.2 Tree-Pattern Graph Kernel

Another popular substructure used to compare graphs is a special type of
subgraph known as a _tree_. There is a notion in which the structure of
organic molecules resemble that of a tree (Yamaguchi et al. (2003), Kriege et
al. (2020)), and based on empirical evidence, it appears that classifying
organic molecules based on trees embedded within the graph may lead to higher
accuracy. In the paper by Mahé and Vert (2009) which this section follows, a
graph kernel is proposed that can measure similarity based on common subtrees,
with a parameter controlling tree complexity. This builds off of the work of
Ramon and Gärtner (2003) with the hope that the kernel provides a way to
capture physicochemical properties of atoms. Mahé and Vert (2009) also provide
recursive algorithms to compute their kernels in practice, which we omit from
this presentation.

We begin by giving the appropriate definitions related to trees from Mahé and
Vert (2009). In this context, there are two important sets associated with
each vertex $v\in V$ in a directed graph $G$. The set of incoming neighbours
$\delta^{-}(v)=\\{u\in V:(u,v)\in E\\}$, and the set of outgoing neighbours
$\delta^{+}(v)=\\{u\in V:(v,u)\in E\\}$. The in-degree of a vertex $v$ is the
quantity $\displaystyle\left\lvert\delta^{-}(v)\right\rvert$ and the out-
degree is the quantity $\displaystyle\left\lvert\delta^{+}(v)\right\rvert$. A
_rooted tree_ $t$ is a directed, connected graph containing no cycles where
every node has in-degree 1, except for one node having in-degree 0, which is
known as the _root_ of the tree. For the remainder of this section, all trees
will be rooted trees. The nodes of $t$ with out-degree 0 are leaf nodes, while
the remaining nodes will be called interior nodes. The depth of a node is
defined as the length of the path from the root to it plus one. A tree where
each leaf node has the same depth $n$ is called a perfectly depth balanced
tree of order $n$, or simply a _balanced tree_.

To define the tree-pattern graph kernel, a way of formalizing the notion of a
tree being embedded in a graph is needed, as well as a way to count how many
times a tree is found within the graph. The following two definitions
accomplish this.

######  Definition 3.3.2.

Let $G=(V,\,E)$ be a graph and let $t=(V_{t},\,E_{t})$ be a tree with vertex
set $V_{t}=\\{\tau_{1},\,\cdots,\tau_{N}\\}$. We say that a $n$-tuple of
vertices $(v_{1},\,\cdots,\,v_{n})\in V^{n}$ is a tree-pattern of $G$ with
respect to $t$ if the following properties hold:

  1. i

) Matching vertex labels: For every $i\in\\{1,\,\cdots,n\\}$,
$l(v_{i})=l(\tau_{i})$.

  2. ii

) Edges correspondence: For every $(\tau_{i},\tau_{j})\in E_{t}$, we have
that $(v_{i},\,v_{j})\in E$ and $l((v_{i},v_{j}))=l(\tau_{i},\tau_{j}))$.
Moreover, if we also have that $(\tau_{i},\tau_{k})\in E_{t}$, then $j\neq k$
if and only if $v_{j}\neq v_{k}$.

The set of all such $n$-tuples will be denoted by $\textup{Pattern}(t,G)$.

An important subtlety of this definition is that a vertex in the graph may be
used multiple times in the tree pattern. In effect, this will allow for two
edges with opposite orientation connecting two vertices (i.e. the edges in
chemical graphs, which represent covalent bonds) to be part of the tree-
pattern.

######  Definition 3.3.3.

The tree-pattern counting function, denoted by
$\varphi_{t}(G)=|\textup{Pattern}(t,G)|$, counts the number of occurrences of
a tree pattern $t$ in $G$.

We may now present the general form of the tree-pattern graph kernel, which is
another instance of the convolution kernel framework (Rupp and Schneider
(2010), Vishwanathan et al. (2010)).

######  Definition 3.3.4.

Let $T$ be a set of trees, $w(t)$ a nonnegative weight function defined on
$T$, and $G_{1}$, $G_{2}$ two labeled graphs. The tree-pattern graph kernel is
defined as

$k_{\textup{TPK}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{t\in
T}w(t)\varphi_{t}(G_{1})\varphi_{t}(G_{2})$ (3.14)

where $\varphi_{t}$ is the tree-pattern counting function.

The two specific forms of tree-pattern graph kernels found in this paper both
take $T$ to be the set of balanced trees of order $h$, which they denote by
$B_{h}$. Where the two kernels differ is how they assign weight to each tree.
Let $\lambda\geq 0$ be a nonnegative hyperparameter. The _size-based balanced
tree-pattern kernel_ $k_{\text{size}}^{h}$ is equation 3.14 with
$w(t)=\lambda^{|t|-h}$, where $|t|$ denotes the number of nodes in the tree:

$k_{\text{size}}^{h}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{t\in
B_{h}}\lambda^{|t|-h}\varphi_{t}(G_{1})\varphi_{t}(G_{2}).$ (3.15)

This kernel may be generalized by enlarging the set of trees under
consideration, such as using the set of trees of depth up to and including
$h$, denoted by $T_{h}$. The _branching-based balanced tree-pattern kernel_
$k_{\text{Branch}}^{h}$ is equation 3.14 with
$w(t)=\lambda^{\textup{branch}(t)}$, where $\textup{branch}(t)$ equals the
number of leaf nodes in $t$ plus one:

$k_{\text{branch}}^{h}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{t\in
B_{h}}\lambda^{\textup{branch}(t)}\varphi_{t}(G_{1})\varphi_{t}(G_{2}).$
(3.16)

A goal of the authors was to control the complexity of tree patterns
incorporated into the kernels. Informally, higher complexity of a tree means
it has a high number of leaf nodes, or a high number of internal nodes. Thus
complexity refers to how “nonlinear” a tree is. When $\lambda>1$, more weight
is placed of complex tree patterns, and when $\lambda<1$, more weight is put
on simpler tree patterns (Rupp and Schneider, 2010). In fact, the authors
remark that as $\lambda$ approaches 0, both kernels approach the value of the
direct product kernel (section 3.2.2). The computational complexity of both of
these tree-pattern kernels is $O(hn^{2}d^{2d})$, where $d$ is the upper bound
on the number of out-degrees of the vertices and $n$ is the number of vertices
in one graph. For chemical compounds in chemoinformatics, Mahé and Vert (2009)
claim that it is almost always the case that $d$ is no greater than 4.
Moreover, due to the computational algorithms introduced, the extension of the
size-based kernel to the set $T_{h}$ carries no extra computational cost.

A phenomenon known as _tottering_ , first noticed with walk-based kernels
(Mahé et al., 2005), can also impact the effectiveness of these tree kernels.
Tottering happens when a tree-pattern contains a vertex from the graph as both
the parent and child of another vertex. Tree-patterns of this type end up
adding noise that tends to obscure the important features, as they vastly
outnumber non-tottering tree patterns as the depth of the tree considered
increases. The way the authors ultimately handle this is to apply a
transformation to the graph that produces a new graph without any tottering
walks (the same transformation mentioned in section 3.3.1). This _non-
tottering kernel_ takes the form

$k_{\text{NT}}(G_{1},G_{2})=\sum_{t\in
T}w(t)\varphi_{t}^{NT}(G_{1}^{\prime})\varphi_{t}^{NT}(G_{2}^{\prime})$ (3.17)

where $\varphi_{t}^{NT}$ is the no-tottering tree-pattern counting function
(see the original paper), and $G_{1}^{\prime},\,G_{2}^{\prime}$ are the
transformed graphs of $G_{1},\,G_{2}$, respectively. This ends up scaling the
computational complexity by

$\frac{(n+m)^{2}}{n^{2}}$ (3.18)

and hence the resulting complexity is $O(h((n+m)d^{d})^{2})$. The average
observed value of equation 3.18 across the experimental data sets used by the
authors was found to be approximately equal to 9.

Mahé and Vert (2009) provide experimental results on their proposed kernels
and their extensions. Classification was performed with SVMs using LIBSVM, as
well as the PyML framework. The graph kernels were computed using ChemCPP. The
first series of experiments were performed on two public data sets, the first
being the common MUTAG data set and the second being the now unavailable NCI
GI50 data set of potential tumor-suppressing compounds. Optimal values of
$\lambda$ appeared to decrease for larger $h$, and the author’s reasoning was
that the number of tree patterns increases exponentially as $h$ increases and
a smaller $\lambda$ allows for less individual influence of each tree. When
$h$ got too large ($h\geq 8$), they noticed convergence issues if $\lambda$
wasn’t small enough. However, when testing the no-tottering extensions, this
problem went away and a small accuracy improvement was seen. It appears that
the optimal value of $h$ is highly dependent on the data set and what
differentiate the molecules. Overall, the kernels performed better than the
walk-based kernels. In Sawada et al. (2014), these kernels are tested against
a few other common kernels. They are mostly comparable to Tanimoto kernels,
but ultimately lose out against them in nearly every test when it comes to
accuracy.

####  3.3.3 Tanimoto Kernels

In the paper by Ralaivola et al. (2005), three closely related kernels are
developed for direct application to problems in chemoinformatics. These
kernels represent molecules using a technique known as _molecular
fingerprinting_ , which encodes certain paths111The authors use a weaker
definition of path than what we presented in section 3.1.1, only requiring
paths to have distinct edges, and not necessarily distinct vertices. We will
keep consistent with this terminology only in this section. found within the
graph in feature vectors. The thinking is that the molecular fingerprinting
representation is more apt for classifying the molecules found in organic
chemistry. The way in which molecules are compared is based on the _Tanimoto
similarity measure_ (Fligner et al., 2002), which is a well-known measure of
chemical distance.

Viewing molecules as labeled graphs, the molecular fingerprinting technique
assigns values to paths emanating from each vertex, and stores them in a
feature vector. The paths are computed via a depth-first search algorithm, for
which not all implementations are created equal (see Ralaivola et al. (2005)
for descriptions of common implementations). A decision needs to be made on
whether cycles are allowed, the maximum depth, and if the same edge can be
visited in two paths with the same starting node after the first point of
divergence. For an example of the latter, if _a_ , _b_ , _c_ , and _d_ label a
square-shaped graph in the clockwise sense, and if the path $(a,b,c,d)$ is
traversed first, then $(a,d,c)$ would not be traversed. All of these
modifications have an effect on both complexity and effectiveness of
fingerprinting. For example, if a molecule contains $n$ atoms and $m$ bonds,
then finding all paths up to depth $d$ with the edge-divergence condition is
$O(nd)$. Without the edge-divergence condition, the complexity becomes
$O(n\alpha^{d})$, where $\alpha$ is the average number of atoms that are
neighbours to any given atom in the graph. The authors note that the value of
$\alpha$ is typically low for organic molecules.

Let $G$ denote the graph of a molecule, and let $\mathcal{P}(d)$ denote the
set of all atom-bond labeled paths of length $d$ in $G$. There are two
approaches to fingerprinting that are considered by Ralaivola et al. (2005).
The first is called the _binary feature map of $G$ for a depth $d$_, and is a
vector-valued function given by

$\varphi_{d}(G)=(\textbf{I}(p\subset G))_{p\in\mathcal{P}(d)}$ (3.19)

where $\textbf{I}(p\subset G)$ is the indicator function that outputs 1 if at
least one depth-first search on $G$ of depth at most $d$ produces the path
$p$, and 0 otherwise. Notice that this disregards multiple separate instances
of $p$ appearing in $G$. The second approach is called the _counting feature
map_ , denoted by $\phi_{d}$, and is given by

$\phi_{d}(G)=(\\#\\{p\subset G\\})_{p\in\mathcal{P}(d)}.$ (3.20)

This feature map does record the number of times a path $p$ appears in the
graph $G$. The notation $\\#\\{p\subset G\\}$ is adapted from Klambauer et al.
(2015). The third approach to fingerprinting is a fixed-length binary feature
map of depth $d$ and length $r$, denoted by $\overline{\varphi}_{d,r}(G)$.
Every $p\in\mathcal{P}(d)$ is mapped to $b$ indices between 1 and r by some
function $f$. The corresponding indices in the vector
$\overline{\varphi}_{d,r}(G)$ are then set to 1 (and remain 1 even if it is
mapped to again). Standard implementations for $f$ fix $b$ to be either 1 or
4, and compute a hash value for each path $p$ as a seed of a random number
generator that chooses $b$ random integers, and reduces them mod $r$. Note
that one may wish to use very large length feature vectors, having a bit
position for each possible path, to eliminate the possibility of hash
collisions that can cause loss of information.

The kernels introduced in this paper are built out of the simpler dot product
kernels on the molecular fingerprint vectors. For example,

$k_{\varphi_{d}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\langle\varphi_{d}(G_{1}),\,\varphi_{d}(G_{2})\rangle=\sum_{p\in\mathcal{P}(d)}\textbf{I}(p\subset
G_{1})\textbf{I}(p\subset G_{2}).$ (3.21)

When using the binary feature map of $G$ of size $r$, this dot product kernel
is denoted by $k_{d,r}$. The dot product kernel using the counting feature map
is not explicitly used. We are now ready to present the main kernel of this
section.

######  Definition 3.3.5.

Let $G_{1},\,G_{2}$ be two graphs, and $d\in\mathbb{N}$ denoting the maximum
search depth being considered. The Tanimoto kernel is defined by

$k_{d}^{t}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{k_{\varphi_{d}}(G_{1},G_{2})}{k_{\varphi_{d}}(G_{1},G_{1})+k_{\varphi_{d}}(G_{2},G_{2})-k_{\varphi_{d}}(G_{1},G_{2})}.$
(3.22)

The Tanimoto kernel essentially computes the ratio between the number of
features extracted from both $G_{1}$ and $G_{2}$, and the total number of
features extracted from $G_{1}$ and $G_{2}$ (without double counting). The
codomain of the kernel is $[0,1]$.

######  Definition 3.3.6.

Let $G_{1},\,G_{2}$ be two graphs, and $d\in\mathbb{N}$. The MinMax kernel is
defined as

$k_{d}^{m}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{\sum_{p\in\mathcal{P}(d)}\min(\\#\\{p\subset
G_{1}\\},\,\\#\\{p\subset
G_{2}\\})}{\sum_{p\in\mathcal{P}(d)}\max(\\#\\{p\subset
G_{1}\\},\,\\#\\{p\subset G_{2}\\})}$ (3.23)

where $\\#\\{p\subset G\\}$ denotes the number of occurrences of the path p in
G.

The biggest difference with the MinMax kernel is that a path appearing
multiple times in a graph is factored into the calculation, where as the
Tanimoto kernel only checks for one instance of a path. In the testing portion
of the paper, this is the reason cited for why the MinMax kernel performs
better on molecules with different sizes.

######  Definition 3.3.7.

Let $G_{1},\,G_{2}$ be two graphs, $r\in\mathbb{N}$ representing feature
vector length, and $d\in\mathbb{N}$ denoting the maximum search depth
considered. Let $c\in(-1,2)$ be a user-defined parameter. The Hybrid kernel is
defined as

$k_{d,r}^{h}(G_{1},\,G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{1}{3}[(2-c)(k_{d,r}^{t}(G_{1},G_{2}))+(1+c)(1-k_{d,r}^{t}(G_{1},G_{2}))].$
(3.24)

The notation $k_{d,r}^{t}$ is used to denote the Tanimoto kernel using the
same base kernel as $k_{d,r}$. Taking one minus this altered Tanimoto kernel
is equivalent to a logical negation. According to the authors of the paper,
$c$ is often taken to be the average density of the feature vectors. This
kernel is designed to be a convex combinations of two altered Tanimoto
kernels.

Next, we discuss the complexity of these kernels following the analysis of
Ralaivola et al. (2005). Using a suffix tree (see Ukkonen (1995) for a
construction algorithm) to store each path emanating from a particular node
and not distinguishing orientation of paths, the complexity of each of the
three kernels acting on two graphs $G_{1}$ and $G_{2}$ is $O(dnm)$, where $n$
and $m$ again represent the number of atoms and bonds of a graph, and $d$ is
the maximum path-length in the search. Note that this can increase depending
on the particulars of the depth-first search algorithm used.

Ralaivola et al. (2005) tested their kernels using the voted perceptron
classifier (Freund and Schapire, 1999) as a model on three data sets to test
the prediction power of mutagenicity, toxicity, and anti-cancer activity. The
first data set used was the MUTAG data set, the second was the PTC data set,
and the third was the NCI GI50 data set. We were only able to find actives
links to retrieve MUTAG. For kernel parameters, the depth was set to 10, $b$
was set to 1, and feature vectors of size of 512 and 1024 (when restricted)
were both used. The tests showed that the Tanimoto and MinMax kernels in
particular were quite effective in classification, with the MinMax kernel
achieving the highest accuracy on the MUTAG data set at the time the paper was
published. Both kernel were above 70% accuracy on the NCI data set, with the
MinMax kernel performing slightly better.

The Tanimoto kernel was also tested against the tree-pattern, marginalized,
and extended marginalized graph kernels in the paper by Sawada et al. (2014).
A section of this paper is dedicated to benchmarking the predictive power of
these four popular kernels on a data set of over 100,000 unique drug-target
interactions using a chemogenomics approach. A pairwise kernel regression
model was employed. The Tanimoto kernel performed best in many of the tests,
and in the others remained very competitive. Detailed descriptions of the
particular classification tasks can be found in the original paper.

###  3.4 Boosting Kernel Performance

Many methods have been proposed to improve both classification accuracy and
runtime performance of preexisting graph kernels. This can be done by
extending existing kernels, for example the extension of the marginalized
graph kernel discussed in section 3.3.1, or by introducing a general scheme
that can be applied to classes of graph kernels. Our focus in this section
will be on the latter.

####  3.4.1 Graph RBF Kernel

While the RBF kernel (section 2.2.2) is commonly viewed as a standalone kernel
with its own feature map and associated RKHS, one could also view it as a
function composed with preexisting kernels. In its original formulation, it
contains the squared Euclidean distance, which can be written as the sum of
dot products, or _linear kernels_. The use of feature maps is how the graph
kernel variant of the RBF kernel is constructed. First, we present the _kernel
metric_ defined in Steinwart and Christmann (2008), which will replace the
Euclidean metric for the general RBF kernel.

######  Definition 3.4.1.

Let $k$ be the kernel of an RKHS on a set $U$ with feature map $\varphi$. The
kernel metric on $U$ is defined as

$d_{k}(x,y)\stackrel{{\scriptstyle\textup{def}}}{{=}}\|\varphi(x)-\varphi(y)\|=\sqrt{k(x,x)-2k(x,y)+k(y,y)}$
(3.25)

for every $x,\,y\in U$.

######  Remark.

The kernel metric is a pseudo-metric on $U$. For it to be a metric, we would
also need for $\varphi$ to be injective.

Now, if $k$ is some fixed graph kernel, then the _graph RBF kernel composed
with k_ is defined as

$k_{\text{GRBF}}(G_{1},G_{2})\stackrel{{\scriptstyle\textup{def}}}{{=}}\exp{\left(-\frac{d_{k}(G_{1},G_{2})^{2}}{2\sigma^{2}}\right)}$
(3.26)

where $G_{1}$ and $G_{2}$ are two graphs, and $\sigma$ is a hyperparameter.
The experiments conducted by Kriege et al. (2020) used cross-validation to
choose $\sigma$ from the set $\\{2^{-7},2^{-6},\cdots,2^{7}\\}$. Overall, they
found that the combination of the graph RBF kernel with other graph kernels
led to a minor accuracy boost, typically a few percentage points. However,
certain kernels received a large performance boost, such as the edge label
histogram kernel (section 3.2.1), which became on-par with more complex
kernels. The general trade-off is that optimizing $\sigma$ can be expensive,
especially as the number of graphs in the data set grows. They recommend that
the RBF kernel be used in conjunction with the vertex and edge label histogram
kernels in any situation, and to avoid it with kernels such as the Weisfeiler-
Lehman family of kernels who already have a hyperparameter to tune, and see
negligible gains. Sugiyama and Borgwardt (2015) composed the graph RBF kernel
with the vertex edge label histogram kernel (definition 3.2.4). They found a
significant increase in classification accuracy with this variant of the
kernel, with up to a 15% increase compared to the original.

####  3.4.2 Weisfeiler-Lehman Kernels

Graphs —like most structures in mathematics —have a very important notion of
_equivalence_ known as _isomorphism_. Informally, two graphs $G_{1}$ and
$G_{2}$ are said to be isomorphic if they have the same _structure_. That is,
if there is a bijection $\varphi:V_{1}\mapsto V_{2}$ with the property that
$(v,v^{\prime})$ is an edge in $G_{1}$ if and only if
$(\varphi(v),\varphi(v^{\prime}))$ is an edge in $G_{2}$, and the
corresponding vertices and edges with respect to $\varphi$ have the same
labels (see Kriege et al. (2020), for example).

Weisfeiler-Lehman kernels (Shervashidze and Borgwardt (2009), Shervashidze et
al. (2011)) are based on a procedure known as the Weisfeiler-Lehman (WL) test
of isomorphism (Weisfeiler and Leman, 1968). This is an algorithm that can be
used to show that two graphs are _not_ isomorphic. No conclusion can be made
if the algorithm terminates naturally. The idea is to iteratively assign
labels to vertices of a graph depending on the previous labels of direct
neighbours. These labels can then be compared between graphs; if they are the
same, repeat the process, if not, the graphs are not isomorphic. The algorithm
naturally terminates after $n$ iterations, where $n$ is the number of vertices
in the graphs. The complexity of this algorithm is $O(hm)$, where $h$ is the
number of iterations, which if leveraged by a graph kernel can offer much
better computational-time scaling with respect to graph size than other
popular kernels. Note that this complexity is dependent on the type of sorting
algorithm implemented within the WL algorithm. The authors achieve this
complexity through the implementation of a _counting sort_. The full WL
algorithm is presented in the paper by Shervashidze et al. (2011).

The general WL kernel is constructed as follows. For each iteration $i$ of the
WL algorithm on a graph $G$, a new label $l_{i}$ for $G$ is created. Let
$G_{i}=(V,E,l_{i})$ be the _WL graph_ at height $i$. When $i=0$, we simply set
$G_{0}=G$ and $l_{0}=l$. If $k$ is a graph kernel, then the _WL kernel with h
iterations and base kernel k_ is defined as

$K_{\textup{WL}}^{h}(G,G^{\prime})\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i=0}^{h}k(G_{i},G_{i}^{\prime}).$
(3.27)

This framework allows for many new kernels to leverage the label refinement
produced by the WL algorithm. The first particular instance presented was
introduced by Shervashidze and Borgwardt (2009), and is called the Weisfeiler-
Lehman _subtree kernel_. It was later shown by Shervashidze et al. (2011) to
be a special case of equation 3.27. The latter formulation is what we will
present. For context, the goal outlined in the original paper by Shervashidze
and Borgwardt (2009) was to develop a fast kernel on labeled graphs utilizing
subtree structure. As the authors note, methods for speeding up computational
time for walk-based and subgraph-based kernels were known, but the same could
not be said for graph kernels using trees. It is known that graph kernels
scale poorly with respect to graph size, which is a limiting factor on any
method. For two graphs $G_{1}$ and $G_{2}$, the _Weisfeiler-Lehman subtree
kernel with h iterations_ , denoted by $k_{WLS}^{h}$, is given by the general
WL kernel in equation 3.27 with the vertex label kernel as its base (section
3.2.1). Hence, it can be written as

$k_{\textup{WLS}}^{h}(G_{1},G_{2})=\sum_{i=0}^{h}\sum_{v_{1}\in
V_{1}}\sum_{v_{2}\in V_{2}}\textbf{I}(l_{i}(v_{1})=l_{i}(v_{2})).$ (3.28)

Note that the value of this kernel is a byproduct of running the algorithm on
two graphs. This implies that complexity of this kernel is $O(hm)$
(Shervashidze et al., 2011, Theorem 5). However, the author provide an
algorithm based on explicit feature map computation that allows for the Gram
matrix of $k_{WLS}^{h}$ on $N$ graphs to be computed in time $O(Nhm+N^{2}hn)$,
where $m$ and $n$ are the maximum number of edges and vertices (respectively)
over the $N$ graphs. See section 3.4.4 for a short discussion regarding
efficiency gains from explicit feature map computation.

The next formulation in Shervashidze et al. (2011) is the _WL edge kernel_ ,
which is again given by the general framework in equation 3.27, but now with
the edge label kernel as its base (section 3.2.1). As the base kernel has
complexity $O(m^{2})$, an upper bound on the complexity of the WL edge kernel
is $O(hm^{2})$. The final WL kernel discussed employs the shortest-path kernel
introduced by Borgwardt and Kriegel (2005) for the base kernel (see section
3.2.3). We briefly review concepts from that section. The first step in
computing the WL-SP kernel is to determine the shortest-paths between any two
vertices of a graph. In Borgwardt and Kriegel (2005), this is done by applying
the Floyd-transformation to a graph, which outputs a matrix whose $ij$th entry
is the distance (measured in terms of number of edges for chemical graphs)
between vertices $v_{i}$ and $v_{j}$. The shortest-path kernel $k_{sp}$ then
counts how many shortest-length paths in each graph share the same path length
and the same labels for the start and end vertices. The complexity for
computing this kernel has an upper bound of $O(hn^{4}+hm)$.

Shervashidze et al. (2011) conducted two sets of experiments. The first is an
empirical assessment of the runtime performance of the WL subtree kernel on
toy data sets using both the pairwise kernel computation scheme, as well as
the “global” $N$-graph implementation. The results showed that “global” option
was far quicker, and appears to scale well with graph size. The second
compares the predicting power of three iterations of WL kernels against
popular graph kernel of the time using standard benchmarking data sets. Many
kernels were used in this experiment, including the generalized random-walk
kernel from Vishwanathan et al. (2010), the graphlet kernel of Shervashidze et
al. (2009). LIBSVM was implemented for SVM classification, and experiments
were repeated 10 times. The height $h$ for the WL subtree kernel was chosen
via 10-fold cross-validation in the range $\\{0,1,\,\cdots,10\\}$. Note that
for WL edge and WL shortest-path kernels, values of 2 and 3 for $h$ were
chosen by the model almost exclusively. The results showed that the WL subtree
kernel was able to handle large graphs with thousands of vertices, and was
very competitive on small data sets too. For the D&D data set containing 1178
protein graphs, the WL kernel took 11 minutes to compute, where as the other
WL-based kernels took anywhere between 23 hours and over a year. In terms of
classification performance, the WL subtree kernel did particularly well on the
NCI1 and NCI109 data sets. The WL shortest-path kernel performed far better
than the other kernels on the smaller ENZYMES data set. Overall, the WL
kernels were either outperforming other kernels, or about on-par with them, in
terms of both classification accuracy and CPU runtime.

Kriege et al. (2016) used the WL kernel framework to develop the WL optimal
assignment kernel (section 3.4.3), which according to the tests done in Kriege
et al. (2020), performs better than the standard WL subtree kernel. Morris et
al. (2021) discusses the applications of the WL algorithm in machine learning,
including the WL kernel and beyond.

####  3.4.3 Optimal Assignment Kernels

So far, all of the graph kernels we have discussed have been instances of the
convolution kernel framework. In this section, we briefly discuss a newer
class of kernels called _optimal assignment kernels_ that leverage a different
framework. While convolution kernels sum over all possible substructures of a
certain type, optimal assignment kernels work by selecting the particular
pairings of substructures that maximized the function value. Many variations
of the optimal assignment kernel have appeared in the literature. After being
introduced by Fröhlich et al. (2005), it was pointed out by Rupp et al. (2007)
that the kernels need not be positive definite. The loss of this property