# INCREASE AND CONQUER: TRAINING GRAPH NEURAL NETWORKS ON GROWING GRAPHS

**Anonymous authors**
Paper under double-blind review

ABSTRACT

Graph neural networks (GNNs) use graph convolutions to exploit network invariances and learn meaningful features from network data. However, on large-scale
graphs convolutions incur in high computational cost, leading to scalability limitations. Leveraging the graphon — the limit object of a graph — in this paper
we consider the problem of learning a graphon neural network (WNN) — the
limit object of a GNN — by training GNNs on graphs sampled Bernoulli from the
graphon. Under smoothness conditions, we show that: (i) the expected distance
between the learning steps on the GNN and on the WNN decreases asymptotically with the size of the graph, and (ii) when training on a sequence of growing
graphs, gradient descent follows the learning direction of the WNN. Inspired by
these results, we propose a novel algorithm to learn GNNs on large-scale graphs
that, starting from a moderate number of nodes, successively increases the size
of the graph during training. This algorithm is benchmarked on both a recommendation system and a decentralized control problem where it is shown to retain
comparable performance to its large-scale counterpart at a reduced computational
cost.

1 INTRODUCTION

Graph Neural Networks (GNNs) are deep convolutional architectures formed by a succession of
layers where each layer composes a graph convolution and a pointwise nonlinearity (Wu et al., 2021;
Zhou et al., 2020). Tailored to network data, GNNs have been used in a variety of applications such
as recommendation systems (Fan et al., 2019; Tan et al., 2020; Ying et al., 2018; Schlichtkrull et al.,
2018; Ruiz et al., 2019a) and Markov chains (Qu et al., 2019; Ruiz et al., 2019b; Li et al., 2015),
and fields such as biology (Fout et al., 2017; Duvenaud et al., 2015; Gilmer et al., 2017; Chen et al.,
2020) and robotics (Qi et al., 2018; Gama & Sojoudi, 2021; Li et al., 2019). Their success in these
fields and applications provides ample empirical evidence of the ability of GNNs to generalize to
unseen data. More recently, their successful performance has also been justified by theoretical works
showing that GNNs are invariant to relabelings (Chen et al., 2019; Keriven & Peyr´e, 2019), stable
to graph perturbations (Gama et al., 2020) and transferable across graphs (Ruiz et al., 2020a).

One of the most important features of a GNN is that, because the linear operation is a graph convolution, its number of parameters does not depend on the number of nodes of the graph. In theory,
this means that GNNs can be trained on graphs of any size. In practice, however, if the graph has
large number of nodes training the GNN is costly because computing graph convolutions involves
large matrix operations. While this issue could be mitigated by transferability — training the GNN
on a smaller graph to execute on the large graph —, this approach does not give any guarantees
on the distance between the optimal solutions on the small and on the large graph. In other words,
when executing the GNN on the target graph we do not know if its error will be dominated by the
transferability error or by the generalization error from training.

In this paper, we address the computational burden of training a GNN on a large graph by progressively increasing the size of the network. We consider the limit problem of learning an “optimal”
neural network for a graphon, which is both a graph limit and a random graph model (Lov´asz,
2012). We postulate that, because sequences of graphs sampled from the graphon converge to it, the
so-called graphon neural network (Ruiz et al., 2020a) can be learned by sampling graphs of growing
size and training a GNN on these graphs (Algorithm 1). We prove that this is true in two steps. In


-----

Theorem 1, we bound the expected distance between the gradient descent steps on the GNN and on
the graphon neural network by a term that decreases asymptotically with the size of the graph. A
consequence of this bias bound is that it allows us to quantify the trade-off between a more accurate
gradient and one that could be obtained with less computational power. We then use this theorem to
prove our main result in Theorem 2, which is stated in simplified form below.

**Theorem (Graphon neural network learning, informal) Let W be a graphon and let {Gn} be a**
sequence of growing graphs sampled from W. Consider the graphon neural network Φ(W) and
assume that it is learned by training the GNN Φ(Gn) with loss function ℓ(yn, Φ(Gn)) on the
sequence **Gn** . Over a finite number of training steps, we obtain
_{_ _}_

_∥∇ℓ(Y, Φ(W)∥≤_ _ϵ with probability 1._

The most important implication of this result is that the learning iterates computed in the sequence of
growing graphs follow the direction of the graphon gradient up to a small ball, which provides theoretical validation to our cost efficient training methodology. We also validate our algorithm in two
numerical experiments. In the first, we learn a GNN-based recommendation system on increasingly
large subnetworks of a movie similarity graph, and compare it with the recommendation system
trained on the full graph. In the second, we consider the problem of flocking and train GNNs to
learn the actions agents need to take to flock. We compare the results obtained when progressively
increasing the number of agents during training and when training directly on the target graph.

2 RELATED WORK

GNNs are data processing architectures that follow from the seminal works in the areas of deep
learning applied to graph theory (Bruna et al., 2013; Defferrard et al., 2016; Gori et al., 2005; Lu &
Getoor, 2003). They have been successfully used in a wide variety of statistical learning problems
(Kipf & Welling, 2016; Scarselli et al., 2018), where their good performance is generally attributed
to the fact that they exploit invariances present in network data (Maron et al., 2019; Gama et al.,
2018; Chami et al., 2021).

More recently, a number of works show that GNNs can be transferred across graphs of different
sizes (Ruiz et al., 2020a; Levie et al., 2019; Keriven et al., 2020). Specifically, (Ruiz et al., 2020a)
leverages graphons to define families of graphs within which GNNs can be transferred with an
error bound that decreases asymptotically with the size of the graph. The papers by Levie et al.
(2019) and Keriven et al. (2020) offer similar results by considering the graph limit to be a generic
topological space and a random graph model respectively. In this paper, we use an extension of the
transferability bound derived in Ruiz et al. (2020a) to propose a novel learning algorithm for GNNs.

3 PRELIMINARY DEFINITIONS

3.1 GRAPH NEURAL NETWORKS

Graph neural networks exploit graph symmetries to extract meaningful information from network
data (Ruiz et al., 2020c; Gama et al., 2020). Graphs are represented as triplets Gn = (V, E, W ),
where V, |V| = n, is the set of nodes, E ⊆V × V is the set of edges and W : E → R is a map
assigning weights to each edge. The graph Gn can also be represented by the graph shift operator
(GSO) S ∈ R[n][×][n], a square matrix that respects the sparsity of the graph. Examples of GSOs include
the adjacency matrix A, the graph Laplacian L = diag(A1) − **A and their normalized counterparts**
Gama et al. (2018). In this paper we consider the graph Gn to be undirected and fix S = A/n.

Graph data is represented in the form of graph signals. A graph signal x = [x1, . . ., xn][T] _∈_ R[n]
is a vector whose i-th component corresponds to the information present at the i-th node of graph
**Gn. A basic data aggregation operation can be defined by applying the GSO S to graph signals x.**
The resulting signal z = Sx is such that the data at node i is a weighted average of the information
in the 1-hop neighborhood of i, zi = _j_ _i_ [[][S][]][ij][x][j][ where][ N][i][ =][ {][j][ |][ [][S][]][ij][ ̸][= 0][}][. Information]

_∈N_
coming from further neighborhoods can be aggregated by successive applications of the GSO, also
called shifts. Using this notion of shift, graph convolutions are defined by weighting the contribu
[P]

tion of each successive application of S to define a polynomial in the GSO. Explicitly, the graph


-----

convolutional filter with coefficients h = [h0, . . ., hK−1] is given by


_K−1_

_hkS[k]x_ (1)

_k=0_

X


**y = h ∗S x =**


where ∗S denotes the convolution operation with GSO S.

Since the adjacency matrix of an undirected graph is always symmetric, the GSO admits an eigendecomposition S = VΛV[H] . The columns of V are the graph eigenvectors and the diagonal elements of Λ are the graph eigenvalues, which take values between −1 and 1 and are ordered as
_−orthonormal basis of1 ≤_ _λ−1 ≤_ _λ−2 ≤ R[n]. . ., we can project (1) onto this basis to obtain the spectral representation of ≤_ 0 ≤ _. . . ≤_ _λ2 ≤_ _λ1 ≤_ 1. Since the eigenvectors of S form an
the graph convolution, which is given by


_K−1_

_hkλ[k]._ (2)

_k=0_

X


_h(λ) =_


Note that (2) only depends on the hk and on the eigenvalues of the GSO. Hence, as a consequence
of the Cayley-Hamilton theorem, convolutional filters may be used to represent any graph filter with
spectral representation h(λ) = f (λ) where f is analytic (Strang, 1976).

Graph neural networks are layered architectures where each layer consists of a graph convolution
followed by a pointwise nonlinearity ρ, and where each layer’s output is the input to the following
layer. At layer l, a GNN can output multiple features x[f]l [,][ 1][ ≤] _[f][ ≤]_ _[F][l][ which we stack in a]_
matrix Xl = [x[1]l _[, . . .,][ x]l[F][l]_ []][ ∈] [R][n][×][F][l] [. Each column of the feature matrix is the value of the graph]
signal at feature f . To map the Fl 1 features coming from layer l 1 into Fl features, Fl 1 _Fl_
convolutions need to be computed, one per input-output feature pair. Stacking their weights in− _−_ _−_ _× K_
matrices Hlk ∈ R[F][l][−][1][×][F][l], we write the l-th layer of the GNN as

_K−1_

**Xl = ρ** **S[k]Xl** 1Hlk _._ (3)

_−_
_k=0_ !

X

In an L-layer GNN, the operation in (3) is cascaded L times to obtain the GNN output Y = XL. At
the first layer, the GNN input is given by X0 = X ∈ R[n][×][F][0] . In this paper we assume F0 = FL = 1
so that Y = y and X = x. A more concise representation of this GNN can be obtained by grouping
all learnable parameters Hlk in a tensor H = {Hlk}l,k and defining the map y = Φ(x; H, S). Due
to the polynomial nature of the graph convolution, the dimensions of the learnable parameter tensor
_H are independent from the size of the graph (K is typically much smaller than n). Ergo, a GNN_
trained on a graph Gn can be deployed on a network Gm with m ̸= n.

3.2 GRAPHON INFORMATION PROCESSING

A graphon is a bounded, symmetric, and measurable function W : [0, 1][2] _→_ [0, 1] which has two
theoretical interpretations — it is both a graph limit and a generative model for graphs. In the first
interpretation, sequences of dense graphs converge to a graphon in the sense that the densities of
adjacency-preserving graph motifs converge to the same densities on the graphon Lov´asz (2012). In
the second, graphs can be generated from a graphon by sampling points ui, uj from the unit interval
and either assigning weight W(ui, uj) to edges (i, j), or sampling edges (i, j) with probability
**W(ui, uj). In this paper, we focus on stochastic graphs Gn where the points ui are defined as**
_ui = (i −_ 1)/n for 1 ≤ _i ≤_ _n and where the adjacency matrix Sn is sampled from W as_

[Sn]ij Bernoulli(W(ui, uj)). (4)
_∼_

Sequences of graphs generated in this way can be shown to converge to W with probability one
(Lov´asz, 2012)[Chapter 11].

In practice, the two theoretical interpretations of a graphon allow thinking of it as an identifying
object for a family of graphs of different sizes that are structurally similar. Hence, given a network
we can use its family’s identifying graphon as a continuous proxy for the graph. This is beneficial
because it is typically easier to operate in continuous than in discrete domains, even more so if the
network is large. We will leverage these ideas to consider graphon data and graphon neural networks
as proxies for graph data and GNNs supported on graphs of arbitrary size.


-----

3.2.1 GRAPHON NEURAL NETWORKS

Graphon data is defined as functions X ∈ _L[2]([0, 1]). Analogously to graph data, graphon data can_
be diffused by application of a linear integral operator parametrized by W and defined as

1
_TWX(v) =_ **W(u, v)X(u)du.** (5)

0

Z

The operator TW is called graphon shift operator (WSO).

The graphon convolution is defined as a weighted sum of successive applications of the WSO.
Explicitly, the graphon convolutional filter with coefficients h = [h0, . . ., hK−1] is given by

_K−1_

_Y = h ∗W X =_ _hk(TW[(][k][)][X][)(][v][)]_ with

_k=0_ (6)

X

1
(TW[(][k][)][X][)(][v][) =] **W(u, v)(TW[(][k][−][1)]X)(u)du**

0

Z

where TW[(0)] = I is the identity (Ruiz et al., 2020b). Since W is bounded and symmetric, TW
is a self-adjoint Hilbert-Schmidt operator (Lax, 2002). Hence, W can be written as W(u, v) =

_i∈Z\{0}_ _[λ][i][ϕ][i][(][u][)][ϕ][i][(][v][)][ where][ λ][i][ are the graphon eigenvalues and][ ϕ][i][ the graphon eigenfunctions.]_
_. . . λThe eigenvalues have magnitude at most one and are ordered asP_ 2 _λ1_ 1. The eigenfunctions form an orthonormal basis of −1 L ≤[2]([0λ−, 1])1 ≤. Projecting the filterλ−2 ≤ _. . . ≤_ 0 ≤
(1) onto this basis, we see that the graphon convolution admits a spectral representation given by ≤ _≤_


_K−1_

_hkλ[k]._ (7)

_k=0_

X


_h(λ) =_


Like its graph counterpart, this spectral representation only depends on the eigenvalues of the
graphon.

Graphon neural networks (WNNs) are the extension of GNNs to graphon data. In the WNN, each
layer consists of a bank of graphon convolutions (6) followed by a nonlinearity ρ. Assuming that
layer l maps Fl 1 features into Fl features, the parameters of the Fl 1 _Fl convolutions (6) can be_
stacked into K matrices− _{Hlk} ∈_ R[F][l][−][1][×][F][l] . This allows writing the− f ×th feature at layer l as

_Fl−1_ _K−1_

_Xl[f]_ [=][ ρ] (TW[(][k][)][X]l[g] 1[)[][H][lk][]][gf] (8)

 _−_ 

_g=1_ _k=1_

X X

for 1 _f_ _Fl. For an L-layer WNN, (8) is repeated for_ 1 _ℓ_ _L. The WNN output is given by_
_≤_ _≤_ _≤_ _≤_
_Y_ _[f]_ = XL[f] [, and][ X]0[g] [is given by the input data][ X] _[g][ for][ 1][ ≤]_ _[g][ ≤]_ _[F][0][. We assume][ F][L][ =][ F][0][ = 1][ so that]_
_XL = Y and X0 = X. A more succinct representation of this WNN is the map Y = Φ(X; H, W),_
where the tensor H = {Hlk}l,k groups the filter coefficients at all layers.

3.2.2 SAMPLING GNNS FROM WNNS

From the representation of the GNN and the WNN as maps Φ(x; H, S) and Φ(X; H, W), we see
that these architectures can share the same filter coefficients H. Since graphs can be obtained from
graphons as in (4), we can similarly use the WNN Φ(X; H, W) to sample GNNs
**yn = Φ(xn; H, Sn) where [Sn]ij ∼** Bernoulli(W(ui, uj)) (9)

[xn]i = X(ui)

i.e., the WNN can be seen as a generative model for GNNs Φ(xn; _, Sn)._
_H_

Conversely, a WNN ca be induced by a GNN. Given a GNN yn = Φ(xn; H, Sn), the WNN induced
by this GNN is defined as


_Yn = Φ(Xn; H, Wn) where Wn(u, v) =_



[Sn]ijI(u ∈ _Ii)I(v ∈_ _Ij)_
_j=1_

X


_i=1_


(10)



[xn]iI(u ∈ _Ii)_
_i=1_

X


_Xn(u) =_


-----

where I denotes the indicator function and the intervals Ii are defined as Ii = [(i − 1)/n, i/n) for
1 ≤ _i ≤_ _n −_ 1 and In = [(n − 1)/n, 1]. The graphon Wn is called the graphon induced by the
_graph Gn and Xn and Yn are called the graphon signals induced by the graph signals xn and yn._

4 GRAPHON EMPIRICAL LEARNING

On graphons, the statistical loss minimization (or statistical learning) problem is given by

minimize Ep(Y,X)[ℓ(Y, Φ(X; H, W))] (11)
_H_

where p(Y, X) is the joint distribution of the data, ℓ is an instantaneous loss function and
**Φ(X; H, W) is a function parametrized by the graphon W and by a set of learnable weights H.**
In this paper, we consider positive loss functions ℓ : R × R → R[+]. The function Φ is parametrized
as a graphon neural network [cf. (8)]. Because the joint probability distribution p(Y, X) is unknown,
we are unable to derive a closed-form solution of (11), but this problem can be approximated by the
empirical risk minimization (ERM) problem over graphons.

4.1 GRAPHON EMPIRICAL RISK MINIMIZATION

Suppose that we have access to samples of the distribution D = {(X _[j], Y_ _[j]) ∼_ _p(X, Y ), j =_
1, . . ., |D|}. Provided that these samples are obtained independently and that |D| is large enough,
the statistical loss in (11) can be approximated as


_|D|_

_ℓ(Y_ _[j], Φ(X_ _[j]; H, W))_ (12)
_j=1_

X


minimize


giving way to the ERM problem (Hastie et al., 2009; Shalev-Shwartz & Ben-David, 2014; Vapnik,
1999; Kearns et al., 1994). To solve this problem using local information, we could adopt some
flavor of gradient descent (Goodfellow et al., 2016). The learning iteration at step k is then

_k+1 =_ _k_ _ηk_ _ℓ(Y_ _[j], Φ(X_ _[j];_ _k, W))_ (13)
_H_ _H_ _−_ _∇H_ _H_

where k = 1, 2, . . . denotes the current iteration and ηk (0, 1) the step size at iteration k.
_∈_

In practice, the gradients on the right hand side of (13) cannot be computed because, being a theoretical limit object, the graphon W is unknown. However, we can leverage the fact that the graphon is a
random graph model to approximate the gradients _ℓ(Y_ _[j], Φ(X_ _[j];_ _, W)) by sampling stochastic_
_∇H_ _H_
graphs Gn with GSO Sn [cf. (4)] and calculating ∇Hℓ(yn[j] _[,][ Φ][(][x][j]n[;][ H][,][ S][n][))][, where][ x][j]n_ [and][ y]n[j] [are]
as in (9). In this case, the graphon empirical learning step in (13) becomes

_k+1 =_ _k_ _ηk_ _ℓ(yn[j]_ _[,][ Φ][(][x]n[j]_ [;][ H][k][,][ S][n][))] (14)
_H_ _H_ _−_ _∇H_

and we have to derive an upper bound for the expected error made when using the gradient calculated
on the graph to approximate the gradient on the graphon. In what follows, we give a closed-form
expression of this bound, and use it as a stepping stone to develop an Algorithm that increases n,
the graph size, over regular intervals during the training process of the GNN. We then prove that our
Algorithm converges to a neighborhood of the optimal solution for the WNN.

4.2 GRADIENT APPROXIMATION

To state our first convergence result, we need following definition.
**Definition 1 (Lipschitz functions). A function f** (u1, u2, . . ., ud) is A-Lipschitz on the variables
_u1, . . ., ud if it satisfies_ _f_ (v1, v2, . . ., vd) _f_ (u1, u2, . . ., ud) _A_ _i=1_
_say that this function is normalized Lipschitz. |_ _−_ _| ≤_ _[|][v][i][ −]_ _[u][i][|][. If][ A][ = 1][, we]_

[P][d]

We also need the following Lipschitz continuity assumptions (AS1–AS4), as well as an assumption
on the size of the graph Sn (AS5).
**AS1. The graphon W and the graphon signals X and Y are normalized Lipschitz.**
**AS2. The convolutional filters h are normalized Lipschitz and non-amplifying, i.e., ∥h(λ)∥** _< 1._


-----

**Algorithm 1 Increase and Conquer: Growing Graph Training**

1: Initialize H0, n0 and sample graph Gn0 from graphon W
2: repeat for epochs 0, 1, . . .
3: **for k =1,..., |D| do**

4: Obtain sample (Y, X) ∼D

5: Construct graph signal yn, xn [cf. (9)]

6:7: **end forTake learning step Hk+1 = Hk −** _ηk∇ℓ(yn, Φ(xn; Hk, Sn))_

8: Increase number of nodes n and sample Sn Bernoulli from graphon W [cf. (4)]

9: until convergence


**AS3. The activation functions and their gradients are normalized Lipschitz, and ρ(0) = 0.**

**AS4. The loss function ℓ** : R _×_ R → R[+] _and its gradient are normalized Lipschitz, and ℓ(x, x) = 0._

**AS5. For a fixed value of ξ ∈** (0, 1), n is such that n − log(2 1n/ξ)/dW > 2/dW where dW denotes
_the maximum degree of the graphon W, i.e., dW = maxv_ 0 **[W][(][u, v][)][du][.]**
R

AS1–AS4 are normalized Lipschitz smoothness conditions which can be relaxed by making the
Lipschitz constant greater than one. AS3 holds for most typical activation functions. AS4 can be
achieved by normalization and holds for most loss functions in a closed set (e.g., the hinge or mean
square losses). AS5 is necessary to guarantee a O( log n/n) rate of convergence of Sn to W

(Chung & Radcliffe, 2011).

p

Under AS1–AS5, the following theorem shows that the expected norm of the difference between the
graphon and graph gradients in (13) and (14) is bounded. The proof is deferred to the appendices.

**Theorem 1. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1 ≤_ _l ≤_ _L −_ 1. Let c ∈ (0, 1] and assume that the graphon
_convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN_
_H_
_sampled from Φ(X; H, W) as in (9). Under assumptions AS1–AS5, it holds that_

log(n[3][/][2])

E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ ] _γc +_ (15)
_∥∇H_ _H_ _−∇H_ _H_ _∥_ _≤_ _O_ _n_
[r] 

_where Yn is the graphon signal induced by [yn]i = Y (ui) [cf. (10)], and γ is a constant that_
_depends on the number of layers L, features F_ _, and filter taps K of the GNN [cf. Definition 7 in the_
_Appendix]._

The bound in Theorem 1 quantifies the maximum distance between the learning steps on the graph
and on the graphon as a function of the size of the graph. This bound is controlled by two terms. On
the one hand, we have the approximation bound, a term that decreases with n. The approximation
bound is related to the approximation of W by Sn. On the other, we have the nontransferable
_bound, which is constant and controlled by c. This term is related to a threshold that we impose on_
the convolutional filter h. Since graphons W have an infinite spectrum that accumulates around zero,
convergence of all spectral components of the filtered data on Wn can only be shown in the limit of n
(Ruiz et al., 2020a). Hence, given that n is finite we can only upper bound distances between spectral
components associated with eigenvalues larger than some c ∈ (0, 1] that we fix. Meanwhile, the
perturbations associated with the other spectral components are controlled by bounding the variation
of the filter response below c.

The bound in Theorem 1 is important because it allows us to quantify the error incurred by taking
gradient steps not on the graphon, but on the graph data. Although we cannot take gradients on the
function that we want to learn [cf. (13)], by measuring the ratio between the norm of the gradient
of the loss on the GNN, and the difference between the two gradients, we can expect to follow the
direction of the gradient on the WNN. This is instrumental for learning a meaningful solution of the
graphon ERM problem (12). Nonetheless, we are only able to decrease this approximation bound
down to the value of the nontransferable bound.


-----

4.3 ALGORITHM CONSTRUCTION

Since the discretization error depends on the number of nodes of the graph, we can iteratively increase the graph size at every epoch. Even if a bias is introduced in the gradient we will be able
to follow the gradient direction of the graphon learning problem (13). At the same time, we need
to keep the computational cost of the iterates under control. Note that the norm of the gradient is
larger at first, but it decreases as we approach the optimal solution. Exploiting this behavior, we may
progressively reduce the bias term as iterations increase. The idea here is to keep the discretization
of the gradient small so as to closely follow the gradient direction of the graphon learning problem,
but without being too conservative as this would incur in high computational cost.

In practice, in the ERM problem the number of nodes that can be added is upper bounded by the
available data, which is defined nodewise on the graph. Hence, we can arbitrarily decide which and
how many nodes to consider for training. The novelty of Algorithm 1 is that we do not use the largest
graph available in the dataset at every epoch to train the GNN. Instead, we set a minimum graph size
_n0 and progressively increase it up to the total number of nodes. The main advantage of Algorithm_
1 is thus that it allows reducing the computational cost of training without compromising GNN
performance. In what follows, we will provide the conditions under which Algorithm 1 converges
to a neighborhood of the optimal solution on the graphon.

4.4 ALGORITHM CONVERGENCE

We have shown that the learning step on the graphon and on the graph are close. Now, it remains
to show the practical implications of Theorem 1 for obtaining the solution of the graphon ERM
problem (12). If the expected difference between the gradients is small, the iterations generated by
(14) will be able to follow the direction of the true gradient on the graphon learning problem. But
because the distance between gradients is inversely proportional to n, we need to strike a balance
between obtaining a good approximation and minimizing the computational cost. In Theorem 2,
we show that Algorithm 1 converges if the rate at which the graph grows is chosen to satisfy the
condition given in (16).
**AS6. The graphon neural network Φ(X;** _, W) is AΦ-Lipschitz, and its gradient_ **Φ(X;** _, W)_
_H_ _∇H_ _H_
_is A∇Φ-Lipschitz, with respect to the parameters H [cf. Definition 1]._
**Theorem 2. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1_ _l_ _L_ 1. Let c (0, 1], ϵ (0, 1 _A_ **Φη), and assume that**
_≤_ _≤_ _−_ _∈_ _∈_ _−_ _∇_
_the graphon convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn)_
_H_
_be a GNN sampled from Φ(X; H, W) as in (9). Consider the iterates generated by equation (16)._
_Under Assumptions AS1-AS6, if at each step k the number of nodes n verifies_

log(n[3][/][2])

_γc + O_ _n_ _<_ [1][ −] _[A][∇]2[ℓ][η][ −]_ _[ϵ]_ _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥_ (16)
[r] 

_then in finite time we will achieve an iterate k[∗]_ _such that the coefficients_ _k∗_ _satisfy_
_H_

E[∥∇Hℓ(Y, Φ(X; Hk∗ _, W))∥] ≤_ 2γc _after at most k[∗]_ = O(1/ϵ) (17)

_where γ is a constant that depends on the number of layers L, features F_ _, and filter taps K of the_
_GNN [cf. Definition 7 in the Appendix]._

Theorem 2 presents the conditions under which Algorithm 1 converges. Intuitively, the condition in
(16) implies that the rate of decrease in the norm of the gradient of the loss function computed on the
graph neural network needs to be slower than roughly n[−][1][/][2]. If the norm of the gradient does not
vary between iterations, the number of nodes does not need to be increased. Note that (115) is equal
to twice the bias term obtained in equation (15). We can keep increasing the number of nodes — and
thus decreasing the bias — until the norm of the GNN gradient is smaller than the nontransferable
constant value. Once this point is attained, there is no gain in decreasing the approximation bound
any further (Ajalloeian & Stich, 2020). Recall that the constant term can be made as small as desired
by tuning c, at the cost of decreasing the approximation bound. To conclude, observe that assuming
smoothness of the GNN is a mild assumption (Scaman & Virmaux, 2018; Jordan & Dimakis, 2020;
Latorre et al., 2020; Fazlyab et al., 2019; Tanielian & Biau, 2021; Du et al., 2019). A characterization
of the Lipschitz constant is an interesting research question but is out of the scope of this work.


-----

5 NUMERICAL RESULTS

5.1 RECOMMENDATION SYSTEM


We consider a social recommendation problem given by movie ratings, for which we use the MovieLens 100k dataset Harper & Konstan (2015). The dataset contains 100, 000 integer ratings between
1 and 5, that were collected between U = 943 users and M = 1682 movies. We consider the problem of predicting the rating that different users would give to the movie “Contact”. To exploit the
social connectivity of the problem, we build the movie similarity graph, by computing the pairwise
correlations between the different movies in the training set Huang et al. (2018). Further details
about training splits, and hyperparameter selection can be found in the supplementary material.

For the experiment, we started the GNNs with 200, and at 300 nodes and added {25, 50, 75, 100}
nodes per epochs. The total number of epochs of the GNN trained on 1000 nodes with which
we compare, is given by the maximum number of epochs that the Algorithm 1 can be run adding
25 nodes per epoch. In Figure 1, we are able to see the empirical manifestation of the benefit of
Algorithm 1. Namely, regardless of the number of nodes added per epoch, in all four cases, using
relative RMSE as a figure of merit, we achieve a comparable performance to that of a GNNs trained
on 1000 nodes for a larger number of epochs. This reinforces the fact that training on a growing
number of graphs can attain similar performance, while requiring a lighter computational cost.

5.2 DECENTRALIZED CONTROL


In this section we consider the problem of coordinating a set of n agents initially flying at random to
avoid collisions, and to fly at the same velocity. Also known as flocking, at each time t agent i knows
its own position ri(t) ∈ R[2], and speed vi(t) ∈ R[2], and reciprocally exchanges it with its neighboring
nodes if a communication link exists between them. Links are govern by physical proximity between
agents forming a time varying graph Gn = (V, E). A communication link exists if the distance
between two agents i, j satisfies rij(t) = _ri(t)_ _rj(t)_ _R = 2m. We assume that at each_
_∥_ _−_ _∥≤_
time t the controller sets an acceleration ui [ 10, 10][2], and that it remains constant for a time
interval Ts = 20ms. The system dynamics are govern by, ∈ _−_ _ri(t + 1) = ui(t)Ts[2][/][2 +][ v][i][(][t][)][T][s]_ [+][ r][i][(][t][)][,]
_vi(t + 1) = ui(t)Ts + vi(t). To avoid the swarm of robots to reach a null common velocity, we_
initialize the velocities at random v(t) = [v1(t), . . ., vn(t)], by uniformly sampling a common
bias velocity vBIAS [ 3, 3], and then adding independent uniform noise [ 3, 3] to each
agent. The initial deployment is randomly selected in a circle always verifying that the minimum ∼U _−_ _U_ _−_
distance between two agents is larger than 0.1m. On the one hand, agents pursue a common average
velocity ¯v = (1/n) _i=1_ _[v][i][(][t][)][, thus minimizing the velocity variation of the team. On the other]_

[P][n]


200.0

10[2]

90.0

80.0

70.0

60.0


200.0

10[2]

90.0

80.0

70.0

60.0

|Col1|Col2|Col3|Col4|Col5|Col6|25 50 75 100|Nodes Nodes Nodes Nodes|Added Added Added Adde|(Starti (Starti (Starti d (Star|ng at 2 ng at 2 ng at 2 ting at|00 No 00 No 00 No 200 N|des) des) des) odes)|
|---|---|---|---|---|---|---|---|---|---|---|---|---|
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||
||||||||||||||

|Col1|Col2|Col3|25 Nodes A 50 Nodes A 75 Nodes A 100 Nodes|dded (Starting a dded (Starting a dded (Starting a Added (Starting|t 300 Nodes) t 300 Nodes) t 300 Nodes) at 300 Nodes)|
|---|---|---|---|---|---|
|||||||
|||||||
|||||||
|||||||
|||||||
|||||||
|||||||
|||||||
|||||||


25 Nodes Added (Starting at 200 Nodes)

50 Nodes Added (Starting at 200 Nodes)

75 Nodes Added (Starting at 200 Nodes)

100 Nodes Added (Starting at 200 Nodes)


25 Nodes Added (Starting at 300 Nodes)

50 Nodes Added (Starting at 300 Nodes)

75 Nodes Added (Starting at 300 Nodes)

100 Nodes Added (Starting at 300 Nodes)


10 11 12 13


Number of Epochs

(a)


Number of Epochs

(b)


Figure 1: Relative RMSE of a recommendation system trained on MovieLens 100k for 20 independent partitions measured over the test set (a) starting with GNNs of 200 nodes, compared to a GNN trained on 1000
nodes for 37 epochs (b) starting with GNNs of 300 nodes, compared to a GNN trained on 1000 nodes for 29
epochs.


-----

20.0

10[1]

9.0
8.0
7.0
6.0

5.0

4.0

3.0

2.0


20.0

10[1]

9.0
8.0
7.0
6.0

5.0

4.0

3.0

2.0

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|Col8|Col9|Col10|Col11|Col12|Col13|Col14|Col15|Col16|Col17|Col18|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
||||||||||2 n|ode|s pe|r ep|och|||||
|||||||||||||||||||
||||||||||5 n 10 Tr Tr|ode nod aine aine|s pe es p d on d on|r ep er e 100 100|och poch nod nod|es ( es (|10 e 30 e|poch poch|s) s)|
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||
|||||||||||||||||||

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|Col8|Col9|2 no|des p|er e|poch|Col14|Col15|Col16|Col17|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
||||||||||||||||||
||||||||||5 no 10 n Train Train|des p odes ed o ed o|er e per n 10 n 10|poch epoc 0 no 0 no|h des des|(9 ep (30 e|och poc|s) hs)|
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||
||||||||||||||||||


8 9 10 11 12 13 14 15 16 17 18

2 nodes per epoch
5 nodes per epoch
10 nodes per epoch
Trained on 100 nodes (10 epochs)
Trained on 100 nodes (30 epochs)

Number of Epochs

(a)


7 8 9 10 11 12 13 14 15 16 17

2 nodes per epoch
5 nodes per epoch
10 nodes per epoch
Trained on 100 nodes (9 epochs)
Trained on 100 nodes (30 epochs)

Number of Epochs

(b)


Figure 2: Velocity variation of the flocking problem for the whole trajectory in the testing set relative to the
centralized controller (a) starting with 10 nodes and (b) starting with 20 nodes.

hand, agents are required to avoid collision. We can thus define the velocity variation of the team
_σv(t) =_ _i=1_ **v(t)**, and the collision avoidance potential

_[∥][v][i][(][t][)][ −]_ [¯] _∥[2]_

2 2
_CA(ri[P], rj[n]) =_ 11/R/∥rCA[2]i(t) − _rj(t)∥CA −[)]_ log(∥ri(t) − _rj(t)∥_ ) otherwise,if ∥ri(t) − _rj(t)(t)∥≤_ _RCA_ (18)


_[−]_ [log(][R][2]

with RCA = 1m. A centralized controller can be obtain by ui(t)[∗] = _n(vi_ **v¯) +**
_n_ _−_ _−_
_j=1_ _i_ _[CA][(][r][i][, r][j][)][ (Tanner et al., 2003).]_

_[∇][r]_
PExploiting the fact that neural networks are universal approximators (Barron, 1993; Hornik, 1991),

_Imitation Learning can be utilized as a framework to train neural networks from information pro-_
vided by an expert (Ross et al., 2011; Ross & Bagnell, 2010). Formally, we have a set of pairs T =
_{xm, u[∗]m[}][, m][ = 1][, . . ., M]_ [, and during training we minimize the mean square error between the]
optimal centralized controller, and the output of our GNN **u[∗]m**
the neighborhood of agent i at time t, the state of the agents ∥ x(t[−]) = [[Φ][(]x[x]([m]t)[;]1[ H], . . ., x[,][ S][)][∥][2]([. Denoting]t)n], xi(t) ∈[ N][i]R[(][t][6][)],
is given by xi(t) = _j:j∈Ni(t)[[][v][i][(][t][)][ −]_ _[v][j][(][t][)][, r][ij][(][t][)][/][∥][r][ij][(][t][)][∥][4][, r][ij][(][t][)][/][∥][r][ij][(][t][)][∥][2][]][. Note that state]_

_xi(t), gets transmitted between agents if a communication link exists between them._

[P]

In Figure 2 we can see the empirical manifestation of the claims we put forward. First and foremost,
we are able to learn a GNN that achieves a comparable performance while taking steps on a smaller
graphs. As seen in Figure 2, GNNs trained with n0 = {10, 20} agents in the first epoch and adding
10 agents per epoch (green line) are able to achieve a similar performance when reaching 100 agents
than the one they would have achieved by training with 100 agents the same number of epochs.
Besides, if we add less nodes per epoch, we are able to achieve a similar performance that we would
have achieved by training on the large network for 30 epochs.


6 CONCLUSIONS

We have introduced a learning procedure for GNNs that progressively grows the size of the graph
while training. Our algorithm requires less computational cost — as the number of nodes in the
graph convolution is smaller — than training on the full graph without compromising performance.
Leveraging transferability results, we bounded the expected difference between the gradient on the
GNN, and the gradient on the WNN. Utilizing this result, we provided the theoretical guarantees
that our Algorithm converges to a neighborhood of a first order stationary point of the WNN in finite
time. We benchmarked our algorithm on a recommendation system and a decentralized control
problem, achieving comparable performance to the one achieve by a GNN trained on the full graph.


-----

REFERENCES

Ahmad Ajalloeian and Sebastian U Stich. Analysis of sgd with biased gradient estimators. arXiv
_preprint arXiv:2008.00051, 2020._

Andrew R Barron. Universal approximation bounds for superpositions of a sigmoidal function.
_IEEE Transactions on Information theory, 39(3):930–945, 1993._

Dimitri P Bertsekas and John N Tsitsiklis. Gradient convergence in gradient methods with errors.
_SIAM Journal on Optimization, 10(3):627–642, 2000._

Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and locally
connected networks on graphs. arXiv preprint arXiv:1312.6203, 2013.

Ines Chami, Sami Abu-El-Haija, Bryan Perozzi, Christopher R´e, and Kevin Murphy. Machine
learning on graphs: A model and comprehensive taxonomy, 2021.

Zhengdao Chen, Soledad Villar, Lei Chen, and Joan Bruna. On the equivalence between graph
isomorphism testing and function approximation with gnns. arXiv preprint arXiv:1905.12560,
2019.

Zhengdao Chen, Lei Chen, Soledad Villar, and Joan Bruna. Can graph neural networks count
substructures? arXiv preprint arXiv:2002.04025, 2020.

Fan Chung and Mary Radcliffe. On the spectra of general random graphs. the electronic journal of
_combinatorics, pp. P215–P215, 2011._

Micha¨el Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on
graphs with fast localized spectral filtering. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon,
and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29. Curran
Associates, Inc., 2016.

Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent finds global
minima of deep neural networks. In International Conference on Machine Learning, pp. 1675–
1685. PMLR, 2019.

Rick Durrett. Probability: Theory and Examples. Cambridge University Press, 2019.

David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alan
Aspuru-Guzik, and Ryan P Adams. Convolutional networks on graphs for learning molecular
fingerprints. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances
_in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015._

Wenqi Fan, Yao Ma, Qing Li, Yuan He, Eric Zhao, Jiliang Tang, and Dawei Yin. Graph neural
networks for social recommendation. In The World Wide Web Conference, pp. 417–426, 2019.

Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George Pappas. Efficient
and accurate estimation of lipschitz constants for deep neural networks. Advances in Neural
_Information Processing Systems, 32:11427–11438, 2019._

Alex Fout, Jonathon Byrd, Basir Shariat, and Asa Ben-Hur. Protein interface prediction using graph
convolutional networks. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30.
Curran Associates, Inc., 2017.

Fernando Gama and Somayeh Sojoudi. Distributed linear-quadratic control with graph neural networks. arXiv preprint arXiv:2103.08417, 2021.

Fernando Gama, Antonio G Marques, Geert Leus, and Alejandro Ribeiro. Convolutional neural
network architectures for signals supported on graphs. IEEE Transactions on Signal Processing,
67(4):1034–1049, 2018.

Fernando Gama, Joan Bruna, and Alejandro Ribeiro. Stability properties of graph neural networks.
_IEEE Transactions on Signal Processing, 68:5680–5695, 2020._


-----

Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural
message passing for quantum chemistry. In Proceedings of the 34th International Conference on
_Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 1263–1272._
PMLR, 06–11 Aug 2017.

Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio. Deep learning, volume 1.
MIT press Cambridge, 2016.

M. Gori, G. Monfardini, and F. Scarselli. A new model for learning in graph domains. Proceedings.
_2005 IEEE International Joint Conference on Neural Networks, 2005., 2:729–734 vol. 2, 2005._

F Maxwell Harper and Joseph A Konstan. The movielens datasets: History and context. Acm
_transactions on interactive intelligent systems (tiis), 5(4):1–19, 2015._

Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data
_mining, inference, and prediction. Springer Science & Business Media, 2009._

Kurt Hornik. Approximation capabilities of multilayer feedforward networks. Neural networks, 4
(2):251–257, 1991.

Weiyu Huang, Antonio G Marques, and Alejandro R Ribeiro. Rating prediction via graph signal
processing. IEEE Transactions on Signal Processing, 66(19):5066–5081, 2018.

Matt Jordan and Alexandros G Dimakis. Exactly computing the local lipschitz constant of relu
networks. arXiv preprint arXiv:2003.01219, 2020.

Michael J Kearns, Umesh Virkumar Vazirani, and Umesh Vazirani. An introduction to computational
_learning theory. 1994._

Nicolas Keriven and Gabriel Peyr´e. Universal invariant and equivariant graph neural networks. In
_Advances in Neural Information Processing Systems (NeurIPS), 2019._

Nicolas Keriven, Alberto Bietti, and Samuel Vaiter. Convergence and stability of graph convolutional networks on large random graphs. arXiv preprint arXiv:2006.01868, 2020.

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. _CoRR,_
abs/1412.6980, 2015.

Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907, 2016.

Fabian Latorre, Paul Rolland, and Volkan Cevher. Lipschitz constant estimation of neural networks
via sparse polynomial optimization. arXiv preprint arXiv:2004.08688, 2020.

P. D. Lax. Functional Analysis. Wiley, 2002.

Ron Levie, Wei Huang, Lorenzo Bucci, Michael M Bronstein, and Gitta Kutyniok. Transferability
of spectral graph convolutional neural networks. arXiv preprint arXiv:1907.12972, 2019.

Qingbiao Li, Fernando Gama, Alejandro Ribeiro, and Amanda Prorok. Graph neural networks for
decentralized multi-robot path planning. arXiv preprint arXiv:1912.06095, 2019.

Yujia Li, Daniel Tarlow, Marc Brockschmidt, and Richard Zemel. Gated graph sequence neural
networks. arXiv preprint arXiv:1511.05493, 2015.

L´aszl´o Lov´asz. Large networks and graph limits, volume 60. American Mathematical Soc., 2012.

Q. Lu and L. Getoor. Link-based classification. In ICML 2003, 2003.

Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman. Invariant and equivariant graph
networks, 2019.

Siyuan Qi, Wenguan Wang, Baoxiong Jia, Jianbing Shen, and Song-Chun Zhu. Learning humanobject interactions by graph parsing neural networks. In Proceedings of the European Conference
_on Computer Vision (ECCV), pp. 401–417, 2018._


-----

Meng Qu, Yoshua Bengio, and Jian Tang. Gmnn: Graph markov neural networks. In International
_conference on machine learning, pp. 5241–5250. PMLR, 2019._

St´ephane Ross and Drew Bagnell. Efficient reductions for imitation learning. In Proceedings of the
_thirteenth international conference on artificial intelligence and statistics, pp. 661–668. JMLR_
Workshop and Conference Proceedings, 2010.

St´ephane Ross, Geoffrey Gordon, and Drew Bagnell. A reduction of imitation learning and structured prediction to no-regret online learning. In Proceedings of the fourteenth international con_ference on artificial intelligence and statistics, pp. 627–635. JMLR Workshop and Conference_
Proceedings, 2011.

Luana Ruiz, Fernando Gama, Antonio Garc´ıa Marques, and Alejandro Ribeiro. Invariancepreserving localized activation functions for graph neural networks. IEEE Transactions on Signal
_Processing, 68:127–141, 2019a._

Luana Ruiz, Fernando Gama, and Alejandro Ribeiro. Gated graph convolutional recurrent neural
networks. In 2019 27th European Signal Processing Conference (EUSIPCO), pp. 1–5. IEEE,
2019b.

Luana Ruiz, Luiz Chamon, and Alejandro Ribeiro. Graphon neural networks and the transferability
of graph neural networks. Advances in Neural Information Processing Systems, 33, 2020a.

Luana Ruiz, Luiz FO Chamon, and Alejandro Ribeiro. The graphon fourier transform. In
_ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing_
_(ICASSP), pp. 5660–5664. IEEE, 2020b._

Luana Ruiz, Fernando Gama, Antonio Garc´ıa Marques, and Alejandro Ribeiro. Invariancepreserving localized activation functions for graph neural networks. IEEE Transactions on Signal
_Processing, 68:127–141, 2020c. doi: 10.1109/TSP.2019.2955832._

Luana Ruiz, Zhiyang Wang, and Alejandro Ribeiro. Graph and graphon neural network stability.
_arXiv preprint arXiv:2010.12529, 2020d._

Kevin Scaman and Aladin Virmaux. Lipschitz regularity of deep neural networks: analysis and
efficient estimation. In Proceedings of the 32nd International Conference on Neural Information
_Processing Systems, pp. 3839–3848, 2018._

Franco Scarselli, Ah Chung Tsoi, and Markus Hagenbuchner. The vapnik–chervonenkis dimension
of graph and recursive neural networks. Neural Networks, 108:248–259, 2018. ISSN 0893-6080.

Michael Schlichtkrull, Thomas N. Kipf, Peter Bloem, Rianne van den Berg, Ivan Titov, and Max
Welling. In The Semantic Web - 15th International Conference, ESWC 2018, Proceedings, pp.
593–607. Springer/Verlag, 2018.

Shai Shalev-Shwartz and Shai Ben-David. Understanding machine learning: From theory to algo_rithms. Cambridge university press, 2014._

Gilbert Strang. Linear Algebra and Its Applications. New York: Academic Press, 1976.

Qiaoyu Tan, Ninghao Liu, Xing Zhao, Hongxia Yang, Jingren Zhou, and Xia Hu. Learning to hash
with graph neural networks for recommender systems. In Proceedings of The Web Conference
_2020, pp. 1988–1998, 2020._

Ugo Tanielian and Gerard Biau. Approximating lipschitz continuous functions with groupsort neural networks. In International Conference on Artificial Intelligence and Statistics, pp. 442–450.
PMLR, 2021.

Herbert G Tanner, Ali Jadbabaie, and George J Pappas. Stable flocking of mobile agents part i:
dynamic topology. In 42nd IEEE International Conference on Decision and Control (IEEE Cat.
_No. 03CH37475), volume 2, pp. 2016–2021. IEEE, 2003._

Vladimir Vapnik. The Nature of Statistical Learning Theory. Springer Science & Business Media,
1999.


-----

Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A
comprehensive survey on graph neural networks. IEEE Transactions on Neural Networks and
_Learning Systems, 32(1):4–24, 2021._

Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec.
Graph convolutional neural networks for web-scale recommender systems. In Proceedings of the
_24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 974–_
983, 2018.

Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang,
Changcheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. AI Open, 1:57–81, 2020. ISSN 2666-6510.


-----

A APPENDIX

B NUMERICAL RESULTS PARAMETERS

All experiments were done on a computer with 32gb of RAM, a CPU Intel Core i9-9900K
@3.60GHz x 16, a GPU GeForce RTX 2080 Ti/PCIe/SSE2, and Ubuntu 18.04.4 LTS.

B.1 RECOMMENDATION SYSTEM

We split the dataset with 90% for the training set, and 10% for the testing set, and we run 20
independent random partitions. For the optimizer, we used 5 samples for the batch size, and ADAM
algorithm Kingma & Ba (2015), with learning rate 0.005, β1 = 0.9, β2 = 0.999, and without
learning rate decay. For the loss, we used the smooth L1 loss. For the GNN, we used ReLU as
non-linearity, we considered F = 32 features, K = 5 filter taps, and L = 1 layers.

[We used the graph neural networks library available online at https://github.com/alelab-upenn/graph-](https://github.com/alelab-upenn/graph-neural-networks/blob/master/examples/movieGNN.py)
[neural-networks/blob/master/examples/movieGNN.py and implemented with PyTorch.](https://github.com/alelab-upenn/graph-neural-networks/blob/master/examples/movieGNN.py)

B.2 DECENTRALIZED CONTROL

We run the system for T = 2s, and used 400 samples for training, 20 for validation, and 20 for the
test set. For the optimizer, we used 20 samples for the batch size, and ADAM algorithm Kingma &
Ba (2015) with learning rate 0.0005, β1 = 0.9, β2 = 0.999, without learning rate decay. We used
a one layer Graph Neural Networks with F = 64 hidden units and K = 3 filter taps, and used the
hyperbolic tangent as non-linearity ρ. We run 10 independent realizations of each experiment.

[We used the graph neural networks library available online at https://github.com/alelab-upenn/graph-](https://github.com/alelab-upenn/graph-neural-networks/blob/master/examples/flockingGNN.py)
[neural-networks/blob/master/examples/flockingGNN.py and implemented with PyTorch.](https://github.com/alelab-upenn/graph-neural-networks/blob/master/examples/flockingGNN.py)


-----

C PROOF OF THEOREM 1

**Definition 2 (Template graphs). Let {ui}i[n]=1** _[be the regular][ n][-partition of][ [0][,][ 1]][, i.e.,]_

_ui =_ _[i][ −]_ [1]


(19)


_for 1 ≤_ _i ≤_ _n. The n-node template graph Gn, whose GSO we denote Sn, is obtained from W as_

[Sn]ij = W(ui, uj) (20)

_for 1 ≤_ _i, j ≤_ _n._

**Definition 3 (Graphon spectral representation of convolutional filter response). As the graphon**
**W is bounded and symmetric, TW is a self adjoint Hilbert-Schmidt operator, which allows to use**
_the operator’s spectral basis W(u, v) =_ _i_ **Z** 0 _[λ][i][ψ][i][(][u][)][ψ][i][(][v][)][. Eigenvalues][ λ][i][ are ordered in]_

_∈_ _{_ _}_
_decreasing order of absolute value i.e.,their only accumulation point is 0 (Lax, 2002, Theorem 3, Chapter 28). Thus, we define the spectral 1 ≥_ _λ1 ≥_ _λ2 ≥· · · ≥_ 0 ≥· · · ≥ _λ−2 ≥_ _λ−1 ≥−1, and_

[P]

_representation of the convolutional filter TH (cf. (6)) as,_


_K−1_

_hkλ[k]_ (21)

_k=0_

X


_h(λ) =_


**Definition 4 (c-band cardinality of W). The c-band cardinality, denoted BW[c]** _[, is the number of]_
_eigenvalues whose absolute value is larger than c._

_BW[c]_ [= #][{][λ][i] [:][ ∥][λ][i][∥≤] _[c][}]_ (22)

**Definition 5 (c-eigenvalue margin of W - Wn). The c-eigenvalue margin of W - Wn is defined as**
_the minimum distance between two different eigenvalues of the integral operator applied to W, and_
_to Wn as follows,_

_δWW[c]_ _n_ [= min]i,j=i[{∥][λ][i][(][T][W][)][ −] _[λ][i][(][T][W][n]_ [)][∥] [:][ ∥][λ][i][(][T][W][n] [)][∥≥] _[c][}]_ (23)
_̸_

**Definition 6 (Graphon Convolutional Filter). Given a graphon W, a graphon signal X, and filter**
_coefficients h = [h0, . . ., hK_ 1] the graphon filter TH : L2([0, 1]) _L2([0, 1]) is defined as,_
_−_ _→_


_K−1_

_hk(TW[(][k][)][X][)(][v][)][.]_ (24)
_k=0_

X


(THX)(v) =


**Proposition 1. Let X ∈** _L2([0, 1]) be a normalized Lipschitz graphon signal, and let Xn be the_
_graphon signal induced by the graph signal xn obtained from X on the template graph Gn [cf._
_Definition 2], i.e., [xn]i = X((i −_ 1)/n) for 1 ≤ _i ≤_ _n. It holds that_

_∥X −_ _Xn∥L2 ≤_ _n[1]_ _[.]_ (25)

_Proof. Let Ii = [(i −_ 1)/n, i/n) for 1 ≤ _i ≤_ _n −_ 1 and In = [(n − 1)/n, 1]. Since the graphon is
normalized Lipschitz, for any u _Ii, 1_ _i_ _n, we have_
_∈_ _≤_ _≤_

_i_

_X(u)_ _Xn(u)_ max _[,]_ (26)
_∥_ _−_ _∥≤_ _n_ _n_ _≤_ _n[1]_ _[.]_
 _[−]_ _[u]_ 

We can then write

_[u][ −]_ _[i][ −]_ [1]

1
_X_ _Xn_ = _X(u)_ _Xn(u)_ _du_ (27)
_∥_ _−_ _∥[2]_ 0 _|_ _−_ _|[2]_
Z

1 1 2 1 2

_du =_ _,_ (28)

_≤_ 0 _n_ _n_
Z    

which completes the proof.


-----

**Proposition 2. Let W : [0, 1][2]** _→_ [0, 1] be a normalized Lipschitz graphon, and let Wn := WGn
_be the graphon induced by the template graph Gn generated from W as in Definition 2. It holds_
_that_
**W** Wn (29)
_∥_ _−_ _∥≤_ _n[2]_ _[.]_

_Proof. Let Ii = [(i −_ 1)/n, i/n) for 1 ≤ _i ≤_ _n −_ 1 and In = [(n − 1)/n, 1]. Since the graphon is
Lipschitz, for any u _Ii, v_ _Ij, 1_ _i, j_ _n, we have_
_∈_ _∈_ _≤_ _≤_

_i_

**W(u, v)** Wn(u, v) max _[,]_ (30)
_∥_ _−_ _∥≤_ _n_ _n_
 _[−]_ _[u]_ 

_j_

+ max _[u][ −]_ _[i][ −]_ [1] _[,]_ (31)

_n_ _n_

 _[−]_ _[v]_ 

(32)

_≤_ _n[1]_ [+ 1]n _[v][= 2][ −]n[j][.][ −]_ [1]


We can then write

1
_∥W −_ Wn∥[2] = 0
Z

1
_≤_ 0
Z

which concludes the proof.


_|W(u, v) −_ Wn(u, v)|[2]dudv (33)

2 2

2 2

_dudv =_ (34)

_n_ _n_

   


**Proposition 3. Consider the L-layer WNN given by Y = Φ(X; H, W), where F0 = FL = 1 and**
_Fℓ_ = F for 1 _ℓ_ _L_ 1. Let c (0, 1] and assume that the graphon convolutions in all layers of
_≤_ _≤_ _−_ _∈_
_this WNN have K filter taps [cf. (6)]. Under Assumptions 1 through 3, the norm of the gradient of_
_the WNN with respect to its parameters H = {Hlk}l,k can be upper bounded by,_

**Φ(X;** _, W)_ _F_ [2][L][√]K. (35)
_∥∇H_ _H_ _∥≤_

_Proof. We will find an upper bound for any element [Hl†k†_ ]g†f † of the tensor H. We start by the
last layer of the WNN, applying the definition given in equation (8),

_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥_ = _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _[X]L[f]_ (36)

_Fl−1_ _K−1_

= _[∇][[][H]l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _[ρ]_ (TW[(][k][)][X]l[g] 1[)[][H][Lk][]][gf] (37)

 _−_ 

_g=1_ _k=1_

X X

 

By Assumption 3, the non-linearity ρ is normalized Lipschitz, i.e.[∇][[][H] _ρ(_ )(u) 1 for all[.] u. Thus,
_∇_ _·_ _≤_
aplying the chain rule for the derivative, and the Cauchy-Schwartz inequality, the right hand side of
the previous expression can be rewritten as,


_Fl−1_

_g=1_

X


_K−1_

(TW[(][k][)][X]l[g] 1[)[][H][Lk][]][gf]
_−_
_k=1_

X


_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥_ = (TW[(][k][)][X]l[g]−1[)[][H][Lk][]][gf]

_g=1_ _k=1_

X X

_[∇][ρ]_  _Fl−1_ _K−1_ 

_l[†]_ _k[†]_ []]g[†] _f_ _[†]_ (TW[(][k][)][X]l[g]−1[)[][H][Lk][]][gf] (38)

_g=1_ _k=1_

X X

_[∇][[][H]_ _Fl−1_ _K−1_

_l[†]_ _k[†]_ []]g[†] _f_ _[†]_ (TW[(][k][)][X]l[g] 1[)[][H][Lk][]][gf] (39)
_≤_ _−_

_g=1_ _k=1_

X X

Note that the a larger bound will occur if l[†] _< L[∇][[][H]_ 1, then by linearity of derivation, and the triangle
_−_
inequality we obtain,


_Fl−1_

_g=1_

X


_K−1_

_TW[(][k][)][(][∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]l[g]−1[)[][H][Lk][]][gf]_
_k=1_

X


(40)


_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥≤_


-----

By Assumption 2, the convolutional filters are non-amplifying, thus it holds that,


_Fl−1_

_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥≤_ _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _[X]l[g]−1_ (41)

_g=1_

X

Now note that as filters are non-amplifying, the maximum difference in the gradient will be attained[∇][[][H]
at the first layer (l = 1) of the WNN. Also note that the derivative of a convolutional filter TH [cf.
Definition 6] at coefficient k[†] = i, is itself a convolutional filter with coefficients hi. The values of
**hi are [hi]j = 1 if j = i and 0 otherwise. Thence,**

_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥≤_ _F_ _[L][−][1]_ (42)

_F_ _[L][−][1]_ _X0_ _._ (43)
_≤_ _∥[h][i][∗]∥[W][X][0]_

To complete the proof note that tensor H has F _[L][−][1]K elements, and each individual gradient is_
upper bounded by (43), and ∥X∥ is normalized by Assumption 1.

**Lemma 1. Let Φ(X; H, W) be a WNN with F0 = FL = 1, and Fl = F for 1 ≤** _l ≤_ _L −_ 1. Let
_c ∈_ (0, 1], and assume that the graphon convolutions in all layers of this WNN have K filter taps

_[cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN sampled from Φ(X;_ _, W) as in (9). Under assumptions_
_H_ _H_
(1),(2),(3), and (5) with probability 1 − _ξ it holds that,_


2 1 +




_n log(_ [2]ξ[n] [)]


**Φ(X;** _, W)_ **Φ(Xn;** _, Wn)_ _LF_ _[L][−][1]_ 1 + _[πB]W[c]_ _n_
_∥_ _H_ _−_ _H_ _∥≤_ _δWW[c]_ _n_



+ [1] (44)

_n_ [+ 4][LF][ L][−][1][c]


_The fixed constants BW[c]_ _[and][ δ]WW[c]_ _n_ _[are the][ c][-band cardinality and the][ c][-eigenvalue margin of][ W]_
_and Wn respectively [cf. Definitions 4,5]._

_Proof. We start by writing the expression on the left hand side, using the definition of WNN [cf._
(8)] we can write,

**Φ(X;** _, W)_ **Φ(Xn;** _, Wn)_ = _XL_ _XnL_ (45)
_∥_ _H_ _−_ _H_ _∥_ _∥_ _−_ _∥_

_FL−1_ _K−1_ _FL−1_ _K−1_

= _ρ_ (TW[(][k][)][X]L[g] 1[)[][H][Lk][]][gf] _ρ_ (TW[(][k][)]n _[X]nL[g]_ 1[)[][H][Lk][]][gf] _._

 _−_  _−_  _−_ 

_g=1_ _k=1_ _g=1_ _k=1_

X X X X

   

Since the non-linearity ρ is normalized Lipschitz by Assumption 3, using the triangle inequality, we
obtain


_FL−1_

_g=1_

X


_K−1_

(TW[(][k][)][X]L[g] 1[)[][H][Lk][]][gf][ −]
_−_
_k=1_

X


_K−1_

(TW[(][k][)]n _[X]nL[g]_ 1[)[][H][Lk][]][gf]
_−_
_k=1_

X


(46)


_XL_ _XnL_
_∥_ _−_ _∥≤_


Using the triangle inequality once again, we split the last inequality into two terms as follows,


_FL−1_

_g=1_

X


_K−1_

_TW[(][k][)][(][X]L[g]_ 1 _nL_ 1[)[][H][Lk][]][gf]
_k=1_ _−_ _[−]_ _[X]_ _[g]_ _−_

X


_XL_ _XnL_ _TW[(][k][)][(][X]L[g]_ 1 _nL_ 1[)[][H][Lk][]][gf]
_∥_ _−_ _∥≤_ _g=1_ _k=1_ _−_ _[−]_ _[X]_ _[g]_ _−_ **[(1)]**

X X

_FL−1_ _K−1_

+ (TW[(][k][)] **Wn** [)][X]L[g] 1[[][H][Lk][]][gf] (47)

_g=1_ _k=1_ _[−]_ _[T][ (][k][)]_ _−_ **[(2)][.]**

X X

Where we have split (47) into terms (1), and (2). On the one hand, by assumption 2, convolutional
filters h are non-amplifying, thus using Cauchy-Schwartz inequality, term (1) can be bounded by,


_FL−1_

_g=1_

X


_K−1_

_TW[(][k][)][(][X]L[g]_ 1 _nL_ 1[)[][H][Lk][]][gf]
_k=1_ _−_ _[−]_ _[X]_ _[g]_ _−_

X


_FL−1_

_g=1_

X


_XL[g]_ _−1_ _[−]_ _[X]nL[g]_ _−1_ _._ (48)


-----

To bound term (2), denoting hLgf the spectral representation of the convolutional filter applied to
_XL[g]_ 1 [at feature][ f][ of layer][ L][ [cf. Definition 3], we will decompose the filter as follows,]
_−_

0 if _λ_ _< c_
_h[≥]Lgf[c]_ [(][λ][)] _|_ _|_ (49)
_hLgf_ (λ) _hLgf_ (c) if _λ_ _c_
 _−_ _|_ _| ≥_

_h[<c]Lgf_ [(][λ][)] _hLgf_ (λ) if |λ| < c (50)
_hLgf_ (c) if _λ_ _c._
 _|_ _| ≥_

Note that hLgf = h[≥]Lgf[c] [+][ h]Lgf[<c] [. Let][ T][ <c][HL]gf [and][ T][ <c][HnL]gf [, be the graphon convoutional filters]
with filter function h[<c]Lgf [on graphons][ W][, and][ W][n][ respectively [cf. Definition 6]. Note that filter]
_h[<c]Lgf_ [, varies only in the interval][ [0][, c][)][, and since filters are normalized Lipschitz by Assumption 2,]
it verifies
_T[[<c]HL]gf_ _[X]L[g]_ _−1_ _[−]_ _[T][ <c][HnL]gf_ _[X]L[g]_ _−1_ _≤_ (hLgf (c) + c) − (hLgf (c) − _c)∥∥XL[g]_ _−1_ (51)

2c _XL[g]_ 1[∥][.] (52)
_≤_ _∥_ _−_

Now we need to upper bound the difference in the high frequencies h[≥]Lgf[c] [. Let][ T][ ≥][H[c]L]gf [and][ T][ ≥][H[c]nL]gf [,]

be the graphon filters with filter function h[≥]Lgf[c] [on graphons][ W][, and][ W][n][ respectively. Let][ S][n]
denote the template graph sampled from the graphon W [cf. definition 2]. We denote Wn, the
induced graphon by template graph Sn as in (10). By introducing T[[≥]H[c]nL]gf [, the graph filter with]

filter function h[≥]Lgf[c] [on graphon][ W][n][, we can use the triangle inequality to obtain,]
_T[[≥]H[c]L]gf_ _[X]L[g]_ _−1_ _[−]_ _[T][ ≥][H[c]nL]gf_ _[X]L[g]_ _−1_ _≤_ _T[[≥]H[c]L]gf_ _[X]L[g]_ _−1_ _[−]_ _[T][ ≥][H[c]nL]gf_ _[X]L[g]_ _−1_ **(2.1)**

+ _T[[≥]H[c]nL]gf_ _[X]L[g]_ _−1_ _[−]_ _[T][ ≥][H[c]nL]gf_ _[X]L[g]_ _−1_ **(2.2).** (53)

Under assumptions 1–5, to bound term (2.1) we can use (Ruiz et al., 2020a, Theorem 1), and to
bound term (2.2) we can use (Ruiz et al., 2020d, Lemma 2). Thus, with probability 1 − _ξ, the_
previous expression can be bounded by,


_n log(_ [2]ξ[n] [)]


2 1 +




2 1 + _n log(_ _ξ_ [)]

**Wn**
_T[[≥]H[c]L]gf_ _[X]L[g]_ _−1_ _[−]_ _[T][ ≥][H[c]nL]gf_ _[X]L[g]_ _−1_ _≤_ 1 + _δ[πB]WW[c]_ _[c]_ _n_   qn  _XL[g]_ _−1_ _._ (54)

Where the fixed constants BW[c] [and][ δ]WW[c] _n_ [are the][ c][-band cardinality and the][ c][-eigenvalue margin]
of W and Wn respectively [cf. Definitions 4,5]. Hence, coming back to (47), we can use (48) to
upper bound (1), and we can use (52), and (54), to upper bound (2) as follows,


_FL−1_

_g=1_

X


_XL[g]_ 1 _nL_ 1 + 2c _XL[g]_ 1[∥]
_−_ _[−]_ _[X]_ _[g]_ _−_ _∥_ _−_

2 1 + _n log(_ [2]ξ[n] [)]

**Wn**
1 + _[πB][c]_  q

_δWW[c]_ _n_ _n_

 


_XL_ _XnL_
_∥_ _−_ _∥≤_


+ 1 + _δ[πB]WW[c]_ **Wnn** _n_ _XL[g]_ _−1_ _._ (55)
 

Now, we arrive at a recursive equation that we can compute for the L layers, with F features per
layer, to obtain,


_XL_ _XnL_ _F0_ _X0_ _Xn0_ + 2LF _[L][−][1]c_ _X0_
_∥_ _−_ _∥≤_ _∥_ _−_ _∥_ _∥_ _∥_


_n log(_ [2]ξ[n] [)]


2 1 +




+ LF _[L][−][1]_ 1 + _[πB]W[c]_ _n_

_δWW[c]_ _n_




_X0_ _._ (56)
_∥_ _∥_


Using Proposition 1, noting that F0 = 1 by construction, and using Assumption 1, concludes the
proof.


-----

**Lemma 2. Let Φ(X; H, W) be a WNN with F0 = FL = 1, and Fl = F for 1 ≤** _l ≤_ _L −_ 1. Let
_c ∈_ (0, 1], and assume that the graphon convolutions in all layers of this WNN have K filter taps

_[cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN sampled from Φ(X;_ _, W) as in (9). Under assumptions_
_H_ _H_
(1),(2),(3), and (5) with probability 1 − _ξ it holds that,_

_∥∇HΦ(X; H, W) −∇HΦ(Xn; H, Wn)∥_


_n log(_ [2]ξ[n] [)]


2 1 +




_KF_ _[L][−][1]_ 2L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_

_δWW[c]_ _n_

 


+ [2][F][ L][−][1][L]


+ 8L[2]F [2][L][−][2]c


_Proof. We will first show that the gradient with respect to any arbitrary element [Hl†k†_ ]g†f † ∈ R
of H can be uniformly bounded. Note that the maximum is attained if l[†] = 1. Without loss of
generality, assuming l[†] _> l −_ 1, we can begin by using the definition given in equation (8) of the
output of the WNN as follows,

_∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W) −∇[Hl†_ _k†_ ]g† _f_ _† Φ(Xn; H, Wn)∥_

= ∥∇[Hl† _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_ (57)

_Fl−1_ _K−1_

= _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _[ρ]_ (TW[(][k][)][X]l[g] 1[)[][H][lk][]][gf]

 _−_ 

_g=1_ _k=1_

X X

_[∇][[][H]_  _Fl−1_ _K−1_ 

_−∇[Hl†_ _k†_ ]g† _f_ _† ρ_  (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_  (58)

_g=1_ _k=1_

X X

  _[.]_

Taking derivatives by applying the chain rule, and applying the triangle inequality it yields,

_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_

_Fl−1_ _K−1_ _Fl−1_ _K−1_

_≤_ _∇ρ_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf] _−∇ρ_ (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_

  _g=1_ _k=1_   _g=1_ _k=1_ 
X X X X

_Fl−1_ _K−1_

_∇[Hl†_ _k†_ ]g† _f_ _†_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf] (59)

 _g=1_ _k=1_ 
X X

_Fl−1_ _K−1_
+ (TW[(][k][)]n _[X]nl[g]_ 1[)[][H][lk][]][gf] (60)

_−_

 _g=1_ _k=1_ 
X X

_[∇][ρ]_ _Fl−1_ _K−1_ _Fl−1_ _K−1_

_∇[Hl†_ _k†_ ]g† _f_ _†_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −∇][[][H]l† _k†_ []]g† _f_ _†_ (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_
 _g=1_ _k=1_ _g=1_ _k=1_ 

X X X X

_[.]_

We can now use Cauchy-Schwartz inequality, Assumptions 3, 4, and Proposition 3 to bound the
terms regarding the gradient of the non-linearity ρ, the loss function ℓ, and the WNN respectively,
as follows,

_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_ (61)

_Fl−1_ _K−1_ _Fl−1_ _K−1_

_≤_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −] (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_

_g=1_ _k=1_ _g=1_ _k=1_

X X X X

_Fl−1_ _K−1_ _[F][ L][−][1][∥][X][0][∥]_

+ _∇[Hl†_ _k†_ ]g† _f_ _†_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −] [(][T][ (]W[k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_ [)]

_g=1_ _k=1_  

X X

_[.]_


-----

We can now apply the triangle inequality on the second term of the previous bound to obtain,

_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_ (62)

_Fl−1_ _K−1_ _Fl−1_ _K−1_

_≤_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −] (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_

_g=1_ _k=1_ _g=1_ _k=1_

X X X X

_Fl−1_ _K−1_ _[F][ L][−][1][∥][X][0][∥]_

+ _∇[Hl†_ _k†_ ]g† _f_ _†_ (TW[(][k][)][)[][H][lk][]][gf][ −] [(][T][ (]W[k][)]n [)[][H][lk][]][gf] [)] _Xnl[g]_ _−1_

_g=1_ _k=1_  

X X

_Fl−1_ _K−1_

+ _g=1_ _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _k=1_ _TW[(][k][)]n_ Xl[g]−1 _[−]_ _[X]nl[g]_ _−1[Hlk]gf_ )

X X

Now note that as we are considering the case in which l _< l_ 1, using Cauchy-Schwartz inequality,

_[∇][[][H]_ _†_ _−_ _[.]_

we can use the same bound for the first and second term of the right hand side of the previous
inequality. Since filters are non-expansive by Assumption 3, it yields

_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_ (63)

_Fl−1_ _K−1_ _Fl−1_ _K−1_

_≤_ 2 (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −] (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_

_g=1_ _k=1_ _g=1_ _k=1_

X X X X

_Fl−1_

_[F][ L][−][1][∥][X][0][∥]_

+ _g=1_ _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ Xl[g]−1 _[−]_ _[X]nl[g]_ _−1_

X

Now notice, that the only term that remains to bound is the exact same bound we obtained in equation

_[∇][[][H]_ _[.]_

(57), but on the previous layer L − 2. Hence, we conclude that by applying the same steps L − 2
times, as the WNN has L layers, we will obtain a bound for any element [Hl†k† ]g†f † of tensor H.

_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥]_ (64)

_Fl−1_ _K−1_ _Fl−1_ _K−1_

_≤_ 2LF _[L][−][2]_ (TW[(][k][)][X]l[g]−1[)[][H][lk][]][gf][ −] (TW[(][k][)]n _[X]nl[g]_ _−1[)[][H][lk][]][gf]_

_g=1_ _k=1_ _g=1_ _k=1_

X X X X

_Fl−1_

_[F][ L][−][1][∥][X][0][∥]_

+ _l[†]_ _k[†]_ []]g[†] _f_ _[†]_ _X1[g]_ 1

_g=1_  _[−]_ _[X]_ _[g]_

X

Note that the derivative of a convolutional filter[∇][[][H] _[.]_ _TH at coefficient k[†]_ = i, is itself a convolutional
filter with coefficients hi [cf. Definition 6]. The values of hi are [hi]j = 1 if j = i and 0 otherwise.
As hi is itself a filter that verifies Assumption 2, as graphons are normalized. Thus, considering l[†] =
0, and using Propositions 1, 2, (Chung & Radcliffe, 2011, Theorem 1) and the triangle inequality,
we obtain,

_n_ _[X][n][0][ −]_ **[h][i][∗][W][X][0]** **W** Wn + Wn **Wn** _X0_ + _Xn0_ _X0_ (65)

_[≤]_ ∥ _−_ _∥_ _∥_ _−_ _∥∥_ _∥_ _∥_ _−_ _∥_

**[h][i][∗][W]** 2 1 + _n log(_ [2]ξ[n] [)]

**Wn**
1 + _[πB][c]_  q  + [1] (66)
_≤_ _δWW[c]_ _n_ _n_ _n_
 

with probability 1 − _ξ. In the previous expression, Wn is the template graphon [cf. Definition 2]._
Now, substituting (64) into (65), and using Lemma 1, with probability 1 − _ξ, it holds that,_


2 1 +




_n log(_ [2]ξ[n] [)]


_∥∇[Hl†_ _k†_ ]g† _f_ _† XL[f]_ _[−∇][[][H]l†_ _k†_ []]g† _f_ _†_ _[X]nL[f]_ _[∥≤][2][L][2][F][ 2][L][−][2]1 +_ _δ[πB]WW[c]_ **W[c]** _nn_


+ [2][F][ L][−][1][L]


+ 8L[2]F [2][L][−][2]c. (67)


To achieve the final result, note that tensor H has KF _[L][−][1]_ elements, and each individual gradient is
upper bounded by (67).


-----

**Lemma 3. Let Φ(X; H, W) be a WNN with F0 = FL = 1, and Fl = F for 1 ≤** _l ≤_ _L −_ 1. Let
_c ∈_ (0, 1], and assume that the graphon convolutions in all layers of this WNN have K filter taps

_[cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN sampled from Φ(X;_ _, W) as in (9). Under Assumptions_
_H_ _H_
(1)–(5) with probability 1 − _ξ it holds that,_

_∥∇Hℓ(Y, Φ(X; H, W)) −∇Hℓ(Yn, Φ(Xn; H, Wn))∥_


_n log(_ [2]ξ[n] [)]


2 1 +




_KF_ _[L][−][1]_ 3L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_

_δWW[c]_ _n_

 


+ [4][F][ L][−][1][L] + 12L[2]F [2][L][−][2]c


_Proof. In order to analyze the norm of the gradient with respect to the tensor H, we can start by_
taking the derivative with respect to a single element of the tensor, [Hl†k† ]g†f † . By deriving the loss
function ℓ using the chain rule it yields,

_∥∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Y, Φ(X; H, W)) −∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Yn, Φ(Xn; H, Wn))∥_

=∥∇ℓ(Y, Φ(X; H, W))∇[Hl† _k†_ ]g† _f_ _† Φ(X; H, W)_

_−∇ℓ(Yn, Φ(Xn; H, Wn))∇[Hl†_ _k†_ ]g† _f_ _† Φ(Xn; H, Wn)∥._ (68)

By Cauchy-Schwartz, and the triangle inequality it holds,

_∥∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Y, Φ(X; H, W)) −∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Yn, Φ(Xn; H, Wn))∥_

_≤∥∇ℓ(Y, Φ(X; H, W)) −∇ℓ(Yn, Φ(Xn; H, Wn))∥∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥_ (69)

+ ∥∇ℓ(Yn, Φ(Xn; H, Wn))∥∥∇[Hl† _k†_ ]g† _f_ _† Φ(X; H, W) −∇[Hl†_ _k†_ ]g† _f_ _† Φ(Xn; H, Wn)∥._

By the triangle inequality and Assumption 4 it follows,

_∥∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Y, Φ(X; H, W)) −∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Yn, Φ(Xn; H, Wn))∥_

_≤∥∇ℓ(Y, Φ(X; H, W)) −∇ℓ(Y, Φ(Xn; H, Wn))∥∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥_ (70)

_∥∇ℓ(Yn, Φ(Xn; H, Wn)) −∇ℓ(Y, Φ(Xn; H, Wn))∥∥∇[Hl†_ _k†_ ]g† _f_ _† Φ(X; H, W)∥_ (71)

+ ∥∇[Hl† _k†_ ]g† _f_ _† Φ(X; H, W) −∇[Hl†_ _k†_ ]g† _f_ _† Φ(Xn; H, Wn)∥_

_≤_ (∥Yn − _Y ∥_ + ∥Φ(Xn; H, Wn)) − **Φ(X; H, W))∥)∥∇[Hl†** _k†_ ]g† _f_ _† Φ(X; H, W)∥_ (72)

+ ∥∇[Hl† _k†_ ]g† _f_ _† Φ(X; H, W) −∇[Hl†_ _k†_ ]g† _f_ _† Φ(Xn; H, Wn)∥._


Now we can use Lemmas 1–2, Propositions 1, and 3, and Assumption 1 to obtain,

_∥∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Y, Φ(X; H, W)) −∇[Hl†_ _k†_ ]g† _f_ _† ℓ(Yn, Φ(Xn; H, Wn))∥_ (73)


_n log(_ [2]ξ[n] [)]


2 1 +




3L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_

_δWW[c]_ _n_




+ [4][F][ L][−][1][L]


+ 12L[2]F [2][L][−][2]c


Noting that tensor H has KF _[L][−][1]_ elements, and each individual term can be bounded by (73), the
desired result is attained.

**Definition 7. We define the constant γ as,**


_KF_ _[L][−][1]L[2]F_ [2][L][−][2], (74)


_γ = 12_


_where K is the number of features, L is the number of layers, and K is the number of filter taps of_
_the GNN._


-----

We will present a more comprehensive statement of Theorem 1, where we include all the smaller
order terms in (15). Notice that the statement of Theorem 1 in the main body of the paper omits
these terms in order to simplify the exposition of the main result. In practice, these smaller order
terms vanish faster as n increases.
**Theorem 1. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1 ≤_ _l ≤_ _L −_ 1. Let c ∈ (0, 1] and assume that the graphon
_convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN_
_H_
_sampled from Φ(X; H, W) as in (9). Under assumptions AS1–AS5, it holds that_
E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ ]
_∥∇H_ _H_ _−∇H_ _H_ _∥_


_n log(2n[3][/][2])_


1 +


1 + _n log(2n_ )

_√KF_ _[L][−][1]_ 6L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_  
_≤_ _δWW[c]_ _n_ p _n_

  

_K_

+ [4][F][ L]n[−][1][L] + 12L[2]F [2][L][−][2]c + [2][F][ 2]√[L]n[√] (75)




_where Yn is the graphon signal induced by [yn]i = Y (ui), ui = (i −_ 1)/n for 1 ≤ _i ≤_ _n [cf. (10)]._
_The fixed constants BW[c]_ _[and][ δ]WW[c]_ _n_ _[are the][ c][-band cardinality and the][ c][-eigenvalue margin of][ W]_
_and Wn respectively [cf. Definitions 4,5 in the supplementary material]._

_Proof of Theorem 1. We begin by considering the event An such that,_

_An =_ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ (76)
_∥∇H_ _H_ _−∇H_ _H_ _∥_



_n log(_ [2]ξ[n] [)]


2 1 +




2 1 + _n log(_ _ξ_ [)]

_√KF_ _[L][−][1]_ 3L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_  q 
_≤_ _δWW[c]_ _n_ _n_

  

+ [4][F][ L][−][1][L] + 12L[2]F [2][L][−][2]c _._

_n_


Thus, by considering the disjoint events An, and A[c]n[, and denoting][ 1][(][·][)][ the indicator function, the]
expectation can be separated as follows,
E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ ]
_∥∇H_ _H_ _−∇H_ _H_ _∥_
= E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ **1(An)]**
_∥∇H_ _H_ _−∇H_ _H_ _∥_
+ E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ **1(A[c]n[)]]** (77)
_∥∇H_ _H_ _−∇H_ _H_ _∥_
We can bound the term regarding A[c]n [using the chain rule, Cauchy-Schwartz inequality, Assumption]
4, and Proposition 3 as follows,
_∥∇Hℓ(Y, Φ(X; H, W)) −∇Hℓ(Yn, Φ(Xn; H, Wn))∥_
_≤∥∇Hℓ(Y, Φ(X; H, W))∥_ + ∥∇Hℓ(Yn, Φ(Xn; H, Wn))∥ (78)
_ℓ(Y, Φ(X;_ _, W))_ **Φ(X;** _, W)_
_≤∥∇_ _H_ _∥∥∇H_ _H_ _∥_
+ ∥∇ℓ(Yn, Φ(Xn; H, Wn))∥∥∇HΦ(Xn; H, Wn)∥ (79)
_≤∥∇HΦ(X; H, W)∥_ + ∥∇HΦ(Xn; H, Wn)∥ (80)

_≤_ 2F [2][L][√]K (81)

Returning to equation (77), we can substitute the bound obtained in equation (81), and by taking
_P_ (An) = 1 _ξ, and using Lemma 3, it yields,_
_−_
E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ ]
_∥∇H_ _H_ _−∇H_ _H_ _∥_


_n log(_ [2]ξ[n] [)]


2 1 +




_KF_ _[L][−][1]_ 3L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_

_δWW[c]_ _n_

 

+ 12L[2]F [2][L][−][2]c + ξ2F [2][L][√]K



_≤_ (1 − _ξ)_

+ [4][F][ L][−][1][L]


(82)


To complete the proof, set ξ = _√1n_ .


-----

E[ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Yn, Φ(Xn;_ _, Wn))_ ]
_∥∇H_ _H_ _−∇H_ _H_ _∥_

E

_Xn;_ _, Wn))_
_H_


Figure 3: In order to satisfy the property that the inner product between the gradient on the GNN
_ℓ(Yn, Φ(Xn;_ _, Wn)) and the gradient on the graphon_ _ℓ(Y, Φ(X;_ _, W)) is positive, we rely on_
_∇H_ _H_ _∇H_ _H_
the condition provided in Theorem 2.

D PROOF OF THEOREM 2

**Definition 8 (Stopping time). We define the stopping time k[∗]** _as,_


_k[∗]_ = min
_k_ 0 _[{∥∇][H][Φ][(][X][;][ H][k][,][ W][n][)][∥≤]_
_≥_

**Definition 9 (Constant ψ).**


_KF_ _[L][−][1]12L[2]F_ [2][L][−][2]c}. (83)


**Lemma 4. Under Assumptions 4, 5, and 6, the gradient of the loss function ℓ** _with respect to the_
_parameters of the GNN H is A∇ℓ-Lipschitz,_

_∥∇Hℓ(Y, Φ(X; A, W)) −∇Hℓ(Y, Φ(X; B, W))∥≤_ _A∇ℓ∥A −B∥_ (84)

_where A_ _ℓ_ = (A Φ + AΦF [2][L][√]K).
_∇_ _∇_

_Proof. To begin with, we can apply the chain rule to obtain,_

_ℓ(Y, Φ(X;_ _, W))_ _ℓ(Y, Φ(X;_ _, W))_
_∥∇H_ _A_ _−∇H_ _B_ _∥_
= _ℓ(Y, Φ(X;_ _, W))_ **Φ(X;** _, W)_ _ℓ(Y, Φ(X;_ _, W))_ **Φ(X;** _, W)_ (85)
_∥∇_ _A_ _∇H_ _A_ _−∇_ _B_ _∇H_ _B_ _∥_

By applying the triangle inequality, and Cauchy-Schwartz it yields,

_ℓ(Y, Φ(X;_ _, W))_ _ℓ(Y, Φ(X;_ _, W))_
_∥∇H_ _A_ _−∇H_ _B_ _∥_
_ℓ(Y, Φ(X;_ _, W))_ **Φ(X;** _, W)_ **Φ(X;** _, W)_
_≤∥∇_ _A_ _∥∥∇H_ _A_ _−∇H_ _B_ _∥_
+ _ℓ(Y, Φ(X;_ _, W))_ _ℓ(Y, Φ(X;_ _, W))_ **Φ(X;** _, W)_ (86)
_∥∇H_ _A_ _−_ _B_ _∥∥∇H_ _B_ _∥_

We can now use Assumptions 1, 4, 5, and 6 as well as Proposition 3, to obtain

_ℓ(Y, Φ(X;_ _, W))_ _ℓ(Y, Φ(X;_ _, W))_ (A Φ + AΦF [2][L][√]K ) (87)
_∥∇H_ _A_ _−∇H_ _B_ _∥≤_ _∇_ _∥A −B∥_

Completing the proof.

**Lemma 5. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1 ≤_ _l ≤_ _L −_ 1. Let c ∈ (0, 1] and assume that the graphon
_convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN_
_H_
_sampled from Φ(X; H, W) as in (9). Under assumptions AS1–AS6, let the following condition be_
_satisfied for every k,_


_n log(2n[3][/][2])_


1 +


1 + _n log(2n_ )

_√KL[L][−][1]_ 6L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_   + [4][F][ L][−][1][L] + 12L[2]F [2][L][−][2]c

_δWW[c]_ _n_ p _n_ _n_

   

+ [2][F][ 2]√[L]n[√]K _<_ [1][ −] _[A]2[∇][ℓ][η][k]_ _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥._ (88)


_then the first k[∗]_ _iterates generated by equation (14), 1(k_ _k[∗])ℓ(Y, Φ(X;_ _k, Wn)) form a positive_
_≤_ _H_
_super-martingale with respect to the filtration Fk generated by the history of the Algorithm up to step_
_k [i.e., {X, Y, Xn, Yn, Wn}k, {X, Y, Xn, Yn, Wn}k−1, . . ., {X, Y, Xn, Yn, Wn}0]. Where k[∗]_ _is_
_the stopping time defined in Definition 8, and 1(·) is the indicator function._


-----

_Proof. To begin with, 1(k < k[∗])ℓ(Y, Φ(X; Hk, W)) ∈Fk, where Fk is the filtration generated by_
the history of the Algorithm up to k. Note that the loss function ℓ is positive by Assumption 4. It
remains to be shown the inequality expression of the super-martingale. For k > k[∗], the inequality
is trivially verified as the indicator function 1(k ≤ _k[∗]) = 0 for k > k[∗]. For k ≤_ _k[∗], as in Bertsekas_
& Tsitsiklis (2000), we define a continuous function g(ϵ) that takes the value of the loss function on
the Graphon data on iteration k + 1 at ϵ = 1, and on iteration k on ϵ = 0 as follows,

_g(ϵ) = ℓ(Y, Φ(X; Hk −_ _ϵηk∇Hℓ(Yn, Φ(Xn; Hk, Wn)), W))._ (89)

Note that function g(ϵ), is evaluated on the graphon data Y, X, W, but the steps are controlled by
the induced graphon data Yn, Xn, Wn. Applying the chain rule, the derivative of g(ϵ) with respect
to ϵ can be obtain as follows,

_∂_

(90)
_∂ϵ_ _[g][(][ϵ][) =]_

_−∇Hℓ(Y, Φ(X; Hk −_ _ϵηk∇Hℓ(Yn, Φ(Hk; Wn; Xn)), W))ηk∇Hℓ(Yn, Φ(Xn; Hk; Wn))._

Now note that the difference in the loss function ℓ between iterations k + 1 and k can be written as
the difference between g(ϵ = 1) and g(ϵ = 0) as follows,

_g(1)_ _g(0) = ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W))._ (91)
_−_ _H_ _−_ _H_

Computing the integration of the derivative of g(ϵ) between [0, 1] it yields

1

_∂_

_ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W)) = g(1)_ _g(0) =_ (92)
_H_ _−_ _H_ _−_ 0 _∂ϵ_ _[g][(][ϵ][)][dϵ]_
Z

1
= 0 _∇Hℓ(Y, Φ(X; Hk −_ _ϵηk∇Hℓ(Yn, Φ(Xn; Hk, Wn)), W))_
Z

(−)ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))dϵ. (93)


Note that the last term of the previous integral does not depend on ϵ. Besides, we can sum and
subtract ∇Hℓ(Y, Φ(Hk, W, X)) inside the integral, to obtain,

_ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W))_
_H_ _−_ _H_
= (−)ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))
1

0 _∇Hℓ(Y, Φ(X; Hk −_ _ϵηk∇Hℓ(Yn, Φ(Xn; Hk, Wn)), W))_

Z

+ ∇Hℓ(Y, Φ(X; Hk, W)) −∇Hℓ(Y, Φ(X; Hk, W))dϵ (94)

1
= −ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))∇Hℓ(Y, Φ(X; Hk, W)) 0 _dϵ_
Z

1
_−_ _ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))_ 0 _∇Hℓ(Y, Φ(Hk −_ _ϵηk∇Hℓ(Yn, Φ(Xn; Hk, Wn)), W, X))_
Z

_−∇Hℓ(Y, Φ(X; Hk, W))dϵ._ (95)


We can now apply the Cauchy-Schwartz inequality to the last term on the previous inequality, and
take the norm of the integral, which is smaller that the integral of the norm to obtain,

_ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W))_
_H_ _−_ _H_
_≤−ηk∇Hℓ(Yn, Φ(Hk; Wn; Xn))∇Hℓ(Y, Φ(Hk, W, X))_

1
+ ηk∥∇Hℓ(Yn, Φ(Hk; Wn; Xn))∥ 0
Z

_ℓ(Y, Φ(_ _k, W, X))_ _[∇][H][ℓ][(][Y,][ Φ][(][H][k][ −]_ _[ϵη][k][∇][H][ℓ][(][Y][n][,][ Φ][(][H][k][;][ W][n][;][ X][n][))][,][ W](96)[, X][))]_
_−∇H_ _H_ _[dϵ.]_


-----

Under Lemma 4, we can take the Lipschitz bound on the gradient on the loss function with respect
to the parameters, using A∇ℓ, to obtain,

_ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W))_
_H_ _−_ _H_
_≤−ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))∇Hℓ(Y, Φ(X; Hk, W))_

1
+ A∇ℓηk∥∇Hℓ(Yn, Φ(Xn; Hk, Wn))| 0 (97)
Z

_≤−ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))∇Hℓ(Y, Φ(X; H[η][k]k[∇], W[H][ℓ][(]))[Y][n][,][ Φ][(][X][n][;][ H][k][,][ W][n][))]_ _[ϵdϵ]_

_k[A][∇][ℓ]_
+ _[η][2]_ 2 _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥[2]._ (98)


Instead of evaluating the internal product between the gradient on the graphon, and induced graphon,
we will use Theorem 1, to bound their expected difference (cf. Figure 3 for intuition). We can add
and subtract the gradient of the loss function on the induced graphon ∇Hℓ(Yn, Φ(Hk; Wn; Xn)),
and use the Cauchy-Schwartz inequality to obtain,

_ℓ(Y, Φ(X;_ _k+1, W))_ _ℓ(Y, Φ(X;_ _k, W))_
_H_ _−_ _H_
_≤−ηk∇Hℓ(Yn, Φ(Xn; Hk, Wn))_
(∇Hℓ(Y, Φ(X; Hk, W)) + ∇Hℓ(Yn, Φ(Xn; Hk, Wn)) −∇Hℓ(Yn, Φ(Xn; Hk, Wn)))

_k[A][∇][Φ]_
+ _[η][2]_ 2 _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥[2]_ (99)

_≤−ηk∥∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥[2]_

+ ηk∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥∥∇Hℓ(Y, Φ(X; Hk, W)) −∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥

_k[A][∇][Φ]_
+ _[η][2]_ 2 _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥[2]._ (100)


We can rearrange the previous expression, to obtain,

_ηk∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥[2]_ 1 − _[A][∇]2[ℓ][η][k]_


_−_ _[∥∇][H][ℓ][(][Y,][ Φ][(][X]∥∇[;][ H]H[k]ℓ[,][ W](Yn[))], Φ[ −∇](Xn[H]; H[ℓ][(][Y]k[n], W[,][ Φ]n[(][X])∥[n][;][ H][k][,][ W][n][))][∥]_ 

_ℓ(Y, Φ(X;_ _k, W))_ _ℓ(Y, Φ(X;_ _k+1, W))._ (101)
_≤_ _H_ _−_ _H_


We can now take the conditional expectation with respect to the filtration Fn to obtain,

E[ℓ(Y, Φ(X; Hk+1, W))|Fk]

_≤_ _ηk∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥[2]_ 1 − _[A][∇]2[ℓ][η][k]_


E _∥∇Hℓ(Y, Φ(X; Hk, W)) −∇Hℓ(Yn, Φ(Xn; Hk, Wn))∥_
_−_  _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥_ 

+ ℓ(Y, Φ(X; _k, W))._ (102)
_H_ _[F][k]_

As step size ηk > 0, and by definition norms are non-negative, using Theorem 1, as condition (88)
holds for k ≤ _k[∗], then_

E[ℓ(Y, Φ(X; Hk+1, W))|Fk] ≤ _ℓ(Y, Φ(X; Hk, W))._ (103)

By definition of super-martingale as in Durrett (2019), we complete the proof.

**Lemma 6. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1 ≤_ _l ≤_ _L −_ 1. Let c ∈ (0, 1] and assume that the graphon
_convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN_
_H_


-----

_sampled from Φ(X; H, W) as in (9). Under assumptions AS1–AS6, for any ϵ ∈_ (0, 1 − _A∇ℓη), if_
_the iterates generated by (14), satisfy,_


_n log(2n[3][/][2])_


1 +


1 + _n log(2n_ )

_√KL[L][−][1]_ 6L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_   + [4][F][ L][−][1][L] + 12L[2]F [2][L][−][2]c

_δWW[c]_ _n_ p _n_ _n_

   

+ [2][F][ 2]√[L]n[√]K _<_ [1][ −] _[A][∇]2[ℓ][η][k][ −]_ _[ϵ]_ _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥._ (104)


_then the expected value of the stopping time k[∗]_ _[cf. Definition 8], is finite, i.e.,_

E[k[∗]] = O(1/ϵ) (105)

_Proof. Given the iterates at k = k[∗], and the initial values at k = 0, we can express the expected_
difference between the loss ℓ, as the summation over the difference of iterates as follows,


E[ℓ(Y, Φ(X; 0, W)) _ℓ(Y, Φ(X;_ _k∗_ _, W))] =_
_H_ _−_ _H_

_k∗_

E _ℓ(Y, Φ(X;_ _k_ 1, W)) _ℓ(Y, Φ(X;_ _k, W)_

"k=1 _H_ _−_ _−_ _H_
X


(106)


Taking the expected value with respect to the final iterate k = k[∗], we get,

E _ℓ(Y, Φ(X;_ _k0_ )) _ℓ(Y, Φ(X;_ _k∗_ ))
_H_ _−_ _H_
 

_k[∗]_
= Ek[∗] E _ℓ(Y, Φ(X;_ _k_ 1, W)) _ℓ(Y, Φ(X;_ _k, W)_ (107)

_H_ _−_ _−_ _H_

  _k=1_  
X

_∞_ _t_

_[k][∗]_

= E _ℓ(Y, Φ(X;_ _k_ 1, W)) _ℓ(Y, Φ(X;_ _k, W)_ _P_ (k[∗] = t)

_H_ _−_ _−_ _H_

_t=0_  _k=1_ 

X X

(108)


Using condition (104), and Lemma 5 for any k ≤ _k[∗], it verifies_

E _ℓ(Y, Φ(X;_ _k_ 1, W)) _ℓ(Y, Φ(X;_ _k, W))_ _η(√_
_H_ _−_ _−_ _H_ _≥_
 

Thus, coming back to (108),


_KF_ _[L][−][1]12L[2]F_ [2][L][−][2]c)[2]ϵ (109)


_KF_ _[L][−][1]12L[2]F_ [2][L][−][2]c)[2]ϵ


_tP_ (k[∗] = t)

_t=0_

X

(110)


E _ℓ(Y, Φ(X;_ _k0_ _, W))_ _ℓ(Y, Φ(X;_ _k∗_ _, W))_ _η(_
_H_ _−_ _H_ _≥_
 

_≥_ _η(_

Note that as the loss function ℓ is non-negative,


_KF_ _[L][−][1]12L[2]F_ [2][L][−][2]c)[2]ϵE[k[∗]] (111)


_ℓ(Y, Φ(X;_ _k0_ _, W))_
_H_

_KF_ _[L][−][1]12L[2]F_ [2][L][−][2]c)[2]ϵ _≥_ E[k[∗]] (112)


_η(_


Thus concluding that k[∗] = O(1/ϵ).

**Theorem 2. Consider the ERM problem in (12) and let Φ(X; H, W) be an L-layer WNN with**
_F0 = FL = 1, and Fl = F for 1 ≤_ _l ≤_ _L −_ 1. Let c ∈ (0, 1] and assume that the graphon
_convolutions in all layers of this WNN have K filter taps [cf. (6)]. Let Φ(xn;_ _, Sn) be a GNN_
_H_
_sampled from Φ(X; H, W) as in (9). Consider the iterates generated by equation (14). Under_


-----

_Assumptions AS1-AS6, for any fixed ϵ ∈_ (0, 1 − _A∇ℓη), if at each step k the number of nodes n is_
_picked such that it verifies_


1 +


_n log(2n[3][/][2])_


1 + _n log(2n_ )

_√KL[L][−][1]_ 6L[2]F [2][L][−][2] 1 + _[πB]W[c]_ _n_   + [4][F][ L][−][1][L] + 12L[2]F [2][L][−][2]c

_δWW[c]_ _n_ p _n_ _n_

   

+ [2][F][ 2]√[L]n[√]K _<_ [1][ −] _[A][∇]2[ℓ][η][k][ −]_ _[ϵ]_ _∥∇Hℓ(Yn, Φ(Xn; Hk, Wn)∥._ (113)


_then in finite time we will achieve an iterate k[∗]_ _such that the coefficients_ _k[∗]_ _satisfy_
_H_


E[∥∇Hℓ(Y, Φ(X; Hk∗ _, W))∥] ≤_ 24

_where A_ _ℓηk = (A_ Φ + AΦF [2][L][√]K).
_∇_ _∇_


_KF_ _[L][−][1]L[2]F_ [2][L][−][2]c _with probability 1_ (114)


_Proof. We can use Lemma 6, to conclude that it must be the case that P_ (k[∗] = ∞) = 0, which
implies that, P (k[∗] _< ∞) = 1. Using stopping time k[∗]_ condition [cf. Definition 8] and the triangle
inequality, it yields,

E[∥∇Hℓ(Y, Φ(X; Hk∗ _, W))∥] ≤∥∇Hℓ(Yn, Φ(Xn; Hk∗_ _, Wn))∥_ (115)
+E[ _ℓ(Yn, Φ(X;_ _k[∗]_ _, Wn))_ _ℓ(Y, Φ(Xn;_ _k[∗]_ _, Wn))_ ]
_∥∇H_ _H_ _−∇H_ _H_ _∥_

Note that the iterates are constructed such that, for every k

E[ _ℓ(Yn, Φ(X;_ _k, Wn))_ _ℓ(Y, Φ(Xn;_ _k, Wn))_ ] _ℓ(Yn, Φ(X;_ _k, Wn))_ _._
_∥∇H_ _H_ _−∇H_ _H_ _∥_ _≤∥∇H_ _H_ _∥_
(116)

Using the stopping time condition, the final result is attained as follows

E[∥∇Hℓ(Y, Φ(X; Hk∗ _, W))∥] ≤2∥∇Hℓ(Yn, Φ(X; Hk∗_ _, Wn))∥_ (117)

_≤24√KF_ _[L][−][1]L[2]F_ [2][L][−][2]c. (118)


-----

