Stable, Efficient, and Flexible Monotone Operator Implicit Graph Neural Networks
Authors: Justin Baker
Qingsong Wang
Bao Wang
Venue: Submitted to ICLR 2023 (2023)
Keywords: implicit graph neural networks
monotone operator
accelerated operator splitting
orthogonal parameterization

Abstract: Implicit graph neural networks (IGNNs) that solve a fixed-point equilibrium equation for representation learning can learn the long-range dependencies (LRD) in the underlying graphs and show remarkable performance for various graph learning tasks. However, the expressivity of IGNNs is limited by the constraints for their well-posedness guarantee. Moreover, when IGNNs become effective for learning LRD, their eigenvalues converge to the value that slows down the convergence, and their performance is unstable across different tasks. In this paper, we provide a new well-posedness condition of IGNNs leveraging monotone operator theory. The new well-posedness characterization informs us to design effective parameterizations to improve the accuracy, efficiency, and stability of IGNNs. Leveraging accelerated operator splitting schemes and graph diffusion convolution, we design efficient and flexible implementations of monotone operator IGNNs that are significantly faster and more accurate than existing IGNNs.

Under review as a conference paper at ICLR 2023
STABLE, EFFICIENT, AND FLEXIBLE MONOTONE OPERATOR IMPLICIT GRAPH NEURAL NETWORKS
Anonymous authors Paper under double-blind review
ABSTRACT
Implicit graph neural networks (IGNNs) that solve a fxxe-point equilibrium equa- tion for representation learning can learn the long-range dependencies (LRD) in the underlying graphs and show remarkable performance for various graph learn- ing tasks. However, the expressivity of IGNNs is limited by the constraints for their well-posedness guarantee. Moreover, when IGNNs hecome effective for learning LRD, their eigenvalues converge the value that slows down the con- vergence, and their performance is unstable across different tasks. In this pa- per, we provide a new well-posedness condition of IGNNs leveraging monotone operator theory. The new well-posedness characterization informs us to design effective parameterizations to improve the accuracy, efficiency, and stability IGNNs. Leveraging accelerated operator splitting schemes and graph diffusion convolution, we design efficient and flexible implementations of monotone oper- ator IGNNs that significantly faster and more accurate than existing IGNNs.
INTRODUCTION
Implicit graph neural networks (IGNNs) that solve a fixed-point equilibrium equation for graph representation learning can learn long-range dependencies (LRD) in the underlying graphs, showing remarkable performance for various tasks [69 39 58 63 22). Let G = (V, E) represent graph, where V is the set of nodes, and E C V V is the set of edges. The connectivity of can be represented by the adjacency matrix A E R with Aij = 1 if there is edge connecting nodes i E V; otherwise Aij = 0. Let X E Rdxn be the initial node features whose i-th column æ; E Rdd is the initial feature of the i-th node. IGNN [39] learns the node representation by finding the fixed point, denoted as Z", of the Picard iteration below
z a (* + (X)), for k 1,2,
where is the nonlinearity (e.g. ReLU), gB is function parameterized by B (e.g. (x(X) Il
Acourseytt)
BXG), matrices W and B Rd are learnable weights, and G is graph-related matrix. In IGNN, Gs chosen as Á = D- 0-112((T (I A)D-1/2 with
105
200 295
600
500
200
202
300
400
500
kco
being the identity matrix and D is the degree ma- trix with D2 1+Σ;= Ag. IGNN constrains W
Acomscytig
5 5
using a tractable projected gradient descent method
Trai
to ensure the well-posedness of Picard iteration at the cost of limiting the expressivity of IGNNs. The prediction of IGNN is given by fe (z"),, a func- tion parameterized by e. IGNNs have several mer- its: 1) The depth of IGNN is adaptive to partic- ular data and tasks rather than fixed. 2) Training
102
202
200
700
200
Goo
Epoch
Epoch
Figure 1: Epoch vs. training. validation, test accuracy of IGNN for classifying directed chains. First row: binary chains of length 100 (left) and 250 (right). Sccond row: three-class chains of length 80 (left) and 100 (right).
IGNNs requires constant memory independent of their depth - leveraging implicit differentiation 666. 5l 13) 3) IGNNS have better potential to capture LRD of the underlying graph compared
to existing GNNs, including GCN 755. GAT 733, SSE (23), and SGC
[79]. The latter GNNs lack
the capability to learn LRD as they suffer from over-smoothing [56
84
62 20 Several methods
have been proposed to alleviate over-smoothing and hence improve learning LRD by adding resid-
ual connections [77 [1] by geometric aggregation [65]. by adding fully-adjacent layer [3, by
improving breadth-wiss backpropagation (59]. and by adding oscillatory layers 227 67

---

Under review as a conference paper at ICLR 2023
Issue 1: Well-posedness of IGNN Limits Its Expressivity. One bottleneck of IGNN is that the magnitude of W's cigenvalues has to be less than one for its well-posedness guarantee; see Sec. 2 for details. This limits the selection of W and thereby limits the expressivity of IGNNs.
Issue 2:
When can IGNNS learn LRD? To understand when IGNN can learn LRD, we run IGNN using the settings in [3 to classify directed chains. Directed chains is synthetic dataset designed
to test the effectiveness of GNNs in learning LRD for node classification 71 39] Fig. plots epoch vs. accuracy of IGNN for the chain classification. Here, each epoch means iterating Equation (I) until convergence and then updating W and B at the end. IGNN can classify the binary chain task perfectly at length 100 but performs near random guesses when the length is 250, as illustrated in [ For the three-class chains. IGNN's s performance is very poor at chain length 100 but performs quite well at length 80. We investigate the results above by studying the dynamics of eigenvalues of the matrix |WD For illustrative purpose, we consider A1 [W]) and )2 WI), the largest and the second largest eigenvalue of [W| in magnitude. Fig. 2 (left) contrasts the evolution of the magni- tude of X I1) and ( of IGNN when classitying nodes on chains with different lengths. We see that the magnitude of both eigenvalues goes to when IGNN becomes accurate. However, Fig. 2 (riihht shows that IGNN takes many more iterations in each epoch when the magnitude of eigenvalues gets close to 1. Indeed, when A >1 ((W)) + 1, the Lipschitz constant of the linear map W ZG + (n(X) is close I, slowing down the convergence of the Picard iterations. The results in Fig. 2 echo our intuition; the representation of a given node aggregates one more hop of infor- mation after each Picard iteration; when the magnitude of eigenvalues gets close to I, Equation (T) converges slowly so that IGNN can capture LRD before fixed point convergence.
We report the classification results of different lengths in Appendix D these results show prevalently that IGNNs suffer from two bottlenecks: 1) An inherent tradeoff between computational efficiency and capability for learning LRD. 2) The performance of IGNNs, based on Picard iteration, unsta- ble in the sense that their performance varies substantially across tasks. In particular, starting from random Gaussian initialization of W the default initialization of W IGNN cannot learn LRD if none of the eigenvalues of W get close to 1 in magnitude.
1.1 OUR CONTRIBUTION
1.00
6.70
150
We develop accurate, stable, and efficient monotone op- erator IGNNs (MIGNNs P
0:50
In particular, we derive new colldposedsess condition for MIGNN leveraging mono- tone operator theory; see See The new well-posedness condition informs us to design 1) a monotone param- eterization of W, whose eigenvalues can take much wider range than that of IGNNs, to boost the expressiv- ity of MIGNNs, addressing Issue 1. And 2) a Cayley transform-based orthogonal paraneterization of W to improve the stability and efficiency of MIGNN for learn- ing LRD, addressing Issue 2: see Sec. 3 Picard itera- tion is inefficient or impossible to find the fixed point of
100
200
40
30
1.00
-
220
Epoch vs. the magnitude of W]) and ()) and the iterations re- quired for each epoch. First row: binary chains, second row: Three-class chains.
Figure 2:
MIGNN with monotone or orthogonal parameterization. As such, we implement MIGNNs leverag- ing muderson-acceleratee operator splitting schemes; see Sec. 4 We verify the efficacy of MIGNN on various benchmark tasks; see Sec. 5
1.2 ADDITIONAL RELATED WORK
We briefly review some representative related works in three directions: deep equilibrium models (DEQs), GNNs, and orthogonal parameterizations for recurrent neural networks (RNNs).
DEQ. IGNN is related to DEQs 7 26 R but the equilibrium equation of IGNN differs from DEQs in that IGNN encodes graph structure. DEQs are class of infinite depth weight-tied feedforward neural networks with forward propagation using root-finding and backpropagation using implicit differentiation. As result, training DEQs only requires constant memory independent of the net- work's depth. Monotone operator theory has been used to guarantee the convergence of DEQS [77) and to improve the robustness of implicit neural networks [44] The convergence of DEQs has also
'The matrix W is obtained by taking the entry-wise absolute value of the matrix W. Starting from here, we use MIGNN stress that the model is based monotone operator theory.
2

---

Under review as a conference paper at ICLR 2023
been considered by constraining the network's weights 449 Jacobian regularization has been used to stabilize the training of DEQs (9. Anderson-accellrrate hyperparameters Graph neural networks. Classical GNNs are defined by stacking explicitly defined graph filter-
Linearized DEQs are studied in [46].
DEQs with learned acceleration
also proposed [T]]
ing layers. Examples include graph convolutional networks (GCNs) [17
24 48), recurrent GNNs
38 30
57 21 GraphSAGE [40], neural graph fingerprints (25 graph isomorphism network
(GIN)
80 message passing neural networks [36). graph attention networks (GATs) [73], GCNs
with convolution kernels learned based on paths (PAN 60] and pathGCN (28). and higher-order
message passing networks [S 14 There are some recent advances in IGNNs: EIGNN removes the nonlinearity in each intermediate iteration and derives closed form of the infinite iterations (58), convergent graph solver (CGS) is an IGNN model with convergence guarantees by constructing the input-dependent linear contracting iterative maps [63] GIND leverages implicit nonlinear diffusion
to access infinite hops of neighbors [22]. In addition to Picard iteration, implicit GNNs have also been defined by parametrizing the diffusion equation on graphs, e.g. [18 19).
Orthogonal parameterization for deep learning. The fixed point iteration Equation I is related the hidden state updates of RNNs [66 29: 2 [50]. Learning LRD is challenging for RNNs due
to exploding and vanishing gradient during backpropagation through time [76] 12 64). Enforcing
orthogonal parameterization for RNNs an effective approach to overcome exploding and vanishing
gradients, benefiting RNNs for learning LRD 55] 78
45 74 61
1.3 NOTATION
We denote scalars by lower- or upper-case letters and vectors/matrices by lower-/upper-catt boldface letters. For vector a, use ||a]|/]| |/]ac so to denotes its f2-/f<c-norm. We use I to denote the identity matrix whose dimension can be inferred from the context. For a matrix A, we denote its transpose as || |//|||||// ol and we denote
A its inverse A-1, its i-th largest eigenvalue in magnitude as (WW).. Given two matrices A and B, we denote their Kronecker/entry-wiss product as A & B/A @ B. and denote A B (A B) if A- B is positive definite (semi-positive definite). We use vec( to denote the vectorization of the matrix A in column-maiji order. The meaning of other notations can be inferred from the context.
,, its Frobenius normm22-nommooonoo
WELL-POSEDNESS OF MIGNN: A MONOTONE OPERATOR PERSPECTIVE
In this section, we characterize the well-posednes of MIGNN leveraging monotone operator theory, see Appendix B for brief review of monotone operator theory. Using the Kronecker produc 3l and vectorization matrix, we can rewrite Equation ( into the following equivalent vectorized form
vec (G @ Wvec(Z (*) + veclgn (X))).
(2)
Gu et al. propose the well poleepossss condition of IGNN as ( & WI) < 1, guaranteeing that 22) can be found by Picard iteration. Selecting G A. all the unique fixed point of Equation
eigenvalues of G are in [-1 with Àn (G) = 1. Therefore, well well-posedness of IGNN is equivalent A1 WT). Then, IGNN parameterizes W by relaxing the well-posednesss condition ((W) < 1 to W xo < 1, which constrains the magnitudes of eigenvalues of W to be less than 1.
to >a WT) < as ( & W = Ar ((G))(((W))
We seek to apply the monotone operator theory to improve the expressivity and efficiency of existing
IGNNs. According to the monotone operator theory (68] 77), finding the fixed point of Equation (2) is equivalent to solving the monotone inclusion problem: find 0 e (F (vec() Z)) with F and g
being two set-valued functions that are given below
F((vo(()) ()) (I G¹ @ W)vee(Z) vec(gn(X)) and af.
(3)
where O denotes the subgradient of a convex closed proper function f that satisfies o = prox with prox7 = arg ming 1 || z|1² + (z)). When o is ReLU, then = proxž for Va > 0 with being the indicator of the positive octant, i.e. f(x) = f{x 2 0}. The above monotone inclusion problem admits a unique solution if the operator F is strongly monotone, i.c. I I-G" -G @ W 2 ml
or,
2 I @ w W (1-m)I.
Therefore, we obtain the following well-posedness condition for MIGNN:
"See Appendix for review of some properties about the Kronecker product.
3

---

Under review as a conference paper at ICLR 2023
Proposition 1 (Well- posedness condition for MIGNN), Let the non-linearity a be ReLU and K (G @ W + G & ). Then the MIGNN model Equation (2) is well-posed as long as K I (1 I some >> 0. As K is symmetric, K N (1 m) is equivalent requiring that each eigenvalue of K is no more than 1 - m.
We provide the proof of Proposition [ the appendix; similarly, the proof of all the subsequent theoretical results are provided in the appendix. The well-posedness condition in Proposition [ allows for more flexible parametrizations than (39] by enabling the real part of eigenvalues of W to be in the range (-o0, 1) and the imaginary part to be arbitrary. Along with providing more flexible well-posedness condition for MIGNN, monotone operator theory guides us in designing efficient algorithms for implementing MIGNN; see Sec.
FLEXIBLE PARAMETERIZATION OF MIGNN
This section presents the monotone and orthogonal parameterizations of W for MIGNN in Equa- tion D. The monotone parameterization can enhance IGNN's expressivity, and the orthogonal parameterization can stabilize and accelerate the training of MIGNNs.
3.1 MONOTONE PARAMETERIZATION
Proposition iiffrrms us to design a more expressive parameterization of W for MIGNN than that used for IGNN leveraging monotone operator theory.
Proposition 2 (Monotone parameterization). Let G = (V, E)) be a graph and let G be L/2 with :: D--DD=/// (D2 A)D 1/2 being the normalized Laplacian, where A is the adjacency matrix and
D is the degree matrix with Dù, = شi=1 Aig. Then the MIGNN model ( o (W z G + 9B gB(X) (X)) is well-posed when the weight matrix W is parameterized as follows
W (1 m) CC F
(4)
where C, F E Rdxd are arbitrary matrices, and m > O.
Remark 1. In monotone parameterization, we first set the graph-related matrix G to be L/2, whose
eigenvalues are in [0, 1]. In contrast, the range of the eigenvalues of used in IGNN, see Sec. is
[-,,1]. Next, we parameterize W in Equation (4. whose eigenvalues have real part in (-00,
m]. Thus, { @ W + G & W 1 (1 m) guaranteeing the well-posedness of MIGNN. Moreover, W = (1-m) -C + F describes all possible W that satisfy W 3 (1-m)
3.2 ORTHOGONAL PARAMETERIZATION
As discussed in Sec. D IGNN learns LRD when X1 ( approaches in magnitude. This of- ten not the case when starting from Gaussian random initialization making IGNN unstable for learning LRD. Inspired by the unitary RNN (S), we propose to use the orthogonal parameteriza- with a learnable scaling factor to stabilize MIGNN in learning LRD. In particular, we parameterize W by the following scaled Cayley map
tion 41 54 53
((>) (I)))II))(+++
(5)
where (6-) is the sigmoid function and 7 E R is a learnable parameter ensuring <(7) E (0, 1). S T is a skew-symmetric matrix with C E dddd being an arbitrary parameterized matrix. is evident that MIGNN with the parameterization Equation ( is well-posed with G being Á defined in Sec.
[ Also, all eigenvalues of (I () + 5)-1 have magnitude 1, see derivation in To effectively learn LRD, MIGNN only requires the scalar Φ(>) to converge to 1.
Appendix
E.3
4 ACCELERATED OPERATOR SPLITTING FOR IMPLEMENTING IGNNS
It is worth noting that monotone and orthogonal parameterizations are beyond the efficient con- vergence regime of the Picard iteration. Thus, we leverage the operator splitting schemes to find the fixed point of the equilibrium equation with monotone orthogonal parameterization. Oper- ator splitting schemes often converge faster than Picard iteration and can guarantee convergence of IGNNs even when Picard iteration fails (68). In particular, for small graphs and tasks where learning LRD is not crucial, we use accereerrteeeeertteedeeeeeeeeedd forward-backe backward splitting (FB) to im- plement MIGNN with monotone parameterization. For tasks that require learning LRD, we employ
4

---

Under review as a conference paper at ICLR 2023
Anderson -cceeeeredd Peaceman- Rachford splitting (PR) with the Neumann series approximation accompanied by diffusion convolution, to implement MIGNN with orthogonal parameterization.
We structure this section follows: In Sec. finding the fixed point of MIGNNs using monotone/orthogoaal parametcrization. In Sec. present backward propagation algorithms for updating the parameters of MIGNN.
present FB (Sec. 4.1.1 MPR (Sec. 4.1.2 for
WC
4.1 FORWARD PROPAGATION FOR FINDING THE FIXED POINT
4.1.1
FB SPLITTING
We can find the fixed point of MIGNN in Equation 2 via FB splitting with iterative scheme
(k+1) (z = prox Z (k) - (z W th G (0=(())())
ar >> 0 is constant. (6)
F.I Note that the Lipschitz con- W [88 Section 5]. Therefore, FB splitting converges to the fixed point if a 2m/ 2m/1II- G @ WI By choosing proper a, FB splitting can converge in the regime that Picard iteration does not. However, when the mono- tone parameterization used || [W]| can be arbitrarily large. Thus a needs be small to guarantee the convergence of FB splitting, in which case the Lipschitz constant is close to 1, and the conver- gence of FB splitting will be significantly slowed. FB splitting appealing for learning with small graphs and tasks where learning LRD is not crucial. In this case, we use monotone parameterization to improve the expressivity of the model, and we denote the MIGNN with monotone parameteri- zation using FB splitting MIGNN-MA Mon. For large graphs and tasks that require learning LRD, FB splitting suffers from slow convergence. Next, we will present PR splitting, which is better for learning large-scale graphs and LRD. Furthermore, we argue that PR splitting is not suitable for implementing MIGNN with monotone parameterization.
We provide a detailed implementation FB splitting in Appendix
stant of the FB iteration is _ := v1 2am + a² I G
4.1.2 PR SPLITTING
PR splitting used in (77] is guaranteed to converge for much broader choice of a and requires fewer iterations than FB splitting. However, each iteration of PR splitting requires inverting large matrices, which is computationally much more expensive and less scalable than FB splitting. PR splitting finds the solution Z' of the MIGNN by letting 2' (T*) where U Raxn is obtained from the fixed-poin exint iteration vec(Ur((+1)) = FPR (Pee((((((())) = CrCo (veo(ve) with Cr and Cg being the Cayley operators (see Appendix B for details) of F and G, respectively. Let r (k) be the shorthand notation of (k)( Then we can formulate the PR splitting as follows
a (++1) ..PR (k) 2V (2 prox (u" (*)) (*) + a veo(ga(X))) - proxf (u (*)) + u (8),
(7)
where the matrix V- = (I aI GT @ W)) and u is the zero vector. With the parametriza- tions discussed in Sec. the linear operator F Equation
3 is strongly monotone and I-LLipschitz where L = 1|1 - G X W||. Therefore. its Cayley operator CJ and hence FPR is contractive with the optimal choice of a being 1/L, see (68] Section 6]. In particular,
suggested to choose 1/(1+ + +(>)) when using orthogonal parametrization W - =I 5) + S) The pseu- docode for the detailed implementation of PR splitting Equation 7 can be found Appendix F.I
Remark 2. Douglas-Rachjan (DR) splitting is another option for solving MIGNN, which is often faster than PR. However, in our case PR contractive, making faster than DR for the same a.
PR splitting also benefits MIGNNs in learning LRD when an orthogonal parameterization is used.
To see this, we have the following Neumann series expansion of (u)
vec (W"DOG)
G &W
V(u (I+a(II G @W) (( I+a
1+1/a
(1+1/a)
where the last equality follows from (ABB)= A B and (A&B)vec B)c vec for AA B and C' that satisfy dimensional consistency. Equation (8) indicates that each node can access information from its ho-hop neighbors in a single PR iteration for MIGNN with orthog- onal parameterization. This cannot be said of monotone parameterization with large [ as
Foor the sake of presentation, we denote Anderson-accellratte FB and PR splitting FB and PR.
5

---

Under review as a conference paper at ICLR 2023
the Neumann series expansion in the last equality of Equation 8 longer applies. Evaluating Stewart algorithm [1]. which con- verts computing V into diagonalizing the matrix G and W, respectively. From Equation 88. we have
Ita (I G¹ &w) "0 ((((*) can be carried out by using
1+1/a
V(vec(U(()) + vec @Qw H @ (Qw (6T )00
where QGT ACT Qot and Qw Qw are the eigen-deccomossiton of GT and of W. respec- tively, and H E Rdxx whose (i. entry is Hij = 1/(1 - (AG))(AAT))))) We provide a proof of Equation ② in Appendix E.44 According to Equation (9). one only needs to calculate the sigen-decomposition of G once prior to training and the cigen-decomposition of W once per epoch. The above matrix inversion procedure echos the idea of EIGNN [58]. MIGNN has multiple layers, with each fixed point iteration representing one layer. In contrast, EIGNN is reducible to a one-layer model; see Appendix A.2 for details EIGNN.
Although PR splitting can capture LRD in a single iteration, computing V in Equation 7 requires computationally prohibitive matrix inversion. We provide two remedies to address this issue for MIGNN using orthogonal parameterization: 1) We use Neumann series expansion to approximate the matrix inversion when orthogonal parameterization is used. 2) We replace the graph-related matrix G with a generalized graph diffusion convolution matrix, e.g. heat kernel or the personalized PageRank [34 B3 Notice that the above two remedies do not work for MIGNN using monotone parameterization since we can no longer use the Neumann series approximation. Therefore, MIGNN with monotone parameterization using PR splitting is not scalable learning large graphs
Neumann series approximation. In the orthogonal parameterization of W we have || ensuring efficient approximation of V in Equation using only few terms of its Neumann series expansion. The K-th order Neumann series expansion of V (vee(U is given by
G &w
|| < 1,
1+1/a
K
( G'
(v((ve(UU(*)))
(10)
(1+ 1/a)
1=0
According Equation 7 the K-th order Neumann series approximated PR iteration function, denoted as F.P.KK can be written as follows
u (k+11) = EPR,K (+) = 2Nx (2 prox (u ()))) u (4)
(*) + vec(ge(X))) - prox (u ) + u (*) (11)
Each node can access information from its K-hop neighbors using the K-th order Neumann series approximated PR iteration, which is more efficient than the existing IGNN. Also, such treatment can significantly accelerate forward propagation. We can intuitively understand this as follows: Each iteration of MIGNN, with K-th order Neumann series approximated PR iteration, aggregates information from K Khhop neighbors, enabling the use of much fewer iterations than that of IGNN, which aggregates one hop per iteration. MIGNN can much smaller A1 ([W) than IGNN to reach the same number of hops, meaning MIGNN converges much faster than IGNN.
MIGNN with diffusion convolution. We can also improve MIGNNs for learning LRD using graph diffusion convolution 34 1. i.e. instead of using Å or L defined in the previous context, we can set G be the combination of higher powers of Á or L, so that each node aggregates features from multi-hop neighbors at each iteration. In particular, any positive integer P, where D is the degree matrix with D = ΣØ Σ (A)) other choices of G can be found in 34 We can show that the eigenvalues of D 1/2 A +.... A are all within [-1, 1]; see E.4 for a proof. As such, the orthogonal parameterization of W still ensures the posl-pose of MIGNN. We write the MIGNN with P-th order diffusion matrix G as follows
let G D-1/2 + D 1/2 for
Z = (WD (A+ A +... A D 1/2 + 9B (X)).
(12)
We can further apply the operator splitting schemes to Equation 12) In particular, we denote the is implemented using P-th
model as MI MIGNN-KDD when W is orthogonal, and Equation
order diffusion and K-th order Neumann series approximated PR iteration.
Now discussion the time complexity of MIGNN-NKDP. The P-th order diffusion matrix only needs to be pre-computed once in preprocessing with time complexity O(Pl)Erl) where n is the number of nodes, and |rl denotes the number of non-zero entries in P.. In each epoch, the

---

Under review as a conference paper at ICLR 2023
parameter K in the K-th order Neumann series affects the training time complexity linearly as Epl) where M denotes the maximal number of iterations, and d is the feature dimension which is much smaller than the number of nodes.
4.1.3 ANDERSON ACCELERATION
We have already seen that the main steps in both FB and PR splitting schemes involve solving iterative equations, e.g. Equations (6) and 77. and we can utilize Anderson acceleration 4 to accelerate the convergence of these iterative equations. We provide the detailed formulation and pseudocode for Anderson-acceless accelerated operator aplitting-basee MIGNNs in Appendix F.3
4.2 BACKWARD PROPAGATION FOR UPDATING MIGNNS
We derive backpropagation for MIGNN based on implicit differentiation 35 7 26 Recall that the vectorized MIGNN vee(Z) a (GT @ Wvec( Z) + voe(gg(X))) has equilibrium point
vec (. J. For any loss function l and any parameter 0 (W or B), we have
ae_ de 00
dved Z-) I (c & ww)
do (G' @ Wvec ( + vecgaa (X))
(13)
00
vec(gn (X)). The values of the first and last term in Equation (13) can be found through automatic differentiation by running one more iteration in the forward pass. Note that the product of the first two terms remains the same for any 0. Hence one only needs compute it once in each backward pass. However, it can still be expensive to find (0t) (Ōvec( (I (GT & W))- Following [77 Theorem 2). the operator splitting methods can be used in the backward pass so that computing (I ( JGT < W)) can be converted into computing V = (I- & w)) which is already calculated in the forward pass; see Appendix F22 Similar to the forward propagation, the backpropagation can also benefit from Anderson acceleration using an iterative formulation, and we provide more details in Appendix F22
where is the Jacobian of a evaluated at GT @ Wvec (*)
5 EXPERIMENTAL RESULTS
In this section, we compare the performance of MIGNN-Mon (MIGNN with monotone parameter- ization implemented via FB splitting) and MIGNN-NKDP (MIGNN with orthogonal parameter- ization implemented via PR splitting accompanied by K-th order Neumann series approximation and P-th order graph diffusion convolution) with IGNN and several other popular GNNs on various graph classification tasks at both node and graph levels. We aim to show that 1) MIGNN-Mon is significantly more expressive than IGNN for both node and graph classifications, and 2) MIGNN- NKDP can learn LRD effectively, efficiently, and stably. The hyperparameters used in each model are provided Appendix K
We conduct all experiments using NVIDIA RTX 3090 graphics cards.
5.1 DIRECTED CHIIN CLASSIFICATION
100
Accuracy &
To show that MIGNNs can capture LRD the underlying graphs, we test them on the syn- thetic chain task using the experimental setup from [58] The chain task dataset comprises of e classes and n. single-linked directed chains, each containing/ nodes. For each chain, only the fea- ture on the first node encodes the label informa- tion. The data is partitioned into training, validation, and test sets of 5%, 10%, and 85%, respectively. We consider binary (c = 2) and three-class classification (c = 3) problems over several different chain lengths. For IGNN, we use the experimental settings used (7I.. For MIGNN, we consider MIGNN-NKDP for this task. Fig. 3 shows the averaged test accuracy over 5 random seeds of different models for classifying directed chains of length ranging from 50 to 300 in an increment of 50 for the binary case and from 40 to 200 in an increment of 20 for the three-class case. For binary classification, MIGNN-N3D3 and MIGNN-N3DS both score perfectly for all random initial- izations of the considered chain lengths. For the three-class task, both MIGNN models achieve high accuracy consistently with the higher order diffusion models, and the higher order diffusion model outperforms the lower order diffusion model on longer chains. In contrast, the accuracy of IGNN is much lower and less stable than that of MIGNNs, and in general, IGNN's performance becomes worse the chain length increases. We provide an ablation study of the impact of the order of
82
CAN
300
150
200
252
150 150
203
Chain Length
Chain Length
Figure 3: The accuracy IGNN and MIGNN dif-
ferent configurations for classifying directed chains of different lengths. Left: binary classification (c = 2). Right: three-class classification (c 3).
7

---

Under review as a conference paper at ICLR 2023
Neumann series approximation and graph diffusion convolution the chain classification accuracy
and computational time in Appendix and H respectively.
We can also set G to be the diffusion matrix in Equation
2 to enhance IGNN's capability in
learning LRD. E.g. we can equip IGNN with a diffusion matrix of order 5, and we denote the resulting model IGNN-D5. Fig. further contrasts the performance of diffusion enhanced models, and we observe that MIGNN is more consistent and more accurate as the chain length increases.
Based the operator split- ting theory, we expect that MIGNNs are more computa- tionally efficient than IGNNs when both models can accu- rately classify the chain nodes. Fig. | compares the accuracy and computational efficiency of MIGNN-N2D5 over IGNN for three-class chain classification. We see that MIGNN-N2DS sta- bly approaches perfect accuracy compared to IGNN, which abruptly changes around epoch 400. When both models accurately classify the chains, MIGNN-N2DS requires fewer iterations and less computational time than IGNN.
Time Elapved (s]
3 6
702
201
1053
116
620
305
10l0
22: 503 8:0
Epoch
Figure 4: The accuracy and efficiency of MIGNN-N2D5 over IGNN for
three class chains, of length 140, classification.
5.2 GRAPH NODE CLASSIFICATION
E0:2T
TT.00
00:05
GENI EL
8K.18
77.01
6D-37
In this subsection, we contrast MIGNN-Mon and MIGNN- NIDI with some existing GNNs for several benchmark graph node classification tasks, including Cora. Citeseer, and Pubmed; each dataset's statistics of nodes/edge/averrge shortest path length between nodes are 2485/5069/5.77,
OCNAGES M4
6-5
7.5. 35
TS.72
DONS pu
85.50
78.24
$7.00
EIGNS 100
85. 80
73:
07-92
M3GNN-N5D1
87.04
74.91
83.55
Table 1: Node classification mean ac-
2120/3679/9.31 and 19717/44324/6.44 respectively. We use the curacy (%) for 10-fold cross- vross-validation.
training procedure outlined in 22 and report the mean accuracy of 10-fold cross validation in Table T The MIGNN-Mon outperforms the implicit model benchmarks IGNN and EIGNN on all three tasks. We provide an ablation study of the impact of the order of Neumann series and graph diffusion convolution for graph node classification in Appendix G and respectively.
5.3 GRAPH CLASSIFICATION
15
Foii-1
Fols-5
In this subsection, we verify that MIGNN-Mon can be more ex- pressive than IGNN for graph classification since the cigenvalues of
10-
Fols-10
monotone parameterization are more flexible than IGNN. We con-
sider five hioinformatics-relatee graph classification benchmarks: MUTAG, PTC, COX2, PROTEINS, and NCII [81), and some de- tails of these datasets are provided in Appendix J The training is performed using 10-fold crossvvaldattin using the experimental setup of 771. The averaged test accuracy and standard deviation across the 10 folds are shown in Table 2 For both IGNN and MIGNN-Mon, we use the hyperparameters outlined in [71]. We present the results for both IGNN and MIGNN-Mon in Table Clearly, MIGNN-Mon outperforms IGNN on all tasks. To verify our theory, we report on the evolution of A ( for three of the ten folds of MUTAG in Fig. 5 For all of the folds ( exceeds one. Table 2 also reports the accuracy of MIGNN-N3DI against several baseline models where it performs better than IGNN and GIND on all tasks and achieves the best accuracy COX2 and PROTEINS tasks among all stud- ied models. These results show that learning LRD effectively is vital for classifying these graphs. We provide ablation study of the impact of the order of Neumann series and graph diffusion convolution on classification accuracy and computational time in Appendix G and H respectively.
200
400
600
300
1000
Epoch
Figure 5: A W) of MIGNN- Mon vs. Epoch on MUTAG.
5.4 LARGER SCALE GRAPH NODE CLASSIFICATION
We further show the advantages of MIGNN-NKDP over IGNN and other GNNs for a larger scale graph node classification task- Amazon co-purchasing dataset, which contains 334863 nodes, 925872 edges, and the diameter of the graph is 44 [22] We provide more details of the Amazon co-purchasing dataset in Appendix g As in (23], we train portions of the graph ranging from 5% to 9%, and test on sets representing 10% of the total graph. We then report both Macro-FI and Micro-F1 consistent with (7T). Fig. contrasts the computational cost of MIGNN-NIDY with
8

---

Under review as a conference paper at ICLR 2023
praplaatog E aodes
198/17.9
34025,5
407412
1135001
WE DA
04. 1...1 1.5
55 PUF.F 4.0
63.2
74.7 4...0.0
54.8 14.5.F.0 0.5
DOCKS 13
85.8
68.6
76.6
T4.4
88 4
5.6
715.2 2:222:00
KANN
76.0
13.4
60.8 0-4
79-7
70.5
3.4
73.5
1.9
0500
92.3 + T.6
68.2
7.2
70.0
5.0
83.55
SIN LM
00.50-
3.3
2-2
SER NRA0
U.1
72.8
17:3
MICON MIDI
86.1
8.1
70.0 184460 6.0
85
2-8
79.0
3-3
75.4
Table 2:
Graph classification mean accuracy (%) = standard deviation for 10-fold cross- voossvvali, We take the results of the baseline models from (22 which are consistent with our reproduced results.
IGNN using 5% of the graph for training. ( of MIGNN-NIDI is much smaller than that of IGNN, implying faster convergence of MIGNN-NIDI than IGNN confirmed by the fact that MIGNN-NIDI saves significantly on the number of iterations and computational time over IGNN.
Time Eispere (s)
0.90
8 Iterations
300
MIGHN-NID1
0.80
0.75
25000 Epoch
Epoch
Epoch
Figure 6: Epoch vs. Az (WD) }}, the time required for each epoch, and iterations required for cach epoch of IGNN and MIGNN-NIDI for the Amazon dataset with 5% training portion.
Fig. contrasts MIGNN-NIDI with baseline mod- els when trained on portions of the graph ranging from 5% to 9%. We see that MIGNN-NIDI out-
Mou-F (N6
performs almost all baseline models over all differ- ent portions of the graph for the training. Though MIGNN-NII NID1 does not outperform IGNN signifi- cantly, MIGNN-NIDI enjoys significant computa- tional advantages over IGNN.
010
170
73.9
Fractoo
Fraction
Figure 7: Fraction vs. Micro-FI (left) and Macro- Fl (right) training accuracy on the Amazon dataset.
5.5 PHYSICAL DIFFUSION IN NETWORKS
0.015
We further consider physical problem of fluid flow in porous media, following [3) The model is a 3D graph whose nodes and edges correspond to pore chambers and throats. We sample train- ing graphs of different sizes between 100 and 500, which are gen- erated to fit into 0.1 m2 cubes. We aim to predict the equilibrium pressures Z inside pore networks G. We train MIGNN to min- imize the mean-squared error (MSE) between the prediction and Z" We utilize the experimental setup of (63 and include their reported results for IGNN, Both IGNN and MIGNN use the same encoder and decoder architecture. Graphs of 50- 200 nodes are sampled in training and 1000 test graphs are generated for pore counts from 200 to 500. Fig. X shows the MSE for the test graphs as the number of nodes (pores) varies from 200 to 500. MIGNN with both monotone and orthogonal parameterizations outperform IGNN by significant margin. For this task of learning physical diffusion in networks, CGS 63 performs better than MIGNN and IGNN in accuracy. As a future direction, we plan to integrate the idea of the learnable graph-related matrix G that is used in CGS with our proposed MIGNN to further improve the performance of MIGNN for learning physical diffusion in networks.
MSB
0.0010
0.005
0.000
700
400
300
Number of Pores
Figure 8: The average MSE of 500 sampled test iterations vs. the num- ber of pores. The error bars repre- sent the standard error of the pre- diction. MIGNN with different pa- rameterizations outperforms IGNN by a significant amount.
CONCLUDING REMARKS
We propose MIGNN based on monotone operator viewpoint of IGNN. In particular, MIGNN can be parameterized more flexibly than the benchmark IGNN. We provide efficient implementations of MIGNN that integrates diffusion convolution leveraging different operator splitting schemes with Anderson acceleration. Numerically, MIGNN remarkably outperforms the existing IGNN in accu- racy, stability, computational efficiency, and learning LRD. As IGNNs are closely related RNNs, an interesting future direction is to explore if the ideas from other RNN architectures 42 59 can be adapted to the improvement of IGNNs.

---

Under review as a conference paper at ICLR 2023
REFERENCES
[1] Sami Abu-El-Haipp Amol Kapoor, Bryan Perozzi, and Joonseok Lee. N-gen: Multi-scale graph convolution for semi-supervised node classification. In uncertainty artificial intelli- gence, pp. 841-851. PMLR, 2020.
[2] Luis B. Almeida. A learning rule for asynchronous perceptrons with feedback in combina- torial environment. Artificial neural networks: concept learning, pp. 102-111, 1990.
[3] Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practi- cal implications. In International Conference on Learning Representations, 2021. URL https: htttps://pperreveww. ftt.:eee/orrru2211111100000000
[4] Donald G. Anderson. Iterative procedures for nonlinear integral equations. Journal of the ACM (JACM), 12(4):547-560, 1965.
[5] Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural net- works. In International Conference on Machine Learning, pp. 1120-1128, 2016.
[6] James Atwood and Don Towsley, Diffussion-convllutional neural networks. In Advances in Neural Information Processing Systems, volume 29, 2016.
[7] Shaojie Bai, J. Zico Kolter, and Vladlen Koltun. Deep equilibrium models. Advances in Neural Information Processing Systems, 32, 2019.
[8] Shaojie Bai, Vladlen Koltun, and J. Zico Kolter. Multiscale deep equilibrium models. In Proceedings of the 34th International Conference on Neural Information Processing Systems, 2020.
[9] Shaojie Bai, Vladlen Koltun, and J. Zico Kolter. Stabilizing equilibrium models by Jacobian regularization. International Conference on Machine Learning, pp. 554-565. PMLR, 2021.
[10] Shaojie Bai, Vladlen Koltun, and J Zico Kolter. Neural deep equilibrium solvers. In Inter- national Conference on Learning Representations, 2022. URL httpps:/
[11] Richard H. Bartels and George W. Stewart. Solution of the matrix equation ax+ xb= [f4]. Communications of the ACM, 15(9):820-83) 826, 1972.
[12] Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE Transactions on Neural Networks, 5(2): 1994.
[13] Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares- López, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentia- tion. arXiv preprint arXiv:2105. 2021.
[14] Cristian Bodnar, Fabrizio Frasca, Nina Otter, Yuguang Wang. Pictro Liò, Guido F Montufar, and Michael Bronstein. Weisfeiler and Lehman go cellular: CW networks. Advances in Neural Information Processing Systems, 34:2625 2021.
[15] Cristian Bodnar, Fabrizio Frasca, Yuguang Wang, Nina Otter, Guido F Montufar, Pietro Liò, and Michael Bronstein. Weisfeiler and Lehman go topological: Message passing simplicial networks. International Conference on Machine Learning, pp. 1026-1037. PMLR, 2021.
[16] Giorgos Bouritsas, Fabrizio Frasca, Stefanos P Zafeiriou, and Michael Bronstein. Improving graph neural network expressivity via subgraph isomorphism counting. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022.
[17] Joan Bruna, Wojeiech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and deep locally connected networks on graphs. In International Conference on Learning Representa- lions, 2014.
10

---

Under review as a conference paper at ICLR 2023
[18] Ben Chamberlain, James Rowbottom, Maria I. Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. GRAND: Graph neural diffusion. In Proceedings of the 38th International Conference Machine Learning. volume 139, pp. 1407-1418. PMLR, 18 24 Jul 2021. URL
m https:://procees: wtt.aea..sssss.sssee/e/eo. v 1 39/ chw.berabbbbeeee a a htmmlll
[19] Benjamin Paul Chamberlain, James Rowbottom, Davide Eynard, Francesco Di Giovanni, Xi- aowen Dong, and Michael M. Bronstein. Beltrami flow and neural diffusion graphs. In Advances in Neural Information Processing Systems, 2021.
[20] Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun. Measuring and relieving the over- smothhin problem for graph neural networks from the topological view. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34 (04), pp. 3438-3445, 2020.
[21] Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. Simple and deep graph convolutional networks. In International Conference on Machine Learning, pp. 1725-1735
PMLR, 2020.
[22] Qi Chen, Yifei Wang, Yisen Wang. Jiansheng Yang, and Zhouchen Lin. Optimization-inducee graph implicit nonlinear diffusion. In International Conference on Machine Learning, pp- 3648 1668.666. PMLR, 2022.
[23] Hanjun Dai, Zornitsa Kozareva, Bo Dai, Alex Smola, and Le Song. Learning steady-states of iterative algorithms over graphs. In International conference on machine learning, pp. 106- 1114. PMLR, 2018
[24] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. Advances in neural information processing systems, 29, 2016.
[25] David Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alan Aspuru- Guzik, and Ryan P. Adams. Convolutional networks on graphs for learning molecular fingerprints. Advances Neural Information Processing Systems, volume 28, 2015.
[26] Laurent EI Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Tsai. Implicit deep learning. SIAM Journal on Mathematics of Data Science, 33)-9930-958, 2021.
[27] Moshe Eliasof, Eldad Haber, and Eran Treister. PDE-GCN: Novel architectures for graph neural networks motivated by partial differential equations. Advances in Neural Information Processing Systems, 34:3836-3848 2021.
[28] Moshe Eliasof, Eldad Haber, and Eran Treister. pathGCN: operators from paths. In International Conference Machine Learning, pp. PMLR, 2022.
Learning general graph spatial
[29] Jeffrey L.. Elman. Finding structure in time. Cognitive Science, 14(2):179-211, 1990.
[30] Claudio Gallicchio and Alessio Micheli. Graph echo state networks. In The 2010 international joint conference on neural networks (IJCNN), pp. 1-8. IEEE, 2010.
[31] Claudio Gallicchio and Alessio Micheli. Fast and deep graph neural networks. In Proceedings of the AAAI conference on artificial intelligence, volume 34 (04), pp. 3898-3905 2020.
[32] Johannes Gasteiger, Aleksandar Bojchevski, and Stephan Günnemann. Predict then propagate: Graph neural networks meet personalized pagerank. In International Conference on Learning Representations, 2018.
[33] Johannes Gasteiger, Aleksandar Bojchevski, and Stephan Günnemann. Combining neural networks with personalized pagerank for classification on graphs. In International Confer- Learning Representations, 2019. URL htttps://opereete
net/forum
(d=H1gL-2A9Yy
[34] Johannes Gasteiger, Stefan Weißenberger, and Stephan Günnemann. Diffusson improves graph learning. In Advances in Neural Information Processing Systems, volume 32, 2019.

---

Under review as a conference paper at ICLR 2023
[35] Jean C. Gilbert. Automatic differentiation and iterative processes. Optinization methods and software, 1(1):13-2,, 1992.
[36] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 1263- 1272. PMLR, 06-10 Aug 2017. URL
https::/prroceddi ngs m press vv00///00/////// m 7a html
[37] Shunwang Gong, Mehdi Bahri, Michael M Bronstein, and Stefanos Zafeiriou. Geometrically principled connections in graph neural networks. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 1 11415-11 1424, 2020.
[38] Marco Gori, Gabriele Monfardini, and Franco Scarselli. A new model for learning in graph domains. In Proceedings. 2005 IEEE international joint conference on neural networks, vol- ume 2, pp. 729-734, 2005.
[39] Fangda Gu, Heng Chang, Wenwu Zhu, Somayeh Sojoudi, and Laurent El Ghaoui Implicit graph neural networks. In Proceedings of the 34th International Conference on Neural Infor- mation Processing Systems, 2020.
[40] Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. In Advances in Neural Information Processing Systems, volume 30, 2017.
[41] Kyle Helfrich, Devin Willmott, and Qiang Ye. Orthogonal recurrent neural networks with scaled cayley transform. In International Conference on Machine Learning, pp. 1969-1978. PMLR, 2018
[42] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, (8):1735-1780, 1997.
[43] Roger A. Horn and Charles R. Johnson. Topics in matrix analysis, 1991. Cambridge University Presss, Cambridge, 37:39, 1991.
[44] Saber Jafarpour, Alexander Davydov, Anton Proskurnikov, and Francesco Bullo. Robust Im- plicit Networks via Non- Kon-Fuclidean Contractions. In Advances Neural Information Process- ing Systems, volume 34, pp. 9857-9868 2021. URL
igs neurips. ineuri cccpaper/2021 http://ee2222a112222001111011110 html
[45] Li Jing, Yichen Shen, Tena Dubcek, John Peurifoy, Scott Skirlo, Yann LeCun, Max Tegmark, and Marin Soljačić. Tunable efficient unitary neural networks (cunn) and their application to mns. In Proceedings of the 34th International Conference on Machine Learning-Vait Volume 70, pp. 1733-1741. JMLR. org, 2017.
[46] Kenji Kawaguchi. On the theory of implicit deep learning: Global convergence with implicit layers. International Conference on Learning Representations, 2021. URL https:// // openreview net. forum? -ZZ
[47] David Kincaid and Ward Cheney. Numerical analysis, brooks. Cole Publishing Company, 20: 10-13, 1991.
[48] Thomas N. Kipf and Max Welling, Semi-supervised classification with graph convolutional networks. In Proceedings of the 5th International Conference on Learning Representatiinss 2017. URL httpttttt//:prreeeeee
n www.eeteeoommmmda--aaaaaaaa.
[49]
Zico Kolter and Gaurav Manek. Learning stable deep dynamics mod- els. In Advances in Neural Information Processing Systems, volume 32, 2019.
httpppps.....pppe/////1111 0a 4bbceda 1 17 226253386bc 232233338cccc65552222222--- www.e.ddara..a.aaap.ppa..ppp.
www.nnnuii.... URL https: http::://prcceedingg wtw.:eewwdddin..eeeiiieeseeeeeee
[50] J. Zico Kolter, David Duvenaud, and Matt Johnson. Deep implicit layers neural odes, deep equilibirum models, and beyond. http::/immlicitt-ayeeee huw.:iit---eyyeeetttttoa.l...... 2020.
tutorial.
12

---

Under review as a conference paper at ICLR 2023
[51] Steven George Krantz and Harold R. Parks. The funciiittunction applications. Springer Science & Business Media, 2002.
function theorem: history, theory, and
[52] Jure Leskovec, Lada A. Adamic, and Bernardo A Huberman. The dynamics of viral marketing. ACM Transactions on the Web (TWEB), 1(1):-ess, 2007.
[53] Mario Lezcano Casado. Trivializations for gradient- based optimization on manifolds. Ad- vances in Neural Information Processing Systems, 32, 2019.
[54] Mario Lezeano-Casado and David Martunez-Rubic Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group. In International Conference on Machine Learning, pp. 3794-3803, PMLR, 2019.
[55] Guohao Li, Matthias Müller, Bernard Ghanem, and Vladlen Koltun. Training graph neural networks with 1000 layers. In International conference on machine learning, pp. 6437-6449. PMLR, 2021.
[56] Qimai Li, Zhichao Han, and Xiao-Ming Wu. Deeper insights into graph convolutional net- works for semi-supervised learning. In Thirty-Second AAAI conference on artificial intelli- gence, 2018.
[57] Yujia Li, Richard Zemel, Mare Brockschmidt, and Daniel Tarlow. Gated graph sequence neural networks. In Proceedings of ICLR' 16, 2016.
[58) Juncheng Liu, Kenji Kawaguchi, Bryan Hooi, Yiwei Wang, and Xiaokui Xiao. Eignn: Efficient dethhite-depth graph neural networks. In Advances in Neural Information Processing Systems, pp. 18762-18773, 2021.
[59] Denis Lukovnikov and Asja Fischer. Improving breadth- backpropagation in graph neural networks helps learning long-range dependencies. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pp. 7180-7191, PMLR, 18-24 Jul 2021. URL httpps://proceeddngg.s ml press/ 39/ wukovni kov html
[60] Zheng Ma, Junyu Xuan, Yu Guang Wang. Ming Li, and Pietro Liò. Path integral based con- volution and pooling for graph neural networks. Advances in Neural Information Processing Systems, 33:16421-16633 2020.
[61] Zakaria Mhammedi, Andrew Hellicar, Ashfaqur Rahman, and James Bailey. Efficient orthogo- nal parametrisation of recurrent neural networks using householder reflections. In Proceedings of the 34th International Conference on Machine Learning- Volume 70, pp. 2401-2409. JMLR.
org, 2017.
[62] Kenta Oono and Taiji Suzuki. Graph neural networks exponentially lose expressive power for node classification. In International Conference on Learning Representations, 2020. URL
net, forum? https:// wtw.:ewuueeetteeteeeeeeea
[63] Junyoung Park, Jinhyun Choo, and Jinkyoo Park. Convergent graph solvers. In Interna- tional Conference on Learning Representations, 2022. URL htt.pps://pperrevie..
net/
[64] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In International Conference on Machine Learning, pp. 1310-1318, 2013.
[65] Honghin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang. Yu Lei, and Bo Yang. Geom-GCN: Geometric graph convolutional networks. In International Conference on Learning Represen- tations, 2020.
wtwpr/eewwwwwweeeeeeeeeeeeeeeeee net wntptt//ooouummm111eeeegg..g.e.
[66] Fernando Pineda. Generalization of back propagation to recurrent and higher order neural networks. In Neural information processing systems, 1987.
[67] T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-couple oscillator networks. In International Conference on Machine Leamn ing, pp. 18888-18909, PMLR, 2022.
13

---

Under review as a conference paper at ICLR 2023
[68] Ernest K. Ryu and Stephen Boyd. Primer monotone operator methods. Appl. Comput. Math, 15(1):3 43, 2016.
[69] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfar- dini. The graph neural network model. IEEE transactions on neural networks, 20(1):61-800 2008.
[70] Nino Shervashidze, SVN Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten Borgwardt. Efficient graphlet kernels for large graph comparison. In Artificial intelligence and statistics, pp. 488-495, PMLR, 2009.
[71] SwiftieH. Implicit graph neural networks. https://github.coe http://ittuub..cooo........t http:///wwiiee//tttta 2020.
[72] Matthew Thorpe, Tan Minh Nguyen, Hedi Xia, Thomas Strohmer, Andrea Bertozzi, Stanley Osher, and Bao Wang. GRAND++: Graph neural diffusion with source term. In Inter- national Conference on Learning Representations, 2022. URL
het:/forumlld-EMMuudddee
[73] Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. In International Conference Learning Rep- resentations, 2018. URL
net wtw.:e/wwwweeee/oooodid--dd----e IXMMpicc2
[74] Eugene Vorontsov, Chiheb Trabelsi, Samuel Kadoury, and Chris Pal. On orthogonality and learning recurrent networks with long term dependencies. In Proceedings of the 34th Interna- tional Conference on Machine Coarning-Volume 70, pp. 3570-3578. JMLR. org, 2017.
[75] Max Welling and Thomas N. Kipf. Semi-supervised classification with graph convolutional networks. International Conference on Learning Representations, 2016.
[76] Paul J. Werbos. Generalization of backpropagation with application to a recurrent gas market model. Neural networks, 1(4):339-333 -363, 1988.
[77] Ezra Winston and J. Zico Kolter. Monotone operator equilibrium networks. In Advances in neural information processing systems, volume 33, pp. 10718-10728, 2020.
[78] Scott Wisdom, Thomas Powers, John Hershey, Jonathan Le Roux, and Les Atlas. Full-capacity unitary recurrent neural networks. Advances in Neural Information Processing Systems, pp- 4880 4888, 2016.
[79] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger. Simplifying graph convolutional networks. In International conference on machine learning, pp. 6861-6871 PMLR, 2019.
[80] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful graph neural networks? In International Conference on Learning Representations, 2019. URL https: //openreview. net hett/oormm22daryys66533
[81] Pinar Yanardag and V. Vishwanathan. Deep graph kernels. Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp- 1365-1374, 2015.
[82] Jaewon Yang and Jure Leskovec. Defining and evaluating network communities based on ground-truth. In Proceedings of the ACM SIGKDD Workshop on Mining Data Semantics, pp- 1-8, 2012.
[83] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen. An end-to-end deep learn- ing architecture for graph classification. Proceedings of the AAAI conference on artificial intelligence, volume 32 (1), 2018.
[84] Lingxiao Zhao and Leman Akoglu. Paimorm: Tackling oversmoothing in gnns. In Inter- national Conference on Learning Representations, 2020. URL het//forum?LL=:xxee111ttt
14

---

Under review as a conference paper at ICLR 2023
A A BRIEF REVIEW OF IGNN AND RELATED MODELS
A.1 IGNN: FORWARD AND BACKWARD PROPAGATION
IGNN employs a projected gradient descent method in the training phase to ensure their pro- posed pell-posedness condition is satisfied. In forward propagation, IGNN finds the equilibrium through direct Picard iteration. During backward propagation, IGNN uses the implicit function the- orem the equilibrium compute the gradient. The computationally expensive terms related to at @vec(2")
(I J J (GT @ W))
(see Section 4.2 for notations) is also computed implicitly through
Picard iteration.
A.2 EIGNN, CGS, AND GIND
EIGNN Efficient infinite-depth graph neural networks (EIGNN) is implicit graph neural net- work model proposed by Liu et [58] whose counterpart explicit GNN is simple graph convo- lution (SGC) [79]. The main update step in EIGNN is given by
((++)) (F)
(14)
where z denotes the hidden feature, G is the normalized augmented adjacency matrix A (See Section [. X is the input feature, g is the weight matrix which is parameterized to guar- antee convergence, and y is constant scalar in (0, 1). Note that, there is non-linearity in the fixed-point Equation 14 and this allows EIGNN to find the equilibrium by the following closed
formula:
lim vec (zttt (I GG @ ((()))) vx)(X).
(15)
k-000
For computation efficiency consideration, the matrix inverse operation is reduced eigenvalue de- composition of G and g(F) where the eigenvalue decomposition G is pre-calculated before training.
CGS
Convergent graph solver (CGS) is an implicit graph neural network proposed by Park et al. where the fixed point equation in use be described as follows
in 663
(b++)) >Z G (X)
(16)
where ( is the hidden feature, y is the contraction factor, G9 E Fax= is the graph-related matrix that is learnable and gB( (x(X) is the input-dependent bias term Similar to the EIGNN case, the linearity in Equation 16 16 the fixed point be found by closed formula.
GIND
The optimization-inducaa graph implicit nonlinear diffusion (GIND) is implicit graph neural network proposed by Chen et al. [22]. GIND involves a fixed point iteration equation of the following form:
Z (k+1) T (Z + G
(17)
where Z-) is the hidden feature, W is the weight matrix, 9B (X) is some input-dependent bias term, and G is a normalization of the adjacency matrix A. The precise definition of G is given 1/2
as G = D 1-2 1///2///2 where D is the degree matrix of the augmented adjacency matrix A + I given as D.
+ E, Aij. The weight matrix W is parameterized so that ||| < 1. Similar to IGNN, the Picard iteration is employed to find the fixed point. The authors claimed that the new fixed-point equation (Equation 17) represents a nonlinear diffusion process with anisotropic properties while IGNN only represents a Tinear isotropic diffusion. However, we observe that GIND is closely related to the following simple variant of IGNN where the main change is to
z (k+1) = (w(-TT)z (K)
G G (
(18)
where the notations are the same as in Equation 17 In fact, once W [|] < and assuming is non-expansive activation function (for example, tanh, ReLU, ELU), then Equation 18 contractive and hence its fixed point exists. Let z" be the fixed-point of Equation (18), then claim that
15

---

Under review as a conference paper at ICLR 2023
- is the fixed point of Equation 117 with the same W, G, and (œ(X) used in
both Equation 18 and Equation 17 This can be seen from the following direct calculation:
2 = =-TGGT
WgB(XiG) = --W'OWW-WW'Z''GG G+W G
=-W (W2G+ Wyn(XiG) G¹
= (W(2+g((X))GG)
B BRIEF REVIEW OF MONOTONE OPERATOR THEORY
B. OPERATORS
In this section, we briefly review the definition and basic theory of monotone operators, more details can be found in [68] We sat T is (set- (set-valued) operator if maps a point in IR to subset of IR and we denote this as T IR Rd We define the graph of an operator as
Gra T {(2z u) T(x)).
Mathematically, operator and its graph are equivalent. In other words, we can view T Rd + Rd as point-to-set mapping and as subset of Rd Rd.
Many notions for functions can be extended to operators. For example, the domain and range of an operator T are defined as
dom T = {a: T(x)) + 0}, range (y lulyy y= E R²)..
If T and S are two operators, we define their composition as
TS(m) T(S(a)),
and their sum as
(T+ = -(xx) + S(x).
Alternately, we can define the operator composition and sum using their graphs,
TS {x x) 3 (z.z) S, (y, T)
T+S= {(x, y + 2) (x. y) T, (z, 8). The identity (Z) and zero (0) operators are defined follows
I ((x, x)) Rd We say an operator T is L-Lipschitz (L > 0) if
((x,0)|xE x (
[[[(x) T(x) < L| yll, Vo, y E dom T,
i.c.
uu v|| s T.|l= vll. (g,) ET. The inverse operator of T is defined
T- = ( (I. (ee.))77 ET.
When 0 T(x), we say that æ is zero T. We write the zero of an operator 7
Zer T {x 0 T(x)) (0).
16

---

Under review as a conference paper at ICLR 2023
B.2 MONOTONE OPERATORS
An operator T on Rd is said to be monotone if
uu U, a y) 2 0, u(,, (U, ET,
where (.)) denotes the inner product between two vectors. Equivalently, we can express mono- tonicity as
(T T - 20 y
Furthermore, we say the operator T is maximal monotone if there is other monotone operator other words, if the monotone operator T is not maximal, then there exists (æ, u) & T s.t. T {(xx u)} is still monotone. A continuous monotone function F IRd + Rd is maximal monotone.
s.t. Gra C Gra properly.
An operator T I2d R is B- B-strongly monotone or B-coercive if B > 0 and
(u U y) > B| y V(x, u) (y. v) T.
We say T is strongly monotone if it is B-strongly monotone for some unspecified constant B (0, o.). In particular, a linear operator F(x) Gæ + h for G E Rdddd and E R4 is maximal monotone if and only if G
(0 stands for the matrix whose entries all zero) and B-strongly monotone if [ (G G ž BI. Similarly, subdifferentiable operator å is maximal monotone if and only if a convex closed proper (CCP) function.
G
An operator T is B-cocoercive or B-inverse strongly monotone if B >> 0 and
(u v, - y) 2
v||3, V(x, u) (y. v) ET.
We say T is cocoercive if it is B-cocoercive for some unspecified constant 3 c0). In particular, if the linear operator F(x) = Gx + h is B-strongly monotone and L-Lipschitz, then F is B cocoercive.
F
C A BRIEF REVIEW OF OPERATOR SPLITTING SCHEMES
In this section, we provide brief review of few celebrated operator splitting schemes for solving ixed-point equilibrium equations.
C.I RESOLVENT AND CAYLEY OPERATORS
The resolvent and Cayley operators of an operator T is defined as, respectively, as follows
RT = (I (I+
and
= 2R
where a > is a constant. The resolvent and Cayley operators are both non-expansive, i.e. they both
have Lipschitz constant L < 1 for any maximal monotone operator T, and the resolvent operator Ry is contractive (i.e. L < 1) for strongly monotone T, the Cayley operator Cy is contractive for strongly monotone and Lipschitz
There are two well-known properties associated with the resolvent operators:
First, when F(xr) = Ga: + h. is linear operator, then
R,((x) (I + aG) (x (ah).
Second, when == 0f for some CCP function f. then the resolvent is given by the follow- ing proximal operator
Ry(xx) prox?(x)) == arg min { -=11= +af(s)}
17

---

Under review as a conference paper at ICLR 2023
C.2 OPERATOR SPLITTING SCHEMES
Operator splitting schemes refer to methods to find zero in a sum of operators (assumed here to be maximal monotone), i.e. find æ s.t.
0 (F + ()))
We present a few popular operator splitting schemes for solving the above monotone inclusion prob- lem.
Forward-backward splitting (FB): Consider the monotone inclusion problem
find 0 E (F + (a))
where F and G are maximal monotone and F is single-valued. Then for any 0, have
0 (F) e (I+ -I--a))aa + (I+ a0) -I- (z) x == -
Therefore, æ is a solution and only if it is a fixed point of Rg aF). Moreover, assume FF is fi-occerrcive, then the Picard iteration using barkudbb splitting can written as
x (k+1) = (x: oFF(()))
which converges a E (0, 23) and Zer(F + 9) +
Peaceman- Rachford splitring (PR): Consider the following monotone inclusion problem
find E (F + (2));
where F and G are maximal monotone. For any > 0, we have
OE (F +)) + az) (I (z)
e 0 6 (Z+))) - Cg( (g)
(Z) -C((2), = (I (
(Z == =o((z)
RxCg(=)=RRg(e) R=(==, =Rg(=)
CxCg(z) =, a = Rg(=).
Therefore, is solution if and only if there is solution of the fixed-point equilibrium equation z = CxCg(z) and æ = Rs(z), which is called Peacemun-Rachforl splitting.
Kouglas-Rachfon splitting (DR): Sometimes the operator CyCg is merely nonexpansive, the Picard iteration with PR given below
=
is not guaranteed to converge. To guarantee convergence, note that for any Va >> 0, have
(F+ GG + (x))) Jg(2).
And the above splitting is called Rouglas-Ruchford wtt.hhhhhhhospppisiiin..... The Picard iteration with DR can be written as follows
x 1++11/ = (1)
x. (k+1) RR (2x (6+1/22) - (b))
(k+1) (A) + x (k+l) 1k+1/2 1/2)
which converges for any if (F 9) + $.
18

---

Under review as a conference paper at ICLR 2023
D PROPERTIES OF KRONECKER PRODUCT
In this section, we collect some Kronecker product results that are used in this paper.
Definition 1. Let A R B R be two matrices. Their Kronecker product AxB E Rxxp defined as follows:
An
B
A B
Ap
Apg B
The following identities about Kronecker product hold:
(A Ø B A & BT AA E B E E
A Ø B = ||||||| B| YAE RPX; B Rrxe
A @ B B x
||||||| & YAER RPXq BE R
(A & B)) = vee (BCA AA E Rd
B E RPXg C Rixr
(A & & A
AA E R B RPxg g C R
A
A C A Rox9 B.CERT
(A+ 3 C - A B C VA, R
(A & (C 8 D) = =C= BD V E R C Raxk D R
Proposition 3 (143) Theorem 4.2.12)). Let A R and B E Rxxm If denote the eigenvalue
(A)) and A(B) = (B) Am (B)) }. then the eigenvalue set of A B is A(A B) (A)(A) ( 1l n....,n,) , m].
sets of A and B as AAA = (A)
E TECHNICAL PROOFS
E.I LIPSCHITZ CONSTANT VS. LARGEST MAGNITUDE OF EIGENVALUE
Let f(Z) = WZG+B be a linear map. With slight abuse notation, we still denote the vectorized version of f as f which reads f(vec(Z) properties of the Kronecker product). The Lipschitz constant Lip (f) of the linear map with
(G' W)ve Z)) vec( (B) (See Appendix D for
respect to the loo vector norm is exactly the co-norm G & I || | W W Recall the
following general result about matrix norm and the largest magnitude of eigenvalue.
Theorem 1 (14 Theorem 4 in Section 4.6)). The largest magnitude of eigenvalue (A) of matrix A satisfies
(A)
inf ||
lll
in which the infimum is taken over all subordinate matrix norms || || including 2-norm and oo- norm.
Meanwhile, note that one has [ flx [ by definition. Hence one has Lip (f) ||G || [|| bo Z AI Aa (W). Note that, when G the normalized adjacency matrix of (G))= = and hence we have (f) 2 A (1)).
undirected graph Á, we have (G (G)
E.2 PROOFS FOR SECTION 2
Proof of Proposition
First recall the operator splitting problem 3 in Section find (F + G))e (2)),
where
& W)veo - vec(yn(X)) and G of, here f is the indicator of the positive octant, i.c. f(x) = 1(x 2 0} for which we have prox equals a, the ReL.U activation function, for all 0. Note that, from the condition K ] (GT @ W
F(veo ( (I-
G & T N (1 I. one has G X W = (1 m and hence
G W mI
19

---

Under review as a conference paper at ICLR 2023
which says F is m-strongly monotone for some m 0. As the function F is a linear and hence continuous function defined on the entire Rdx", it is then automatically maximal monotone once it is monotone. Since is CCP function, its subdifferential operator ς 0f is maximal mono- tone. In particular, the linear map F is ingle-valued, can apply the FB splitting scheme in Appendix C.2 as the following: for any a > O, we have
E (F -vvccv vec(2)
RR (I aF)c()c())).
+ vec Z) = proxř (vec(Z) - a (vec G @ Wvec vec(gg(Xi))
vec Z (vec( Z) a veel Z) -G
W (Z) vec(gs(X))))
When a = 1 in the last above, we recover the MIGNN model
2
vec Z) = oG @ Wvec + vec(gB(X))
This shows the equivalence between finding a fixed point of MIGNN model 2 and finding a zero of the operator splitting problem Therefore, when K 3 (1 - I. the linear map F is strongly monotone and Lipschitz, the monotone splitting problem and hence the MIGNN model is well- sposed, see Appendix C.2
E.3 PROOFS FOR SECTION
The following properties of the Cayley map used in this paper.
Proposition 4. Let S be a skew-symmeetic matrix. Then its image under the Cayley map Cay(S) (I - SI + S) is an orthogonal matrix, and hence the magnitude of all its eigenvalues
Proof. To verify that the Cayley map is well-defined, it suffices to show that - is not an eigenvalue of S. This can be derived from the general fact that each eigenvalue of any skew-symmetric matrix is purely imaginary. To see this, let À be an eigenvalue of S with corresponding eigenvector v where both À and possibly contain complex numbers. Let v and S denote the conjugate transpose of the vector v and the matrix S respectively. We then have
SSS Sv o" (Ae) A/vle-
where l/c denotes the Euclidean norm for complex vector. At the same time, one has
SS O) H -xp12
where Ā denotes the complex conjugate of X. Hence ->, that is À is purely imaginary. This concludes the proof that (I SI + S) -1 is well-defined.
Note that (I - SI (1- (I-))++)) Therefore, (I - 5) S) is (real) orthogonal.
S)
=
In the last part, we present a short proof that the magnitude of all eigenvalues of a (real) orthogonal matrix O equals Let Ao be an eigenvalue O and ur is its eigenvector. Then we have
Doo |w ( = (Ow) H (Ow) = w HHO = (Ow) (on))((Ou))) (Ow) = T Ow = ww
Hence,
!ool 1.
Proof of Proposition
Since the normalized Laplacian is symmetric, we have
K = L @ W L ( - 1 L 0 (
2 (2
The property Kronecker product (Theorem 3) tells us that the eigenvalues K are the products of the eigenvalues of L and ( (W + W )) Therefore, the MIGNN model satisfies the well- posedness condition Proposition Tonce
A
2 L
23
W + WT) <
20

---

Under review as a conference paper at ICLR 2023
for all eigenvalues from L and (W + WT) Notice that 2,L is positive semi-definite and all its eigenvalues are within [0, 1/. Therefore, W guarantees the well-posedness of MIGNN as long as all eigenvalues satisfy
A 6 2 ( WT) S
When W (1 --)m CC F we have 2W W -1- - As CC positive semi-definite, all eigenvalues of
1 (W T) are no more than - (1-m).
E.4 PROOFS FOR SECTION
The following result about Kronecker product is adapted from 58 which we include here for com- pleteness.
Proof of Formula 9 used in Section f Since G' is symmetric, it admits eigen-decomposition GT = QαT AgT Qo- where QGT is orthogonal and hence satisfies QGT = Qa . As W is
diagonalizable, it admits eigen-decomposition W = AwwAQQw Then we can write
G @ W = =Qar AcT Q @ AwwAwQw Let n = (G) and = dim W), we have
[Qa @ QwllAGt & AwliQu @ Qw'
I In I, I [Qc I, L @ 1Qw @wll (Q QwliI, @ 1m @ Qw'l Therefore, for some matrix B E Rdxn
V (vve(())
1
I nd
G & W)
(vec(U))
1 1
1+0 a
1
Ind
a
G @ w)
(vec(U))
I+a
I+a
[QaT @ Qwl
Ind 1 Act & Aw
Q6 8 qn'l) (vec(U))
1+a
(Qat 8 @Qw]
Ind
a - @ Aw 1+a
[Q @ Qw'l (vec(U))
1+a a
Note that [nd AGT Aww] is diagonal matrix whose inverse is given by the diagonal matrix Diag(vec(H)) where the entires of H is given Hij 1/ (1 (Aw)a(AA )ss)
Here the notation Diag(u) denotes the diagonal matrix that has its diagonal for any vector From this we have,
1 (1QGT Qw] Diag(vec H)) [Q Qw') (vec(U))
((ve((U))
1+a
1
(QQT @ Qw Diag(vec (H)) UQ&T)
++
1
vec QQ @ Qw vec H (QWUQGT))
+
1
(Qw H
I+a
where ( denotes entry-wise multiplication.
For the reader's convenience, we present the following fact that implies D-1/2 + A +.... .... 1/2 D-122 has its eigenvalues within [-1,1] which is used in MIGNN with diffusion convolution
A D P (Equation 12}.
Proposition 5. Let S E RR be non-singular symmetric matrix and let D be the degree matrix defined as the diagonal matrix where D,s = Σn [Syl- Since is non-singular, D D-1/2 is well- defined. Then the normalization $ D-1/2 1/2 of S has its eigenvalues with (-1,1] 1].
21

---

Under review as a conference paper at ICLR 2023
Proof. Note that, the normalization S satisfies
D 1/2 1/2 D SD
1/2 = =,
that is, S is symmetric. To complete the proof, it then suffices to show that both 1 S and I I- S are positive semi-definite. Indeed, from the construction, both symmetric matrices are diagonal dominant, and their diagonal entries are positive, hence they are positive semi-definite by Gershgorin's Circle Theorem. Meanwhile, for any vector V E R", we have
S and + S
2 T (I Si = 0 (I + D DDDD//DDDD1-- 11/2
= D (DD 1/2.0
= D 1-1 T (D+
2 0
This shows that I + S is positive semi-definite. Similarly, one can derive that - S is positive semi-definite from D positive semi-definite.
F MIGNN VIA ANDERSON-ACCELERATTH OPERATOR SPLITTING SCHEMES
In this section, we present the pseudocodes of Anderson-aacclleaate MIGNN operator splitting schemes discussed Section 4
F.1 PSEUDOCODE FOR MIGNN WITH OPERATOR SPLITTING SCHEMES
FB Splitting. The detail of the FB splitting scheme iteration function Equation 6 of solving MIGNN is presented in Algorithm [
Algorithm fooward-MIGNN
Z 0: while em
( = (1 Z aW ZG + (g)(X)
(+) = prox}(22 (+)
|
Z Z
end while return
PR splitting, The details of the PR splitting scheme encoded in the iteration function Equation of solving MIGNN is presented in Algorithm 2
Algorithm 2 PR-forward-MIGNN
2, u ver(U)= = em: 1; V= aI+a G w))-
while en > edo
z (1/2) = -proxf (u)
u (1/2) = (1/2) -
(+) V(u (12 a vec(ggBX)))
u (+) = (+)+) u (12))
[ult)-ull
Cm = Nutt
2, u := end while
return prox.(u) (u)
F.2 MORE DETAILS ON BACKWARD PROPAGATION
In the backward propagation, the following result from (77] allows convert the computing of the inverse Jacobian term (I J & W))- to the (transpose of) matrix inverse term V (I GT WI which is already calculated in the forward pass.
22

---

Under review as a conference paper at ICLR 2023
Proposition 6 (Adapted from 177, Theorem 31). Let vec ) be the fixed point of the MIGNN
model 2 and J is the Jacobian a of the non-linearity at the GT @ Wvee + For any E IR" the solution u" of the equation
u = (I - J(G w))
is given by
u' = 2 (G W )
where u is a solution of the operator splitring problem 0 E (F + (t). with operators defined as
F(u)= = (I- W (ü), G(ü)= = Dü
(19)
where D is the diagonal matrix defined by J (I D) (where D if J 0).
Note that, since the non-linearity a is applied entry-wise, the Jacobian J is a diagonal matrix, and its diagonal entries consist of the vectorization of the Jacobian @r(W az z Therefore, the Jacobian J and hence D can be efficiently computed. We provide the pseudo-codes of FB and PR splitting schemes for the backward propagation described in the above proposition as Algorithm and Algorithm 4 respectively and their Anderson-acceeleaatee version can be found in Algorithm and Algorithm
FB backward propagation We now present the pseudo-code of FB splitting method (Algo- rithm 33) for the backward propagation with the procedure described in Proposition
Algorithm 3 FBB-bckward-MIGGNN
vec(U) (U) == 0:
er : 1:
Svec
while CTT > do
u = o vee W UG
+am,
if Da <,
0
if Da
en |
|x(+]]|
u == u
end while
Set U = (n)
return v veo (ec(WTul ue
Let u (k) be the intermediate variable, the procedure of applying FB splitting on monotone splitting problem 19 can be summarized as finding the fxee-point u of the following iteration function
u (++1) = BFB (u (k)) = (I + a ((1
+ +oW v).
(20)
PR backward propagation We now present the pseudo-code of PR splitting method (Algo- rithm 4 for the backward propagation with the procedure described in Proposition
Algorithm bR-backwardd-GGNN
= 0; u vec(U) = 0: em: 1;
imme(Z-)
V == (I+a(I - @ WW)
while er > ddo
u (1/2)
Bettn I+a(I+Das
if Da <
if Di
y (1/2) = (1/21
u (+) = V y (1/2)
y = 2u y (1/2)
| -yl
y. = y (+)
end while
if Da < 00
Compute u where un =
0
if Di := 00
Set U = vec 1 (u)
return v + vee W TUG)
23

---

Under review as a conference paper at ICLR 2023
Let (8) be the intermediate variable, the procedure of applying PR splitting Equation be summarized as first finding the fixed-point y of the following iteration function
199 can
y (8++) = B (g
2V (2(1 aD) aD) (k) + (v))
y( -2 -((++ - +a0) +y (k)
(21)
and then the final solution of the operator splitting problem is a 'y" +av).
F.3 ANDERSON ACCELERATION
We first introduce the general Anderson acceleration scheme. Let f IR R" be function the Lipschitz constant L(f) < 1. Therefore, the function f admits unique fixed point and can be obtained through Picard iteration. Let h(x)) f(x) x be the residual function. Let x(() be the initial guess, 3 (0, 1) be relaxation parameter, and m 1 be an integer parameter. Then the Anderson acceleration update x
x (k+1) = (1 B)
m m1
+ 3 Σ
7
(k)
(a(l-m)tt)
(22)
1=0
1=0
where the coefficients (A) = (70 (k) following:
..... mm
are determined by a least-square problem as the
min
Σ h (k--+))))))))) ))
s.t.
7 = 1.
(70... 7m)
1=0
Note that, when 1, the trivial weight (A) = (O....... 0, 1) recovers Picard iteration. Therefore, when the Picard iteration converges, the Anderson acceleration also converges and typically faster.
In Algorithm 5
we present the FB MIGNN forward propagation with Anderson acceleration on the FB iteration function FFB which is introduced in Section 4 and recalled here:
Z (k++)) = F Z (()
proxf (z(h) a (z WZ G (C-g((X)) (x())
Algorithm MIGNN--TB-Forwart FB MIGNN forward propagation
Input: initial point Z (0) = 0, FB damping parameter o, relaxation parameter 3, max storage size m 2 1.
Compute F(0) for k 1,... K do
F.B Z (0)).
H.
F(O)
Z (0)
Set m. min(m, Compute (+)) = r:P
PB
Z
H = (()) (k)
Update H = (H (k mk-
m.
H (k)
Determine (k) = 70
(4)
Tm
that solves
me
min Tmm....Tm
HHY| s.l. Σ = 1.
1=0
Set
mi
mg
z (k+1) =
Σ
%
(k)
p.
((k ((k- (1--)))) + (1 - Σ
%
(k)
(ik (+)
end for
return Z (1+1)
In Algorithm 6 we present the PR MIGNN forward propagation with Anderson acceleration on the
PR iteration function FPR which is introduced in Section 4 and recalled here:
(k+1) P.R PPP((u ()) 2V (u u (8) veo(ga(())) proxx (x) (b) (23)
24

---

Under review as a conference paper at ICLR 2023
Algorithm 6 MIGNN-PR-forward, PR MIGNN forward propagation
Input:
initial point ()) = vec : 0, PR damping parameter a, AA relaxation parameter
8. max storage size m 1.
Compute f(0) FFR (u (1))
(0) = f(() =
for 1,..., K do
Set my min (m, k) Compute (k) F (u
h (k) f(h) r
Update H (h' (k mm) ..... (k))
Determine (()) 70 ...., 7im
(k)
(8)
T
that solves
min Tmp
H s.t. Σ 76 1.
1=0
Set
mg
me
w (1+1)
B
Ys
(1)
FPR (2 ((k- ((--)))) + (1 (3)
Σ
(k)
(k-m))++)
0
1=0
end for
Set (k++)) = vec (z (x++1))
return prox3 prox3(U (U((-1))
The FB iteration function for the backpropagation BB is introduced in Appendix F.2 and recalled here:
(++1)
B. a ( (I + aD) ((1 - a) + +oW v).
(24)
We now present the Anderson-aaccelerated FB MIGNN backward propagation as Algorithm 7
Algorithm MIGGNN-FB-Backwar FB MIGNN backward propagation
Input: initial point (0) == vec(U) = 0, t = ōvec Z.j. ation parameter 3, max storage size 2 1.
PR damping parameter o, AA relax-
Compute (()) := B (u (0)), h (0) = f (0) for I..... K do
Set := min
Compute (M) BFB (R) (k) (k)
u
Update H (h (mm2))
h (8))
T
Determine (k) = 770 .eany 7mn
(k)
that solves
Hyll s.t. Σ 7/ 1.
min
(mp
1=0
Set
mi
u (++1) = 3
T BFB (2 (k-mm))+ + (1 - (3)
(k)
Σg
Y u (+==m)++))
(k)
==0
2=0
end for
Set (k+1) = vec vec return + vec(W] U[++1)
25

---

Under review as a conference paper at ICLR 2023
The PR iteration function for the backpropagation B is introduced in Appendix F.2 and recalled here: let y (k) be the intermediate variable,
y " ())) = (2(1+ aD) (t) av) - 2(I aD) (y av) yh) (25)
and then the final solution of the operator splitting problem is ū = (I aD)) (y" + tav). We now present the
PR MIGNN backward propagation Algorithm
Algorithm 8 MIGNN-PR- PR MIGNN backward propagation
Input: initial point y := 0, v == Ovect
88
PR damping parameter a, AA relaxation parameter
8, max storage size m 2 1.
Compute f(0)= : BP (y (0)), h0)) = f()) (y for .....: K do
Set := min (mi k)
Compute (k) B.P (y (k) h (k) f(k) y
Update H (h (z-
h (k)
(R) Determine (k) = ( 10 ...., Ymn
T
that solves
min Tmp
H Hyl s.t.
7/- 1.
1=0
Set
mi
mg
y (K+1)
B
(k) B ( (k-mm))+t)) + (1 - 3)
Σ
?
(k)
((k ((-mm)++))
0=0
==0
end for
y +att+ 1+a(1+Di))
Compute u (k+1) where u,
(k+1)
D;; < 00
=
0
if = 00
Set (k+1) == vec vec (m
return + vec(W W U(++1)
G EFFECTS OF THE ORDER OF NEUMANN SERIES EXPANSION
In this section, we perform ablation studies on the effects of the order of the Neumann series for approximating matrix (I a(I G & W))-1 in MIGNN-NKDP with fixed P P= 1. We study
the performance of MIGNN-NKDP for synthetic directed chain classification, benchmark graph node and graph classification.
G.1 DIRECTED CHAIN CLASSIFICATION
Examining the Neumann series expansion for the synthetic chain classification task demonstrates the trade-off between accuracy and time complexity. We train MIGNN-NKDI for three-class clas- sification, where the order K ranges from 1 to 5 in increments of 1. Fig. plots the resulting test accuracy, number of iterations, and time elapsed for each training epoch.
We make three observations as the order of the Neumann series increases. First the accuracy in- creases with respect to the order with diminishing returns. Second the number of iterations increases relative to the order up 3. Finally the time elapsed also increases with respect to the order up and 5 which are similar. These observations underscore the trade-off between accuracy and time complexity as the order increases.
26

---

Under review as a conference paper at ICLR 2023
Iterations
40.0
-0
Time (s)
0.3
37.5
20
35.0
0.1
100
200
400
660
800
1000
200
400
600
800
1000
200
4015
600
800
1000
Epoch
Epoch
Epoch
Figure 9: Comparison of Neumann expansion for accuracy, number of iterations, and elapsed time using three-class chain classifications with chain length 140.
G.2 NODE CLASSIFICATION
The graph node classification tasks also highlight the trade-off between accuracy and time complex- ity. We train MIGNN-NKDI using 10-fold cross validation on Cora, Citeseer and Pubmed. We consider K in the range from 5, incrementing along with their standard deviations are reported in Table
The mean test accuracy and time elapsed
Deteets
Cora Acarary)
Cor (1 Fimeo
Choeer Actumayi
Chorerr chorer)llree)
I'lmall (Aceimes)
Pnbrred (imel
U1-154
1--11
13377
1.11
MIGNS NID1
88.1
0.016
MIGNS- 9301
85.8
5.55
0.511
0-022
73.88
5.2
0-232
0.022
83.3
0:78
0.211
0.017
MIGNN-N501
67.0
0.895
0.064
74.88 2.3
0.282
0.015
83.8 0.66
0.272 0.024
Table 3: Graph node classification mean accuracy (%) +: standard deviation for 10-fold vros-validatiin.
For node classification
a very clear trend across all datasets. Both the accuracy and time clapsed increase with the order of the Neumann expansion. However, the accuracy scales with diminishing returns; notice N4 and N5 have the same accuracy for both Citeseer and Pubmed.
G.3 GRAPH CLASSIFICATION
In this subsection, apply MIGNN-NKDI to classify the MUTAG dataset, where K ranges from to incrementing by 1. Fig. 10 plots the test accuracy, the number iterations, and the time elapsed for training one fold of the 10-fold cross validation.
90
30
Time Ellapsed (s)
H503
Itorations
Accuracy
30
70
20
Epoch
Epoch
Epoch
Figure 10: Comparison of Neumann expansion for accuracy, number of iterations and elapsed time using the first fold the MUTAG graph data set.
Unlike the directed chains and node classification tasks, the graph classification does not show sig. nificant improvements from higher order Neumann expansion on this fold. However, from Table 2 observe that over 10-fold cross validation diffusion improves the results. Although the accuracy and iteration count remain similar among all orders, the time elapsed still scales with the order.
H EFFECTS OF THE ORDER OF GRAPH DIFFUSION CONVOLUTION
In this section, we use MIGNN-NKDP with fixed K 1 and varying order of graph diffusion P study the effects of the order of graph diffusion convolution. We report the performance of MIGNN benchmarking on synthetic directed chain classification, benchmark graph node and graph classification tasks.
H.1 DIRECTED CHAIN CLASSIFICATION
The three-clase chain classification task benefits tremendously for high orders of diffusion. We train MIGNN-NIDP on chain lengths of 140, where P ranges from to 5 incrementing by Fig. plots the test accuracy, number of iterations, and time elapsed for each training epoch.
27

---

Under review as a conference paper at ICLR 2023
70
Iteraaions
0.14
Accuracy
00.12
Time
15
20.10
40
30
10
0.08
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Epoch
Epoch
Epoch
Figure 11:
Comparison of graph diffusion convolution for accuracy, numher iterations and clapsed time
using three-class chains of length 140.
For diffusion convolution we make two observations. First, the accuracy scales with the order of diffusion with remarkable gap between D3 and D4, Second, the iteration count and time elapsed remain relatively constant among all orders, with DI standing out as the least among all others.
Our theory informs of the following: 1) Accuracy scaling occurs when the introduced P-hop edges contain relevant information for the task 2) Time elapsed scales relative to the number of edges in the higher order graph diffusion matrix. Our observations support our theory and strongly suggest using diffusion as an inexpensive improvement to simple learning tasks.
H.2 NODE CLASSIFICATION
this subsection, we study the effects of the order diffusion convolution on the node classification tasks outlined in the citation datasets (Cora, Citeseer, Pubmed), We consider MIGNN-NIDE with reports the test accuracy and the time elapsed for each epoch for different MIGNN models. We observe that diffusion does provide any benefit for graph node classification.
P ranging from to with an increment of 1. Table Treporrr
(Accuracy)
(Accuracy)
(Timc)
(Accuracy)
NIGNN--NIDI
3 1.81
0.036
09.3
0.141
0.017
80.9
3.87
0.151
0.015
MIGNN-NIDE
86.5
1,30
0.367
0.032
69.3
6.5
0.146
0.021
77.8
4.82
0.410
0.040
MIGNN-NIDS
83.7 + 1.33
0.766
+ 0.021
68.0
+ 7.0
0.164 = 0.057
83.3 - 0.76
6.25 1.14
Table 4: Graph node classification mean accuracy (%) + standard deviation for 10-fold cross- validaaaii.
H.3 GRAPH CLASSIFICATION
We further apply MIGNN-NIDP to classify the MUTAG dataset, where P ranges from I to incrementing by 1. Fig. 12 plots the test accuracy, number of iterations, and time elapsed for each epoch for different MIGNN models.
-0
RDDS
10D4
HIDS
Accuracy
Itorations
Time Ellapsed
85
0.6
60
20
40
Epoch
Epoch
Epoch
Figure 12:
Comparison of diffusion convolution for accuracy, number of iterations and clapsed time using the first fold the MUTAG graph data set.
We observe that higher order diffusion convolution has little impact the time complexity when each connected subgraph is small relative to the underlying graph.
28

---

Under review as a conference paper at ICLR 2023
I MORE DISCUSSION ON WHEN IGNNS BECOME EXPRESSIVE FOR LEARNING LRD
In this section, we further confirm the interconnection between the accuracy of IGNN for classifying directed chains and the eigenvalues W.. The accuracy and number of iterations of IGNN and the dynamics of the two leading eigenvalues are plotted in Figs. 13 and 14 respectively, for the binary and three-class cases. These results confirm the phenomena we have discussed Sec.
1.0
100
MMVVV
Accuracy(%)
80
Az(lW)
Iterations
200
0.6
60
100
Chain Length
100
100
200
300
400
500
0
100
200
300
400
500
100
200
300
400
500
Epoch
Epoch
Epoch
(a)
(b)
(C)
100
300
Accuracy(%%
Iterations
Ax(W
200
100
Chain Length
-
0.4
100
200
300
400
500
100
200
300
400
500
100
705
300
400
500
Epoch
Epoch
Epoch
(d)
12
Length
Accuracy
As(/W))
# Iterations
10
+ 200
Aa((W
0.2
50
48
100
200
100
400
500
100
200
300
400
500
100
200
300
400
500
Epoch
Epoch
Epoch
(g)
(h)
(i)
10
Length
Accuracy
52
0.3
(s(WW Ao((W
Iterations
250
50
0.1
100
200
360
400
500
100
200
300
400
500
100
200
300
400
500
Epoch
Epoch
Epoch
(j)
(k)
(l)
0.4
12.5
Length
Accuracy(%)
50
Au(IW
A,(WW
* Iterations
10.0
+ 300
0.2
7.5
88
0.1
5.0
0.0
2.5
100
200
300
400
500
100
200
300
400
500
200
300
400
500
Epoch (m)
Epoch (n)
Epach (0)
Figure 13: first column, training, test, and validation accuracies of IGNN are depicted for several varying chain lengths. In the second column, the corresponding top two eigenvalues are plotted. The third column depicts the number of Picard iterations for each chain length. When IGNN becomes accurate for chain
classification, the corresponding A ([W)) becomes close I and requires substantially more iterations for the Picard iteration to converge.
29

---

Under review as a conference paper at ICLR 2023
0.5
42
Chain Length
Accuracy(50
# Iterations
40
0.3
Ag([W]
10
38
0.2
36
0.1
200
300
500
Epoch
Epoch
Epoch
(a)
(b)
(c)
100
300
Accuracy
0.8
# Iterations
80
Ax((W)
200
60
0.6
100
Chain Length
Train - Validate - Test
0.4
120
400
500
100
200
300
400
500
700
300
400
500
Epoch
Epoch
Epoch
(d)
(e)
(f)
0.6
40
12.5
Length
Accuracy(%8
36
0.4
(1()
Ax([W
lterations
140
36
7.5
0.2
5.0
34
100
200
300
400
500
100
200
300
400
500
2.50
Epoch
Epoch
Epoch (i)
(g)
(h)
40
0.4
Length
Accuracy
Ad((W)
Itteraionns
10
160
28
0.3
36
0.2
34
0.1
100
200
360
400
500
100
200
300
400
500
100
200
300
400
500
Epoch
Epoch
Epoch (1)
(j)
(k)
0.5
10
Chain Length
Accuracy(9)
38
0.4
Iterations
+ 180
0.3
0.2
32
100
200
300
400
500
200
300
400
500
100
200
300
400
500
Epoch (m)
Epoch (n)
Epoch (0)
Figure 14: The first column shows the training, test, and validation accuracies of IGNN for several chain lengths of three classes. column, we plot the number of Picard iterations for each chain length. As the maximum eigenvalue of the sys- tem approaches 1, IGNN becomes more accurate for chain classification at the cost of a significantly increased number of training iterations.
second column, we plot the corresponding top two eigenvalues. In the third
J DETAILS ABOUT DATASETS
Synthetic chains dataset. To evaluate the LRD learning ability of models, construct synthetic chains dataset as in Gu et al. [39]. Both binary classification and multiclass classification are con- sidered. Let be the number of classes, that there are types of chains. The label information is only encoded as a one-hot vector in the first c-dimensions of the node feature of the starting nodes of each chain. With classes, ne chains for each class, and nodes in each chain, the chain dataset has x x nodes in total.
Bioinformatics datasets. MUTAG is a dataset of 188 mutagenic aromatic and heteroaromatic nitro compounds. PTC is dataset of 344 chemical compounds that report carcinogenicity for male
30

---

Under review as a conference paper at ICLR 2023
and female rats. COX2 is dataset of 467 cyclooxygenase-2 (COX-2) inhibitors. PROTEINS is dataset of 1113 secondary structure elements (SSEs). NCII is a public dataset from the National Cancer Institute (NCI) and is subset of balanced datasets of chemical compounds screened for the ability to suppress or inhibit the growth of a panel of human tumor cell lines.
Amazon product co-purchasing network. This dataset contains 334863 nodes (representing goods), 925872 edges, and 58 label types. An edge is formed between two nodes if the represented goods have been purchased together [52].
Pore networks. The pore network is simulated dataset that models fluid flow in porous media. Each porous network is randomly generated inside a cubic domain of width 0.1m by Delaunay Voronoi tessellation. The prediction of equilibrium pressure in a pore network under physical diffusion is introduced as GNN task in [63]. The GNN model prediction accuracy is compared with the ground truth obtained through solving the diffusion equation directly, see 63 Appendix CJ for more details.
Citation dataset. Cora and Citeseer are large citation datasets that deseribe the presence of specific words in publications. Pubmed is large citation dataset that contains information about papers classified for studying one of the three diabetes. The following table adapted from describes the statistics of the three datasets.
Dataset
Type
Classes
Features
Nodes
Edges
Label rate
Avg. SP
Cora
Citation
2879
2810
7981
0.047
5.27
Citeseer
Citation
6
3703
2110
3668
0.036
9.31
Pubmed
Citation
500
19717
44324
0.003
6.34
Table 5: Dataset statisties. The shortest path length is denoted by Avg. SP.
K DETAILS ABOUT HYPERPARAMETERS
The default parameter settings for MIGNN are the following. For the fixed-point schemes 0.9, 3 = 0.9, the default maximum iteration is 300, the tolerance is le-6, and convergence is mea- sured in the -nom of the difference between two consecutive fixed point iterations. The learnable parameter is initialized to y = 1.0.
Synthetic chains dataset. For both binary and three-class classification we use the parameters outlined by IGNN 77I. In both classification tasks we make the same modifications. We set the clipping and dropout to
Citation dataset. In the citation datasets follow the training procedure used by GIND 22 For the Cora dataset we set the weighted-decay to le-4 and the fixed-point tolerance to 1e-3. For all three models we set the fixed-point
0.5, and the number of hidden layers to 64.
Bioinformatics datasets. In the bioinformatics datasets we follow the training procedure used by IGNN 77 On the C12 dataset for MIGNN-Mon, we use a 0.5. On all other datasets we extend the number of training epochs to 500.
Amazon product co-purchasing network. In the Amazon product co-purchasing dataset we fol- low the training procedure used by IGNN 71
Pore networks. In the physical diffusion pore networks we follow the training procedure used by CGS [63] and the default parameter values for MIGNN.
31


Reviews:
[{"review_id": "rZm8_-QzK2I", "reviewer": "Reviewer_h7JA", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes to improve the stability, accuracy, and efficiency of implicit graph neural networks (IGNN) by new parameterizations (i.e., the Cayley transform-based orthogonal parameterization and monotone parameterization) and advanced solvers (i.e., operator splitting, Anderson acceleration). Theoretical justification for the well-posedness for the proposed MIGNN is presented, and empirical experiments demonstrate the performance of MIGNN.", "summary_of_the_review": "The paper introduces interesting observations and ideas to improve IGNN. However, the contribution and significance are unclear, and the evaluation is not convincing enough to justify the effectiveness of the proposed algorithm.\n\n## After rebuttal\nThe revision significantly improves the paper. I am willing to increase my score.", "strength_and_weaknesses": "# Strength\n\n1. The paper makes a nice observation that the accuracy of IGNN seems to correlate with the eigenvalue of $W$ in the implicit model. This is an interesting observation and motivates the work in the paper.\n\n2. The paper generalizes the well-posedness condition for IGNN and proposes more flexible constructions based on monotone operator theory. New parameterizations are proposed to define a new model MIGNN.\n\n3. The paper comprehensively discusses the technical details of related works, which is helpful to understand the idea in the paper. \n\n4. Extensive experiments are presented with discussions.\n\n# Weakness\n\n1. The contribution and significance of the proposed ideas are unclear. For instance, many techniques in MIGNN have been extensively studied in the literature and even in the study of GNNs, such as Anderson acceleration, Neumann series approximation, high-order graph diffusion, eigendecomposition, etc. Moreover, monotone operator theory and operator splitting are not new (even in the context of deep learning). Therefore, the proposed idea seems an incremental combination of many techniques but the main contribution is unclear.\n\n2. Among all the techniques in the proposed MIGNN, the flexible parameterization in Section 3 seems to be new and interesting. However, the paper fails to justify the effectiveness of this key innovation due to the lack of a comprehensive ablation study. The improvements over IGNN may be due to other existing techniques that can be trivially applied in IGNN as well such as Anderson acceleration, Neumann series approximation, high-order graph diffusion, eigendecomposition, etc. In fact, Figure 3 partially confirms this since the high-order diffusion in IGNN-D5 significantly improves the performance of IGNN. \n\n3. The comparison with IGNN is unfair. In fact, each forward iteration of MIGNN has more computation and graph aggregations than each forward iteration of IGNN due to the Neumann series approximation and diffusion convolutions. For instance, it is mentioned in the paper that \"Each node can access information from its K-hop neighbors using the K-th order Neumann series approximated PR iteration\". Therefore, the concept of iteration needs to be clearly and fairly defined when comparing their stability, accuracy, and efficiency. \n\n4. There is a lack of computation complexity analysis. In fact, the computation cost of MIGNN is pretty high, and it is unclear how it outperforms IGNN or EIGNN theoretically and empirically.\n\n5. The motivation for monotone parameterization in Section 3.2 is unclear. Why do you define G as L/2? What are the intuition and advantages？\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "g9y9dooYm12", "reviewer": "Reviewer_Ti7k", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This work proposes to use monotone and orthogonal operator to build implicit graph neural network. The adopted approach is well build upon monotone operator theory. The new operators allows more stable (via orthogonal operator) or expressive (via monotone operator) implicit GNN models with provable convergence. The discussion is super extensive. The experiments on node-level and graph-level tasks show the superiority of the method.  ", "summary_of_the_review": "The paper is strong in its technical contributions. The empirical contributions are about the acceptance bar. More extensive studies are suggested. The exposition of main text can be improved by showing the results more directly instead of detouring the discussion. ", "strength_and_weaknesses": "Strengths\n1. This work is strong in its math-grounded modeling technique. In particular, it connects GNN modeling with monotone operator theory, which is exciting. Although such connection has been leveraged to study deep equilibrium model [1], I think it is still novel to apply such connection to study GNN models. \n\n2. The derivation and the discussion are super extensive, covering different types of operators to match different properties, different operator splitting iterative methods, acceleration ideas, back-propagation ideas.   \n\n3. The experiments make sense and show the benefits of the model.\n\nWeaknesses\n1. The main concern of this paper is about its exposition. As it tries to discuss too many things, many details in the main text are missed and I have to check the appendix even if I am fairly familiar with monotone operator theory. Some discussions may be not necessary, kind of detour, and may distract the readers. \n\nFor example, in Sec. 2, the authors advertise a lot about the new well-posed condition, where they claim that Eq. (2) is well posed even if $W$ has eigenvalue less than -1. However, this result is not very practically useful, because in sec. 3, their finally adopted $W$ still have to have eigenvalues with absolute values less than 1 to guarantee convergence. Moreover, in Remark 2, the way to make W symmetric will naturally remove the asymmetric part of W in the monotone operator ($F- F^T$). \n\n2. Some ablation studies are missing. For example, the authors also adopt diffusion convolution in Eq.(12), where the diffusion convolution automatically has the benefit in capturing long range dependence. So, it is unclear that in the node classification tasks, whether the benefit comes from diffusion convolution or implicit operators. Here,  I also do not know what \"N3D5\", \"N5D5\" used in the model name in the experiment means, which are not defined. I assume they are related to the diffusion convolution. Please clarify. Also, provide the study of using different convolutions to show the benefits indeed come from implicit modeling. \n\n3. Node-level experiments are too few. Actually, Long range dependence is not that crucial for graph-level tasks, because one can adopt transformer models and run on dense graphs by adding jump-hop edges. This is also why the proposed models are hard to outperform many baselines in graph-level tasks (the used baselines and datasets in this work are actually not SOTA and large). LRD is more important for node-level tasks. However, this work only provides one dataset for such evaluation. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "FooUNaGqqO", "reviewer": "Reviewer_FLVV", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper presents a method to stabilize the training of implicit graph neural networks (gnns) and in less required iterations. \nThe proposition is based on learning a monotone operator and employing splitting techniques for a better solution of the iterative problem.\n\nSeveral experiments are conducted and show the improvement of the proposed method compared to previous implicit gnns.", "summary_of_the_review": "The paper suggest an interesting method to improve implicit graph neural networks and the results compared to such methods seem promising. However there is a lack of discussion with other relevant works, and the experimental part is limited and hard to compare with other methods.", "strength_and_weaknesses": "Pros:\n1. The authors provide a sufficient background of implicit gnns\n2. The paper is well motivated and mostly easy to follow\n3. There are several interesting experiments \n\nCons:\n1. The paper lacks discussion of relevant prior work. For example, the use of Cayley parameterization was previously proposed in \"CayleyNets: Graph Convolutional Neural Networks with Complex Rational Spectral Filters\". In the sense of learning long-range dependencies several works proposed various approaches, for instance learning multiple hop filters \"Path Integral Based Convolution and Pooling for Graph Neural Networks\" or dedicated oscillatory layers in \"PDE-GCN: Novel Architectures for Graph Neural Networks Motivated by Partial Differential Equations\" and in \"Graph-Coupled Oscillator Networks\".\n\n2.While the experiments do show the significance of the proposed method compared to existing implicit gnns, it is hard to evaluate the performance of the proposed method compared to recent work in the field of gnns due to two reasons.\nA. The authors do not compare where possible to recent methods. For example table 1 lacks many works that perform significantly better. See for example \"Improving graph neural network expressivity via subgraph isomorphism counting.\", \"Weisfeiler and lehman go topological: Message passing simplicial networks.\" , \"Weisfeiler and lehman go cellular: Cw networks.\" , \nB. Besides this experiment, all other datasets are not directly comparable with other works. This is not bad, but to really assess the performance of this method I think that it needs to include experiments that are more directly comparable like node classification on Cora,Citeseer,Pubmed,OGBN-ARXIV and others. \n\n3.Regarding equation 1. This formulation of adding the initial features looks similar to the operation in \"Simple and Deep Convolutional Graph Neural Networks\" which is not discussed in this paper. Can the authors explain?\n\n4. Regarding the claim at the bottom of the first page. I am not convinced that it is true that IGNN does not suffer from the very same problem. What is the difference here? Why should IGNNs that smooth the node features not oversmooth if the propagation is the same?\n\n5.I understand that the proposed method requires less iterations, but I do not know how much time each iteration of the proposed method takes. I believe that if the authors add run times their results will be more convincing.\n\n6.Proposition 2 is not clear - where does the final eqation W = ... comes from?\n\n7.One of the hyperparameters of this method is the order of diffuion. However in different experiments the authors use different orders. It is hard to deduce what is the actual influence of the order this way. I think that the paper can benefit from an ablation study of the various hyperparameters.\n\n8.What are the chosen hyperparameters in your experiments? How did you choose them?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "6nbQDuYfhF", "reviewer": "Reviewer_Y8Ag", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposed a variant of Implicit Graph Neural Network (IGNN) based on the monotone operator theory. First, this paper formalized the convergence of IGNN as a monotone inclusion problem and derived the sufficient condition for well-posedness. Then, two types of parameterization for the weight matrix were proposed based on this condition: the orthogonal parametrization and the monotone parametrization. In addition, to speed up the convergence of IGNN, this paper proposed two acceleration methods based on the forward-backward splitting scheme and the Pearceman-Rachford splitting scheme, respectively. Finally, this paper applied the proposed model to the chain classification task on synthesis data and to the graph and node prediction tasks on real data to assess the prediction accuracy and computational efficiency.", "summary_of_the_review": "Theoretical part is clearly written and sound. Also, it is accessible to those who are not familar with the monotone operator theory.\nIf I understand correctly, most improvements are the application of existing techniques for computing fixed points by iterative algorithms, except for the monotone parameterization. Also, I have questions about the design of experiment settings and the interpretation of the experiment results.", "strength_and_weaknesses": "【Strengths】\n- Provides basics of monotone operator theory, which enables the paper to be more accessible to those unfamiliar with monotone operator theory.\n\n\n【Weaknesses】\n- Novelty of the proposed method is limited because the methodology primarily relied on the combination of existing methods.\n- Interpretation of $\\lambda_1(W)$, which may be related to IGNN convergence and prediction accuracy, is somewhat questionable.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject