IGNN-Solver: A Graph Neural Solver for
Implicit Graph Neural Networks

Similar to traditional explicit GNNs, IGNNs [1, 2, 3, 4] also involve an aggregation process, with the distinction that the depth of layers (iteration step $k$) is infinite. The aggregation process in IGNNs is typically defined as

where $b_{\Omega}$ represents affine transformation parameterized by $\bm{\Omega}$, and the weight matrices ${\bm{W}}$ and $\bm{\Omega}$ are globally shared at each iteration step. The IGNN model $f_{\theta}$ is formally described by

where the representation, given as the “internal state” ${\bm{Z}}^{\star}$, is obtained as the fixed-point solution of the equilibrium (2). The final representation theoretically encompasses information from all neighbors in the graph. Therefore, IGNNs capture LRD within the graph, offering better performance compared to GNNs with finite iterations [1, 2]. Another notable advantage of this framework is its memory efficiency, as it only needs to retain the current state ${\bm{Z}}$ without requiring additional intermediate representations.

Notably, to obtain a valid representation for any given input ${\bm{X}}$ from IGNNs, we need to solve the unique internal state ${\bm{Z}}^{\star}$ in (3). However, there may not always exist a well-defined solution ${\bm{Z}}^{\star}$ for some inputs ${\bm{X}}$ [1]. To guarantee the existence and uniqueness of the solution to (3), we introduce the concept of well-posedness for IGNNs $f_{\theta}$ with activation function $\sigma$ [1]. Specifically, the weight matrix ${\bm{W}}\in$ $\mathbb{R}^{m\times m}$ and the adjacency matrix $\widehat{{\bm{A}}}\in\mathbb{R}^{n\times n}$ are said to be well-posed for $\sigma$ if the solution ${\bm{Z}}^{\star}\in\mathbb{R}^{m\times n}$ of  (3) exists and is unique for any $b_{\Omega}({\bm{X}})\in\mathbb{R}^{m\times n}$.

Based on the monotone operator theorem [39, 34], a sufficient condition for the well-posedness for IGNNs is defined as follows. Let the nonlinearity $\sigma$ be a non-expansive activation function (e.g., Sigmoid, Tanh, ReLU, etc.), and let $\widehat{{\bm{W}}}=\widehat{{\bm{A}}}^{\top}\otimes{\bm{W}}$. ¹¹1We represent the Kronecker product of matrices ${\bm{A}}$ and ${\bm{B}}$ as ${\bm{A}}\otimes{\bm{B}}$. Then, the IGNN model defined in (3) is well-posed as long as ${\bm{I}}-\widehat{{\bm{W}}}\succeq m{\bm{I}}$ for some $m>0$. ²²2 We represent ${\bm{A}}\succeq{\bm{B}}$ if ${\bm{A}}-{\bm{B}}$ is semi-positive definite. For the case of non-symmetric matrices $\widehat{{\bm{W}}}$, positive definiteness can be defined as the positive definiteness of the symmetric component, i.e., ${\bm{I}}-\widehat{{\bm{W}}}\succeq m{\bm{I}}$ if and only if ${\bm{I}}-\frac{\widehat{{\bm{W}}}+\widehat{{\bm{W}}}^{T}}{2}\succeq m{\bm{I}}$.

It establishes that if $\widehat{{\bm{W}}}$ is symmetric and ${\bm{I}}-\widehat{{\bm{W}}}\succeq m{\bm{I}}$ for some $m>0$, then the existence and uniqueness of the equilibrium point ${\bm{Z}}^{\star}$ are guaranteed, ensuring the model is well-posed. This follows from the sufficient condition that ${\bm{I}}-\widehat{{\bm{W}}}$ is strongly monotone.

To satisfy the well-posedness of IGNNs, a commonly used method [40, 3], is to orthogonally reparameterize ${\bm{W}}$ in the IGNN layer $f_{\theta}$ using the following Cayley map, which is inspired from the unitary operator [41], as

where $\kappa\in(0,1)$ and ${\bm{S}}={\bm{C}}-{\bm{C}}^{\top}$ is a skew-symmetric matrix with ${\bm{C}}$ as an arbitrary parameterized matrix. From (4), it follows that ${\bm{W}}$ is positive definite and its real eigenvalues are in $(0,1)$. Moreover, based on the definition of $\widehat{{\bm{A}}}$, we can derive the eigenvalues of $\widehat{{\bm{W}}}$ are in $(0,1)$, which ensures that ${\bm{I}}-\widehat{{\bm{W}}}\succeq m{\bm{I}}$, thereby ensuring the well-posedness of IGNN. We refer interested readers to [3] for detailed proof.

To accelerate the convergence of predicted values, we propose constructing an initializer

where $h_{\phi}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d^{\prime}}$, for a rapid and reasonable input-based initial estimate value, rather than simply setting the initial values to $\mathbf{0}$ or random, where $\phi$ are learnable parameters of the initializer. We set the intermediate dimension $d^{\prime}$ to be significantly smaller than the feature dimension $d$ to reduce the training overhead of the initializer.

In the original AA solver, the parameter $\beta$ is predetermined and fixed [13], whereas the parameters $\alpha^{[k]}$s are determined by solving numerous linear equations using the least squares method. Although this method operates adequately, its efficiency and adaptability in optimizing IGNN models are limited. Therefore, we propose introducing a tiny and learnable graph neural network as

where $\,s_{\xi}:(\mathbb{R}^{(m_{k}+1)\times d^{\prime}}\times\mathbb{R}^{n})\rightarrow(\mathbb{R}^{(m_{k}+1)}\times\mathbb{R}^{1})$, to predict the two parameters instead of setting them as the least squares solution on the past residuals ${\bm{G}}$.

Firstly, the number of edges in graph $\widehat{{\bm{A}}}$ is usually large in practice, as it is affected by the scale of the input (for example, in the ogbn-products dataset, the number of edges in the graph is close to 62M). To maintain the lightweight nature of tiny-GNN and reduce the computational cost of prediction, we propose introducing graph sparsification [43, 44, 45] for a light graph:

where $\widehat{{\bm{A}}}_{s}$ represents the sparse subgraph of $\widehat{{\bm{A}}}$, obtained via a graph sparsification operator $\varsigma$, with the parameter $\beta$ controls the sparsity of $\widehat{{\bm{A}}}_{s}$. We choose the RPI-Graph approach [17] in this paper, which is a plug-and-play graph sparsification method based on the principle of relevant information. It is worth noting that this subgraph is used exclusively in the tiny-GNN predictor $s_{\xi}$ and does not alter the original graph of the IGNN model $f_{\theta}$.

Besides, considering that $n$ (the number of nodes) is relatively large (for example, $n$ is approximately 169K in the ogbn-arxiv dataset), it is not appropriate to map $s_{\xi}$ from a high-dimensional space to a very low-dimensional one in (9) directly [25, 46, 47] owing to the curse of dimensionality [19]. Therefore, to maintain $s_{\xi}$ fast and compact, we recommend compress ${\bm{G}}^{[k]}_{i}$ to form a smaller yet still representative version $\bar{{\bm{G}}}^{[k]}_{i}$. Specifically, We map it from $\mathbb{R}^{(m_{k}+1)\times d^{\prime}}$ to an appropriate space $\mathbb{R}^{(m_{k}+1)\times p}$ by multiple layer perception (MLP). Then, the nearest $m_{k}+1$ sets of fixed points are merged into one feature matrix:

where $\mathop{||}_{i=k-m_{k}}^{k}$ represents the concatenation operation that stacks the compressed storage matrix $\bar{{\bm{G}}}^{[k]}_{i}\in\mathbb{R}^{(m_{k}+1)\times p}$ in the $k$-th iteration.

Once the smaller yet representative storage matrix $\bar{{\bm{G}}}$ is obtained as described above, we treat it as a $p$-channel matrix and employ a tiny graph neural network layer:

to predict the relative weights $\alpha^{[k]}$ assigned to past $m_{k}+1$ residual, as well as the AA mixing coefficient $\beta^{[k]}$ at the $k$-th iteration. Notably, we normalize the original output $\widehat{\alpha}^{[k]}$ using the function

to ensure that the constraint $\mathbf{1}^{\top}\alpha^{[k]}=1$ is satisfied.

In this way, $s_{\xi}$ shall autonomously learn and adjust these parameters $\alpha^{[k]}$ and $\beta^{[k]}$ based on previous solver steps, while also receiving gradients from subsequent iterations. We provide a comparison of different choices for $s_{\xi}$ in Section V-D, demonstrating that IGNN-Solver improves the convergence path. We also introduce the training procedure below.

The evaluation is conducted on nine different-field, real-world benchmarks for node classification tasks, including four citation datasets (Citeseer, ACM, CoraFull, ogbn-arxiv), three social interaction datasets (BlogCatalog, Flickr, Reddit), and two product network datasets (Amazon-all, ogbn-products); as well as five bioinformatics benchmarks (MUTAG, PTC_MR, COX2, PROTEINS, NCI1) for graph classification tasks. Further details can be found in Appendix A. Notably, to demonstrate the scalability of the proposed model on larger datasets, we use four large-scale datasets: Amazon-all, Reddit, ogbn-arxiv, and ogbn-products, two of which are from the Open Graph Benchmark (OGB) [22].

We compare IGNN-Solver with a comprehensive set of baselines in graph node classification tasks, including three implicit GNNs, i.e., MIGNN [3], EIGNN [2], IGNN [1], and nine explicit/traditional GNNs, i.e., AM-GCN [25], GCN [12], GAT [23], SGC [48], APPNP [49], JKNet [50], DEMO-Net [51], GCNII [52], ACM-GCN [53]. For the sake of fairness, we use the same hidden layer configuration for all baseline models. For example, in the Flickr dataset, the dimension of the hidden layers is uniformly set to $512$ for all methods (including ours, see in Table IV). In addition, for EIGNN, the arbitrary gamma is set to $0.8$, which is consistent with the original paper. For AM-GCN, the number of nearest neighbors in the KNN graph is set from $5,6,7$. For GCN, the hidden layer setup is the same as that for others. For GAT, three attention headers are used for each layer. For SGC, the power of self-loops in the graph adjacency matrix is set to $3$. For APPNP, we set the number of iterations to $3$ and teleport probability $\alpha$ to $0.5$. For JKNet, we set layers to $1$ in small datasets and set to $8$ in large ones. For DEMO-Net, the regularization parameter is set as $5e-4$, the learning rate is set as $5e-3$, and the hash dimension is set as $256$. For GCNII, we set $\alpha_{\ell}=0.1$ and $L_{2}$ regularization to $5e-4$, consistent with the original paper. For ACM-GCN, we set layers to $1$ in small datasets and set to $4$ in large ones. For our IGNN-Solver and its variants, we employ the same settings as the basic IGNN model. Furthermore, we adjust these parameters affecting the convergence and performance of IGNN through a combined approach of hierarchical grid search and manual tuning. The Adam optimizer [54] is used for optimization.

To demonstrate the superiority of the IGNN-Solver over IGNNs in terms of both performance and efficiency, we analyze the movement of the speed/accuracy Pareto curve across various datasets, rather than concentrating on a single point on the curve, which is depicted in Figures  2 and 3. All experiments are conducted five times, and the best results are reported. Furthermore, we present the node classification performance results for each method on various datasets, which are shown in Tables I and II. All experiments are conducted five times, and we report the average results along with the standard deviation. Additional details on the training procedure and the hyperparameters used for each task are provided in Appendices A and B.

From Figure 2 (the speed/accuracy Pareto curve on large-scale datasets: Amazon-all, Reddit, ogbn-arxiv and ogbn-products), we observe that:

1. 1.

   Comprehensive performance and efficiency: IGNN-Solver generally outperforms all other methods, particularly in large datasets with a greater graph radius. This advantage is more pronounced in such cases. The improved performance is attributed to IGNN-Solver’s significant enhancement in the convergence path and the better initialization of the fixed-point equation through the initializer;

2. 2.

   Rapid inference advantage: IGNN-Solver demonstrates a notable speed advantage during inference. For instance, on the large-scale dataset Amazon-all, our IGNN-Solver achieves $8\times$ faster inference speed than IGNN while maintaining the same performance, and there is consistently at least $1.5\times$ acceleration and a similar pattern has been observed as well on other datasets.

From Table I, we observe that implicit GNNs with infinite depth generally outperform shallow explicit GNNs in most cases. Furthermore, IGNN-Solver consistently shows higher accuracy percentages across most datasets compared to other state-of-the-art explicit GNNs like DEMO-Net [51], GCNII [52] etc, indicating its superior performance.

We further evaluate IGNN-Solver on several small-scale graph node classification tasks, including Citeseer, ACM, CoraFull, BlogCatalog, and Flickr. The hyperparameters used are outlined in Table IV, Section A-B. The mean accuracy and standard deviation are reported in Table II, while the speed/accuracy Pareto curve is presented in Figure 3. In general, learning LRD is not crucial for these tasks, as the graph diameters are relatively small [3]. However, as shown in Table II, even for these small-scale node classification tasks, IGNN-Solver outperforms IGNNs and many explicit GNNs, as well as other advanced implicit models. As illustrated in Figure 3, the solver continues to improve the convergence path of IGNN, even in small-scale tasks. These results, along with those in Section V-B, confirm the expressive power of the IGNN-Solver using learnable graph neural networks, even surpassing that of many explicit and enhanced implicit GNNs.

We extend our evaluation to graph classification tasks to further demonstrate the effectiveness and versatility of IGNN-Solver. We conduct experiments on several widely used bioinformatics benchmark datasets for graph classification, including MUTAG, PTC_MR, COX2, PROTEINS, and NCI1. These datasets vary in size and complexity, as detailed in Appendix A. We train the model using 10-fold cross-validation using the experimental setup of [1]. For each dataset, we employ standard 10-fold cross-validation to evaluate the test performance of the IGNN-Solver, with the baseline and their results drawn from [3, 62]. The average prediction accuracy and standard deviations are in Table III.

The results of the graph classification experiments, summarized in Table III, show that our IGNN-Solver almost outperforms the baseline models across all datasets, clearly demonstrating its ability to effectively capture graph structures and significantly improve classification accuracy. These findings confirm the effectiveness of the IGNN-Solver in graph classification tasks, further supporting its potential for deployment in real-world applications where graph data is prevalent.

To investigate the proportion of training overhead attributed to the solver throughout the entire training process, we depict in Figure 4 the percentage of training computations on IGNN-Solver relative to the total time overhead on IGNN. Our approach is not only effective but also incurs minimal training overhead for the solver: the solver module is relatively small and requires only about $1\%$ of the original training time required by the IGNN model. Besides, this proportion will be lower for large-scale data. For instance, on the Amazon dataset, IGNN necessitates a total runtime of $3$ hours, whereas the solver only requires approximately $1.6$ minutes, indicating that our solver maintains its lightweight nature across any-scale datasets.

We present additional evidence on the convergence and generalizability of the neural solvers in Figure 5, where the ordinate represents the relative residual $\frac{||f(z)-z||_{F}}{||z||_{F}}$, with $\|\cdot\|_{F}$ denoting the Frobenius norm. We compare the convergence of a pre-trained IGNN model under two conditions: 1) canonical iteration in IGNN [1], and 2) IGNN-Solver with $K=6$ unrolling steps.

From Figure 5, we observe that the IGNN-Solver, trained with $K$ unrolling steps, is capable of generalizing beyond $K$, with both solvers ultimately reaching a stable state, demonstrating good convergence. Moreover, thanks to the fixed-point solution of the graph neural parameterization, IGNN-Solver continuously enhances the convergence of the canonical iterators while being more computationally efficient. Specifically, each step of the neural solver is cheaper than the standard iterator step. This explains the observed improvement in inference efficiency of at least $1.5\times$ in Section V-B.

We provide additional evidence on the convergence and generalization of IGNN-Solver by exploring alternatives to our proposed graph neural (GN) layers, including neural network (NN), temporal convolutional network (TCN) [63], and gated recurrent unit (GRU) [64]. The relative residual curves on the Citeseer dataset, observed during both the warm-up and final inference stages, are presented in Figure 6. Interestingly, despite their larger number of parameters and more complex architectures, TCN and GRU consistently underperform compared to NN. Notably, IGNN-Solver (GN) significantly improves the convergence path of all solvers and exhibits the fastest rate of residual reduction, establishing itself as the most effective predictor in IGNN-Solver due to its ability to effectively exploit rich graph information.

In this section, we analyze the effect of different loss components on the IGNN-Solver. The performance of IGNN-Solver and its variants on the Citeseer dataset is shown in Figure 7, with similar patterns observed across other datasets and solvers. For the main convergence loss $\mathcal{L}_{\text{conv}}$ in fixed-point iterations (16), we design two contrasting schemes concerning its weight $\lambda_{1}^{[k]}$s: setting them to a constant value (i.e., $\lambda_{1}^{[k]}=\lambda_{1}$ for all $k$), or setting all values except the $K$-th term to zero (i.e., $\lambda_{1}^{[k]}=\lambda_{1}$ if $k=K$ else $0$). Both variant solvers still perform well, yet our suggested monotonically increasing scheme (i.e., emphasizing the later steps to a greater extent) exhibits the best performance.

Furthermore, removing either the initializer or the alpha loss $\mathcal{L}_{\alpha}$ negatively impacts performance, with the former having a more detrimental impact on the solver. Nonetheless, all these ablation configurations significantly outperform methods without a solver, underscoring the advantages of employing a custom learnable solver for the IGNN model.