Machine Learning (2021) 110:1585–1607
https://doi.org/10.1007/s10994-021-05991-y

Learning subtree pattern importance for Weisfeiler‑Lehman
based graph kernels
Dai Hai Nguyen1,2 · Canh Hao Nguyen2 · Hiroshi Mamitsuka2,3Received: 22 November 2020 / Revised: 23 March 2021 / Accepted: 3 May 2021 /Published online: 13 June 2021
© The Author(s), under exclusive licence to Springer Science+Business Media LLC, part of Springer Nature 2021

Abstract
Graph is an usual representation of relational data, which are ubiquitous in many domains
such as molecules, biological and social networks. A popular approach to learning withgraph structured data is to make use of graph kernels, which measure the similarity
between graphs and are plugged into a kernel machine such as a support vector machine.
Weisfeiler-Lehman (WL) based graph kernels, which employ WL labeling scheme to
extract subtree patterns and perform node embedding, are demonstrated to achieve great
performance while being efficiently computable. However, one of the main drawbacks of ageneral kernel is the decoupling of kernel construction and learning process. For moleculargraphs, usual kernels such as WL subtree, based on substructures of the molecules, consider all available substructures having the same importance, which might not be suitablein practice. In this paper, we propose a method to learn the weights of subtree patterns in
the framework of WWL kernels, the state of the art method for graph classification task
(Togninalli et al., in: Advances in Neural Information Processing Systems, pp. 6439–6449,
2019). To overcome the computational issue on large scale data sets, we present an efficient learning algorithm and also derive a generalization gap bound to show its convergence. Finally, through experiments on synthetic and real-world data sets, we demonstrate
the effectiveness of our proposed method for learning the weights of subtree patterns.
Keywords Graph kernel · Optimal transport · Weisfeiler Lehman scheme

Editors: Annalisa Appice, Sergio Escalera, Jose A. Gamez, Heike Trautmann.
* Dai Hai Nguyen
hai@k.u-tokyo.ac.jp
Canh Hao Nguyen
canhhao@kuicr.kyoto-u.ac.jp
Hiroshi Mamitsuka
mami@kuicr.kyoto-u.ac.jp
1

Graduate School of Frontier Sciences, The University of Tokyo, 5‑1‑5 Kashiwa‑no‑ha, Kashiwa,
Chiba 277‑8561, Japan

2

Bioinformatics Center, Institute for Chemical Research, Kyoto University, Uji 611‑0011, Japan3

Department of Computer Science, Alato University, Espoo 02150, Finland

13

Vol.:(0123456789)

1586

Machine Learning (2021) 110:1585–1607

1 Introduction
Graphs are natural data structures, which appear in various domains such as bioinformatics (Sharan and Ideker 2006), cheminformatics (Trinajstic 2018), social network analysis (Scott 2011) and so on, where nodes (vertices) represent objects and edges represent
the relations between them. A popular approach to learning with graph structured data is
to make use of graph kernels. Essentially, a graph kernel is a measure of the similarity
between two graphs and must satisfy two fundamental requirements of being a valid kernel: (1) symmetric and (2) positive semi-definite (PSD). Furthermore, the requirements ofdesigning a graph kernel are: it should capture the semantic inherent in the graph structures(e.g. substructures of different levels), and it must be efficiently computable (Vishwanathan
et al. 2010).
A number of graph kernels have been proposed in literature such as random walk
(Kashima 2003), shortest path (Borgwardt and Kriegel 2005), Weisfeiler-Lehman (WL)
subtree (Shervashidze and Borgwardt 2009) kernels, just to name a few. Most of them arebased on R-Convolution framework (Haussler 1999), which decomposes two graphs into
substructures and adds up the similarities between their substructures to compute kernel
values. Different graph kernels are defined under different ways of decomposition (typesof substructures). For instance, the substructures can be random walks (Kashima 2003),
shortest paths (Borgwardt and Kriegel 2005) or subtree patterns (Vishwanathan et al.
2010). Among these, WL subtree kernels have been shown to achieve great prediction performance while being efficiently computable. The key point is that it simply employs a WL
based color refinement scheme to embed each node in a given graph into a vector of WL
labels, which correspond to subtree patterns of the graphs. Then, the kernel between two
graphs is defined as the sum of all pairwise similarities between any two node embeddingsof the two graphs.
Following a different approach, WL based optimal assignment (WL-OA) kernel (Kriege
et al. 2016) assigns one node embedding of one graph to one embedding of the other suchthat the total similarities between assigned node embeddings is maximized. This is also
known as optimal assignment problem in combinatorial mathematics (Munkres 1957). In asimilar vein, Wasserstein WL (WWL, Togninalli et al. 2019) uses optimal transport (OT),also known as Wasserstein distance (Villani 2008), for measuring the distance between two
graphs based on their WL node embeddings. The distance is then converted into a similarity matrix through Laplacian kernel. Furthermore, both of these similarity matrices areshown to be valid kernels due to the hierachy property of WL labels (see Kriege et al. 2016
and Togninalli et al. 2019 for more details).
One of the main drawbacks of these kernels is that they are predefined feature extraction without learning the importance of substructures to the problem. This results in thedecoupling of data representation and learning process. In these kernels, substructures
are given the same weights. However, for the problems such as molecule classification, it
is known that only subparts of the molecules are responsible for their properties. Therefore, we wish to be able to give weights to their substructures to have higher classification
performance and model interpretation. Based on this motivation, we propose a model tolearn the weights of subtree patterns (extracted by WL labeling scheme). Our work extends
WWL kernels (Togninalli et al. 2019) by formulating an OT based distance as a parametricfunction of subtree pattern weights before converting into kernels. We also propose an efficient stochastic learning algorithm to estimate the weights and derive a generalization gap
bound for the algorithm. Finally, through experiments on synthetic and four real-world data13

Machine Learning (2021) 110:1585–1607

1587

sets, we show that learning important subtree patterns by our proposed method can lead tomore accurate predictive performance and extract important patterns which enhance theclassification results.
The remainder of the paper is organized as follows: in Sect. 2, we review graph kernels which are based on WL labeling scheme, including WL subtree, WL-OA and WWL
kernels. In Sect. 3, we present our method that parameterizes the Wasserstein distance
between two graphs with WL labeling scheme as a function of subtree patterns and present the stochastic algorithm for learning parameters of the function. In Sect. 4, we derivea generalization bound for the learning algorithm. In Sect. 5, experimental results on thesynthetic and real-world data sets are provided. Finally, we conclude by summarizing thiswork and discussing possible extensions in Sect. 6.

2 Related work
In this paper, we consider the binary classification problem for graph structured data:
given a collection of labeled graphs (gi , yi ), i = 1, .., n (where n is the number of examples) drawn from an unknown joint distribution P over G × {−1, 1}, where G is a space ofgraphs. We wish to learn a classifier h ∶ G → {−1, 1}, which is based on a similarity function K ∶ G × G → [−1, 1]. If K is symmetric and positive semi-definite (PSD), it is called avalid kernel.
There are a number of proposed kernels on graphs, see (Kashima 2003; Borgwardt
and Kriegel 2005; Shervashidze and Borgwardt 2009). In general, they are defined basedon R-Convolution framework (Haussler 1999), that is, each graph g ∈ G is decomposed
into substructures, and a kernel value K(g, g� ) is defined as a sum of pairwise similarities
between their substructures. In fact, many graph kernels can be considered as instances of
the R-Convolution framework under different decomposition into substructures. The substructures can be random walks (Kashima 2003), shortest paths (Borgwardt and Kriegel2005) or circle subtrees (Shervashidze and Borgwardt 2009). Among these, WeisfeilerLehman (WL) subtree kernels (Shervashidze and Borgwardt 2009) and its variants havebeen shown to achieve great performance for the graph classification tasks. In this work,
we focus on WL based kernels and will review them in the following subsections.

2.1 Weisfeiler‑Lehman (WL) scheme for node embeddings
Weisfeiler-Lehman (WL) subtree kernels (Shervashidze and Borgwardt 2009) are based on
an iterative colour refinement (also known as WL labeling scheme) and have been shownto achieve great performance for graph classification task. For each node of a given graph,
the WL labeling scheme creates a sequence of ordered strings by the aggregation of thelabels of the node and its neighbors; these strings are then hashed or indexed to produce
compressed updated node labels or new indices. If the iteration of the scheme is increased,these obtained labels represent increasingly broader neighborhood of each node. More specifically, for a graph G = (V, E) with initial labels 𝓁0 (v) for v ∈ V and let H be the number
of WL iterations, we can define a sequence of refined labels (𝓁0 , 𝓁1 , ..., 𝓁H ), where 𝓁h+1 isobtained from 𝓁h by the following procedure: 𝓁h+1 (v) = hash(𝓁h (v), Nh (v)), where Nh (v)denotes a lexicographically sorted sequence of labels of v’s neighbors at iteration h and thehash function is to create a updated compressed node label for v. We use perfect hashing

13

1588

Machine Learning (2021) 110:1585–1607

Fig. 1  a Illustration of Weisfeiler-Lehman (WL) node embedding scheme with two iterations ( H = 2). b
Illustration of different WL based graph kernels: WL Subtree kernel is computed by considering all pairwise similarities between node embeddings of two graphs; WL optimal assignment kernel takes the similarities of only aligned node embeddings into account; For Wasserstein WL kernel, one node embedding of
one graph can be coupled with multiple node embeddings of the other graph for computing the kernel

for the hash function, as in Shervashidze and Borgwardt (2009), ensuring two nodes at
iteration h + 1 have the same label if and only if their label and those of their neighbors at
iteration h are the same. Throughout the rest of paper, we denote the set of WL labels at
Σh. The WL labeiteration h as Σh, for h = 1, ..., H , and the set of all WL labels as Σ = ∪H
h=1ling scheme is illustrated in Fig. 1a.
Based on WL labeling scheme, a WL node embedding scheme was proposed to generate node embeddings from node labels of the graphs (Shervashidze and Borgwardt 2009).

Definition 1 (WL based node embedding scheme, Shervashidze and Borgwardt 2009). Let
G = (V, E) and let H be the number of WL iterations. For every h ∈ {1, ..., H}, we define
the node embedding xh (v) of a node v ∈ V and graph embedding 𝐗h (G) of G at iteration has follows:

13

1589

Machine Learning (2021) 110:1585–1607

xh (v) = 𝓁 h (v), 𝐗h (G) = [xh (v1 ), ..., xh (vnV )]T(1)

Then the graph embedding of G can be defined as:f H (G) ∶ G ↦ ΣnV ×H]T[G ↦ 𝐗1 (G)T , ..., 𝐗H (G)T
where nV is the number of nodes in G. With the WL node embedding scheme above, weare ready to introduce some notations which will be used throughout the rest of paper.
(f h (G), f h (G� )) be the Hamming distance matrix where each entry is
Notations Let DHamhthe normalized Hamming distance between the corresponding node embeddings of G andG′ at iteration h, defined as:

1 ∑ discd (u, v), didisc (u, v) =h i=1 ihdhHam (u, v) ={0 if ui = vi
1 otherwise(2)

where u and v denote two node embeddings, dhHam (u, v) denotes the normalized Hamming
distance between u and v at iteration h and didisc (u, v) denotes the discrete distance betweenu and v at iteration i. Similarly, let DDisc(𝐗h (G), 𝐗h (G� )) be the discrete distance matrixhwhere each entry is the discrete distance between the corresponding node embeddings of]ij ∈ [0, 1] and [DDisc]ij ∈ {0, 1}. We also
G and G′ at iteration h. It is easy to see that [DHamhhdefine a base kernel which corresponds to the averaged number of feature shared by twonode embeddings as:

1∑
𝟙(ui = vi )h i=1hkh (u, v) =(3)

It is easy to see that kh (u, v) = 1 − dhHam (u, v). Thanks to WL node embedding scheme, agraph can be represented as a point cloud or a set of node embeddings. Measuring thesimilarity (or dissimilarity) between two graphs boils down to measuring the similarity (ordissimilarity) between two sets of embeddings. Each node embedding captures information
about the neighborhood of the corresponding node (or a rooted subtree pattern). In the following subsections, we will review different WL based graph kernels derived from different ways of comparing their sets of WL node embeddings.

2.2 WL subtree kernels
WL subtree kernels (Shervashidze and Borgwardt 2009) simply employ the aforementioned WL labeling scheme to extract subtree patterns, which represent the neighborhood of each node in the graph up to a given distance (or number of hops H). Essentially
WL subtree kernel counts the number of common WL labels. In the context of WL node
embedding scheme, it can be computed by summing all pairwise similarities between node
embeddings of two graphs. More formally, for two graphs G and G′ with two sets of node
embeddings f H (G) and f H (G� ), respectively, the WL subtree kernel value between themcan be defined as:∑∑
kH (𝐱, 𝐲)
𝐊WL (G, G� ) =
(4)
H
H
�
𝐱∈f (G) 𝐲∈f (G )

13

1590

Machine Learning (2021) 110:1585–1607

It is obvious that as the base kernel kH (𝐱, 𝐲) (defined in Eq. (3)) is equal to the number of
WL labels (subtree patterns) shared by two node embeddings, the kernel 𝐊WL (G, G� ) isequal to the total number of WL subtree patterns shared by two graphs.

2.3 WL‑based optimal assignment kernels
Optimal assignments are natural measures of similarity between two sets of points. In particular, for two sets of points, the goal is to assign one point of one set to another pointof the other set (one-to-one correspondence) such that the sum of similarities between
assigned points is maximized. Finding such an optimal alignment or bijection is also
known as the well-studied assignment problem in combinatorial optimization (Munkres
1957). However, a challenge is how to design a valid kernel based on optimal assignments.
Kriege et al. (2016) introduced a restricted class of kernels, called strong kernels, thatguarantees the construction of valid optimal assignment based kernels. An important result
is that the strong kernels give rise to hierarchies defined on the domain of kernels. Basedon this, the authors proposed WL-based optimal assignment (WL-OA) kernels with the
WL node embedding scheme. WL-OA kernels employ the base kernel, defined in Eq. (3),which satisfies the requirement of being a strong kernel as the sequence of refined WLlabels (𝓁0 , 𝓁1 , ..., 𝓁H ) gives rise to a family of nested subsets, which can be represented bya hierarchy. Consequently, WL-OA kernels are valid. Formally, given two graphs G andG′ with two sets of node embeddings f H (G) and f H (G� ), respectively, the WL-OA kernel
value between them is defined as:∑
𝐊OA (G, G� ) =max
kH (𝐱, 𝐲)(5)
B∈B(f H (G),f H (G� ))(x,y)∈B

where B(f H (G), f H (G� )) is the set of all bijections between two sets f H (G) and f H (G� ). To
apply this kernel to graphs of different number of nodes, we can fill up the graph withsmaller number of nodes, says f H (G� ), by new node embeddings z with k(𝐱, 𝐳) = 0 for all𝐱 ∈ f H (G) without changing the result. It is worth noting that the WL-OA kernels take thesimilarities of aligned node embeddings into account, while the WL subtree kernels consider all pairwise similarities.

2.4 Wasserstein WL kernels
Optimal transport (OT), also known as Wasserstein distance function (Villani 2008), hasgained much attraction in machine learning community as a powerful tool for the comparsion of two probability distributions. The naive computation of this distance between two
discrete measures, e.g. point clouds, involves solving transport problem. Formally, let
X = {𝐱1 , ..., 𝐱m } and Y = {𝐲1 , ..., 𝐲n } be two sets of points, where m and n denote the size of
two sets X and Y, respectively; 𝐩 ∈ ℝm and 𝐪 ∈ ℝn are two discrete probability distributions over X and Y, respectively. We use dij to denote the distance between 𝐱i and 𝐲j , e.g. the‖‖2squared Euclidean distance dij = ‖𝐱i − 𝐲j ‖ . The Wasserstein distance is formulated as a‖‖2linear program over the transportation matrix (or joint probability) 𝐏 ∈ ℝm×n:

13

1591

Machine Learning (2021) 110:1585–1607

W1 (X, Y, d) = min𝐏ij

subject to

m
∑

m n
∑
∑

𝐏ij dij

i=1 j=1

𝐏ij = 𝐪j , ∀j ∈ [1, n]

i=1
n

∑

(6)

𝐏ij = 𝐩i , ∀i ∈ [1, m]

j=1With WL node embedding scheme to generate node embeddings for graphs, Togninalli
et al. (2019) evaluated the pairwise Wasserstein distance between graphs with the normalised Hamming (2) as the ground distance. Then, Wasserstein WL (WWL) kernel is defined
as an instance of Laplacian kernels, see Eq. (7).

𝐃WWL (G, G� ) = W1 (f H (G), f H (G� ), dHHam )𝐊WWL (G, G� ) = e−𝛾𝐃

WWL (G,G� )(7)where 𝛾 is a hyperparameter.
In general cases, it is not necessarily possible to derive a valid kernel from the Wasserstein distance. However, thanks to the special property of WL labels and normalized Hamming distance (2), 𝐃WWL was shown to be conditionally negative definite (CND), resulting
in the validity of 𝐊WWL , by proving the following lemma (see Togninalli et al. 2019 for its
proof):
Lemma 1 If a transportation matrix PH is optimal solution of (6) with the ground distancedHHam (2) between node embeddings at iteration H, then we have the two following claims:
1. PH is also optimal solution of (6) with the discrete distance dHDisc between H-iteration
values.
Ham
between
2. PH is also optimal solution of (6) with the normalised Hamming distance dH−1
node embeddings at iteration H − 1.
Let P∗ be the optimal solution of (6) for DHam
. From the above lemma, it is also the optiH,.
The
Wasserstein
distance between two graphs G and G′
h
=
1,
..,
H
mal solution for DDischin (7) can be simplified as follows:

1 ∑W (𝐗h (G), 𝐗h (G� ), dhDisc )H h=1 1H

𝐃WWL (G, G� ) =(8)

The Eq. (8) is a sum of OT distances with the discrete distances as ground metrics, which
are CND. Therefore, the sum is also CND, leading to the validity of the similarity matrix
𝐊WWL.

13

Machine Learning (2021) 110:1585–1607

1592

3 Incorporating subtree pattern importance into WL based graphkernels
One of the main limitations of kernels is the decoupling of data representation and learning
process, that is, the kernel must be predefined prior to learning, leading to limited predictive performance. Furthermore, in prediction tasks for molecular data, the output might bedetermined by the presence of a few important substructures, while these kernels containall substructures with equal weights. Motivated by this drawback, in this paper we address
the problem of incorporating subtree pattern weights for WWL kernel (Togninalli et al.
2019). To this end, we aim to learn new kernels from a parametric form of Wasserstein
distance taking into account subtree pattern weights (8), and learn these weights from data
optimally for the task.

3.1 Parametric form of Wasserstein distance with subtree pattern weights
To derive a parametric form of the distance function (8), we rely on the following simple
observation:
Lemma 2 Let S be a set of elements, X = {𝐱1 , ..., 𝐱m } and Y = {𝐲1 , ..., 𝐲n } ( X, Y ⊆ S ) be two
multiset of S of m and n samples, respectively, the Wasserstein distance between them withthe discrete distance as the ground metric is determined by:∑min(𝜇X (v), 𝜇Y (v))
W1 (X, Y) = 1 −(9)v∈Swhere 𝜇X (v) denotes the mass density function of the multiset X with v ∈ S.
Applying this lemma to Eq. (8), we have:1 ∑∑min(𝜇𝐗h (G) (v), 𝜇𝐗h (G� ) (v))
H h=1 v∈Σh
H

𝐃WWL (G, G� ) = 1 −(10)

Our idea is to give each substructure or WL label v a nonnegative weight wv ∈ ℝ≥0 for its
importance to the problem, so the parametric form of Eq. (10) is defined as follows:��1 �⟨𝐰h , 𝐳h G, G� ⟩
H h=1
H

b−

(11)

(
)
h
where 𝐰h, 𝐳h G, G� ∈ ℝ|Σ | are the vectors of entries wv and min(𝜇𝐗h (G) (v), 𝜇𝐗h (G� ) (v)),
h ; b is a constant to ensure that the value of parametric function isrespectively, for v ∈ Σnonnegative. In vector form, this can be expressed as:��
d𝐖 (G, G� ) = b − ⟨𝐖, 𝐙 G, G� ⟩
(12)
� �
⎧ �
�
�T �T
�T
�
⎪ 𝐙 G, G� = 1 𝐳1 G, G� , ..., 𝐳H G, G�
H
where ⎨�
�T
⎪ 𝐖 = 𝐰T1 , ..., 𝐰TH⎩
The parametric form (12) is a linear
functionwith respect to the parameter vector)(
𝐖 ∈ ℝd (d = |Σ1 | + ... + |ΣH |) and 𝐙 G, G′ is considered as a feature vector of a pair of

13

1593

Machine Learning (2021) 110:1585–1607

Fig. 2  Illustration of Hinge loss max(0, 𝛼1 − x) and max(0, x − 𝛼2 ); and their smooth versions: V1 (left) andV2 (right), respectively, with 𝛼1 = 3, 𝛼2 = 4 and 𝜎 = 1.5

graphs G and G′. Once the parameters are estimated, we can derive a similarity matrix
through the Laplacian kernel as in Eq. (7). More importantly, as 𝐖 is nonnegative, it iseasy to see that d𝐖 is a CND function, and thus the derived similarity matrix is valid.

3.2 Formulation of learning subtree pattern weights 𝐖
We aim to learn the parameters W in Eq. (12) using the notions of metric learning (Kulis2012). That is two input graphs with the same labels are encouraged to be closer while thetwo with different labels become far away from each other. In other words, within classdistances should be small, while between class distances should be large. For this purpose,as a loss function for a graph pair g and g′, we can use the following two hinge loss function: max(0, 𝛼1 − d𝐖 (g, g� )) if g and g′ are with different labels and max(0, d𝐖 (g, g� ) − 𝛼2 )otherwise, for learning subtree pattern weights, where 𝛼1 and 𝛼2 are constants (𝛼1 ≥ 𝛼2).
The former yields a penalty if g and g′ of different labels are closer than 𝛼1 while the latteryields a penalty when g and g′ of the same label are more distant than 𝛼2. Instead of using
these functions in the optimization problem, we use their smooth versions: V1 and V2 (seeFig. 2), which offer useful properties for deriving a generalization bound for the problem
in Sect. 4. The derivation of these functions is based on the connection between the strong
convexity of a function and Lipschitz continuous gradient of its Fenchel dual (see moredetails in Nesterov (2005)).
More concretely, let Dn = {z1 = (g1 , y1 ), ..., zn = (gn , yn )} where gi ∈ G andyi ∈ Y = {−1, 1}, for i = 1, ..., n, we formulate a constrained minimization problem forlearning subtree pattern weights as follows:

minimize𝐖

1 ∑∑𝓁(𝐖, zi , zj )n2 i=1 j=1
n

n

(13)

subject to
𝐖∈C
]
[
T
where C = {𝐖 ∈ ℝd ∶ 𝐖 = 𝐰T1 , ..., 𝐰TH , ||𝐰h − 𝐜h ||2 ≤ 𝜖h , 1 ≤ h ≤ H} (𝐜h and 𝜖h areconstant vectors and{scalars); 𝓁 is a continuously differentiable function and defined as folV1 (𝐖, gi , gj ) if yi ≠ yj
lows: 𝓁(𝐖, zi , zj ) =V2 (𝐖, gi , gj ) otherwise13

Machine Learning (2021) 110:1585–1607

1594

V1 (𝐖, gi , gj )⎧0
if d𝐖 (gi , gj ) ≥ 𝛼1
⎪ 𝛼 − 𝜎 − d (g , g ) if d (g , g ) ≤ 𝛼 − 𝜎=⎨ 1 2
𝐖 i j
𝐖 i j
1
�
�
⎪ 1 d𝐖 (gi , gj ) − 𝛼1 2 if 𝛼1 − 𝜎 < d𝐖 (gi , gj ) < 𝛼1
⎩ 2𝜎

(14)

V2 (𝐖, gi , gj )⎧0
if d𝐖 (gi , gj ) ≤ 𝛼2
⎪ d (g , g ) − 𝛼 − 𝜎 if d (g , g ) ≥ 𝛼 + 𝜎=⎨ 𝐖 i j
2
𝐖 i j
2
2�
�
⎪ 1 d𝐖 (gi , gj ) − 𝛼2 2 if 𝛼2 < d𝐖 (gi , gj ) < 𝛼2 + 𝜎
⎩ 2𝜎

(15)

)
(
(
)
and 𝛼1, 𝛼2 𝛼1 ≥ 𝛼2 and 𝜎 are constants; d𝐖 (gi , g(j ) is )calculated from 𝐙 gi , gj as in Eq.
(12). To reduce notation, we use 𝐙i,j rather than 𝐙 gi , gj in the rest of the paper.
2‖ ‖Lemma 3 Let 𝛽 = max1≤i<j≤n ‖𝐙ij ‖ , L = 𝛽𝜎 , and M = max(b − 𝜎2 − 𝛼2 , 𝛼1 − 𝜎2 ), the loss
‖ ‖2function 𝓁 defined in problem (13) is L-lipschitz, 𝛽-smooth and upper bounded by M.

Proof First, it is obvious to see that as 0 ≤ d𝐖 ≤ b, we have the following bounds:
0 ≤ V1 ≤ 𝛼1 − 𝜎2 and 0 ≤ V2 ≤ b − 𝜎2 − 𝛼2. Therefore, 𝓁 is upper bounded by M.
We derive the first and second order derivatives of the smooth function V1 as follows:∇V1 (𝐖, gi , gj )⎧0
if d𝐖 (gi , gj ) ≥ 𝛼1⎪
if d𝐖 (gi , gj ) ≤ 𝛼1 − 𝜎= ⎨ 𝐙�ij�⎪ 1 𝛼1 − d𝐖 (gi , gj ) 𝐙ij if 𝛼1 − 𝜎 < d𝐖 (gi , gj ) < 𝛼1
⎩𝜎

(16)

∇2 V1 (𝐖, gi , gj )⎧0
if d𝐖 (gi , gj ) ≥ 𝛼1⎪
0
if d𝐖 (gi , gj ) ≤ 𝛼1 − 𝜎=⎨⎪ 1 𝐙ij 𝐙Tij if 𝛼1 − 𝜎 < d𝐖 (gi , gj ) < 𝛼1
⎩𝜎

(17)

In order to prove the function V1 is L-Lipschitz and 𝛽-smooth, it is sufficient to show that‖‖the norm of its derivative is always less than L: ‖∇V1 (𝐖, gi , gj )‖ ≤ L and the spectral‖‖2
norm (or the maximum eigen value) of its second order derivative is always less than 𝛽 :||∇2 V1 (𝐖, gi , gj )|| ≤ 𝛽, ∀gi , gj ∈ G . Indeed, from Eqs. (16) and (17), we have:‖‖ 1‖‖ 1 ‖ ‖2‖‖‖ ‖‖V1 (𝐖, gi , gj )‖ ≤ ‖𝐙ij ‖ and ‖∇2 V1 (𝐖, gi , gj )‖ ≤ 𝜎 ‖𝐙ij 𝐙Tij ‖ = 𝜎 ‖𝐙ij ‖ . Also, we can
‖
‖
‖
‖
‖
‖
‖
‖
‖ ‖2
2
2
‖ ‖
‖ ‖
bound ‖𝐙ij ‖ by the inequality ‖𝐙ij ‖ ≤ 𝛽 . Thus V1 is 𝛽-smooth and L-Lipschitz. Similarly,
‖ ‖2
‖ ‖2
we can also show that V2 is 𝛽-smooth and L-Lipschitz. The lemma is proven.□

13

1595

Machine Learning (2021) 110:1585–1607

3.3 A stochastic learning algorithm for constrained optimization
The constrained optimization problem (13) is convex and thus guarantees to find its global
optimum. Standard methods such as projected gradient descent can be used to solve theproblem (13). However, for large scale data sets, solving the problem (13), involving n2
terms with d parameters, might be computationally expensive. For instance, the data set
PROTEIN (see Table 1) has n2 > 106 pairs of examples and the weight vector size ofd > 105 with the number of WL iterations H = 5. In this subsection, we present an efficientstochastic learning algorithm for dealing with this issue.

Let 𝐖(t) denote the weight at iteration t. The weight is initialized by a vector of ones:
𝐖(0) = 𝟏d , which is also the case of WWL kernel without learning subtree pattern importance. At each iteration t, we randomly pick up a pair of examples ( z = (g, y), z� = (g� , y� ))from the training data set Dn and compute the gradient grad(t) corresponding to( this )pair. In
fact this step can be done efficiently due to the sparsity of the feature vector 𝐙 g, g′ in Eq.
(12). Then, we update the current solution 𝐖(t) to 𝐖(t+1) by the following rule:{
where projC [𝐖]h =𝐖(t+1) = projC [𝐖(t) − 𝜇grad(t) ]if ||𝐰h − 𝐜h ||2 ≤ 𝜖h
𝜖h
(𝐰
−
𝐜
)
+
𝐜
otherwisehhh||𝐰h −𝐜h ||2𝐰hmaps a point 𝐰h(1 ≤ h ≤ H ) back to the bounded feasible region. The procedure is illustrated in
Algorithm 1.

4 Theoretical guarantees: a bound on generalization gap
In this section, we provide a bound on the generalization gap of the proposed stochastic
learning algorithm for solving the problem (13). The gap is defined as the expected difference between the generalization error R(.) and empirical error RDn (⋅). In order to derive thegeneralization bound, we first provide basic setup and notations; then prove that our learning algorithm has a uniform stability, which is established in Theorem 5 using Lemma 3

13

Machine Learning (2021) 110:1585–1607

1596

and Lemma 4; finally derive our generalization bound, which is established in Theorem 9
using the McDiarmid inequality (see Theorem 6).

4.1 Basic setup and notations
Generalization Error. Let 𝐖n be the parameters of the parametric function (12) obtainedby training on the data set Dn using Algorithm 1. Then the generalization error (or risk)R(𝐖n ) with a loss function 𝓁 is defined as:

R(𝐖n ) = 𝐄z,z� [𝓁(𝐖n , z, z� )]where 𝐄z,z� [𝓁(⋅, ⋅, ⋅)] denotes the expectation of function 𝓁 when z and z′ are sampledaccording the distribution P .
Empirical Error. The empirical error RDn (𝐖n ) is defined on the training data set Dn as :RDn (𝐖n ) =1 ∑ ∑𝓁(𝐖n , zi , zj )n2 z ∈D z ∈D
i

n j

n

Expected Generalization Gap. As Algorithm 1 is based on a randomized procedure, we use
the definition of the expected generalization gap as follows:

𝕂n = 𝐄SGD [R(𝐖n ) − RDn (𝐖n )]where 𝐄SGD denotes the expectation taken over the inherent randomness of the stochastic
algorithm.

4.2 Uniform stability of the stochastic learning algorithm
Intuitively, a learning algorithm is said to have a uniform stability if its output is stable
under a small modification of the training data set. For a randomized learning algorithm,
the uniform stability property is defined as follows:

Definition 2 (Uniform Stability of the randomized algorithm) A randomized algorithm 𝔸
is 𝛽n-uniformly stable with respect to a loss function 𝓁 , if the following inequality holds:∀(Dn , k), sup |𝐄SGD [𝓁(𝐖n , z, z� )] − 𝐄SGD [𝓁(𝐖n,k , z, z� )]| ≤ 𝛽nz,z�where Dn,k is the new data set obtained from Dn by replacing zk ∈ Dn with a new exampleẑ k sampled from P ; 𝐖n and 𝐖n,k are the outputs of 𝔸 trained on two data sets Dn and Dn,k ,respectively.
In order to prove that Algorithm 1 has the uniform stability property, we need the following lemma (its proof is placed in the "appendix" section):
Lemma 4 Let the loss function 𝓁 defined in the problem (13) be 𝛽-smooth and L-Lipschitz;
and 𝐖(T)be the parameters of the parametric form (12) trained on Dn and Dn,k ,
𝐖(T)
n
n,k
respectively, using Algorithm 1 with the number of iterations T and learning rate 𝜇. Then,the expected difference in the model parameters is upper bounded by:

13

1597

Machine Learning (2021) 110:1585–1607

]
[
4
‖
‖
− 𝐖(T)
𝐄SGD ‖𝐖(T)‖ ≤ 𝜇TL
n
n,k
‖
‖2
n

(18)

Using Lemma 4 and L-Lipschitz property of function 𝓁 (see Lemma 3), we can now
prove the stability of Algorithm 1.
Theorem 5 [Uniform Stability of Algorithm 1] Let the loss function 𝓁 defined in the problem (13) be 𝛽-smooth and L-Lipschitz. Then Algorithm 1 with the fixed learning rate 𝜇 is kn-uniformly stable where kn = 4n 𝜇TL2.
Proof We have the following inequalities:

4
2‖|𝐄SGD [𝓁(𝐖n , z, z� )] − 𝐄SGD [𝓁(𝐖n,k , z, z� )]| ≤ L𝐄SGD ‖‖𝐖n − 𝐖n,k ‖2 ≤ n 𝜇TL(19)where the first and second inequalities are obtained by the L-Lipschitz property of 𝓁 andLemma 4, respectively. This completes the proof. □

4.3 Bound on generalization gap
Using the property of uniform stability in the previous subsection, we can derive the generalization bound which is done by the McDiarmid inequality (McDiarmid 1989).
Theorem 6 [McDiarmid inequality (McDiarmid 1989)] Let X1 , X2 , ..., Xn be n independent
random variables taking values in X and let Z = f (X1 , ..., Xn ). If, for each 1 ≤ i ≤ n, there
exists a constant ci such thatsupx1 ,...,xn ,x� |f (x1 , ..., xn ) − f (x1 , ..., xi� , ..., xn )| ≤ ci , ∀1 ≤ i ≤ n,i��2then for any 𝜖 > 0, Pr[�Z − 𝐄[Z]� ≥ 𝜖] ≤ 2exp ∑−2𝜖
n
2
c
i=1 i

[ ]
To derive the bound on R(Wn ), we replace Z by Kn in Theorem 6 and bound 𝐄SGD 𝕂n
and |𝕂n − 𝕂n,k | which are established by the following lemmas (see their proofs in the"appendix’ section).
Lemma 7 For the loss function satisfying a uniform stability in kn, we have the following
inequality:(20)

𝐄Dn [𝕂n ] ≤ 2kn

Lemma 8 For the loss function satisfying a uniform stability in kn and upper bounded byM, we have the following inequality:|𝕂n − 𝕂n,k | ≤ 2kn +

2M
n

(21)

Now we can derive the generalization bound for R(𝐖n ) in the following theorem:

13

Machine Learning (2021) 110:1585–1607

1598

Theorem 9 Let Dn be a training data set with n samples, 𝐖n be the solution obtained byminimizing the optimization problem (13) using Algorithm 1 with uniform stability kn. Thenthe following inequality holds for probability of at least 1 − 𝛿 (0 ≤ 𝛿 ≤ 1):
√
[
]
2
2
(22)
𝐄SGD R(𝐖n ) − RDn (𝐖n ) ≤ 2kn + (nkn + M)
logn𝛿
Proof Applying McDiarmid’s concentration inequality (6) by replacing Z with 𝕂n, we
have:

⎛
⎞
� �
⎜
⎟
−2𝜖 2
Pr(𝕂n − 𝐄Dn 𝕂n ≥ 𝜖) ≤ 2exp⎜ �
�2 ⎟⎜ n 2k + 2M ⎟
n
⎝
⎠
n

(
By fixing 𝛿 = 2exp
proof of Theorem 9. 

)

−2𝜖 2

)2
(
n 2kn + 2Mn

, we get 𝜖 = (nkn + M)√

2
log 𝛿1nwhich completes the□

The generalization bound is meaningful if the bound converge to 0 as n → ∞. Ourderived generalization bound converges to 0 as kn decays with O( 1n ). This confirms Algorithm 1 converges.

5 Experiments
In this section, we demonstrate the benefit of learning subtree pattern weights by experiments on both synthetic and real-world data. We performed classification experiments using
the C-SVM implementation LIBSVM (Chang and Lin 2011). The necessary parameters of
SVM were selected by cross-validation on the training set. These are the regularization
parameter C ∈ {10−3 , 10−2 , ..., 102 , 103 } and kernel parameter 𝛾 ∈ {0.0001, 0.001, 0.01}.
For learning the weights 𝐖 of subtree patterns (or WL labels) in Algorithm 1, the learning rate 𝜇 and maximum number of iterations T were set as 0.0001 and 500, respectively;𝜖h ∈ {0.1, 0.5, 1.0} were selected by cross validation based on the training set and 𝐜h was
fixed as a vector of ones, i.e. 𝐜h = 𝟏|Σh | for h = 1, .., H ; the hyperparameters 𝛼1, 𝛼2 and 𝜎
were empirically determined as 1.0, 0.5 and 0.1, respectively. All kernels were implemented in Python 3.0 and experiments were conducted on an Intel Core i9 at 2.3 Ghz with
64GB of RAM using a single processor only. The source code can be accessed throughhttps://​github.​com/​haidn​guyen​0909/​weigh​tedWWL.

5.1 Synthetic data
We designed eight substructures, shown in Fig. 3, in which substructures indexed by 1, 2,
5 and 6 are assumed to be indicative to positive class (+ 1) as they have a pattern ’1–0(−2)–0’ in common. The others are indicative to negative class (− 1). Our synthetic data set
consists of eight groups of graphs, each corresponds to one of these eight substructuresby randomly adding noisy nodes and edges. We used groups corresponding substructures
1, 2, 3 and 4 as training data and the others as testing data. We constructed two kernels
for graphs: WWL and the proposed method with number of WL iterations H = 2, then13

Machine Learning (2021) 110:1585–1607

1599

Fig. 3  Designed substructures 1–8: substructures 1, 2, 5 and 6 are indicative to positive class as they contain pattern ’1–0(− 2)–0’ (emphasized in yellow). The rest are indicative to negative class. Graphs are generated from substructures by adding random nodes and noisy edges (20 examples for each substructure).
Graphs generated from the substructures 1, 2, 3 and 4 are used for training. Graph generated from the substructures 5, 6, 7 and 8 are used for testingFig. 4  Examples of selected subtree patterns of the first level h = 1 (left) and their weights learned by the
proposed learning algorithm (right)

used SVM for classification. We reported mean accuracy obtained by ten synthetic data
sets generated in this way.
We observed that WWL obtained mean accuracies of 82.4%, while the proposed
method achieved significantly higher accuracy of 95%. It is noted that the testing examples
were confusing the classifier. For instance, the substructure 5 has the same similarity with
substructures 2 (indicative to positive class) and 4 (indicative to negative class) according to the WWL kernel, making it hard for the classifier to distinguish graphs containing
the substructure 5. In contrast, this confusion can be alleviated by learning subtree patternweights. In particular, the pattern ’1–0(− 2)–0’ present in the substructure 5 is assigned
a high weight by the proposed method (see Fig. 4). Therefore, graphs generated from thesubstructure 5 are more likely to be classified as positive.

13

Machine Learning (2021) 110:1585–1607

1600
Table 1  Statistics of datasets used in experiments
Datasets#Graphs#Classes

Avg. card (V)Avg. card (E)#labels

MUTAG​

188

2 (125 vs. 63)

17.9

39.6

7

PTC-MRPROTEIN
NCI1

344
1113
4110

2 (192 vs. 152)
2 (663 vs. 450)
2 (2053 vs. 2057 )

25.6
39.1
N/A

51.9
145.63
N/A

19
3
N/A

5.2 Real‑world data
In this subsection we present an experimental evaluation of the proposed method on realworld data. We report experimental results on four benchmark bioinformatics data sets,
involving node-labeled graphs, particularly, MUTAG, PTC-MR, PROTEIN and NCI1.
The MUTAG dataset consists of graph structures of 188 chemical compounds which areeither mutagenic aromatic or heteroromatic nitro compounds and nodes can have 7 discrete labels. The PTC-MR dataset consists of 344 chemical compounds which are known
to cause or not cause cancer in rats and mice, and nodes can have 19 discrete labels. The
PROTEIN dataset consists of relations between secondary structure elements represented
by nodes and neighborhood in the amino-acid sequence or in 3D space by edges, and nodes
can have 3 discrete labels. The NCI1 dataset is a balanced subset of chemical compounds
screened for their ability to inhibit the growth of a panel of human tumor cell lines, andnodes can have 37 discrete labels. Some statistics of these data sets are shown in Table 1.
We compared the proposed method to several state-of-the-art graph kernels. Due to thelarge number of graph kernels in the literature, we selected representatives of the majorfamilies of graph kernels. In particular, for the family of walk based kernels, we compared
the proposed method to the fast random walk kernel (Kashima 2003) that essentially counts
the common labeled walks. For the family of path based graph kernels, we compared to theshortest path kernel (Borgwardt and Kriegel 2005). For the family of WL based graph kernels, we compared to WL subtree (Vishwanathan et al. 2010), WL-OA (Kriege et al. 2016)
and WWL kernels (Togninalli et al. 2019). The WL based kernels have been shown to besuperior to previous approaches.
We report mean predictive accuracies and standard deviations obtained by 10-foldcross-validation repeated 10 times with random fold assignments. Within each fold, thenumber of hops H ∈ {1, 2, ..., 6} was selected by cross validation based on the training
set. The results evaluated by classification accuracy are summarised in Table 2. We used
one-sided paired t-test to verify if the accuracy differences between two methods on data
sets are statistically significant. We empirically observed that random walk and shortest
path kernels were less competitive to WL-based kernels on four data sets. On three datasets MUTAG, PROTEIN and NCI1, the proposed method was comparable with WL-OA
while it improved its unweighted version WWL by 1.4%, 1.5% and 0.8%, respectively(the calculated p values were 0.061, 0.0087 and 0.055, respectively, smaller than the significance level of 𝛼 = 0.1). On PTC-MR, the proposed method improved WWL by 0.6%
while outperforming WL-OA by nearly 5% (the calculated p values were 0.0035 and <0.001, respectively). In all these data sets, random walk, shortest path and WL subtree kernels were dominated by the rest in large margins. Furthermore, we also investigated someselected subtree patterns at the first and second levels (h = 1, 2) along with their weights
learned by the proposed algorithm from two data sets: MUTAG and PTC-MR (see Figs. 5

13

Machine Learning (2021) 110:1585–1607

1601

Fig. 5  Examples of selected subtree patterns and their weights of the first level h = 1 a and the second levelh = 2 b extracted from MUTAG by Algorithm 1

and 6, respectively). Interestingly, the weights were found different over substructures,indicating their different degrees of importance in the prediction task. These experimental results confirmed the effectiveness of learning important subtree patterns for WWL
kernels.

5.3 Computational efficiency of the proposed stochastic algorithm
In this subsection, we evaluate the computational efficiency of the proposed Algorithm 1.
We empirically compared two variants: batch and stochastic (Algorithm 1), for solving theminimization problem (13) in terms of running time. The first variant considers all pairs ofgraphs for every step of projected gradient descent. The second variant considers one pairof graphs at a time to take a single step.
First we assessed the running time of two variants on randomly generated graphs (asdescribed in Subsect. 5.1) with respect to two parameters: number of graphs N and number of WL iterations H. We varied N in range {50, 100, 200, 400, 600, 800, 1000} and H
in range {1, 2, 3, 4, 5, 6, 7}. For each individual experiment, we fixed one parameter at its
default value and varied the other. The default values were 100 for N and 2 for H. We
report CPU running times in seconds in Fig. 7. Empirically, we observed that the running
time of full batch variant increased quickly when increasing the number of graphs N andthe number of WL iterations H. In contrast, the stochastic variant scaled much better withmuch lower running times, indicating that Algorithm 1 has high scalability for large scaledata sets. The computational efficiency of Algorithm 1 can be explained by the fact that13

Machine Learning (2021) 110:1585–1607

1602

Fig. 6  Examples of selected subtree patterns and their weights of the first level h = 1 a and the second levelh = 2 b extracted from PTC-MR by Algorithm 1
Table 2  Classification accuracies and standard deviation on real-world graph data sets: MUTAG, PTC-MR,
PROTEIN and NCI1
Kernels

MUTAG​

PTC-MRPROTEIN

NCI1

Random Walk (Kashima 2003)

85.06 ± 0.18 55.74 ± 3.64 71.11 ± 0.83 62.88 ± 0.22

Shortest path (Borgwardt and Kriegel 2005)
WL Subtree (Shervashidze and Borgwardt
2009)
WL-OA (Kriege et al. 2016)
WWL (Togninalli et al. 2019)
Proposed

85.49 ± 0.59 53.29 ± 0.92 73.03 ± 1.13 61.36 ± 0.19
85.61 ± 0.85 61.89 ± 1.97 72.5 ± 0.32 85.61 ± 0.13

Best results are shown in bold

13

88.17 ± 1.98 60.49 ± 1.39 75.89 ± 0.41 86.17 ± 0.35
86.95 ± 1.35 64.86 ± 1.57 74.25 ± 0.74 85.69 ± 0.28
88.37 ± 1.82 65.44 ± 0.97 75.73 ± 0.57 86.45 ± 0.11

1603

Machine Learning (2021) 110:1585–1607

Fig. 7  Running time in seconds on synthetic data sets of two variants: full batch and stochastic algorithms,
for learning the subtree pattern importance W (in the optimization problem (13)) (Default values: dataset
size N = 100, WL iteration H = 2)
Table 3  Running time in seconds of two variants: full batch and stochastic on real-world data sets:
MUTAG, PTC-MR, PROTEIN and NCI1
Variants/data sets

MUTAG​

PTC-MRPROTEIN

NCI1

Full batch

1’ 32”

3’ 48”

1 h 20’40”

5 h 5’

Stochastic

4h 29’

8h 50’> 3 days> 1 week

computing gradient of the loss function for a graph pair g and g′ involves a few substruc′tures) WL labels) v shared by g and g in Eq. (12), i.e., sparsity of the feature vector( (or𝐙 g, g′ .
Second we assessed the running time of two variants on real-world data sets: MUTAG,
PTC-MR, PROTEIN and NCI1. We reported the running time of two variants to finish theentire classification tasks, including learning subtree pattern weights, computing kernels
and doing classification, in Table 3. The running time were taken average by 10-fold crossvalidation. We empirically observed that the full batch variant was slow when running oneven small data sets MUTAG and PTC-MR, taking in approximately 4 h and 9 h, respectively. But the stochastic variant was much faster, taking only less than 4 min on the twodata sets. We also observed that the stochastic variant could easily scale up to data setswith thousands of graphs. Particularly, on data sets PROTEIN and NCI1, the tasks were
performed in nearly 1h 30’ and 5 h, respectively. However, it was impossible for the fullbatch variant to finish the tasks in less than 3 days. These evidence showed that the proposed stochastic variant is highly scalable.

6 Conclusion and discussion
In this work, we proposed to learn the weights of substructures of graphs, particularly,
subtree patterns (extracted by WL labeling scheme), to overcome the limitations of current graph kernels. We considered the problem of incorporating subtree pattern weights for
WWL kernel (Togninalli et al. 2019) by formulating the parametric form of Wasserstein

13

1604

Machine Learning (2021) 110:1585–1607

distance taking into account subtree pattern weights, and learning these weights from data
optimally for the tasks.
Our proposed method has several advantages. First, it can learn the importance of subtree patterns specifically for the tasks through their weights in the parametric distance function. Second, the kernels converted from the learned parametric function of subtree patternweights are valid. Third, the efficient stochastic algorithm for learning the weights has highscalability for large scale data sets, and its theoretical guarantees are provided.
Although we considered WWL kernel for extracting important subtree patterns, aninteresting and worthwhile extension of our work would be to apply this idea to other WL
based graph kernels such as WL subtree and WL-OA kernels. The improvements of the
optimization algorithm for learning subtree pattern weights in terms of both convergenceand efficiency would also be our future work.

A AppendicesA.1 Proof of Lemma 4.
Proof We prove the lemma by following the notion of using the same randomness for two
data set Dn and Dn,k as in Hardt et al. (2016). Particularly, we supply the sample sequences
S = {p1 = (zi1 , zj1 ), ..., pT = (ziT , zjT )} to two identical learning algorithms except that forsome t (1 ≤ t ≤ T ), if pt contains zk ( pt = (zk , zjt ) or (zit , zk )), we replace it with p̂ t = (̂zk , zjt )or (zit , ẑ k ). So, there are two cases to consider:
Case 1: At step t, Algorithm 1 picks a pair of samples (z, z′) that contain no zk ( z ≠ zk
and z′ ≠ zk ) and this case occurs with probability (1 − 1n )2. Then, we have:
)]‖2
(
[
‖
[ (t)
( (t)
)]
‖ (t+1)(t+1) ‖2
(t)
(t)
�
�
‖
,
z,
z
proj
−
𝜇∇𝓁
𝐖
𝐖
−
𝜇∇𝓁
𝐖
,
z,
z
−
proj
𝐖
‖𝐖n − 𝐖n,k ‖ = ‖
C
C
‖
n
n
n,k
n,k
‖
‖2 ‖
‖
‖2
(
)‖2
)
‖ (t)
( (t)
(t)
�
t
�
‖
≤‖
‖𝐖n − 𝜇∇𝓁 𝐖n , z, z − 𝐖n,k + 𝜇∇𝓁 𝐖n,k , z, z ‖
‖
‖2
(
(
)( (
))
)
‖ (t)
(t) ‖2
(t)
(t)
(t)
�
= ‖𝐖n − 𝐖n,k ‖ − 2𝜇 𝐖n − 𝐖n,k ∇𝓁 𝐖n , z, z − ∇𝓁 𝐖(t), z, z�
n,k
‖
‖2
(
)‖2
‖
(t)
(t)
�
� ‖
+ 𝜇2 ‖‖∇𝓁(𝐖n , z, z ) − ∇𝓁 𝐖n,k , z, z ‖
‖2
‖
(
)
(
)‖2
‖ ( (t))2‖ (t)(t) ‖2−𝜇 ‖≤ ‖𝐖n − 𝐖n,k ‖ − 𝜇∇𝓁 𝐖n , z, z� − ∇𝓁 𝐖(t), z, z� ‖
‖
‖
n,k
‖
‖2
𝛽
‖2
‖

�
The second line is obtained by the fact that ��projC [𝐮] − projC [𝐯]�2 ≤ ‖𝐮 − 𝐯‖2 for 𝐮, 𝐯 in
the domain. The last line is obtained by the 𝛽-smoothness of function 𝓁 (see Lemma 3). So
we have the following inequality (by selecting 𝜇 ≤ 𝛽2):‖ (t+1)‖‖(t+1) ‖− 𝐖(t)‖𝐖n − 𝐖n,k ‖ ≤ ‖𝐖(t)n,k ‖‖‖2 ‖ n‖2(23)Case 2: At step t, Algorithm 1 picks a pair of samples (z, z′) that contain zk ( z = zk andz� = zk ) and this case occurs with probability 1 − (1 − 1n )2. Then we have:

13

1605

Machine Learning (2021) 110:1585–1607‖ (t+1)‖‖(t+1) ‖+ 2𝜇L
− 𝐖(t)‖𝐖n − 𝐖n,k ‖ ≤ ‖𝐖(t)n,k ‖‖‖2 ‖ n‖2(24)

The above inequality holds as the norm of gradient of loss function 𝓁 is upper bounded by
L. From two inequalities (23) and (24), and considering the probabilities of two cases, wehave the following inequality:

)
]
[
] (
[
1 2
‖
‖
(t) ‖
(t+1) ‖−𝐖
𝐄SGD ‖𝐖(t)𝐄SGD ‖𝐖n(t+1) − 𝐖n,k‖
‖ ≤ 1−n,k ‖2
‖ n
‖
‖2
n

(25)

(
(
) )(
[
]
)
1 2
‖
‖
+ 1− 1−− 𝐖(t)
𝐄SGD ‖𝐖(t)‖ + 2𝜇L
n
n,k
‖
‖2
n

(26)

(
[
]
(
) )
1 2
‖
(t) ‖
−
𝐖
= 𝐄SGD ‖𝐖(t)
+
2𝜇L
1
−
1
−
n,k ‖
‖ n
‖2
n

(27)

By the above recursive formula, we obtain the following:
(
[
]
(
) )
1 2
4
‖ (T)
(T) ‖
𝐄SGD ‖𝐖n − 𝐖n,k ‖ ≤ 2 1 − 1 −𝜇TL ≤ 𝜇TL
‖
‖2
n
n

(28)

which completes the proof. □

A.2 Proof of Lemma 7
Proof By the definition of 𝕂n as in (4.1), we have the following:
𝐄Dn [𝕂n ] = 𝐄Dn 𝐄SGD [R(𝐖n ) − RDn (𝐖n )]1 ∑= 𝐄Dn 𝐄SGD 𝐄z,z� 𝓁(𝐖n , z, z� ) − 𝐄Dn 𝐄SGD 2𝓁(𝐖n , zi , zj )n i,j=1
[
]
n
n
]
1 ∑∑[
�
�
�
𝓁(𝐖n , z, z ) − 𝓁(𝐖n , zk , z ) + 𝓁(𝐖n , zk , z ) − 𝓁(𝐖n , zk , zj )= 𝐄Dn ,z,z� 𝐄SGD 2
n k=1 j=1
[
]
n
n
]
1 ∑∑[
�
�
�
𝓁(𝐖n , z, z ) − 𝓁(𝐖n , zk , z ) + 𝓁(𝐖n , zk , z ) − 𝓁(𝐖n , zk , zj )= 𝐄Dn ,z,z� 𝐄SGD 2
n k=1 j=1
[
]
n
n
]
1 ∑∑[
�
�
𝓁(𝐖n , z, z ) − 𝓁(𝐖n , zk , z )= 𝐄Dn ,z,z� 𝐄SGD 2
n k=1 j=1⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟n[+ 𝐄Dn ,z,z� 𝐄SGD

(a)

]

n
n
1 ∑∑[

�

]

𝓁(𝐖n , zk , z ) − 𝓁(𝐖n , zk , zj )n2 k=1 j=1⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟(b)

We first process the part (a) which is equivalent to the following:13

Machine Learning (2021) 110:1585–1607

1606

(a) =

n
n
[
] 1∑
[
]
1∑
𝐄Dn ,z,z� 𝐄SGD 𝓁(𝐖n , z, z� ) −𝐄Dn ,z,z� 𝐄SGD 𝓁(𝐖n , zk , z� )
n k=1
n k=1

=

n
n
[
[
]
] 1∑
1∑
𝐄Dn ,z,z� 𝐄SGD 𝓁(𝐖n , z, z� ) −𝐄Dn ,̂zk ,z� 𝐄SGD 𝓁(𝐖n , zk , z� )
n k=1
n k=1

=

n
n
[
[
]
] 1∑
1∑
𝓁(𝐖n,k , z, z� )𝐄Dn ,z,z� 𝐄SGD 𝓁(𝐖n , z, z� ) −
𝐄
�𝐄
n k=1
n k=1 Dn ,z,z SGD

≤

n
[
]
1∑
𝐄
𝓁(𝐖n , z, z� ) − 𝓁(𝐖n,k , z, z� ) ≤ kn
�𝐄
n k=1 Dn ,z,z SGD

Similarly, we also prove that ( b) ≤ kn which completes the proof. □

A.3 Proof of Lemma 8
Proof By the definition of 𝕂n as in (4.1), we have:|𝕂n − 𝕂n,k | = |𝐄SGD [R(𝐖n ) − RDn (𝐖n )] − 𝐄SGD [R(𝐖n,k ) − RDn,k (𝐖n,k )]|≤ |𝐄SGD R(𝐖n ) − 𝐄SGD R(𝐖n,k )| + |𝐄SGD RDn (𝐖n ) − 𝐄SGD RDn,k (𝐖n,k )|⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟
(c)

(d)

We bound the two terms (c) and (d) as follows:(c) = |𝐄SGD 𝐄z,z� [𝓁(𝐖n , z, z� ) − 𝓁(𝐖n,k , z, z� )]|≤ 𝐄z,z� |𝐄SGD [𝓁(𝐖n , z, z� ) − 𝓁(𝐖n,k , z, z� )]| ≤ kn
1 ∑ ∑
1 ∑ ∑
(d) = |𝐄SGD 2𝓁(𝐖n , zi , zj ) − 𝐄SGD 2𝓁(𝐖n,k , zi , zj )|n z ∈D z ∈D
n z ∈D z ∈D
i

n j

n

i

n,k

j

n,k

1 ∑
𝐄 [𝓁(𝐖n , zi , zj ) − 𝓁(𝐖n,k , zi , zj )]=| 2
n i≠k,j≠k SGD
1 ∑+ 2
𝐄 [𝓁(𝐖n , zi , zk ) − 𝓁(𝐖n,k , zi , zk )]n i≠k SGD1 ∑+ 2
𝐄 [𝓁(𝐖n , zk , zj ) − 𝓁(𝐖n,k , zk , zj )]n j≠k SGD≤

The two above inequalities complete the proof. 2(n − 1)(n − 1)2
2M
kn +
Mkn +
n
n2
n2
□

Acknowledgements D. H. N. has been supported in part by Otsuka Toshimi scholarship and JSPS Research
Fellowship for Young Scientists (DC2) with KAKENHI [grant number 19J14714]. C. H. N. has been supported in part by MEXT KAKENHI [grant number 18K11434]. H. M. has been supported in part by JST
ACCEL [grant number JPMJAC1503], MEXT KAKENHI [grant numbers 16H02868, 19H04169], FiDiPro
by Tekes (currently Business Finland) and AIPSE program by Academy of Finland.

13

Machine Learning (2021) 110:1585–1607

1607

References
Borgwardt, K.M., & Kriegel, H.P. (2005). Shortest-path kernels on graphs. In: Fifth IEEE International
Conference on Data Mining (ICDM’05), pp. 8–pp. IEEE
Chang, C. C., & Lin, C. J. (2011). Libsvm: A library for support vector machines. ACM Transactions onIntelligent Systems and Technology (TIST), 2(3), 1–27.
Hardt, M., Recht, B., & Singer, Y. (2016). Train faster, generalize better: Stability of stochastic gradientdescent. In: International Conference on Machine Learning, pp. 1225–1234. PMLR
Haussler, D. (1999). Convolution kernels on discrete structures. Tech. rep., Technical report, Department ofComputer Science, University of California.
Kashima, H., Tsuda, K., & Inokuchi, A. (2003). Marginalized kernels between labeled graphs. In: Proceedings of the 20th International Conference on Machine Learning (ICML-03), pp. 321–328
Kriege, N.M., Giscard, P.L., Wilson, R. (2016). On valid optimal assignment kernels and applications tograph classification. In: Advances in Neural Information Processing Systems, pp. 1623–1631
Kulis, B., et al. (2012). Metric learning: A survey. Foundations and Trends in Machine Learning, 5(4),287–364.
McDiarmid, C. (1989). On the method of bounded differences. Surveys in Combinatorics, 141(1), 148–188.
Munkres, J. (1957). Algorithms for the assignment and transportation problems. Journal of the Society for
Industrial and Applied Mathematics, 5(1), 32–38.
Nesterov, Y. (2005). Smooth minimization of non-smooth functions. Mathematical Programming, 103(1),127–152.
Scott, J. (2011). Social network analysis: Developments, advances, and prospects. Social Network Analysis
and Mining, 1(1), 21–26.
Sharan, R., & Ideker, T. (2006). Modeling cellular machinery through biological network comparison.
Nature Biotechnology, 24(4), 427–433.
Shervashidze, N., & Borgwardt, K. (2009). Fast subtree kernels on graphs. In: Advances in neural information processing systems, pp. 1660–1668
Togninalli, M., Ghisu, E., Llinares-López, F., & Rieck, B. (2019). Borgwardt, K.: Wasserstein weisfeilerlehman graph kernels. In: Advances in Neural Information Processing Systems, pp. 6439–6449
Trinajstic, N. (2018). Chemical graph theory. London: Routledge.
Villani, C. (2008). Optimal transport: Old and new (Vol. 338). NewYork: Springer Science & BusinessMedia.
Vishwanathan, S. V. N., Schraudolph, N. N., Kondor, R., & Borgwardt, K. M. (2010). Graph kernels. The
Journal of Machine Learning Research, 11, 1201–1242.
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

13