FACULTY OF COMPUTER SCIENCE, ELECTRONICS AND TELECOMMUNICATIONS
INSTITUTE OF COMPUTER SCIENCE
Master’s Thesis

Application of Graph Neural Networks and graph descriptors for graph
classification

Zastosowanie grafowych sieci neuronowych oraz deskryptorów grafowych w
klasyfikacji grafów

Author: | Jakub Adamczyk
---|---
Field of study: | Computer Science - Data Science
Supervisor: | Wojciech Czech, PhD

Krakow, 2022

## Acknowledgements

I would like to sincerely thank my supervisor, Wojciech Czech, for his
guidance. I am thankful for his suggestions, including the proposed subject of
this thesis. His hard work and vast knowledge of graph methods have been a
great source of inspiration and motivation for my work.

I also want to thank AGH University of Science and Technology, and ACK
Cyfronet in particular, for computational resources, which aided my research.

In particular, I would like to thank my grandparents and parents for teaching
me discipline and diligence, owing to which I could finish my thesis.

## Abstract

Graph classification is an important area in both modern research and
industry. Multiple applications, especially in chemistry and novel drug
discovery, encourage rapid development of machine learning models in this
area. To keep up with the pace of new research, proper experimental design,
fair evaluation, and independent benchmarks are essential. Design of strong
baselines is an indispensable element of such works.

In this thesis, we explore multiple approaches to graph classification. We
focus on Graph Neural Networks (GNNs), which emerged as a de facto standard
deep learning technique for graph representation learning. Classical
approaches, such as graph descriptors and molecular fingerprints, are also
addressed. We design fair evaluation experimental protocol and choose proper
datasets collection. This allows us to perform numerous experiments and
rigorously analyze modern approaches. We arrive to many conclusions, which
shed new light on performance and quality of novel algorithms.

We investigate application of Jumping Knowledge GNN architecture to graph
classification, which proves to be an efficient tool for improving base graph
neural network architectures. Multiple improvements to baseline models are
also proposed and experimentally verified, which constitutes an important
contribution to the field of fair model comparison.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Research goals
    3. 1.3 Contributions
    4. 1.4 Thesis overview
    5. 1.5 Notation and definitions
  2. 2 Literature review
    1. 2.1 Graph neural networks
      1. 2.1.1 Introductory works and general architecture
      2. 2.1.2 GNN taxonomy
      3. 2.1.3 Convolutional architectures
      4. 2.1.4 Attentional architectures
      5. 2.1.5 Readout
    2. 2.2 Graph descriptors
    3. 2.3 Evaluation methodology
  3. 3 Theoretical background
    1. 3.1 Graph neural networks
      1. 3.1.1 Graph convolutional layers
      2. 3.1.2 Readout layer
      3. 3.1.3 Skip connections
      4. 3.1.4 Baseline architectures
    2. 3.2 Graph descriptors
    3. 3.3 Evaluation methodology
  4. 4 Experiments, results, and discussion
    1. 4.1 Experimental setup
      1. 4.1.1 Datasets
      2. 4.1.2 Data preprocessing
      3. 4.1.3 Architectures
      4. 4.1.4 Model training, selection, and evaluation
    2. 4.2 The results and discussion
      1. 4.2.1 Overall models performance comparison
      2. 4.2.2 Baselines performance
      3. 4.2.3 GNN models performance
      4. 4.2.4 LDP models performance
      5. 4.2.5 Fingerprint models performance
      6. 4.2.6 Learning curves analysis
      7. 4.2.7 Hyperparameters analysis
  5. 5 Thesis summary
    1. 5.1 Key findings and conclusions
    2. 5.2 Future work
  6. A Full results tables
    1. A.1 HIV dataset
    2. A.2 BACE dataset
    3. A.3 BBBP dataset
    4. A.4 Lipophilicity dataset

##  1 Introduction

Many practical applications require solving problems represented as tasks on
graphs. They arise in many areas, such as theoretical computer science,
bioinformatics, linguistics, and social sciences. Graph structure is a very
abstract and general representation and, as such, can serve as a framework to
formulate many distinct problems in a unified way. Graph theory, pioneered by
Euler in the 18th century, has very strong mathematical foundations. A
multitude of practical applications, combined with a strong theoretical
background, resulted in the dynamic development of this branch of mathematics
with many sub-disciplines. One of them, particularly interesting from a
computational perspective, is graph comparison and matching, which aims to
infer the properties of a graph and use them to answer questions such as:

  * •

How can we define the similarity between graphs?

  * •

How can we represent graphs in metric spaces?

  * •

How can we use graph representations combined with their similarity for
predicting the properties of other graphs?

Tasks of this discipline can also be formulated in terms of machine learning,
which has long traditions in computer science and computational learning
theory. This resulted in the creation of graph representation learning, a
field of study which aims to learn graph properties from the gathered data and
use them for supervised or unsupervised learning tasks on graphs [53]. The
tools used include graph descriptors combined with classification algorithms,
graph kernels directly defining graph similarity, and graph neural networks.

###  1.1 Motivation

Most real-world graph datasets are far too large for manual analysis, and
instead require scalable algorithms. They can be divided into three groups
[53, 113, 123]: node-level, edge-level, and graph-level. We can perform both
supervised and unsupervised tasks on each level, e.g. node classification,
node clustering, edge regression, graph classification. Graph-level tasks are
arguably the most well studied in mathematics, as they use many graph features
very well studied theoretically, e.g. algebraic connectivity [107], efficiency
[69] or chemical indices [109, 50, 91]. Graph kernel methods are also used at
this level, as they describe whole graph structures to compare them to each
other. There are many practical applications of this type of graph task in
both supervised (when we know target labels) and unsupervised (we do not)
settings:

  1. 1.

Graph classification: supervised learning, when the target is a class from a
finite, unordered, discrete set. Applications include, e.g. predicting
mutagenicity [64], classifying molecules with activity against HIV, [63] or
detecting malicious software based on analysis of its control flow graph [77].

  2. 2.

Graph regression: supervised learning, when a target label is a continuous
number. Applications are, e.g. inferring solubility of molecules [86] or
predicting quantum properties of compounds, such as atomization energy [47].

  3. 3.

Graph clustering: unsupervised learning, when we aim to group graphs into a
given number of clusters. The example application is clustering substrate-
reaction metabolic networks for constructing phylogenetic trees [32].

Out of those three groups, graph classification and regression are more
popular. It is because a plethora of problems from bioinformatics and related
domains (pharmacology, molecular biology, biochemistry, quantum physics, and
chemistry) are naturally expressed as supervised learning on molecular graphs.
Those are often binary classification (with two classes, whether a molecule
has or does not have a given property) or regression (magnitude of physical
property, e.g. solubility). There are multiple datasets available in this
area, and the interest in the scientific community as well as in the private
sector is growing. In silico experiments using machine learning-based
predictions are becoming an important tool, e.g. in novel drug discovery [51].

Most works on machine learning on graphs focus on undirected, unweighted,
homogenous simple graphs, possibly with attributed nodes (vertices), and do
not consider edge attributes, as this type of graphs is the most common in
general. Directed graphs have asymmetric adjacency matrices, which require
dedicated GNN architectures to handle them, as it results in undesirable
numerical properties, such as potentially complex eigenvalues. Heterogeneous
graphs (with nodes with different sets of attributes) can be represented as
homogenous graphs with high-dimensional node feature vectors, or we can work
on them directly with dedicated GNN architectures. They are much more common
in the context of recommendation systems using graph-based approach, where
users and items are distinct types of nodes. Edge attribute information can be
incorporated with Levi graph transformation [49, 71] and regular GNN
architectures can be applied then, or we can use dedicated architectures.
Multigraphs (graphs with more than 1 edge between any pair of nodes) are much
less frequent than simple graphs and not applicable in many areas (e.g. in
molecular structures). However, many of the ideas in this thesis apply to
graphs of any type.

Traditional machine learning methods, including deep learning, have been used
for structured data, such as tabular datasets, audio, or images. This data has
a well-defined metric space, which makes them easier to analyze with both
manual feature engineering and deep learning. Graphs are combinatorial
structures and, as such, require dedicated algorithms and feature engineering.
However, all approaches work under the hood in the same way: they extract a
graph embedding, i.e. a vector of typically continuous numbers, which
represents a graph as a point in some high-dimensional metric space. Given
this representation, the problem is just a tabular classification. Therefore,
the most common approaches really differ in the way they create the graph
embedding.

The traditional approach calculates features selected by the researcher and
uses them as a direct input for the classifier. Another common name for those
features is graph invariants, as the vast majority of them depend solely on
the structure of the graph, not on any particular labeling or ordering of the
nodes or edges in the graph, or their permutation (they are isomorphism
invariant). They are most often based either on graph theory (e.g. node
centrality, clustering coefficient), or algebraic graph descriptors [33, 53])
or on domain knowledge, e.g. using molecular fingerprints in chemistry (e.g.
Extended Connectivity Fingerprints (ECFP) [93]). They can describe both the
local and global structure of the graph, but typically they require domain
knowledge and choosing the right set of descriptors is not straightforward
[104]. Some of them are learned from data, e.g. ECFP, and some are
predetermined by expert knowledge, e.g. MACCS keys [37].

Historically, graph kernels [46] were often used to capture the similarity
between graphs using their structural properties (e.g. using random walks or
subtrees [45], or isomorphism testing heuristics [97]). They are an example of
shallow learning, as they learn from data, but do not learn representations in
an end-to-end fashion, like deep learning. While they are highly expressive,
they come at a high computational cost, both in terms of their calculation
(typically at least $O(n^{3})$) and later use for similarity search (at least
$O(n^{2})$), making their usage infeasible for large datasets.

Deep learning applied to non-euclidean domains, known as geometric deep
learning [24, 23], has seen a surge in popularity in recent years. A large
part of that revolution is graph neural networks (GNNs), designed specifically
to learn node and graph embeddings in an end-to-end fashion. An important
feature of those architectures is preserving input invariance, i.e. the order
of nodes does not matter, which was not the case for traditional deep learning
models [94, 113]. The key advantages of GNNs are relatively low computational
cost (compared to graph kernels) and end-to-end training (learning from data
for efficient graph representation tailored to a given task). However, they
are based on the message-passing paradigm [94, 53, 47], in which neighboring
nodes exchange information. This means that capturing more global features and
explicit structural graph features is problematic, requiring sophisticated
architectures, deep networks, and solving problems arising from them (e.g.
exploding/vanishing gradients, overfitting, oversmoothing [75, 117]).

###  1.2 Research goals

The aim of this work is to investigate techniques for graph classification, on
tasks requiring capturing both local and global features of graphs. To this
end, we focus on modern, general-purpose GNN architectures (not limited to any
particular domain), and on graph descriptors, such as molecular fingerprints.
Those approaches extract both local and global features, exploiting structural
graph knowledge as well as available node features. The specific goals of the
research project presented in this thesis are as follows:

  * •

Investigate GNN architectures for graph classification and their building
blocks, and how the architectural choices (e.g. depth, sophistication,
information aggregation strategies) affect the quality of models.

  * •

Compare the efficiency of GNNs and graph descriptor-based models, especially
on tasks requiring extraction of hierarchical information, such as molecular
property prediction.

  * •

Test solutions on a varied collection of state-of-the-art benchmark datasets,
and compare models in a unified way, following experimental procedures
designed with recommended fair evaluation and fair comparison techniques.

###  1.3 Contributions

Our contributions can be summarized as follows:

  1. 1.

We designed and created a fair evaluation benchmarking environment for GNNs. A
collection of challenging datasets with varied properties has also been
selected.

  2. 2.

We performed numerous experiments with GNN architectures, obtaining one of the
largest sets of results in the literature. Furthermore, we analyzed those
results from different perspectives, giving many insights into GNN
architectures.

  3. 3.

In addition to GNN architectures, as a part of fair comparison, we propose 4
new baselines for graph classification:

     * •

Simple Graph Convolution (SGC), applied to this task for the first time

     * •

Local Degree Profile (LDP) with additional graph descriptors

     * •

Molecular Fingerprints (MFPs) and Deep Multisets, combined with node features
embedding

  4. 4.

We performed experiments with graph descriptor-based LDP model and with
molecular fingerprints, showing that the latter outperform GNNs in most cases.
We also proposed using a new model, based on concatenation of multiple
fingerprints, obtaining a more robust classifier.

###  1.4 Thesis overview

This thesis is organized in the following way. Firstly, in Section 1.5 we
introduce the notation and definitions from graph theory used throughout this
text.

Further, in Section 2 we present the review of the literature relevant to the
thesis. In particular, it is divided into three subsections, following the
main research goals of the work. In Section 2.1 we review GNN architectures,
their particular building blocks, and strategies for using them for graph
classification. In Section 2.2 we describe the available graph descriptors
from various domains. Finally, in Section 2.3 we introduce evaluation
strategies for ML in graph classification.

Section 3 presents the detailed theoretical background of the methods relevant
to performed experiments. In particular, Section 3.1 focuses on building
blocks of GNN architectures, which are explained and compared based on their
theoretical properties and the capability of capturing graph properties. In
Section 3.2 we describe graph descriptors, algorithms for their calculation,
their interpretation, strengths, and weaknesses. We outline the details of
evaluation methodologies, good practices of fair evaluation, and their
applications to the graph domain in Section 3.3.

In Section 4 we provide descriptions of experiments preparation, results, and
their discussion. Experimental setup, including architectures used, datasets,
hyperparameter tuning, and evaluation methodology are in Section 4.1. We
present results and discuss them in Section 4.2.

The thesis summary, drawn conclusions, possible future work, and final remarks
are presented in Section 5.

###  1.5 Notation and definitions

In this section, we recall the basic mathematical notation and definitions
used further in this work. Like most works on GNNs, we consider primarily
undirected, unweighted, homogenous simple graphs, possibly with attributed
vertices (nodes). We do not consider edge attributes for GNNs, since they
require graph transforms or dedicated architectures; however, graph
descriptors may make use of this information. We do not consider directed or
heterogeneous graphs or multigraphs, as they require dedicated GNN
architectures and the problems requiring them are more specific and less
popular. However, most of the ideas in this thesis also apply to them. Because
of this, we will assume in the rest of this thesis that the graphs possess
appropriate properties of undirected, unweighted, homogenous simple graphs,
e.g. symmetric adjacency matrix.

First, we present the mathematical notation used. The special sets of numbers
are represented with a blackboard bold font, e.g. $\mathbb{R}$ is a set of
real numbers.

Vectors are represented with lowercase Latin or Greek alphabet letters.
Matrices and tensors are written as uppercase letters. For example: $x$, $A$,
$\Theta$. Scalars at particular indices are denoted as $x_{1}$, $A_{1,2}$.
This notation is also often used to denote different vectors/matrices, e.g.
$G_{A}$ vs. $G_{B}$, however, the meaning should be always obvious from the
context in those cases. For notational convenience, commas in indices will be
sometimes omitted, e.g. $A_{ij}$. In the case of matrices, if the entire
dimension is selected, a colon is used, e.g. $A_{i,:}$ means selecting the
entire $i$-th row. Vectors are indexed with integers from 1 to $n$, and
matrices are indexed similarly and are of shape $(m,n)$ ($m$ rows and $n$
columns), unless otherwise specified. Tensors are indexed similarly to
matrices. Shapes of vectors, matrices and tensors are written in upper index
as an abbreviation, e.g. $x^{n}\Leftrightarrow x\in\mathbb{R}$,
$A^{m,n}\Leftrightarrow A\in\mathbb{R}^{m\times n}$. Series of matrices (e.g.
weight matrices in neural network layers) are marked with index in parentheses
in upper index, e.g. $W^{(i)}$, $H^{(i+1)}$. The identity matrix is denoted as
$I$. A diagonal matrix with vector $x$ on diagonal is marked as $diag(x)$.

The dot product between vectors is marked as $\cdot$, e.g. $z=x\cdot y$.
Scalar, matrix, or tensor multiplication is marked as two letters next to each
other, e.g. $C=AB$. Elementwise vector or matrix multiplication, also known as
Hadamard product, is denoted with $\odot$, e.g. $C=A\odot B$. Matrix
multiplication with itself $k$ times, also known as matrix $k$-th power, is
written as $A^{k}$, e.g. $AAA=A^{3}$. Matrix or vector concatenation is marked
$\operatorname*{\scalerel*{\|}{\sum}}$ operator, with vectors/matrices being
concatenated in squared brackets, e.g.
$c=\left[a\operatorname*{\scalerel*{\|}{\sum}}b\right]$. The concatenation
dimension for matrices should be obvious from the context, and most often is
the horizontal dimension (increasing the number of columns).

Functions are marked with lowercase Latin or Greek alphabet letters, e.g.
$f(x)$. It is assumed that arguments are vectors, unless otherwise specified.
In particular, $\sigma(x)$ denotes the sigmoid function, not just any function
in general, like in some works on GNNs:

$\sigma(x)=\frac{1}{1+e^{-x}}$

Below, we present the basic graph theory definitions and notation.

######  Definition 1.1 (Undirected graph).

An undirected graph $G$ is defined as an ordered pair $G=(V_{G},E_{G})$, where
$V_{G}$ is a set of vertices and $E_{G}$ is a set of edges. An edge
$e_{ub}=\\{u,v\\}\in E_{G}$ is an unordered pair of vertices ($u\in V_{G}$,
$v\in V_{G}$). Two vertices $u$ and $v$ are adjacent ($u\sim v$) if they are
joined by an edge.

If the context of a graph is obvious, the indices are dropped for all
definitions concerning it, e.g. $G=(V,E)$. Edges are written as $(u,v)$ and it
is assumed that the order of vertices in the pair does not matter. Vertices
are also called nodes, which is a more standard naming in GNN literature.
Graphs are interchangeably called networks. Following standard notation from
GNN literature, $|V|$ will be also written as $n$ if the context is obvious.
For example, we use $O(n^{2})$ instead of $O(|V|^{2})$, as typically the graph
is sparse, i.e. $|V|>>|E|$, and the complexity bound depends only on $|V|$,
while in computer science the main complexity variable is traditionally
denoted as $n$.

######  Definition 1.2 (Simple graph, multigraph).

A graph is simple if there is at most one edge between any two vertices $u$
and $v$. Otherwise, if there is at least one pair of vertices $u$ and $v$
connected with two or more edges, the graph is a multigraph.

######  Definition 1.3 (Attributed graph).

A graph is attributed if its vertices or edges have feature vectors
(attributes) attached to them. Elements of those vectors are typically
numerical, either real numbers or integers. Graphs with all vectors with the
same attributes are homogenous, while graphs with at least two distinct types
of feature vectors are heterogeneous.

######  Definition 1.4 (Walk).

A walk of length $k$ from vertex $u$ to vertex $v$ is a sequence of $k$ edges
connecting $u$ and $c$. The closed walk is a walk where $u=v$.

######  Definition 1.5 (Path).

A path $p_{G}(u,v)$ between vertices $u$ and $v$ is a walk without repeated
vertices.

######  Definition 1.6 (Vertices distance).

The distance $dist_{G}(u,v)$ between vertices $u$ and $v$ is a length of the
shortest path between $u$ and $v$. If a path does not exist, then
$dist_{G}(u,v)=\infty$.

######  Definition 1.7 (Vertex neighborhood).

The (direct) neighborhood $N(v)$ of vertex $v$ is a set of vertices adjacent
to $v$. $k$-th order ($k$-hop) neighborhood $N^{k}(v)$ of a vertex $v$ is a
set of vertices of distance to $v$ lesser or equal to $k$.

######  Definition 1.8 (Vertex degree).

The degree $d(v)$ of vertex $v$ is the number of edges incident to $v$, equal
to the size of $N(v)$ for simple graphs.

######  Definition 1.9 (Permutation invariance).

A function $f$ of a vector argument $x$ is permutation invariant if the value
of $f$ does not change when the order of vector elements changes.

Unordered sets, such as vertices and edges, may be represented with indexable
structures, e.g. with vectors or matrices. However, they are still permutation
invariant. Only for implementational purposes, the indexing of nodes is fixed
after loading the graph as a matrix.

######  Definition 1.10 (Graph isomorphism).

Let $G$ and $H$ be simple graphs. A graph isomorphism between $G$ and $H$ is a
bijection: $\alpha:V_{G}\rightarrow V_{H}$ such that:

$(u,v)\in E_{G}\Leftrightarrow(\alpha(u),\alpha(v))\in E_{H}$

Isomorphic graphs $G$ and $H$ are denoted as $G\simeq H$. Isomorphism
preserves the structural information of the graph represented by its edges. It
constitutes an equivalence relation on graphs that partitions the set of
graphs into equivalence classes.

######  Definition 1.11 (Adjacency matrix).

Let $G$ be a simple graph. Adjacency matrix $A_{G}$ is a matrix for which

$A_{u,v}=\begin{cases}1,&\mbox{if }{u,v}\in E_{G}\\\ 0,&\mbox{if }{u,v}\notin
E_{G}\\\ \end{cases}$

Adjacency matrix is symmetric for undirected graphs. For graphs without self-
loops, the diagonal entries are zeros. For weighted graphs, a non-zero entry
$A_{G}(u,v)$ represents the weight of the edge between $u$ and $v$. Those
matrices are typically (very) sparse.

######  Definition 1.12 (Degree matrix).

For a simple graph $G=(V,E)$ with $n$ vertices, the degree matrix
$D\in\mathbb{R}^{n\times n}$ is a diagonal matrix:

$D_{ii}=d(v_{i})$

######  Definition 1.13 (Laplace matrix).

For a simple undirected graph $G=(V,E)$ with $n$ vertices, its Laplace matrix
(Laplacian) $L\in\mathbb{R}^{n\times n}$ is defined as:

$L=D-A$

The elements are defined as:

$L_{ij}=\begin{cases}d(v_{i}),&\mbox{if }i=j\\\ -1,&\mbox{if }(v_{i},v_{j})\in
E\mbox{ and }i\neq j\\\ 0,&\mbox{otherwise}\end{cases}$

$L$ is a positive semidefinite matrix. It is (very) sparse for sparse graphs,
which is often the case, as are all other Laplacian matrices (see below).

######  Definition 1.14 (Symmetric normalized Laplacian).

The symmetric normalized Laplacian $L^{sym}\in\mathbb{R}^{n\times n}$ is
defined as:

$L^{sym}=D^{-\frac{1}{2}}LD^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$

The elements are defined as:

$L^{sym}_{i}j=\begin{cases}1,&\mbox{if }i=j\mbox{ and }d(v_{i})\neq 0\\\
-\frac{1}{\sqrt{d(v_{i})d(v_{j})}},&\mbox{if }(v_{i},v_{j})\in E\mbox{ and
}i\neq j\\\ 0,&\mbox{otherwise}\end{cases}$

$L^{sym}$ is a positive semidefinite matrix.

######  Definition 1.15 (Vertex descriptor).

A vertex descriptor is a scalar or vector attached to the given element,
calculated using its structural properties (e.g. neighbors) or attributes
vector.

######  Definition 1.16 (Graph descriptor).

A graph descriptor (graph invariant, graph feature) is a scalar or vector
attached to the graph. It is calculated using some properties of the graph
(e.g. its structure) or using a distribution of vertex (node) descriptors
(e.g. histogram of values of a particular descriptor for all vertices of the
graph).

######  Definition 1.17 (Graph embedding).

A graph embedding is a vector attached to the graph, representing it as a
point in a metric space.

##  2 Literature review

In this chapter, we present an overview of the literature relevant to the
thesis. We review works on GNN, graph descriptors, and the evaluation
methodologies for ML models.

###  2.1 Graph neural networks

This section describes GNN architectures. First, we provide a historical
perspective and works describing the general framework of GNNs. Then, the
taxonomy of different architectures used in the literature is reviewed. Later,
we present the related works on convolutional and attentional architectures,
as well as on readout layer.

####  2.1.1 Introductory works and general architecture

Applying neural networks to graph data requires the models to be able to work
with graph matrix representation (e.g. adjacency matrix). It also requires the
model to be invariant to the input permutation, so that the nodes of the graph
can be seen by the network in any order. Most GNN models aim to learn a node
embedding for each node, which (possibly after many layers, which may also be
of many types) are aggregated for the graph classification as a graph
embedding. Due to this formulation, most works on GNNs focus on node
classification, but those architectures can be easily applied to graph
classification, as it just requires one additional step of node embedding
aggregation and classifying the graph embedding, typically with a multilayer
perceptron (MLP).

The first works for applying neural networks to graphs and introducing the GNN
concept are [95, 48, 94]. Those works are a series by the same authors, with
the most mature model presented in [94], which is sometimes called vanilla GNN
in literature. This model is limited to undirected homogenous graphs, but can
use both node and edge features. It introduces the concept of learning node
embeddings through message passing, where each node exchanges its information
(messages) with its neighbors. Information from neighbors is aggregated and
used to update the node embedding, and this process is repeated multiple
times, extending the radius of information exchange for each node. First, the
authors propose first using a local transition function to learn a node
embedding, and then using a local output function to learn an output state of
the node using the input node features and its embedding. Both functions are
parametrized and are MLPs, and their parameters are learned through
backpropagation. This is also called an encoding network, as this encodes node
states and is a recurrent neural network (RNN). Therefore, this type of
network is an example of RecGNN (Recurrent Graph Neural Network) [113].
Computing the embedding of all nodes is done through an iterative scheme with
Banach’s fixed point theorem. While this model allows learning with graph
data, it has significant limitations. The representation is quite limited, as
it does not learn hierarchical representations. Moreover, the learning scheme
is inefficient, requiring both obtaining a fixed point solution and
backpropagation in each iteration.

The general framework of neural message passing for graphs presented in [94]
has been generalized in [47] and is known as a message-passing paradigm [82]
or as a Message-Passing Neural Network (MPNN) model (also called neighborhood
aggregation strategy [115]). It is used by a vast majority of GNNs. The single
layer model in this framework represents message passing between a node and
its neighborhood, therefore in $k$ layer GNN each node captures the
information of its $k$-hop neighborhood. Formally, GNNs are parametrized with
two functions reflecting how the node neighborhood function is aggregated
($a_{v}=AGGREGATE(N(v))$) and how this information is used, along with its
embedding, to obtain the new, updated embedding
($h_{v}=COMBINE(a_{v},h_{v}^{prev}$). Additionally, for graph classification
[47] a layer with $h_{G}=READOUT(H_{v})$ function is used (often called global
pooling layer [82, 113, 40] or graph-level pooling layer [115]), which
aggregates the node embeddings matrix and outputs the graph embedding vector.
Different GNN architectures can therefore be defined in terms of choice of
those functions and their properties. For example, the difference between
families of spatial convolutional and attentional architectures lies in a
different weighting of neighbors in $AGGREGATE$ function. Such layers are
known as message-passing layers [40] or computational/propagation modules
[123]. In GNN literature, layers are also often called modules, blocks,
operators or message-passing functions.

####  2.1.2 GNN taxonomy

In recent years, numerous surveys on GNNs have been conducted. Most of them
propose a particular taxonomy, since the classification of architectures into
a few distinct families makes the comparison of technical details easier.

In [24] convolutional GNNs are reviewed and divided into two groups, depending
on the approach to deriving and implementing the convolution operator: spatial
and spectral. The work [55] contains an overview of many graph representation
learning approaches, including traditional “shallow” approaches, but also
features GNNs, divided into two groups: neighborhood aggregation
(convolutional) and message-passing (vanilla GNN, MPNN).

GNNs can also be decomposed into basic building blocks, as proposed in [18].
GNN architectures are analyzed as compositions of 6 basic building blocks,
called Graph Networks (GNs), and proposes classifying the GNN architectures in
this framework.

The authors of [123] observe that GNNs can be classified on many axes: working
on directed vs. undirected graphs, homogenous vs. heterogeneous graphs, static
vs. dynamic (changing in time) graphs, and the choice of computational
modules, which is further divided into computational modules, sampling and
pooling methods. General GNN architectures are divided by message passing
architecture into two main groups: convolutional and recurrent. Convolutional
aggregation is further defined as either spectral or spatial, where spatial
convolution can be either “basic” convolution or attentional. Recurrent
operators can be convergence-based or gate-based (using a gating mechanism
known from, e.g. Long-Short Term Memory (LSTM) [60] or Gated Recurrent Unit
(GRU) [31] networks).

Alternatively, as described in the work [113], GNNs are divided into four
categories: recurrent (RecGNNs), convolutional (ConvGNNs), graph autoencoders
(GAEs), and spatial-temporal (STGNNs). Convolutional architectures are
classified as either spatial or spectral. Basic attention-based architectures
are treated as a variant of spatial convolutional architecture, but it is
noted that there are interesting attentional GNNs that do not fit this
classification. GAEs are divided by their application, as either used to learn
network embeddings or for graph generation, and are built in practice using
either recurrent or convolutional blocks.

The work [122] provides an overview of deep learning on graphs, but focuses on
GNNs specifically. It categorizes deep learning methods on graphs into five
categories: graph recurrent neural networks (Graph RNNs), graph convolutional
neural networks (GCNs), graph autoencoders (GAEs), graph reinforcement
learning (Graph RL), and graph adversarial methods. Recurrent architectures
are divided into node-level or graph-level, based on the task they were
created to perform (e.g. node classification vs. graph generation).
Convolutional layers are further specified as either spatial or spectral and
are also differentiated by the convolution filter type (e.g. first-order vs.
polynomial), readout mechanism, or computational complexity. Attention-based
modules are included in the convolutional approach, but are detailed
separately as an improvement.

While [67] focuses on exploring the relationship between GNNs and Neural-
Symbolic Computing (NSC), authors in their analysis divide GNNs into three
groups: convolutional, message-passing (vanilla GNN, MPNN), and attentional
(using attention mechanism).

In the work [23], not only a division of GNNs into groups is proposed, but the
authors also organize them into a strict hierarchy: convolutional,
attentional, and message-passing. Each group is a subset of the previous one,
however, the authors specify that using a more general group of GNNs is not
always advantageous, as simpler architectures are naturally regularized (and
therefore prevent overfitting).

The authors of [28] propose a taxonomy for graph encoding methods: shallow
embedding, graph regularization, graph auto-encoding, and neighborhood
aggregation (GNNs). GNNs are further classified as either global filtering
(spectral) or using local message-passing, which can be general frameworks,
spectrum-free methods, or spatial. Graph convolutional methods are categorized
separately as spectrum-based, spectrum-free, spatial, or attentional.
Spectrum-based methods compute the full graph spectrum and spectrum-free
methods rely on its approximation, typically using some polynomial base.

In turn, the goal of [112] is applying GNNs to graphs for natural language
processing (NLP). The authors present different GNNs as different
implementations of graph filter learning and divide them into four groups:
spectral convolutional, spatial convolutional, attentional, and recurrent.
Pooling layers are also considered and divided into global and hierarchical.

As the focus of this work is graph classification, layered feedforward
architectures are of the most interest, as they can easily learn using
supervised learning with labels for the whole graphs. While GAEs can be used
for graph classification, they do not learn end-to-end for this specific task
(since they optimize the encoding error during training, not misclassification
error), therefore their performance is expected to be worse, and for this
reason, they will not be considered. We will analyze the following
categorization of architectures, inspired by the most of the works described
above: recurrent GNNs (RecGNNs), convolutional GNNs (ConvGNNs), attentional
GNNs (AttGNNs); ConvGNNs can be either spectral or spatial. We will focus on
convolutional and attentional approaches. Recurrent architectures are mostly
used in the pioneering works on GNNs, have high computational costs, and in
general have been outperformed by convolutional and attentive architectures
for node and graph classification [113, 123], therefore they will not be
considered further. We treat the attentional architectures as a separate
group, next to the convolutional layers, as the attention mechanism is a
sophisticated weighting scheme and many recent works focus on it specifically
[102, 22, 119, 100, 66]. However, it should be stressed that attentional
architectures also use graph convolution, i.e. they simply implement
attention-based graph convolution.

####  2.1.3 Convolutional architectures

Convolutional GNN (ConvGNNs) architectures have been the most popular ones in
recent years, as they are efficient, convenient to compose, and their
similarities to CNNs can be exploited [113]. They are divided into spectral-
based or spatial-based, depending on how the graph convolution is defined, as
for graphs the definition is not as obvious as for images. Spectral approaches
stem from graph signal processing, where the convolution operation is
interpreted as smoothing the graph signal, i.e. removing the noise from it
[98]. It operates on the eigendecomposition of the Laplacian and uses graph
Fourier transform to transform the signal, convolving it in the spectral
(frequency) domain with learnable filters. Spectral approaches are sometimes
considered more suited for the transductive setting, where the test set is
known and consists of a set number of cases; this is a common case for node
classification, where a given subset of nodes in each graph has unknown
labels. Spatial ConvGNNs are more similar to CNNs, as they define graph
convolutions based on a node’s spatial relations, which can be also seen as
direct information propagation with vertex neighborhood [113]. While this
approach aims to have both lower computational cost (not requiring the
eigendecomposition) and better generalization properties (not depending on a
particular eigenbasis), defining the convolution for differently sized
neighborhoods and maintaining the local invariance property is not
straightforward. They aim to perform well in inductive settings, i.e. to learn
general rules and apply them to test data having the same distribution as the
training data.

One of the first works on ConvGNNs is [25], which introduced the Spectral
Network architecture and basic ideas of spectral analysis of graphs: Laplacian
eigendecomposition, treating convolution as signal smoothing, and learnable
localized filters with multiple channels. A spatial interpretation is also
provided. This pioneering work has four major limitations [113, 82], all
related to computing the exact eigendecomposition of graphs. Firstly, any
change to the graph changes the eigenbasis, which heavily affects stability
and therefore the generalization properties. This can also be interpreted as
inherent overfitting. The learned filters are domain-dependent, so transfer
learning similar to CNNs is not possible. Additionally, the filters may be
non-spatially localized. Lastly, the eigendecomposition itself is
computationally costly, having $O(n^{3})$ complexity.

Further works on spectral ConvGNNs focused on simplifying the ideas from [25].
The approximations both reduce the computational complexity and regularize the
model. The work [56] introduced the general idea of approximating
convolutional filters with truncated expansion using $K$-th order Chebyshev
polynomial basis. ChebNet [35] is an implementation of this idea for ConvGNNs.
This convolution has the advantage of being $K$-localized, since it is a
$K$-th order polynomial in the Laplacian, therefore the filters can extract
local features independently of graph size [113]. It also does not require
computing the eigenvectors of the Laplacian, only the largest eigenvalue,
which greatly reduces the computational load. Additionally, pooling is used,
employing the Graclus algorithm [36] for graph coarsening and node ordering
algorithm to transform the graph to a 1D signal. Another parametrization, with
Cayley polynomials (which are parametric rational complex functions), which
aims to capture narrow frequency bands in the spectral domain, is presented in
[72]. An even stronger simplification is presented in well-known work [65] as
Graph Convolutional Network (GCN). This architecture limits the layer-wise
convolution to just $K=1$, since previous architectures had the tendency to
overfit on local neighborhood structures for graphs with wide degree
distributions [82]. For this reason, this first-order approximation became a
popular technique for ConvGNNs. Additionally, the largest eigenvalue is
approximated as $\lambda_{max}=2$ and free parameters are reduced to just one
vector. Those vectors are composed into a matrix to allow multichannel input
and output. Numerical instability (vanishing/exploding gradients) of deep
multilayer GNNs is also addressed through the introduction of “renormalization
trick”, in which self-loops are added to the graph (i.e. identity matrix is
added to adjacency matrix). The resulting model is well regularized,
localized, and greatly simplified. It can be also treated as a spatial method
[82, 113], which in each layer aggregates the information from the direct node
neighborhood.

Further works on spectral ConvGNNs mostly focused on expanding the structure
learned in the convolution layer. While all previous models use the original
graph structure to propagate information, Adaptive GCN (AGCN) [76]
additionally learns the underlying implicit relations that may exist in the
graph. Additional “residual” edges added this way shorten the paths in the
graph between nodes important to each other, increasing the flow of
information between them. To this end, AGCN computes “residual” the Laplacian
matrix using the learned adjacency matrix, which is in turn computed via a
learned metric. Mahalanobis distance and Gaussian kernel are used to compare
nodes and, as a very generic framework, create a large space of learnable
solutions. However, this results in a dense adjacency matrix, which greatly
increases the computational cost of processing all layers in the graph.

Spatial approaches define convolutions directly using the topological
structure of the graph, and have been largely developed in parallel to
recurrent and convolutional spectral approaches. Neural Network for Graphs
(NN4G) [83], proposed in parallel with vanilla GNN, learns graph dependencies
through a compositional layered architecture, where the convolution directly
sums up a node’s neighborhood information. Residual and skip connections are
used to use information from different layers. The downside of this early
approach is the usage of an unnormalized adjacency matrix (which may cause
stability problems). Neural Fingerprints (Neural FPs) [38] architecture takes
a similar approach, as it also sums up the embeddings of the node and its
entire neighborhood. However, it also multiplies that sum with a weight
matrix, where each node degree has a separate matrix. Unfortunately, this
parametrization makes it infeasible to apply this model for large-scale graphs
with many node degrees.

Diffusion Convolutional Neural Network (DCNN) [14] interprets graph
convolution as a diffusion process, where information is transferred between
nodes with a certain transition probability, so that information distribution
can reach an equilibrium. It is efficiently implemented as a probability
tensor power series, where the probability transition matrix is computed using
degree and adjacency matrices, so the graph topological structure drives the
probability distribution. Hidden representation (embedding) matrices from each
diffusion step are concatenated as the final model outputs. For graph
classification, the average across nodes is taken, and the averaged vectors
are concatenated instead. Diffusion Graph Convolution (DGC) [79] takes a
similar approach, but sums up the outputs of each diffusion step instead of
concatenating them. In both of those models, the power of a transition
probability matrix is used, therefore very little information is contributed
by distant neighbors. While this may be beneficial for applications where
direct neighbors are exponentially more important (e.g. social networks), this
means that more global structure is not considered, which in turn may be
harmful in applications where long dependency chains exist (e.g. biology).

Some spatial convolutional models use the neighborhood sampling approach,
where a subset of $k$ neighboring nodes is taken into consideration as a
convolutional patch. In those approaches node ordering is often introduced,
where nodes (either all from the graph or in each neighborhood separately) are
ranked by their importance, and only a subset of fixed size $k$ is used.

PATCHY-SAN [85] takes a 4-step approach: node selection, neighborhood
assembly, graph normalization, and convolution with traditional CNN
architecture. First, all nodes in the graph are ordered using some labeling
procedure (e.g. node degree, centrality, or Weisfeiler-Lehman color [113]),
and $w$ nodes are selected using stride $s$ (padding with zeros if required).
For each one of the selected nodes, the neighborhood is assembled using
breadth-first search, so that each node has exactly $k$ neighbors; note that
BFS node selection order is known, as the order of nodes is fixed. Graph
normalization is a mathematically involved step, in which the graph is
transformed into a regular Euclidean grid structure. The idea is that nodes
from different graphs that have similar structural roles should be assigned
similar positions. After this step, a regular CNN architecture can be used, as
the normalized neighborhoods act as receptive fields and node and edge
attributes are regarded as channels. This architecture inherently leverages
both node and edge attributes and is able to apply regular CNN architectures
to graph data. However, it is computationally expensive (optimal neighborhood
normalization is even NP-hard, so approximations have to be used [85]), the
multistep method is complicated, and the method does not really work on graph
data - it merely transforms it into an easier problem in Euclidean domain,
which may lose some important information and makes layered composition much
harder.

GraphSAGE [54] proposes two innovations: a general inductive framework for
different ConvGNNs and a neighborhood sampling method. Node neighbors are
uniformly sampled to create fixed-size sets, which keeps computational cost
constant, which is especially advantageous for graphs with varying node
degrees, e.g. social networks. However, this comes at the cost of lowered
accuracy. Neighbors are then aggregated using a chosen aggregating function,
leading to different GraphSAGE subtypes. Mean aggregation can be interpreted
as an approximation of a transductive GCN approach, or as a form of skip
connection [82]. LSTM aggregator has a larger expressive capability but is not
permutation invariant and node neighbors are permuted to simulate this. The
pooling aggregator uses a fully connected layer with max pooling of outputs.
In addition to those mechanisms, an unsupervised loss function is proposed to
encourage nearby nodes to have similar representations and distant nodes to
have different representations.

A quite special type of spatial ConvGNN is Graph Isomorphism Network (GIN)
[115], which is constructed using very different principles compared to many
other architectures. Authors show the connection between GNNs and Weisfeiler-
Lehman isomorphism test, showing that typical GNNs (e.g. GCN, GraphSAGE) are
at most as powerful as WL graph isomorphism test. They also provide examples
of graph structures that cannot be distinguished by such architectures and
propose GIN model to solve this problem, which is proved to be as powerful as
WL test. The main requirement for this is using injective $AGGREGATE$ and
$COMBINE$ functions, which are summation and 2-layer MLP, respectively. For
graph classification, the choice of readout function is also studied in this
work, and the sum function is proved to be more powerful in distinguishing
graph structure than mean and max.

####  2.1.4 Attentional architectures

Attention mechanism [15, 30, 101] can be considered an additional learned
weighting scheme, using which the more “important” input features receive
higher weights. This contrasts with plain convolutional architectures, which
assume either the identical contribution of neighbors to the central node
(e.g. GraphSAGE) or with predetermined weights (e.g. GCN) [113]. A special
type of attention is self-attention [30], most often used in GNNs, where input
attends to itself, identifying the most important parts. Attentional
architectures are inherently spatial, as they operate on a node neighborhood
to assess the importance of neighbors.

Graph Attention Network (GAT) [102] applies the self-attention mechanism (a
single layer feedforward neural network) to node neighbors as an attention
layer to weight them. Additionally, multihead attention [101] stabilizes the
learning process [82] and increases the model capacity [113]; $k$ independent
attention mechanisms are used in parallel and their outputs are concatenated
or averaged to obtain the node embedding. This method deals with the problem
of different node degrees through learned weights, so that wide node
distributions are not problematic. It is also efficient in terms of
parallelization capabilities and can be easily applied to inductive learning
problems. The downside is that weighting may sometimes result in highly
varying predictions and high standard deviations [96]. The authors of [22]
observe that applying linear layers one after another in GAT means that the
ranking of attended nodes is unconditioned on the query node. To deal with
this problem, they propose GATv2 architecture, in which every node can attend
to any other node. They show that this way the model is more resistant to
structural noise (e.g. small changes in graph structure) and requires fewer
attention heads to perform as good or better than plain GAT. Since this method
requires just a simple reordering of operations, it does not increase the
computational cost and retains the advantages of GAT.

Gated Attention Network (GaAN) [119] also uses the multihead attention
mechanism, but in contrast to GAT it uses key-value attention and dot product
attention. Moreover, while in GAT heads have uniform weights, in GaAN they are
weighted using an additional soft gate, called gated attention aggregation,
which is implemented as a convolutional network. This allows differentiating
more, and less important subspaces created by attention heads, which may be
especially different between nodes, e.g. with different structural roles in
the graph.

####  2.1.5 Readout

Readout layer aggregates the node representations flatly into a single graph
embedding vector. It is also called global pooling [40], global aggregation,
[115] or direct pooling [123].

The most basic and, at the same time the most popular, readout techniques are
mean, max, and sum [113, 40], which are applied as the last layer before MLP.
As shown in [115], the sum function is the most powerful of the three, as
there exist graph structural patterns that mean and max aggregators cannot
differentiate, but the sum can. Alternatively, a weighted average can be used,
with weights computed using attention, called global (soft) attention [78,
40].

SortPool [121] proposes sorting the node embeddings according to the
structural roles of the nodes. WL coloring is used to this end as a graph
labeling technique since it provides an ordering consistent among graphs, i.e.
if nodes have similar structural properties in different graphs, they will be
assigned similar relative positions. Such structure can be fed to
1-dimensional CNNs and max pooling layers, and their output is classified with
MLPs.

Sequence-to-sequence models have been adapted to sets [103] and are known as
Set2Set [40, 123]. They can be used for global pooling, as a graph is an
unordered set of nodes. It uses an LSTM-based method and special attention,
called content-based attention, which has the property that the vector
retrieved from memory using the attention will not change if the memory is
randomly shuffled.

###  2.2 Graph descriptors

Graph descriptors, also called graph invariants or features, have been
developed historically to capture the topological structure of the graph.
Since they are hand-crafted and based on mathematical properties or domain-
specific observations, they are highly interpretable and their advantages and
disadvantages in particular situations are known.

A distinction can be made based on which graph element is described: the whole
graph or its element - vertex, edge, or a pair of such. Graph invariants
describe the graph as a whole, typically reflecting a single topological
property. For this reason, they cannot be used to reconstruct the graph
univocally, but are also largely independent (uncorrelated), meaning that
multiple can be used together to form a good feature set for graph
classification. They directly describe the global structure of the graph.
Examples include density, diameter, radius, efficiency, graph clustering
coefficient, and multiple features derived from graph spectrum, such as
spectral radius or algebraic connectivity [17, 108, 107, 32, 49].

Vertex descriptors determine vertex importance or its topological role in the
graph. Edge descriptors are similar to vertex descriptors, but they describe
edges. Pair descriptors operate on two vertices or edges. They often describe
the similarity or dissimilarity of the given objects. For graph
classification, the values of the given element descriptors are aggregated,
similarly to the node feature aggregation in GNNs. Typically, a histogram of
values with a given number of bins (which is a hyperparameter) is calculated.
Element descriptors include:

  * •

vertices: degree, clustering coefficient, betweenness centrality, PageRank,
communicability betweenness [32, 105, 21]

  * •

edges: edge connectivity, range of edge, edge betweenness [32, 42]

  * •

pairs: shortest path length, longest path length, commute time [32, 90]

Many complex graph descriptors, called molecular fingerprints, were created in
the field of cheminformatics (area between computer science and computational
chemistry). They are domain-specific, as they have been designed precisely to
capture the physiochemical properties of molecules. Their calculation
typically requires chemical graphs, which are node-attributed with atom types
and edge-attributes with bond types. They have been used mostly in scenarios
requiring graph classification, such as virtual screening (high-throughput
screening (HTS), discarding inactive molecules in large molecular databases)
or structure-activity relationship modeling (SAR modeling, predicting
bioactivity of molecules).

From a mathematical point of view, they are often similar to bag-of-words,
where a number of topological substructures are first counted and then hashed
into the fixed-length vector (representing a dictionary/hashmap). This makes
those descriptors lie somewhere between node and graph descriptors. Arguably
the most famous descriptor is an Extended Connectivity Fingerprint (ECFP)
[93], also known as circular or Morgan fingerprint. It relies on counting atom
types in node neighborhoods with increasing radii (in a fashion similar to
Breadth-First Search (BFS)). RDKit fingerprint [9], introduced by RDKit open-
source library for cheminformatics, has been inspired by the proprietary
Daylight fingerprint. It identifies all subgraphs with a given size in a
molecule, taking into account atom and bond types, counts them, and hashes
them into fixed-size vectors. MACCS keys [37] is a unique fingerprint because
its 166 features are not learned from data, but instead were chosen by domain
experts as useful in molecular property prediction. They are mostly node
feature counts or subgraph counts.

###  2.3 Evaluation methodology

Proper model selection and assessment are crucial for selecting the classifier
that will really perform the best for future data. Additionally, a clear and
concise description of the model is required to understand the formal basis
for proposed algorithms. Ablation studies and the analysis of particular
elements are also required to assess the broader applicability of the model
and its sensitivity to domains, dataset sizes, and other variables. This makes
the task of proper model evaluation a hard and important problem in machine
learning in general.

In work [57], the descriptions of typical model validation methods are
provided. In particular, authors describe 3 groups of methods: holdout, cross-
validation (CV), and bootstrapping. In holdout, a single test set it is used
to assess the model performance and is used when the data is sufficiently
large. Cross-validation uses multiple test folds, but each sample from the
dataset is used only once as a test sample. This method is more useful for
small datasets, and it also provides means to calculate the standard deviation
of the test score. Bootstrapping methods rely on creating multiple (tens or
hundreds) test sets, drawn with replacement from the dataset, and as such,
they require sophisticated improvements to not be overly optimistic.
Typically, authors precisely differentiate test and validation sets, as
evaluation methods can be used in many configurations, still the final results
should always be reported on the test set.

The authors of [68] focus on the design of experiments and their evaluation in
machine learning. In particular, they highlight the importance of train-test
division and the importance of baselines. The final model metrics should be
compared to simple baselines, and other algorithms from relevant literature
and also take into consideration class imbalance, which is common for many
domains, including molecular property prediction.

The rapid development of GNN algorithms raised a concern about the quality of
ML practices and evaluation methodologies specifically in this domain. The
authors of [96] analyze architectures for node classification on standard
datasets and also introduce 4 new datasets. The last part is important, as new
models may be optimized toward standard benchmark datasets, and testing their
performance on previously unseen datasets from a similar domain leads to a
challenging and fair comparison. Indeed, the authors find out, e.g. that GAT
has a considerably high variance on 2 new datasets, a behavior that has not
been seen before on the standard datasets. The authors also follow good
practices from [80], where for fair comparison the model architectures and
many parameters are fixed, hyperparameter grids are provided, and the same
training, validation, and model selection procedure is applied to all compared
algorithms. For each dataset, 100 random training/validation/test splits are
used, where each test result is obtained as an average over 20 random
initializations for that split. 100 results are averaged as a final result for
the dataset. For model comparison, instead of performing paired t-tests (which
compare one model to another), a model ranking scheme is proposed to compare
every model to every other model.

Additionally, the work [39] concerns the fair evaluation of graph
classification and is, therefore, the most relevant to this thesis. The
authors compare 5 GNN architectures and a strong baseline, which turns out to
be the best model in many cases. Adding node degree as an additional node
feature is also considered, giving better results on social datasets. The
following evaluation scheme is proposed: 10-fold cross-validation for testing,
where for each test fold an average of 3 random initializations is used as a
result. Inner validation with holdout is used for hyperparameter tuning.

The works above raise a question of proper evaluation practices in GNNs and
the importance of comparing the results between works. To assess the
possibility of referring to the results from other works, we perform the
overview of evaluation methodology in the studies about graph classification
with GNNs.

Multiple works report only the performance of models on node classification
tasks and not for graph classification. It is nevertheless important to check
their methodology for a comprehensive overview. GCN [65] uses a holdout for
testing, repeated 10 times, and holdout for validation. AGCN [76] uses holdout
for testing and a 5-fold CV for validation. DCNN [14] contains results both
for node and graph classification, and for node classification, holdout is
used for both validation and testing. GraphSAGE [54] and GAT [102] use holdout
for validation and testing.

Other models provide results for graph classification and typically use the
same or similar datasets from biology and chemistry. Neural FPs [38] use
holdout for testing and validation. DCNN [14] for graph classification uses
holdout for validation and testing, and it should be noted that holdout parts
are unusually large, both consisting of 1/3 of the dataset. PATCHY-SAN [85]
has an unclear description, but as we understand it, the hyperparameters on
the same datasets were used for reporting results, meaning that no test set
was really used, and only the validation results are provided, using 10-fold
CV. GIN [115] only provides 10-fold validation results and does not use the
test set.

Some GNNs are evaluated only on non-standard tasks, which makes the direct
comparison of those approaches with other works impossible. DGCN [120]
evaluates the model on datasets from image segmentation, using holdout for
testing and validation. DGC [79] is used for spatio-temporal time series
forecasting, with a holdout for testing and validation.

Nevertheless, using sophisticated evaluation methodologies is problematic from
the implementational and computational points of view. Retraining neural
networks multiple times can be prohibitively costly, even on small datasets.
Using larger and more challenging datasets makes evaluation easier, as the
holdout method on a properly challenging test set is often enough. Due to the
rising popularity of GNNs, such benchmarks have been created, e.g. MoleculeNet
[114] or Open Graph Benchmark [63]. Standardized benchmarks also often feature
precalculated training-validation-test splits for the holdout method to
compare algorithms in a controlled environment.

##  3 Theoretical background

In this chapter, we introduce the theory behind the main elements of the
thesis. We explain the mathematical foundations of graph neural networks,
graph descriptors, and graph classification algorithms. Lastly, we describe
approaches to model evaluation.

###  3.1 Graph neural networks

The goal of a GNN is to learn a function assigning each graph a label from a
set of classes. This is typically realized as a two or three-step process,
highly resembling how classic convolutional neural networks (CNNs) work:

  1. 1.

Node embedding - each node from the graph is represented as a vector, i.e.
embedding

  2. 2.

Graph embedding (readout) - node embeddings are combined with a readout
function, creating an embedding vector for the whole graph

  3. 3.

Graph classification - a multilayer perceptron (MLP) is used to classify the
graph vector

This is a typical end-to-end deep learning process, as backpropagation is used
to learn the weights of both MLP and embedding parts. The 3-step process stems
from adapting GNNs for node classification to graph classification tasks,
since node classification works very similarly, albeit without the readout
step (only steps 1 and 3).

Figure 1: Visualization of general convolutional GNN architecture. Source:
[81]

The architecture of a typical GNN is similar to CNNs, consisting of stacked
graph convolutional layers, each passing its output to the next one, as shown
in Figure 1. At the end, the readout function aggregates the output of the
last layers to a vector.

GNNs use both graph topology and node features to learn node embeddings. Graph
structure is represented as an adjacency matrix $A$, while node features are
stacked as horizontal vectors into the matrix $X$. Convolutional layers output
transformed feature vectors, i.e. $H^{(i)}$ is the output of $i$-th layer,
where $H^{(0)}=X$. A layer is therefore a function $H=f(X,A)$.

In the rest of this section, we emphasize the most important GNN architectures
by writing them in bold.

####  3.1.1 Graph convolutional layers

First, we introduce spectral ConvGNNs, as they have a solid theoretical
foundation in graph signal processing and were pioneering works in ConvGNNs.
This section closely follows multiple related works [25, 65, 55, 113] and the
works about respective architectures.

For degree matrix $D$ and adjacency matrix $A$, we can represent the graph
with symmetric normalized Laplacian matrix
$L=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$. As this matrix is positive
semidefinite, its eigendecomposition is $L=U\Lambda U^{T}$, where $U$ denotes
the matrix forming an orthonormal space. The graph Fourier transform of a
signal $x$ (e.g. a vector of a single graph feature, $x\in\mathbb{R}^{N}$) is
$\mathcal{F}(x)=U^{T}x=\hat{x}$, and the inverse Fourier transform is
$\mathcal{F}^{-1}(\hat{x})=U\hat{x}$.

Spectral graph convolution is then defined for convolutional filter
$g_{\theta}=\text{diag}(\theta)$ as:

$g_{\theta}\star x=Ug_{\theta}U^{T}x.$ (1)

The choice of $g_{\theta}$ can vary. Most naturally it can be parametrized
with eigenvalues of the Laplacian $\Lambda$, so we have $g_{\theta}(\Lambda)$.
Using the exact definition has two considerable downsides. First, computing
the eigendecomposition has $O(n^{3})$ complexity, which is prohibitively high
in practice. Additionally, any change to the graph changes the eigenbasis, and
it is heavily domain-dependent, meaning that it will not generalize well to
graphs with even moderately different structures than seen during training;
this is similar to overfitting.

In the work [56] it is shown that $g_{\theta}(\Lambda)$ can be approximated by
truncated series of Chebyshev polynomials $T_{k}(x)$ up to $K$-th order. This
results in convolution, which is used in ChebNet [35]:

$\begin{split}&g_{\theta}\star
x\approx\sum_{k=0}^{K}\theta_{l}T_{k}(\tilde{L})x\\\
&\tilde{L}=\frac{2}{\lambda_{max}}L-I,\end{split}$ (2)

where $\lambda_{max}$ is the largest eigenvalue. The convolutional filter
$\theta$ is now a vector of Chebyshev coefficients, which are recursively
defined as $T_{0}(x)=1$, $T_{1}(x)=x$, $T_{k}(x)=2xT_{k-1}(x)-T_{k-2}(x)$.
This formulation is $K$-localized, i.e. for each node it depends only on its
$K$-th order neighborhood. This is a desirable feature, as this means that
filters can extract local features independently of graph size [113]. Also,
for this reason, this method is also a precursor of spatial methods, which are
characterized by their locality.

However, the previous approach still overfits, especially on local
neighborhood structures for graphs with very wide node degree distributions,
common e.g. in social networks. It is also hard to stack layers and is not
fully localized. Graph Convolutional Network (GCN) [65], arguably the most
commonly known and used ConvGNN (or GNN in general), proposes a few
simplifications.

Firstly, we set $K=1$ and learn a function linear on the graph Laplacian
spectrum in each layer. This limits the range of possible filters, reducing
overfitting, but removes the need for explicit parametrization of
$g_{\theta}$, e.g. by the Chebyshev polynomials, and allows easy stacking of
multiple layers. Furthermore, we approximate $\lambda_{max}\approx 2$, as we
expect other parameters to adapt to this during training. Setting $K=1$ leaves
us with 2 filter parameters $\theta_{0}$ and $\theta_{1}$, and to reduce
overfitting even further, we restrain the number of parameters by setting
$\theta=\theta_{0}=-\theta_{1}$. After introducing all those simplifications,
we have:

$g_{\theta}\star
x\approx\theta\left(I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}\right)x.$ (3)

Since $I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ has eigenvalues in range $[0,2]$,
stacking multiple such layers could lead to numerical instabilities and
exploding/vanishing gradients. For this reason, we add self-loops to the graph
(an edge from each vertex to itself) in the so-called renormalization trick,
replacing $I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ by
$\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}$, where
$\tilde{A}=A+I$ and $\tilde{D}_{ii}=\sum_{j}\tilde{A}_{ij}$. This graph
convolution can easily be formulated for $D$-dimensional node features (i.e.
$D$ channels):

$Z=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}X\Theta,$ (4)

where $X\in\mathbb{R}^{N\times D}$ is the graph feature matrix and
$\Theta\in\mathbb{R}^{D\times F}$ is a matrix of filter parameters. Since the
filters are learnable weights, $W$ is often used instead of $\Theta$ to mark
this matrix. The full layer-wise propagation rule is therefore defined as:

$H^{(k+1)}=\sigma\left(\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}H^{(k)}W^{(k)}\right),$
(5)

where $\sigma$ denotes an activation function, most often $ReLU(x)=\max(0,x)$.

It is also a first-order approximation of ChebNet, which is a localized
spectral filter, therefore making GCN $1$-localized, i.e. each node aggregates
messages from its direct neighborhood. For this reason, this method can also
be interpreted as a spatial convolution. This notion is further reinforced by
observing that GCN has node-wise formulation [65, 53]:

$h^{(k+1)}_{u}=\sigma\left(W^{(k)}\sum_{v\in
N(u)\cup\\{u\\}}\frac{h_{v}^{(k)}}{\sqrt{d(u)d(v)}}\right)$ (6)

where $d(v)=|N(v)|$ is the degree of node $v$. This weighted average
normalizes the influence of neighbors using their degrees. It can be
interpreted as smoothing the hidden representations locally, along the edges
of the graph, encouraging similar predictions for close nodes. Additionally,
as noted in [65], this can be loosely interpreted as a differentiable and
parametrized generalization of the Weisfeiler-Lehman algorithm with $k=1$.

In [75] it is shown that the convolution formulation of GCN is basically a
form of Laplacian smoothing. This also explains why a small number of layers
(e.g. 1-3 in [65]) often works best for ConvGNNs - having too many layers
causes oversmoothing, a negative phenomenon when embeddings of neighboring
nodes become indistinguishable.

One of the latest developments of the line of reasoning presented above is
Simplified Graph Convolutional network (SGC) [110]. It is based on the
hypothesis that nonlinearity $\sigma$ between layers is not critical, but
rather the results of GCN stem simply from local averaging. Denoting
$S=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}$, in GCN we have:

$H^{(k+1)}=\sigma\left(\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}H^{(k)}W^{(k)}\right)=\sigma\left(SH^{(k)}W^{(k)}\right).$
(7)

In SGC, on the other hand, we remove the $\sigma$ function, and the expression
for $K$-layer GNN becomes:

$H^{(K)}=S^{K}XW,$ (8)

where matrix $S$ is raised to the $K$-th power. This model has the receptive
field of $K$ hops, same as similar GCN, but has only one weight matrix instead
of $K$, which makes it extremely fast. The authors of SGC also show that SGC
corresponds to a fixed low-pass filter on the graph spectral domain. While SGC
is aimed at node classification, it can be used for graph classification as
well as a simple baseline model. For a visual comparison between GCN and SGC,
see Figure 2.

Figure 2: The visualization of GCN and SGC for node classification. Source:
[110]

The majority of modern ConvGNNs are spatial approaches, since their spatial
locality typically corresponds to lower computational cost and better
predictive performance due to capturing the properties of node neighborhoods.
Two arguably most well-known architectures in this group are GraphSAGE and
Graph Isomorphism Network (GIN).

GraphSAGE [54] proposes two changes compared to GCN: neighborhood sampling and
parameterizable aggregation function. The first technique enforces using the
same, constant neighborhood size for all nodes, no matter their degree, which
makes the computational cost predictable and constant, which is not the case
for other GNNs. This can be used for any GNN to reduce the computational cost,
but is reported [54] to reduce the performance quality as a consequence. The
propagation scheme in GraphSAGE is highly parameterizable and explicitly
models the $AGGREGATE$ and $COMBINE$ functions. Each node first aggregates the
representation of neighboring nodes, but not including itself, as
$h_{N(v)}^{k}=AGGREGATE(N(v))$, where the aggregating function can be any
function invariant to input permutation (i.e. symmetric) and with sufficiently
high representational capacity, for example:

  1. 1.

Mean - we take the average value of each feature (or hidden dimension), which
is similar to GCN. However, in GraphSAGE we do not include the node itself in
this averaging. While not presented in the paper, it is possible to use
another permutation invariant function, such as sum or max, which has been
implemented, e.g. as a part of the PyTorch Geometric [40] library.

$h_{N(v)}^{(k+1)}=\text{mean}\left(\\{h_{u}^{(k)},\forall u\in N(v)\\}\right)$
(9)

  2. 2.

LSTM - it is a trainable aggregator, representing a more complex function,
with high expressive capability. Since LSTMs process input sequentially, the
neighbors are randomly shuffled before aggregation, to simulate permutation
invariance.

$h_{N(v)}^{(k+1)}=\text{LSTM}\left(\text{permute}\left(\\{h_{u}^{(k)},\forall
u\in N(v)\\}\right)\right)$ (10)

  3. 3.

Pooling - neighboring vectors are first transformed by MLP (typically a
1-layer network with ReLU activation), and then pooled with max function. This
approach aims to capture the different aspects of neighboring nodes and
implicitly weight them [116].

$h_{N(v)}^{(k+1)}=\max\left(\\{\text{MLP}(h_{u}^{(k)}),\forall u\in
N(v)\\}\right)$ (11)

The $COMBINE$ function concatenates the node vector with its aggregated
neighbors vector and feeds the result to a single layer MLP, creating the
output hidden representation for that node:

$h_{v}^{(k+1)}=\sigma\left(W^{(k)}\left[h_{v}^{(k)}\operatorname*{\scalerel*{\|}{\sum}}h_{N(v)}^{(k)}\right]\right),$
(12)

where $\sigma$ is almost always the ReLU function. The most important
difference between GCN and the mean variant of GraphSAGE (the most similar
one) is the unique treatment of the node itself in the latter - in GCN it is
averaged just like its neighborhood. In GraphSAGE this can be interpreted as a
type of skip (residual) connection, enabling putting more weight on the node
itself, compared to its neighborhood. Formally, the difference is:

$\begin{split}&\text{GCN:}\quad
h_{v}^{(k+1)}=\sigma\left(W^{(k)}\cdot\text{mean}\left(\\{h_{u}^{(k)},\forall
u\in N(v)\cup\\{v\\}\\}\right)\right)\\\ &\text{GraphSAGE:}\quad
h_{v}^{(k+1)}=\sigma\left(W^{(k)}\left[h_{v}^{(k)}\operatorname*{\scalerel*{\|}{\sum}}\text{mean}\left(\\{h_{u}^{(k)},\forall
u\in N(v)\\}\right)\right]\right).\end{split}$ (13)

Graph Isomorphism Network (GIN) [115] formally explores the connection between
GNNs and Weisfeiler-Lehman (WL) graph isomorphism test [97]. Authors show that
generally, GNNs are at most as powerful as WL test in distinguishing graph
structures and show examples of structures that previous architectures such as
GCN or GraphSAGE are unable to distinguish. However, those GNNs have
considerable advantages over WL test: it focuses only on graph topology, while
GNNs make use of graph features, and additionally, they learn not only to
distinguish graphs structures but to embed them, so that the distance between
similar graphs in the embedding space will be small. The GIN architecture aims
to be as powerful as WL test at discriminating between graph structures, and
at the same time provide advantages of GNNs, in particular spatial ConvGNNs.
Authors prove that this can be achieved in a very simple way: neighborhood
aggregation and readout functions need to be injective. As MLPs are universal
function approximators [62, 61] and can therefore approximate injective
functions, they can be used here. The authors show that single-layer MLP has
insufficient capabilities in this regard and that 2-layer MLP is required.
Neighborhood vectors are summed before transforming them with MLP. The
equation for update becomes:

$h_{v}^{(k+1)}=\sigma\left(\text{MLP}^{(k)}\left(\left(1+\varepsilon^{(k)}\right)\cdot
h_{v}^{(k)}+\sum_{u\in N(v)}h_{u}^{(k)}\right)\right).$ (14)

$\varepsilon$ can be either a learnable parameter, a fixed constant, or set to
0, each option less powerful than the previous one. However, the authors show
that setting $\varepsilon=0$, while making the network less powerful, makes it
generalize better. With this simplification, the update can be efficiently
represented in matrix form as:

$H_{v}^{(k)}=\sigma\left(\text{MLP}^{(k)}((A+I)\cdot X)\right).$ (15)

The attention mechanism [15, 30, 101] is basically a learnable weighting
scheme, allowing nodes to “attend” to some neighbors, which are deemed more
interesting and informative, more than others. Those networks are inherently
spatial, as nodes attend to their neighbors. This can be thought of as a
generalization of static weight in GCN $\frac{1}{\sqrt{d(u)d(v)}}$, which acts
as a normalizing constant there; in attentional GNNs this is learned, rather
than set explicitly.

Arguably the most well-known attention-based GNN is Graph Attention Network
(GAT) [102]. In this model, every two neighboring nodes attend to each other,
using the arbitrary attention function $attention(u,v)$. Additionally, the
node attends to itself (similar to a self-loop). To obtain sufficient
expressive power, features are linearly transformed into higher-level space
before computing attention scores, and attention scores are normalized with
the softmax function afterward to form a probability distribution for each
node. For clarity, we omit the $(k)$ mark below, but all relevant equations
concern a single layer, e.g. $W$ means $W^{(k)}$.

$\begin{split}&\tilde{h}_{u}=Wh_{u}\quad\quad\tilde{h}_{v}=Wh_{v}\\\
&e_{uv}=attention\left(\tilde{h}_{u},\tilde{h}_{v}\right)\\\
&\alpha_{uv}=\text{softmax}_{v}(e_{uv})=\frac{\exp(e_{uv})}{\sum_{v\in
N(u)\cup\\{u\\}}\exp(e_{uv})}\end{split}$ (16)

GAT model is agnostic to the choice of the attention mechanism. The authors of
[102] chose Bahdanau attention [15], also known as additive score function, in
the form of:

$attention_{add}(\tilde{h}_{u},\tilde{h}_{v})=\text{LeakyReLU}\left(a^{T}\left[\tilde{h}_{u}\operatorname*{\scalerel*{\|}{\sum}}\tilde{h}_{v}\right]\right),$
(17)

where $a$ is a vector of attention weights (for the given layer) and LeakyReLU
is a smooth approximation of ReLU, with negative input slope $\alpha$ (in
[102] $\alpha=0.2$ is used).

$\text{LeakyReLU}=\begin{cases}x&x\geq 0\\\ \alpha*x&x<0\end{cases}$ (18)

Normalized attention scores are applied to linearly transformed node features,
and the final update equation for the node in GAT becomes:

$\begin{split}&h_{u}^{(k+1)}=\sigma\left(\sum_{v\in
N(u)\cup\\{u\\}}\alpha_{uv}W^{(k)}h_{v}^{(k)}\right)\\\
&\alpha_{uv}=\frac{\exp\left(\text{LeakyReLU}\left(a^{T}\left[W^{(k)}h_{u}^{(k)}\operatorname*{\scalerel*{\|}{\sum}}W^{(k)}h_{v}^{(k)}\right]\right)\right)}{\sum_{v\in
N(u)\cup\\{u\\}}\exp\left(\text{LeakyReLU}\left(a^{T}\left[W^{(k)}h_{u}^{(k)}\operatorname*{\scalerel*{\|}{\sum}}W^{(k)}h_{v}^{(k)}\right]\right)\right)}\end{split}$
(19)

To stabilize the learning process, multihead attention [101] is used in GAT.
Multiple attention mechanisms, called heads in this context, are used in
parallel and can be either concatenated or averaged. In GAT concatenation is
used in all layers except the last one, where averaging is used, in order to
keep the proper dimensionality (since otherwise, $T$ attention heads with
concatenation result in $T\cdot D$ output dimensionality). This does not
necessarily result in increased computational cost, since when using multiple
attention heads the output dimensionality stays the same, while accuracy is
often higher [101].

In Equation 20, following [102], we mark concatenation of series of vectors
with $\operatorname*{\scalerel*{\|}{\sum}}$ operator.

$\begin{split}&h_{u}^{(k+1)}=\operatorname*{\scalerel*{\|}{\sum}}_{t=1}^{T}\sigma\left(\sum_{v\in
N(u)\cup\\{u\\}}\alpha_{uv}W^{(k)}h_{v}^{(k)}\right)\\\
&h_{u}^{(K)}=\sigma\left(\frac{1}{T}\sum_{t=1}^{T}\sum_{v\in
N(u)\cup\\{u\\}}\alpha_{uv}W^{(k)}h_{v}^{(k)}\right)\end{split}$ (20)

Figure 3: Left: GAT attention mechanism. Right: illustration of multihead
attention for $K=3$. Source: [102]

####  3.1.2 Readout layer

Convolutional layers exchange information between neighboring nodes. While
this propagates information in the graph (using its structure), it still
requires gathering global information for graph classification. This is
achieved through $READOUT$ function, gathering information from all node
embeddings and flattening them into a single graph embedding vector, which can
be then classified. This operation is also sometimes called global pooling
[40].

Most often, simple and permutation invariant functions are used: mean, max and
sum. The selections found in the literature often seem arbitrary, however, the
authors of GIN [115] analyze this and prove that sum is the most powerful
function of the three, being able to distinguish graph structures that are
indistinguishable for mean and max functions. Some works suggest that using
more than one function is beneficial, as their results can be concatenated,
and more information is extracted, e.g. [27, 70] use mean and max, while [44]
uses all three.

More advanced readout schemes also exist, which aim to aggregate the
information from all nodes in a more sophisticated way, weighting their
influence. The examples are SortPool [121], which sorts the nodes and
transforms the output using 1D CNN (due to the loss of injectivity because of
the sorting), global attention [78] using the attention mechanism for
weighting nodes, and Set2Set (S2S) [103] based on the sequence-to-sequence
model for learning on sets. Those techniques are more computationally involved
and are outside the scope of this thesis.

####  3.1.3 Skip connections

Using skip connections (also called shortcut connections), i.e. passing the
input, or output of previous layers, directly to further layers, has been a
useful tool in CNNs since the introduction of residual connections in ResNet
[59]. In the context of GNNs skip connections are especially useful for tasks
that show homophily, i.e. where each node is strongly related to the features
of its local neighborhood [53]. They are also used to prevent oversmoothing
and vanishing gradients [82].

Residual connection is the simplest form of skip connection, where layer
implements $f(x)+x$ instead of plain $f(x)$. The authors of [59] hypothesize
that training networks with access to identity mapping makes it easier to
train the network. This notion has been formalized and explained in [74],
where it is shown that using residual connections greatly simplifies the
surface of the loss function and therefore greatly reduces the number of local
minima, making training easier. In GNNs this approach has been used e.g. in
GraphSAGE [54] and ResGCN [73].

However, many layers may contain different information useful for final
classification. For example, early layers contain more “raw” features, since
they cover a small graph radius, using just the node and its direct
neighborhood. On the other hand, layer layers are more smoothed and cover
larger subgraphs, with larger radii. We may want to capture structure
information at many layers of smoothing. Simple residual connections cannot
solve this problem, since even when a particular layer uses a skip connection
and sends an unchanged $x$ signal further, all subsequent layers must use this
output, i.e. any later representation cannot decide to use skipped or non-
skipped outputs from previous layers.

The problems outlined above have been studied in [116], where Jumping
Knowledge (JK) model is proposed. Authors show that GNNs have problems with
varying neighborhood sizes and differentiating nodes with varying structural
roles in the graph. They propose a more adaptive scheme, where networks can
learn representations of different order for particular graph substructures.
This is done through a type of skip connection called jump connection, where
the output of each layer is passed directly to the end, just before the
readout. This way we have access to the rich information about different graph
substructures, and can also adapt the effective size of the neighborhood for
each node. Three schemes of aggregation are proposed:

  1. 1.

Max pooling - for each node, we perform global max pooling, i.e. we gather its
representations at different layers and select the maximal value for each
feature among all layers. This should maximize the usage of both local
features from early layers and more globally aggregated information from later
layers. This scheme is node adaptive and has a very low computational cost.
Note that other global pooling functions can also be used in a very similar
way.

  2. 2.

Concatenation - a very straightforward, but not injective, way of combining
vectors is concatenating them. This approach is not node adaptive, i.e. it
does not differentiate between individual node structural roles well. However,
it optimizes for combining the subgraph features, which may work better for
small and more regular graphs. It has higher computational cost than max
variant, since later MLP network has to work with larger input.

  3. 3.

LSTM pooling - a bidirectional LSTM with an attention mechanism is used to
learn attention scores, and features at different levels are combined using
those weights. It is a node adaptive method, as attention scores are different
for each node. This approach is particularly effective for large, complex
graphs, but is susceptible to overfitting. The computational and memory costs
are high, since a complex bidirectional RNN with backpropagation through time
has to be computed.

Figure 4: The visualization of Jumping Knowledge architecture. “N.A.” stands
for neighborhood aggregation, i.e. graph convolution. For graph
classification, readout has to be applied after the last layer shown. Source:
[116]

Since jumping connections excel at effectively combining representations at
different scales, they are of great importance for graph classification, since
it relies on extracting hierarchical information about the graph
substructures. It allows merging of both more local and global information for
each node.

GIN [115] is based on differentiating graph substructures, and for this
reason, should greatly benefit from jumping connections. Summation is used as
a readout for each layer, and JK concatenation scheme is used.

####  3.1.4 Baseline architectures

Having proper baselines is very important for model assessment. Model
performance does not exist in a void, it should be always compared to both
simple and strong baselines. This way we can understand how much improvement
we actually obtained using new methods.

Baselines for graph classification have been created over the years. They are
based on the observation that graph data is very unique, consisting of both
node features and graph structural information, and therefore requires
dedicated algorithms. While they are simple, they often achieve results
comparable with GNNs, or even outperform them. This signifies the importance
of datasets sufficiently large for GNNs, proper architectures able to extract
task-relevant features, and hyperparameter tuning, to be able to outperform
baselines.

In [39], two approaches have been proposed. Molecular Fingerprints (MFPs) is a
simple, feature-based baseline aimed at chemical graph classification. It is
based on molecular fingerprints, a group of methods used in chemistry to
vectorize molecules; those classical methods, however, use hand-crafted
features, while this baseline is learnable end-to-end since it is similar to a
very simple GNN. The algorithm applies global sum pooling and then applies a
single layer MLP with ReLU activation. This works particularly well for
molecular datasets since they typically use one-hot encoded features, where
each feature is an atom type. The sum of such features is injective, and the
sum is simply the count of particular atoms in a whole molecule graph. This
approach should not be mistaken with molecular fingerprints for feature
extraction, which are graph descriptors and may use other classifiers than
MLP. Another proposed technique, inspired by deep learning on sets, is called
Deep Multisets. It is based on learning permutation-invariant functions on
node sets and is aimed at datasets with no node features (e.g. social
datasets) or with continuous features (e.g. protein datasets). First, node
vectors are transformed with a single layer MLP, followed by global sum
pooling, and another single-layer MLP for classification. Those approaches do
not leverage graph topology at all. For this reason, if GNNs are unable to
match or outperform those baselines, it means that either the structural
information is not very important for the problem, or that the architecture is
not able to use it properly. As shown in [39], most architectures struggle to
match the results of those methods on popular, small datasets.

###  3.2 Graph descriptors

A great advantage of GNNs is using node features and smoothing neighborhood
representations. However, they do not make explicit use of graph topology,
whereas graph kernels use graph topology exclusively and are able to achieve a
good performance, sometimes surpassing GNNs [39, 26]. To explicitly describe
the graph structure, node, edge, and graph-level descriptors can be used, or
domain-specific descriptors such as molecular fingerprints. For graph
classification, node and edge descriptors are aggregated, typically as
histograms. LDP [26], described in Section 3.1.4, is a descriptor-based
method. It also shows that features can be calculated not only for nodes
themselves but also for their neighborhoods, better describing their
structural roles in the graph. The authors of LDP also mention that they
experimented with a few other node-level and edge-level descriptors, but did
not get a consistent improvement (e.g. closeness centrality or Fiedler
vector). A particularly important group of features are centrality measures,
which can describe the importance of specific nodes, e.g. the existence of
large central stems in molecules, from which the branching parts start.

Local Degree Profile (LDP) [26] has been designed as a baseline for GNNs to be
measured against. It is a graph descriptor-based method, leveraging only
simple topology-based features describing the local neighborhood of vertices.
For each node, 5 simple structural features are computed: degree, minimal
neighbor degree, maximal neighbor degree, average neighbor degree, and
standard deviation of neighbors’ degrees. Those features summarize the degree
information about the node and its neighborhood. Each feature is then
aggregated with either a histogram (number of bins is a hyperparameter) or an
empirical distribution function (EDF). This creates a graph-level summary,
which is then classified with an SVM (of course, any other algorithm can be
used here). The authors also experiment with adding all-pairs shortest paths
aggregation to the final vector, getting about 2% improvement on chemical
datasets. Open source code for LDP [4] also shows that the authors added the
sum of neighbors’ degrees as a feature, which was not described in the paper.
This method is shown to be able to closely match or outperform (sometimes
quite significantly) graph kernels and GNNs.

Clustering coefficient [106] measures neighborhood connectivity, indicating
how close is it to a complete graph. For a single node $v$ it is defined as a
proportion of the number of edges in its neighborhood to the total possible
number of edges that could exist in that neighborhood:

$Clust_{coef}(v)=\frac{2*|\\{e_{ij}\\}|}{d(v)(d(v)-1)},$ (21)

where $i,j\in N(v)$, $d(v)=|N(v)|$, $e_{ij}=(i,j)\in E_{G}$ are edges in the
neighborhood of $v$. Alternatively, we can define it as a fraction of possible
triangles through a node [32]:

$Clust_{coef}(v)=\frac{2T(v)}{|N(v)|(|N(v)|-1)},$ (22)

where $T(v)$ is the number of triangles through a node $v$.

Betweenness centrality [41] measures how much influence a node has over the
flow of information in a graph. It finds “bridge” or “central” nodes. For node
$v$ is defined as a fraction of all-pairs shortest paths that go through $v$:

$C_{B}(v)=\sum_{s,t\in V}\frac{\sigma(s,t|v)}{\sigma(s,t)},$ (23)

where $\sigma(s,t)$ is the number of shortest paths between nodes $s$ and $t$,
and $\sigma(s,t)$ is the number of those paths that pass through node $v$.

Closeness centrality [43] is a measure of node centrality. It measures the
closeness of a node to all other nodes in terms of shortest paths lengths and
is defined as:

$C_{C}=\frac{|V|-1}{\sum_{u\in V}dist(u,v)}.$ (24)

PageRank [21] is a well-known descriptor measuring the influence of nodes in
the network. Originally made for ranking websites, it assigns high influence
to nodes connected to other important nodes, measuring the “core” elements of
the graph. For undirected graphs, it is defined as:

$PR(v)=\frac{1-d}{|V|}+d\sum_{u\in N(v)}\frac{PR(v)}{d(v)},$ (25)

where $d$ is a damping factor, a constant of value typically around $0.85$.

Figure 5: General scheme of hashed binary fingerprints. Subgraphs are
identified, hashed and the corresponding bit in a bit vector dictionary is
set. Source: [16]

Extended Connectivity Fingerprints (ECFP) [93], also called circular
fingerprints or Morgan fingerprints, is a molecule descriptor calculated using
hashing graph substructures. Hydrogen atoms and bonds to hydrogens are ignored
during calculation. Firstly, each atom in a molecule gets a unique integer
identifier, using e.g. number of “heavy” (non-hydrogen) neighbors, the atomic
number, or atomic mass. Then a given number of rounds of iterative updating
happen (see Figure 7 for visualization), up to a given radius $k$, which is
typically 2 or 3 (resulting in ECFP4 and ECFP6 fingerprint, respectively,
since traditionally diameter is used instead of radius in the name). In each
round, an extended neighborhood is considered, and all unique substructures
are given distinct identifiers. This process is quite similar to the
Weisfeiler-Lehman isomorphism test since it is a variant of Morgan’s
algorithm, which was created to solve the molecule isomorphism problem
(similar to the graph isomorphism problem of WL-test). Finally, all
identifiers are hashed into a fixed-length vector. There are two variants of
ECFP fingerprints: bit vector (binary), where only the existence or absence of
a feature at a given position is saved, and count vector, where all
occurrences are accounted for. Note that due to hash collisions this mapping
is not unequivocal and reversible, but this is not a problem in practice for
graph classification. This algorithm is able to detect large and complex
subgraph structures, especially in ECFP6 variant (with radius 3).

Figure 6: Visualization of ECFP fingerprint. Each iteration increases the
explored subgraph radius by 1, up to a given maximum. Source: [93]

Figure 7: Features discovered by ECFP algorithm for butyramide molecule.
Source: [93]

RDKit fingerprints [9, 10, 11] are inspired by publicly available descriptions
of a proprietary Daylight fingerprint. They are hashed substructure
fingerprints, similar to Morgan fingerprints - first a set of subgraph
structures is identified, then they are hashed into a fixed-length vector. It
can be a bit or count vector, similarly to ECFP. First, all subgraphs in a
molecule are identified, enumerating them as paths up to a given length,
typically between 5 and 7. Each subgraph is hashed and a relevant bit in a
dictionary is set (or count is increased). Individual subgraphs are hashed
based on their bonds, using atomic numbers, atom degrees, bond types, and atom
and bond aromaticity. In comparison to atom-centric ECFP fingerprints, the
RDKit fingerprints generate subgraphs using walks, which enables them to
discover and remember structures with more complex topology, such as path-like
subgraphs.

MACCS keys [37], also known as MDL keys [2] (after the company that developed
them), is a set of hand-crafted binary features. They consist of 166
descriptors that have been selected using expert knowledge, which are detected
in a molecule, creating a fixed-length bit vector (0 for the absence of a
feature, 1 for existence). Those are structural keys, i.e. they are not hashed
like ECFP fingerprints. Features are mostly defined using SMARTS patterns [3],
which are similar to regular expressions (regexes) but targeted towards SMILES
strings. Calculating them is available as a part of multiple open source
packages, which also list SMARTS patterns, e.g. in RDKit [5]. Apart from a few
atom counts (e.g. iodine I or sulfur S), a vast majority are groups and
hierarchical structures, e.g. carbonyl groups C=O (common in organic compounds
and typically a part of larger functional groups) or aromatic rings
(participating in many organic reactions). This means that this descriptor can
detect predetermined hierarchical structures that can be problematic to catch
for other methods, since the chances of learning topologically complex or
large substructures are low for limited datasets. For examples of some MACCS
keys and their SMARTS patterns, see Figure 8.

Figure 8: Visualization of selected MACCS keys using SMARTSviewer. Above each
diagram, the first line reports the index of the respective MACCS key and its
SMARTS pattern. Source: [29]

Using any descriptor-based method basically performs feature extraction for a
graph, generating a vector - graph embedding. While in GNNs the classification
of those embeddings is done through MLP, because of gradient flow and end-to-
end learning, in descriptor-based methods it is not required. We can use any
type of classifier, and in practice Random Forest (RF) is a very popular
choice. It typically gives very good results and is not very sensitive to
hyperparameters choice, which reduces the computational cost of hyperparameter
optimization. It also handles high-dimensional vectors generated by Morgan or
RDKit fingerprints well. It is very scalable, as all trees can be trained in
parallel, which is a considerable advantage over another popular choice,
kernel SVM, which has computational complexity between $O(^{2})$ and
$O(n^{3})$. Similarly to neural networks, class weighting scheme can be used
in RF, in two modes: globally, when class weight is computed from its
proportion in the whole datasets, or adaptively per subsample, when we compute
class weights separately in each tree based on the proportion of classes in
its bootstrap sample. The minimal number of samples to perform a split can be
tuned to prevent overfitting.

###  3.3 Evaluation methodology

Every proper evaluation of the machine learning model requires dividing the
available dataset into three parts: training, validation, and testing. We
train the model on the training set, making the labels available during
training. A validation set is used for model selection, checking the model’s
performance with given hyperparameters, and labels are not used to train the
model, but rather to steer the process of tuning and prevent overfitting to
the training data. The test set is used only once, after the entire training
and model selection, to assess the performance of the final model. The matter
of how to choose those sets is not straightforward, since while properties of
the particular dataset influence this choice, we should use the same
benchmarking process for all datasets for a fair comparison.

The two most popular approaches are the holdout method and cross-validation
(CV) [57]. In the holdout method, we simply divide the available dataset
randomly into training/validation/test parts, typically in proportions
$80\%/10\%/10\%$. This works well for large datasets, where the number of test
samples is large enough to get a proper evaluation of the model. Cross-
validation (CV) divides the dataset randomly into $K$ equal-sized, non-
overlapping parts called folds. A fold acts as a test set, and all the other
parts are merged into a training set, and this process is repeated for all
folds. This way we have $K$ test scores. A validation set is created each time
from the training set, and we can use either holdout or another CV, where
$K^{\prime}$ folds are validation sets. The latter approach is called nested
CV. Cross-validation has a high computational cost, as $K$ models have to be
trained (or $K*K^{\prime}$ for nested CV), and for this reason, it is used
only for the small datasets. It is also recommended in those cases since the
single test set would be too small and would give biased results. In addition,
using CV allows for computing the standard deviations of the test results
providing a practical measure of the model stability.

Figure 9: Visualization of nested cross-validation procedure. The source of
complexity, nested loop over folds, is clearly visible. This is 5 outer folds,
5 inner folds variant, resulting in training 25 models. Source: [1]

In the literature on graph classification, both approaches were used. The
reason for using cross-validation is that older graph classification datasets
are quite small, having from as few as 180 samples in the case of the MUTAG
[34] dataset, up to about 1000 samples, e.g. in the case of the PROTEINS [20]
dataset. However, in some works only validation scores are provided, e.g. in
GIN paper [115] (10-fold CV). Such results are always biased, giving more
optimistic results than on a proper test set. This is a major problem in
comparing the results between works, as the authors of [39] have highlighted.

Small datasets cause many issues. Having to use cross-validation, or nested
cross-validation, may increase the computational cost so much that the
training and evaluation procedure takes as long as for large datasets. At the
same time, the small sizes of single validation and test sets make the model
comparison hard. For example, if the test set consists of 100 samples (which
for many older datasets is quite normal), each sample changes accuracy by
$1\%$. For a perspective, in graph classification, the gain of $2-3\%$
accuracy is often a considerable improvement. Smaller training sets mean easy
overfitting, penalizing more sophisticated and often very interesting
architectures. They also require shorter training, which makes the initial
random initialization of weight much more important. Many of those problems
have been described in recent works [114, 39, 63]. The authors of [39] propose
a pipeline for evaluation, using static training/validation/test splits, CV
for testing, and holdout for validation. Additionally, in order to minimize
the effect of unlucky random initialization, they retrain the model 3 times
for testing, and for each test fold report the average of those three runs. We
argue that this tries to fix the symptoms, not the problem - too small
datasets. Large graph classification datasets, consisting of thousands or tens
of thousands of graphs, are currently available, e.g. many datasets from the
TUDatasets collection [84], molecular benchmark MoleculeNet [114], and Open
Graph Benchmark (OGB) [63]. Those datasets often come with pre-selected
validation and test sets, offering a unified, fair comparison between models.
They typically use the holdout method, since they are large enough.

Holdout and CV are general machine learning tools for model assessment. Those
methods, using randomized division into subsets, assume the i.i.d. property,
i.e. that samples (data points) are independently and identically distributed.
Since they have been created for tabular data, which typically has this
property, it works well in general. However, this often does not work that
well for molecular property prediction, which is probably the most common task
in graph classification, as shown by the authors of MoleculeNet [114]. Purely
random division, used e.g. in MUTAG or PROTEINS datasets, may cause very
similar molecules to be in both training and validation/test set, which makes
the task easier.

In real-world applications, when the laboratory experiments are used to verify
the ground truth for new molecules, a time split is used - a dataset created
up to a given date is a training set, and newer data serves as a test set. The
validation set is typically the part of the newest data from the training set.
This also causes a phenomenon known as a distribution shift, where newer data
follows a slightly different distribution compared to the historical data.
This is unavoidable in real-world datasets, and as such should be reflected in
model evaluation. However, this approach requires time labels to be available
for all molecules, which is almost never the case. In MoleculeNet, which is
already a large and varied collection, only a single dataset has this
information available.

Instead, we can use a scaffold split [19] (specifically Bemis-Murcko
scaffolds), which takes into consideration the structure of molecules. A
molecular scaffold reduces the chemical structure of a compound to its core
components, essentially by removing all side chains and only keeping ring
systems and parts that link together ring systems. Data is divided into
equivalence classes, i.e. groups of samples, where each group has the same
scaffold. Training, validation, and test sets are created by merging those
groups so that all molecules with a particular scaffold are in only one
subset. This creates an out-of-distribution prediction problem, where the
model is forced to learn the underlying properties of molecules to achieve
good test performance because the distribution of training and test data is
sufficiently different.

Figure 10: Comparison of train-validation-test split methods. Random split
disregards the nature of the data. Scaffold split relies on chemical structure
of molecules. Stratified random sampling for regression and time split are
shown, but are not used in this thesis. Source: [114]

Graph classification datasets are often highly imbalanced, i.e. the number of
samples from different classes is very different. It is not uncommon for
binary classification to have $90\%$ of the samples from the negative class
(e.g. non-active compound) and just $10\%$ from the positive class (e.g.
compound active against HIV). Additionally, we are often interested in the
minority class. This requires the appropriate design of algorithm,
training/validation/test splitting, and metrics used for evaluation. Neural
networks of any kind naturally lend themselves to imbalanced problems, as they
minimize loss function, which can be weighted by class proportion, so that we
put more weight on classes with fewer samples. Random sampling would use the
majority class, threatening that validation or test sets could even not
contain the minority class at all. For this reason, stratified sampling is
used in those cases. It ensures that while samples are random, the class
distribution is preserved.

The model selection is typically done through loss minimization since during
training neural networks minimize it and as such it’s the most natural way.
However, on test set we are interested in measuring the actual performance of
the model, e.g. how good is it in classifying previously unseen samples.
Different metrics measure the quality of the model in different ways. In
particular, we are interested in metrics that work well in the case of
imbalanced datasets. Many such metrics can be used, but the most common is the
Area Under Receiver Operating Characteristic Curve (AUROC, ROC AUC, or AUC).
For binary classification with positive class $1$ and negative class $0$, we
denote:

  * •

$P$ \- number of positive class samples

  * •

$N$ \- number of negative class samples

  * •

$TP$ \- true positives, the number of samples from positive class with model
prediction $1$

  * •

$TN$ \- true negatives, the number of samples from positive class with model
prediction $0$

  * •

$FP$ \- false positives, the number of samples from the positive class
incorrectly classified as $0$ by the model

  * •

$FN$ \- false negatives, the number of samples from the negative class
incorrectly classified as $1$ by the model

The AUROC is defined as an area under the ROC curve, which is created by
plotting the recall, or true positive rate (TPR), against the false positive
rate (FPR), at various positive class probability thresholds. It measures the
probability that a randomly chosen positive class sample will be assigned a
higher probability (of being from the positive class) than a randomly chosen
sample from the negative class. TPR and FPR are defined as:

$\begin{split}&TPR=\frac{TP}{P}=\frac{TP}{TP+FN}\\\
&FPR=\frac{FP}{N}=\frac{FP}{FP+TN}\end{split}$ (26)

The AUROC metric has many desirable properties for measuring efficiency of
graph classification models on imbalanced datasets. It is a single number,
which makes it a convenient measure to compare models, compared to e.g. ROC or
PRC (Precision-Recall) curves. Additionally, it is sensitive to the class
imbalance that may be present in the dataset, compared to e.g. accuracy. For
example, for a dataset with 90% negative class and 10% positive class,
accuracy will be 90%, which does not account for class imbalance in any way,
while AUROC will be low if the model does not discriminate the positive class
well. At the same time, it is not biased towards models that perform well on
the minority class at the expense of the majority class [58].

Calculating the AUROC value requires outputs of positive class probability,
not only a binary 0/1 decisions from the classifiers. This is not a major
limitation, however, since most algorithms support outputting probabilities,
e.g. neural networks, Random Forests or SVMs (through Platt scaling).
Moreover, this requires the classifier to not only output good class values,
but also to be sure of its predictions, i.e. output low probabilities for
negative class samples and high probabilities for positive class samples.

##  4 Experiments, results, and discussion

In this section, we present the experiment design, results, and discussion.
Experimental setup, including datasets, models and their configurations, is
described in Section 4.1. Results are presented and discussed in Section 4.2.

###  4.1 Experimental setup

####  4.1.1 Datasets

Because of the problems with small datasets described in Section 3.3, we
decided to use datasets from MoleculeNet [114]. The datasets from MoleculeNet
have reasonable size and have predetermined training, validation and test
split, available through the OGB interface. Those splits are often calculated
using scaffold splitting, making them more challenging. Additionally, we limit
ourselves to the datasets having at least 1500 graphs. This allows us to use
the holdout method, without incurring the large computational cost of cross-
validation. We select a varied collection of datasets. Their basic statistics
are summarized in Table 1. The descriptions below closely follow MoleculeNet
paper [114]. The visualizations of class frequency (for classification
datasets) and target value distribution (for regression dataset) are provided
in Figure 11 and Figure 12.

Category | Dataset | # tasks | # graphs | Task | Split | Metric
---|---|---|---|---|---|---
Physical
chemistry | Lipophilicity | 1 | 4200 | Regression | Random | RMSE
Biophysics | HIV | 1 | 41127 | Classification | Scaffold | ROC-AUC
BACE | 1 | 1513 | Classification | Scaffold | ROC-AUC
Physiology | BBBP | 1 | 2039 | Classification | Scaffold | ROC-AUC
Table 1: Datasets used and their basic statistics, taken from [114].

Lipophilicity dataset aims to predict lipophilicity - the ability of a
chemical compound to dissolve in fats and oils. It is important for drugs, as
it affects solubility and membrane permeability. The measure to predict is the
octanol/water distribution coefficient $\log{D}$ at pH $7.4$. It uses random
split, as it is a regression dataset.

HIV dataset presents the task of predicting whether the molecules are active
against HIV. Molecules inhibiting HIV replication are crucial for developing
anti-HIV drugs. Scaffold split is used, as we are interested in discovering
new HIV inhibitors, and therefore have to evaluate previously unseen molecular
structures. It is by a large margin the largest and the most imbalanced
dataset used.

BACE dataset provides binary binding/non-binding results for a set of
inhibitors of human $\beta$-secretase 1 (BACE1). Drugs that block this enzyme
(BACE inhibitors), as they could prevent the buildup of beta-amyloid peptides
in neurons, which are the major drivers of Alzheimer’s disease pathogenesis
[89]. Scaffold split is used, for the similar reason as in HIV dataset.

BBBP dataset concerns prediction of the blood-brain barrier penetration (BBBP)
ability of small molecules. The membrane separating circulating blood and
brain extracellular fluid blocks most drugs, hormones, and neurotransmitters.
Any drugs that aim to directly affect the central nervous system have to
penetrate that barrier. Scaffold splitting is also used, similarly to two
previous datasets.

Figure 11: Positive and negative class percentages in classification datasets
used.

Figure 12: Regression values distribution in the Lipophilicity dataset.

####  4.1.2 Data preprocessing

Since GNNs work by neighborhood aggregation, they operate on node features.
For this reason, choosing proper featurization is essential for GNN
performance. Fortunately, for chemical datasets many atomic features have been
studied and have become a standard in molecular graph classification, in part
because of the OGB initiative [63]. They are one-hot encoded, since
neighborhood aggregation works feature-wise, and also because GIN makes such
assumption about node features [115]. Those features can be easily calculated
by standard RDKit software [92] based on SMILES strings. The full list of
features is available online as part of OGB package [7] and they have been
summarized in Table 2.

Feature |  Values
---|---
Atomic number |  Integer in $[1,118]$ or misc
Chirality |  Tetrahedral-CW, tetrahedral-CCW, unspecified or other
Degree |  Integer in $[0,10]$ or misc
Formal charge |  Integer in $[-5,5]$ or misc
Number of connected hydrogen atoms |  Integer in $[0,8]$ or misc
Number of radical electrons |  Integer in $[0,4]$ or misc
Hybridization |  SP, SP2, SP3, SP3D, SP3D2 or misc
Aromaticity |  True or false
Part of a ring |  True or false
Table 2: Features extracted in OGB for molecular datasets, taken from [7].

Nevertheless, using those raw features would result in huge dimensionality and
extremely sparse features. In many datasets there would be potentially many
constant features, as, e.g. molecules with very high degree atoms (9 or 10)
are rare. For this reason, we can instead embed those features, similarly to
embedding input words in NLP models [6, 8]. This results in smaller, more
dense input vectors. Features can be embedded independently and summed,
creating a fixed length, continuous representation [6]. In PyTorch, the
embedding class also linearly transforms features, which allows creating a
weighted summation, where weights are learned through backpropagation.

For all GNN models, we use embedded atomic features. Embedding dimensionality
is the same as the number of convolutional channels for a given model.
Molecular Fingerprints and Deep Multisets have been proposed for datasets with
pre-calculated features [39], so we check raw features, but we also propose
using embedded features and evaluate this approach, for the reasons outlined
above. While this approach results in lower interpretability in case of MFPs,
it may improve the classification results.

We do not perform any other kind of preprocessing. We use dataset splits as
provided by OGB data loaders, without filtering any molecules.

####  4.1.3 Architectures

As baselines, we use:

  * •

Molecular Fingerprints (MFPs)[39]

  * •

Deep Multisets [39]

  * •

SGC

For MFPs and Deep Multisets, we treat layer sizes in MLPs and learning rate as
hyperparameters. We check MFPs and Deep Multisets both without and with atom
embeddings, as separate models. For SGC, tune use the same hyperparameters as
for other GNNs (see below).

We focus on four most well-known GNN architectures: GCN, GraphSAGE, GIN, and
GAT. They are varied, since they follow different design principles. Firstly,
we perform a set of basic experiments, checking “flat” architectures, i.e.
just $K$ layers stacked one after another, then the readout function and MLP.
The MLP is always a single linear layer, with either softmax activation
(classification) or no activation (regression), which is typical for GNNs.
Additionally, for all GNN-based approaches we tune the following
hyperparameters: the number of convolutional layers, size (number of channels)
of convolutional layers, readout function, dropout in convolutional layers,
learning rate. For all other settings, we use PyTorch Geometric model
defaults, e.g. ReLU activation and no normalization between layers (we found
during initial experiments that adding batch normalization either does not
improve performance or harms it).

We also check the Jumping Knowledge architectures. For each one of the four
basic GNN architectures, we aggregate the outputs of layers with JK: max,
concatenation and LSTM. We do not treat the aggregation type as a
hyperparameter, but rather check 12 different architectures, e.g. GCN + JK
max, GCN + JK concatenation, GCN + JK LSTM. For LSTM, we use default PyTorch
Geometric model defaults.

We evaluate Local Degree Profile (LDP) as a purely topological graph
descriptor approach. We check four variants:

  1. 1.

Original - 5 node features proposed as a baseline

  2. 2.

Extended - with neighbors degree sum and shortest path descriptors

  3. 3.

With additional features - we add clustering coefficient, betweenness
centrality, closeness centrality and PageRank to the extended variant features

  4. 4.

With additional features and feature selection - we take the extended variant
and treat using (or not) each one of additional 4 features as a binary
hyperparameter

We treat the number of bins and whether to normalize the histograms as
hyperparameters.

For fingerprints-based approaches, we check:

  1. 1.

Morgan fingerprints, hyperparameters: radius, number of bits per radius. The
total number of bits passed to RDKit for fingerprint calculation is simply a
multiplication of those two. We choose this approach instead of optimizing the
number of bits directly because this approach gave better results in
preliminary experiments.

  2. 2.

RDKit fingerprints, hyperparameters: max path length, number of bits

  3. 3.

MACCS keys, no additional hyperparameters

  4. 4.

Concatenation of Morgan, RDKit and MACCS fingerprints, simultaneously tuning
hyperparameters for all of them

Inspired by OGB leaderboard architectures and motivated by preliminary
experiments, we apply subsampling the negative class in the training data for
LDP and fingerprints on HIV dataset, which has very significant negative class
dominance. We treat the percentage of retained negative class samples as a
hyperparameter. For both LDP and fingerprints, we use Random Forest as a
classifier/regressor, and we treat minimal number of samples for split as a
hyperparameter. For classification, we use entropy minimization as a split
criterion, and for regression we use squared error minimization. For
classification tasks, we also treat the class weighting approach (none,
balanced, balanced subsample) as a hyperparameter.

For both LDP and fingerprints, we also drop constant features and heavily
correlated features (Pearson correlation coefficient 0.95 or higher). Feature
can be constant, e.g. for MACCS keys, when a particular feature is not present
in the dataset, which can happen, since they are predetermined. For
concatenation of fingerprints, we can get potentially many heavily correlated
features, since particular fingerprints can extract the same substructures.

We summarize hyperparameter optimization configurations in Section 4.1.4. Note
that for GNNs on a large HIV dataset, we check larger convolution sizes.

Out of all models in this section, our original propositions are:

  * •

SGC - as a baseline in the context of graph classification

  * •

Molecular Fingerprints and Deep Multisets - the combination with node features
embedding

  * •

LDP - using well selected additional features, and combining it with feature
selection

  * •

using concatenation of multiple fingerprints - while used before for HIV
dataset, there has been no through study of this method on multiple datasets

####  4.1.4 Model training, selection, and evaluation

To implement our methods, we use:

  * •

loading MoleculeNet datasets: OGB [63]

  * •

GNNs: PyTorch [87] and PyTorch Geometric [40]

  * •

LDP and graph descriptors: PyTorch Geometric [40], NetworkX [52], and Scikit-
learn [88]

  * •

fingerprint-based methods: RDKit [92] and Scikit-learn [88]

  * •

hyperparameter optimization: Optuna [13]

To train GNNs, we optimize cross-entropy loss (classification) or mean squared
error (regression), using the Adam optimizer. We train for 3000 epochs max,
with early stopping, if we do not get better validation result (AUROC for
classification, RMSE for regression) for set number of rounds, we end the
training. The patience parameter is:

  * •

GNN models - 100 rounds

  * •

MFPs and Deep Multisets - 500 rounds, we find that improves results

The batch size is:

  * •

GNN models - 64, a typical value in the literature

  * •

MFPs and Deep Multisets - full batch training for all datasets except HIV,
where we use 4096; we find that this greatly stabilizes training

Those settings have been inspired by the literature, especially by [39], and
validated during preliminary experiments. For heavily imbalanced
classification datasets (HIV and BBBP) we use class weighting, as implemented
in the Scikit-learn library [12].

Since the size of many hyperparameter grids are large, and training the neural
network-based models (Molecular Fingerprints, Deep Multisets, GNNs) takes very
long time, we cannot use grid search for hyperparameter optimization.
Additionally, many parameters are inherently continuous or ordered (e.g.
learning rate, number of layers) and discretizing them may have a negative
effect. Instead of using simple random search, we use the Bayesian
optimization, which is commonly used in recent works because of the excellent
performance and being able to detect good sets of parameters in just a few
iterations. Specifically, we use the Tree Parzen Estimator (TPE) sampling,
implemented in the Optuna [13] library, which allows extremely efficient
optimization. In the literature, using about 20-30 iterations seems the most
popular, but to ensure good exploration of hyperparameter space, we use 40
trials.

For the models using Random Forest (LDP, fingerprints), the model training is
extremely fast (seconds or at most minutes). Moreover, all parameters are
discrete or can be easily discretized. For this reason, we use grid search in
most cases, and fingerprints parameters are all categorically distributed
(since we check all possibilities anyway, Optuna will not benefit from
specifying more specific distribution type). The only exception is the LDP
with additional features and feature selection, since the feature selection
has an exponential complexity ($2^{4}$ possibilities for 4 features). This, in
conjunction with other hyperparameters, results in a very large grid, so in
this case, we use Bayesian optimization with 100 iterations.

In all cases of hyperparameter selection, we select the model with highest
AUROC or lowest RMSE, depending on the dataset task.

For LDP and fingerprints models using RF, we optimize hyperparameters for
feature extraction and for classifier / regressor in a 4-step process:

  1. 1.

We use RF with default Scikit-learn parameters (apart from split criterion) to
choose optimal feature extraction (LDP or fingerprint) hyperparameters. For
HIV dataset, negative class sampling rate is also tuned here.

  2. 2.

The features are calculated for the entire dataset using selected
hyperparameters, creating a regular tabular dataset.

  3. 3.

We tune RF hyperparameters on transformed dataset (without modifying training
/ validation / test splits).

  4. 4.

The final classifier is trained on the training dataset, using 1000 trees and
determined RF hyperparameters.

Hyperparameter ranges are summarized in Table 3.

Architecture |  Hyperparameters
---|---
GNN |  Number of channels: [64, 512] (HIV) or [64, 256] (other datasets), discrete integer distribution with step 64
Number of convolutional layers: [1, 5], discrete integer distribution with
step 1
Readout: ["mean", "max", "sum"], categorical distribution
Dropout: [0.0, 0.3], discrete real distribution with step 0.1
Learning rate: [1e-5, 1e-3], continuous real distribution

Only for GAT architectures:
Number of attention heads: [2, 8], discrete integer distribution with step 2
Molecular Fingerprints, Deep
Multisets |  MLP layer size: [64, 512], discrete integer distribution with step 64
Learning rate: [1e-5, 1e-3], continuous real distribution

LDP |  Number of bins: [10, 100], discrete integer distribution with step 10
Normalize distributions: [false, true], categorical distribution
Only for additional features and feature selection variant:
Add clustering coefficient: [false, true], categorical distribution
Add betweenness centrality: [false, true], categorical distribution
Add closeness centrality: [false, true], categorical distribution
Add PageRank: [false, true], categorical distribution

Only for HIV:
Negative class sample size: [0.3, 0.5], discrete real distribution with step
0.1
Fingerprints |  For Morgan fingerprints:
Radius: [2, 3], categorical distribution
Bits per radius: [1024, 2048, 4096], categorical distribution

For RDKit fingerprints:
Max path length: [5, 6, 7], categorical distribution
Number of bits: [1024, 2048, 4096], categorical distribution

Only for HIV:
Negative class sample size: [0.3, 0.4, 0.5], categorical distribution
Random
Forest |  Minimal samples for split: [2, 10], discrete integer distribution with step 2

Only for classification datasets:
Class weight: [none, "balanced", "balanced_subsample"], categorical
distribution
Table 3: Hyperparameters ranges and distributions.

###  4.2 The results and discussion

Since the tables with full results are very large, they are presented in
Appendix A. Here, we refer to selected, interesting subsets of results for
easier analysis. They have been broken down into separate sections, as they