11institutetext: LAMSADE, Université Paris-Dauphine, PSL
<EMAIL_ADDRESS>
# Scaling up graph homomorphism
for classification via sampling

Paul Beaujean  Florian Sikora  Florian Yger

###### Abstract

Feature generation is an open topic of investigation in graph machine
learning. In this paper, we study the use of graph homomorphism density
features as a scalable alternative to homomorphism numbers which retain
similar theoretical properties and ability to take into account inductive
bias. For this, we propose a high-performance implementation of a simple
sampling algorithm which computes additive approximations of homomorphism
densities. In the context of graph machine learning, we demonstrate in
experiments that simple linear models trained on sample homomorphism densities
can achieve performance comparable to graph neural networks on standard graph
classification datasets. Finally, we show in experiments on synthetic data
that this algorithm scales to very large graphs when implemented with Bloom
filters.

###### Keywords:

Graph embedding Graph homomorphism Subgraph counting.

##  1 Introduction

Statistical learning was first developed to study the performance of
predictive models on Euclidean data where each observation is a point in some
vector space. In domains where the focus is relationships between entities,
data is often represented by graphs of varying size and connectivity. This
mismatch between models operating on vectors and graph data makes it difficult
to adapt existing machine learning techniques to certain application domains.
In particular in chemistry and biology, many structures such as molecules
[15], protein interactions, interactions between genes, and even brain regions
[16] are naturally modeled as graphs.

#### Graph embeddings

There are two main ways to find vector representations for non-Euclidean data:
either computing an embedding of the data into a vector space directly, or
learning this embedding [4]. For example, the bag of words model is a standard
method in natural language processing which maps a sentence to the number of
occurrences of each word that appears in it. This usually sparse high-
dimensional vector is then amenable to further preprocessing or directly used
as an input to machine learning models. On the other hand, richer embeddings
that are obtained by learning such as word2vec [30] or Glove [35] often
outperform direct embeddings and improve the accuracy of several downstream
machine learning tasks. However, their learning procedures require
significantly more computational resources than what is required to compute
simple embeddings like the bag of words model and to train simple models on
the resulting vectors.

There is ongoing research directed towards finding invariant vector
representations of graphs and studying their properties [19]. Among computed
embedding methods, a popular one is to use for each graph the histogram of
colors obtained in the canonical coloring given by the Weisfeiler–Leman color
refinement algorithm [40]. Such histograms of colors serve as vector
representations which are guaranteed to correctly detect isomorphic graphs but
can fail to tell apart two non-isomorphic graphs. While color refinement
provides global information about the graph, other embedding approaches aim to
represent a graph via its substructures. In particular, subgraphs have been a
main point of interest in the design of features used in machine learning
tasks on graphs. Shervashidze et al. have proposed kernel methods based on
graphlet statistics [38] both with exact counting algorithms and sampling
alternatives which amount to studying the distribution of small induced
subgraphs. NT and Maehara have studied the use of graph homomorphism numbers
as graph embeddings [33] which can be seen as counting a specified collection
of partial subgraphs. In the graph kernel literature, statistical information
regarding random walks or specific subgraphs such as shortest paths and
subtrees has been the basis of several methods which can be computed
efficiently [25]. In the field of learned graph embeddings, graph neural
networks remain the most popular technique with a large part of the literature
being focused on finding the best variants of message-passing graph neural
networks [37].

In this work, we contribute to the study of feature generation for graph data,
i.e. we directly compute invariant graph embeddings that will be given as
inputs to standard machine learning operating on Euclidean data [6]. Good
embeddings should be fast to compute, have high representational power, and
ideally allow the user to factor in inductive bias.

#### Limitations of homomorphism numbers

NT and Maehara [33] leverage results by Lovász [28] on graph homomorphisms to
promote the use of homomorphism numbers as graph embeddings that carry
desirable theoretical guarantees. In particular they introduce a class of
embeddings parameterized by a collection of pattern graphs or motifs
$\mathcal{F}$. Counting graph homomorphisms from the patterns $\mathcal{F}$
corresponds to counting partial subgraphs of that collection.

Graph homomorphism numbers are closely related to several classes of existing
graph machine learning techniques. Given an appropriate choice of
$\mathcal{F}$ they are equivalent to color refinement algorithms [9] and to a
large class of graph neural networks [32]. Moreover, NT and Maehara [33] prove
universality theorems which state that if a function operating on graphs is
invariant over graphs that have the same homomorphism numbers from
$\mathcal{F}$, then it can be accurately approximated by a polynomial of
homomorphism numbers from $\mathcal{F}$.

To highlight the practical use cases of these embeddings, NT and Maehara focus
on graph homomorphisms that can be counted efficiently. In particular they
develop custom counting methods for tree homomorphisms and cycle homomorphisms
which are reminiscent of existing kernel methods. They also give an
implementation of the algorithm proposed by Díaz et al. [10] which is
currently the best algorithm for counting homomorphisms from any arbitrary
motif $F$. In particular this algorithm is polynomial if $F$ has bounded
treewidth, that is $F$ is close to being a tree.

However, under commonly used computational complexity assumptions, it is not
possible to design an algorithm with better worst-case time complexity than
the one given by Díaz et al. In practice, this algorithm is fast but cannot be
used realistically for large graphs such as those found in large-scale social
network datasets [26].

#### Our contributions

To remedy the problem of computing homomorphism numbers on large graphs we
attempt to find an alternative that would retain its properties and remain
fast to compute at scale.

An homomorphism density is a normalized version of an homomorphism number
which, instead of the count, represents the frequency at which a given
subgraph $F$ with $k$ nodes appears in a graph $G$ with $n$ nodes. Unlike
homomorphism numbers which can grow as large as $n^{k}$, homomorphism
densities are always contained in the (0,1) interval. For this reason, while
approximating the number of homomorphisms remains a difficult problem from the
point of view of computational complexity, it is still possible to give an
additive approximation of homomorphism densities in polynomial time. Note that
this only goes one way: an approximate number would give an approximate
density but an additive approximation of the homomorphism density does not
allow to recover a good approximation of the number of homomorphisms.

The main contributions of this paper can be summarized as follows. (1) We
demonstrate in graph classification benchmarks on standard datasets that
sample homomorphism densities retain the same representational power as
homomorphism numbers and simple models using those as features are comparable
in performance to popular graph neural network architectures. (2) We show that
a high-performance implementation of a simple sampling algorithm which
computes an $\varepsilon$-additive approximation of homomorphism densities
scales to very large graphs in synthetic experiments.

##  2 An alternative to homomorphism numbers

In this section we discuss existing theoretical results on graph homomorphisms
that justify the relevance of approximate homomorphism densities as graph
embeddings.

###  2.1 The representational power of homomorphism numbers

In the context of graph theory, graph morphisms are functions which map the
node set $V(F)$ to $V(G)$. In particular, graph homomorphisms are morphisms
that preserve adjacency, e.g. $\forall\,uv\in E(F),f(u)f(v)\in E(G)$. A graph
homomorphism that also preserves non-adjacency is called a graph isomorphism.
We write $G_{1}\cong G_{2}$ when there exists a graph isomorphism between
$G_{1}$ and $G_{2}$.

The most important fact about graph homomorphisms is that any undirected graph
$G$ can be uniquely determined (up to isomorphism) by the set of homomorphisms
from pattern graphs to $G$. Following NT and Maehara [33], we use the
notation:

$\operatorname{Hom}(F,G)=\left\\{f:V(F)\to V(G)\mid\forall\,uv\in
E(F),f(u)f(v)\in E(G)\right\\}$

for the set of homomorphisms from a pattern graph $F$ to our target graph $G$.
The cardinality of this set, $\hom(F,G)=\left|\operatorname{Hom}(F,G)\right|$,
is called the homomorphism number of $G$ from $F$. As it turns out, to
identify a graph $G$, it is enough to know its homomorphism numbers. Since we
frequently consider several homomorphism numbers at once it is convenient to
gather them in a vector of homomorphism numbers as follows:
$\hom_{\mathcal{F}}(G)=(\hom(F,G))_{F\in\mathcal{F}}$ where $\mathcal{F}$ is a
set of pattern graphs. This notation allows us to restate the property more
formally as follows:

###### Theorem 2.1

([27]) Given two undirected graphs $G_{1}$ and $G_{2}$ with at most $n$ nodes.
Denoting by $\mathcal{G}_{n}$ the set of all simple graphs with at most $n$
nodes, we have:

$G_{1}\cong
G_{2}\iff\hom_{\mathcal{G}_{n}}(G_{1})=\hom_{\mathcal{G}_{n}}(G_{2}).$

In particular, this theorem can be interpreted as saying that
$\hom_{\mathcal{G}}(G)$ is a canonical vector representation of $G$ that is
invariant to labeling, where $\mathcal{G}$ is the set of all simple graphs,
i.e. graphs without loops or multiple edges. It may appear that we have traded
a finite graph $G$ for a vector of infinite length $\hom_{\mathcal{G}}(G)$.
However Theorem 2.1 already indicates that to tell apart two graphs with $n$
nodes it is enough to look only at the coordinates of the vector that
correspond to pattern graphs with at most $n$ nodes.

A natural question to ask is whether it is possible in practice to consider
significantly less than $\sum_{k=1}^{n}2^{k(k-1)/2}$ coordinates to separate
graphs. In the context of machine learning, this would mean that a short
vector of homomorphism numbers could suffice, e.g. for graph classification.
NT and Maehara [33] provide a preliminary answer to this question and
demonstrate that this approach is adequate even with a small number of small
pattern graphs or a few well-chosen patterns graphs.

###  2.2 The computational complexity of counting graph homomorphisms

Even for a single pattern $F$ computing the value $\hom(F,G)$ which is the
number of homomorphisms from $F$ to $G$ is $\\#P$-hard in general [10] where
$\\#P$-hard is the analogue of $NP$-hardness for counting problems. This
computational complexity is inherited from the hardness of deciding the
existence of a specific subgraph in a graph, e.g. finding a large clique. In
simpler terms, finding the difference between 0 homomorphism from some $F$ to
$G$ and more than 1 homomorphisms is already hard for some $F$.

For many pattern graphs however, it is still possible to compute homomorphism
numbers in polynomial time. For example, if we consider pattern graphs with no
more than $k$ nodes, we can obtain a naive polynomial-time algorithm by
listing out the $O(n^{k})$ subsets of $k$ nodes in $G$ and checking which of
these subsets correspond to $F$. This simple approach is often combined with
clever heuristics in the graphlet literature to obtain high-performance
software to compute subgraph statistics [24]. For specific families of pattern
graphs like cycles or trees, there are known polynomial-time algorithms that
are significantly more efficient. However, in the case of an arbitrary pattern
$F$ the best option is to use the algorithm of Díaz et al. [10]. This
algorithm relies on a tree decomposition of the pattern graph $F$ and has a
worst-case time complexity of $O(\operatorname{tw}(F)\cdot
kn^{\operatorname{tw}(F)+1})$ where $G$ and $F$ have respectively $n$ and $k$
nodes and the treewidth $\operatorname{tw}(F)$ is the size of the smallest
tree decomposition of $F$ [36]. This graph parameter, which is $NP$-hard to
compute, measures how “tree-like” a graph is: lower values signify being very
close to trees (which have a treewidth of 1) such as cycles which have
treewidth 2 while at the extreme opposite the complete graph $K_{n}$ has a
treewidth of $n-1$. Unfortunately, this means that even computing the
homomorphism number from a relatively small pattern graph e.g. $\hom(K_{5},G)$
with the algorithm of Díaz et al. is almost as impractical as listing out
every subset of 5 nodes.

The algorithm of Díaz et al. means that computing homomorphism numbers is
tractable when the treewidth parameter is bounded. A natural question is to
ask whether there are other cases or parameters that make the problem easier.
It turns out that not only this is the only parameter that can make the
problem easier [18] but also that even computing a multiplicative
approximation of homomorphism numbers cannot be done efficiently unless the
treewidth of $F$ is bounded [5]. On a side note, as the number of
homomorphisms can range from 0 to $n^{k}$ it is difficult to define what would
be an efficient additive approximation of homomorphism numbers.

This leaves us with two remaining possibilities for computing homomorphism
numbers: either we limit ourselves to patterns with low treewidth or we rely
on graphs for which a given heuristic would perform well in practice. However
there is one blind spot of computational complexity that we can exploit if we
switch from homomorphism numbers to homomorphism densities.

###  2.3 Efficient sampling for additive approximation of homomorphism density

Homomorphism densities are normalized homomorphism numbers. Formally, if we
write $n=|V(G)|$ and $k=|V(F)|$, we can define the homomorphism density from
$F$ to $G$ as follows:

$t(F,G)=\hom(F,G)/n^{k}.$ (1)

While $\hom(F,G)$ is a magnitude, that is the number of graph homomorphisms
from $F$ to $G$, $t(F,G)$ is a proportion. Alternatively, the homomorphism
density from $F$ to $G$ is the probability that a morphism $f:V(F)\to V(G)$
drawn uniformly at random is a graph homomorphism, i.e.
$t(F,G)=\operatorname{\mathbb{P}}_{f\sim V(F)\to
V(G)}\\{f\in\operatorname{Hom}(F,G)\\}$. Another way to interpret it is the
probability that a random morphism will preserve all edges of the pattern
graph $F$ into the target graph $G$.

Like in Theorem 2.1 where we use homomorphism numbers to distinguish between
two graphs, we can use homomorphism densities $t_{\mathcal{F}}$ as a canonical
graph embedding:

###### Theorem 2.2

(Isomorphism via homomorphism densities) Given two undirected graphs $G_{1}$
and $G_{2}$ with the same number of nodes $n$, and $\mathcal{G}_{n}$ the set
of all simple graphs with at most $n$ nodes, we have:

$G_{1}\cong G_{2}\iff t_{\mathcal{G}_{n}}(G_{1})=t_{\mathcal{G}_{n}}(G_{2}).$

The above can be readily derived from Theorem 2.1 by dividing each equality of
homomorphism numbers $\hom(F_{i},G_{1})=\hom(F_{i},G_{2})$ on both sides by
$n^{|V(F_{i})|}$ for every simple graph $F_{i}\in\mathcal{G}_{n}$.

From section 2.2 it is clear that we should not expect homomorphism densities
to be easier to compute exactly as multiplying by $n^{k}$ would suffice to
recover the corresponding homomorphism densities. Similarly, if we could
obtain a multiplicative approximation of $t(F,G)$ we would obtain the same
multiplicative approximation of $\hom(F,G)$ but that is ruled out by hardness
results.

However, since $t(F,G)$ is by definition the probability of a polynomial-time
testable property, it is possible to turn sampling into a polynomial-time
$\varepsilon$-additive approximation algorithm. Such result would be allowed
despite the previously mentioned hardness results because there is no
guarantee that an additive approximation of $t(F,G)$ can be processed into any
meaningful approximation of $\hom(F,G)$. In particular, unless the precision
is small enough $\varepsilon<n^{k}$ there is no way to distinguish between 0
or 1 homomorphism. Furthermore such additive precision would lead to an
exponential running time. In other terms, an efficient additive approximation
algorithm for homomorphism densities would not be guaranteed to detect very
small densities but would produce accurate estimates for larger densities.

##  3 Implementing a high-performance sampling algorithm for approximate
homomorphism density

In this section we describe an efficient parallel implementation of a sampling
algorithm which gives a $\varepsilon$-additive approximation of homomorphism
densities and we give an analysis of its time complexity. This algorithm
relies on standard sampling arguments similar to those used in sketching
algorithms for graph data [1].

As earlier, we say that $G$ and $F$ have respectively $n$ and $k$ nodes.
Notice that we can sample a morphism $f:V(F)\to V(G)$ uniformly at random from
the space of all functions $V(F)\to V(G)$ by drawing $k$ integers
independently at random from the uniform discrete distribution
$\mathcal{U}(0,n-1)$. Here we represent $f$ by an array of size $k$ where
function application $f(u)$ is realized by accessing the array $f[u]$.
Furthermore, sampling $N$ such morphisms can be done by requesting a
$N$-by-$k$ array instead. Sampling uniform integers is implemented in every
high-performance pseudo-random number generator software suite such as the
randint routine provided by the NumPy library [21]. We are left with computing
the sample mean for the probability of $f$ to be a homomorphism.

However, for the sample mean to be close to the actual mean, we need to make
sure to compute it from a large enough number of samples. In essence,
estimating an homomorphism density $t(F,G)$ is identical to estimating the
unknown bias of a coin, i.e. a Bernoulli distribution of parameter
$p\in[0,1]$. Here, a random morphism corresponds to a coin flip which lands on
heads when the morphism is an homomorphism, while the density $t(F,G)$
corresponds to the unknown parameter $p$ of the biased coin.

Fortunately, it is well known that we can to use Chernoff bounds [7] to derive
a number of samples that is sufficient to reach an $\varepsilon$-additive
approximation of $p$. This leads to the following folklore result:

###### Theorem 3.1

(Sampling lemma) Let $X_{1},\ldots,X_{N}$ be $N$ i.i.d samples of a Bernoulli
distribution with parameter $p\in[0,1]$. The following implication holds:

$N\geq\dfrac{1}{2\varepsilon^{2}}\log\dfrac{2}{\delta}\implies\operatorname{\mathbb{P}}\left(\left|\frac{1}{N}\sum_{i=1}^{N}X_{i}-p\right|>\varepsilon\right)\leq\delta$

where $\varepsilon>0$ is the additive precision of our sampling and
$1-\delta\in(0,1)$ the degree of confidence in our estimate.

0: $G$ an undirected graph on $n$ nodes, $F$ a pattern graph on $k$ nodes and
$l$ edges, $\varepsilon>0$ the requested additive precision,
$1-\delta\in(0,1)$ the desired confidence.

0: $\bar{t}$ such that
$\operatorname{\mathbb{P}}(|t(F,G)-\bar{t}|>\varepsilon)\leq\delta$

1: $N\leftarrow O(\varepsilon^{-2}\log\delta^{-1})$

2: for $i=1$ to $N$ do

3: $f_{i}\sim\left(\mathcal{U}(0,n-1)\right)_{[k]}$

4: end for

5: $\bar{t}\leftarrow\frac{1}{N}$$\sum_{i=1}^{N}\prod_{uv\in
E(F)}\mathds{1}_{E(G)}(f_{i}(u)f_{i}(v))$

6: return $\bar{t}$

Algorithm 1 Sample homomorphism density

Let us take a quick look at the worst-case time complexity of Algorithm 1. The
first step (line 2) amounts to drawing $N$ times $k$ random integers in
$\\{0,\ldots,n-1\\}$ which incurs a running time of $O(Nk\log n)$. Note that
this procedure is embarrassingly parallel and can be split into at least $N$
independent threads. The second step (line 5) is the bulk of the computational
content: for each morphism $f_{i}$ the algorithm must query $l$ edges of $G$
through an adequate data structure. As soon as an edge query fails, i.e. the
morphism fails to preserve that particular edge, the indicator function
returns 0 and short-circuits the product. The worst case happens on actual
homomorphisms which must query all $l$ edges. The running time of this step is
then $O(Nl)$ which amounts to a total time complexity of $O\left((k\log
n+l)\cdot\varepsilon^{-2}\log\delta^{-1}\right)$. While this analysis appears
to price the random number generation (line 2) higher than the homomorphism
tests (line 5), in practice generators benefits from advanced implementations
that can leverage current hardware architectures efficiently. For this reason
we focus on optimizing edge queries with adequate data structures described
later in this section.

We summarize some properties of our sample homomorphism density algorithm. (1)
It can be computed with close to no regard to the size of $G$ as the
dependency is in $O(\log n)$. Furthermore, $G$ does not need to fit in memory
as we only need: (a) to access its number of nodes to initialize our random
number generator and (b) to be able to query random edges. (2) The running
time depends linearly on the number of edges in the pattern graph $l=|E(F)|$
which invites us to consider small pattern graphs. (3) Higher precision
(smaller $\varepsilon$) is expensive and prohibits the use of the sample
homomorphism density to reliably decide the existence of large cliques. We
stress again that additive approximation algorithms generally cannot provide
useful answers to NP-complete decision problems. (4) A positive sample
homomorphism density $\bar{t}>0$ implies that the real homomorphism density is
positive $t(F,G)>0$, but no quantitative comparison can be given without
knowledge of $\varepsilon$ and $\delta$. (5) Higher degrees of confidence are
extremely cheap, making it possible to add extra nines to a $\delta=99.9\%$
without jeopardizing the running time. (6) The free choice of $F$ allows us to
emulate existing techniques such as statistics on paths or walks [25] but also
to consider more complex substructures.

#### Data structures for approximate membership

The problem of set membership is one of the fundamental building blocks in the
study of data structures. A wealth of data structures and algorithms provide
various trade-offs and guarantees going from sequential search in an unsorted
array, to tries, and hash tables. Of particular note is perfect hashing which
allows for constant-time set membership queries at the cost of linear space
(the space required to represent a perfect hashing function is proportional to
the size of the set). In some way, this space complexity is optimal and cannot
be improved upon.

Surprisingly, a simple probabilistic data structure, the Bloom filter [2],
allows to implement set membership close to this optimal space complexity.
Each element of the set is represented by no more than a dozen of bits and the
membership query is realized in constant time by applying a few hash functions
to the element to be queried. Bloom filters are implemented via a array of
bits which are filled at indices corresponding to the outputs of the different
hash functions. The randomness of this data structure is solely contained in
the choice of hash functions and is perfectly deterministic once the hash
functions have been chosen.

Bloom filters alleviate the need to store the elements themselves unlike e.g.
hash tables and other traditional data structures. This convenience in the
specific context of the set membership problem comes at the cost of a false
positive rate. Indeed, when a Bloom filter returns that an element is not in
the set, that answer is always truthful. However a positive answer does not
guarantee that the element is actually in the set.

With regards to computing the sample homomorphism density, each graph $G$ can
thus be compressed into a Bloom filter containing information about its edges.
In practice this means that each edge, usually represented by a pair of
unsigned integers each coded over 32, 64, or 128 bits is compressed down to a
constant dozen of bits, regardless of the number of nodes of $G$. Furthermore,
each homomorphism test corresponds to computing grouped queries which must all
return true. Since the queries are independent, this implies a lower false
positive rate for the grouped query than its constituent queries. Finally, the
Bloom filter associated with $G$ is constructed once as a read-only data
structure which can be queried repeatedly in parallel when computing the
sample homomorphism density of $G$ with regards to queries representing the
edge set of some pattern graph $F$.

##  4 Sample homomorphism densities in graph classification tasks

We conduct numerical experiments to validate the relevance of sample
homomorphism densities as computed graph embeddings for graph classification.
We have implemented Algorithm 1 and our experimental methodology in the Python
3 language leveraging the NetworkX library for graphs [20], the NumPy library
for array processing [21], and the scikit-learn library for machine learning
models [34]. We have split our code into two parts, a Python library on one
hand and a set of Python scripts that depend on that library on the other.
Source code for both parts is provided as a supplementary document. All our
experiments are run on a Linux computer with an 8-core Intel Xeon Skylake
processor clocked at 2.2 GHz with hyper-threading and 14 GB of memory. We have
not used a GPU for computing embeddings or training models.

We start by describing the settings of our experiments as well as our choice
of parameters and datasets. Then, we discuss experimental results that
demonstrate the relevance of sample homomorphism densities in the context of
graph classification.

###  4.1 Sampling in practice

The guarantees provided by the Chernoff bounds implicitly require the use of
an additive precision of $\varepsilon=O(n^{-k})$ to detect the presence of a
single partial subgraph, i.e. in the case where $t(F,G)=1/n^{k}$. However, we
instead sample at fixed levels of precision, with low values of
$\varepsilon\in\\{0.1,0.05,0.01\\}$. This choice introduces two biases. First,
a low precision can be seen as a preprocessing method which filters out low-
frequency patterns from the feature vectors, leaving subgraphs that occur at a
high frequency in the target graph, e.g. $O(10^{-1})$ or $O(10^{-2})$, and for
which we can guarantee a good additive approximation. Second, fixed precision
implies low precision on larger pattern graphs and higher precision on smaller
ones which in some extreme cases can lead to sampling being equivalent to
exhaustive search, e.g. the target graph is relatively small and the pattern
graph has 2 or 3 nodes. This means that only highly frequent larger pattern
graphs can appear in the feature vectors while smaller pattern graphs are
detected even at lower frequencies.

As mentioned above, the cost of achieving high confidence is extremely low
compared to that of precision. However, exploratory testing reveals that
confidence has virtually no impact on graph classification tasks so we keep it
fixed at $\delta=95\%$.

###  4.2 Pattern graphs and weighted homomorphisms

NT and Maehara [33] focus their study on homomorphism numbers for specific
families of patterns, such as trees and cycles. This choice is partly
motivated by the lower computational complexity of computing homomorphism
numbers for these families with the $O(\operatorname{tw}(F)\cdot
kn^{\operatorname{tw}(F)+1})$ algorithm of Díaz et al. [10], as trees and
cycles have respectively a treewidth of 1 and 2. On the other hand, since
Algorithm 1 does not depend on low treewidth or on a modest size of $G$, we
choose to consider graph patterns among small connected graphs of the Atlas of
Graphs provided by the NetworkX library [20]. The first 10 graphs of this list
correspond to all connected undirected graphs with at most 4 nodes. The next
10 graphs are connected undirected graphs with 5 nodes (out of 21 such graphs)
and contain several graphs of treewidth 3. Each component of the feature
vectors thus corresponds to a graph in that list. To study the impact of
information from larger subgraphs we vary our family of pattern graphs
$\mathcal{F}$ with increasing numbers of subgraphs with the following values:
10, 15, 20.

We consider three variants of homomorphism densities, the original unweighted
variant, a weighted variant which assigns to each morphism the product of the
numerical attributes of the target nodes, and finally a second weighted
variant which disregards node attributes if they exist and instead weighs each
node by its degree.

###  4.3 Model evaluation methodology

We follow standard model evaluation methodology with a 10-fold cross-
validation procedure closely following the study of Errica et al. [13]. Each
dataset is first split in 10 blocks, with an inner holdout procedure splitting
each block with a 4:1 training/validation ratio. Each block is validated
independently and once a set of hyper-parameters has been selected, the
corresponding model is trained and scored 3 times to smooth out the randomness
that could be caused by different initial conditions during the training
procedure.

However, since our underlying features are the result of a sampling procedure,
we also average the entire 10-fold cross-validation procedure over 10
independent samples obtained for each dataset. Reported cross-validation
scores, e.g. mean over 10 blocks of the test accuracy of the best models, are
reported as averages over these 10 independent samples, and the standard
deviation reported is the mean over 10 samples of their corresponding standard
deviations.

To perform sensitivity analysis over the different parameters of the sampling
procedure mentioned above, we hold every parameter constant except the topic
of interest. All accessible random seeds are accounted for and controlled,
both in the model evaluation procedure as well as in the feature generation
code. Nevertheless, subroutines present in the scikit-learn library [34] carry
randomness that we cannot control which leads to an unavoidable variance in
our experiment results. From exploratory testing, this variance is only
visible in the third decimal of our reported scores which is why we decide to
not perform an in-depth analysis of its effect.

###  4.4 Models and hyperparameters

The features we consider simply concatenate the number of nodes of a given
graph with an approximate homomorphism density vector. These features are then
given as inputs to multiple machine learning models for classification.
Because we study the properties of homomorphism density features and are not
attempting at challenging the state of the art in graph classification, we
focus on simple models that are well understood. The simplest model we
consider is a logistic regression classifier with $\ell_{1}$, $\ell_{2}$, or
no regularization. We select the regularization parameter $C$ among the
values: $10^{-4},10^{-2},10,10^{4}$. Training is done using the liblinear
solver through the interface provided by scikit-learn. We also study several
non-linear models which we describe in further details in the supplementary
material.

###  4.5 Datasets

To compare our method with existing research, we use publicly available
datasets from the TUDataset collection [31] which can be split in two groups.
First, we consider datasets of graphs obtained from biochemistry applications
such as MUTAG [8]. Note that the MUTAG dataset, widely used in the graph
machine learning literature, only contains 188 graphs with 17.9 nodes on
average. This small dataset size artificially inflates the variance of cross-
validation procedures (a single error costs 5.3% test accuracy). We also
consider datasets of with a larger number of graphs such as NCI1 [39] (4110
graphs, 29.9 nodes on average, 0.2% per error), ENZYMES (600 graphs, 36.2
nodes on average, 1.7% per error) and PROTEINS [3] (1113 graphs, 39.1 nodes on
average, 0.9% per error), or DD [11] (1178 graphs, 284.3 nodes on average,
0.8% per error). A second group consists of datasets obtained from online
social networks, with the COLLAB (5000 graphs, 74.5 nodes on average, 0.2% per
error), REDDIT-BINARY (2000 graphs, 429.6 nodes on average, 0.5% per error),
and IMDB-BINARY datasets (1000 graphs, 19.7 nodes on average, 1.0% per error)
[41].

| MUTAG | NCI1 | PROTEINS | DD | ENZYMES | REDDIT-B | COLLAB | IMDB-B
---|---|---|---|---|---|---|---|---
SGHD-$\mathcal{A}_{10}$ $\varepsilon=0.1$ | $83.6\pm 8.7$ | $62.6\pm 2.9$ | $72.0\pm 4.1$ | $76.2\pm 3.2$ | $21.2\pm 4.1$ | $73.6\pm 3.2$ | $68.1\pm 2.0$ | $63.3\pm 3.7$
SGHD-$\mathcal{A}_{10}$ $\varepsilon=0.01$ | $86.3\pm 7.9$ | $62.7\pm 3.1$ | $72.3\pm 3.6$ | $76.1\pm 3.2$ | $26.3\pm 4.5$ | $73.8\pm 3.2$ | $67.7\pm 2.1$ | $68.7\pm 2.4$
GHC-$\mathcal{T}_{13}$ (NT&M.) | $88.2\pm 7.4$ | $65.4\pm 2.5$ | $70.6\pm 4.7$ | $75.3\pm 3.6$ | $21.7\pm 3.6$ | $84.6\pm 2.3^{*}$ | $62.1\pm 1.9$ | $69.7\pm 4.4$
GHC-$\mathcal{C}_{7}$ (NT&M.) | - | - | - | $76.1\pm 3.9$ | $29.5\pm 3.2$ | - | - | -
DGCNN (Errica+) | - | $76.4\pm 1.7$ | $72.9\pm 3.5$ | $76.6\pm 4.3$ | $38.9\pm 5.7$ | $77.1\pm 2.9$ | $57.4\pm 1.9$ | $53.3\pm 5.0$
DiffPool (Errica+) | - | $76.9\pm 1.9$ | $73.7\pm 3.5$ | $75.0\pm 3.5$ | $59.5\pm 5.6$ | $76.6\pm 2.4$ | $67.7\pm 1.9$ | $68.3\pm 6.1$
ECC (Errica+) | - | $76.2\pm 1.4$ | $72.3\pm 3.4$ | $72.6\pm 4.1$ | $29.5\pm 8.2$ | OOR | OOR | $67.8\pm 4.8$
GIN (Errica+) | - | $80.0\pm 1.4$ | $73.3\pm 4.0$ | $75.3\pm 2.9$ | $59.6\pm 4.5$ | $\textbf{87.0}\pm 4.4$ | $\textbf{75.9}\pm 1.9$ | $66.8\pm 3.9$
GraphSAGE (Errica+) | - | $76.0\pm 1.8$ | $73.0\pm 4.5$ | $72.9\pm 2.0$ | $58.2\pm 6.0$ | $86.1\pm 2.0$ | $71.6\pm 1.5$ | $69.9\pm 4.6$
No Topology (Errica+) | - | $69.8\pm 2.2$ | $\textbf{75.8}\pm 3.7$ | $\textbf{78.4}\pm 4.5$ | $\textbf{65.2}\pm 6.4$ | $72.1\pm 7.8$ | $55.0\pm 1.9$ | $50.7\pm 2.4$
GIN (NT&M.) | $75.5\pm 7.6$ | $76.7\pm 1.2$ | $73.0\pm 1.8$ | - | - | $74.1\pm 2.3$ | $\textbf{75.9}\pm 0.8$ | $70.7\pm 1.1$
GNTK (NT&M.) | $\textbf{89.5}\pm 7.1$ | $\textbf{83.8}\pm 1.6$ | - | - | - | - | - | $\textbf{75.6}\pm 3.9$

Table 1: Test accuracy scores. OOR: out of resources during training, “-”:
experiment not provided in the corresponding study.

###  4.6 Results and discussion

We compare the test accuracy of simple linear models trained on sample
homomorphism density features with recent studies in graph classification
realized with a sound methodology. Numerical results are summarized in Table
1.

The first category of models corresponds to those trained on embeddings
computed from the graph topology without taking into account node features.
Among the numerous variants of our proposed method we report scores of the
simplest model: SGHD-$\mathcal{A}_{10}$. The Sample Graph Homomorphism Density
$\mathcal{A}_{10}$ model is a logistic regression classifier trained on sample
homomorphism densities computed from 10 motifs of the Atlas of graphs with
coarse or medium fixed precision. Note that features are given as inputs to
the model without preprocessing. Directly related to our method is the study
of homomorphism numbers for graph classification of NT and Maehara [33] which
follows a similar model evaluation methodology except for the lack of repeated
scoring of selected models. We consider two of their SVC models based on their
Graph Homomorphism Convolution framework: GHC-$\mathcal{T}_{13}$ which
computes homomorphism numbers from the first 13 trees, and
GHC-$\mathcal{C}_{7}$ from the first 7 cycles.

On the other hand, models in the second category learn embeddings from node
features with implicit access to graph topology. Since our model evaluation
methodology follows closely the work of Errica et al. [13] we have high
confidence in comparing our test accuracy scores with the ones they report.
The study of Errica et al. gives a fair and reproducible assessment of the
performance of popular graph neural network architectures (DGCNN, DiffPool,
ECC, GIN, GraphSAGE) which aggregate node features via a message-passing
training procedure in comparison to a baseline (“No Topology”) that disregards
the graph topology information and only aggregates node features. Furthermore
we also include a different implementation of GIN together with the Graph
Neural Tangent Kernel method as tested in [33].

On most datasets we observe that our SGHD-$\mathcal{A}_{10}$ models achieve
similar or better performance than the GHC models based on homomorphism
numbers. Furthermore, test accuracy scores are very close to those of multiple
GNN architectures which were trained over a maximum of 2 days.

However, it is clear that the dataset NCI1 cannot be properly learned with
graph homomorphism techniques, and that sample homomorphism densities lose
some information compared to homomorphism numbers. Similarly, the ENZYMES
dataset with 6 classes cannot be learned from the graph homomorphism features
utilized here although increasing accuracy or specializing the patterns
towards cycles does improve performance slightly. We note that two GNN
architectures also fail to classify the ENZYMES dataset.

We notice that our SGHD models lose to GHC models on 2 datasets, with an 11
points difference on IMDB-B. It is interesting to note that in the
supplementary material to [33] we find that GHC-$\mathcal{T}_{13}$ achieves a
very similar $73.8\%\pm 2.8\%$ test accuracy when trained after min/max or
max/abs preprocessing while it achieves over 80% accuracy with other scaling
techniques (standard, quantile, power). This hints at the possibility that
additional preprocessing or models with higher capacity may be able to extract
more information from sample homomorphism densities.

The following observations are taken from experimental results detailed in our
supplementary material. We observe a slight increase in test accuracy from
considering weighted homomorphisms for SGHD on biochemistry datasets in every
combination of parameters. On social datasets, which do not have node
features, there is a strong positive effect of sampling at higher precision.
This effect is compounded if we consider higher precision together with more
graph patterns.

Finally, we acknowledge that while comparable to some learned embeddings, our
models using computed embeddings are less competitive than the best methods
that rely primarily, or even exclusively on node features. This latter fact
constitutes a challenge to the field of graph machine learning.

##  5 Testing the scalability of sample homomorphism densities

As mentioned previously, standard benchmark datasets in graph classification
such as those featured in the TUDataset collection [31] as well as more recent
large-scale datasets such as the Open Graph Benchmark [23] or the datasets
introduced by Dwivedi et al. [12] only contain graphs which have at most a few
hundred nodes. While these modern datasets include significantly more graphs
per dataset, each graph remains small. For example, the recent OGB Large Scale
Challenge [22] introduces a new dataset for graph classification which is made
of 3.8M graphs each containing on average 14.5 edges. Our method on the other
hand is designed to be used on large graphs and because of the lack of
corresponding dataset, we settle for experiments conducted on synthetic
graphs.

We describe experiments which attempt to compare two currently available
approaches to compute graph homomorphism features. The first method is the C++
implementation of the algorithm of Díaz et al. [10] given by the homlib
library of NT and Maehara [33]. This algorithm includes a tree decomposition
routine together with a dynamic programming algorithm to compute homomorphism
numbers. Obtaining homomorphism densities is done by dividing by $n^{k}$ which
is the cardinality of all morphisms from a pattern graph $F$ with $k$ nodes to
a target graph $G$ with $n$ nodes. The second method is our implementation of
Algorithm 1 with two variants: exact edge membership queries to an adjacency
list and approximate queries to a Bloom filter representing the edge set of
$G$. Our implementation is written in the Rust programming language [29] and
will be made available as a library at a later date. Our experiments use
Erdős-Rényi random graphs $G(n,p)$ to control the presence of subgraphs via
the edge density $p$. This allows us to manipulate the target homomomorphism
density for many small pattern graphs. Furthermore, this method allows us to
generate families of graphs of varying sizes which retain common properties.
In Figure 1 we report the running time of algorithms computing homomorphism
densities from three small clique patterns: $K_{3}$, $K_{4}$, and $K_{5}$ to
random graphs that are slightly above the connectivity threshold while
remaining relatively sparse with an average degree of $\log^{2}n$.

$10^{2}$$10^{3}$$10^{4}$$10^{5}$$10^{2}$$10^{3}$$10^{4}$$10^{5}$$10^{6}$$n$Running
time (ms)homlib implementation of [10] Algorithm 1
($\varepsilon=$5\text{\times}{10}^{-3}$$) with adjacency listAlgorithm 1
($\varepsilon=$5\text{\times}{10}^{-3}$$) with Bloom filterAlgorithm 1
($\varepsilon=$1\text{\times}{10}^{-2}$$) with Bloom filter Figure 1: Running
time of homomorphism density algorithms w.r.t. $K_{3}\to G(n,\log^{2}n/n)$

The homlib library implements graphs via adjacency lists and as such
demonstrates great performance on smaller graphs. However, its running time is
highly dependent on the homomorphism density of the target graph. Consider for
example a random graph $G(1000,$9\text{\times}{10}^{-3}$)$ (disconnected
w.h.p.) which contains very few triangles, for which homlib computes an exact
triangle density of $7.9\text{\times}{10}^{-7}$ in only 8ms. The running time
drastically increases to 102ms for $G(1000,$4\text{\times}{10}^{-2}$)$
(connected w.h.p.) which has a notably larger density of
$6.4\text{\times}{10}^{-5}$. More importantly, the running time of homlib
scales linearly with the size of $G$ as shown in Figure 1. Surprisingly,
increasing the treewidth parameter of the pattern graph $F$, from triangle to
$K_{5}$ the complete graph over 5 nodes, incurs at most a linear increase in
the running time which does not match the worst-case time complexity.

For the sample homomorphism density implementation of Algorithm 1 we set a
requested additive error of $\pm$1\text{\times}{10}^{-2}$$ (this is not 1%
error) and $\pm$5\text{\times}{10}^{-3}$$ with 95% confidence, i.e. no more
than 5% of the sample densities exceed the error bound. In practice it is
clear this theoretical error bound is extremely conservative. For example if
we consider a graph $G(400,0.05)$ with approximately $4000$ edges, the exact
triangle homomorphism density is $1.26\text{\times}{10}^{-4}$ which should be
well below a “detection threshold” of $\pm$1\text{\times}{10}^{-2}$$. However,
the sample homomorphism density correctly identifies the order of magnitude
with a sample density of $1.94\text{\times}{10}^{-4}$.

In the case of our implementation of Algorithm 1 with Bloom filters we select
a 1% false positive rate. Compared to the implementation with adjacency lists,
the Bloom filter variant appears to be agnostic to the size of $G$.
Surprisingly, we observe a behavior that is comparable to homlib in the sense
that a low homomorphism density is associated with lesser running time. This
can be explained by the fact that the less homomorphisms the higher the
frequency of negative membership queries which short-circuits the entire group
query. As noted in the analysis of its worst-case time complexity, the most
important factor influencing the running time of Algorithm 1 is the requested
additive error $\varepsilon$. Figure 1 reveals that a fixed additive error of
$\varepsilon=\pm$1\text{\times}{10}^{-2}$$ results in running times ranging
from 30ms to 60ms while raising the precision to
$\varepsilon=\pm$5\text{\times}{10}^{-3}$$ increases the running time to
values ranging from 100ms to 200ms. We underline that this additional cost is
needed when the target homomorphism density becomes small. For example, the
triangle homomorphism density of a random graph $G(n,\log^{2}n)$ with
$n=${10}^{5}$$ is $2.1\text{\times}{10}^{-8}$ which cannot be detected (the
sample density is 0) with a precision of
$\varepsilon=\pm$1\text{\times}{10}^{-2}$$. However, a slightly better
precision with $\varepsilon=\pm$5\text{\times}{10}^{-3}$$ gives a coarse
approximation of this value which is “only” off by 2 orders of magnitude. This
detection threshold can also be used as an implicit filter discarding low-
frequency patterns. Finally, we observe like with homlib that the size of $F$
has very little impact on the running time.

We have conducted additional experiments with cuckoo filters [14] and XOR
filters [17] which reproduce the known tradeoffs between fast filter
construction and fast filter querying. However the cuckoo filter
implementation we have tested was slower approximately by 10% to 20% on both
construction and querying than Bloom filters while XOR filters were up to 30%
faster than Bloom filters at the cost of more than double the construction
time.

The above experimental results provide evidences that our implementation of
Algorithm 1 with Bloom filters would allow to scale a standard training
pipeline including feature generation from existing datasets with graphs of
${10}^{2}$ nodes to currently unavailable datasets with graphs of ${10}^{5}$
nodes with no significant increase of the total running time, whereas using
homlib would multiply the running time by at least 3 orders of magnitude.

##  6 Conclusion

In this paper, we have provided evidences that homomorphism numbers do not
scale to large graphs due to limitations inherent to computational complexity
but also from observations of the performance of practical implementations of
state-of-art homomorphism counting algorithms. On the other hand we have shown
that homomorphism numbers and homomorphism densities have equal
representational power from a theoretical point of view. This is also true in
practice when comparing homomorphism numbers and approximate homomorphism
densities as features for graph classification where simple models using them
as input attain comparable performance to popular graph neural network
architectures. Finally we have demonstrated in synthetic experiments that our
high-performance implementation of a simple algorithm for approximate
homomorphism densities using Bloom filters is highly scalable.

As future works, we are interested in studying the families of pattern graphs
used to compute homomorphism densities and would like to explore whether these
families could be learned. A second direction would be to explore model
explainability and the relationship between learned model weights and
homomorphism information. Furthermore, graph homomorphisms are known to form
an algebra [28], where the product of homomorphism numbers of $G_{1}$ and
$G_{2}$ corresponds to the homomorphism of a graph product between $G_{1}$ and
$G_{2}$. This invites further investigation into whether polynomial kernels or
other non-linear models implicitly capture unseen homomorphism information as
well as studying its potential impact on classification performance. Finally,
given the scalability of our algorithm to large graphs, we we would like to
construct datasets for graph classification that contain sufficiently large
graphs in order to fairly compare the performance of commonly used models on
large-scale graph data.

## References

  * [1] K. J. Ahn, S. Guha, and A. McGregor.  Graph sketches: sparsification, spanners, and subgraphs.  In PODS, 2012.
  * [2] B. H. Bloom.  Space/time trade-offs in hash coding with allowable errors.  Communications of the ACM, 13(7):422–426, 1970.
  * [3] K. M. Borgwardt, C. S. Ong, S. Schönauer, S. Vishwanathan, A. J. Smola, and H.-P. Kriegel.  Protein function prediction via graph kernels.  Bioinformatics, 2005.
  * [4] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst.  Geometric deep learning: going beyond Euclidean data.  IEEE Signal Processing Magazine, 2017.
  * [5] A. Bulatov and S. Živný.  Approximate counting CSP seen from the other side.  ACM Transactions on Computation Theory, 12:1–19, 05 2020.
  * [6] H. Cai, V. W. Zheng, and K. C.-C. Chang.  A comprehensive survey of graph embedding: Problems, techniques, and applications.  IEEE Trans. Knowl. Data Eng., 2018.
  * [7] H. Chernoff.  A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations.  The Annals of Mathematical Statistics, 1952.
  * [8] A. K. Debnath, R. L. Lopez de Compadre, G. Debnath, A. J. Shusterman, and C. Hansch.  Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity.  J. Med. Chem, 1991.
  * [9] H. Dell, M. Grohe, and G. Rattan.  Lovász meets Weisfeiler and Leman.  In ICALP, 2018.
  * [10] J. Díaz, M. Serna, and D. M. Thilikos.  Counting H-colorings of partial k-trees.  Theor. Comput. Sci., 2002.
  * [11] P. D. Dobson and A. J. Doig.  Distinguishing enzyme structures from non-enzymes without alignments.  J. Mol. Biol, 2003.
  * [12] V. P. Dwivedi, C. K. Joshi, T. Laurent, Y. Bengio, and X. Bresson.  Benchmarking graph neural networks.  arXiv preprint arXiv:2003.00982, 2020.
  * [13] F. Errica, M. Podda, D. Bacciu, and A. Micheli.  A fair comparison of graph neural networks for graph classification.  In ICLR, 2019.
  * [14] B. Fan, D. G. Andersen, M. Kaminsky, and M. D. Mitzenmacher.  Cuckoo filter: practically better than Bloom.  In Proceedings of the 10th ACM International on Conference on emerging Networking Experiments and Technologies, pages 75–88, 2014.
  * [15] B. Gaüzère, L. Brun, and D. Villemin.  Graph kernels in chemoinformatics.  In Quantitative Graph Theory: Mathematical Foundations and Applications, chapter 15. Taylor & Francis Group, 2014.
  * [16] J. Gonzalez-Astudillo, T. Cattai, G. Bassignana, M. Corsi, and F. D. V. Fallani.  Network-based brain computer interfaces: principles and applications.  J. Neural Eng., 2020.
  * [17] T. M. Graf and D. Lemire.  Xor filters: Faster and smaller than Bloom and cuckoo filters.  Journal of Experimental Algorithmics (JEA), 25:1–16, 2020.
  * [18] M. Grohe.  The complexity of homomorphism and constraint satisfaction problems seen from the other side.  J. ACM, 54(1), Mar. 2007.
  * [19] M. Grohe.  word2vec, node2vec, graph2vec, x2vec: Towards a theory of vector embeddings of structured data.  In PODS, 2020.
  * [20] A. Hagberg, P. Swart, and D. S Chult.  Exploring network structure, dynamics, and function using NetworkX.  In Python in Science Conference (SciPy), 2008.
  * [21] C. Harris et al.  Array programming with NumPy.  Nature, 2020.
  * [22] W. Hu, M. Fey, H. Ren, M. Nakata, Y. Dong, and J. Leskovec.  OGB-LSC: a large-scale challenge for machine learning on graphs.  arXiv preprint arXiv:2103.09430, 2021.
  * [23] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec.  Open Graph Benchmark: datasets for machine learning on graphs.  arXiv preprint arXiv:2005.00687, 2020.
  * [24] K. Jamshidi, R. Mahadasa, and K. Vora.  Peregrine: a pattern-aware graph mining system.  In European Conference on Computer Systems, 2020.
  * [25] N. M. Kriege, F. D. Johansson, and C. Morris.  A survey on graph kernels.  Applied Network Science, 2020.
  * [26] J. Leskovec and A. Krevl.  SNAP Datasets: Stanford large network dataset collection.  http://snap.stanford.edu/data, June 2014.
  * [27] L. Lovász.  Operations with structures.  Acta Mathematica Academiae Scientiarum Hungarica, 1967.
  * [28] L. Lovász.  Large networks and graph limits.  American Mathematical Soc., 2012.
  * [29] N. D. Matsakis and F. S. Klock.  The Rust language.  ACM SIGAda Ada Letters, 34(3):103–104, 2014.
  * [30] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean.  Distributed representations of words and phrases and their compositionality.  In NIPS, 2013.
  * [31] C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel, and M. Neumann.  TUDataset: A collection of benchmark datasets for learning with graphs.  In ICML Workshop on Graph Representation Learning and Beyond (GRL+), 2020.
  * [32] C. Morris, M. Ritzert, M. Fey, W. L. Hamilton, J. E. Lenssen, G. Rattan, and M. Grohe.  Weisfeiler and Leman go neural: Higher-order graph neural networks.  In AAAI, 2019.
  * [33] H. NT and T. Maehara.  Graph homomorphism convolution.  In ICML, 2020.
  * [34] F. Pedregosa et al.  Scikit-learn: Machine learning in Python.  JMLR, 2011.
  * [35] J. Pennington, R. Socher, and C. D. Manning.  Glove: Global vectors for word representation.  In EMNLP, 2014.
  * [36] N. Robertson and P. D. Seymour.  Graph minors. II. Algorithmic aspects of tree-width.  Journal of algorithms, 1986.
  * [37] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini.  The graph neural network model.  IEEE Trans. Neural Netw, 2008.
  * [38] N. Shervashidze, S. Vishwanathan, T. Petri, K. Mehlhorn, and K. Borgwardt.  Efficient graphlet kernels for large graph comparison.  In AIStat, 2009.
  * [39] N. Wale, I. A. Watson, and G. Karypis.  Comparison of descriptor spaces for chemical compound retrieval and classification.  Knowl. Inf. Syst., 2008.
  * [40] B. Weisfeiler and A. Leman.  The reduction of a graph to canonical form and the algebgra which appears therein.  NTI, Series 2, 1968.
  * [41] P. Yanardag and S. Vishwanathan.  Deep graph kernels.  In KDD, 2015.


#  EPick: Multi-Class Attention-based U-shaped Neural Network for Earthquake
Detection and Seismic Phase Picking

Wei Li, Megha Chakraborty, Darius Fenner, Johannes Faber, Kai Zhou, Georg
Rümpker, Horst Stoecker and Nishtha Srivastava Wei Li and Nishtha Srivastava
are with Frankfurt Institute for Advanced Studies, 60438 Frankfurt am Main,
Germany (e-mail: srivastava@fias.uni-frankfurt.de.)Georg Rümpker and Megha
Chakraborty are with Frankfurt Institute for Advanced Studies, 60438 Frankfurt
am Main, Germany, and also with Department of Geosciences, Goethe-University,
60438 Frankfurt am Main, Germany.Darius Fenner is with Frankfurt Institute for
Advanced Studies, 60438 Frankfurt am Main, Germany, and also with Johannes
Gutenberg-Universität Mainz, 55122 Mainz, Germany.Johannes Faber and Kai Zhou
are with Frankfurt Institute for Advanced Studies, 60438 Frankfurt am Main,
Germany, and also with Institute for Theoretical Physics, Goethe Universität,
60438 Frankfurt am Main, Germany.Horst Stoecker is with Frankfurt Institute
for Advanced Studies, 60438 Frankfurt am Main, Germany, with Xidian-FIAS
International Joint Research Center, and also with Institute for Theoretical
Physics, Goethe Universität, 60438 Frankfurt am Main, Germany, and also with
GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt, Germany.

###### Abstract

Earthquake detection and seismic phase picking not only play a crucial role in
travel-time estimation of body-waves (P and S waves) but also in the
localisation of the epicenter of the corresponding event. Generally, manual
phase picking is a trustworthy and the optimum method to determine the phase-
arrival time, however, its capacity is restricted by available resources and
time. Moreover, noisy seismic data renders an additional critical challenge
for fast and accurate phase picking. In this study, a deep learning-based
model, EPick, is proposed which benefits both from U-shaped neural network
(also called UNet) and attention mechanism, as a strong alternative for
seismic event detection and phase picking. On one hand, the utilization of
UNet structure enables addressing different levels of deep features. On the
other hand, attention mechanism promotes the decoder in the UNet structure to
focus on the efficient exploitation of the low-resolution features learned
from the encoder part to achieve precise phase picking. Extensive experimental
results demonstrate that EPick achieves better performance over the benchmark
method, and show the model’s robustness when tested on a different seismic
dataset.

###### Index Terms:

Seismic phase picking, U-shape neural network, attention mechanism.

##  I Introduction

In Earthquake monitoring, event detection and phase picking are the first
major steps for travel time estimation of P- and S-wave arrivals. This
information is further used for event localization and in seismic tomography.
The highly complex earthquake waveform pattern is defined by the convolution
of source (source time function), path (attenuation, scattering, dispersion
etc.) and site (changes in waveform due to local site conditions) effects.
Nonetheless, overlapping of different phases and poor signal-to-noise ratio
(SNR) often makes it difficult to define a clear seismic onset [1]. With the
abundant and ever increasing database of seismic event recordings, manually
performing the task of earthquake signal detection and seismic phase picking
is not only challenging but also very time-consuming and labor-intensive.

To achieve reliable automatic phase picking, a wide spectrum of traditional
automatic pickers have been developed such as short-term average/long-term
average (STA/LTA) [2], auto regression-Akaike information criterion (AR-AIC)
pickers [3] and higher-order statistics [4]. In particular, the first two
approaches require intensive human involvement. For example, STA/LTA requires
experts to carefully set up parameters whereas the STA/LTA ratio is sensitive
to the detection threshold, and the AR-AIC picker mainly depends on the
STA/LTA trigger. Besides, they can not take advantages of the prior knowledge
of previous picks due to the fact that each measurement in these two methods
is treated individually. While these automatic picking algorithms are useful
in seismological studies, their performances are less accurate than manual
picks, which makes them inefficient in seismic monitoring [5]. Meanwhile, when
applied to real-time seismic data, the accuracy of traditional automatic
pickers may not be satisfactory, particularly for noisy data.

The increasing number of seismic sensors deployed for earthquake monitoring
produces a huge amount of seismic data every day, making data flow and
processing along with defining the manual features for traditional automated
method more difficult and time consuming. Therefore, earthquake monitoring has
an increasing need for more efficient and robust tools to process the large
volumes of seismic data. Deep learning, as an effective method, has achieved
widespread success in a broad range of applications [6, 7, 8]. Inspired by the
success of those applications, phase picking attracts a new wave of deep
learning applications in seismology. In contrast to traditional automated
methods where only a limited set of defined features of seismograms are used,
with the help of deep learning, more abundant features can be successfully
extracted from seismic data and further be used for phase picking. Recent
years have witnessed remarkable achievements of the application of deep
learning in seismic processing tasks especially for earthquake detection and
seismic phase picking [9, 10, 11, 12, 13, 14, 15, 16, 17]. For example, Perol
et al. [9] proposed ConvNetQuake based on convolutional neural network (CNN)
[18] to perform event detection and clustering based on three component
seismic waveform. Their proposed approach is mainly compared with the
similarity search methods and proved to obtain superior performance over
conventional methods in less computational time. Zhu et al. [14] proposed and
trained a deep learning-based phase picker for the Wenchuan earthquake data,
and then fine-tuned the trained model for the Oklahoma dataset. In [17], with
the advent of transfer learning, a trained deep neural network (DNN) phase
picker on local seismic data is retrained to boost performance.

Despite significant progress in deep learning-empowered phase picking
techniques, it still remains a major challenge to develop a comprehensive
approach that not only achieves the satisfied accuracy in seismological
studies, but also widely adapts to small-magnitude seismic events. Moreover,
different forms including Gaussian, box, and triangle are used to label the
phase arrival time when training the neural network in recent works [11] [16]
[19]. For example, Zhu et al. applied a mask with the shape of a Gaussian
distribution around the manual picks for the purpose that the ground truth
arrival times will be centred on the manual picks with some uncertainty.
However, it has the potential to introduce a bias in the training process.

Figure 1: Example of a seismic data [20] where the plots in the figure from
top to bottom represent the three-component seismic recordings (including
east-west, north-south, and vertical directions), respectively. The red and
blue vertical lines denotes the first arrival times of picked P and S phase
from the earthquake catalog.

In this work, a new deep learning-based model, EPick, is proposed for
earthquake signal detection and phase picking. In contrast to the previous
methods of setting labels with a time window, in this paper, only sample
labeling will be adopted to train the neural network. Fig. 1 gives an example
of a seismic data marked with phase arrival time. In summary, the main
contributions of this work are as follows:

  * •

The UNet associated with skip connection is used to perform earthquake signal
detection and seismic phase arrival time picking where the ground truth
arrival time is labeled by one sample, not a time window.

  * •

Considering the heterogeneities of low-level detailed feature maps from the
encoder, multi-class attention mechanism is incorporated into vanilla UNet
structure to effectively explore different scale-level features to improve
phase picking performance.

  * •

Given a global seismic dataset, extensive experiments are carried out to
certify the model performance. Meanwhile, in order to evaluate the model’s
generalisability, the proposed model is tested on a new seismic dataset and
compared with several previous methods.

The structure of the paper is organized as follows. Section II introduces the
basic knowledge related to our method, involving UNet and attention mechanism.
The proposed method details are presented in Section III. Section IV gives a
general view of the dataset used in the article. The model performance of the
proposed approach is analysed in Section V. Section VI makes a conclusion for
the paper.

##  II Related Works

U-shaped Neural Network (UNet). UNet was originally proposed to perform
biomedical image segmentation [21]. It mainly consists of two parts: a
contracting path and a expansive path, which shares similar spirit with the
encoder-decoder architectures.

  * •

The encoder comprises several convolution modules to encode the input with
feature representations of multiple different levels, where each convolution
module is followed by a maxpool downsampling operation.

  * •

The decoder is composed of several deconvolutional modules to perform
upsampling associated with concatenation operations. It aims to semantically
project the higher resolution features extracted by the encoder into upsampled
feature space to do a dense classification.

In the last few years, UNet has been used in the area of seismology for phase
picking. For instance, the architecture of PhaseNet proposed in [11] is a
modified UNet that utilizes a fully convolutional encoder-decoder network
associated with skip connections, for the picking of P and S phases. Zhao et
al. [22] developed a U-shaped neural network for distinct phase detection (P
and S phases) and arrival time picking. In [23], the RLU-Net network is
introduced to pick the first arrival of micro-seismic signals with low SNR.

Attention Mechanism. In the field of natural language processing (NLP),
attention mechanism is proposed to enhance the performance over the encoder-
decoder driven machine translation system [24]. Specifically, an attention
mechanism helps a neural networks focus only on those useful aspects of input
data that boost prediction performance. It succeeds in splitting complex tasks
into small regions of attention and then processing those small tasks
sequentially. Recently, this mechanism including its variants has achieved
superior performance in wide applications, including computer vision, speech
processing, etc. Two most common attention techniques used in attention
mechanisms are presented.

  * •

Self-Attention is an attention mechanism that relates different positions of a
single sequence to compute a representation of the sequence [24]. It enables
the input interacting with itself, calculating the attention scores, and
finally achieving the aggregated output by using these interactions and
attention scores.

  * •

Multi-Head Attention is a mechanism that can be utilized to improve the
performance of self-attention layer [25]. Unlike self-attention where the
attention is only computed once, in the the multi-head mechanism, the process
of the scaled dot-product attention runs multiple times in parallel. As
illustrated in [24], those independent attention outputs are simply
concatenated and linearly transformed into the expected dimensions.

Despite the common use of UNet in seismology, there are still bottlenecks
(i.e., insufficient utilization of information flow) left for UNet, which
impede the potential applications of the raw UNet architecture. To tackle the
inadequate utilization of different level features, in this work, attention
mechanism is incorporated into raw UNet structure, given the fact that
features generated at different stages often possess different levels of
discrimination. Moreover, in contrast to previous methods, when performing
model training, the target phase is located without using a time interval.

Figure 2: EPick model architecture. The model input is the 1-minute three-
component seismic waveforms with a sampling rate 100 Hz. The output is
composed of three probabilities corresponding to noise, P pick and S pick. The
colored blocks represent different layers of the neural network, where the
extracted features will go through sets of transformations. The whole
architecture can be broadly divided into an encoder and a decoder. The encoder
consists of several down-sampling operations (including a 1D convolution block
and stride). The decoder contains several up-sampling processes conducted by
deconvolution layers to recover the feature length of the previous stage. The
skip connection after attention mechanism at each stage enables the network to
concatenate the low-level features to the right block. Here, the skip
connection directly copies the output of attention mechanism (denoted as the
blue modules) to the right layer. The last layer adopts the classical softmax
function to calculate the probabilities of the above-mentioned three classes.

##  III Methodology

In this work, we propose EPick for earthquake detection and phase picking
combined with attention mechanism applied to reduce long-range dependencies
between small-scale features and effectively use low-level features in UNet.
This section presents a detailed description of the EPick architecture.

The model architecture is illustrated in Fig. 2, which is based on the
fundamental UNet structure, combined with an attention mechanism. Similar to
[21], EPick can be roughly decomposed as an encoder network followed by a
decoder network. However, different from [21], the exploitation of the
attention mechanism helps EPick sufficiently use those low-level features
extracted in the encoder part.

The fixed-length three-component waveform is represented by a three-channel
1-dimensional vector, which is fed as the input to EPick. EPick classifies
noise, first P-arrival and first S-arrival. One data example is illustrated in
Fig. 1. Here each channel of the seismic waveform is sampled at 100 Hz of 60
seconds time length. We label the output of our model as $Y_{i}\quad(i=0,1,2)$
to indicate the three class labels of “noise”, “first P-arrival” and “first
S-arrival”. It should be noted that the “noise” here refers to data samples
that are not first arrivals of P or S waves in the metadata, which would be
identified as class $Y_{0}$. Class $Y_{1}$ and class $Y_{2}$ correspond to the
provided first arrivals of P phase and S phase in the metadata.

In the training process, these discrete classes will be embedded as one-hot
encoding that maps a variable to a binary vector, where the length of the
binary vector equals to the number of categories. For instance, for the label
$i$ ( $i$ is the ground truth class), its one-hot encoding label can be
represented as the vector whose element being 1 at index $i$, and 0 for
remainders. Therefore, in this work, the one-hot encoding for the defined
categories is formed as follows:

$Y_{i}=\begin{cases}Y_{0}:[1,0,0]\\\ Y_{1}:[0,1,0]\\\
Y_{2}:[0,0,1]\end{cases}$ (1)

When seismic data is fed into the neural network, it undergoes several down-
sampling and up-sampling modules. Each module is comprised of 1D convolutional
block associated with rectified linear unit (ReLU) activation function. The
goal of down-sampling operation is to learn scale-invariant features from
seismic data. These extracted features are forwarded to the up-sampling stage
such that at each moment, one can obtain the corresponding probability
distribution over three classes. A skip connection involving an attention
block at each stage enables the concatenation of the features learned from the
encoder part to the up-sampling stage. Similar to [11], the deconvolution
operation [26] is used to achieve the up-sampling with the aim to recover
previous feature size. Besides, padding operations are also performed before
and after convolutions to keep the output size the same as the input size.

To be able to output the probability distribution over three classes at each
moment for each sample, a softmax function [27] is used as the final layer in
the network. The fundamental process of the softmax module is to convert the
output representation of the decoder part into probability within interval
$(0,1)$. In this work, the detection process is regarded as multiple
classification problem. Hence, the processed results of EPick on each data
sample can be mapped into the probability by using softmax function as below:

$Y_{i}=\frac{e^{z_{i}(x)}}{\sum_{k=1}^{3}e^{z_{k}(x)}}$ (2)

where $i=1,2,3$ denotes the noise, first P-arrival and first S-arrival; and
$z(x)$ represents the processed results of EPick before using softmax
function.

In this study, the cross-entropy [28] between the predicted label and the
ground truth label is utilized to compute the loss.

$L(Y,Y^{\prime})=-\sum^{3}_{i=1}\sum^{n}_{j=1}Y_{ij}\cdot\log(Y^{\prime}_{ij})$
(3)

where $n$ denotes the sampling number.

Moreover, in order to avoid overfitting, weight decay [29], a popular
regularization techniques, is adopted in the model training. Here, $L_{2}$
regularization [30], the default implementation of weight decay, is added to
the loss function. Eventually, the total loss function is given as

$\mathcal{L}=L(Y,Y^{\prime})+w_{d}*L_{2}$ (4)

where $w_{d}$ denotes the weight decay factor.

The data used in this work has a class imbalance problem [31], where the class
“noise” far exceeds the other two classes, which poses a challenge to build a
reliable classification model statistically. To deal with this issue, a weight
is introduced for each class which helps the model place more emphasis on the
minority classes during the model training process. It aims at ensuring a
classifier that is capable of learning equally from all categories.

##  IV Data and Model Training

For this study, several subsets of STanford EArthquake Dataset (STEAD) [20]
are used to train and test our EPick neural network. STEAD is a high-quality
and global seismic dataset involving non-earthquake and earthquake signals of
labeled P- and S-arrivals. Besides, it has been used in previous papers to
develop and compare detection and phase picking approaches for local
earthquakes.

Figure 3: Frequency-magnitude distributions of earthquakes events used in the
analysis.

In our study, about 100,000 samples including earthquake and noise waveforms
are used to train the proposed model. Here, the seismic event data is selected
based on the unique source id in the dataset [20]. It should be noted that in
contrast to [11][16][19], in this article, we do not use time offsets to
generate the arrival times of P and S phases, especially in the model training
phase. The ratio between training data and testing data is around $4:1$. To
evaluate the earthquake detection and phase picking performance of EPick,
25,000 test seismic waveforms (including 20,000 earthquake and 5000 pure-noise
examples) are utilized, wherein a comparison is made between EPick and the
baseline method UNet without attention modules.

Figure. 3 shows the earthquake magnitude distribution of the training and
testing dataset, respectively. Note that for the comparison, the baseline
model used here are trained on the same training set and then is applied to a
common test set from STEAD.

Fig. 4 gives examples of training instance visualizations with their
corresponding label including earthquake and non-earthquake signals. The three
plots at the bottom of each subfigure represent the normalized three component
seismograms. The upper three plots denote the label information and then each
sample is labeled with a one-hot encoding vector.

(a) Earthquake signal.

(b) Non-earthquake signal.

Figure 4: Training examples [20]. The three plots at the bottom of each
subfigure represents the normalized three component seismograms. The upper
three plots denote the label information and then show each sample encoded in
a one-hot vector.

In this work, the proposed model is constructed and implemented in TensorFlow
[32], and is trained and tested on a Nvidia A100 GPU. The Adam optimizer is
used for optimizing the method. In the model training process, we stop the
model training if the training accuracy does not improve for around 10,000
successive steps.

##  V Results and Model Generalization

In this section, the results of extensive evaluations are presented to denote
the performance of the proposed model. Please note that in the model testing
phase, a time interval is employed to locate the predicted phase pick.
Moreover, due to the interference of P-wave signals with the later arriving
shear waves, it is more difficult to pick first S arrival than first P arrival
[33], hence, a larger time interval is chosen for first S arrival picking. The
time interval for first P arrival is $\pm$ 0.1s from the true label while for
first S arrival it is $\pm$ 0.2s.

The proposed model is also compared with previous methods on the INSTANCE
dataset [34], which is chosen as an unbiased way to evaluate and compare the
performance of EPick with previous methods.

###  V-A Metrics

In order to evaluate the performance, a confusion matrix [35], a powerful
analytical tool in deep learning and data science, is used. The confusion
matrix is capable of displaying detailed information about how a deep learning
classifier has performed with respect to the target classes in the dataset.
Table I gives a visualization of the confusion matrix.

TABLE I: Definition of confusion matrix. | True values
---|---
| Positive | Negative
Predicted values | Positive | True Positive (TP) | False Positive (FP)
Negative | False Negative (FN) | True Negative (TN)

The confusion matrix shows the examples that have been correctly classified
against those examples that have been incorrectly classified. Given a
confusion matrix, different numbers of true-positive samples (TP), true-
negative samples (TN), false-positive samples (FP), and false-negative (FN)
samples can be easily seen. Then, the following metrics can be defined using
the confusion matrix.

Accuracy refers to the ratio between the sum number of true positives and true
negatives and the number of all instances. It is used to measure the ratio of
the correctly predicted samples to all samples. The accuracy metric can be
formulated as below:

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$ (5)

Precision is computed as the fraction between the amount of positive samples
correctly identified to the total number of samples correctly or incorrectly
categorized as positive. It aims to measure the model’s accuracy in
identifying a sample as positive.

$Precision=\frac{TP}{TP+FP}$ (6)

Recall describes how many of the real positive instances a model is able to
identify correctly. In summary, it is the ratio of correctly predicted
positive cases to the all observations in true positive examples.

$Recall=\frac{TP}{TP+FN}$ (7)

F1-score. When handling imbalanced datasets where one set of classes dominates
over another set of classes, accuracy is not the best option to evaluate the
model performance. In contrast to accuracy, the F1-score is usually more
effective, since it is able to punish those extreme values more by using
harmonic mean in place of arithmetic mean. Besides, it succeeds in capturing
both the trends of Precision and Recall into one single value. The following
equation is used to formulate the F1-score through combining both of the
metrics we have reviewed.

$F1=2*\frac{Precision*Recall}{Precision+Recall}$ (8)

Picking Error refers to the time residuals $t$ (also called “time difference”)
between picks of the deep learning model and ground truth.

(a) Unet.

(b) EPick.

Figure 5: Confusion matrix for earthquake event detection in the testing data for UNet and EPick. In the testing, EPick misclassifies less noise as earthquake events. TABLE II: Accuracy, Precision, recall, and F1-score for non-earthquake and earthquake signals detection. Model | Accuracy | Precision | Recall | F1-score
---|---|---|---|---
UNet[22] | 0.9756 | 0.9704 | 1.0 | 0.9850
EPick | 0.9766 | 0.9716 | 1.0 | 0.9856
TABLE III: Evaluation metrics on the test dataset.

|  | Mean (s) | St. Dev. (s) | Precision | Recall | F1 Score
---|---|---|---|---|---|---
P phase | EPick | 0.0095 | 0.1649 | 0.97 | 0.96 | 0.96
| UNet[22] | 0.0094 | 0.1769 | 0.97 | 0.97 | 0.97
S phase | EPick | 0.0036 | 0.1488 | 0.95 | 0.95 | 0.95
| UNet[22] | 0.0117 | 0.2627 | 0.94 | 0.94 | 0.94

  * •

Mean and standard deviation (st. dev.) of the differences of model predicted
arrivals minus manual picked arrivals in seconds

Figure 6: Distribution of time residuals (also called time difference, $\Delta
t$) of UNet and EPick on the test data set.

(a) Magnitude distribution of seismic data.

(b) Seismic data visualization marked with different time intervals.

Figure 7: Magnitude distribution and seismic signal visualization of the used
subset of INSTANCE dataset [34]. (a) Frequency-magnitude distribution of the
test data. (b) 60s seismic waveform is sampled at 100Hz. The evaluated models
are using different lengths of the same seismic event marked with different
dashed lines.

###  V-B Results

First, Figure 5 shows the confusion matrix for the involved models: UNet and
EPick. Furthermore, the overall detection performance for non-earthquake and
earthquake events is shown in Table II. As seen from Table II, EPick improves
slightly in all categories over UNet (without using attention mechanism).
Compared to earthquake event detection, phase picking is much more difficult.

Table III lists the experimental results of both EPick and UNet tested in this
study for seismic phase arrival time picking. The means and standard
deviations of arrival-time residuals between ground truth picks and model
detected picks are computed for model performance evaluation. The before
mentioned statistics suggest that EPick achieves a better performance over
UNet when dealing with seismic phase picking.

Figure. 6 displays the distribution of picking error between the predicted and
ground truth first P and first S picks for EPick and UNet, respectively.
Similar to [11], we also find that the picking error distribution of first P
picks is narrower than first S picks, which is in line with [33].

###  V-C Model generalization

In deep learning, the model generalization describes how well a trained model
is performing on unseen data, which is regarded as one of the most important
criteria in practical applications. To investigate the generalization
abilities on a separate dataset, the trained model is tested on a subset of
the INSTANCE dataset [34], where seismic waveform data are collected from weak
and strong motion stations that have been extracted from the Italian EIDA
node. Consider that the time interval of $\pm$ 0.1s is used in PhaseNet model
[11] for first P and S arrival picking and in ARRU [19], the standard
deviation of the target Gaussian function is 0.2s for the P picks and 0.3s for
the S picks. In this study, we use the same time interval $\pm$ 0.1s in our
model to evaluate the model’s generalization on the INSTANCE dataset.

Figure 7 (a) shows the frequency-magnitude distribution of the used subset of
INSTANCE data, and (b) gives an example of one seismic event signal of
different lengths, where dashed lines represent the length of the input fed
into different models. 20s data is used in ARRU [19] and 30s data is used in
PhaseNet [11]111Here, the used PhaseNet model [11] and ARRU model [19] are
directly cited from their saved trained models without training in this paper,
respectively.. Hence, for this comparison, only the seismic data whose P and S
arrival times are within 20s are considered.

TABLE IV: Experiment results on INSTANCE [34].

|  | Mean (s) | St. Dev. (s)
---|---|---|---
P phase | EPick | -0.885 | 3.947
| PhaseNet[11] | -0.188 | 2.075
| ARRU[19] | 0.354 | 1.181
S phase | EPick | 0.116 | 1.895
| PhaseNet[11] | -0.178 | 2.234
| ARRU[19] | 0.319 | 3.130

  * •

Mean and standard deviation (st. dev.) of the differences of model predicted
arrivals minus manual picked arrivals in seconds

The model comparison are shown in Table IV. The results demonstrate that
compared to PhaseNet [11] and ARRU [19], our model achieves better performance
for first S-phase arrival picking on the INSTANCE dataset. On the first
P-phase arrival picking, the error is comparatively large. The reason might be
that a large weight is assigned to the first S-phase arrival picking during
the model training. We also observed that by using the similar weights for
both first P-phase and S-phase, the performance of S-phase reduced without any
significant improvement in P-phase picking error. Additional reason might be
that PhaseNet and ARRU were trained on data labeled in window format, while
our model used single sample labeled data. Furthermore, only earthquake
signals are used in the training process of PhaseNet and ARRU. Meanwhile, our
proposed model is also trained on non-earthquake signals, which leads to a
more robust performance in the case of noisy seismic data.

##  VI Conclusion

In this study, we investigate the combination of U-shaped neural network and
attention mechanism involving self-attention and multi-head attention for
seismic phase picking. To fully leverage the power of attention mechanism,
EPick is proposed, which not only completes the task of seismic event
detection, but also well utilizes the low-level extracted features by using
the UNet architectural design to achieve phase arrival time detection. As an
alternative framework for seismic phase picking, EPick achieves superior
performance than the baseline method for first S-phase arrival picking,
whereas the performance is more robust in the case of P-phase. The
experimental results well-demonstrate the generalization ability of the
proposed model in S-phase picking. This model can be used in tasks that
require fast seismic data processing, as well as dealing with big data. EPick
can further be developed in monitoring real-time seismic signals.

## Acknowledgment

This research is supported by the “KI-Nachwuchswissenschaftlerinnen” - grant
SAI 01IS20059 by the Bundesministerium für Bildung und Forschung - BMBF.
Calculations were performed at the Frankfurt Institute for Advanced Studies’
new GPU cluster, funded by BMBF for the project Seismologie und Artifizielle
Intelligenz (SAI). We also thank Dr Jan Steinheimer and Jonas Köhler for their
kind suggestions.

## References

  * [1] D. Seidl and W. Stammler, “Restoration of broad-band seismograms. (part 1),” _Geophysics_ , pp. 114–122, 1984.
  * [2] R. V. Allen, “Automatic earthquake recognition and timing from single traces,” _Bulletin of the Seismological Society of America_ , vol. 68, no. 5, pp. 1521–1532, 1978.
  * [3] R. Sleeman and T. Van Eck, “Robust automatic p-phase picking: an on-line implementation in the analysis of broadband seismogram recordings,” _Physics of the earth and planetary interiors_ , vol. 113, no. 1-4, pp. 265–275, 1999.
  * [4] Z. Ross, M. White, F. Vernon, and Y. Ben-Zion, “An improved algorithm for real-time s-wave picking with application to the (augmented) anza network in southern california,” _Bulletin of the Seismological Society of America_ , vol. 106, no. 5, pp. 2013–2022, 2016.
  * [5] J. Wang, Z. Xiao, C. Liu, D. Zhao, and Z. Yao, “Deep learning for picking seismic arrival times,” _Journal of Geophysical Research: Solid Earth_ , vol. 124, no. 7, pp. 6612–6624, 2019.
  * [6] G. E. Hinton, S. Osindero, and Y.-W. Teh, “A fast learning algorithm for deep belief nets,” _Neural computation_ , vol. 18, no. 7, pp. 1527–1554, 2006\.
  * [7] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [8] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot _et al._ , “Mastering the game of go with deep neural networks and tree search,” _nature_ , vol. 529, no. 7587, pp. 484–489, 2016.
  * [9] T. Perol, M. Gharbi, and M. Denolle, “Convolutional neural network for earthquake detection and location,” _Science Advances_ , vol. 4, no. 2, p. e1700578, 2018.
  * [10] Z. E. Ross, M.-A. Meier, and E. Hauksson, “P wave arrival picking and first-motion polarity determination with deep learning,” _Journal of Geophysical Research: Solid Earth_ , vol. 123, no. 6, pp. 5120–5129, 2018.
  * [11] W. Zhu and G. C. Beroza, “Phasenet: a deep-neural-network-based seismic arrival-time picking method,” _Geophysical Journal International_ , vol. 216, no. 1, pp. 261–273, 2019.
  * [12] Y. Zhou, H. Yue, Q. Kong, and S. Zhou, “Hybrid event detection and phase-picking algorithm using convolutional and recurrent neural networks,” _Seismological Research Letters_ , vol. 90, no. 3, pp. 1079–1087, 2019.
  * [13] E. Pardo, C. Garfias, and N. Malpica, “Seismic phase picking using convolutional networks,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 57, no. 9, pp. 7086–7092, 2019.
  * [14] L. Zhu, Z. Peng, J. McClellan, C. Li, D. Yao, Z. Li, and L. Fang, “Deep learning for seismic phase detection and picking in the aftershock zone of 2008 mw7. 9 wenchuan earthquake,” _Physics of the Earth and Planetary Interiors_ , vol. 293, p. 106261, 2019.
  * [15] Y. Chen, “Automatic microseismic event picking via unsupervised machine learning,” _Geophysical Journal International_ , vol. 222, no. 3, pp. 1750–1764, 2020.
  * [16] S. M. Mousavi, W. L. Ellsworth, W. Zhu, L. Y. Chuang, and G. C. Beroza, “Earthquake transformer–an attentive deep-learning model for simultaneous earthquake detection and phase picking,” _Nature communications_ , vol. 11, no. 1, pp. 1–12, 2020.
  * [17] C. Chai, M. Maceira, H. J. Santos-Villalobos, S. V. Venkatakrishnan, M. Schoenball, W. Zhu, G. C. Beroza, C. Thurber, and E. C. Team, “Using a deep neural network and transfer learning to bridge scales for seismic phase picking,” _Geophysical Research Letters_ , vol. 47, no. 16, p. e2020GL088651, 2020.
  * [18] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” _Advances in neural information processing systems_ , vol. 25, pp. 1097–1105, 2012.
  * [19] W.-Y. Liao, E.-J. Lee, D. Mu, P. Chen, and R.-J. Rau, “Arru phase picker: Attention recurrent-residual u-net for picking seismic p-and s-phase arrivals,” _Seismological Research Letters_ , 2021.
  * [20] S. M. Mousavi, Y. Sheng, W. Zhu, and G. C. Beroza, “Stanford earthquake dataset (stead): A global data set of seismic signals for ai,” _IEEE Access_ , vol. 7, pp. 179 464–179 476, 2019.
  * [21] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [22] M. Zhao, S. Chen, L. Fang, and A. Y. David, “Earthquake phase arrival auto-picking based on u-shaped convolutional neural network,” _Chinese J. Geophys_ , vol. 62, no. 8, pp. 3034–3042, 2019.
  * [23] J. Zhang and G. Sheng, “First arrival picking of microseismic signals based on nested u-net and wasserstein generative adversarial network,” _Journal of Petroleum Science and Engineering_ , vol. 195, p. 107527, 2020.
  * [24] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” _arXiv preprint arXiv:1706.03762_ , 2017.
  * [25] K. Han, Y. Wang, H. Chen, X. Chen, J. Guo, Z. Liu, Y. Tang, A. Xiao, C. Xu, Y. Xu _et al._ , “A survey on visual transformer,” _arXiv preprint arXiv:2012.12556_ , 2020.
  * [26] H. Noh, S. Hong, and B. Han, “Learning deconvolution network for semantic segmentation,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1520–1528.
  * [27] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016, http://www.deeplearningbook.org.
  * [28] K. P. Murphy, _Machine learning: a probabilistic perspective_. MIT press, 2012.
  * [29] Z. Xie, I. Sato, and M. Sugiyama, “Stable weight decay regularization,” _arXiv preprint arXiv:2011.11152_ , 2020.
  * [30] C. Cortes, M. Mohri, and A. Rostamizadeh, “L2 regularization for learning kernels,” _arXiv preprint arXiv:1205.2653_ , 2012.
  * [31] H. Kaur, H. S. Pannu, and A. K. Malhi, “A systematic review on imbalanced data challenges in machine learning: Applications and solutions,” _ACM Computing Surveys (CSUR)_ , vol. 52, no. 4, pp. 1–36, 2019.
  * [32] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/
  * [33] T. Diehl, N. Deichmann, E. Kissling, and S. Husen, “Automatic s-wave picker for local earthquake tomography,” _Bulletin of the Seismological Society of America_ , vol. 99, no. 3, pp. 1906–1920, 2009.
  * [34] A. Michelini, S. Cianetti, S. Gaviano, C. Giunchi, D. Jozinovic, and V. Lauciani, “Instance–the italian seismic dataset for machine learning,” _Earth System Science Data Discussions_ , pp. 1–47, 2021.
  * [35] C. Sammut and G. I. Webb, _Encyclopedia of machine learning_. Springer Science & Business Media, 2011.