Year: 2020

Abstract:
We present Iterated Graph Neural Network System (IGNNS), a new framework of Graph Neural Networks (GNNs), which can deal with undirected graph and directed graph in a unified way. The core component of IGNNS is the Iterated Function System (IFS), which is an important research field in fractal geometry. The key idea of IGNNS is to use a pair of affine transformations to characterize the process of message passing between graph nodes and assign an adjoint probability vector to them to form an IFS layer with probability. After embedding in the latent space, the node features are sent to IFS layer for iterating, and then obtain the high-level representation of graph nodes. We also analyze the geometric properties of IGNNS from the perspective of dynamical system. We prove that if the IFS induced by IGNNS is contractive, then the fractal representation of graph nodes converges to the fractal set of IFS in Hausdorff distance and the ergodic representation of that converges to a constant matrix in Frobenius norm. We have carried out a series of semi-supervised node classification experiments on citation network datasets such as citeser, Cora and PubMed. The experimental results show that the performance of our method is obviously better than the related methods.

heading: 1 INTRODUCTION
text: GNN (Scarselli et al., 2009) has been proved to be effective in processing graph structured data, and has been widely used in natural language processing, computer vision, data mining, social network and biochemistry. In recent years, GNN has developed a variety of architectures, such as GCN (Kipf & Welling, 2017), GraphSAGE (Hamilton et al., 2017), GAT (Veličković et al., 2018), DGI (Veličković et al., 2019), GIN (Xu et al., 2019), GCNII (Ming Chen et al., 2020) and GEN (Li et al., 2020). These architectures have a common feature, that is, the representation of each node is updated using messages from its neighbors but without distinguishing the direction (or angle) of message passing between two nodes. Recent studies have shown that considering directed message passing between nodes can improve the performance of GNN and achieve success in related fields. For example, DimeNet (Klicpera et al., 2020) considers the spatial direction from one atom to another and can learn both molecular properties and atomic forces. R-GCN (Schlichtkrull et al., 2018) and Bi-GCN (Marcheggiani & Titov, 2017; Fu et al., 2019) are models for directed graph, applied in the field of natural language processing. We note that the above direction based model does not consider the bidirectional mixed passing of messages.
But in real life, message passing is interactive in different directions. For example, node A obtains a message from node B. After processing the message, node A not only passes it to the next node C, but also feeds back to node B. Suppose there are only two directions for message passing, forward and backward, represented by 0 or 1, respectively. The symbol space of the first generation message passing path is {0, 1} = {0, 1}1, and that of the second generation message passing path is {00, 01, 10, 11} = {0, 1}2. Generally, the symbol space of the n-th generation message passing path is {0, 1}n and the size of the symbol space is 2n. This means that the scope of message passing spreads with exponent 2. However, in Bi-GCN (similar to Bi-LSTM) and R-GCN architectures, the symbol space is {{0}n, {1}n}, and its size is 2, which indicates that a lot of information will be lost in the process of message passing (see Appendix A).
How to characterize the above message passing patterns? We use two mappings to represent message passing process in two directions. Then the interactive passing of messages in different directions is equivalent to the composite operation of corresponding mappings. In addition, the direction of
message passing is often random, so we endow the two mappings with an adjoint probability vector to reflect the randomness. Because the symbol space of the iterative path of the Iterated Function System (IFS) with two mappings is also {0, 1}n and the mapping is selected with a certain probability, the iterative process of IFS is similar to the message passing process. In other words, the above message passing pattern can be described perfectly by IFS with probabilities. We naturally present the Iterative Graph Neural Network System (IGNNS), whose core layer is constructed by IFS. Figure 1 describes the differences in message passing patterns among GCN, Bi-GCN and IGNNS. At the same time, we regard undirected graph as a directed graph with equal probability of bidirectional message passing (see Figure 1(a)), so the IGNNS architecture can handle directed graph and undirected graph in a unified way.

heading: 2 PRELIMINARIES
text: A graph G = (V,E) is defined by its note set V = {v1, v2, ..., vN} and edge set E = {(vi, vj)|vi, vj ∈ V }. Let A ∈ RN denote the adjacency matrix of G, providing with relational information between nodes. A[i, j] denote i, jth element of A, A[i, :] means the ith row, and A[:, j] means the jth column. In this paper, we assume that all nodes of G are self adjacent, that is A[i, i] = 1, i = 1, 2, ..., N . let D = diag(d1, d2, ..., dN ) be the degree matrix of A, where di = ∑N j=1 A[i, j].
Neighborhood Normalization. There are two ways to normalize A. One approach is the following mean-pooling employed by Hamilton et al. (2017) and Veličković et al. (2019) for inductive learning:
Amp = D −1A.
Another approach is the following symmetric normalization employed by Kipf & Welling (2017):
Asym = D − 12AD− 1 2 .
Iterated Function System. A mapping f : RN → RN is said to be a contractive mapping on RN if there exists a constant 0 < c < 1 such that ‖f(x1)− f(x2)‖2 < c‖x1 − x2‖2 for all x1, x2 ∈ RN . An iterated function system (Hutchinson, 1981) is defined by
IFS = {RN ; f1, f2, ..., fm;p}, where each fi : RN → RN is a contractive mapping and p = (p1, p2, ..., pm) is an adjoint probability vector meaning that fi is selected by probability pi for each iteration. Hutchinson (1981) showed that there exists a unique nonempty compact set F such that
F = m⋃ i=1 fi(F).
A =
 1 1 0 01 1 0 10 1 1 1
0 1 1 1
 , A0 =

1 2 1 2 0 0 0 12 0 1 2 0 0 12 1 2
0 0 0 1
 and A1 =

1 0 0 0 1 2 1 2 0 0
0 12 1 2 0 0 13 1 3 1 3
 .
We call F the fractal set or invariant set of IFS. More conclusions on IFS can be found in the Appendix D. It is well known that there exists a unique probability measure µ with support F satisfying the equation
µ = m∑ i=1 piµ ◦ f−1i . (1)
The probability measure µ in (1) is called the self-similar measure of IFS with probability vector p.

heading: 3 IGNNS ARCHITECTURE
text: In this section, we will introduce the architecture of the IGNNS according to the input layer, IFS layer, representation layer and output layer, which is described in Figure 2.

heading: 3.1 INPUT LAYER
text: Given a graph structure data X ∈ RN×F of G = (V,E), called as the feature matrix of node set V . A row of X represents the F -dimensional feature vector of a node in V . Let W int ∈ RF×H be a learnable parameter matrix, where H is the dimension of the latent space. Then XW int ∈ RN×H . The output of input layer is defined by
X int = σ(XW int) ∈ RN×H , (2)
where σ(·) is the activation function. Generally, ReLU(x) = max(0, x) is used as the nonlinear activation function. Here, each column of X int is regarded as a point in RN , so X int is the set of H points in RN and arranged in a certain order. The vector composed of the ith component of these points (the ith row of X int) is a feature representation of the ith node of graph G.

heading: 3.2 IFS LAYER
text: Let A be the adjacency matrix of G. Let triu(A) denote the upper triangular matrix of A and tril(A) denote the lower triangular matrix of A. The symmetric normalization of triu(A) and tril(A) are
A0 = D − 12 0 triu(A)D − 12 0 and A1 = D − 12 1 tril(A)D − 12 1 ,
where D0 and D1 are degree matrices of triu(A) and tril(A) respectively. Sometimes, we use the mean-pooling of triu(A) and tril(A), i.e. A0 = D−10 triu(A),A1 = D −1 1 tril(A). Let f0, f1 be the two affine transformations on RN , induced by A0,A1 respectively, defined as follows: f0 : x→ A0x+ b0, x ∈ RN , b0 ∈ R, f1 : x→ A1x+ b1, x ∈ RN , b1 ∈ R,
where b0 and b1 are learnable biases, namely add constants b0 and b1 to each component of A0x and A1x respectively. Constructing iterated function system
IFS = {RN ; f0, f1;p}, where p = (p0, p1) is a learnable adjoint probability vector, satisfying p0 > 0, p1 > 0 and p0+p1 = 1. Using the symbol space Ωm = {0, 1}m, then for each i = (i1, i2, ..., im) ∈ Ωm the length of i is m, denoted as |i| = m, and defining pi = pi1pi2 · · · pim and fi = fi1 ◦ fi2 ◦ · · · ◦ fim . Let n be the number of iterations of IFS. For IGNNS, n is a preset parameter. The iterative process of IFS is described as follows:
The first iteration (|i| = 1). The result of the first iteration is denoted by H(1) = {f0(X int), f1(X int)} = {Hi}|i|=1,
where Hi = fi(X int), ∀i ∈ Ω1. Since IFS selects the iteration branch fi with probability pi, the mathematical expectation of H(1) is computed by
E1 = p0f0(X int) + p1f1(X int) = p0H0 + p1H1 = ∑ |i|=1 piHi.
If choose to use bias in iterations, then H0 = A0X int + b0,H1 = A1X int + b1, where b0 and b1 are learnable H-dimensional vectors.
The second iteration (|i| = 2). Using the results of the first iteration as the input of the second iteration, then the result of the second iteration is denoted by
H(2) = {f0(f0(X int)), f0(f1(X int)), f1(f0(X int)), f1(f1(X int))} = {f00(X int), f01(X int), f10(X int), f11(X int)} = {Hi}|i|=2,
where Hi = fi(X int), ∀i ∈ Ω2. Note that IFS selects the iteration path fi with probability pi, then the mathematical expectation of H(2) is computed by
E2 = ∑ |i|=2 pifi(X int) = ∑ |i|=2 piHi.
We expand the expression of E2 and perceive its powerful feature representation ability. First,
H00 = f0(H0) = A0(A0X int + b0) + b0,H01 = f0(H1) = A0(A1X int + b1) + b0,
H10 = f1(H0) = A1(A0X int + b0) + b1,H11 = f1(H1) = A1(A1X int + b1) + b1.
Then E2 = p00H00 + p01H01 + p10H10 + p11H11
= (p00A00 + p01A01 + p10A10 + p11A11)X int
+ (p00A0b0 + p01A0b1 + p10A1b0 + p11A1b1) + (p00b0 + p01b0 + p10b1 + p11b1),
where Ai = Ai1Ai2 , ∀i = (i1, i2) ∈ Ω2. The n-th iteration (|i| = n). Inductively, we have
H(n) = {Hi}|i|=n, Hi = fi(X int), En = ∑ |i|=n piHi. (3)
Note that Hi ∈ RN×H and each column of Hi is a point in RN , so we regard it as a subset of RN withH elements. Thus H(n) is a subset of RN withH×2n elements (including duplicate elements). Because of Theorem 4.1, we call H(n) the fractal representation with depth n of notes.

heading: 3.3 REPRESENTATION LAYER
text: After n iterations of IFS layer, the dynamic trajectory of IFS is obtained:
O = {E1,E2, ...,En}.
In general, the global representation R of notes is obtained by time average or concatenation operations on O.
R = 1
n n∑ i=1 Ei ∈ RN×H or R =‖ni=1 Ei ∈ RN×nH ,
where ‖ is the concatenation operator. Because of Theorem 4.2, we call R the ergodic representation of notes. In practice, we adopt weighted time average or weighted concatenation. According to the Theorem C.1, we use heuristic weights (Here, we understand it as the average expansion factor of the distance between two points after affine transformation). Let r = √ ln(N) + γ, where γ ≈ 0.577215664 is the Euler constant. Suppose r = (r1, r2, ..., rn) is a learnable n-dimensional vector with initial value ri = ( 1 r )i−1 . Then the ergodic representation of notes is
R = n∑
i=1
riEi ∈ RN×H or R =‖ni=1 riEi ∈ RN×nH .

heading: 3.4 OUTPUT LAYER
text: Let W out ∈ RH×P be a learnable parameter matrix, where P is the dimension of the output layer (such as the number of class labels). If R is generated byO concatenation, then let W out ∈ RnH×P . There are two ways to construct output layer, one is to use a Single-Layer Perception (SLP) as output, that is
O = RW out + bout;
the other is to use f0, f1 for Mixed Propagation (MP), that is, let R0 = f0(RW out) and R1 = f1(RW out), where the biases of f0, f1 are removed, then the output
O = p0R0 + p1R1 + bout.
Where the bias bout ∈ RP is an optional learnable parameter vector.

heading: 3.5 INITIALIZATION OF LEARNABLE VARIABLES
text: The learnable parameters of IGNNS include input layer matrix W int ∈ RF×H , adjoint probability vector p = (p0, p1) ∈ R2 of IFS, biases b0, b1 ∈ RH of IFS layer, weight coefficient r = (r1, r2, ..., rn) of representation layer, matrix W out ∈ RH×P of output layer and bias bout ∈ RP of output layer. Among them, W int and W out are the required learnable parameters, using the initialization described in Glorot & Bengio (2010); b0, b1 and bout are optional learnable parameters with an initial value 0; p is a optional learnable parameter, for undirected graph, setting p0 ∈ [0.5− 0.1, 0.5 + 0.1], and for directed graph, setting (For the reasons, see Appendix G)
p0 = detD1
detD0 + detD1 , p1 = detD0 detD0 + detD1 ;
r is a optional learnable parameter, and its initial value as defined in 3.3. Let n be the number of iterations of IFS. We regard n as the depth of IGNNS. Thus, IGNNS is denoted as
O = IGNNS(X,A;W int, n,p, b0, b1, r,W out, bout) or simply O = IGNNS(X, IFS),
where IFS is induced by relational matrix A. The output of IGNNS can be used as the input of downstream tasks, and can also be connected to other network architectures.

heading: 3.6 THEORETICAL TIME COMPLEXITY OF IGNNS
text: Let n,N,H, P be defined as above. Let T (·) denote the number of calculations of an object. For input layer,
T (input layer) = NFH +NH = O(NFH).
For IFS layer, during the iterative calculation, we will store the previous calculation results, thus T (H(1)) = 2N2H , T (H(i)) = 2×T (H(i−1)), i = 2, 3, ..., n. It follows that T (H(i)) = 2iN2H, i = 1, 2, ..., n. Similarly, T ({pi||i| = i}) = 2i, i = 1, 2, ..., n. Complete the above calculation, it is easy to see that T (Ei) = 2iNH . Thus
T (IFS layer) = n∑
i=1
( T (H(i)) + T ({pi||i| = i}) + T (Ei) ) = O(2nN2H).
It is easy to verify that T (representation layer) = O(nNH). We assume that the output layer is constructed by mixed propagation, then T (output layer) = O(N2P ) if W out ∈ RH×P , and T (output layer) = O(N2P + nNHP ) if W out ∈ RnH×P . Then
T (IGNNS) = O(2nN2H +N2P ) or O(2nN2H +N2P + nNHP ). In practice, for large graphs, 2nN2H N2P nNHP , thus 2nN2H is the main factor affecting the time complexity of IGNNS. Furthermore, for large graphs of the same size, n is the main important factor affecting time complexity. For citation network datasets such as citeser, Cora and PubMed, we suggest that n ≤ 8 (see Appendix B).

heading: 4 GEOMETRIC PROPERTIES OF IGNNS
text: The discussion here assumes that affine f0, f1 are contractive. Otherwise, let f0 : x → 1 ‖A0‖F+1A0x + b0 and f1 : x → 1
‖A1‖F+1A1x + b1. In practice, IGNNS does not use contractive affine. If contractive affine is used in IGNNS, it can be seen from the following theorems that the characterization ability of IGNNS decreases with the increase of IFS iterations, which is similar to the performance of Graph Convolution Network (GCN). Such a phenomenon is called over-smoothing (Li et al., 2018b; Xu et al., 2019; Chen et al., 2020), which suggests that as the number of layers increases, the representations of the nodes in GCN are inclined to converge to a certain value and thus become indistinguishable. In order to overcome the over-smoothing problem of deep GNN, people need to use some new methods such as Skip connection (Xu et al., 2018), Drop edge (Rong et al., 2020), Residual connection (Klicpera et al., 2019a), Identity mapping (Chen et al., 2020), Generalized message aggregation functions (Li et al., 2020) and so on. Generally speaking, deep network may lead to the decrease of generalization performance. To analyze which type of deep GNN would achieve better generalization performance, Xu et al. (2020) proposes a guiding theoretical framework.
Theorem 4.1 (Fractal generation) Let H(n) = {Hi}|i|=n, which is a subset of RN with H × 2n elements (including duplicate elements), then
dH(H(n),F)→ 0, n→∞, where dH is the Hausdorff distance defined on H(RN ), the set of all nonempty compact subsets of RN , and F is the fractal set of IFS in IGNNS. In other words, as the number of iterations increases, H(n) will be independent of node feature X , only related to the graph structure described by A.
Let T be the Hutchinson operator on H(RN ), defined as T (B) = f0(B) ⋃ f1(B),∀B ∈ H(RN ). Then the updated rule of H(n) satisfies H(n) = T (H(n−1)) = · · · = Tn(H(0)),
where H(0) = {X int} is a subset of RN with H elements. In fractal geometry, H(n) is used to draw the fractal image on the plane. First, taking initial value H(0) = {x0}, where x0 is a point in plane. For enough n, printing all the points of H(n) on the screen to obtain the approximate fractal image.
Theorem 4.2 (Ergodic property) Let En = ∑ |i|=n piHi be the mathematical expectation of H(n) = {Hi}|i|=n, then En converges to a constant matrix E ∈ RN×H in Frobenius norm, i.e. limn→∞En = E, where E[i, :] = (ei, ei, ..., ei)> ∈ RH and ei ∈ R is a constant, i = 1, 2, ..., N . Furthermore, the time average of the dynamic trajectory O of IFS satisfies
lim n→∞
1
n n∑ i=1 Ei = lim n→∞ En = E
and series ∑∞
i=1 riEi ∈ RN×H converges in Frobenius norm.
Theorem 4.2 shows that as long as the number of iterations is large enough, the embeddings of nodes will be close to linear correlation, and the representation ability of IGNNS will decline. However, in the framework of IGNNS, because the spectral radius ρ(A0) = ρ(A1) = 1, IFS is not contractive in general, and IGNNS still has the ability of depth feature representation.

heading: 5 EXPERIMENTS

heading: 5.1 EXPERIMENTAL TASK: SEMI-SUPERVISED NODE CLASSIFICATION
text: Let Z = softmax(O), where Z ∈ RN and softmax(·) is the softmax activation function, defined as softmax(xi) = 1Z exp(xi) with Z = ∑ i exp(xi), is applied row-wise. For semi-supervised multiclass classification, we employ the following cross-entropy to evaluate error over all labeled examples: L = − ∑ l∈YL ∑P i=1 Y [l, i] lnZ[l, i],where YL is the set of node indices that have labels with P classes, Y [l, :] is a one-hot vector of size P representing the class of node l and Z[l, :] is the row l of the matrix Z.

heading: 5.2 EXPERIMENTAL SETUP
text: Datasets. In our experiment, we use three standard citation network benchmark datasets for evaluation, including Cora, Citeseer, Pubmed and apply the standard fixed training/validation/testing split (Yang et al., 2016; Kipf & Welling, 2017; Veličković et al., 2018) on above datasets, with 20 nodes per class for training, 500 nodes for validation and 1,000 nodes for testing. In these citation networks, papers are represented as nodes, and citations of one paper by another are denoted as edges. Node features are the bag-of-words vector of papers, and node label is the only one academic topic of a paper. See Table 1 for more details.
Parameter Setting. Random seed for Tensorflow and Numpy is set to 1234. ReLU (Nair & Hinton, 2010) is used as the activation function in input layer and output layer. Dropout (Srivastava et al., 2014) is applied to input layer, IFS layer and output layer. In representation layer, we adopt weighted time average to get the global representation of notes. In output layer, we adopt the method of mixed propagation to get the output of IGNNS. We use the AdamOptimizer (P.Kingma & Ba, 2015) during training. More details of hyper-parameters are shown in Table 2. During training stage, we select the best model to maximize the accuracy of the validation set and use early stopping with a patience of 100 epochs.

heading: 5.3 EXPERIMENTAL RESULT
text: We compare with those models that strictly follow the standard of experiment setup of semisupervised node classification, i.e. the standard fixed training/validation/testing split (Yang et al., 2016; Kipf & Welling, 2017) is applied on dataset. For baselines, we include recent deep GNN models such as JKNet (Xu et al., 2018), APPNP (Klicpera et al., 2019a), Attention-based models such as GAT (Veličković et al., 2018), AGNN (Thekumparampil et al., 2018) and H-GAT (Gulcehre et al., 2019) , and other models such as TAGCN (Du et al., 2017) and N-GCN (Abu-El-Haija et al., 2018). We also include three state-of-the-art shallow GNN models: Planetoid (Yang et al., 2016), GCN (Kipf & Welling, 2017) and DGCN (Zhuang & Ma, 2018). The detailed results are shown in Table 3.
We can see from Table 3 that the improved performance of model IGNNS in dataset Cora and Citeseer is much higher than that in dataset Pubmed. To understand why this happens, we analyze
the characteristics of these citation networks. We consider two statistical properties of networks, one is the Network Density d(G), which is defined as d(G) = 2LN(N−1) , where N is the number of nodes and L is the number of edges, and the other is the Average Clustering Coefficient C, which is defined as C = 1N ∑ i∈V Ci, where V is the set of nodes, Ci = 2ei ki(ki−1) , ki is the number of the neighbors of node vi and ei is the number of undirected edges between ki neighbors. The small Network Density means the strong global sparsity of the network, and the small Average Clustering Coefficient means the strong sparsity of the neighbors of nodes. The calculation results of the statistical characteristics of the network are shown in Table 4. We can see from Table 4 that Pubmed is more sparse than Cora and Citeseer. The performance of IGNNS benefits from the bidirectional mixed propagation of information between nodes, but this sparsity weakens the gain of IGNNS.

heading: 5.4 PERFORMANCE OF COMPLETELY LINEAR IGNNS
text: In Nonlinear IGNNS, we use the nonlinear activation function ReLU(x), learn adjoint probability vector p = (p0, p1) by p0 ← ReLU(p0)+0.1ReLU(p0)+ReLU(p1)+0.2 , p1 ← ReLU(p1)+0.1 ReLU(p0)+ReLU(p1)+0.2 and learn the rep-
resentation layer coefficient r = (r1, r2, ..., rn) by ri ← ReLU(ri) with initial value ri = ( 1 r )i−1 where r = √ ln(N) + 0.577215664. In this experiment, to get a completely linear IGNNS, we let all the activation functions be the identity function, i.e. σ(x) = x, and let the adjoint probability vector p = (p0, p1) and the representation layer coefficient r = (r1, r2, ..., rn) be hyperparameters
without learning. For Citeseer, we let p0 = 0.6, and for Pubmed, we use bias for IFS layer. Except for the above changes, experimental task and the other experimental setups remain unchanged as showed in section 5.1 and 5.2 respectively. We can see from Table 5 that the performance of completely linear IGNNS is better than that of baseline model GCN. Compared with other models, completely linear IGNNS is still competitive. This is due to the fact that the IFS can extract more features than spectral filters. For more discussion, Let A be the normalization adjacency matrix of graph G, and let A0 and A1 be defined as in section 3.2. We further assume that the dimension of the hidden space is equal to 1. This means that the input of GNN (GCN or IGNNS) is a point x0 = X
int = XW int ∈ RN×1. Let n be the depth of GNN, for IGNNS, equal to the number of iterations of IFS. For convenience, we ignore the activation function and parameter matrix. Let f(x) = Ax+ b, f0(x) = A0x+ b0 and f1(x) = A1x+ b1. For GCN, the message passing results are {f(x0)}, {f2(x0)}, ..., {fn(x0)}. We see that each iteration only gets one value, i.e. |{fn(x0)}| = 1. For IGNNS, the message passing results are as follows:
{f0(x0), f1(x0)}, {f0(f0(x0)), f0(f1(x0)), f1(f0(x0)), f1(f1(x0))}, ..., {fi(x0)}|i|=n.
If f0 and f1 satisfy separation condition, i.e. f0(x) = f1(y) implies x = y, then |{fi(x0)}|i|=n| = 2n. This means that IGNNS can extract more information than GCN. Even if f0 and f1 are contractive mappings, by Theorem 4.1, we have {fi(x0)}|i|=n → F, where F is the fractal set of IFS induced by f0 and f1. Generally speaking, F is a uncountable compact set. This means that when n is large, the features may still be distinguishable.

heading: 6 CONCLUSION
text: In this paper, we propose a new framework of graph neural networks, IGNNS, which give a connection between Graph Neural Networks and Iterated Function System. We use IFS to simulate the bidirectional message passing process of graph neural network, and obtain the fractal representation and ergodic representation of graph nodes, which are very helpful for downstream tasks. The experiments show that we have achieved good results in semi-supervised node classification task. Interesting directions for future work include pruning the iterative path space {0, 1}n to reduce the computational complexity, coding graph structured data with IFS, and establishing more interesting connections between IFS and graph neural networks.

heading: A THE FRACTAL REPRESENTATION OF GRAPH G WITH ONLY ONE SELF
text: ADJACENT NODE v.
For the sake of discussion, we assume that the dimension of the hidden space is equal to 1. We assume that messages are sent from node v, propagate in two directions (clockwise and anticlockwise), and are finally received by node v. The received messages in the clockwise direction become one-third of the original, and the received messages in the anticlockwise direction become one-third of the original plus a constant of 23 . It is expressed by mathematical formula as follows
f0(x) = x
3 , f1(x) =
x 3 + 2 3 , x ∈ R.
For Bi-GCN, messages are delivered independently in both directions. In other words, there are two independent channels, and the message passing (transmitting or receiving) can only be carried out by their own channels. Let x0 ∈ R be the initial message. In the clockwise direction, after n passes, the received messages are
x0 3 , x0 32 , ..., x0 3n → 0.
In the anticlockwise direction, the received messages are
x0 3 + 2 3 , x0 32 + 2 32 + 2 3 , ..., x0 3n + n∑ i=1 2 3i → 1.
For IGNNS, the two channels have a connection point at node v. First, node v sends the message x0 in both directions, and the connection point of node v will receive two messages {f0(x0), f1(x0)}. In the second launch, any message (f0(x0) or f1(x0)) can be sent in both directions, so the received messages are
{f0(f0(x0)), f0(f1(x0)), f1(f0(x0)), f1(f1(x0))}.
In summary, After n passes, the received messages are
H(n) = {fi(x0)}|i|=n → C,
where C is the famous Cantor Set. This means that we have not only received one message, but 2n, since f0 and f1 satisfy the separation condition. We can see from Figure 3 that Bi-GCN gets boundary messages and IGNNS gets all messages. Let p = (p0, p1) be the adjoint probability vector, then the mathematical expectation En of H(n) is ∑ |i|=n pifi(x0). We interpret En as the average of all messages received. Now the question is, fractal representation gets enough messages, but is there redundancy in these messages? How to select the valid message from the fractal representation becomes the focus of our research in the next stage.

heading: B ANALYSIS OF TIME COMPLEXITY ON CORA, CITESEER AND PUBMED.
text: From section 3.6, we have known that the time complexity in Experiment 5.1 isO(2nN2H+N2P ), where n is the number of iterations of IFS (the depth of IGNNS), N is the number of nodes, H is the dimension of the latent space and P is the dimension of the output layer. In this section, we compare the real running time (100 epochs) of IGNNS on Cora, Citeseer and Pubmed. Let H = 8, then 2nN2H is the main factor affecting the time complexity of IGNNS. Let the depth of IGNNS
go from 1 to 8, and compare the real training time on Cora, Citeseer and Pubmed. We can see from Figure 4 that the actual results are basically consistent with the theoretical analysis results.

heading: C FROBENIUS NORM OF MATRIX
text: Theorem C.1 Let A ∈ RN×N be the adjacency matrix of a graph with no weights, i.e. Ai,j = 1 if there exists an edge i → j in the graph and Ai,j = 0 otherwise, and A1 = D−11 tril(A) (or A1 = D − 12 1 tril(A)D − 12 1 ) as defined in IGNNS, then
‖A1‖F ≥ √ ln(N) + γ,
where γ ≈ 0.577215664 is the Euler constant.
Proof. Case 1: A1 = D−11 tril(A). Let tril(A) = (aij)N×N , D1 = diag(d1, d2, ..., dN ) be the degree matrix of tril(A) and A1 = (bij)N×N . Note that tril(A) is an lower triangular matrix, then di = ∑N j=1 aij = ∑i j=1 aij . Since aij ∈ {0, 1}, we have i ≥ di. Computing the Frobenius Norm of A1 as follows:
‖A1‖2F = N∑ i=1 N∑ j=1 b2ij = N∑ i=1 i∑ j=1 b2ij = N∑ i=1 i∑ j=1 ( aij di )2 . (4)
∀i ∈ {1, 2, ..., N}, note that di elements in {aij}ij=1 are 1 and the rest are 0. It follows that
i∑ j=1 ( aij di )2 = ( 1 di )2 × i∑ j=1 a2ij = ( 1 di )2 × di = 1 di ≥ 1 i . (5)
It follows from (4) and (5) that
‖A1‖2F ≥ N∑ i=1 1 i .
So ‖A1‖F ≥ √∑N i=1 1 i ≈ √ ln(N) + γ, where γ ≈ 0.577215664 is the Euler constant.
Case 2: A1 = D − 12 1 tril(A)D − 12 1 . Computing the Frobenius Norm of A1 as follows:
‖A1‖2F = N∑ i=1 N∑ j=1 b2ij = N∑ i=1 i∑ j=1 b2ij = N∑ i=1 i∑ j=1 a2ij didj . (6)
∀i ∈ {1, 2, ..., N}, it follows from j ≥ dj that
i∑ j=1 a2ij didj ≥ 1 di i∑ j=1 1 j · a2ij . (7)
Note that di elements in {aij}ij=1 are 1 and the rest are 0. It follows from Rearrangement inequality that
i∑ j=1 1 j · a2ij ≥
1
i− di + 1 · a2i(i−di+1) +
1
i− di + 2 · a2i(i−di+2) + · · ·+
1 i · a2ii, (8)
where ai(i−di+1) = ai(i−di+2) = · · · = aii = 1. Thus
i∑ j=1 1 j · a2ij ≥ 1 i · 12 + 1 i · 12 + · · ·+ 1 i · 12 = di × 1 i . (9)
It follows from (6), (7) and (9) that
‖A1‖2F ≥ N∑ i=1 1 i .
Which completes the proof.
D INTRODUCTION TO ITERATED FUNCTION SYSTEM
In order to prove Theorem 4.1 and Theorem 4.2, we will briefly introduce the relevant conclusions on IFS in this section, and we will not give the proof here. More details of IFS Theory can be found in Hutchinson (1981); Elton (1987); Barnsley (1988); Falconer (1990); Massopust (2017).
D.1 FRACTAL SPACE
Let (X; d) be a complete metric space. Let H(X) denote a set consisting of all nonempty compact subsets of X . Hausdorff distanceon dH onH(X) defined by
dH(A,B) = max{max a∈A min b∈B d(a, b),max b∈B min a∈A d(a, b)},∀A,B ∈ H(X). (10)
Theorem D.1 (H(X); dH) is a complete metric space.
We call (H(X); dH) a Fractal space. Let {fi}ni=1 be a set of mappings on (X; d). Hutchinson operator T : (H(X); dH)→ (H(X); dH) defined as
T (B) = n⋃ i=1 fi(B),∀B ∈ H(X). (11)
Theorem D.2 If {fi}ni=1 is a set of contractive mappings on (X; d), then Hutchinson operator T is a contractive mapping on (H(X); dH).
D.2 MARKOV OPERATOR OF IFS
Let (X; d) be a complete metric space. Let M(X) be the set of all probability measures on X . Let C(X) be the set of all continuous functions mapping X to R. We say that f ∈ L̃ip1, if |f(x)− f(y)| ≤ d(x, y),∀x, y ∈ X . It is easy to see that if f ∈ L̃ip1 then f ∈ C(X). Hutchinson metric dM onM defined as
dM (µ, ν) = sup {∣∣∣∣∫ X fdµ− ∫ X fdν ∣∣∣∣ |f ∈ L̃ip1} ,∀µ, ν ∈M(X). (12) Theorem D.3 (M(X); dM ) is a complete metric space.
Let IFS = {X; f1, f2, ..., fn;p}, the Markov operator M :M(X)→M(X) of IFS defined as
Mµ = n∑ i=1 piµ ◦ f−1i , µ ∈M(X). (13)
Theorem D.4 Markov operator M of IFS is a contractive mapping on space (M(X); dM ).
Let measure sequence {µi} ⊂ M and µ ∈M, we call {µi} weakly convergent to µ if the following equation holds:
lim i→∞ ∫ X fdµi = ∫ X fdµ,∀f ∈ C(X), (14)
denoted as µi w−→ µ as i→∞.
Theorem D.5 If µi dM−−→ µ as i→∞, then µi w−→ µ as i→∞.

heading: E THE PROOF OF THEOREM 4.1
text: Proof. By Theorem D.1, (H(RN ); dH) is a complete metric space. Let T be the Hutchinson operator on H(RN ), defined as T (B) = f0(B) ⋃ f1(B),∀B ∈ H(RN ). By Theorem D.2, T is a contractive mapping on (H(RN ); dH). It follows from the Banach fixed point theorem that there exits a unique compact subset F ofH(RN ) such that
F = T (F) = f0(F) ⋃ f1(F),
which implies that F is the fractal set of IFS. Further more, ∀B ∈ H(RN ), we have Tn(B) dH−−→ F. The above convergence is independent of the choice of initial value. Thus, let H(0) = {X int} = {x1, x2, ..., xH}, xi ∈ RN , we have
H(n) = T (H(n−1)) = · · · = Tn(H(0)) dH−−→ F. The above result indicates that when n is large enough, H(n) is close to the fractal set F of IFS in the sense of Hausdorff distance, and has nothing to do with the choice of initial value H(0).

heading: F THE PROOF OF THEOREM 4.2
text: Proof. It suffices to prove that ∀j ∈ {1, 2, ...,H} , En[:, j], the j column of En, satisfies
lim n→∞ En[:, j] =  e1 e2 ... eN  N×1 .
For this purpose, let xj be the j column of X int as defined in (2), then xj is a point in RN . Define a Dirac measure as follows:
δx(B) = { 1 x ∈ B, 0 other.
(15)
It is easy to see that δx ∈ M(RN ). The Markov operator M :M(RN )→M(RN ) of IFS, defined as
Mµ = 1∑ i=0 piµ ◦ f−1i , µ ∈M(R N ). (16)
Now take µ0 = δxj , and the results of iterative calculation are as follows:
µ1 = Mµ0 = 1∑ i=0 piµ0 ◦ f−1i = 1∑ i=0 piδxj ◦ f−1i = 1∑ i=0 piδfi(xj) = ∑ |i|=1 piδfi(xj).
µ2 = M 2µ0 = Mµ1 = 1∑ i=0 piµ1 ◦ f−1i = 1∑ i=0 pi( ∑ |i|=1 piδfi(xj)) ◦ f −1 i
= p0p0δf0(f0(xj)) + p0p1δf0(f1(xj)) + p1p0δf1(f0(xj)) + p1p1δf1(f1(xj)) = ∑ |i|=2 piδfi(xj)
Inductively, we have
µn = M nµ0 = ∑ |i|=n piδfi(xj). (17)
By Theorem D.4, it follows from the Banach fixed point theorem that there exits a unique probability measure µ∗ such that
µn dM−−→ µ∗.
The above µ∗ is actually the self-similar measure of IFS. By Theorem D.5, we have µn w−→ µ∗, i.e.
lim n→∞
∫ Fdµn = ∫ Fdµ∗, ∀F ∈ C(RN ).
It follows from (17) and (3) that
∫ Fdµ∗ = lim
n→∞
∫ Fdµn = lim
n→∞
∫ Fd( ∑ |i|=n piδfi(xj))
= lim n→∞ ∑ |i|=n piF (fi(xj))
= lim n→∞ ∑ |i|=n piF (Hi[:, j]), ∀F ∈ C(RN ). (18)
In (18), ∀i ∈ {1, 2, ..., N}, take the continuous function Fi to satisfy
Fi(t) = ti, ∀t = (t1, t2, ..., tN )> ∈ RN .
It follows from (18) that ∀i ∈ {1, 2, ..., N},
lim n→∞ En[i, j] = lim n→∞ ∑ |i|=n piHi[i, j] = lim n→∞ ∑ |i|=n piFi(Hi[:, j])
= ∫ Fi(t)dµ∗(t) = ∫ tidµ∗(t) = ∫ F tidµ∗(t) = ei. (19)
Which combined with simple mathematical analysis technology completes the proof.

heading: G HOW TO SET THE INITIAL VALUE OF ADJOINT PROBABILITY VECTOR?
text: The geometric meaning of matrix determinant detA is the expansion factor of graph volume under linear transformation A. Let F be the fractal set of IFS as defined in IGNNS. Then
F = f0(F) ⋃ f1(F).
Note that
detAi = volume(fi(F))
volume(F) , i = 0, 1.
It can be seen that if the value of detAi is large, it reflects that fi(F) has a large share in F. therefore, when selecting the iterative function, fi should have a greater probability of being selected. So set
p0 = detA0
detA0 + detA1 , p1 = detA1 detA0 + detA1 ;
Note that A0, A1 are triangular matrixes and the diagonals of triu(A), tril(A) are equal to 1. It follows from A0 = D−10 tril(A) (or A0 = D − 12 0 tril(A)D − 12 0 ) that detA0 = 1 detD0
. Similarly, detA1 =
1 detD1 . It follows that
p0 = detD1
detD0 + detD1 , p1 = detD0 detD0 + detD1 ,
where D0 and D1 are degree matrices of triu(A) and tril(A) respectively.

Review summaries:
This work proposes a new graph neural network architecture with modified rules for message passing, Iterated Graph Neural Network System (IGNNS). The paper then provides a theoretical analysis of the proposed architecture by connecting it with Iterated Function System (IFS), an important research field in fractal geometry. This paper further demonstrates empirically that the proposed architecture outperforms related models on citation network datasets.

This paper proposes a new framework of GNN which can deal with undirected and directed graphs in a unified way. The authors argue that the size of the symbol space for a message passing path with length n is 2^n, while previous architectures only have constant size. Motivated by this observation, the authors borrow ideas from Iterated Function System to augment the symbol space. 

The paper proposes a new definition of GNNs designed to cope with bi-directional message-passing processes.  To do so, a new symbols space, different from the one adopted by Bidirectional GCN,  is considered, together with an iterated function system. These lead to an architecture composed of 4 steps: an input layer that acts as a classic FC layer; an IFS layer that applies the iterated function system considering the adjacency matrix; a layer to concatenate or sum the expected values of each iteration; and an output layer that combines the results using the functions of the IFS and a new learnable weight matrix.