CONVERGENT GRAPH SOLVERS
Year: 2022

Abstract:
We propose the convergent graph solver (CGS)1, a deep learning method that learns iterative mappings to predict the properties of a graph system at its stationary state (fixed point) with guaranteed convergence. The forward propagation of CGS proceeds in three steps: (1) constructing the input-dependent linear contracting iterative maps, (2) computing the fixed points of the iterative maps, and (3) decoding the fixed points to estimate the properties. The contractivity of the constructed linear maps guarantees the existence and uniqueness of the fixed points following the Banach fixed point theorem. To train CGS efficiently, we also derive a tractable analytical expression for its gradient by leveraging the implicit function theorem. We evaluate the performance of CGS by applying it to various network-analytic and graph benchmark problems. The results indicate that CGS has competitive capabilities for predicting the stationary properties of graph systems, irrespective of whether the target systems are linear or non-linear. CGS also shows high performance for graph classification problems where the existence or the meaning of a fixed point is hard to be clearly defined, which highlights the potential of CGS as a general graph neural network architecture.

heading: 1 INTRODUCTION
text: Our world is replete with networked systems, where their overall properties emerge from complex interactions among the system entities. Such networked systems attain their unique properties from their stationary states; hence, finding these stationary properties is a common goal for many problems that arise in the science and engineering field. Examples of such problems include the minimization of the molecule’s potential energy that finds the stationary positions of the atoms to compute the potential energy of the molecule (Moloi & Ali, 2005), the PageRank algorithm that finds the stationary importance of online web pages to compute recommendation scores (Brin & Page, 1998), and the network analysis of fluid flow in porous media that finds the stationary pressures in pore networks to compute the macroscopic properties of the media (Gostick et al., 2016). In these network-analytic problems, the network is often represented as a graph, and the stationary states are often computed through a problem-specific iterative method derived analytically from the prior knowledge of the target system. By applying the iterative map repeatedly, these methods compute the stationary states (i.e., fixed points) and the associated properties of the target system. However, analytically deriving such iterative maps for highly complex systems typically requires tremendous efforts and time.
Instead of deriving/designing the problem-specific iterative methods, researchers employ deep learning approaches to learn the iterative methods (Dai et al., 2018; Gu et al., 2020; Hsieh et al., 2019; Huang et al., 2020). These approaches learn an iterative map directly without using domain-specific knowledge, but using only the input and output data. By applying the learned iterative map, these approaches (approximately) compute the fixed points and predict the properties of a target system. Graph Neural Networks (GNNs) have been widely used to construct iterative maps (Dai et al., 2018; Gu et al., 2020; Scarselli et al., 2008; Alet et al., 2019). However, when applying this approach, the existence of the fixed points is seldom guaranteed. As a result, the number of iterative steps is often required to be specified, as a hyperparameter, to ensure the termination of the iterations. This may cause premature termination or inefficient backward propagation if the number of iterations is inadequately small or too large.
1The code is available at https://github.com/Junyoungpark/CGS.
In this study, we propose a convergent graph solver (CGS), a deep learning method that can predict the solution of a target graph analytical problem using only the input and output data, and without requiring the prior knowledge of existing solvers or intermediate solutions. The forward propagation of CGS is designed to proceed in the following three steps:
• Constructing the input-dependent linear-contracting iterative maps. CGS uses the input graph, which dictates the specification of the target network-analytic problem, to construct a set of linear contracting maps. This procedure formulates/set up the internal problem to be solved by considering the problem conditions and contexts (i.e., boundary conditions or initial conditions in PDE domains – the physical network problems). Furthermore, the input-dependent linear map can produce any size of transition map flexibly depending on the input size graph; thus helping the trained model to generalize over unseen problems with different sizes (size transferability).
• Computing the fixed points via iterative methods. CGS constructs a set of linear contracting maps, each of which is guaranteed to have a unique fixed point that embeds the important features for conducting various end tasks. Thus, CGS computes the unique solutions of the constructed linear maps via iterative methods (or direct inversion) with convergence guarantee.
• Decoding the fixed points to estimate the properties. By using a separate decoder architecture, we compute the fixed points in the latent space while expecting them to be an effective representation that can improve the predictive performance of the model. This enables CGS to be used not only for finding the real fixed points (or its transportation) if they exist, but also for computing the "virtual fixed point" as a representation learning method in general prediction tasks.
The parameters of CGS are optimized with the gradient computed based on the implicit function theorem, which requires O(1) memory usage when computing the gradient along with the iterative steps. CGS is different from the other studies that solve the constrained optimization (with convergence guarantee) (Gu et al., 2020; Scarselli et al., 2008; Tiezzi et al., 2020) in that it does not impose any restriction when optimizing the parameters. Instead, CGS is inherently structured to have the unique fixed points owing to the uses of the linear map. Note that the structural restriction is only in the form of an iterative map; we can flexibly generate the coefficients of the linear map using any network (i.e., GNN) and utilize multiple linear maps to boost the representability.
We evaluate the performance of CGS using two types of paradigmatic network-analytic problems: physical diffusion in networks and the Markov decision process, where the true labels indeed exist and can be computed analytically using linear and non-linear iterative methods, respectively. We also employ CGS to solve various graph classification benchmark tasks to show that CGS can serve as a general (implicit) layer like other GNN networks when conducting typical graph classification tasks. In these experiments, we seek to compute a virtual fixed point that can serve as the best hidden representation of the input when predicting the output. The results show that CGS can be (1) an effective solver for graph network problems or (2) an effective general computational layer for processing graph-structured data.

heading: 2 RELATED WORK
text: Convergent neural models. Previous studies that have attempted to achieve the convergence property of neural network embedding (e.g., hidden vectors of MLP and hidden states of recurrent neural networks) can be grouped into two categories: soft and hard approaches. Soft approaches typically attain the desired convergence properties by augmenting the loss functions (Erichson et al., 2019; Miyato et al., 2018). Although the soft approaches are network-architecture agnostic, these methods cannot guarantee the convergence of the learned mappings. On the other hand, hard approaches seek to guarantee the convergence of the iterative maps by restricting their parameters in certain ranges (Gu et al., 2020; Tiezzi et al., 2020; Miller & Hardt, 2018; Kolter & Manek, 2019). This is achieved by projecting the parameters of the models into stable regions. However, such projection may lead to non-optimal performances since it is performed after the gradient update, i.e., the projection is disentangled from training objective.
Implicit deep models. The forward propagation of CGS, which solves fixed point iterations, is closely related to deep implicit models. Instead of defining the computational procedures (e.g., the depth of layer in neural network) explicitly, these models use implicit layers, which accept inputdependent equations and compute the solutions of the input equations, for the forward propagation.
For instance, neural ordinary differential equations (NODE) (Chen et al., 2018; Massaroli et al., 2020) solve ODE until the solver tolerance is satisfied or the integration domain is covered, the optimization layers (Gould et al., 2016; Amos & Kolter, 2017) solve the optimization problem until the duality gap converges, and the fixed point models (Bai et al., 2019; Winston & Kolter, 2020) solve the network-generated fixed point iterations until some numerical solver satisfy the convergence condition. Implicit models use these (intermediate) solutions to conduct various end tasks (e.g., regressions, or classifications). In this way, implicit models can impose desired behavioral characteristics into layers as inductive bias and hence, often show superior parameter/memory efficiency and predictability.
Fixed points of graph convolution. Methods that find the fixed points of graph convolutions have been suggested in various contexts of graph-related tasks. Several works have utilized GNN along with RNN-like connections to (approximately) find the fixed points of graph convolutions (Liao et al., 2018; Dai et al., 2018; Li et al., 2015; Scarselli et al., 2008). Some have suggested to constrain the parameter space of GNN so that the trained GNN becomes a non-expansive map, thus producing the fixed point (Gu et al., 2020; Tiezzi et al., 2020). Others have proposed to apply an additional GNN layer on the embedded graph and penalize the difference between the output of the additional GNN layer and the embedded graph to guide the GNN to find the fixed points (Scarselli et al., 2008; Yang et al., 2021). It has been shown that regularizing GNN to find its fixed points improves the predictive performance (Tiezzi et al., 2020; Yang et al., 2021).
Comparison between CGS and existing approaches. Combining the ideas of (1) computing/using the fixed points of the graph convolution (as a representation learning) and (2) utilizing the implicit differentiation of the model (as a training method) has been proposed by numerous studies (Scarselli et al., 2008; Liao et al., 2018; Johnson et al., 2020; Dai et al., 2018; Bai et al., 2019; Gallicchio & Micheli, 2020; Bai et al., 2019; Tiezzi et al., 2020; Gu et al., 2020). Majority of those studies assume that the convolution operators (iterative map) induce a convergent sequence of the (hidden) representation; however, this assumption typically does not hold unless certain conditions hold for the convolution operator (iterative map). When neither the convergence nor uniqueness holds, it can possibly introduce biases in the gradient computed by implicit differentiation (Liao et al., 2018; Blondel et al., 2021). Hence, to impose convergence, some methods restrict the learned convolutions to be strictly contractive (i.e., the hidden solutions are convergent) by projecting the learned parameters into a certain region and solving the constraint training problems, respectively (Gu et al., 2020; Tiezzi et al., 2020). Unlike these methods restricting the parameters of the graph convolutions directly, CGS guarantees the convergence and the uniqueness of the fixed point by imposing the structural inductive bias on the iterative map, i.e., using the contractive linear map that has a unique fixed point, thus alleviating the need to solve the constrained parameter optimization (Tiezzi et al., 2020).

heading: 3 PROBLEM DESCRIPTION
text: The objective of many network-analytic problems can be described as:
Find a solution vector Y ∗ from a graph G that represents the target network system.
In this section, we briefly explain a general iterative scheme to compute Y ∗ from G. The problem specification G = (V,E) is a directed graph that is composed of a set of nodes V and a set of edges E. We define the ith node as vi and the edge from vi to vj as eij . The general scheme of the iterative methods is given as follows:
H [0] = f(G), (1) H [n] = T (H [n−1];G), n = 1, 2, ... (2) Y [n] = g(H [n];G), n = 0, 1, ... (3)
where f(·) is the problem-specific initialization scheme that transforms G into the initial hidden embedding H [0], T is the problem-specific iterative map that updates the hidden embedding H [n] from the previous embedding H [n−1], and g(·) is the problem-specific decoding function that predicts the intermediate solution Y [n].
T is designed such that the fixed point iteration (equation 2) converges to the unique fixed point H∗:
lim n→∞
H [n] = H∗ s.t. H∗ = T (H∗,G) (4)
The solution Y ∗ is then obtained by decoding H∗, i.e., g(H∗;G) ≜ Y ∗. In many real-world network-analytic problems, we can obtain G and its corresponding solution Y ∗, but not T , H [n] and Y [n]. Therefore, we aim to learn a mapping T from G to Y ∗ without using T , H [n] and Y [n].

heading: 3.1 EXAMPLE: GRAPH VALUE ITERATION
text: Let us consider a finite Markov decision process (MDP), whose goal is to find the state values through the iterative applications of the Bellman optimality backup operator (Bellman, 1954). We assume that the state transition is deterministic.
We define G = (V,E), where V and E are the set of states and transitions of MDP respectively. vi corresponds to the ith state of the MDP, and eij corresponds to the state transition from state i to j. eij exists only if the state transition from state i to j is allowed in the MDP. The features of eij are the corresponding state transition rewards. The objective is to find the state values V ∗. In this setting, the Bellman optimal backup operator T is defined as follows:
V [n] i = T (V [n−1];G) ≜ max j∈N (i) (rij + αV [n−1] j ) (5)
where V [n]i is the state value of vi estimated at n-th iteration, N (i) is the set of states that can be reached from the ith state via one step transition, rij is the immediate reward related to eij , and α is the discount rate of the MDP.
In this graph value iteration (GVI) problem, the initial values V [0] are set as zeros (i.e., f(·) = 0), then T is applied until V [n] converges, and V [n] is decoded with the identity mapping (i.e., g(·) is identity and H [n] ≜ V [n]). We will show how CGS constructs the transition map T from the input graph and predicts the converged state values V ∗ without using equation 5 in the following sections.

heading: 4 CONVERGENT GRAPH SOLVERS
text: We propose CGS that predicts the solution Y ∗ from the given input graph G in three steps: (1) constructing linear iterative maps Tθ from G, (2) computing the unique fixed points H∗ of Tθ via iterative methods, and (3) decoding H∗ to produce Y ∗, as shown in Figure 1.

heading: 4.1 CONSTRUCTING LINEAR ITERATIVE MAPS
text: CGS first constructs an input-dependent contracting linear map Tθ(· ;G) such that the repeated application of Tθ(· ;G) always produces the unique fixed point H∗ (i.e. limn→∞ H [n] = H∗) that embeds the essential characteristics of the input graph G for conducting end tasks. In other words, CGS learns to construct iterative maps that are tailored to each input graph G, from which the unique fixed points of the target system are guaranteed to be computed and used for conducting end tasks.
To impose the contraction property on the linear map, CGS utilizes the following iterative map:
Tθ(H [n];G) ≜ γAθ(G)H [n] +Bθ(G) (6)
where γ is the contraction factor, Aθ(G) ∈ Rp×p is the input-dependent transition parameter, Bθ(G) ∈ Rp is the input-dependent bias parameter, and p is the number of nodes in graph. To construct an input-dependent iterative map that preserves the structural constraints, which is required to guarantee the existence and uniqueness of a fixed point, we employ GNN-based parametergenerating network fθ(·). The parameter generation procedure for Tθ starts by encoding G = (V,E): V′,E′ = fθ(V,E) (7) where V′ and E′ are the set of updated node embeddings v′i ∈ Rq and edge embeddings e′ij ∈ R respectively. CGS then constructs Aθ(G) by computing the (i, j)th element of Aθ(G) as follows:
[Aθ(G)]i,j =
{ σ(e′ij)
d(i) if eij exists, 0 otherwise.
(8)
where σ(x) is a differentiable bounded function that projects x into the range [0, 1] (e.g., Sigmoid function), and d(i) is the outward degree of vi (i.e.„ the number of outward edges of vi). Bθ(G) is simply constructed by vectorizing the updated node embeddings as [Bθ(G)]i,: = v′i.
Theorem 1. The existence and uniqueness of H∗ induced by Tθ(H [n];G). The proposed scheme for constructing Aθ(G), along with the bounded contraction factor 0 < γ < 1, is sufficient for Tθ(H [n];G) to be γ-contracting and, as a result, Tθ has unique fixed point H∗. Refer to Appendix A.1 for the proof.
Multi-head extension. Tθ can be considered as a graph convolution layer defined in equation 6. Thus, CGS can be easily extended to multiple convolutions in order to model a more complex iterative map. To achieve such multi-head extension with M graph convolutions, one can design fθ to produce a set of transition parameters [A1, ...,Am, ...,AM ] and a set of bias parameters [B1, ...,Bm, ...,BM ] for Tm(H [n]m ;G) ≜ γAmH [n]m +Bm for m = 1, ...,M .

heading: 4.2 COMPUTING FIXED POINTS
text: The fixed point H∗m of the constructed iterative map Tm(H [n] m ;G) ≜ γAmH [n]m + Bm satisfies H∗m = γAmH ∗ m +Bm for m = 1, ...,M . Due to the linearity, we can compute the fixed point of Tm via matrix inversion: H∗m = (I − γAm)−1Bm (9)
where I ∈ Rp×p is the identity matrix. The existence of (I − γAm)−1 is assured from the fact that Tm is contracting (see Appendix A.2). The matrix inversions can be found by applying various automatic differentiation tools while maintaining its differentiability. However, the computational complexity of the matrix inversion scales O(p3), which can limit this approach from being scaled to large scale problems favorably.
To scale CGS to larger graph inputs, we compute the fixed point of Tm by repeatedly applying the iterative map Tm starting from an arbitrary initial hidden state H [0]m ∈ Rp until the hidden embedding converges, i.e., limn→∞ H [n] m = H∗m.
One can choose a way to compute the fixed point between inversion and iterative methods depending on the size of the transition matrix Am and its sparsity because these factors can result in different computational speed and accuracy. In general, for small-sized problems, inversion methods can be favorable; while for large-sized problems, iterative methods are more efficient. (See Appendix D.6)

heading: 4.3 DECODING FIXED POINTS
text: The final step of CGS is to aggregate the fixed points of multiple iterative maps and decode the aggregated fixed points to produce Y ∗. The entire decoding step is given as follows:
H∗ = [H∗m|| ... ||H∗M ] (10) Y ∗ = gϕ(H
∗;G) (11) where H∗ is the aggregated fixed points, M is the number of heads, and gϕ(·) is the decoder which is analogous to the decoding function of the network-analytic problems (equation 3).

heading: 5 TRAINING CGS
text: To train CGS with gradient descent, we need to calculate the partial derivatives of the scalar-valued loss L with respect to the parameters of Tθ. To do so, we express the partial derivatives using chain rule taking H∗ as the intermediate variable
∂L ∂(·) = ∂L ∂H∗ ∂H∗ ∂(·) (12)
where (·) denotes the parameters of Aθ(G) or Bθ(G). Here, ∂L∂H∗ is readily computable via an automatic differentiation package. However, computing ∂H ∗
∂(·) is less straightforward since H ∗ and
(·) are implicitly related via equation 6. One possible option to compute the partial derivatives is to backpropagate through the iteration steps.
Although this approach can be easily employed using most automatic differentiation tools, it entails extensive memory usage. Instead, exploiting the stationarity of H∗, we can derive an analytical expression for the partial derivative using the implicit function theorem as follows:
∂H∗
∂(·) = −
( ∂g(H∗,A,B)
∂H∗
)−1 ∂g(H∗,A,B)
∂(·) (13)
where g(H,A,B) = H − (γAH + B). Here, we omit the input-dependency of A and B for notational brevity. We compute the inverse terms via an iterative method. This option allows one to train CGS with constant memory consumption over the iterative steps. The derivation of the partial derivatives and the software implementation of equation 13 are provided in Appendix A.3 and B respectively.

heading: 6 EXPERIMENTS
text: We first evaluate the performance of CGS for two types of network-analytic problems: (1) the stationary state of physical diffusion in networks where the true solutions can be computed from linear iterative maps, and (2) the state values via GVI where the true solutions can be calculated from non-linear iterative maps. We then assess the capabilities of CGS as a general GNN layer by applying CGS to solve several graph property prediction benchmarks problems.

heading: 6.1 PHYSICAL DIFFUSION IN NETWORKS
text: Diffusion of fluid, heat, and other physical quantities are omnipresent in the science and engineering applications. Mathematically, physical diffusion in networks (e.g. pipe/pore networks) is often described by a graph. The stationary state of the graph can be expressed as∑
j∈N (i)
kij(pi − pj) = 0, ∀vi ∈ V \ ∂(V), (14)
pi = p b i , ∀vi ∈ ∂(V), (15)
where pi and pj are the potentials at vi and vj respectively, kij is the conductance of the edge the connects the two nodes, and pbi is the prescribed potentials at the boundary nodes that belong to the set ∂(V). Equation 14 specializes to a particular diffusion problem according to how p is prescribed (e.g. pressure for fluid flow, and temperature for heat transfer).
Figure 4: Pore network graph
As an example of physical diffusion, we consider fluid flow in porous media – particularly, finding the fluid pressures of a pore network that is in equilibrium state (the solutions of equations 14 and 15). We model the pore network as a 3D graph whose nodes and edges correspond to pore chambers and throats respectively as shown in figure 4. We assume linear diffusion such that p (i.e., Y ∗) can be computed using a linear iterative map (Gostick et al., 2016).
We employ CGS to predict the equilibrium pressures Y ∗ inside pore networks G. The node features are the Cartesian coordinates, volume, diameter, and boundary indicator of its corresponding pore. The boundary pressure is also as a node feature if the node corresponds to a boundary pore. The edge features are the cylinder volume, diameter, and length of its corre-
sponding throat. We sample training graphs such that the graphs fit into 0.1 m3 cubes. The training graphs, which have 50–200 nodes, are then randomly generated (See Appendix C.1). We train CGS such that it minimizes the mean-squared error (MSE) between the predicted ones and Y ∗.
To investigate the effectiveness of the multi-head extension, we train CGS(4), CGS(8) and CGS(16), where CGS(m) denotes CGS with m heads. For the baseline models, we use implicit GNNs, IGNN, (Gu et al., 2020), SSE (Dai et al., 2018), and n-layer GNN models GNN(n). IGNN and SSE find the fixed points in the forward propagation step. We utilize the same GNN architecture as the encoders for all baselines except for SSE. Please refer to Appendix C.1 for the details about the data generation, network architectures, and training schemes.
All CGS models show better generalization capabilities than the baselines in predicting Y ∗ as shown in Figure 2, even though all models show similar prediction errors during training (See Appendix C.1.2). CGSs with higher m show superior prediction results for the test cases. This difference evinces that the use of multi-head extension (i.e., multiple linear iterative maps) is advantageous due to the increased expressivity. Also, when comparing CGS(8) and GNN(1), which utilize the same encoder architecture and thus has the same number of GNN layers as shown in Figure 3, CGS(8) shows better prediction performance. This is because the 1-hop aggregation cannot provide enough information to compute the equilibrium pressure. This result indicates that CGS successfully accommodates the long-range patterns in graphs without adopting additional graph convolution layers.

heading: 6.2 GRAPH VALUE ITERATION
text: We investigate the performance of CGS on graph value iteration (GVI) problems, where the iterative map (equation 5) is non-linear as explained in section 3.1. The goal of the experiments is to show that CGS can estimate the state values that are computed from the nonlinear iterative map accurately, even when using the set of learned linear iterative maps, as shown in figure 5.
We train three CGS models, CGS(16), CGS(32), and CGS(64), and two baseline models, SSE and IGNN, on randomly generated MDP graphs. Each MDP graph has ns nodes and each node has na edges (i.e. the MDP has ns distinct states and na possible actions from each state). We sample ns and na from the discrete uniform distributions with lower and upper bounds of (20, 50) and (5, 10), respectively. We evaluate the predictive performance of the trained CGS on randomly generated GVI problems to verify the generalization capability of CGS for different ns and na. We refer to Appendix C.2 for the details of the model architectures and training schemes.
Table 1 summarizes the evaluation results of CGS and the baseline models. We report the meanabsolute-percentage error (MAPE) between the predicted state values and their true values, and the accuracy between the derived and optimal policies (in %), following Deac et al. (2020).
All the CGS models show reliable value and policy predictions for the in-training cases (ns = 20) as well as for the out-of-training cases. In general, CGS with many heads shows better prediction results than the model with smaller number of heads. The two baseline models shows significant performance drops as ns and, especially, na increase. (i.e. the number of incoming edges per node becomes larger). The results show that CGS can solve the unseen problems well (i.e., generalization). From the results, we can conclude constructing the transition map adaptively by utilizing the input graph information is a crucial factor for achieving better generalization and higher prediction accuracy of CGS, as compared to other baseline models (SSE and IGNN) that uses fixed transition maps.
Ablation studies. We conduct various ablation studies and confirm the following: • Appendix D.1: the flexibility of fθ(·) design is beneficial to attain higher predictive performances, • Appendix D.2: the input-dependency of both Aθ(G) and Bθ(G) is essential for generalization, • Appendix D.3: the contraction factor γ balances the predictability and computational time, • Appendix D.4: the linear Tθ(·) has sufficient expressivity, compared to the nonlinear extension,
while being robust to the hyperparameters,
• Appendix D.5: a relatively small of training samples (≥ 512) is enough to attain higher predictive performances, compared to IGNN.
• Appendix D.6: the iterative fixed point computation scales better than the direct method in terms of memory usage.
From the results, we confirmed the proposed design of Tθ(·) is effective and robust to the hyperparameters in solving GVI , and the proposed training scheme – iterative fixed point computation and
computing the gradient via the implicit function theorem – is practically suitable in terms of memory usage. Due to the page limit, we refer to D for the details and results of the ablation studies.

heading: 6.3 GRAPH CLASSIFICATION
text: We show that CGS can also perform general graph classification tasks accurately, where the existence or the meaning of a fixed point is hard to be clearly defined, although it is originally designed to predict quantities related to the fixed points.
We assess the graph classification performance of CGS on six graph classification benchmarks: two social-network datasets (IMDB-Binary, IMDB-Multi) and four bioinformatics datasets (MUTAG, PROTEINS, PTC, NCI1). Since the social-network datasets do not have node features, they are generated based on the node degrees following Xu et al. (2018). Also, the edge features are initialized with one vectors for all datasets. To conduct the graph classification tasks, we perform the sum readout over the outputs of CGS (i.e., summing all fixed points), and then utilize additional MLP to predict the graph labels from the readout value. We perform 10-fold cross validation and report the average and standard deviation of its accuracy for each validation fold, following the evaluation scheme of Niepert et al. (2016). We refer to Appendix C.3.1 for the details of the network architecture, training, and hyperparamter searchings.
The results in Table 2 show that the classification performance of CGS is comparable to those of other methods using fixed point iteration (LP-GNN and IGNN). Note that we reproduced the results of IGNN as the original paper used a different performance metric (IGNN (ours)). We provide the additional benchmark results in Appendix C.3.2 to compare with IGNN following their test metric. In general, CGS shows better performance than LP-GNN and IGNN, which also find the fixed points of graph convolutions on the social-network datasets where the node features are not given. From the results, CGS can be interpreted as that CGS finds "virtual fixed points" that contain the most relevant information to classify the graph labels. These results indicate that CGS has a potential as an general graph convolution layer.

heading: 7 CONCLUSION
text: We propose the convergent graph solver (CGS) as a new learning-based iterative method to compute the stationary properties of network-analytic problems. CGS generates contracting input-dependent linear iterative maps, finds the fixed points of the maps, and finally decodes the fixed points to predict the solution of the network-analytic problems. Through various network-analytic problems, we show that CGS has competitive capabilities in predicting the outputs (properties) of complex target networked systems in comparison with the other GNNs. We also show that CGS can solve general graph benchmark problems effectively, showing the potential that CGS can be used as a general graph implicit layer for processing graph structured data.

heading: 8 ETHIC STATEMENTS AND REPRODUCIBILITY
text: Ethics statement We propose a deep learning method that learns iterative mappings of the target problems. As we discussed in section 1, various types of engineering, science, and societal problems are formulated in the form of fixed point-finding problems. On the bright side, we expect the proposed method to expedite scientific/engineering discoveries by serving as a fast simulation. On the other hand, as our method rooted in the idea of finding fixed points, it can be used to analyze networks and finding an adversarial or weak point of a network that can change the results of many network-analytic algorithms that are ranging from everyday usages, such as recommendation engines of commercial services, and maybe life-critical usages.
Reproducibility As machine learning researchers, we consider the reproducibility of numerical results as one of the top priorities. Thus, we put a significant amount of effort into pursuing the reproducibility of our experimental results. As such, we set and tracked the random seed used for our experiments and confirmed the experiments were reproducible.

heading: Convergent Graph Solvers

heading: Supplementary Material

heading: Table of Contents

heading: A Proofs and derivations 14
text: A.1 Existence of converged embedding . . . . . . . . . . . . . . . . . . . . . . . . 14 A.2 Existence of inverse matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.3 The complete derivation of the partial derivatives . . . . . . . . . . . . . . . . . 14

heading: B Software Implementation 16

heading: C Details on experiments 18
text: C.1 Physical diffusion experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 18 C.2 Graph value iteration experiments . . . . . . . . . . . . . . . . . . . . . . . . . 19 C.3 Graph classification experiments . . . . . . . . . . . . . . . . . . . . . . . . . 19

heading: D Ablation studies 21
text: D.1 Effect of fθ(·) architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 D.2 Design of the input-dependent transition maps . . . . . . . . . . . . . . . . . . 21 D.3 Effect of contraction parameter γ . . . . . . . . . . . . . . . . . . . . . . . . . 21 D.4 Comparisons to the non-linear iterative maps . . . . . . . . . . . . . . . . . . . 22 D.5 Sample efficiency of CGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 D.6 Runtime comparisons of direct inversion and iterative methods . . . . . . . . . . 23

heading: A PROOFS AND DERIVATIONS
text: In this section, we prove the existence of the converged hidden embedding H∗ and (I − γAθ)−1 and give the complement derivation of the partial derivatives that is used to train CGS.

heading: A.1 EXISTENCE OF CONVERGED EMBEDDING
text: The proposed transition map is defined as follows:
Tθ(H [n];G) = γAθ(G)H [n] +Bθ(G) (A.1)
where 0.0 < γ < 1.0, Aθ(G) ∈ Rp×p with ||Aθ(G)|| ≤ 1.0, and Bθ(G) ∈ Rp×q . We first show that the proposed transition map is contracting. Then, the existence of unique fixed point can be directly obtained by applying the Banach fixed point theorem (Banach, 1922).
Lemma 1. Tθ(H [n];G) is a γ-contraction mapping.
Proof. The proof is trivial. Consider the following equations:
||Tθ(H [n+1];G)− Tθ(H [n];G)|| = ||γAθ(G)H [n+1] +Bθ(G)− γAθ(G)H [n] −Bθ(G)|| = γ||Aθ(G)(H [n+1] −H [n])|| ≤ γ||Aθ(G)|| × ||(H [n+1] −H [n])|| = γ||(H [n+1] −H [n])||
The inequality holds by the property of the spectral norm. Therefore, Tθ(H [n];G) is γ-contracting.

heading: A.2 EXISTENCE OF INVERSE MATRIX
text: CGS can find the fixed point H∗ by pre-multiplying the inverse matrix of (I − γAθ) to Bθ, which is given as follows:
H∗ = (I − γAθ)−1Bθ (A.2)
The following lemma shows the existence of the inverse matrix. Lemma 2. With matrix Aθ, which is constructed using equation 6, and 0.0 < γ < 1.0, (I − γAθ) is invertible.
Proof. According to the invertible matrix theorem, (I − γAθ) being invertible is equivalent to (I − γAθ)x = 0 only having the trivial solution x = 0. By rewriting the equation as x = γAθx and defining the γ-contraction map F(x) = γAθx, we can get the following result:
lim n→∞
Fn(x) = 0, for any x. (A.3)
From the Banach fixed point theorem, we can conclude that the unique solution x is 0. Therefore, (I − γAθ) is invertible.

heading: A.3 THE COMPLETE DERIVATION OF THE PARTIAL DERIVATIVES
text: The following equality shows the relationship of the gradient of the scalar-valued loss L, the fixed point H∗, and the parameters of the fixed point equation Aθ(G), Bθ(G):
∂L ∂(·) = ∂L ∂H∗ ∂H∗ ∂(·) (A.4)
where (·) denotes Aθ(G) or Bθ(G). Here, ∂L∂H∗ are readily computable via automatic differentiation packages. However, computing ∂H ∗
∂(·) is less straightforward since H ∗ and (·) are implicitly related
via equation 6.
We reformulate the transition map (equation 6) in a root-finding form g(H,A,B) = H − (γAH + B). Here, we omit the input-dependency of A and B for the brevity of the notations. For the fixed point H∗, the following equality holds:
∂g(H∗(A,B),A,B)
∂A = 0 (A.5)
We can expand the derivative by using chain rule as follows:
∂g(H∗(A,B),A,B)
∂A =
∂g(H∗,A,B)
∂A +
∂g(H∗,A,B) ∂H∗ ∂H∗(A,B) A = 0 (A.6)
By rearranging the terms, we can get a closed-form expression of ∂H ∗
∂A as follows:
∂H∗
∂A = −
( ∂g(H∗,A,B)
∂H∗
)−1 ∂g(H∗,A,B)
∂A (A.7)
Here, we can compute ∂g(H ∗,A,B) ∂A easily by using either automatic differentiation tool or manual gradient calculation. However, in practice, directly computing (∂g(H ∗,A,B)
∂H∗ ) −1 can be problematic
since it involves the inversion of Jacobian matrix. Instead, we can (1) construct a linear system whose solution is (∂g(H
∗,A,B) ∂H∗ ) −1 and solve the system via some matrix decomposition or (2) solve another fixed point iteration which converges to ( g(H
∗,A,B) ∂H∗ ) −1. We provide the pseudocode that shows how to compute the inverse Jacobian with the second option in section B.
The partial derivative with respect to B can be computed through the similar procedure and it is given as follows:
∂H∗
∂B = −
( ∂g(H∗,A,B)
∂H∗
)−1 ∂g(H∗,A,B)
∂B (A.8)

heading: B SOFTWARE IMPLEMENTATION
text: In this section, we provide a Pytorch style pseudocode of CGS which computes the derivatives via the backward fixed point iteration.
1 import torch 2 import torch.nn as nn 3
4 5 class CGP(nn.Module): 6 """ 7 Convergent Graph Propagation 8 """ 9
10 def __init__(self, 11 gamma: float, 12 activation: str, 13 tol: float = 1e-6, 14 max_iter: int = 50): 15 16 super(CGP, self).__init__() 17 self.gamma = gamma 18 self.tol = tol 19 self.max_iter = max_iter 20 self.act = getattr(nn, activation)() 21 22 self.frd_itr = None # forward iteration steps 23 24 def forward(self, A, b): 25 """ 26 :param A: A matrix [#.heads x #. edges x #. edges] 27 :param b: b matrix [#.heads x #. nodes x 1] 28 :return: z: Fixed points [#. heads x #. nodes] 29 """ 30 31 z, self.frd_itr = self.solve_fp_eq(A, b, 32 self.gamma, 33 self.act, 34 self.max_iter, 35 self.tol) 36 37 # re-engage autograd and add the gradient hook 38 z = self.act(self.gamma * torch.bmm(A, z) + b) 39 40 if z.requires_grad: 41 y0 = self.gamma * torch.bmm(A, z) + b 42 y0 = y0.detach().requires_grad_() 43 z_next = self.act(y0).sum() 44 z_next.backward() 45 dphi = y0.grad 46 J = self.gamma * (dphi * A).transpose(2, 1) 47 48 def modify_grad(grad): 49 y, bwd_itr = self.solve_fp_eq(J, 50 grad, 51 1.0, 52 nn.Identity(), 53 self.max_iter, 54 self.tol) 55 56 return y 57 58 z.register_hook(modify_grad) 59 z = z.squeeze(dim=-1) # drop dummy dimension 60 return z
61 62 @staticmethod 63 @torch.no_grad() 64 def solve_fp_eq(A, b, 65 gamma: float, 66 act: nn.Module, 67 max_itr: int, 68 tol: float): 69 """ 70 Find the fixed point of x = act(gamma * A * x + b) 71 """ 72 73 x = torch.zeros_like(b, device=b.device) 74 itr = 0 75 while itr < max_itr: 76 x_next = act(gamma * torch.bmm(A, x) + b) 77 g = x - x_next 78 if torch.norm(g) < tol: 79 break 80 x = x_next 81 itr += 1 82 return x, itr
Listing 1: CGS pseudocode

heading: C DETAILS ON EXPERIMENTS
text: In this section, we provide the details of experiments including generation schemes and hyperparatemers of models. We run all experiments on a single desktop equipped with a NVIDIA Titan X GPU and AMD Threadripper 2990WX CPU.

heading: C.1 PHYSICAL DIFFUSION EXPERIMENTS

heading: C.1.1 DATA GENERATION
text: In this section, we provide the details of porous network problems.
Diffusion Equation For fluid flow in porous media described by Darcy’s law, equation 14 is specific to ∑
j∈N (i)
π
8µ r4ij lij (pi − pj) = 0, ∀vi ∈ V \ ∂(V), (A.9)
where µ is the dynamic viscosity of the fluid, and rij and lij are the radius and length of the cylindrical throat between the i-th and j-th pore chambers respectively.
Data generation We generate random pore networks inside a cubic domain of width 0.1m using Voronoi tessellation. We sample the pore diameters from the uniform distribution of U(9.9× 10−3 m, 10.1×10−3 m) and assume that the fluid is water with µ = 10−3 Nsm−2 under the temperature of 298K(25 °C). The boundary conditions are the atmospheric pressure (101, 325Pa) on the front surface of the cube, zero pressure on the back surface, and no-flux conditions on all other surfaces. We simulate the flow using OpenPNM (Gostick et al., 2016). We normalize the pressure values (targets) by dividing the pressure by the maximum pressures of the pores. Note that this normalization is always viable because the maximum pressure is the prescribed boundary pressure on the front surface due to the physical nature of the diffusion problem.

heading: C.1.2 DETAILS OF CGS AND BASELINES
text: In this section, we provide the details of CGS and the baseline models. For brevity, we refer an MLP with hidden neurons n1, n2, ... nl for each hidden layer as MLP(n1, n2, ..., nl).

heading: Network architectures
text: • CGS(m): fθ is a single layer attention variant of graph network (GN) (Battaglia et al., 2018) whose edge, attention, and node function are MLP(64). The output dimensions of the edge and node function are determined by the number of heads m. ρ(·) is the summation. gθ is MLP(64, 32). All hidden activations are LeakyReLU. We set γ as 0.5.
• IGNN: fθ, gθ is the same as the one of CGS(8).
• SSE: We modify the original SSE implementation (Dai et al., 2018) so that the model can take the edge feature as an additional input. As gθ, we use the same architecture to the one of CGS(m).
• GNN(n): It is the plain GNN architecture having the stacks of n different GNN layers as fθ. For each GNN layer, we utilize the same GN layer architecture to the one of CGS(m). gθ is the same as the one of CGS(m).
Training details We train all models with the Adam optimizer (Kingma & Ba, 2014), whose learning rate is initialized as 0.001 and scheduled by the cosine annealing method (Loshchilov & Hutter, 2016). The loss function is the mean-squared error (MSE) between the model predictions and the ground truth pressures. The training graphs were generated on-fly as specified in the data generation paragraph. We used 32 training graphs per gradient update. On every 32 gradient update, we sample the new training graph. We train 1000 gradient steps for all models.
Training curves The training curves of the CGS models and baselines are provided in figure 6.

heading: C.2 GRAPH VALUE ITERATION EXPERIMENTS
text: In this section, we provide the details of data generation, the CGS models and baseline architecture and their training schemes.

heading: C.2.1 DETAILS OF GVI DATA GENERATION
text: Data generation We generate the MDP graph by randomly sampling na out-warding edges for all nodes. The discount factor of MDP α is 0.9. Rewards are sampled from the uniform distribution whose upper and lower bounds are 1.0 and -1.0 respectively. The true state-values (labels) are computed by iteratively employing the exact analytical Bellman operator (equation 5) until the state-values converge (i.e., value iteration). The convergence tolerance is 0.001.

heading: C.2.2 DETAILS OF CGS AND BASELINES

heading: Network architectures
text: • CGS(m): fθ is a three layer attention GN network as whose edge, attention, and node function are MLP(128). The output dimensions of the edge and node function are determined by the number of heads m. ρ(·) is the summation. gθ is MLP(64, 32). All hidden activations are LeakyReLU. We set γ as 0.5.
• IGNN: fθ, gθ is the same as the one of CGS(32). • SSE: We modify the original SSE implementation (Dai et al., 2018) so that the model can take the
edge feature as an additional input. As gθ, we use the same architecture to the one of CGS(m).
Training details We train all models with the Adam optimizer whose learning rate is initialized as 0.001 and scheduled by the cosine annealing method. The loss function is MSE between the model predictions and the ground truth state-values. The training graphs are generated on-fly as specified in the data generation paragraph. We use 64 training graphs per gradient update. On every 32 gradient graphs, we sample the new training graph. We train 5000 gradient steps for all models.

heading: C.3 GRAPH CLASSIFICATION EXPERIMENTS

heading: C.3.1 EXPERIMENT DETAILS AND HYPERPARAMETERS
text: In this section, we explain the network architecture and training details of the six graph benchmark problems. Across all the benchmark dataset, we use the dataset implementation of DGL (Wang et al., 2019) and cross-validation indices generated with Scipy (Virtanen et al., 2020). We set the contraction factor γ as 0.5. We train all models with the Adam optimizer whose learning rate is initialized as 0.001 and scheduled by the cosine annealing method for 500 (100 for NCI1 dataset
due to the large datset size) epochs with 128 mini-batch size. We set the random seed of Scipy, Pytorch Paszke et al. (2019), and DGL as 0.
Hyperparameter tuning Due to our limited computational resources, we search at most 10 different pairs of hyperparameters for each dataset. To find the initial hyperparemters, we first tune CGS(4) on MUTAG, which is the smallest dataset. For each benchmark dataset, we start hyperparmeter tunings from the hyperparameters that are used for MUTAG-CGS(4) and tune the activation functions of fθ and gθ, the number of GN layers in fθ, the number of layers of gθ, and the dropout rate of gθ.

heading: Network architectures
text: • IMDB-Binary: fθ is a two layer attention GN network as whose edge, attention, and node function are MLP(128). The output dimensions of the edge and node function are determined by the number of heads m. The hidden dimensions of edge, and node are 64. ρ(·) is the summation. gθ is MLP(64, 32). All hidden activations are Swish Ramachandran et al. (2017).
• IMBD-Multi: The same as IMDB-Binary. • MUTAG: The same as IMDB-Binary. All hidden activations are LeakyReLU. • PROTEINS: The same as the MUTAG. Apply dropout Srivastava et al. (2014) with probability
0.2 after the activation functions of gθ. • NCI1: The same as IMDB-Binary.

heading: C.3.2 EXTENDED BENCHMARK RESULTS
text: We provide the benchmark results of CGS with the test metric that is used from (Gu et al., 2020). The test metric is to calculate the average maximal test accuracies over the test folds. CGS shows better predictive accuracies compared to IGNN.

heading: D ABLATION STUDIES
text: In this section, we provide the results of the ablation studies in GVI. The ablation studies were done to understand
• the effect of fθ(·) architecture, • the effect of transition map design, • the effect of contraction factor γ, • the effect of employing the non-linear transition maps, • the sample efficiency of CGS.
We also provide the runtime comparisons of direct inversion and iterative methods to find the fixed points of hidden embedding.
D.1 EFFECT OF fθ(·) ARCHITECTURE
As the Aθ(G) and Bθ(G) generation schemes of CGS allow the employment of the arbitrary architecture of GNN as fθ(·), CGS has different predictive performances depending on the architectural selection of fθ(·). Here, we investigate the effect of the number of GN layers in fθ(·) to the predictive performance of CGS for the GVI problems. The number of the independent GN layers controls the range of information when CGS constructs the transition maps. That is, a larger number of GN layers allows the information to be gathered from far neighborhoods while constructing Aθ(G) and Bθ(G). Table 4 shows the predictive performances of CGS models with the different number of GN layers in fθ(·). The model with n GN layers is referred to as n-CGS(16). In general, the model with a larger number of GN layers performs better. These results highlight that allowing the flexibility in fθ(·) can be practically beneficial when we derive the model with guaranteed convergence.

heading: D.2 DESIGN OF THE INPUT-DEPENDENT TRANSITION MAPS
text: The transition map Tθ(·) of CGS is given as follows:
Tθ(H [n];G) = γAθ(G)H [n] +Bθ(G) (A.10)
By restricting Aθ(G) and Bθ(G) to be input independent, we can consider simpler variants of CGS. We can consider three variants which has (1) input-independent A and B; (2) input-independent A and input-dependent B; (3) input-dependent A and input-independent B. Out of three variants, (1) cannot differentiate different G and (2) is similar to Gu et al. (2020). Therefore, we provide the GVI results of (3) - Fixed B CGS - in here. From Table 5, we can confirm that having input-dependent A and B in Tθ(·) shows best performances.
D.3 EFFECT OF CONTRACTION PARAMETER γ
The contraction parameter γ, as a hyperparameter, changes the fixed point (because it changes the linear map) and controls the rate of convergent speed as well. We investigate the effect of γ to the
predictive performance of CGS. We test three γ = 0.3, 0.5, 0.7 and the experimental results are given in Table 6. The CGS model with smaller γ tends to have higher prediction errors compared to the larger γ. However, using larger γ enlargers the number of iterative steps. In this regard, we use γ = 0.5 to balance the computational speed and the predictive performance of CGS.

heading: D.4 COMPARISONS TO THE NON-LINEAR ITERATIVE MAPS
text: A natural question to the linear iterative map of CGS is "can we achieve a performance gain if we employ non-linear contracting iterative maps?" To answer the question, we provide the extended experiment results.
The analysis of the existence and uniqueness of fixed points still holds when CGS employs componentwise non-expansive (CONE) activation (e.g., ReLU, LeakyReLU, Tanh, Swish, Mish) to the outputs of Tθ(·) given as follows:
Tθ(H [n];G) = ϕ(γAθ(G)H [n] +Bθ(G)) (A.11) where ϕ(·) is a CONE activation. The gradient of loss w.r.t Aθ(G) and Bθ(G) can be computed similarly to the non-linear activation cases as described in Appendix B.
We compare the GVI results of linear CGS to the non-linear CGS utilizing LeakyReLU as ϕ(·). Table 7 shows the GVI experiment results. The linear CGS, CGS(m), and non-linear CGS, nl-CGS(m), shows similar predictive performance on GVI experiments in general. From these experiments, we can observe that the non-linear extension of CGS does not give significant performance gain.
Furthermore, the training of non-linear CGSs can be challenging as (1) they exhibit higher variance in loss while training, and (2) the choice of non-linearity can severely change the performance of the entire model. For instance, the rectifying units such as ReLU and LeakyReLU can result in the premature termination of the iterative schemes of CGS when Bθ(G) has large negative values. Bounded activation such as Tanh limits the range of hidden embeddings to a certain range.

heading: D.5 SAMPLE EFFICIENCY OF CGS
text: We assumed the training graph and corresponding labels are easily sampled while we solving physical diffusion problem and GVI. For some practical cases, this assumptions cannot be made. Hence, we investigate the sample efficiency of CGS and our closest baseline IGNN(Gu et al., 2020).
We prepare 2048 training and 2048 test graphs and their labels. We then train CGS and IGNN by using the first 128, 256, 512, 1024, and 2048 training graphs. As shown in figure 7, we can confirm that the trade-off between the sample efficiency and expressivities. IGNN only utilizes the input-dependent bias terms in the transition maps; thus, such structural assumptions can serve as an effective regularizer when the training samples are limited. However, IGNN shows more minor improvements along with the number of training samples. On the other hand, CGS performs worse than IGNN when the training graphs are limited, but it starts to outperform IGNN as more training samples are used. Finally, when the number of training samples increases (≥ 512 graphs), CGS outperforms IGNN significantly.

heading: D.6 RUNTIME COMPARISONS OF DIRECT INVERSION AND ITERATIVE METHODS
text: In this paragraph, we provide the experimental results that shows the runtime of CGS(2) models which solve the fixed point equation via direct inversion and iterative methods. For all size of GVI graphs, we test the models 100 times with na = 5. As shown in figure 8, for small graphs ns ≤ 4000, solving the fixed point equation with the direct inversion is faster than solving it with the iterative scheme. However, the direct inversion scales worse than iterative method.

Review summaries:
This paper proposes a method for learning how to predict the properties of graphs. These properties are fixed points (e.g. the PageRank vector or Equilibrium Values of a Diffusion Process).Wrt existing methods (SSE and IGNN), CGS learn multiple transition matrices instead of using a unique one. The underlying idea is to: (a) learn several linear operators whose iterations lead to a fixed point, (b) aggregate/combine the fixed points and (c) make the final prediction. With these elements to hand it is possible to find fixed points for several tasks. The proposed approach is tested in PHYSICAL DIFFUSION, GRAPH VALUE ITERATION and GRAPH CLASSIFICATION. The results show that CGSs are competitive with the state of the art and sometimes generalize better. 

The paper introduces the equilibrium GNN-based model with a linear transition map. The transition map is made contracting to ensure that the fixed point exists and is unique. The paper provides the extensive comparisons and ablation studies for the proposed model. 

This paper proposes a graph neural network termed "Convergent Graph Solvers" that computes the stationary distribution of stochastic matrices. This is achieved by an implicit definition of a linear system and by computing the fixed point of that system. The paper shows on a range of problems that the proposed approach works well and that it is particularly fast for smaller graphs.