Year: 2021

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.

heading: 1 INTRODUCTION
text: 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 & Peyré, 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 (Lovász, 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.

heading: 2 RELATED WORK
text: 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.

heading: 3 PRELIMINARY DEFINITIONS

heading: 3.1 GRAPH NEURAL NETWORKS
text: 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 ∈ Rn×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 ∈ Rn 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∈Ni [S]ijxj where Ni = {j | [S]ij 6= 0}. Information 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 contribution 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
y = h ∗S x = K−1∑ k=0 hkS kx (1)
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ΛVH . 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 −1 ≤ λ−1 ≤ λ−2 ≤ . . . ≤ 0 ≤ . . . ≤ λ2 ≤ λ1 ≤ 1. Since the eigenvectors of S form an orthonormal basis of Rn, we can project (1) onto this basis to obtain the spectral representation of the graph convolution, which is given by
h(λ) = K−1∑ k=0 hkλ k. (2)
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 xfl , 1 ≤ f ≤ Fl which we stack in a matrix Xl = [x1l , . . . ,x Fl l ] ∈ Rn×Fl . 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 ∈ RFl−1×Fl , we write the l-th layer of the GNN as
Xl = ρ ( K−1∑ k=0 SkXl−1Hlk ) . (3)
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 ∈ Rn×F0 . 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 6= n.

heading: 3.2 GRAPHON INFORMATION PROCESSING
text: 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 Lovász (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 (Lovász, 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.

heading: 3.2.1 GRAPHON NEURAL NETWORKS
text: Graphon data is defined as functions X ∈ L2([0, 1]). Analogously to graph data, graphon data can be diffused by application of a linear integral operator parametrized by W and defined as
TWX(v) = ∫ 1 0 W(u, v)X(u)du. (5)
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
Y = h ∗W X = K−1∑ k=0 hk(T (k) W X)(v) with
(T (k) W X)(v) = ∫ 1 0 W(u, v)(T (k−1) W X)(u)du
(6)
where T (0)W = 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 as −1 ≤ λ−1 ≤ λ−2 ≤ . . . ≤ 0 ≤ . . . λ2 ≤ λ1 ≤ 1. The eigenfunctions form an orthonormal basis of L2([0, 1]). Projecting the filter (1) onto this basis, we see that the graphon convolution admits a spectral representation given by
h(λ) = K−1∑ k=0 hkλ k. (7)
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} ∈ RFl−1×Fl . This allows writing the f th feature at layer l as
Xfl = ρ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf  (8) for 1 ≤ f ≤ Fl. For an L-layer WNN, (8) is repeated for 1 ≤ ` ≤ L. The WNN output is given by Y f = XfL, and X g 0 is given by the input data X
g for 1 ≤ g ≤ F0. We assume FL = F0 = 1 so that XL = Y and X0 = X . A more succinct representation of this WNN is the map Y = Φ(X;H,W), where the tensorH = {Hlk}l,k groups the filter coefficients at all layers.

heading: 3.2.2 SAMPLING GNNS FROM WNNS
text: 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)) [xn]i = X(ui)
(9)
i.e., the WNN can be seen as a generative model for GNNs Φ(xn;H,Sn). 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) = n∑ i=1 n∑ j=1 [Sn]ijI(u ∈ Ii)I(v ∈ Ij)
Xn(u) = n∑ i=1 [xn]iI(u ∈ Ii) (10)
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.

heading: 4 GRAPHON EMPIRICAL LEARNING
text: On graphons, the statistical loss minimization (or statistical learning) problem is given by
minimize H Ep(Y,X)[`(Y,Φ(X;H,W))] (11)
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.

heading: 4.1 GRAPHON EMPIRICAL RISK MINIMIZATION
text: Suppose that we have access to samples of the distribution D = {(Xj , 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
minimize H |D|∑ j=1 `(Y j ,Φ(Xj ;H,W)) (12)
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
Hk+1 = Hk − ηk∇H`(Y j ,Φ(Xj ;Hk,W)) (13)
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∇H`(Y j ,Φ(Xj ;H,W)) by sampling stochastic graphs Gn with GSO Sn [cf. (4)] and calculating ∇H`(yjn,Φ(xjn;H,Sn)), where xjn and yjn are as in (9). In this case, the graphon empirical learning step in (13) becomes
Hk+1 = Hk − ηk∇H`(yjn,Φ(xjn;Hk,Sn)) (14)
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.

heading: 4.2 GRADIENT APPROXIMATION
text: 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 ∑d i=1 |vi − ui|. If A = 1, we say that this function is normalized Lipschitz.
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: InitializeH0, 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: Take learning stepHk+1 = Hk − ηk∇`(yn,Φ(xn;Hk,Sn)) 7: end for 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(2n/ξ)/dW > 2/dW where dW denotes the maximum degree of the graphon W, i.e., dW = maxv ∫ 1 0 W(u, v)du.
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).
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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under assumptions AS1–AS5, it holds that
E[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖] ≤ γc+O (√ log(n3/2)
n
) (15)
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.

heading: 4.3 ALGORITHM CONSTRUCTION
text: 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.

heading: 4.4 ALGORITHM CONVERGENCE
text: 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;H,W) isAΦ-Lipschitz, and its gradient∇HΦ(X;H,W) is A∇Φ-Lipschitz, with respect to the parametersH [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;H,Sn) 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
γc+O (√ log(n3/2)
n
) <
1−A∇`η − 2 ‖∇H`(Yn,Φ(Xn;Hk,Wn)‖ (16)
then in finite time we will achieve an iterate k∗ such that the coefficientsHk∗ satisfy
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.

heading: 5 NUMERICAL RESULTS

heading: 5.1 RECOMMENDATION SYSTEM
text: 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.

heading: 5.2 DECENTRALIZED CONTROL
text: 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) ∈ R2, and speed vi(t) ∈ R2, 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)T 2s /2 + vi(t)Ts + ri(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 ∼ U [−3, 3], and then adding independent uniform noise U [−3, 3] to each agent. The initial deployment is randomly selected in a circle always verifying that the minimum distance between two agents is larger than 0.1m. On the one hand, agents pursue a common average velocity v̄ = (1/n) ∑n i=1 vi(t), thus minimizing the velocity variation of the team. On the other
hand, agents are required to avoid collision. We can thus define the velocity variation of the team σv(t) = ∑n i=1 ‖vi(t)− v̄(t)‖2, and the collision avoidance potential
CA(ri, rj) = { 1/‖ri(t)− rj(t)‖2 − log(‖ri(t)− rj(t)‖2) if ‖ri(t)− rj(t)(t)‖ ≤ RCA 1/R2CA − log(R2CA) otherwise, (18)
with RCA = 1m. A centralized controller can be obtain by ui(t)∗ = −n(vi − v̄) +∑n j=1∇riCA(ri, rj) (Tanner et al., 2003).
Exploiting 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 provided 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 −Φ(xm;H,S)‖2. Denoting Ni(t) the neighborhood of agent i at time t, the state of the agents x(t) = [x(t)1, . . . , x(t)n], xi(t) ∈ R6, is given by xi(t) = ∑ j:j∈Ni(t)[vi(t) − vj(t), rij(t)/‖rij(t)‖
4, rij(t)/‖rij(t)‖2]. Note that state xi(t), gets transmitted between agents if a communication link exists between them.
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.

heading: 6 CONCLUSIONS
text: 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.

heading: A APPENDIX

heading: B NUMERICAL RESULTS PARAMETERS
text: 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/graphneural-networks/blob/master/examples/movieGNN.py and implemented with PyTorch.
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/graphneural-networks/blob/master/examples/flockingGNN.py and implemented with PyTorch.

heading: C PROOF OF THEOREM 1
text: Definition 2 (Template graphs). Let {ui}ni=1 be the regular n-partition of [0, 1], i.e.,
ui = i− 1 n
(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., 1 ≥ λ1 ≥ λ2 ≥ · · · ≥ 0 ≥ · · · ≥ λ−2 ≥ λ−1 ≥ −1, and their only accumulation point is 0 (Lax, 2002, Theorem 3, Chapter 28). Thus, we define the spectral representation of the convolutional filter TH (cf. (6)) as,
h(λ) = K−1∑ k=0 hkλ k (21)
Definition 4 (c-band cardinality of W). The c-band cardinality, denoted BcW, is the number of eigenvalues whose absolute value is larger than c.
BcW = #{λ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,
δcWWn = mini,j 6=i {‖λi(TW)− λi(TWn)‖ : ‖λi(TWn)‖ ≥ 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,
(THX)(v) = K−1∑ k=0 hk(T (k) W X)(v). (24)
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 ≤ 1
n . (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
‖X(u)−Xn(u)‖ ≤ max (∣∣∣∣u− i− 1n ∣∣∣∣ , ∣∣∣∣ in − u ∣∣∣∣) ≤ 1n. (26)
We can then write
‖X −Xn‖2 = ∫ 1
0
|X(u)−Xn(u)|2du (27)
≤ ∫ 1
0
( 1
n
)2 du = ( 1
n
)2 , (28)
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‖ ≤ 2
n . (29)
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
‖W(u, v)−Wn(u, v)‖ ≤ max (∣∣∣∣u− i− 1n ∣∣∣∣ , ∣∣∣∣ in − u ∣∣∣∣) (30)
+ max (∣∣∣∣v − j − 1n ∣∣∣∣ , ∣∣∣∣ jn − v ∣∣∣∣) (31) ≤ 1 n + 1 n = 2 n . (32)
We can then write
‖W −Wn‖2 = ∫ 1
0
|W(u, v)−Wn(u, v)|2dudv (33)
≤ ∫ 1
0
( 2
n
)2 dudv = ( 2
n
)2 (34)
which concludes the proof.
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 parametersH = {Hlk}l,k can be upper bounded by,
‖∇HΦ(X;H,W)‖ ≤ F 2L √ K. (35)
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),
‖∇[H l†k† ]g†f† Φ(X;H,W)‖ = ∥∥∥∥∇[Hl†k† ]g†f†XfL ∥∥∥∥ (36) =
∥∥∥∥∇[Hl†k† ]g†f†ρ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[HLk]gf ∥∥∥∥. (37) By Assumption 3, the non-linearity ρ is normalized Lipschitz , i.e. ∇ρ(·)(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,
‖∇[H l†k† ]g†f† Φ(X;H,W)‖ = ∥∥∥∥∇ρ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[HLk]gf ∥∥∥∥ ∥∥∥∥∇[Hl†k† ]g†f† Fl−1∑
g=1 K−1∑ k=1 (T (k) W X g l−1)[HLk]gf ∥∥∥∥ (38) ≤ ∥∥∥∥∇[Hl†k† ]g†f† Fl−1∑
g=1 K−1∑ k=1 (T (k) W X g l−1)[HLk]gf ∥∥∥∥ (39) Note that the a larger bound will occur if l† < L− 1, then by linearity of derivation, and the triangle inequality we obtain,
‖∇[H l†k† ]g†f† Φ(X;H,W)‖ ≤ Fl−1∑ g=1 ∥∥∥∥K−1∑ k=1 T (k) W (∇[Hl†k† ]g†f†X g l−1)[HLk]gf ∥∥∥∥ (40)
By Assumption 2, the convolutional filters are non-amplifying, thus it holds that,
‖∇[H l†k† ]g†f† Φ(X;H,W)‖ ≤ Fl−1∑ g=1 ∥∥∥∥∇[Hl†k† ]g†f†Xgl−1 ∥∥∥∥ (41)
Now note that as filters are non-amplifying, the maximum difference in the gradient will be attained 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,
‖∇[H l†k† ]g†f† Φ(X;H,W)‖ ≤ FL−1 ∥∥∥∥hi∗WX0∥∥∥∥ (42)
≤ FL−1‖X0‖. (43)
To complete the proof note that tensor H has FL−1K 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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under assumptions (1),(2),(3), and (5) with probability 1− ξ it holds that,
‖Φ(X;H,W)−Φ(Xn;H,Wn)‖ ≤LFL−1 (
1 + πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 1
n + 4LFL−1c (44)
The fixed constants BcW and δ c WWn 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;H,W)−Φ(Xn;H,Wn)‖ = ‖XL −XnL‖ (45)
= ∥∥∥∥∥∥ρ FL−1∑
g=1 K−1∑ k=1 (T (k) W X g L−1)[HLk]gf
− ρ FL−1∑
g=1 K−1∑ k=1 (T (k) Wn XgnL−1)[HLk]gf ∥∥∥∥∥∥ . Since the non-linearity ρ is normalized Lipschitz by Assumption 3, using the triangle inequality, we obtain
‖XL −XnL‖ ≤ FL−1∑ g=1 ∥∥∥∥∥ K−1∑ k=1 (T (k) W X g L−1)[HLk]gf − K−1∑ k=1 (T (k) Wn XgnL−1)[HLk]gf ∥∥∥∥∥ . (46) Using the triangle inequality once again, we split the last inequality into two terms as follows,
‖XL −XnL‖ ≤ FL−1∑ g=1 ∥∥∥∥∥ K−1∑ k=1 T (k) W (X g L−1 −X g nL−1)[HLk]gf ∥∥∥∥∥ (1) +
FL−1∑ g=1 ∥∥∥∥∥ K−1∑ k=1 (T (k) W − T (k) Wn )XgL−1[HLk]gf ∥∥∥∥∥ (2). (47) 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 ∥∥∥∥∥ K−1∑ k=1 T (k) W (X g L−1 −X g nL−1)[HLk]gf ∥∥∥∥∥ ≤ FL−1∑ g=1 ∥∥XgL−1 −XgnL−1∥∥ . (48)
To bound term (2), denoting hLgf the spectral representation of the convolutional filter applied to XgL−1 at feature f of layer L [cf. Definition 3], we will decompose the filter as follows,
h≥cLgf (λ) { 0 if |λ| < c hLgf (λ)− hLgf (c) if |λ| ≥ c
(49)
h<cLgf (λ) { hLgf (λ) if |λ| < c hLgf (c) if |λ| ≥ c.
(50)
Note that hLgf = h ≥c Lgf + h <c Lgf . Let T <c [HL]gf and T<c[HnL]gf , be the graphon convoutional filters with filter function h<cLgf on graphons W, and Wn respectively [cf. Definition 6]. Note that filter h<cLgf , varies only in the interval [0, c), and since filters are normalized Lipschitz by Assumption 2, it verifies∥∥∥T<c[HL]gfXgL−1 − T<c[HnL]gfXgL−1∥∥∥ ≤ ∥∥(hLgf (c) + c)− (hLgf (c)− c)‖‖XgL−1∥∥ (51)
≤ 2c‖XgL−1‖. (52)
Now we need to upper bound the difference in the high frequencies h≥cLgf . Let T ≥c [HL]gf and T≥c[HnL]gf , be the graphon filters with filter function h≥cLgf on graphons W, and Wn respectively. Let Sn 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≥c[HnL]gf , the graph filter with filter function h≥cLgf on graphon Wn, we can use the triangle inequality to obtain,∥∥∥T≥c[HL]gfXgL−1 − T≥c[HnL]gfXgL−1∥∥∥ ≤ ∥∥∥T≥c[HL]gfXgL−1 − T≥c[HnL]gfXgL−1∥∥∥ (2.1)
+ ∥∥∥T≥c[HnL]gfXgL−1 − T≥c[HnL]gfXgL−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,
∥∥∥T≥c[HL]gfXgL−1 − T≥c[HnL]gfXgL−1∥∥∥ ≤ (1 + πBcWnδcWWn )2(1 +√n log( 2nξ )) n
∥∥XgL−1∥∥ . (54) Where the fixed constants BcW and δ c WWn
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,
‖XL −XnL‖ ≤ FL−1∑ g=1 ∥∥XgL−1 −XgnL−1∥∥+ 2c‖XgL−1‖ + ( 1 +
πBcWn δcWWn )2(1 +√n log( 2nξ )) n
∥∥XgL−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‖+ 2LFL−1c‖X0‖
+ LFL−1 (
1 + πBcWn δcWWn )2(1 +√n log( 2nξ )) 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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under assumptions (1),(2),(3), and (5) with probability 1− ξ it holds that,
‖∇HΦ(X;H,W)−∇HΦ(Xn;H,Wn)‖
≤ √ KFL−1 ( 2L2F 2L−2 ( 1 +
πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 2FL−1L
n + 8L2F 2L−2c
) .
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,
‖∇[H l†k† ]g†f† Φ(X;H,W)−∇[H l†k† ]g†f† Φ(Xn;H,Wn)‖
= ‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ (57)
= ∥∥∥∥∇[Hl†k† ]g†f†ρ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf  −∇[H
l†k† ]g†f† ρ Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ∥∥∥∥. (58) Taking derivatives by applying the chain rule, and applying the triangle inequality it yields,
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖
≤ ∥∥∥∥(∇ρ( Fl−1∑
g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf ) −∇ρ ( Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ))
∇[H l†k† ]g†f† ( Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf )∥∥∥∥ (59) +
∥∥∥∥∇ρ( Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ) (60)
( ∇[H
l†k† ]g†f† Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf −∇[Hl†k† ]g†f† Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf )∥∥∥∥. 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,
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ (61)
≤ ∥∥∥∥ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf − Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ∥∥∥∥FL−1‖X0‖ +
∥∥∥∥ Fl−1∑ g=1 ∇[H l†k† ]g†f† K−1∑ k=1 ( (T (k) W X g l−1)[Hlk]gf − (T (k) Wn Xgnl−1)[Hlk]gf ) )∥∥∥∥.
We can now apply the triangle inequality on the second term of the previous bound to obtain,
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ (62)
≤ ∥∥∥∥ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf − Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ∥∥∥∥FL−1‖X0‖ +
∥∥∥∥ Fl−1∑ g=1 ∇[H l†k† ]g†f† K−1∑ k=1 ( (T (k) W )[Hlk]gf − (T (k) Wn )[Hlk]gf ) ) Xgnl−1 ∥∥∥∥ +
Fl−1∑ g=1 ∥∥∥∥∇[Hl†k† ]g†f† K−1∑ k=1 T (k) Wn ( Xgl−1 −X g nl−1 ) [Hlk]gf ) ∥∥∥∥. Now note that as we are considering the case in which l† < l−1, using Cauchy-Schwartz inequality, 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
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ (63)
≤ 2 ∥∥∥∥ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf − Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ∥∥∥∥FL−1‖X0‖ +
Fl−1∑ g=1 ∥∥∥∥∇[Hl†k† ]g†f†(Xgl−1 −Xgnl−1) ∥∥∥∥.
Now notice, that the only term that remains to bound is the exact same bound we obtained in equation (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 tensorH.
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ (64)
≤ 2LFL−2 ∥∥∥∥ Fl−1∑ g=1 K−1∑ k=1 (T (k) W X g l−1)[Hlk]gf − Fl−1∑ g=1 K−1∑ k=1 (T (k) Wn Xgnl−1)[Hlk]gf ∥∥∥∥FL−1‖X0‖ +
Fl−1∑ g=1 ∥∥∥∥∇[Hl†k† ]g†f†(Xg1 −Xg1) ∥∥∥∥.
Note that the derivative of a convolutional filter 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,∥∥∥∥hi∗WnXn0 − hi∗WX0∥∥∥∥ ≤ (‖W −Wn‖+ ‖Wn −Wn‖)‖X0‖+ ‖Xn0 −X0‖ (65)
≤ (
1 + πBcWn δcWWn )2(1 +√n log( 2nξ )) n + 1 n (66)
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,
‖∇[H l†k† ]g†f† XfL −∇[Hl†k† ]g†f†X f nL‖ ≤2L
2F 2L−2 (
1 + πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 2FL−1L
n + 8L2F 2L−2c. (67)
To achieve the final result, note that tensorH has KFL−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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under Assumptions (1)–(5) with probability 1− ξ it holds that,
‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖
≤ √ KFL−1 ( 3L2F 2L−2 ( 1 +
πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 4FL−1L
n + 12L2F 2L−2c ) 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,
‖∇[H l†k† ]g†f† `(Y,Φ(X;H,W))−∇[H l†k† ]g†f† `(Yn,Φ(Xn;H,Wn))‖ =‖∇`(Y,Φ(X;H,W))∇[H
l†k† ]g†f† Φ(X;H,W)
−∇`(Yn,Φ(Xn;H,Wn))∇[H l†k† ]g†f† Φ(Xn;H,Wn)‖. (68)
By Cauchy-Schwartz, and the triangle inequality it holds,
‖∇[H l†k† ]g†f† `(Y,Φ(X;H,W))−∇[H l†k† ]g†f† `(Yn,Φ(Xn;H,Wn))‖ ≤‖∇`(Y,Φ(X;H,W))−∇`(Yn,Φ(Xn;H,Wn))‖‖∇[H
l†k† ]g†f† Φ(X;H,W)‖ (69)
+ ‖∇`(Yn,Φ(Xn;H,Wn))‖‖∇[H l†k† ]g†f† Φ(X;H,W)−∇[H l†k† ]g†f† Φ(Xn;H,Wn)‖.
By the triangle inequality and Assumption 4 it follows,
‖∇[H l†k† ]g†f† `(Y,Φ(X;H,W))−∇[H l†k† ]g†f† `(Yn,Φ(Xn;H,Wn))‖ ≤ ‖∇`(Y,Φ(X;H,W))−∇`(Y,Φ(Xn;H,Wn))‖‖∇[H
l†k† ]g†f† Φ(X;H,W)‖ (70)
‖∇`(Yn,Φ(Xn;H,Wn))−∇`(Y,Φ(Xn;H,Wn))‖‖∇[H l†k† ]g†f† Φ(X;H,W)‖ (71) + ‖∇[H
l†k† ]g†f† Φ(X;H,W)−∇[H l†k† ]g†f† Φ(Xn;H,Wn)‖
≤ (‖Yn − Y ‖+ ‖Φ(Xn;H,Wn))−Φ(X;H,W))‖)‖∇[H l†k† ]g†f† Φ(X;H,W)‖ (72) + ‖∇[H
l†k† ]g†f† Φ(X;H,W)−∇[H l†k† ]g†f† Φ(Xn;H,Wn)‖.
Now we can use Lemmas 1–2, Propositions 1, and 3, and Assumption 1 to obtain,
‖∇[H l†k† ]g†f† `(Y,Φ(X;H,W))−∇[H l†k† ]g†f† `(Yn,Φ(Xn;H,Wn))‖ (73)
≤ ( 3L2F 2L−2 (
1 + πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 4FL−1L
n + 12L2F 2L−2c ) Noting that tensor H has KFL−1 elements, and each individual term can be bounded by (73), the desired result is attained.
Definition 7. We define the constant γ as,
γ = 12 √ KFL−1L2F 2L−2, (74)
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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under assumptions AS1–AS5, it holds that
E[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖]
≤ √ KFL−1 ( 6L2F 2L−2 ( 1 +
πBcWn δcWWn )(1 +√n log(2n3/2)) n
+ 4FL−1L
n + 12L2F 2L−2c
) + 2F 2L √ K√
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 BcW and δ c WWn
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 = ( ‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖ (76)
≤ √ KFL−1 ( 3L2F 2L−2 ( 1 +
πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 4FL−1L
n + 12L2F 2L−2c
)) .
Thus, by considering the disjoint events An, and Acn, and denoting 1(·) the indicator function, the expectation can be separated as follows,
E[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖] = E[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖1(An)]
+ E[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖1(Acn)] (77) We can bound the term regardingAcn 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;H,W))‖‖∇HΦ(X;H,W)‖
+ ‖∇`(Yn,Φ(Xn;H,Wn))‖‖∇HΦ(Xn;H,Wn)‖ (79) ≤ ‖∇HΦ(X;H,W)‖+ ‖∇HΦ(Xn;H,Wn)‖ (80) ≤ 2F 2L √ 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[‖∇H`(Y,Φ(X;H,W))−∇H`(Yn,Φ(Xn;H,Wn))‖]
≤ (1− ξ) √ KFL−1 ( 3L2F 2L−2 ( 1 +
πBcWn δcWWn )2(1 +√n log( 2nξ )) n
+ 4FL−1L
n + 12L2F 2L−2c
) + ξ2F 2L √ K (82)
To complete the proof, set ξ = 1√ n .

heading: D PROOF OF THEOREM 2
text: Definition 8 (Stopping time). We define the stopping time k∗ as,
k∗ = min k≥0 {‖∇HΦ(X;Hk,Wn)‖ ≤
√ KFL−112L2F 2L−2c}. (83)
Definition 9 (Constant ψ). Lemma 4. Under Assumptions 4, 5, and 6, the gradient of the loss function ` with respect to the parameters of the GNNH is A∇`-Lipschitz,
‖∇H`(Y,Φ(X;A,W))−∇H`(Y,Φ(X;B,W))‖ ≤ A∇`‖A − B‖ (84)
where A∇` = (A∇Φ +AΦF 2L √ K).
Proof. To begin with, we can apply the chain rule to obtain,
‖∇H`(Y,Φ(X;A,W))−∇H`(Y,Φ(X;B,W))‖ = ‖∇`(Y,Φ(X;A,W))∇HΦ(X;A,W)−∇`(Y,Φ(X;B,W))∇HΦ(X;B,W)‖ (85)
By applying the triangle inequality, and Cauchy-Schwartz it yields,
‖∇H`(Y,Φ(X;A,W))−∇H`(Y,Φ(X;B,W))‖ ≤ ‖∇`(Y,Φ(X;A,W))‖‖∇HΦ(X;A,W)−∇HΦ(X;B,W)‖ + ‖∇H`(Y,Φ(X;A,W))− `(Y,Φ(X;B,W))‖‖∇HΦ(X;B,W)‖ (86)
We can now use Assumptions 1, 4, 5, and 6 as well as Proposition 3, to obtain
‖∇H`(Y,Φ(X;A,W))−∇H`(Y,Φ(X;B,W))‖ ≤ (A∇Φ +AΦF 2L √ K‖A − B‖) (87)
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;H,Sn) be a GNN sampled from Φ(X;H,W) as in (9). Under assumptions AS1–AS6, let the following condition be satisfied for every k,
√ KLL−1 ( 6L2F 2L−2 ( 1 +
πBcWn δcWWn )(1 +√n log(2n3/2)) n + 4FL−1L n + 12L2F 2L−2c
) + 2F 2L √ K√
n < 1−A∇`ηk 2 ‖∇H`(Yn,Φ(Xn;Hk,Wn)‖. (88)
then the first k∗ iterates generated by equation (14), 1(k ≤ k∗)`(Y,Φ(X;Hk,Wn)) form a positive super-martingale with respect to the filtrationFk 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,
∂
∂ g( ) = (90) −∇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;Hk+1,W))− `(Y,Φ(X;Hk,W)). (91)
Computing the integration of the derivative of g( ) between [0, 1] it yields
`(Y,Φ(X;Hk+1,W))− `(Y,Φ(X;Hk,W)) = g(1)− g(0) = ∫ 1
0
∂
∂ g( )d (92)
= ∫ 1 0 ∇H`(Y,Φ(X;Hk − ηk∇H`(Yn,Φ(Xn;Hk,Wn)),W))
(−)η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;Hk+1,W))− `(Y,Φ(X;Hk,W)) = (−)ηk∇H`(Yn,Φ(Xn;Hk,Wn))∫ 1 0 ∇H`(Y,Φ(X;Hk − ηk∇H`(Yn,Φ(Xn;Hk,Wn)),W))
+∇H`(Y,Φ(X;Hk,W))−∇H`(Y,Φ(X;Hk,W))d (94) = −ηk∇H`(Yn,Φ(Xn;Hk,Wn))∇H`(Y,Φ(X;Hk,W)) ∫ 1
0
d
− ηk∇H`(Yn,Φ(Xn;Hk,Wn)) ∫ 1
0
∇H`(Y,Φ(Hk − ηk∇H`(Yn,Φ(Xn;Hk,Wn)),W, X))
−∇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;Hk+1,W))− `(Y,Φ(X;Hk,W)) ≤ −ηk∇H`(Yn,Φ(Hk; Wn;Xn))∇H`(Y,Φ(Hk,W, X))
+ ηk‖∇H`(Yn,Φ(Hk; Wn;Xn))‖ ∫ 1
0 ∥∥∥∥∇H`(Y,Φ(Hk − ηk∇H`(Yn,Φ(Hk; Wn;Xn)),W, X)) −∇H`(Y,Φ(Hk,W, X)) ∥∥∥∥d . (96)
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;Hk+1,W))− `(Y,Φ(X;Hk,W)) ≤ −ηk∇H`(Yn,Φ(Xn;Hk,Wn))∇H`(Y,Φ(X;Hk,W))
+A∇`ηk‖∇H`(Yn,Φ(Xn;Hk,Wn))| ∫ 1
0 ∥∥∥∥ηk∇H`(Yn,Φ(Xn;Hk,Wn))∥∥∥∥ d (97) ≤ −ηk∇H`(Yn,Φ(Xn;Hk,Wn))∇H`(Y,Φ(X;Hk,W))
+ η2kA∇`
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;Hk+1,W))− `(Y,Φ(X;Hk,W)) ≤ −ηk∇H`(Yn,Φ(Xn;Hk,Wn)) (∇H`(Y,Φ(X;Hk,W)) +∇H`(Yn,Φ(Xn;Hk,Wn))−∇H`(Yn,Φ(Xn;Hk,Wn)))
+ η2kA∇Φ
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))‖
+ η2kA∇Φ
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∇`ηk 2
− ‖∇H`(Y,Φ(X;Hk,W))−∇H`(Yn,Φ(Xn;Hk,Wn))‖ ‖∇H`(Yn,Φ(Xn;Hk,Wn)‖ ) ≤ `(Y,Φ(X;Hk,W))− `(Y,Φ(X;Hk+1,W)). (101)
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∇`ηk 2
− E [ ‖∇H`(Y,Φ(X;Hk,W))−∇H`(Yn,Φ(Xn;Hk,Wn))‖
‖∇H`(Yn,Φ(Xn;Hk,Wn)‖ ∣∣∣∣Fk]) + `(Y,Φ(X;Hk,W)). (102)
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;H,Sn) be a GNN
sampled from Φ(X;H,W) as in (9). Under assumptions AS1–AS6, for any ∈ (0, 1 − A∇`η), if the iterates generated by (14), satisfy,
√ KLL−1 ( 6L2F 2L−2 ( 1 +
πBcWn δcWWn )(1 +√n log(2n3/2)) n + 4FL−1L n + 12L2F 2L−2c
) + 2F 2L √ K√
n < 1−A∇`ηk − 2 ‖∇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;H0,W))− `(Y,Φ(X;Hk∗ ,W))] =
E [ k∗∑ k=1 `(Y,Φ(X;Hk−1,W))− `(Y,Φ(X;Hk,W) ] (106)
Taking the expected value with respect to the final iterate k = k∗, we get, E [ `(Y,Φ(X;Hk0))− `(Y,Φ(X;Hk∗)) ] = Ek∗
[ E [ k∗∑ k=1 `(Y,Φ(X;Hk−1,W))− `(Y,Φ(X;Hk,W) ]∣∣∣∣k∗] (107)
= ∞∑ t=0 E [ t∑ k=1 `(Y,Φ(X;Hk−1,W))− `(Y,Φ(X;Hk,W) ] P (k∗ = t)
(108)
Using condition (104), and Lemma 5 for any k ≤ k∗, it verifies E [ `(Y,Φ(X;Hk−1,W))− `(Y,Φ(X;Hk,W)) ] ≥ η( √ KFL−112L2F 2L−2c)2 (109)
Thus, coming back to (108), E [ `(Y,Φ(X;Hk0 ,W))− `(Y,Φ(X;Hk∗ ,W)) ] ≥ η( √ KFL−112L2F 2L−2c)2 ∞∑ t=0 tP (k∗ = t)
(110)
≥ η( √ KFL−112L2F 2L−2c)2 E[k∗] (111)
Note that as the loss function ` is non-negative, E [ `(Y,Φ(X;Hk0 ,W)) ] η( √ KFL−112L2F 2L−2c)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;H,Sn) be a GNN 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
√ KLL−1 ( 6L2F 2L−2 ( 1 +
πBcWn δcWWn )(1 +√n log(2n3/2)) n + 4FL−1L n + 12L2F 2L−2c
) + 2F 2L √ K√
n < 1−A∇`ηk − 2 ‖∇H`(Yn,Φ(Xn;Hk,Wn)‖. (113)
then in finite time we will achieve an iterate k∗ such that the coefficientsHk∗ satisfy
E[‖∇H`(Y,Φ(X;Hk∗ ,W))‖] ≤ 24 √ KFL−1L2F 2L−2c with probability 1 (114)
where A∇`ηk = (A∇Φ +AΦF 2L √ K).
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[‖∇H`(Yn,Φ(X;Hk∗ ,Wn))−∇H`(Y,Φ(Xn;Hk∗ ,Wn))‖]
Note that the iterates are constructed such that, for every k
E[‖∇H`(Yn,Φ(X;Hk,Wn))−∇H`(Y,Φ(Xn;Hk,Wn))‖] ≤ ‖∇H`(Yn,Φ(X;Hk,Wn))‖. (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 √ KFL−1L2F 2L−2c. (118)

Review summaries:
This paper proposed an effective and scalable algorithm to train the graph neural network. Their method leverages the framework of graphon neural networks to enlarge the training size of GNN during the training. Specifically, a growing size subgraph is sampled from graphon by a Bernoulli distribution and fed to a graph neural network. A Theoretical guarantee is given for the convergence of the algorithm by proving the absolute value of gradient will decrease to a small value surely.

In this paper, the authors propose a method that learns a large graph neural network (GNN) starting from a relatively smaller GNN which would increase its number of nodes step by step (epoch by epoch). Most importantly, the paper gives mathematical proof to show that the gradient descent steps on GNN are close to the graphon neural network (WNN) which is the limit object of the GNN. This paper also demonstrates its proposed method on a recommendation system and a decentralized control problem and shows reduced computational cost.

This paper shows that a graph neural network (GNN) trained on an increasingly larger graph behaves like a graphon neural network (WNN), the "limit object of a GNN", which arises when considering a continuous graphon instead of a discrete graph (the graphon can be seen as a generative model from which the discrete graph is sampled).   This observation is translated into two theorems that bound the expected distance between the GNN and the WNN as a function of the number of nodes used to train the GNN.   The authors then propose an algorithm to train a GNN by gradually adding nodes to the training data, leading to a computational advantage compared to training the GNN on the full dataset from the start.