Matérn Gaussian Processes on Graphs

Viacheslav Borovitskiy$*1,5$ Peter Mostowsky$1$ Iskander Azangulov$*1$ Marc
Peter Deisenroth$3$ Alexander Terenin$*2$ Nicolas Durrande$4$

$1$St. Petersburg State University $2$Imperial College London $3$Centre for
Artificial Intelligence, University College London $4$Secondmind $5$St.
Petersburg Department of Steklov Mathematical Institute of Russian Academy of
Sciences

###### Abstract

Gaussian processes are a versatile framework for learning unknown functions in
a manner that permits one to utilize prior information about their properties.
Although many different Gaussian process models are readily available when the
input space is Euclidean, the choice is much more limited for Gaussian
processes whose input space is an undirected graph. In this work, we leverage
the stochastic partial differential equation characterization of Matérn
Gaussian processes—a widely-used model class in the Euclidean setting—to study
their analog for undirected graphs. We show that the resulting Gaussian
processes inherit various attractive properties of their Euclidean and
Riemannian analogs and provide techniques that allow them to be trained using
standard methods, such as inducing points. This enables graph Matérn Gaussian
processes to be employed in mini-batch and non-conjugate settings, thereby
making them more accessible to practitioners and easier to deploy within
larger learning frameworks.

##  1 Introduction

Gaussian process (GP) models have become ubiquitous in machine learning, and
have been shown to be a data efficient approach in a wide variety of
applications [29]. Key elements behind the success of GP models include their
ability to assess and propagate uncertainty, as well as encode different kinds
of prior information about the function they seek to approximate. For example,
by choosing different covariance kernels, one can encode different degrees of
differentiability, or specific patterns, such as periodicity and symmetry.
Although the input and output spaces of GP models are typically subsets of
$\mathbb{R}$ or $\mathbb{R}^{d}$, this is by no means a restriction of the GP
framework, and it is possible to define models for other types of input or
output spaces [23, 24, 6].

$*$Equal contribution. Code available at: http://github.com/spbu-math-
cs/Graph-Gaussian-Processes. Correspondence to:
<EMAIL_ADDRESS>
In many applications, such as predicting street congestion within a road
network, using kernels based on the Euclidean distance between two locations
in the city does not make much sense. In particular, locations that are
spatially close may have different traffic patterns, for instance, if two
nearby roads are disconnected or if traffic is present only in one direction
of travel. Here, it is much more natural for the model to account directly for
a distance based on the graph structure. In this work, we study GPs whose
inputs or outputs are indexed by the vertices of an undirected graph, where
each edge between adjacent nodes is assigned a positive weight.

Gaussian Markov random fields (GMRF) [30] provide a canonical framework for
such settings. A GMRF builds a graph GP by introducing a Markov structure on
the graph’s vertices, and results in models and algorithms that are
computationally efficient. These constructions are well-defined and effective,
but they require Markovian assumptions, which limit model flexibility.
Although it may be tempting to replace the Euclidean distance, which can
typically be found in the expression of stationary kernels, by the graph
distance, this typically does not result in a well-defined covariance kernel
[13].111Approaches of this type can still be employed, provided one introduces
additional restrictions on the class of graphs and/or on allowable ranges for
parameters such as length scale. See [3] for an example.

As a consequence, working with GPs on graphs requires one to define bespoke
kernels, which, for graphs with finite sets of nodes, can be viewed as
parameterized structured covariance matrices that encode dependence between
vertices. A few covariance structures dedicated to graph GPs have been
explored in the literature, such as the diffusion kernel or random walk
kernels [20, 41], but the available choices are limited compared to typical
Euclidean input spaces and this results in impaired modeling abilities.

In this work, we study graph analogs of kernels from the Matérn family, which
are among the most commonly used kernels for Euclidean input spaces [29, 35].
These can be used as GP input covariances or GP output cross-covariances. Our
approach is related to [44, 23, 6], where GPs with Matérn kernels are defined
on Euclidean and Riemannian manifolds via their stochastic partial
differential equation (SPDE) representation. For the graph Matérn GPs with
integer smoothness parameters we obtain sparse precision matrices that can be
exploited for improving the computational speed. For example, they can benefit
from the well-established GMRF framework [30] or from recent advances in non-
conjugate GP inference on graphs [10]. As an alternative, we present a Fourier
feature approach to building Matérn kernels on graphs with its own set of
advantages, such as hyperparameter optimization without incurring the cost of
computing a matrix inverse at each optimization step.

We also discuss important properties of the graph Matérn kernels, such as
their convergence to the Euclidean and Riemannian Matérn kernels when the
graph becomes more and more dense. In particular, this allows graph Matérn
kernels to be used as a substitute for manifold Matérn kernels [6] when the
manifold in question is obtained as intermediate output of an upstream
manifold learning algorithm in the form of a graph.

##  2 Gaussian processes

Let $X$ be a set. A random function $f:X\rightarrow\mathbb{R}$ is a Gaussian
process $f\sim\operatorname{GP}(\mu,k)$ with mean function $\mu(\cdot)$ and
kernel $k(\cdot,\cdot)$ if, for any finite set of points $\boldsymbol{x}\in
X^{n}$, the random vector $f(\boldsymbol{x})$ is multivariate Gaussian with
mean vector $\boldsymbol{\mu}=\mu(\boldsymbol{x})$ and covariance matrix
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}=k(\boldsymbol{x},\boldsymbol{x})$.
Without loss of generality, we assume that the prior mean $\mu$ is zero.

For a given set of training data $(x_{i},y_{i})$, we define the model
$y_{i}=f(x_{i})+\varepsilon_{i}$ where $f\sim\operatorname{GP}(0,k)$ and
$\varepsilon_{i}\sim\operatorname{N}(0,\sigma^{2})$. The posterior of $f$
given the observations is another GP. Its conditional mean and covariance are

$\displaystyle\mu_{\mid\boldsymbol{y}}(\cdot)$
$\displaystyle=\mathbf{K}_{\boldsymbol{\cdot}\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I})^{-1}\boldsymbol{y}$
(1) $\displaystyle k_{\mid\boldsymbol{y}}(\cdot,\cdot^{\prime})$
$\displaystyle=k(\cdot,\cdot^{\prime})-\mathbf{K}_{\boldsymbol{\cdot}\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{K}_{\boldsymbol{x}\boldsymbol{\cdot}^{\prime}}$
(2)

which uniquely characterize the posterior distribution [29]. Following [45],
posterior sample paths can be written as

$f(\cdot)\mid\boldsymbol{y}=f(\cdot)\\!+\\!\mathbf{K}_{\boldsymbol{\cdot}\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}\\!+\\!\sigma^{2}\mathbf{I})^{-1}(\boldsymbol{y}\\!-\\!f({\boldsymbol{x}})\\!-\\!\boldsymbol{\varepsilon})$
(3)

where $(\cdot)$ is an arbitrary set of locations.

###  2.1 The Matérn kernel

When $X=\mathbb{R}^{d}$ and $\tau=x-x^{\prime}$, Matérn kernels are defined as

$\displaystyle
k_{\nu}(x,x^{\prime})=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\bigg{(}$
$\displaystyle\sqrt{2\nu}\frac{\left\lVert\tau\right\rVert}{\kappa}\bigg{)}^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{\left\lVert\tau\right\rVert}{\kappa}\right)$
(4)

where $K_{\nu}$ is the modified Bessel function of the second kind [15]. The
parameters $\sigma^{2}$, $\kappa$ and $\nu$ are positive scalars that have a
natural interpretation: $\sigma^{2}$ is the variance of the GP, the length-
scale $\kappa$ controls how distances are measured in the input space, and
$\nu$ determines mean-square differentiability of the GP [29]. As
$\nu\rightarrow\infty$, the Matérn kernel converges to the widely-used squared
exponential kernel

$k_{\infty}(x,x^{\prime})=\sigma^{2}\exp\bigg{(}-\frac{\left\lVert\tau\right\rVert^{2}}{2\kappa^{2}}\bigg{)}.$
(5)

Euclidean Matérn kernels are known to possess favorable asymptotic properties
in the large-data regime—see [34, 19] for details.

In our setting, an important property of Matérn kernels is their connection to
stochastic partial differential equations (SPDEs). [44] has shown that Matérn
GPs on $X=\mathbb{R}^{d}$ satisfies the SPDE

$\left(\frac{2\nu}{\kappa^{2}}-\Delta\right)^{\frac{\nu}{2}+\frac{d}{4}}f=\mathcal{W}$
(6)

for $\nu<\infty$, where $\Delta$ is the Laplacian and $\mathcal{W}$ is
Gaussian white noise [22] re-normalized by a certain constant—see [23] or [6]
for details. Similarly, the limiting squared exponential GP satisfies

$e^{-\frac{\kappa^{2}}{4}\Delta}f=\mathcal{W}$ (7)

where $e^{-\frac{\kappa^{2}}{4}\Delta}$ is the (rescaled) heat semigroup [12,
16]. These equations have been studied as a means to extend Matérn Gaussian
processes to Riemannian manifolds, such as the sphere and torus [23, 6]. Since
these spaces can be discretized to form a graph, we explore the relationship
between these settings in the sequel.

###  2.2 Gaussian processes on graphs

A number of approaches have been proposed to define GPs over a weighted
undirected graph $G=(V,E)$. The aim is to define a GP indexed by the vertices
$V$, which reflects the notion of _closeness_ induced by the edges $E$ and
their associated weights. This has been studied by a variety of authors,
including [20, 30], and others, in areas such as GMRFs and diffusion kernels.
A number of variations, such as introducing dependence on attributes contained
in nodes [26], are possible. We will interpret the ideas presented in this
work from multiple viewpoints in order to synthesize these different
perspectives.

A number of authors, such as [40, 46] have studied graph-structured multi-
output Gaussian processes $f:\mathbb{R}\rightarrow\mathbb{R}^{|V|}$, where
$|V|$ is the number of nodes for which the outputs dependencies should reflect
the graph structure. Although such settings may appear different from the
scope outlined above, the problem can be cast into the proposed graph GP
framework by constructing a Gaussian process $f:\mathbb{R}\times
V\rightarrow\mathbb{R}$ through partial function application.

It is also possible to consider Gaussian processes
$f:\mathcal{G}\rightarrow\mathbb{R}$ where $\mathcal{G}$ is an appropriately
defined _space of graphs_. Here, each individual input to the GP is an _entire
graph_ rather than just a node on a fixed graph. This setting departs
significantly from the preceding ones, and we do not study it in this work: a
recent survey on the topic is given by [21].

##  3 Matérn GPs on graphs

We now define the _Matérn_ family of Gaussian processes on graphs by
generalizing their SPDE characterization to the graph setting. This will
entail introducing appropriate notions for the left-hand-side and right-hand-
side of equations (6) and (7). Note that since a graph is a finite set, such a
Gaussian process can be viewed as a multivariate Gaussian whose indices are
the graph’s nodes. Throughout this work, we refer to these as Gaussian
processes to emphasize that their covariance reflects the structure of the
graph.

Let $G$ be a weighted undirected graph whose weights are non-negative—for an
unweighted graph, assume all weights are equal to one. Denote its adjacency
matrix by $\mathbf{W}$, its diagonal degree matrix by $\mathbf{D}$ with
$\mathbf{D}_{ii}=\sum_{j}W_{ij}$, and define the _graph Laplacian_ by

$\mathbf{\Delta}=\mathbf{D}-\mathbf{W}.$ (8)

The graph Laplacian is a symmetric, positive semi-definite matrix, which we
view as a linear operator acting on a $|V|$-dimensional real space. Note that
this operator should be viewed as an analog of $-\Delta$ in the Euclidean or
Riemannian setting.222This sign convention ensures the positive semi-
definiteness of $\mathbf{\Delta}$ and is equivalent to the analyst’s (rather
than geometer’s) convention for studying the continuous Laplacian.
Notwithstanding the different sign convention, the graph and Euclidean
Laplacian are intimately linked. For example, the Laplacian of a graph given
by a grid corresponds to the finite-difference approximation of the Euclidean
Laplacian, and diffusions on graphs involve the graph Laplacian in the same
way diffusion in continuous spaces involve the classical Laplacian. If $G$ is
connected, then $\mathbf{\Delta}$ has a nullspace of dimension $1$, following
the case of compact manifolds [32]. Since $\mathbf{\Delta}$ is symmetric
positive semi-definite, it admits an eigenvalue decomposition
$\mathbf{\Delta}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$ where
$\mathbf{\Lambda}$ is diagonal with non-negative entries and $\mathbf{U}$ is
orthogonal.

To construct graph-theoretic analogs of the differential operators in (6) and
(7), we introduce a notion of _functional calculus_ for $\mathbf{\Delta}$. Let
$\Phi:\mathbb{R}\rightarrow\mathbb{R}$ be a function. For the diagonal matrix
$\mathbf{\Lambda}$, let $\Phi(\mathbf{\Lambda})$ be a diagonal matrix defined
by applying $\Phi$ to the diagonal of $\Lambda$ element-wise. Define the
matrix $\Phi(\mathbf{\Delta})$ to be

$\Phi(\mathbf{\Delta})=\mathbf{U}\Phi(\mathbf{\Lambda})\mathbf{U}^{T}.$ (9)

Although such a definition may seem arbitrary, this generalization of $\Phi$
to square matrices can be interpreted intuitively as plugging
$\mathbf{\Delta}$ into the Taylor expansion of $\Phi$, which immediately boils
down to (9) due to the orthogonality of $\mathbf{U}$. Taking $\Phi$ to be one
of

$\displaystyle\Phi(\lambda)$
$\displaystyle=\left(\frac{2\nu}{\kappa^{2}}+\lambda\right)^{\frac{\nu}{2}}$
$\displaystyle\Phi(\lambda)$ $\displaystyle=e^{\frac{\kappa^{2}}{4}\lambda}$
(10)

gives the operators on the left-hand side of (6) and (7), respectively.333We
use
${\Phi(\lambda)=\big{(}\frac{2\nu}{\kappa^{2}}+\lambda\big{)}^{\frac{\nu}{2}}}$
and ${\Phi(\lambda)=e^{\frac{\kappa^{2}}{4}\lambda}}$ instead of
${\Phi(\lambda)=\big{(}\frac{2\nu}{\kappa^{2}}-\lambda\big{)}^{\frac{\nu}{2}}}$
and ${\Phi(\lambda)=e^{-\frac{\kappa^{2}}{4}\lambda}}$ because of the
different sign convention for the graph Laplacian. Note that the term $d/4$
present in the Euclidean and Riemannian cases to ensure regularity is not
needed here. This will result in a slightly different scaling for graph Matérn
kernels, which for Markovian cases are indexed by integers rather than the
half-integers.

Figure 1: Here we illustrate prior variance for a complete graph, a star
graph, and two randomly generated graphs. For the complete graph, the variance
is uniform because of its symmetry (reordering the nodes always gives the same
graph). In the star graph, the center node has much lower variance, this
mirrors the behavior of random walk for which the return time is lower for the
center node than for the other nodes. The random graphs illustrate the idea
that variance depends on graph structure rather than on the degrees of
individual nodes.

Replacing Gaussian white noise process with a standard Gaussian
$\bc{W}\sim\operatorname{N}(\boldsymbol{0},\mathbf{I})$ gives the equations

$\displaystyle\left(\frac{2\nu}{\kappa^{2}}+\mathbf{\Delta}\right)^{\frac{\nu}{2}}\boldsymbol{f}$
$\displaystyle=\bc{W},$ $\displaystyle
e^{\frac{\kappa^{2}}{4}\mathbf{\Delta}}\boldsymbol{f}$ $\displaystyle=\bc{W},$
(11)

where we associate the vector of coefficients $\boldsymbol{f}$ with the
function $f:V\rightarrow\mathbb{R}$, which gives the graph Gaussian process of
interest. This gives

$\displaystyle\boldsymbol{f}$
$\displaystyle\sim\operatorname{N}\left(\boldsymbol{0},\left({\textstyle\frac{2\nu}{\kappa^{2}}}+\mathbf{\Delta}\right)^{-\nu}\right),$
$\displaystyle\boldsymbol{f}$
$\displaystyle\sim\operatorname{N}\left(\boldsymbol{0},e^{-\frac{\kappa^{2}}{2}\mathbf{\Delta}}\right)$
(12)

as _graph Matérn_ and _graph squared exponential_ Gaussian processes,
respectively. Note that the former covariance is _not_ obtained by adding
$\frac{2\nu}{\kappa^{2}}$ to the entries of the matrix $\mathbf{\Delta}$, but
rather by adding this number to its _eigenvalues_ , in the sense of Equation
(9). By construction their covariance matrices are positive semi-definite, and
we refer to them as the _graph Matérn_ and _graph diffusion_ kernels,
respectively. These possess a number of key properties, which we enumerate
below.

#### Sparsity.

For sparse graphs, $\mathbf{\Delta}$ is sparse. Hence, for sufficiently small
integers $\nu$ and most graphs, the precision matrices
$\left(\frac{2\nu}{\kappa^{2}}+\mathbf{\Delta}\right)^{\nu}$ are sparse. This
can be exploited to significantly accelerate their computational efficiency,
and is utilized extensively by prior work on GMRFs [30].

#### Non-uniform variance.

The prior variance for the introduced kernels varies along vertices. This
behavior is similar to the variance of the Matérn kernel on certain non-
homogeneous Riemannian manifolds [6]. Figure 1 illustrates the prior variance
of the graph Matérn kernel on a set of graphs. It can be seen that the
kernel’s variance is not simply a function of degree and depends in a complex
manner on the graph in question.

[38] study a similar phenomenon in the context of _random walk kernels_. They
show that the variance is determined by the return time of a certain random
walk defined on the graph.

#### Use with symmetric normalized graph Laplacian.

The above kernels are defined in terms of the graph Laplacian
$\mathbf{\Delta}$. In some applications, it might be preferable to instead
work with the _symmetric normalized graph Laplacian_
$\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}$. Doing so in above
expressions yields

$\displaystyle\boldsymbol{f}$
$\displaystyle\sim\operatorname{N}\left(\boldsymbol{0},\left({\textstyle\frac{2\nu}{\kappa^{2}}}+\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}\right)^{-\nu}\right)$
(13) $\displaystyle\boldsymbol{f}$
$\displaystyle\sim\operatorname{N}\left(\boldsymbol{0},e^{-\frac{\kappa^{2}}{2}\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}}\right)$
(14)

which we call the _symmetric normalized graph Matérn GP_ and the _symmetric
normalized graph squared exponential GP_ , respectively. Whether one should
use the symmetric normalized or the standard Laplacian generally depends on
the particular application—see [42, 412] for a discussion in the related
setting of spectral clustering.

#### Connection with the graph diffusion equation.

The graph diffusion kernel $e^{-\frac{\kappa^{2}}{2}\mathbf{\Delta}}$ is the
Green’s function of the graph diffusion equation. That is, if
$\boldsymbol{\phi}:[0,\infty)\times V\to\mathbb{R}$ solves the differential
equation

$\displaystyle\left.\frac{\mathop{}\\!\mathrm{d}{\boldsymbol{\phi}}}{\mathop{}\\!\mathrm{d}{t}}\right.+\mathbf{\Delta}\boldsymbol{\phi}=0$
$\displaystyle\left.\boldsymbol{\phi}\right|_{t=0}=\boldsymbol{v}$ (15)

then
$\left.\boldsymbol{\phi}\right|_{t=\tau}=e^{-\tau\mathbf{\Delta}}\boldsymbol{v}$.
This equation describes heat transfer along the graph. If $\mathbf{\Delta}$ is
replaced with the symmetric normalized graph Laplacian
$\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}$, then the value
$\left.\boldsymbol{\phi}\right|_{t=\tau}$ can be interpreted as the
unnormalized density of a continuous-time random walk which moves along the
graph.

#### Limits and connection with random walks.

Figure 2: The Matérn kernel $k_{3/2}(x,\cdot)$ defined on increasingly fine
triangular meshes (three left columns), and Matérn kernel $k_{1/2}(x,\cdot)$
on the sphere (rightmost). The true kernel is computed using 512 analytically
known Laplace-Beltrami eigenpairs via the technique in [6]. The discrepancy in
the smoothness parameter is due to a different parameterization of the kernel,
which accounts for dimensionality of the Riemannian manifold. The point $x$ is
chosen to correspond to the north pole. For graphs defined by the meshes,
kernel values are projected onto the sphere via piecewise-constant
interpolation based on spherical Voronoi diagram induced by the graph nodes.
The graph kernel converges to the true kernel on the right.

We detail two different ways in which the graph diffusion kernel arises as a
limit of a sequence of kernels.

First, mirroring the Euclidean case, the graph Matérn kernel converges to the
graph diffusion kernel as $\nu\rightarrow\infty$, when scaled appropriately.
To see this, note that the eigenvectors of both kernels are identical, and
consider the eigenvalues

$(2\nu/\kappa^{2})^{\nu}\left(\frac{2\nu}{\kappa^{2}}+\lambda\right)^{-\nu}\xrightarrow[]{\nu\rightarrow\infty}e^{-\frac{\kappa^{2}}{2}\lambda}.$
(16)

This implies that the corresponding matrices converge, and hence that the
Gaussian processes converge in distribution.

Second, the graph diffusion kernel arises as a limit of the _random walk
kernel_ of [32]. This kernel is defined as

$\big{(}\mathbf{I}-(1-\alpha)\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}\big{)}^{p},$
(17)

which arises by symmetrizing the $p$-step transition matrix of a lazy random
walk on the graph. This kernel looks superficially similar to the Matérn
kernel defined previously, but contains a number of important differences. In
particular, the power $p$ is positive rather than negative, the laziness
parameter $\alpha$ is restricted to be in $[0,1)$, and the Laplacian is
subtracted rather than added. Nonetheless, this kernel also converges to the
diffusion kernel, in the limit $p\rightarrow\infty$ with
$\alpha=1-\frac{\kappa^{2}}{2p}$ where $\kappa$ is a fixed constant.
Specifically,

$\displaystyle\begin{gathered}\big{(}\mathbf{I}-(1-\alpha)\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}\big{)}^{p}\\\
\qquad\qquad\xrightarrow[]{p\rightarrow\infty}e^{-\frac{\kappa^{2}}{2}\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}}\end{gathered}$
(20)

which yields the normalized graph diffusion kernel.

#### Relationship with Matérn Gaussian processes on Riemannian manifolds.

One of the most common ways that linear (S)PDEs in the Euclidean and
Riemannian settings are solved is via the _finite-element method_. Roughly
speaking, this technique expresses the (S)PDE solution using a truncated basis
expansion, whose coefficients are obtained by solving a sparse linear system.
A typical choice for the finite element basis are _piecewise polynomial_
functions constructed using the nodes and edges of a graph, which corresponds
to a mesh discretization of the manifold.

By their construction, the basis weights of a finite-element-discretized
Riemannian Matérn Gaussian process define a GMRF approximation to the Matérn
Gaussian process. The precise form of the GMRF approximation obtained depends
on the choice of finite element space. Under appropriate regularity
conditions, [23] have proven convergence of the finite element discretizations
to the solution of the true SPDE solution.

Not all kernels that arise this way are graph Matérn kernels. Fortunately, a
number of other convergence results [4, 7] are available. One can show under
appropriate regularity conditions and proper weighting of edges that the
eigenvalues and eigenvectors of the graph Laplacian converge to the
eigenvalues and eigenfunctions of the Laplace–Beltrami operator. Using these
ideas, [31] provide a theoretical framework for studying how graph Matérn
kernels converge to their Riemannian counterparts. This is illustrated in
Figure 2, which shows convergence of a graph Matérn kernel to its Riemannian
limit for a sphere.

Graph Matérn kernels can also be used with graphs that are obtained as
intermediate steps of manifold learning algorithms. This enables them to be
applied to data that is concentrated on an unknown manifold in a way that
captures its geometry—a common situation in practice. Provided the manifold
learning algorithm is well-posed and possesses appropriate guarantees, the
resulting graph kernels will converge to the manifold Matérn kernels studied
by [6].

#### Use as an output cross-covariance.

The above section describes Matérn Gaussian processes to model functions
$f:V\rightarrow\mathbb{R}$. We now briefly recall how to lift this
construction to functions
$\boldsymbol{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{|V|}$, where the output
cross-covariance is induced by the graph $G$. The simplest way to do so is to
apply the _intrinsic coregionalization model_ [2, 14]. Here, one views a
multi-output function $\boldsymbol{f}(x)$ as a single-output function with an
additional input $f(x,i)$ such that $\boldsymbol{f}(x)_{i}=f(x,i)$. One then
uses a separable kernel

$k((x,i),(x^{\prime},i^{\prime}))=k_{\mathbb{R}^{d}}(x,x^{\prime})k_{G}(i,i^{\prime}),$
(21)

where $k_{G}(i,i^{\prime})$ is the graph kernel. This defines model that
respects the graph structure over the output space.

###  3.1 Scalable Training

Here we discuss how to train Matérn Gaussian processes in a scalable manner.
The techniques we present fall in one of two broad areas: (1) graph Fourier
feature approximations, and (2) efficient computation via sparse linear
algebra. All of the techniques described here apply even if the graph is too
large to store its full covariance matrix in memory.

Graph Fourier feature methods approximate the kernel matrix using a truncated
eigenvalue expansion. Here, we first obtain the $\ell$ smallest eigenvalues
and eigenvectors of $\mathbf{\Delta}$ using a routine designed for efficiently
computing the eigenvalues of sparse matrices, such as the Lanczos algorithm.
Since the mappings

$\displaystyle\lambda$
$\displaystyle\mapsto{\frac{2\nu}{\kappa^{2}}+\lambda}^{-\nu}$
$\displaystyle\lambda$ $\displaystyle\mapsto e^{-\frac{\kappa^{2}}{4}\lambda}$
(22)

are decreasing functions on $\mathbb{R}^{+}$, we thus obtain the $\ell$
largest eigenvalues of the Matérn or diffusion kernels. Once this is obtained,
scalable training proceeds by relating the approximate GP with a Bayesian
linear regression model, mirroring Euclidean Fourier features [28]. Here,
reduced computational costs come from the orthogonality of features, i.e. the
eigenvectors of the kernel matrix. This mirrors ideas in [11, 33, 6] where
Gaussian processes are specified using Karhunen–Loéve type decompositions—for
further details, see Appendix B. The main difficulty with this approach is
that it can exhibit _variance starvation_ issues, which causes approximation
quality to deteriorate in the large-data regime [43].

A different way to train graph GPs scalably is by restricting $\nu$ to small
integer values, for which in most graphs the inverse kernel matrices, termed
_precision_ matrices

$\left(\frac{2\nu}{\kappa^{2}}+\mathbf{\Delta}\right)^{\nu}$ (23)

are sparse. The prior is then a _Gaussian Markov random field_ for which a
variety of scalable training procedures are already available off-the-shelf
[30]. This approach generally avoids variance starvation, at a cost often no
worse than that of graph Fourier features. The main limitation is that the
computational gains diminish when $\nu$ increases, and that the numerical
routines needed for scalable training may not be well-supported by standard
Gaussian process packages. This support is rapidly improving: as shown by [10,
1], it is possible to leverage certain sparse computations in automatic
differentiation frameworks, such as TensorFlow, to accelerate variational GP
models and their sparse counterparts.

(a) Mean

(b) Standard deviation

Figure 3: Traffic flow speed (mph) interpolation over a graph of San Jose
highways. Colored circles represent graph nodes, black lines represent edges
and colored circles with a white point in the middle represent training data.
The standard deviation color bar is clipped: the $10\%$ most variable points
are painted in yellow.

(a) Mean

(b) Standard deviation

Figure 4: Traffic flow speed (mph) interpolation over a particular road
junction on a graph of San Jose highways. Colored circles represent graph
nodes, black lines represent edges and colored circles with a white point in
the middle represent training data.

#### Non-conjugate learning via doubly stochastic variational inference.

(a) Ground truth

(b) Prediction errors. Red nodes are misclassified.

Figure 5: Topic classification for the _Cora_ citation network, using the
citation graph.

In non-conjugate settings, such as classification, the posterior distribution
is no longer Gaussian, and the training techniques from the preceding sections
do not apply. Here we briefly sketch how to obtain analogous techniques using
variational inference and inducing points.

Let $\boldsymbol{z}\subset V$ be a set of variational _inducing inputs_ [18],
let $\boldsymbol{x}\subset V$ be the training locations, let
$\boldsymbol{u}=\boldsymbol{f}(\boldsymbol{z})$. We approximate the
distribution $p(\boldsymbol{u},f\mid\boldsymbol{y})$ with a variational family
$q(\boldsymbol{u},f)$ of the form
$q(\boldsymbol{u},f)=p(f\mid\boldsymbol{u})q(\boldsymbol{u})$, and train the
model by minimizing the Kullback–Leibler divergence

$\displaystyle D_{\operatorname{KL}}(q(\boldsymbol{u},f)\mid\mid
p(\boldsymbol{u},f\mid\boldsymbol{y}))=\mathop{\mathchoice{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}missing}_{q(\boldsymbol{u},f)}\\!\ln\frac{q(\boldsymbol{u},f)}{p(\boldsymbol{u},f\mid\boldsymbol{y})}$
(24) $\displaystyle=D_{\operatorname{KL}}(q(\boldsymbol{u})\mid\mid
p(\boldsymbol{u}))-\sum_{i=1}^{N}\mathop{\mathchoice{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}{\mathop{\mathbb{E}}}missing}_{q(\boldsymbol{u},f)}\ln
p(y_{i}\mid f_{i})$ (25) $\displaystyle\quad-\ln p(\boldsymbol{y})$

where we henceforth ignore the term $\ln p(\boldsymbol{y})$ because it does
not depend on $q$.

Following [36], we choose
$q(\boldsymbol{u})=\operatorname{N}(\boldsymbol{\mu},\mathbf{\Sigma})$ to be a
free-form multivariate Gaussian. Minimizing the Kullback–Leibler divergence
therefore entails finding the closest _GP posterior_ to the true (non-
Gaussian) posterior. The inducing points $\boldsymbol{z}$ can be chosen to be
a subset of the data, or even the entire graph, provided one works with sparse
precision matrices.

To minimize the Kullback–Leibler divergence, one employs a stochastic
optimization algorithm, such as ADAM. At each step, this entails sampling (1)
a mini-batch of data, (2) the variational posterior from $q(\boldsymbol{u},f)$
so that one can calculate the expression inside the expectation. For
classification, this expression reduces to a cross-entropy loss. This yields
the correct loss in expectation, thereby defining a doubly stochastic
variational inference [37] algorithm for scalable training.

#### Variational inference with interdomain inducing variables.

An alternative to using $\boldsymbol{u}=\boldsymbol{f}(\boldsymbol{z})$ as
inducing variables is to consider _interdomain inducing variables_ such as
$u_{i}=\langle f,\psi_{i}\rangle$. Provided that the $\psi_{i}$ are chosen
wisely, this approach can yield orders of magnitude speed ups [17]. For the
problem at hand, it is natural to use the kernel eigenvectors as inducing
functions $\phi_{i}$ [8]. It is worth noting that this construction is exactly
the discrete counterpart to the _variational inference with spherical
harmonics_ method developed by [11] in the continuous case.

##  4 Experiments

This section is dedicated to examples that illustrate the relevance of the
proposed approach. Full experimental details are provided in Appendix A.

###  4.1 Probabilistic graph interpolation of traffic data

To illustrate the flexibility of graph GPs in settings where other GP models
may not be easy to apply, we consider probabilistic interpolation of traffic
congestion data of highways in the city of San Jose, California. The graph of
the road network is obtained from OpenStreetMap (OSM, OpenStreetMap), and it
is pre-processed to associate to each edge a weight equal to its inverse
weight and by removing hanging nodes. We obtain traffic congestion data from
the _California Performance Measurement System_ database [9], which provides
traffic flow speed in miles per hour at different locations, which we add as
additional nodes onto the graph. Note that both travel directions on a highway
are considered separate roads and thus form separate edges—this can be seen in
Figure 4.

This process yields a graph with $1,016$ nodes and $1,173$ edges. Traffic flow
speed is available at $325$ nodes. We use $250$ of these nodes, chosen at
random, as training data, and the remainder as test data. The data is
normalized by subtracting the mean over the training set and dividing the
result by its standard deviation. For output values, we examine traffic
congestion on Monday at 17:30 in the afternoon.

We train the Matérn graph GP model, whose kernel is approximated using $500$
eigenpairs of the graph Laplacian. During training, we optimize the kernel
hyperparameters $\kappa$ (length scale), $\sigma^{2}$ (variance) and $\nu$
(smoothness), as well as the Gaussian likelihood error variance
$\varsigma^{2}$.

Figure 3 shows the predicted mean traffic on the entire graph, along with the
standard deviation. This illustrates that global dependencies are captured by
the model in different regions of the graph. In particular, the posterior GP’s
standard deviation increases for nodes for which traffic sensors are further
away.

To illustrate local dependencies, Figure 4 details the model prediction at a
particular road junction. It can be seen that despite being very close in the
Euclidean distance, two locations on opposite sides of the road can have
associated predictions that differ considerably. This shows that the graph GP
incorporates the structure of the graph into its predictions, and allows them
to differ when the graph distance is large.

###  4.2 Multi-class classification in a scientific citation network

To illustrate graph GPs in a more complex non-conjugate setting, we consider
multi-class classification on the largest connected component of the _Cora_
citation network, considered purely as a graph, i.e. without abstracts. This
yields a data set consisting of $2,485$ nodes corresponding to scientific
publications, with $5,069$ edges corresponding to citations. We choose $140$
nodes at random for training. Each node is labeled according to one of seven
classes, which are different scientific topics. We employ only the graph and
labels as part of the classification problem.

For the model, we use a multi-output Matérn graph GP where each output
dimension is given using a separable kernel, yielding a GP
$f:V\rightarrow\mathbb{R}^{7}$. We employ a categorical likelihood, with class
probabilities given by applying a robust max function (Appendix A) to the GP.
We approximate the graph kernel using $500$ eigenpairs of the graph Laplacian,
in order to allow its pointwise evaluation, and train the model using
variational inference in GPflow [25, 39].

Performance is illustrated in Figure 5. The fraction of misclassified nodes is
small, showing that in this data set, a paper’s topic can be accurately
predicted from its works cited alone even in a low-data regime.

###  4.3 Comparison with other kernels

| Lap. | Matérn kernel | Diffusion kernel | Random walk kern. | Inverse cosine kern.
---|---|---|---|---|---
Traffic MSE | $\mathbf{\Delta}$ | 1.37 (0.08) | 1.38 (0.08) | — | —
$\mathbf{\widetilde{\Delta}}$ | 1.48 (0.11) | 1.49 (0.11) | 1.52 (0.12) | 1.92 (0.11)
Citation accuracy | $\mathbf{\Delta}$ | 0.79 (0.01) | 0.79 (0.02) | — | —
$\mathbf{\widetilde{\Delta}}$ | 0.79 (0.01) | 0.78 (0.01) | 0.78 (0.02) | 0.48 (0.02)
Table 1: Mean and standard deviation (in parentheses) of model performance for
different graph GP kernels. Note that the standard deviation of the ground
truth traffic flow speed across the entire graph is approximately $17.1$ mph.
We denote the unnormalized Laplacian by $\mathbf{\Delta}$, and the symmetric
normalized Laplacian by
$\mathbf{\widetilde{\Delta}}=\mathbf{D}^{-1/2}\mathbf{\Delta}\mathbf{D}^{-1/2}$.

We now evaluate the performance of GPs induced by the graph Matérn on both of
the examples considered. For comparison, we consider the graph diffusion
kernel, which is the limit of the Matérn kernel as $\nu\rightarrow\infty$, as
well as the random walk and inverse cosine kernels defined in [32]. We
consider variants employing both normalized and unnormalized graph Laplacians.
We evaluate performance via mean squared error, using $250$ training nodes
with $75$ test nodes for the traffic example, and $140$ training nodes with
$1000$ test nodes for the citation network example. We repeat this $10$ times
to obtain average mean squared error and accuracy. Results are presented in
Table 1. This shows that use of the more flexible Matérn class results in
similar or slightly improved performance compared to alternatives on these
data sets.

##  5 Conclusion

In this work, we study graph Matérn Gaussian processes, which are defined as
analogous to the Matérn stochastic partial differential equations used in the
Euclidean and Riemannian cases. We discuss a number of their properties, and
provide scalable training algorithms by means of Gaussian Markov random
fields, graph Fourier features and sparse GPs. We demonstrate their
effectiveness on a simple graph interpolation problem, and on a multi-class
graph classification problem. Compared to other graph kernels, the graph
Matérn class is flexible, interpretable, and shown to perform well, mirroring
the behavior of Euclidean and Riemannian Matérn kernels. We hope these
techniques inspire new use of Gaussian processes by machine learning
practitioners.

## Acknowledgments

We are grateful to Mohammed Charrout for finding a bug in our software, which
enabled us to fix incorrect results in Table 1 that were present in prior
drafts of this manuscript. This research was partially supported by “Native
towns”, a social investment program of PJSC Gazprom Neft and by the Ministry
of Science and Higher Education of the Russian Federation, agreements No
075-15-2019-1619 and No 075-15-2019-1620.

## References

  * [1] Vincent Adam et al.  “Doubly sparse variational Gaussian processes”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 2874–2884
  * [2] Mauricio A Alvarez, Lorenzo Rosasco and Neil D Lawrence  “Kernels for vector-valued functions: A review”  In _arXiv:1106.6251_ , 2011
  * [3] Ethan Anderes, Jesper Møller and Jakob G. Rasmussen  “Isotropic covariance functions on graphs and their edges”  In _Annals of Statistics_ 48.4 Institute of Mathematical Statistics, 2020, pp. 2478–2503
  * [4] Mikhail Belkin and Partha Niyogi  “Convergence of Laplacian eigenmaps”  In _Advances in Neural Information Processing Systems_ , 2007
  * [5] Geoff Boeing  “OSMnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks”  In _Computers, Environment and Urban Systems_ 65 Elsevier, 2017, pp. 126–139
  * [6] Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky and Marc Peter Deisenroth  “Matérn Gaussian processes on Riemannian manifolds”  In _Advances in Neural Information Processing Systems_ , 2020
  * [7] Dmitri Burago, Sergei Ivanov and Yaroslav Kurylev  “A graph discretization of the Laplace–Beltrami operator”  In _Journal of Spectral Theory_ 4.4 European Mathematical Society Publishing House, 2014, pp. 675–714
  * [8] David R. Burt, Carl E. Rasmussen and Mark Wilk  “Variational orthogonal features”  In _arXiv:2006.13170_ , 2020
  * [9] Chao Chen et al.  “Freeway performance measurement system: mining loop detector data”  In _Transportation Research Record_ 1748.1, 2001, pp. 96–102
  * [10] Nicolas Durrande et al.  “Banded matrix operators for Gaussian Markov models in the automatic differentiation era”  In _International Conference on Artificial Intelligence and Statistics_ , 2019
  * [11] Vincent Dutordoir, Nicolas Durrande and James Hensman  “Sparse Gaussian processes with spherical harmonic features”  In _International Conference on Machine Learning_ , 2020
  * [12] Lawrence C Evans  “Partial Differential Equations”  American Mathematical Society, 2010
  * [13] Aasa Feragen, Francois Lauze and Soren Hauberg  “Geodesic exponential kernels: When curvature and linearity conflict”  In _Conference on Computer Vision and Pattern Recognition_ , 2015
  * [14] Pierre Goovaerts  “Geostatistics for Natural Resources Evaluation”  Oxford University Press, 1997
  * [15] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik  “Table of Integrals, Series, and Products”  Academic Press, 2014
  * [16] Alexander Grigoryan  “Heat kernel and analysis on manifolds”  American Mathematical Society, 2009
  * [17] James Hensman, Nicolas Durrande and Arno Solin  “Variational Fourier features for Gaussian processes”  In _Journal of Machine Learning Research_ 18, 2017, pp. 151–202
  * [18] James Hensman, Nicolo Fusi and Neil D. Lawrence  “Gaussian processes for big data”  In _Uncertainty in Artificial Intelligence_ , 2013
  * [19] C.. Kaufman and B.. Shaby  “The role of the range parameter for estimation and prediction in geostatistics”  In _Biometrika_ 100.2 Oxford University Press, 2013, pp. 473–484
  * [20] R. Kondor and J. Lafferty  “Diffusion kernels on graphs and other discrete input spaces”  In _International Conference on Machine Learning_ , 2002
  * [21] Nils M Kriege, Fredrik D Johansson and Christopher Morris  “A survey on graph kernels”  In _Applied Network Science_ 5.1, 2020, pp. 1–42
  * [22] Mikhail Lifshits  “Lectures on Gaussian Processes”  Springer, 2012
  * [23] Finn Lindgren, Håvard Rue and Johan Lindström  “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73.4, 2011, pp. 423–498
  * [24] Anton Mallasto and Aasa Feragen  “Wrapped Gaussian process regression on Riemannian manifolds”  In _Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 5580–5588
  * [25] Alexander G. G. Matthews et al.  “GPflow: A Gaussian process library using TensorFlow”  In _Journal of Machine Learning Research_ 18.40, 2017, pp. 1–6
  * [26] Yin Cheng Ng, Nicolò Colombo and Ricardo Silva  “Bayesian semi-supervised learning with graph Gaussian processes”  In _Advances in Neural Information Processing Systems_ , 2018, pp. 1683–1694
  * [27] OpenStreetMap contributors  “Planet dump retrieved from https://planet.osm.org ”, 2017
  * [28] Ali Rahimi and Benjamin Recht  “Random features for large-scale kernel machines”  In _Advances in Neural Information Processing Systems_ , 2008
  * [29] Carl Edward Rasmussen and Christopher KI Williams  “Gaussian Processes for Machine Learning”  MIT Press, 2006
  * [30] Havard Rue and Leonhard Held  “Gaussian Markov Random Fields: Theory and Applications”  CRC Press, 2005
  * [31] Daniel Sanz-Alonso and Ruiyi Yang  “The SPDE approach to Matérn fields: Graph representations”  In _arXiv:2004.08000_ , 2020
  * [32] Alexander J Smola and Risi Kondor  “Kernels and regularization on graphs”  In _Learning Theory and Kernel Machines_ Springer, 2003, pp. 144–158
  * [33] Arno Solin and Simo Särkkä  “Hilbert space methods for reduced-rank Gaussian process regression”  In _Statistics and Computing_ 30.2 Springer, 2020, pp. 419–446
  * [34] Michael Stein  “Asymptotics for spatial processes”  In _Handbook of spatial statistics_ CRC Press, 2010, pp. 79–88
  * [35] Michael L Stein  “Interpolation of Spatial Data: Some Theory for Kriging”  Springer Science & Business Media, 1999
  * [36] Michalis Titsias  “Variational learning of inducing variables in sparse Gaussian processes”  In _International Conference on Artificial Intelligence and Statistics_ , 2009
  * [37] Michalis Titsias and Miguel Lázaro-Gredilla  “Doubly stochastic variational Bayes for non-conjugate inference”  In _International Conference on Machine learning_ , 2014, pp. 1971–1979
  * [38] Matthew J Urry and Peter Sollich  “Random walk kernels and learning curves for Gaussian process regression on random graphs”  In _Journal of Machine Learning Research_ 14.1, 2013, pp. 1801–1835
  * [39] Mark van der Wilk et al.  “A Framework for interdomain and multioutput Gaussian processes”  In _arXiv:2003.01115_ , 2020
  * [40] Arun Venkitaraman, Saikat Chatterjee and Peter Handel  “Gaussian processes over graphs”  In _International Conference on Acoustics, Speech and Signal Processing_ , 2020
  * [41] S Vichy N Vishwanathan, Nicol N Schraudolph, Risi Kondor and Karsten M Borgwardt  “Graph kernels”  In _Journal of Machine Learning Research_ 11, 2010, pp. 1201–1242
  * [42] Ulrike Von Luxburg  “A tutorial on spectral clustering”  In _Statistics and Computing_ 17.4 Springer, 2007, pp. 395–416
  * [43] Zi Wang, Clement Gehring, Pushmeet Kohli and Stefanie Jegelka  “Batched large-scale bayesian optimization in high-dimensional spaces”  In _International Conference on Artificial Intelligence and Statistics_ , 2018
  * [44] Peter Whittle  “Stochastic processes in several dimensions”  In _Bulletin of the International Statistical Institute_ 40.2, 1963, pp. 974–994
  * [45] James T. Wilson et al.  “Efficiently sampling functions from Gaussian process posteriors”  In _International Conference on Machine Learning_ , 2020
  * [46] Yin-Cong Zhi, Yin Cheng Ng and Xiaowen Dong  “Gaussian processes on graphs via spectral kernel learning”  In _arXiv:2006.07361_ , 2020

##  Appendix A Appendix: experimental details

### Probabilistic graph interpolation of traffic data

Here we consider the problem of interpolating traffic congestion data over a
road network, consisting of San Jose highways. We build the road network, a
weighted graph, from OpenStreetMap data [27] obtained using the OSMnx library
[5].

We obtain traffic congestion data from the _California Performance Measurement
Systems_ database [9]. This system maps sensor location and time & date to
(among other factors) traffic flow speed in miles per hour measured by the
given sensor at a particular time and date. We chose the specific date and
time to be 17:30 on Monday, the 2th of January, 2017. We bind the traffic
congestion data to the graph by adding additional nodes that subdivide
existing edges of the graph at the location of the measurement points. The
edge weights are assigned to be equal to the inverse travel distances between
the corresponding nodes. When we encounter two nodes connected by multiple
edges, we remove all but one edge between them, and set its weight to be the
inverse mean travel distance between the pair of nodes. After this, the graph
is processed by (1) twice subsequently removing hanging nodes, and (2)
removing two other nodes located particularly far away from data. We end up
with $325$ labeled nodes on a sparse connected graph with $1016$ nodes and
$1173$ edges.

For the model, we employ a GP with a graph Matérn kernel and the graph
diffusion kernel. Both kernels are approximated using $500$ eigenpairs of the
graph Laplacian. The eigenpairs required to build graph kernels are computed
using TensorFlow’s eigenvalue decomposition (tf.linalg.eigh) routine and are
thus exact up to a numerical error arising from the floating point arithmetic.
We train by optimzing the kernel hyperparameters $\kappa$ (length scale),
$\sigma^{2}$ (variance) and for the case of the graph Matérn kernel, $\nu$
(smoothness), as well as the Gaussian likelihood variance $\varsigma^{2}$. We
use the following values for initialization: $\nu=3/2$ for the graph Matérn
kernel, $\kappa=3$, $\sigma^{2}=1$, $\varsigma^{2}=0.01$. We randomly choose
$250$ nodes with traffic flow rates as training data. The likelihood is
optimized for $20000$ iterations with the ADAM optimizer, with learning rate
set to the TensorFlow default value of $0.001$.

We repeat the experiment 10 times. The average MSE over the 10 experiments and
its standard deviation are presented in Table 1 for the graph Matérn kernel
and for the graph diffusion kernel. The smoothness parameter $\nu$ after
optimization varies across the runs. Its mean is $1.8$ and the standard
deviation is $0.9$. One of these ten runs is visualized in Figures 3 and 4.

### Multi-class classification in a scientific citation network

Here we consider multi-class classification on the largest connected component
of the _Cora_ citation network.

There are $2708$ nodes corresponding to scientific publications in the _Cora_
citation network. Usually, each node of this graph is attributed by a
$0$/$1$-valued word vector indicating the absence/presence of the
corresponding word from the dictionary of $1433$ unique words in the
abstract—we discard this portion of data to focus on other problem
aspects.444It’s also possible to use a separable kernel given by the product
of a graph Matérn kernel and a Euclidean kernel operating on word vectors, but
this model is unlikely to perform well due to usual GP difficulties in high
dimension. Incorporating this portion of the data effectively therefore
necessitates additional modeling considerations, so we do not focus on it
here. Two nodes are connected by an edge if one paper cites another. The edges
are unweighted and undirected, and there are $5429$ edges overall. The nodes
are labeled into one of seven classes. We concentrate on the largest connected
component of the Cora graph, which has $2485$ nodes and $5069$ edges.

We consider a Gaussian process $V\rightarrow\mathbb{R}^{7}$ with independent
components as a latent process for the classification task. The process is
then combined with the robust max function $\boldsymbol{\varphi}(\cdot)$ to
yield a vector of probabilities over the classes. Let
$\boldsymbol{f}=(f_{1},..,f_{7})$ denote the values of $7$-dimensional latent
process. Then we have

$\displaystyle\boldsymbol{\varphi}(\boldsymbol{f})$
$\displaystyle=(\varphi_{1}(\boldsymbol{f}),..,\varphi_{7}(\boldsymbol{f}))$
$\displaystyle\varphi_{c}(\boldsymbol{f})$
$\displaystyle=\begin{cases}1-\varepsilon,&\mathrel{\text{missing}}{if}c=\operatorname*{\arg\max}_{c}f_{c}\\\
\varepsilon/6,&\mathrel{\text{missing}}{otherwise},\\\ \end{cases}$ (26)

where $\varepsilon$ is fixed to be $10^{-3}$. For the likelihood, we employ
the categorical distribution. We use the graph Matérn kernel and the graph
diffusion kernel for the prior. Both kernels are approximated using $500$
eigenpairs of the graph Laplacian. The eigenpairs required are computed using
TensorFlow eigenvalue decomposition (tf.linalg.eigh) and are thus exact up to
a numerical error arising from use of floating point arithmetic. We use $140$
random nodes as a training set and predict labels for the remaining $2345$
nodes.

For the variational GP, we employ the standard SVGP model of GPflow [25, 39].
We fix the the coordinates of inducing points, which are located on the graph,
to the data locations, and do not optimize them. We restrict the covariance of
the inducing distribution to a diagonal form, since this results in the best
performance in this setting. The inducing mean is initialized to be zero and
the inducing covariance is initialized to be the identity. The prior
hyperparameters are initialized to be $\nu=3$ for the Matérn kernel,
$\kappa=5$, $\sigma^{2}=1$. The variational objective is optimized for $20000$
iterations with the ADAM optimizer whose learning rate set to $0.001$.
Convergence is usually achieved in a significantly lower number of iterations.
The SVGP batch size is equal to the number of training data points, and
whitening is enabled.

We repeat the experiment 10 times. Each time, we additionally randomly choose
a test set of $1000$ nodes. The average (over the 10 experiments) accuracy
(fraction of the correctly classified points) and its standard deviation are
presented in Table 1 for the graph Matérn kernel and for the graph diffusion
kernel. The value of the smoothness parameter $\nu$ after optimization varies
considerably across runs. The mean $\nu$ is $3.2$ and the standard deviation
is $2.7$. The results of one of said ten runs are visualized on Figures 6, 7,
8. Note that Figures 6 and 8 repeat Figure 5 (a), (b) at a higher resolution.

Figure 6: Topic classification for the _Cora_ citation network, using the
citation graph. Ground truth.

Figure 7: Topic classification for the _Cora_ citation network, using the
citation graph. Prediction.

Figure 8: Topic classification for the _Cora_ citation network, using the
citation graph. Prediction errors. Red nodes are misclassified.

##  Appendix B Appendix: graph Fourier features

Assume that $\mathbf{\Delta}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$ with
$\mathbf{U}$ orthogonal and $\mathbf{\Lambda}$ diagonal. For a set of of
locations $\boldsymbol{x}\in V^{n}$ define $\mathbf{U}(\boldsymbol{x})$ to be
the $n\times d$ submatrix of $\mathbf{U}$ with rows corresponding to the
elements of $\boldsymbol{x}$. Define also $\Psi(\lambda)=\Phi(\lambda)^{-2}$,
where $\Phi$ is given in (10). Then, if $f$ is a graph Matérn process and
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$ is the covariance matrix of
$f(\boldsymbol{x})$, we have

$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}=\mathbf{U}(\boldsymbol{x})\Psi(\mathbf{\Lambda})\mathbf{U}(\boldsymbol{x})^{T}\approx\widetilde{\mathbf{U}}(\boldsymbol{x})\Psi(\widetilde{\mathbf{\Lambda}})\widetilde{\mathbf{U}}(\boldsymbol{x})^{T},$
(27)

where $\widetilde{\mathbf{U}}(\boldsymbol{x})$ and
$\widetilde{\mathbf{\Lambda}}$ are the approximations of
$\widetilde{\mathbf{U}}(\boldsymbol{x})$ and $\widetilde{\mathbf{\Lambda}}$
with elements in the last $d-l$ rows of matrices
$\widetilde{\mathbf{U}}(\boldsymbol{x})$ and $\widetilde{\mathbf{\Lambda}}$
set to zero. This approximation arises from finding only the $l$ smallest
eigenpairs, which can be done more efficiently than the full
eigendecomposition by employing, for instance, the Lanczos algorithm.

With this, using the Woodbury matrix identity, we get

$\displaystyle\left(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I}\right)^{-1}$
$\displaystyle=\sigma^{-2}\mathbf{I}-\sigma^{-4}\mathbf{U}(\boldsymbol{x})\big{(}\Psi(\mathbf{\Lambda})^{-1}+\overbracket{\sigma^{-2}\mathbf{U}(\boldsymbol{x})^{T}\mathbf{U}(\boldsymbol{x})}^{\sigma^{-2}\mathbf{I}}\big{)}^{-1}\mathbf{U}(\boldsymbol{x})^{T},$
(28)
$\displaystyle\left(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I}\right)^{-1}$
$\displaystyle\approx\sigma^{-2}\mathbf{I}-\sigma^{-4}\widetilde{\mathbf{U}}(\boldsymbol{x})\big{(}\Psi(\widetilde{\mathbf{\Lambda}})^{-1}+\underbracket{\sigma^{-2}\widetilde{\mathbf{U}}(\boldsymbol{x})^{T}\widetilde{\mathbf{U}}(\boldsymbol{x})}_{\text{diagonal}}\big{)}^{-1}\widetilde{\mathbf{U}}(\boldsymbol{x})^{T},$
(29)

where a diagonal matrix is inside the brackets. It is clear that computing
$\left(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I}\right)^{-1}\boldsymbol{y}$
when the (approximation of) eigenpairs are known can be done without incurring
the $O(n^{3})$ cost, resulting in efficient computation.

The resulting kernel can be (approximately) expressed by

$k(i,j)=\sum_{s=1}^{d}\Psi(\lambda_{s})\boldsymbol{u_{s}}(i)\boldsymbol{u_{s}}(j)\approx\sum_{s=1}^{l}\Psi(\lambda_{s})\boldsymbol{u_{s}}(i)\boldsymbol{u_{s}}(j),$
(30)

where $\lambda_{s}$ are eigenvalues of $\mathbf{\Delta}$, while
$\boldsymbol{u_{s}}(i)$ and $\boldsymbol{u_{s}}(j)$ are the $i$-th and $j$-th
component of the eigenvector $\boldsymbol{u_{s}}$ corresponding to
$\lambda_{s}$. For Matérn kernels we have
$\Psi(\lambda)=\left(\frac{2\nu}{\kappa^{2}}+\lambda\right)^{-\nu}$ from which
it is clear that differentiating $k$ with respect to the smoothness parameter
$\nu$ is straightforward. This allows to optimize the marginal likelihood over
the smoothness parameter $\nu$ in the same way this optimization over the
length scale parameter $\kappa$ is usually done.